Cu(100) surface alloy formation

Surface Science 459 (2000) 365–389 www.elsevier.nl/locate/susc

Atomistic modeling of Pd/Cu(100) surface alloy formation Jorge E. Garce´s a, Hugo O. Mosca b,c, Guillermo H. Bozzolo d,e, * a Centro Ato´mico Bariloche, Comisio´n Nacional de Energı´a Ato´mica, 8400 Bariloche, Argentina b Comisio´n Nacional de Energı´a Ato´mica, 1429 Buenos Aires, Argentina c Departamento de Fı´sica, Facultad de Ciencias Exactas y Naturales, Universidad de Buenos Aires, 1428 Buenos Aires, Argentina d Ohio Aerospace Institute, 22800 Cedar Point Rd., Cleveland, OH 44142, USA e National Aeronautics and Space Administration, Glenn Research Center, Cleveland, OH 44135, USA Received 21 June 1999; accepted for publication 10 March 2000

Abstract A straightforward modeling approach, using the Bozzolo–Ferrante–Smith (BFS) method for alloys, is introduced and applied to the formation process of Pd/Cu(100) surface alloys. Ranging from the deposition of one single Pd atom to the formation of a c(2×2) structure, the proposed methodology helps explain all the experimentally observed features for coverages of up to 0.5 ML. In excellent agreement with all the known experimental facts, the approach sheds light onto the exchange mechanism between adatoms and substrate atoms, interdiffusion to subsurface layers, the formation of ordered structures and Cu islands and the c(2×2) structure at 0.5 ML Pd coverage. © 2000 Elsevier Science B.V. All rights reserved. Keywords: Adatoms; Alloys; Compound formation; Computer simulations; Copper; Low index single crystal surfaces; Palladium; Semi-empirical models and model calculations; Single crystal surfaces; Surface defects; Surface structure, morphology, roughness, and topography

1. Introduction The rapid growth in surface analysis techniques has resulted in a correspondingly growing interest in the phenomenon of surface alloy formation [1,2]. In some systems, the observed structures are sufficiently simple to be accurately described both theoretically and experimentally [1], while in other cases, the available experimental techniques are not accurate enough to provide a detailed and definite description of the system at hand. This is particularly true in those cases where the surface alloys formed exhibit a very complex * Corresponding author. Fax: +1-440-962-3075. E-mail address: [email protected] (G.H. Bozzolo)

structure involving many layers below the surface, thus constituting a severe test for both experimental and theoretical methods. One such case is the Pd/Cu(100) surface alloy, where a rather complex behavior is observed with increasing Pd coverage [4–19]. The growth and structure of Pd surface alloys formed on Cu(100) have been the subject of numerous recent studies, which have resulted in an abundant wealth of information related to each step of the alloy formation process. Experimental confirmation of a Pd/Cu(100) surface alloy was obtained by low-energy electron diffraction (LEED) analysis by Wu et al. [3], which established the preference for the formation of a surface alloy over the formation of an ordered c(2×2) overlayer. This was followed by experimental [4– 12] and theoretical studies [13], which raised ques-

0039-6028/00/$ - see front matter © 2000 Elsevier Science B.V. All rights reserved. PII: S0 0 39 - 6 0 28 ( 00 ) 0 04 7 3 -8

366

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389

Fig. 1. (a) Pure, equilibrium crystal (reference atom denoted by the arrow), (b) a reference atom (denoted by the arrow) in the alloy to be studied (atoms of other species denoted with other shading) and (c) the same reference atom in a monatomic crystal, with the identical structure of the alloy to be studied, but with all the atoms of the same atomic species as the reference atom, for the calculation of the strain energy term for the reference atom. The strain energy is the difference in energy of the reference atom between (c) and (a).

Fig. 2. (a) Reference atom (denoted by an arrow) in the actual alloy environment and (b) the reference atom surrounded by a chemical environment equivalent to that in (a) but with the different neighboring atoms occupying equilibrium lattice sites corresponding to the ground state of the reference atom.

tions regarding the Pd coverage corresponding to the optimum c(2×2) structure [3,6–10] and the presence of Pd below the surface layer [6,10,14,15] as well as its influence on the formation of the c(2×2) structure. Moreover, complex structures involving atoms in the second layer appear at a higher Pd coverage [12,15,16 ]. A study of the annealing properties of Pd/Cu(100) surface alloys by Anderson et al. [17] first indicated the possibility of a loss of some Pd from the surface alloy and its replacement with Cu, resulting in a less ordered alloy surface or smaller islands. These results were interpreted as evidence for Pd dissolution into the bulk and the formation of an ordered

underlayer [17]. From a kinetics standpoint, the deposition of an overlayer of Pd at low temperatures, allowing it to alloy with the underlayer at higher temperatures [11], yields clear evidence for two mechanisms for diffusion: ordinary diffusion, resulting in the formation of Pd islands (at low temperature), and exchange jumps, resulting in c(2×2) nuclei. Better c(2×2) ordering can thus be obtained by increasing the heating rate. Recent scanning tunneling microscopy (STM ) studies by Murray et al. [4,5] have provided a firm experimental foundation for all the steps in the surface alloy formation process with increasing Pd coverage, addressing issues like the initial nucleation and growth leading to the formation of the c(2×2) alloy [4], including the mechanism by which Pd goes to subsurface layers, as well as the formation of the c(2×2)-p4g second-layer phase [5]. In spite of the lack of consensus on some of the details concerning the structure of Pd/Cu(100), it constitutes a rich example of surface alloying. It contains several features that, from a theoretical standpoint, constitute a perfect workbench for testing the validity and scope of theoretical techniques. The purpose of this work is to introduce a simple methodology designed to provide detailed insight on the different aspects of the surface alloying process, without any constraint that could arise from any predetermined knowledge of the system under study. This is accomplished by defin-

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389

367

Table 1 Experimental results used for the determination of the ECT and BFS parametersa Experimental results

Pd Cu

ECT parameters

˚) Lattice parameter (A

Cohesive energy (eV )

Bulk modulus (GPa)

p

˚ −1) a (A

˚ −1) l (A

˚) l (A

3.890 3.615

3.94 3.50

195.77 142.08

8 6

3.612 2.935

0.666 0.765

0.237 0.272

BFS parameters D =−0.0495 PdCu

D =−0.0431 CuPd

a The last four columns display the resulting Equivalent Crystal Theory (ECT ) [27] parameters: p is related to the principal quantum number, n, for the atomic species considered ( p=2n−2), a parameterizes the electron density in the overlap region between two neighboring atoms, l is a screening factor for atoms at distances greater than nearest-neighbor distance, and l is a scaling length needed to fit the lattice parameter dependence of the energy of formation with the universal binding energy relationship of Rose et al. [28]. The last row displays the BFS parameters D and D used in this work. PdCu CuPd

Fig. 3. Schematic representation of the BFS contributions to the total energy of formation. The left-hand side represents the reference atom (denoted by an arrow) in an alloy. The different terms on the right-hand side indicate the strain energy (atoms in their actual positions but of the same atomic species as the reference atom), the chemical energy term (atoms in ideal lattice sites) and the reference chemical energy (same as before, but with the atoms retaining the original identity of the reference atoms).

ing a simple operating procedure based on the creation of ‘catalogues’ of possible configurations (i.e. atomic distributions in a given computational cell ) suited to answer specific questions regarding the process of surface alloy formation. The energy analysis of the configurations included in these catalogues is then performed by the Bozzolo– Ferrante–Smith (BFS ) method for alloys [20,21], which is particularly suited for treating problems of surface structure and composition. As the results obtained warrant drawing specific conclusions, we single out features of the system for which corroboration with available experimental evidence can be

made. If successful, we will have laid the necessary framework for similar applications to other systems for which no previous knowledge exists, or for those cases where experimental results are not conclusive enough. As an initial application of this procedure, in this work, we concentrate on the early stages of Pd/Cu(100) surface alloy formation, comparing our results with those found experimentally. At higher coverages, a much more complex surface alloy forms, for which there is still some ambiguity between the experimental interpretation and the existing structural theoretical modeling [5,7,8,10].

368

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389

Fig. 4. Side view of a Cu slab with a (100) surface, showing different stages of interdiffusion of a Pd adatom: (a) the Pd adatom in the overlayer [Pd(O)], (b) occupying a surface site with the ejected Cu atom at nearest-neighbor distance ([Pd(S )+Cu(O)] ), (c) n same, but with the Cu(O) atom decoupled from the Pd(S) atom ([Pd(S )+Cu(O)] ), (d ) occupying a site one layer below the surface f layer and the Cu atom at second neighbor distance ([Pd(1b)+Cu(O)] ), (e) with the Cu atom decoupled from Pd(1b) n ([Pd(S)+Cu(O)] ) and (f ) with Pd in the second layer below the sur-face and the Cu atom in the overlayer ([Pd(2b)+Cu(O)] ). f f

As such, we will limit our analysis only to those coverages for which there is ample experimental evidence that provides a solid framework for our theoretical descriptions.

2. BFS method Since its inception a few years ago, the BFS method has been applied to a variety of problems, ranging from bulk properties of solid solution fcc alloys [20] and the defect structure in ordered bcc alloys [22,23] to more specific applications including detailed studies of the structure and composition of alloy surfaces and surface alloys [24]. In what follows, we provide a brief description of the operational equations of BFS. The reader

is encouraged to seek further details in previous papers where a detailed presentation of the foundation of the method, its basis in perturbation theory and a discussion of the approximations made are clearly shown [20–26 ]. The BFS method provides a simple algorithm for the calculation of the energy of formation of an arbitrary alloy (the difference between the energy of the alloy and that of its individual constituents). In BFS, the energy of formation DH is written as the superposition of elemental contributions e of all the atoms in the alloy i DH=∑ (E∞−E )=∑ e (1) i i i i i where E ∞ is the energy of atom i in the alloy, and i E is the corresponding value in a pure equilibrium i

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389

369

computing the difference, e , for each atom in the i alloy, a two-step approach is introduced for such a calculation in order to identify contributions to the energy due to structural and compositional effects. Therefore, e is broken up into two separate i contributions: a strain energy (eS) and a chemical i energy (e =eC−eC0). While there is a certain level i i i of arbitrariness in how this separation is implemented, it is only meaningful when a good representation of the actual process is obtained by properly linking both contributions. This is achieved by recoupling the strain and chemical contributions by means of a coupling function, g , i properly defined to provide the correct asymptotic behavior of the individual components. Each individual contribution, e , can therefore be written as i Fig. 5. Energy spectrum of the configurations shown in Fig. 4. DE is the difference between the energy of formation per adatom (in eV/atom) of a configuration and the lowest energy state.

monatomic crystal. In principle, the calculation of DH would simply imply computing the energy of each atom in its equilibrium pure crystal and then its energy in the alloy. In BFS, beyond directly

e =eS +g(eC −eC0 ). i i i i

(2)

The BFS strain energy contribution, eS , is i defined as the contribution to the energy of formation from an atom in an alloy computed as if all the surrounding atoms were of the same atomic species, while maintaining the original structure of the alloy. To visualize this concept, Fig. 1a represents the atom in question (identified with an

Fig. 6. Cu(100) surface with a step (denoted with slightly darker disks) and a Pd adatom (a) occupying a site in the overlayer [Pd(O)], (b) inserted in a surface site with the ejected Cu atom at a nearest-neighbor site on the overlayer ([Pd(S )+Cu(O)] ), (c) with the n ejected Cu atom in a site away both from the Pd(S ) atom and the step and (d) attached to the step.

370

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389

Fig. 7. Configurations with two Pd atoms corresponding to 0.033 ML Pd coverage. Pd atoms are denoted with black disks, whether they occupy overlayer or surface sites. Ejected Cu atoms are denoted with dark gray disks, while surface Cu atoms are indicated with light gray disks. The subindices denote the relative location between the last two listed atoms. The relative location between the last atom in the first term and the last atom in the second term of the expression is denoted with a subindex for the whole expression. (a–l ) Top view of a Cu(100) surface showing Pd and Cu atoms in the surface and overlayer. (m) Side view of a Cu(100) slab, showing a Pd atom in the first layer below the surface and one in the surface layer, with the ejected Cu atom occupying an overlayer site.

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389

Fig. 8. Energy spectrum of the N=2 configurations shown in Fig. 7. The energies DE (in eV/atom) are relative to the lowest energy state.

arrow) in an equilibrium position in its groundstate crystal (arbitrarily represented by a simple cubic lattice). Fig. 1b shows the same atom in the alloy being studied (also arbitrarily represented by a different crystallographic symmetry). Two things can be different between the reference crystal and the alloy. First, atoms of other species may occupy neighboring sites in the crystal and, second, the crystal lattice may not be equivalent in size or structure to that of the ground-state crystal of the reference atom. In Fig. 1b, the different atomic species are denoted with different symbols from that used for the reference atom, and the differences in size and/or structure are denoted with a schematically different atomic distribution as compared to the ground-state crystal shown in Fig. 1a. The BFS strain energy accounts for the change in energy due only to the change in geometrical environment of the crystal lattice (from Fig. 1a to Fig. 1b), ignoring the additional degree of freedom introduced by the varying atomic species in the alloy. In this context, Fig. 1c shows the environment ‘seen’ by the reference atom when computing its BFS strain energy contribution. The neighboring atoms conserve the sites in the actual alloy (Fig. 1b), but their chemical identity has changed to that of the reference atom (Fig. 1b), thus simpli-

371

fying the calculation to that of a single-element crystal. The BFS strain energy term represents the change in energy of the reference atom in going from the configuration denoted in Fig. 1a–c. In this sense, the BFS strain energy differs from the commonly defined strain energy in that the actual chemical environment is replaced by that of a monatomic crystal. Its calculation is then straightforward, even amenable to first-principles techniques. In our work, we use Equivalent Crystal Theory ( ECT ) [26,27] for its computation, due to its proven ability to provide accurate and computationally economical answers to most general situations. In all cases considered in this work, a rigorous application of ECT is reduced to that of its two leading terms, which describe average density contributions and bond-compression anisotropies. We neglect the three- and four-body terms dealing with the bond angle and face-diagonal anisotropies. The chemical environment of atom i is considered in the computation of eC , the first term in the i total BFS chemical energy contribution, where the surrounding atoms maintain their identity but are forced to occupy equilibrium lattice sites corresponding to the reference atom i. Following the convention introduced in Fig. 1, Fig. 2a shows the reference atom in the actual alloy (similar to Fig. 1b), while Fig. 2b indicates the atomic distribution used in computing the BFS chemical energy, eC (note that the lattice used in Fig. 2b i corresponds to that of the ground-state crystal of the reference atom, as shown in Fig. 1a). The total BFS chemical energy is then the difference between the energy of the reference atom in Fig. 2b, eC , i and its energy in its ground-state crystal (Fig. 1a). Building on the concepts of ECT, a straightforward approach for the calculation of the chemical energy is defined, properly parameterizing the interaction between dissimilar atoms. The second contribution to the chemical energy ( Fig. 1a) is included in order to completely free the chemical energy from structural defects, taking into account the possibility that the reference atom is not in a full-coordination environment. This is accomplished by recomputing the contribution, eC , i defined before, but once again assuming that all atoms are of the same species as the reference atom.

372

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389

Fig. 9. Configurations with three Pd atoms corresponding to 0.05 ML Pd coverage. Pd atoms are denoted with black disks, whether they occupy overlayer or surface sites. Ejected Cu atoms are denoted with dark gray disks and surface Cu atoms are indicated with light gray disks. The symbol Cu(i) is used to denote a Cu island. (a–l.1) Top view of a Cu(100) surface showing Pd and Cu atoms in the surface and overlayer. ( l.2) Side view of a Cu(100) slab, showing a Pd atom in the second layer below the surface and one in the surface layer, with the ejected Cu atom occupying an overlayer site.

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389

373

next-nearest neighbors, respectively, and where p, l, a and l are ECT parameters that describe element i (see Ref. [24] for definitions and details), r denotes the distance between the reference atom and its neighbors, S(r) describes a screening function [27] and the sum runs over nearest- and nextnearest neighbors. This equation determines the lattice parameter of a perfect equivalent crystal where the reference atom i has the same energy as it has in the geometrical environment of the alloy under study. R and R denote the nearest- and 1 2 next-nearest-neighbor distances in this equivalent crystal. Once the lattice parameter of the (strain) equivalent crystal, aS, is determined, the BFS strain energy contribution is computed using the universal binding energy relation of Rose et al. [28], which contains all the relevant information concerning a single-component system: Fig. 10. Energy spectrum of the N=3 configurations shown in Fig. 9. The energies DE (in eV/atom) are relative to the lowest energy state.

As mentioned above, the BFS strain and chemical energy contributions take into account different effects, i.e. geometry and composition, computing them as isolated effects. A coupling function, g , i restores the relationship between the two terms. This factor is defined in such a way as to properly consider the asymptotic behavior of the chemical energy, where chemical effects are negligible for large separations between dissimilar atoms. Within the framework of this discussion, the total BFS contribution, e , of each atom in the alloy can be i graphically depicted by the combination of strain and chemical effects shown in Fig. 3. In what follows, we provide the basic operational equations needed to compute each one of the terms introduced above. The BFS strain energy contribution, eS , is obtained by solving the ECT i perturbation equation NRpi e−aiR1 +MRpi e−(ai+1/li)R2 =∑ rpi e−(ai+S(rj))rj j 2 1 j (3) where N and M are the number of nearest- and

eS =Ei (1−(1+aS1 )e−aS1i ) (4) i C i where Ei is the cohesive energy of atom i and C where the scaled lattice parameter aS1 is given by i aS1 =q(aS −ai )/l (5) i i e i where q is the ratio between the equilibrium Wigner–Seitz radius and the equilibrium lattice parameter ai . e The BFS chemical energy is obtained by a similar procedure. As opposed to the strain energy term, the surrounding atoms retain their chemical identity, but are forced to be in equilibrium lattice sites of an equilibrium (otherwise monatomic) crystal i. The BFS equation for the chemical energy is given by NRpi e−aiR1 +MRpi e−(ai+1/li)R2 2 1 =∑ (N rpi e−aikr1 +M rpi e−(aik+1/li)r2 ) (6) ik 1 ik 2 k where N and M are the number of nearest- and ik ik next-nearest neighbors of species k of atom i. The chemical environment surrounding atom i is reflected in the parameters a , given by ik a =a +D (7) ik i ki where the BFS parameters D (a perturbation on the single-element ECT parameter a ) describe the i

374

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389

Fig. 11. Configurations with four Pd atoms corresponding to 0.067 ML Pd coverage. Pd atoms are denoted with black disks, whether they occupy overlayer or surface sites. Ejected Cu atoms are denoted with dark gray disks, and surface Cu atoms are indicated with light gray disks. Cu(i) denotes a Cu island, and 3Pd(2×2) denotes a set of three Cu atoms in surface sites following a c(2×2) pattern. (a–l.1) Top view of a Cu(100) surface showing Pd and Cu atoms in the surface and overlayer. ( l.2) Side view of a Cu(100) slab, showing a Pd atom in the second layer below the surface and one in the surface layer, with the ejected Cu atom occupying an overlayer site. (m–p) Top view of a Cu(100) surface showing Pd and Cu atoms in the surface and overlayer. (q–t) Top and side views of a Cu(100) slab, showing a Pd atom in the layers below the surface.

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389

375

Fig. 11. (continued ).

changes of the wave function in the overlap region between atoms i and k. Once Eq. (6) is solved for the equivalent chemical lattice parameter aC , the i BFS chemical energy is then (8) eC =c Ei (1−(1+aC1 )e−aCi * ) i i C i where g =1 if aC1> 0 and g =−1 if aC1<0. The i i i i scaled chemical lattice parameter is given by aC1 =q(aC −ai )/l . (9) i i e i Finally, as mentioned above, the BFS chemical

and strain energy contributions are linked by a coupling function, g , which describes the influence i of the geometrical distribution of the surrounding atoms in relation to the chemical effects and is given by g =e−aSi * (10) i where the scaled lattice parameter, aS1 , is defined i in Eq. (5). In this work, we used the BFS interaction parameters, D, determined following the procedure

376

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389

Fig. 12. Energy spectrum of the N=4 configurations shown in Fig. 11. The energies DE (in eV/atom) are relative to the lowest energy state.

outlined in Ref. [20]. The pure element parameters a , E , l, a, l and the BFS parameters D and e c PdCu D used in this study are listed in Table 1. CuPd 3. Results and discussion The objective of this paper is to introduce the BFS-based methodology for a detailed study of the most important features of the Pd/Cu(100) surface alloy formation. The methodology assumes no a priori information on the system at hand, and as such, it is designed to provide such information based solely on the results obtained. The only input necessary consists of the basic parameterization of the participating elements and lattice structures needed. We will begin this section by listing the most salient features of the Pd/Cu(100) system based

on recent experimental data [5,29,30], with the sole purpose of tracing a path for our discussion of the results within the framework of what is known experimentally. It should be noted that none of this information is used in the application of the method, and its presence at the beginning of this section serves the goal of orienting and organizing the results as they appear. The characteristic features of the low coverage Pd/Cu(100) surface alloy are: (f1) Initial growth proceeds via the formation of chains of Pd atoms alloyed into the surface along the [010] and [001] directions, with the Pd atoms preferentially occupying second neighbor sites, as opposed to a random distribution. (f2) No evidence for Pd diffusion in the surface, indicating their stability in the Cu substrate at room temperature. (f3) Ejected Cu atoms form islands on the surface or migrate to a nearby step, and as the Pd coverage increases, the mobility of Cu atoms is lowered, as manifested by the nucleation of new islands in areas between existing ones. (f4) Cu islands are observed to nucleate on top of alloyed areas, resulting in subsurface Pd. (f5) Formation of a c(2×2) structure as alloyed chains converge with increasing coverage. (f6) The long-range order of the c(2×2) structure is determined by an interplay between defects within the c(2×2) at lower coverages, and the initiation of second-layer growth at increased coverages, before the c(2×2) structure is fully completed. (f7) Coalescence of well-ordered Cu islands above the c(2×2) alloy surface [29]. (f8) Interdiffusion of Pd atoms to the second layer [30]. In what follows, and ignoring the previously listed knowledge, we will implement the BFS-based methodology and correlate the conclusions drawn from the results with the above-mentioned features. Starting with a single Pd atom deposited on a Cu(100) surface, this section is devoted to the analysis of the evolution of the system with increasing Pd coverage. In all cases, we introduce a stepby-step methodology, drawing conclusions for the behavior of the system by examining an appropri-

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389

377

Fig. 13. Configurations with five Pd atoms corresponding to 0.083 ML Pd coverage. Cu(i) denotes a Cu island in the overlayer and nPd(2×2) denotes a set of n Pd atoms located in surface sites following a c(2×2) pattern. (a–k.1) Top view of a Cu(100) surface showing Pd and Cu atoms in the surface and overlayer. (k.2) Side view of a Cu(100) slab, showing a Pd atom in the second layer below the surface and one in the surface layer, with the ejected Cu atom occupying an overlayer site.

378

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389

Fig. 14. Energy spectrum of the N=5 configurations shown in Fig. 13. The energies DE (in eV/atom) are relative to the lowest energy state.

ately defined set of possible configurations, which include all the necessary atomic distributions for addressing specific issues. These configurations consist of cells with 60 surface atoms and several planes in the direction perpendicular to the surface. We then compute the BFS energy of these states and plot them in the form of an energy spectrum. While it is true that the system will most likely reach those states of lowest energy, it is also true that other metastable states (close in energy to the lowest energy state) have a role in determining the behavior of the system, i.e. the closer these states are to the lowest energy state, the greater the likelihood that these states have of appearing in the actual system. As will be seen later, this approach will enable us to have a deeper understanding of the system.

Obviously, the first question is whether Pd atoms deposited on Cu(100) do or do not penetrate in the surface layer, and if they do, what happens to the ejected Cu atoms. Fig. 4 shows the set of configurations, designed to study the different options for the Pd and Cu atom and the ability of the Pd atom to interdiffuse in the Cu slab. The configurations include the Pd atoms in the overlayer (O), in the surface (S), or in the first (1b) or second (2b) layer below the surface plane, with the ejected Cu atom either close ( Fig. 4b and d) or far from the Pd atom. For all the examples studied in this paper, it is useful to visualize the results in the form of an energy spectrum, with each energy level in the spectrum indicating the energy of a given configuration, as shown in Fig. 5. Each level is labeled by its difference in energy, DE, with respect to the lowest state (the most likely to be found experimentally) and by a shorthand notation indicating the type of structure: X(O) denotes an X atom in the overlayer, X(S) in the surface plane, X(1b) one plane below the surface, X(2b) two planes below, etc. A subindex (n or f ) denotes whether the ejected Cu atom remains close (n) or far ( f ) from the Pd atom where it first inserts itself in the surface plane. A subindex [hkl ] indicates the direction on the (100) surface chosen by the atoms preceding the subindex. Fig. 5 clearly indicates a strong preference for a [Pd(S)+Cu(O)] state, as measured by the gap n between this and any other state in the spectrum. At this level of coverage, the position of the Pd(1b) or the Pd(2b) state indicates that it is highly unlikely that Pd will penetrate the Cu slab beyond the surface plane. The tightly ‘bound’ ground state ( lowest energy) hints at the possibility that this particular arrangement could constitute the building block of more complicated structures during the growth process. We conclude the analysis for N=1 (where N is the number of Pd atoms included in the configurations) by commenting on the behavior of the ejected Cu atoms by Pd substitutions. While we do not account for diffusion (i.e. all calculations are static, with atoms located at fixed sites), it is still possible to understand the behavior of Cu atoms and the ability to diffuse along the surface

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389

379

Fig. 15. Configurations with six Pd atoms corresponding to 0.1 ML Pd coverage. (a–c) Top view of a Cu(100) surface showing Pd and Cu atoms in the surface and overlayer. Cu(i) denotes a Cu island in the overlayer, and nPd(2×2) denotes a set of n Pd atoms located in surface sites following a c(2×2) pattern. (d, e) Side view of a Cu(100) slab, showing (d ) one Pd atom in the first layer below the surface and (b) one Pd atom two layers below the surface.

plane, by generating a small catalogue of possible configurations that schematically account for what is likely to be found in reality. This catalogue, shown in Fig. 6, includes configurations with a step on the surface (dark gray disks on the left

Fig. 16. Energy spectrum of the N=6 configurations shown in Fig. 15. The energies DE (in eV/atom) are relative to the lowest energy state. Only the lowest-lying energy levels are shown. The highest energy state (not shown) corresponds to all six Pd atoms in the overlayer, far from each other, with an energy gap of 1.27 eV/atom.

side of the figures), so that the ejected Cu atom has the possibility of migrating to a nearby step. Fig. 6a shows a single Pd atom in the overlayer, occupying a fourfold hollow site. Fig. 6b–d shows the Pd atom in a surface site and the displaced Cu atom in overlayer sites at the nearest-neighbor (NN ) distance of Pd ( Fig. 6b), somewhere between the Pd atom and the step (Fig. 6c) and attached to the Cu step ( Fig. 6d ). The last three configurations are also labeled with the gain in energy with respect to the configuration shown in Fig. 6a. Any configuration with the Pd atom in the surface plane [noted from now on as Pd(S)] lowers the energy. Moreover, attachment of the ejected Cu atom to either the step or the Pd(S) atom is, in that order, energetically favorable. This result is consistent with the experimental observation (f3) that Cu atoms tend to migrate to nearby islands at low coverage (note that with the limited size of the cell chosen for our calculations, one Pd atom corresponds to 0.017 ML Pd coverage) [5]. The fact that the migration of Cu atoms to

380

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389

Fig. 17. Configurations with seven Pd atoms corresponding to 0.117 ML Pd coverage. Cu(i) denotes a Cu island in the overlayer, and nPd(2×2) denotes a set of n Pd atoms located in surface sites following a c(2×2) pattern. (a, b) Top and side view of a Cu(100) surface showing Pd and Cu atoms in the overlayer, surface plane, and planes below the surface. (c, d ) Top view of a Cu(100) slab.

nearby steps is favored in this particular system could result in a substantial reduction in the number of configurations to be included in catalogues for higher coverage. However, the opposite is true if we were to attach to the catalogues the necessary generality so that they could be easily used for the analysis of any other system. Dealing with the presence of steps introduces, from the modeling standpoint, the need to deal with the specific preparation of the surface, as the presence of steps and terraces and their influence on the surface alloy formation are direct consequences of such a process. To avoid these unnecessary complications, and in order to understand the type of Cu growth in the overlayer with increasing Pd coverage, we will ignore the presence of steps in what follows and perform our calculations, assuming that Cu atoms do not diffuse on the Cu surface. However, the analysis of higher coverage cases, as will be discussed below, will have to be performed within the framework provided by this choice. This can be seen by examining the case of two Pd atoms. A substantial set of possible configurations is shown in Fig. 7, with the corresponding energy spectrum (once again, the energies shown

are referenced to the lowest energy state) shown in Fig. 8. The lowest energy state corresponds to one where both Pd atoms locate themselves initiating a chain in surface sites along the [010] direction, and the ejected Cu atoms closely attached to them [experimental fact (f1)]. The preference for a specific direction is clearly seen by comparing this state with the next one up in energy, where the only structural difference is the two Pd(S) atoms locating themselves at NN sites. Moreover, the structure of this ground state can be seen as a particular coupling of two N=1 states (as in Fig. 4b). We will see that this coupling scheme gives rise to the formation of the c(2×2) structure at higher coverages [experimental fact (f5)]. Also, as seen for N=1, the Pd(1b) state is too high in energy and therefore is unlikely to be found. The exchange mechanism by which Pd atoms choose surface sites and the way it develops for higher coverage situations can be partially explained in terms of the strain effective coordination, i.e. the environment ‘seen’ by the Pd atom in the host lattice from a BFS standpoint. We define the strain effective coordination of an atom X in a host lattice Y as the number of nearest neighbors

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389

Fig. 18. Energy spectrum of the N=7 configurations shown in Fig. 17. The energies DE (in eV/atom) are relative to the lowest energy state. Only the lowest-lying energy levels are shown. The highest energy state (not shown) corresponds to all seven Pd atoms in the overlayer, far from each other, with an energy gap of 1.26 eV/atom.

of the X atom in Y sites, needed to simulate the equilibrium conditions of an X atom in a pure crystal of its own species [i.e. the X atom has zero BFS strain energy, as defined in Eq. (4)]. As described by BFS, such a number can be computed by equating the BFS strain energy defect [righthand side of Eq. (3)] describing the X atom in the Y lattice with the ideal situation where X is located in an equilibrium X crystal: (11) n rpX e−aXrY =NrpX e−aXrX X eff Y where r (r ) is the equilibrium nearest-neighbor X Y distance in a crystal X (Y ), and p , a ( p , a ) are X X Y Y the ECT parameters of element X (Y ). For the case in hand, N=12 and n =10.69. This means eff that (beyond the fact that this is not an integer) a Pd atom needs 10.69 nearest neighbors in a Cu lattice to have zero BFS strain energy. This concept also helps to explain the attachment of the ejected Cu atom to the Pd(S) atom: by remaining in this site, the Pd atom increases its strain effective coordination, approaching the equilibrium value (n =12) as it increases the number of Cu nearest eff neighbors from four to nine (assuming that the ejected Cu atom is also a nearest neighbor). When attaching two N=1 states as shown in Fig. 7i and k, it is reasonable to expect that the growth pattern will be consistent with an increase of coordination leading to the ideal value noted above. This concept, while simple and intuitive, does not account for chemical interactions between atoms of dis-

381

similar species. Therefore, while it might help explain some noticeable effects (i.e. the choice of Pd for surface sites against overlayer sites), it only provides a partial explanation for more subtle effects. For example, the difference in energy between the configurations shown in Fig. 7i and k is not consistent with the difference in strain coordination (for the Pd atoms): while the Pd atoms in Fig. 7i have nine and 10 nearest neighbors respectively, those in Fig. 7k have 10 nearest neighbors each. This would indicate a preference for the configuration shown in Fig. 7k. It is therefore clear that a more detailed description of the chemical interactions is needed in order to fully account for the observed behavior based solely on the strain coordination concept. We conclude our analysis of the N=2 case by noting that the spectrum shown in Fig. 8 still indicates that interdiffusion of Pd in the Cu slab is energetically unfavorable [i.e. Pd(1b) and Pd(2b) states are still high in energy and therefore less likely to appear] related to experimental fact (f2), in that Pd atoms are highly stable in the Cu substrate. The case N=3 is also interesting as it shows the first hint of additional features observed experimentally. A small sample of some of the many accessible configurations is shown in Fig. 9 and the corresponding spectrum in Fig. 10. Once again, we note that the lowest energy state (to which all others are referenced ) is formed by the combination of Pd atoms in surface sites forming a (still incomplete) c(2×2)-like structure, with a Cu island firmly attached to such a structure [experimental facts (f3) and (f4)]. An additional feature is observed after comparing the spectrum shown in Fig. 10 with those for lower Pd coverage based on the fact that we can establish a correspondence between the size of the energy gap for the Pd(1b) state and the ground state with the probability of a Pd atom penetrating the surface layers. The states including interdiffusion in the Cu slab are now closer to the lowest energy state than for lower coverages, indicating an increasing likelihood for Pd penetration in the Cu lattice. Although this gap is still somewhat large, the fact that it is much smaller than it is for the N=1 and N=2 cases can be taken as a strong indication that Pd interdiffusion is not only to be expected but also

382

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389

Fig. 19. Configurations with eight Pd atoms corresponding to 0.133 ML Pd coverage. The energies of formation per adatom (in eV/atom) are also shown. (i–l ) Top view of four configurations with eight Pd atoms corresponding to 0.133 ML Pd coverage and their corresponding side views.

to become more prevalent for increasing Pd coverage. We conclude the discussion for N=3 by noting once again that the ground-state configuration can be seen as the coupling of three N=1 states in a specific orientation, due to an improvement of the strain effective coordination of the Pd(S) atoms.

The analysis of the N=4 case is a natural continuation of the features already identified for a lower coverage. Fig. 11 shows a small sample of some of the relevant configurations, and Fig. 12 shows the corresponding energy spectrum. All the features noted at a lower coverage are found in this case. In particular, it is seen that the lowest

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389

383

Fig. 19. (continued ).

energy state corresponds to a c(2×2) structure with a Cu island firmly attached to it [experimental facts (f3) and (f4)]. The trend for interdiffusion hinted at lower coverages becomes more noticeable as the states characterized by Pd(1b) and Pd(2b) atoms are in direct competition with alternative shapes of the Pd(S) atoms, and the ensuing Cu island [Pd(1b) and Pd(2b) states are clustered just above the ground state]. Moreover, the gap between the Pd(1b) and Pd(2b) states and the ground state is significantly smaller than that for lower values of N, indicating the possibility that interdiffusion of Pd in the Cu slab increases with coverage [experimental fact (f6)]. For N=5 and higher, no new features are observed. If any, closure of the c(2×2) arrangement at N=4 helps us focus our attention once again on the way additional Pd atoms join the pre-existing structure. Fig. 13 shows the most rele-

vant N=5 configurations and Fig. 14 the corresponding energy spectrum. Of particular interest are two slightly different N=5 configurations, where the difference resides in the way the fifth Pd(S ) atom joins the N=4 ground state (Fig. 13h and i). This specific choice of surface sites for the extra Pd atom indicates that a clear direction for growth is seen [experimental fact (f1)], which maximizes the strain effective coordination of the intervening atoms. The symmetry broken by adding a fifth atom to the c(2×2) arrangement results in a slightly higher energy for the whole configuration (−0.11897 eV/atom for N=5 and −0.17375 eV/atom for N=4). However, the narrowing of the gap observed between the ground state and Pd(1b) and Pd(2b) configurations continues, confirming the relationship between interdiffusion and coverage [experimental facts (f6) and (f8)].

384

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389

Fig. 20. Energy spectrum of the N=8 configurations shown in Fig. 19. Only those configurations with Pd atoms below the surface layer are explicitly labeled. The energies DE (in eV/atom) are relative to the lowest energy state.

Figs. 15–18 highlight the main features of the N=6 (Figs. 15 and 16) and N=7 (Figs. 17 and 18) cases, indicating the most relevant part of the energy spectrum and the lowest energy configurations. In the case N=8 (corresponding to 0.13 ML

Pd coverage), symmetry is restored, and the ground state shows features similar to those observed for N=4. Several significant configurations with N=8 and the corresponding energy spectrum are shown in Figs. 19 and 20. However,

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389

385

Fig. 21. Cu island on a Cu(100) surface (a) alloyed with a c(2×2) Pd structure in the surface layer, (b) somewhere between the c(2×2) Pd structure and a terrace and (c) adjoining a terrace. Cu atoms in the overlayer are indicated by dark gray disks, Pd atoms with black disks. The energies of (a) and (b) are referred to the lowest energy configuration (c).

by now, the gap between Pd(S) and Pd(1b) states is nearly washed out, indicating a higher likelihood for Pd interdiffusion. It should be noted that, while favored, it is still reduced to individual atoms only: states with more than 2 Pd atoms in layers below the surface are less energetically favorable, hinting at the possibility that new ordering patterns might be found at higher coverages, where Pd atoms in subsurface layers participate. The understanding of the process of Pd interdiffusion is not easily achieved by means of models or concepts that isolate particular effects in an attempt to identify the driving force for a particular mechanism. For example, the phenomenon of surface alloying of immiscible metals can be satisfactorily explained by means of atomic mismatch effects [31]. One single concept suffices to understand the behavior of the atoms in the surface: the

competition between the elastic energy of the adatoms and those same atoms in the bulk, even in the absence of any surface energy effects, provides a simple and clear explanation of the observed phenomenon. In the case of Pd interdiffusion, however, simple and intuitive atomic size-effect concepts are not necessarily sufficient, and a general description, which involves (BFS ) strain energy and (BFS ) chemical energy effects as well as the interplay between them, is necessary to reproduce the observed behavior. One final feature can be recognized, in spite of the fact that we do not account for the motion of Cu(O) atoms along the Cu surface. In the presence of terraces and steps, three distinct regimes can be modeled where the Cu islands are either attached to the step, nucleated on top of the underlying Pd c(2×2) phase, or randomly located at any other

386

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389

Fig. 22. Summary of the energy spectra N Pd atoms (N=1,…,8). The lowest energy state is indicated by a thick line, corresponding to the configurations shown in previous figures. The thick dashed line indicates the first configuration that includes a Pd atom below the surface layer [Pd(1b) or Pd(2b)]. The dotted lines indicate the first configuration that includes Pd atoms in the overlayer [the fact that such states do not appear for N=1, 6 and 8 is because configurations with Pd(O) were not included in the catalogue of selected distributions (N=6, 8) or because they exceed the highest energy shown (0.60 eV/atom)]. The lowest energy configuration for each value of N is also shown.

point on the plain Cu surface. The calculations show that attaching the ejected Cu atoms to steps is always energetically favorable. However, such a preference is less significant for increasing coverage as the difference in energy between these states and those where the Cu island nucleates on Pd or form somewhere else on the surface is less marked. Fig. 21 shows some significant results for N=4. This lowering of the ‘energy barrier’ between these two distinct behaviors is consistent with the experimental finding of a decrease in mobility of Cu islands with increasing coverage [5]. A similar situation is found for other values of N: considering the same three island locations depicted in Fig. 21, it is found that there is a trend toward a lowering of the ‘energy barrier’ (i.e. the energy

of the configuration where the island is located somewhere between the step and the Pd atoms), as well as a less marked difference between the energy of the other two configurations, where the island is attached either to the step or to the Pd cluster. For example, the gap between these last two states decreases from 0.37 eV/atom for N=1 to 0.20 eV/atom for N=4, while the difference in energy between the energy barrier and the configuration where the Cu island is attached to the Cu step decreases from 0.49 eV/atom for N=1 to 0.32 eV/atom for N=4. Fig. 22 summarizes the results found so far for N=1 through N=8, by showing the relevant portions of the different energy spectra shown in Figs. 5, 8, 10, 12, 14, 16, 18 and 20, indicating the

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389

closing gap between the c(2×2) configurations and those with Pd(1b) and Pd(2b). We now increase the coverage to 0.5 ML Pd. It can be easily shown that no new features appear for intermediate coverages, beyond those already found for a low coverage. The results shown in Fig. 23 prove that up to this coverage, the c(2×2) structure is still energetically favored against other ordered structures ( Fig. 23a and b) or island shape (Fig. 23a and c). While the number of different

Fig. 23. Configurations with 30 Pd atoms corresponding to 0.5 ML Pd coverage: a 30-atom Cu island (a) alloyed with a c(2×2) Pd structure in the surface layer, (b) alloyed with a p(2×2) Pd structure in the surface layer and (c) same as (a), but with the 30-atom Cu island broken into two disjointed 15-atom Cu islands. The notation labeling each figure denotes the number and type of atoms on the left- and right-hand side of the cell (top line) and the constitution of the surface layer. Also, the energy of formation per adatom (in eV/atom) is included.

387

accessible configurations is now much greater than that for lower coverages, the 0.5 ML case can be summarized by considering three basic configurations. Most of the remaining possibilities can be easily ruled out as they would entail growth patterns for which there is no evidence at low coverage (i.e. a cluster of Pd atoms in the surface layer, etc.). The lowest energy configuration for this coverage is, as expected, the one where all the Pd atoms form a c(2×2) structure and a Cu island forms in the overlayer [experimental fact (f7)]. However, as noted earlier, increasing Pd coverage leads to the possibility of increased interdiffusion of Pd atoms to subsurface layers. Within the restrictions imposed in our calculation (mainly the lack of individual atomic relaxations), a somewhat surprising result is found when considering the energetics of one Pd atom located in subsurface layers for an otherwise ‘perfect’ configuration such as that shown in Fig. 23a. In Fig. 24, two particular cases are shown: one of the Pd atoms in Fig. 23a is now allowed to occupy a site in the first layer below the surface (Fig. 24a) or in the second subsurface layer (Fig. 24b). Surprisingly, this second configuration is slightly lower in energy, reversing the trend observed for lower coverages where Pd(1b) states were always favored over Pd(2b) states. However, the energy difference between these two configurations is exceedingly small, leading us to consider the possibility that slight relaxations of the Pd(1b) atom in Fig. 24a might suffice to make Pd(1b) states once again energetically favorable against Pd(2b) states. This latest result indicates that, at this level of coverage, a proper treatment of individual atomic relaxations, not included in the current analysis, is necessary in order to provide a complete description of the ensuing structures. However, within the framework established for this study, it is also clear that in spite of the simplicity of the methodology and the straightforward analysis of the results, a great deal of information can be obtained, allowing us to identify every single feature of relevance during the surface alloy formation process and, as is the case for 0.5 ML, also give us an indication of what restrictions should be lifted in order to provide a complete description of the system.

388

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389

Fig. 24. Side view of configurations with 30 Pd atoms corresponding to 0.5 ML Pd coverage. A 30-atom Cu island is located in the overlayer (top row of gray disks), while 29 Pd atoms are located in the surface plane in the c(2×2) structure denoted by alternating Pd (black disks) and Cu (shaded disks) atoms. One Pd atom is located in (a) the first plane below the surface and (b) two planes below the surface aligned with a Pd atom in the surface layer. The energy of formation per adatom (in eV/atom) is also noted.

4. Summary

References

Summarizing, all the essential features observed experimentally up to 0.5 ML Pd coverage [5] are properly described in this simple calculation: the alloying of Pd atoms into the surface in the [010] and [001] directions, the nucleation of Cu islands on top of alloyed areas resulting in subsurface Pd, the decreased mobility of Cu islands with increasing coverage, the formation of a c(2×2) phase as the chains converge, and the interplay between the c(2×2) phase and the initiation of second-layer growth at increasing coverages before the c(2×2) is fully completed. We also showed that a simple methodology (i.e. the generation of catalogues and the analysis of results via energy spectra), tied with a powerful but computationally simple energy method (the BFS method for alloys), provides a straightforward and detailed description of the surface alloy formation process.

[1] U. Bardi, Rep. Prog. Phys. 57 (1994) 939. [2] A. Christensen, A.V. Ruban, P. Stoltze, K.W. Jacobsen, H.L. Skriver, J.K. Norskov, Phys. Rev. B 56 (1997) 5822. [3] S.C. Wu, S.H. Lu, Z.Q. Wang, C.K.C. Lok, J. Quinn, Y.S. Li, D. Tian, F. Jona, P.M. Marcus, Phys. Rev. B 38 (1988) 5363. [4] P.W. Murray, I. Stengsgaard, E. Laegsgaard, F. Besenbacher, Phys. Rev. B 52 (1995) R14404. [5] P.W. Murray, I. Stensgaard, E. Laegsgaard, F. Besenbacher, Surf. Sci. 365 (1996) 591. [6 ] G.W. Graham, P.J. Schmitz, P.A. Thiel, Phys. Rev. B 41 (1990) 3353. [7] T.D. Pope, K. Griffiths, P.R. Norton, Surf. Sci. 306 (1994) 294. [8] T.D. Pope, M. Vos, H.T. Tang, K. Griffiths, I.V. Mitchell, P.R. Norton, W. Liu, Y.S. Li, A.R. Mitchell, Z.-J. Tian, J.E. Black, Surf. Sci. 337 (1995) 79. [9] G.W. Graham, Surf. Sci. 171 (1986) L432. [10] Y.G. Shen, J. Yao, D.J. O’Connor, B.V. King, R.J. MacDonald, J. Phys.: Condens. Matter 8 (1996) 4903. [11] T.D. Pope, K. Griffiths, V.P. Zhdanov, P.R. Norton, Phys. Rev. B 50 (1994) 18553. [12] J. Yao, Y.G. Shen, D.G. O’Connor, B.V. King, J. Vac. Sci. Technol. A 13 (1995) 1443. [13] R. Tetot, J. Kudrnovsky, A. Pasturel, V. Drchal, P. Weinberger, Phys. Rev. B 51 (1995) 17910. [14] M. Valden, J. Aaltonen, M. Pessa, M. Gleeson, C.J. Barnes, Chem. Phys. Lett. 228 (1994) 519. [15] T.D. Pope, M. Vos, H.T. Hang, K. Griffiths, I.V. Mitchell, P.R. Norton, W. Liu, Y.S. Li, K.A.R. Mitchell, Z.-J. Tian, J.E. Black, Surf. Sci. 337 (1995) 79. [16 ] J. Yao, Y.G. Shen, D.J. O’Connor, B.V. King, Surf. Sci. 359 (1996) 65. [17] G.W. Anderson, T.D. Pope, K.O. Jensen, K. Griffiths, P.R. Norton, P.J. Schultz, Phys. Rev. B 48 (1993) 15283. [18] T.D. Pope, G.W. Anderson, K. Griffiths, P.R. Norton, G.W. Graham, Phys. Rev. B 44 (1991) 11518.

Acknowledgements Fruitful discussions with N. Bozzolo are gratefully acknowledged. G.B. would like to thank the Centro Ato´mico Bariloche, CNEA (Argentina) for their hospitality. J.G. acknowledges the support of the NASA/OAI Collaborative Fellowship Program. This work was partially funded by the HITEMP, PPM and HOTPC programs at NASA Glenn Research Center.

J.E. Garce´s et al. / Surface Science 459 (2000) 365–389 [19] J. Kudrnovsky, S.K. Bose, V. Drchal, Phys. Rev. Lett. 69 (1992) 308. [20] G. Bozzolo, J. Ferrante, J. Computer-Aided Mater. Design 2 (1995) 113. [21] G. Bozzolo, R.D. Noebe, J. Ferrante, C. Amador, J. Comput.-Aided Mater. Design 6 (1999) 1. [22] G. Bozzolo, J. Ferrante, R.D. Noebe, C. Amador, Scripta Mater. 33 (1995) 1907. [23] G. Bozzolo, J. Ferrante, R.D. Noebe, C. Amador, Scripta Mater. 36 (1997) 813. [24] G. Bozzolo, R. Ibanez-Meier, J. Ferrante, Phys. Rev. B 51 (1995) 7207.

389

[25] G. Bozzolo, R.D. Noebe, J. Ferrante, A. Garg, F. Honecy, C. Amador, J. Comput.-Aided Mater. Design 6 (1999) 33. [26 ] G. Bozzolo, J. Ferrante, A.M. Rodriguez, J. Comput.Aided Mater. Design 1 (1993) 285. [27] J.R. Smith, T. Perry, A. Banerjea, J. Ferrante, G. Bozzolo, Phys. Rev. B 44 (1991) 6444. [28] J.H. Rose, J.R. Smith, J. Ferrante, Phys. Rev. B 28 (1983) 1835. [29] A.R. Koymen, K.H. Lee, G. Yang, K.O. Jensen, A.H. Weiss, Phys. Rev. B 48 (1993) 2020. [30] J.B. Hannon, H. Ibach, P. Stoltze, Surf. Sci. 355 (1996) 63. [31] J. Tersoff, Phys. Rev. Lett. 74 (1995) 434.