. 

cr2 SIMULATION STUDIES FOR SURFACES 

AND MATERIALS STRENGTH 

Semi-Annual Report 

for 

Cooperative Agreement NCC2-297 

for the period 

November 1, 1986 - April 30, 1987 

Submitted to 

National Aeronautics and Space Administration 

Ames Research Center 

Moffett Field, California 94305 

Computational Chemistry Branch 

Dr. David Cooper, Chief and Technical Monitor 

Thermosciences Division 

Dr. Jim Arnold, Chief 

Prepared by 

ELORET INSTITUTE 

3788 Fabian Way 

Phone: 408 730-8422 and 415 493-4710 

Palo Alto, CA 94303 

K. Heinemann, President and Grant Administrator 

Timur Halicioglu, Principal Investigator 

(tiASA-Ca-180539) ISIEULATICI SIUCIES FOR N87-22277 

SUBPACES AND HAIERIALS STEENG1IP Semiannual 

iieport, 1 Hov. 1466 - 30 Apr- 1987 (Bloret 

Corp.) 13 p Avail: "XIS EIC AOz/HF A01 Unclas 

CSCL 20K G3/39 007U105 I 
Three different investigations based on computer simulations were carried out 

during this period. Simulation calculations were performed using model potentials 

with two- and three-body interactions which were represented by the Mie and the 

Axilrod-Teller potentials, respectively. In the first part, to analyze stability criteria 

for simple crystals, a parametrical investigation was carried out and stability regions 

for five different crystalline structures were determined with respect to potential 

energy parameters. To simulate some of the basic surface processes such as diffusion 

and nucleation, in the second and third parts, we considered silicon surfaces as a 

specific example. In these studies diffusion of adatoms on stepped surfaces were 

simulated and energetics associated with the formation of kink sites were calculated. 

1. Stability Diagrams for Simple Crystals 

In simulation calculations structural stability of the model system is an important aspect. Often, it is desired that the structure of the system under consideration 

corresponds to the energetically most favorable configuration. For a given potential 

energy function, finding the most stable structure among different possible configurations is a cumbersome and time consuming procedure. This problem emerges, in 

particular, for potential energy functions involving three-body interactions which 

furnish stability regions for a large number of crystalline structures ranging from 

the close-packed (HCP or FCC) to more open (e.g., diamond cubic or p-tin) structures. The significance of three-body interactions in calculating the energy and 

structure-related properties of simple crystals were shown in our earlier investigat ions. 

In the present study, stability regions for five different crystalline structures 

(HCP, FCC, BCC, @-tin and diamond cubic) were calculated for a wide range of 

potential parameters. Calculations were carried out in a dimensionless form using 

the Mie potential and the Axilrod Teller function to represent twe and three-body 

interactions, respectively. The results are presented in a graphical form in Figure 

1. This representation is found to be very useful in deciding parameters for various 

specific systems in which two-body and three-body interactions are operational. 

The stability map shown in Figure 1 represents clearly the effect of the three-body 

interaction on the stability of crystals. A potential energy function based only 

on a two-body potential is able to provide stable regions only for HCP and FCC 

structures. In analyzing Figure 1 one should remember that every point on these 

maps corresponds to the energetically most stable structures among the five different crystals taken into consideration in this work. The general tendency in the 

stability of these crystals as a function of 2’ is changing as HCP-tFCCdBCC 

for small n, as HCP-FCC-t P-tin-tdiamond cubic for intermediate values of n 

1 
and as HCP-tFCC +diamond cubic for larger n values. Because of the structural 

similarities, the difference between energies calculated for HCP and FCC structures 

turned out to be very small. Within our accuracy limits the HCP and FCC structures were considered energetically degenerate for regions denoted by HF which 

occupies large areas in the stability maps for n 5 6 (see Figure 1). 

Boundaries between the stability regions should not be taken as sharp limits. 

Due to the discrete nature of the calculations all boundaries must be regarded 

rather as somewhat broader transition regions. The map for n = 5 is interesting; 

for example for Z* = 0.825 and rn = 9 all five crystalline structures seem to have a 

common region. The @-tin structure has a narrow stability region from n = 5 to 

n = 8. 

2. Adatom Diffusion and Ledge Interactions on the Si(ll1) Surface 

Knowledge concerning the diffusion of atoms on surfaces containing ledges is 

essential for the proper understanding of crystal growth phenomena at an atomic 

level. Early studies on metal surfaces revealed the importance of adatom-ledge interactions for surface diffusion when attempting to evaluate different growth mechanisms. Based on calculations of surface structure and energetics for single atoms, 

the motion of particles across an otherwise flat ideal Si(ll1) surface will occur by 

activated jumps between hole sites. These jumps will carry the atom through the 

cradle sites so that the minimum energy path for the motion of an atom from the 

hole to the cradle site will control this diffusion process. Figure 2 defines these and 

other possible adatom positions. Of course, the actual macroscopic process of surface diffusion depends on the detailed structure of the surface so that diffusion via 

this simple unit process is expected to be somewhat modified by (7 x 7), (2 x 1) or 

(1 x 1) surface reconstruction patterns and by the presence of ledges on the surface. 

In the present study, the ”Constrained Statics” method is utilized to find the 

minimum energy configuration of the surface when one adatom is constrained to 

relax only in a specific plane while all other atoms relax in all three dimensions. 

Many similar computations are performed for the specific atom in different initial 

positions and minimum energy trajectory which takes the atom from one position 

to another is thus mapped out. The result for the hole to cradle trajectory is shown 

in Figure 3. The moving atom stayed in the vertical plane containing the two sites, 

as expected, so that it passed through the bridge site with an interaction energy of 

1.24 eV (see Figure 2), which is clearly an inflection point on the energy trajectory. 

A small additional energy hump (0.12 eV) exists between the bridge site and the 

metastable cradle site (0.99 eV). Upon analysis of overall structure when the adatom 

was near the energy maximum (1.36 eV), the energy hump was seen to be caused 

by the development of a significant local reconstruction underneath the cradle site. 

2 
The subsequent cradle to hole site energy barrier, to complete the overall hole to 

hole motion, is E, = 0.36 eV. 

These are the absolute minimum energy changes that can be associated with 

the jump. In practice, the large local atomic reconstruction at the cradle site is not 

likely to happen during an adatom jump so that the hole site to hole site adatom 

movement could exhibit an average maximum activation energy somewhat greater 

than 1.36 eV. Experimental estimates for the Si( 111) surface diffusion activation 

barrier range from 0.2 to 2.0 eV due to the very indirect nature of the data which 

must be analyzed and interpreted, and because of the uncertainty of the local surface 

structure and composition. 

Adatom-ledge Interactions 

Although it is important to study all adatom-ledge interactions, only the hole 

adatom-[211] ledge interactions are considered here. Figure 4 gives side views of the 

12111 lower and upper ledges, respectively, displaying the adatom positions tested. 

Simulation calculations were performed with ledges at 25.45A spacings. Figure 4c 

plots the hole formation energies, Ef, using circles for the lower ledge and squares 

for the upper ledge. The data points are vertically coincident with their spatial 

positions in Figures 4a and 4b. The horizontal dashed line shows the formation 

energy for the adatom at a hole site on the flat ideal surface, labeled as Ef, indicating 

the energy that would obtain at an infinite distance from the ledge. 

From our earlier study, we found that the hole adatom puts the surface further 

into compression, and that a 0.5% in-plane tensile strain caused the Ef of the hole 

adatom to decrease (become more favorable) by 0.08 eV. In Figure 4, where the 

upper ledge causes a major reduction of the surface compression, the hole Ef should 

be much lower near the edge, as Figure 4c clearly confirms. In fact, the stress effect 

is larger on the terrace above the ledge than on the lower terrace, and the Ef 

is also more negative above the ledge. The increase in the adatom energies right 

next to the ledge is arguably due to the interference of short range reconstruction 

displacements of the adatoms and ledge which begin to dominate the longer range 

stress interactions. 

The lower [211] ledge has a much smaller effect on the surface compression, but 

it still decreases it somewhat, particularly on the terrace above the ledge. However, 

in the explicit calculation of Figure 4, the hole adatom is not at all favored on the 

upper site of this ledge, and it is repulsed significantly from the terrace below the 

ledge. It may be that the shear fields may offset the otherwise favorable effect of the 

normal stresses above the ledge. The repulsion below the ledge does not have any 

obvious explanation in terms of simple stress interactions. As argued above, the 

local reconstructions interfering with one another may be the source of repulsion. 

Or, the stresses in the underlying puckered layers, which have not been extensively 

analyzed thus far, may be involved in the repulsion effect. 

3 
The only firm conclusion to be drawn from these results is that the surface 

behavior is not always so easily cast into a simple top layer stress interaction model. 

In the case of the upper ledge, where the stress effect is huge, the stress analysis 

is sufficient. For the lower ledge, where the surface compression is not a dominant 

feature, there is no simple correlation between the stress effects of the separate 

defects and their interaction energy. 

Perhaps the most important point to notice about Figure 4 is that the adsorption energies of adatoms near the ledge are symmetrical, an assumption that has 

been used in most computations of adatom diffusion to ledges via both the lower 

and upper terraces. Of course, these numbers are expected to change when more 

adatoms and surface vacancies are present and when the terraces reconstruct into 

distinctive patterns. 

3. Kink-Site formation Energies on the Si(ll1) Surface 

A semiempirical potential energy function (PEF) has been developed for Si 

comprising both two-body and three-body effective potentials that has significant 

application to a variaty of phenomena ranging from cluster stabilities to crystal 

growth considerations. The parameters contained in this PEF were obtained from 

fits to various experimental properties. Using this PEF, as described in our earlier 

reports, more than 25 different materials properties of Si bulk, clusters and surfaces 

were successfully reproduced. The calculated energy quantities are all within *20% 

while structural properties, such as bond distances for small clusters and the nearest 

neighbor distance in relaxed crystals, were reproduced to within 7%. However, 

quantities related to the second derivatives of energy, such as bulk modulus, elastic 

constants and vibrational properties were not predicted so well and only qualitative 

agreement is obtained for such parameters. On the other hand, surface energy 

values for the Si(lOO), (110) and (111) surfaces were predicted with good precision 

and the PEF produces the correct surface reconstruction geometry for the (2 x 1) 

Si(100) surface and for the (2 x 1) Si(ll1) surface. 

In the crystal growth, as well as in the surface reconstruction processes, energy 

values associated with the kink formation are very important. It would be very 

nice if one could make first principles ab initio calculations of surface quantities like 

ledges and kinks; however, such calculations are only possible for very small systems 

containing only a few atoms. Since long-range surface stress effects due to interacting ledges on Si(ll1) surfaces have been calculated and found to be important to 

make meaningful ledge or kink calculations with any PEF, about 500-1000 atoms 

must be taken into account. It is for this reason that one is presently limited to the 

use of semiempirical potential functions. Of the alternate semiempirical Si PEF’s 

available today none has the track record of that proposed here. Each undoubtedly 

has its strengths and weaknesses but insufficient experience is available to catalogue 

4 
their applicability in a given situation. Perhaps it is best to look at each as the 

analogue of a particular experimental technique that provides a measurement with 

a certain error bar. Refinement of the technique generally reduces the error bar. 

Such semiempirical PEF's are able to make a determination of ledge and kink energies which is more than present experimental techniques can do. Granted, one 

does not really know the magnitude of the error bar generated by using a particular 

PEF so the number should be treated with caution until substantiated in some alternate fashion; however, it is the best path presently available to determine some 

quantities important to crystal growth modeling. 

Using the present PEF, most recently, the structure and energetics for the 

dominant ledges on the Si(ll1) surface have been calculated. These ledges are free of 

kinks and, therefore, are not the equilibrium configuration except at 0"K. However, 

the properties of four different kind of [2ii] and [211] ledges that can be found on 

the Si(ll1) provide a necessary baseline and are an important structural basis for 

kink formation studies. These four structures are presented in Figure 5. The ledge 

energy per atom for these four structures at large ledge spacing are X = 0.24 and 

0.01 eV for [21T] lower and upper ledges, respectively, while X = 0.192 and 0.01 eV 

for (2111 lower and upper ledges, respectively. This report presents calculations of 

the kink structure and the kink energetics on these ledges. These kink properties 

are also expected to depend upon the state of surface reconstruction, either by a 

direct structural connection or indirectly through the surface stresses. Although the 

surface reconstruction problem should be attacked first and a more thorough study 

of ledges in contact with such reconstructed surfaces is needed before the nature of 

the kink sites can be studied in detail, there is pedagogical value to be gained from 

a study of kink formation on these 0°K ledges. Because of these limitations, only 

a preliminary study of a few kink sites will be presented here. 

Results for isolated kink formation on both the (ZiI] and [Zll] upper ledges 

and the [211] lower ledge at 0°K are summarized in Figure 6. Although three 

examples do not constitute a statistically significant trend, the main feature of 

interest is that the kink formation energy is low when the reconstructed ledge energy 

is high. The ledge energy is low when the reconstruction produces a large tensile 

stress component and reduces the net local stress to very small values. Thus, little 

additional stress relaxation can occur via kink formation. For those ledges that do 

not reconstruct strongly, kink formation has the opportunity of strongly relaxing 

the local compressive stress. 

Large values of Ek will lead to a small kink density and thus slow ledge movement during crystal growth at a given driving force. Small values of Ek lead to high 

kink densities and thus rapid ledge movement during growth at such driving force. 

Thus, although the [Zii] and [Zll] ledges will bound the terraces on the Si(ll1) 

surface at equilibrium for energetic reasons, they will be the slowest growing ledges 

and will thus also bound the terrgces during crystal growth for kinetic reasons. In 

5 
actual practice, the foregoing conclusion is always valid, however, for doped materials interaction energies must be properly modified to obtain correct kink energy 

values. 

6 
Figure Captions 

1) Stability regions for five simple crystalline structures, HCP (H), FCC (F), BCC 

(B), ,&tin (S) and Diamond cubic (D). Parameters m and n are the repulsive and 

attractive exponents of the Mie potential, respectively, while 2' represents the 

three- body intensity parameter. 

2) Top view of postulated adatom adsorption sites on the Si(ll1) surface 

3) Minimum energy path for the motion of an adatom from the hole site to the 

cradle site. Energies are in eV measured relative to the hole site adatom formation 

energy. 

4) Formation energies for hole adatoms near the (a) lower and (b) upper 121 l] ledges. 

Adatom positions are shown with filled triangles in (a and b). Energies are plotted 

directly below in (c), using circles for the lower ledge and squares for the upper 

ledge. 

5) Relaxed structures of the basic (2111 and [Zii] single height ledge on Si(ll1). 

6) Structures and energies for selected single height ledges (Aoo) and kink sites (Ek) 

on Si(ll1). 

7 
P L. 

X n 

0. P L. 

I 

-. * 

L. 

I 

bn . . . . . . , . , . N 

d P 
I 

A Hole 

0 Cradle 

-1 Bridge 

0 Point 

FIGURE 2 

FIGURE 3 
A 

B 

0 

0 o0 o0 o0 00 

0 

0 

0 0 

0 0 0 0 

0 

0 

0 0 

0 

0 .. 0 0 0 

0 0 0 0 

0 0 0 

0 00 

0 e 0 

C 

FIGURE 4 
- 

SIDE VIEW TOP MtW ' (tii) 

LEDGE 1 

II OI 

~00000 

...... * 000000 

6..... LOWER 

000 0 0 0 .O .O 

00000 

0.00. 

FIGURE 5 
i Ledge 1 ! 

I i i 

! (zii) I 

UPPER 

(211) 

UPPER 

I 

0.01 

0.01 

I (zii) I i 

LOWER 1 0.24 i 

- 

Top View of Structure Ex I 

I 

0.75 

Y~OOOC) 

...... 

000 

e..... a a o 0 o t. 0 

0.73 

FIGURE 6