- Email: [email protected]
MBE growth mechanisms; studies by Monte-Carlo simulation
ELSEVIER
Thin Solid Films 267 ( 1995) 3746
MBE growth mechanisms; studies by Monte-Carlo simulation
’
H. Sitter Johannes
Kepler Universitiir Linz, Institute oj’Experimenta1 Physics, A-4040 Linz, Ausrria
-. Abstract A hlonte-Carlo-based method for simulation of epitaxial crystal growth is reported. The improvement of this concept in comparison with based on the solid-on-solid model is the possibility of simulating growth without the assumption of a discrete crystal lattice structure. The used continuous space model requires only the assessment of the Hamiltonian. In the first step the model was applied for studying growth in a two-dimensional cross-section through the substrate and the epilayer with a width of up to 80 atoms. Two different potentials can be chosen for the simulation: (a) Lennard-Jones and (b) directional Lennard-Jones to describe covalent systems. In addition to the selection of the potential, we can set the following parameters for different kinds of atoms: radii of the potential minimum, binding energy between the atoms for all possible pairings, growth temperature and substrate structure. Island and layer growth in the nucleation stage were investigated by selecting different binding energies of the epilayer atoms to the substrate atoms. The epitaxial growth of strained layers was simulated to study the incorporation of dislocations at the interface. The continuous space model also allowed calculation of the rhombic distortion of a strained layer grown on a tilted substrate with atomic steps. As a further extension of the two-dimensional model the re-evaporation of particles was allowed which is essential for the simulation of atomic-layer epitaxy. The obtained results for the growth rate in atomic layer epitaxy as a function of substrate temperature can be compared with experimental findings. In .I second attempt we extended the continuous space model to three dimensions based on the same assumption as tested in the twodimensional model. The three-dimensional model allowed the same simulations for islands and layer growth but also the selection of vicinal surfaces of the substrates with respect to the usually used ( 111 )-oriented surfaces. The simulation showed a clear stabilisation of the growth direct,on by the atomic steps on the substrate and even in lattice -matched structures the formation of twins was suppressed, the simulations
Keywc~rd~: Adsorption; Deposition process; Growth mechanism; Nucleation
1. Introduction Acvances
in crystal
growth
technologies,
in particular
in
beam epitaxy (MBE) [ 1 ] and atomic layer epitaxy (ALE) [ 21 have stimulated the search for realistic Irheoretical models of crystal growth. So far, phenomenological theories using thermodynamics were used to describe processes such as condensation and re-evaporation or misfit relaxation and incorporation of dislocations. With the ard of modern computer technology it became possible to realixe ab-initio calculations with realistic interaction potentials. These microscopic models increase the understanding of the physics of crystal growth and epitaxy. In principle, there are two possible simulation methods. First. one can solve the equation of motion of the molecules to track the system. This method is calledmoleculardynamics (MD). Second, it is possible to calculate the properties of a system at equilibrium from its distribution function in the the field of molecular
’ Presented as an invited talk at the Workshop on MBE-Growth and Technology, Warsaw, Poland, December 1994.
Physics
0040-6090/95/$09.50 0 1995 Elsevier Science S.A. All rights reserved .SSDIOO40-6000(95)06593-8
phase space. Because this is done with the use of random sampling, it is called the Monte-Carlo method (MC). The name Monte-Carlo method is used for all models requiring random numbers, but in this work we restrict ourselves to algorithms similar to the so-called important sampling method as described by Metropolis et al. [ 3 I. For both, MD and MC, one only needs to know or assume the Hamiltonian of the system. The main difference between MD and MC is the fact that MD can describe dynamic systems, while by MC no real time behaviour can be calculated because the algorithm relies on thermodynamic equilibrium. However, the MC calculations are performed in such a way that the system is disturbed by insertion of a new particle, e.g. an adatom which can be incorporated into the growing crystal, and then equilibrated until the change in the potential energy by more MC steps is negligible. This procedure can be repeated several times, and in this way dynamic processes with an artificial time scale can also be simulated by MC methods. On the other hand, the solution of the coupled equations of motion for any particle of the system in MD restricts the number of particles
38
H. Sitter / Thin Solid Films 267 (I 995) 3746
and also the range of real time simulation because of limitations in CPU time. Especially if the interaction potential is a complicated function, the MC method is favourable, because only one scalar number, the change of the potential energy, has to be calculated. In a huge number of computer simulations the growth of lattice-matched crystals is described by an Ising model [47]. All these models use a predefined lattice structure. Although these models are quite simple, they have been used successfully to study the surface roughening with increasing temperature [ 51, the transition between island and layer growth [ 61 and propagation of spiral dislocations [ 71. More sophisticated simulations based on Ising models allow the simulation of lattice-mismatched structures [ 8,9]. However, calculations using continuous space simulations are quite rare [ 10-17 1. Due to the complexity of continuous space models, their application has to be restricted to very small three-dimensional crystals or two-dimensional crosssections. Continuous space models using MD and MC methods provide some advantages in comparison with discrete space models, because the assumption of the Hamiltonian is sufficient to describe the physical behaviour of the system and the models are somewhat closer to reality. The most often used potential for the simulation of crystal growth is of the Lennard-Jones (LJ) type, which has been used for MC [ 10,15-171 and MD [ 111 simulations. For real semiconductor structures, a covalent bonding mechanism should be taken into account. This was achieved by a directional LJ (DLJ) potential, where the particles are defined by their locations and their orientation [ 12,15-171. By defining for each atom three bonds in an angle of 120” [ 121 or by using four bonds in an two-dimensional projection [ 161, this model was used to predict the stability range in strained layer superlattices. Another way to describe covalent bonds is to set up an interaction potential consisting of a LJ term together with a three-particle part, which favours a tetrahedral atomic configuration. With this model the epitaxial growth of silicon on silicon was simulated by MD calculations [ 131. At present, there is no general rule for the choice of the simulation model. All the models mentioned above have advantages and disadvantages. With the model presented in this paper it is principally possible to tackle the problem of lattice mismatch and miscut, substrate orientation, substratetemperature dependence, surface morphology, incorporation of dislocations and the flip from island to layer growth for heteroepitaxy, including the growth of binary and ternary compounds on binary substrates. In Section 2 we will discuss the Hamiltonian and certain restrictions for the simulation due to limitations in CPU time. In Section 3 we present some of the already achieved results. The most important results of our calculations are: the incorporation of misfit dislocations in a GaAs/CdTe heterostructure, and the change from island to layer growth as a function of binding energies between substrate and epilayer atoms in the two-and three-dimensional model. The prediction of the stability range in an ALE
process could be calculated in the two-dimensional and is comparable with experimental results.
2. Description
model
of the model
2. I. Two-dimensional
model
We use a continuous-space MC method to simulate epitaxial crystal growth in a two-dimensional cross-section through the substrate and the epilayer parallel to a (001) plane of L zinc-blende structure. In this cross-section the four tetrahedra: bonds of each atom are oriented in the crystal lattice as shown in Fig. 1 (a). It is evident that the projection of the tetragonal zinc-blende lattice into the plane gives a four-fold symmetry. This restriction to two dimensions was applied to the simulation to save CPU time during the first calculations in order to test the model. In our model each particle is represented by two twodimensional vectors in the plane of the cross-section, namely its location r and its orientation e with (e 1= 1. The pairwise interaction potential is characterised by three parameters: the binding energy E,, the binding length r,,, and the parameter p which describes the covalent character of the bond. The potential energy of two interacting particles is given by Eqn. (1).
‘\ ~100) (001)
Surface
(4
atom 2 Fig. 1. (a) (001) surface of a zinc-blende structure with the usual crystallographic directions and the corresponding bonds. (b) Schematic drawing of two particles at a distance of rj2, with their orientations e,, eZ and the definition of the angles cp,. cpz,
H. Sifter/ Thin Solid Films 267 (1995) 3746
39
r= Irl
r=r,-r2
r.e, cos
((0”)
=7
f(cp,)=cos2 o= 2$,6
Fig. Pairwise binding energy El2 (r. 9, 0) in units of EO, plotted in Cartesim coordinates for three different values of p.
I
/co~“-(1-P)~-“-Pf((p,>f((p2)m-6 r
a2
r
CUlOff
0 r 2 rculoff
(1) The pairwise interaction energy is a function of the variables r, 2, ‘p, and (p2which are shown in Fig. 1 (b) . The distance r12, the angular dependencef, the parameter (Y,and the cutoff radius rcut_Off are taken as
(2cp,)withnE(1,2}
(2)
rCuloff= 2Jr,,
A graphical representation of the potential energy for three different values of the parameter /3 is given in Fig. 2. The interaction potential can be tuned by the parameter /3 from the pure LJ potential with p = 0 (Fig. 2(a)) to the pure directional LJ potential with /3 = 1 (Fig. 2(c) ) . Since different species of atoms are allowed, E,, req,,, and j3 are specific for the atom species involved in the interaction. Therefore, we denote these parameters by E,,,j, requ,;, and PO to describe an interaction of an atom of kind i with an atom of kindj. To simulate the growth of compound epilayers on compound substrates our model allows one to select I?,,,,, rqu.,, and p,, for n different kinds of atoms. That means, in principle, 3n( n + 1) /2 parameters describing the interaction of atoms can be set at will. To save CPU time we set rq “,,, = ( rCq“,,, + r,.,,) /2 and pi,= pii = p,j= p. From now on the symbols 0, X and * denote the three kinds of atoms according to the representation of atoms in the lattice picture. By definition, 0 indicates a substrate atom while X and * indicate an atom of the epilayer. To simulate the growth of a crystal we use a substrate width of 80 atoms and periodic-boundary conditions in the horizontal direction. We start by placing some layers of substrate atoms in the distance of rqu OO. The substrate is equilibrated by means of a conventional’important sampling MC method [ 31; this means that one particle is selected at random and moved to a position within a rectangle with a horizontal length of 2.5 rqu,OO and a height of roqu,OO,The orientation e of the particle is chosen at random also. This move and rotation is accepted with the probability p = exp ( - AHIKJ). This procedure, commonly called the MC step, is repeated until an acceptable approach to the equilibrium is achieved. After the substrate is equilibrated, another particle is placed near the surface at a random position and the system is equilibrated again. After some particles are deposited, expectation values of the position of each particle are calculated. The difference between the actual and the average position is used as an indication of the stability of the lattice. This procedure is repeated until the required number of particles is deposited.
2.2. Model for atomic-layer epitaxy
Ultra-high vacuum (UHV) ALE of a binary compound on an elemental substrate is characterised by the following conditions.
H. Sitter/Thin
Solid Films 267 (1995) 3746
Fig. 3. Two-dimensional lattice used for the simulations of UHV-ALE.
(i> UHV ALE is only possible for binary systems where the bonds between the different atomic species are stronger than the bonds between the same atomic species, which is true for many II-VI compounds. (ii) In contrast to MBE growth, in an UHV ALE growth process a shutter sequence is applied to guarantee that only one species of atoms is present at the growing surface at any instant. As a consequence, in a certain temperature range only one monolayer of a single type of atom should be found on the surface after each type of evaporation pulse. For the simulation of the ALE process we use a ( 11) surface, corresponding to a cross-section along a ( 100) plane in a real crystal as seen in Fig. 3, in this case atoms 1, 3,5 and so on are atoms of the same species. Additional atoms of the same kind in position 2 and 4 are bound very weakly and can easily re-evaporate if the substrate temperature is high enough. If we deposit a different kind of atom in position 2 or 4 then two strong bonds to the underlying monolayer will give a new stable surface. If a monolayer is not complete and an atom of a wrong species would be deposited in this monolayer only two weak bonds would be available for this atom, so it is very improbable that the atom is incorporated in the wrong layer. Notice that an UHV ALE process consists of a series of different thermodynamic growth processes. For the simulation we split the UHV ALE process into these thermodynamic processes and use the result of one process as the initial setup of the following process. In addition to our two-dimensional model described in Section 2.1 we have to change the particle species to be deposited after each growth cycle. This is equivalent to the experimental situation in an UHV ALE process. Within each of these ALE growth cycles, atoms of one species sufficient for a coverage of four layers were deposited. By using different values of r,,, for substrate and epilayer atoms we simulate lattice-mismatched ALE. 2.3. Three-dimensional
model
Also in the three-dimensional case we use a LJ potential to describe the interaction between the adatoms and the substrate. So, we can choose the binding length and the binding energy between particles of species i and j as described in Section 2.1.
Fig. 4. Representation of the angles cpand 6 in the xy plane (a) and in the perpendicular plane which contains the vector r12 (h).
To simulate a cubic primitive structure we introduced a directional LJ potential with four-fold symmetry, which seemed to be the easiest case with respect to vectorisation in the computing. For a cubic symmetry we used the following analytical expression for the binding energy between atoms of species i and j.
r-‘2-f(~01,91).~(~*P2,~292).~i,.r-6 =Eo,!y.
(3)
ff;
with f((p,6)
= [cos2 (29) +cos2 &sin2
(2~)]
.cos2 (26)
(4)
The definition of cp and 8 follows the usual convention for polar coordinates. The meaning of the angles is depicted in Fig. 4 for the case where the distance between the atoms is described by the vector ti2 and the orientation of the atoms is represented by the vectors e, and e2, respectively. The parameter Q is defined as in the two-dimensional case in Eq. (2). In principle, it is possible to simulate with our three-dimensional model the growth of compound epilayers on compound substrates or to study the incorporation of impurities in a host lattice, because we can set the pairwise binding energies and length for n kinds of atoms at will. To simulate the growth of a crystal, we started by placing three layers of 30 X 30 substrate atoms according to their interaction potential, which means in a closed packed hexagonal structure for the spherical LJ potential and in a cubic primitive structure for the fourfold symmetric directional LJ potential. 2.4. Presentation of the results We present the results in three different ways. As long as we are working in a two-dimensional cross-section, it is possible to show a picture of the whole calculated crystal lattice.
H. Sitter/Thin Solid Films 267 (I 995) 37-46
However, for simplicity we show typical parts of the lattice picture, while the system taken into account is much larger. In the case of the three-dimensional model we can show the lattice picture in a perspective view or in a top view, where usually the atoms of the top layer of the substrate are shown as light grey circles. The adatoms are shaded with increasing darkness depending on their distance from the substrate. So, we usually have to restrict our presentations of the threedimensional model calculations to only some monolayers which are cut from the whole lattice and depicted in different projections. To demonstrate the lattice quality we plot the normalised pair density as a function of the pair distance r, which is defined as the number of pairs with a pair distance within r and r t A r. The resulting figure shows a series of peaks, each corresponding to a translation vector in the two-dimensional or three-dimensional lattice. The horizontal axis is always scaled in units of reqU,oo, which gives us the possibility of identifying the lattice misfit and the long range order in the lattice.
3, Applications
41
epilayer (E,, x 0 = 2.0E0, x x, Eo,oo = EO,x x ) and exact lattice match, the result is shown in Fig. 7. The lattice pictures Figs. 6(a) and 7(a) show clearly either island or layer growth. Figs. 6(b) and 7(b) show the pair density for island and layer growth. The insets define the unit cells with the base
+
t
4
1
Fig. 5. Typical part of the simulated lattice of an epilayer with a misfit of 14% simulating the situation present in CdTe on GaAs. The other parameters arep=O.S,E~~.,,= Ea.xo = E 0.30 and kaT= 0.04 Eo.,,. The resulting dislocations are marked by arrows.
and discussion
By setting the parameter /3 between 0 and 1, we could choose between LJ and DLJ potential. Whenever we set /3 = 0 we obtained a hexagonal lattice structure, while for /3 = 0.5 the results showed a quadratic lattice. It is obvious that there must exist a transition from a hexagonal to a quadratic lattice, but this transition was not investigated in detail. 3.1. Results
of the
two-dimensional
model
1’0 test our model we simulated the growth of a structure with a mismatch of 14% using DLJ with p= 0.5 and rc’,qu,>x = 1II 4 ~equ.00 The binding energies were chosen as E,, x < = EC, x 0 = EO,>:x and T= 0.04 E,,, x x kB ‘. The result is depicted in Fig. 5, which shows a typical section of the lattice but the simulation was calculated for a substrate width of 80 atoms. The resulting dislocations are marked by arrows. It is evident that the mean distance between two dislocations is seven substrate atoms. It is thus in excellent agreement with the TEM observations of CdTe/GaAs heterostructures [ 181. In the next set of simulations we investigated the transition from island to layer growth. The substrate width was kept constant during all simulations with 80 atoms and the temperature was set to T=O.O4 E,,, x k; I. The ratio of the bindmg energy between substrate and epilayer atoms E,,, 0 to the binding energy between epilayer atoms E,,. x was varied in seven steps from 0.5: 1 to 2: 1. For all the simulated systems we found a strong influence of the ratio of the binding energies on the nucleation stage. Out of 48 simulations we show two typical results. For a weak binding between substrate and epilayer E,,, o = OSE,, x x, E,,oO = E,_ x x and a lattice-mismatched (5%) system, we obtained the results shown in Fig. 6. For a strong bond between substrate and
Fig. 6. Result of a simulation to demonstrate island growth by selecting a weak bond between epilayer and substrate. The parameters are p=O, r,,.xx =0.95 reqU.oo, E0.x~=0.5 Eo.xx; Eo.~~=Eo,.,: (a) typical part of the lattice, (b) pair density vs. distance in units of T,.~,“.~~,
#x
XYXXXrYXXXXXXYXXX~XXxx~~~~xxxx~ XXX~XxrX 9000000000000000000000000000000000000 0000000000000000000000000000000000000~00000000
pooooooooooOOOOOoOOOoooooooooooooooooooooooooooo
o-
--
(a) :
Y
4
Fig. 7. Result of a simulation to demonstrate layer growth by selecting a strong bond between epilayer and substrate. The parameters are /3=0, f&O = I$. x x : (a) typical part of the r,,,. x x = re,.Oov &x0=2 Eo.xx; lattice, (b) pair density vs. distance in units of rcqU.oO.
42
H. Sitter/Thin
Solid Films 267 (1995) 3746
vectors for the identification of the translation vectors. Within an island one can find pairs for almost all translation vectors while in a single layer there exist only pairs for translation vectors parallel to the interface. Therefore, we find peaks for the translation vectors (11)) (21), (22)) and (31) in Fig. 6(b), whereas in Fig. 7(b) the (50), (60), and (70) peaks are clearly resolved. The lattice mismatch of 5% is evident by a shift of all peaks in Fig. 6(b) in comparison with the scale, which is given in lattice constants of the substrate.
The macroscopic parameter which can be compared with the experiment is the growth rate given in monolayers per growth cycle. For all the simulations, we calculated this
(a)
3.2. Results of ALE simulation To investigate the influence of the substrate temperature on the ALE process, we simulate a system with a width of 80 atoms consisting of an elemental substrate and a binary epilayer. The lattice pictures presented in the following figures show only typical parts of the whole calculated systems. The strong bonds between unlike epilayer atoms and the strong bonds to the substrate were obtained by setting
(b)
E 0,x x =0.2Eo.x * E O,** =0.2&,x
(d)
*
E o,x o = 0.2Eo.00 E 0,x * = Eo,oo
(o)
E o.*o--E 0.00 At first, we deposited an amount of atoms of one kind, which corresponds to 4 monolayers (MLs) . After equilibration of the system, we changed the atom species and repeated the same procedure to finish one digital cycle (1 ALE growth cycle). For each simulation we performed 2 ALE growth cycles. The temperature was chosen in the range from T= 0.01-0.035 Eo,Oo kI; ‘. The lattice pictures of the ALE simulation for different substrate temperatures are shown in Fig. 8. At the lowest temperature, no layer growth was obtained (see Fig. 8(a) ) . With increasing temperature, the self-regulatory process became increasingly pronounced (see Figs. 8 (b)-8 (d) ) . If the temperature was chosen too high, some parts of the completed monolayers began to evaporate (see Figs. 8 (e) and 8(f)). This temperature dependence of the ALE process is also evident in the pair densities as shown in Fig. 9; parts (a) and (b) correspond to the lattice pictures of Figs. 8(a) and 8 (d) . In Fig. 9(a) only peaks with low indices can be found together with a homogeneous density for larger distances which shows that only short-range order is present. Contrary to the inferior lattice quality represented in Figs. 8(a) and 9(a) the lattice shown in Fig. 8(d) has a pair density with very pronounced peaks, as seen in Fig. 9(b) . Because the lattice consists of only four monolayers on a ( 11) substrate surface are seen, high-indexed peaks as (40) and (50) and SO on cannot appear in the pair density. However, a peak at 5r eqU,ooarises due to the existence of pairs with a (43) lattice vector.
Fig. 8. Cross-sections of the resulting lattices from ALE growth simulations for different temperatures given in units of !I,., x k, ‘: (a) T= 0.01, (b) T=0.017, (c) T=0.020, (d) T=0.022, (e) T=0.027 and (f) T=0.035.
Fig. 9. (a) Pair density for the lattice plotted in Fig. 8(a). for the lattice plotted in Fig. 8(d),
(b) Pair density
H. Sitter/Thin
Solid Films 267 (1995) 37-46
43
3.3. Results of the three-dimensional model
00..
0.00
..!
0,Ol
-!
0.02 Growth
temperature
0.03
ON
LEg,,,okBm’l
Fig. IO. Calculated growth rate in monolayer (ML) per cycle in the ALE simulation as a function of temperature. parameter by dividing the number of deposited particles, given In monolayers, by the number of growth cycles. The calculated growth rate as a function of temperature is shown in Fig. 10. It is clear that there exists a temperature range from 0.02 to 0.03 E o,OOk; I, where a constant growth rate of 1 ML cycle- ’ can be obtained, reflecting the self-regulatory regime of ALE growth. The lower limit of the plateau is given by the incomplete re-evaporation of weakly bound atoms. The high temperature limit of the plateau is determined by the evaporation of the compound. The temperature dependence of the growth rate obtained by the MC simulation is in excellent agreement with the experimental results reported in Ref. [ 191. In another set of computer simulations we investigate the influence of a tilted substrate surface in connection with a large lattice mismatch between substrate and epilayer. This is achieved by setting r,,,, f . = rqU,x x = 0.9r,,,,, at T= 0.102E o,OOk; ‘, where we obtained layer-by-layer growth
for lattice-matched systems on flat substrates. No significant influence on the growth rate is found; however, a significant change of the lattice structure occurred. This can be seen in the lattice picture shown in Fig. 11 (a), where the atomic rows of the epilayer are tilted by about 5” with respect to the atomic rows Iof the substrate which stays quadratic. The epilayer is strained to meet the horizontal lattice constant of the substrate and all additional lattice plane is inserted close to the step in the substrate. As a consequence, the vertical lattice constant of the epilayer shrinks. As a result, the quadratic lattice is strained to a rhombic structure. This can be seen even better in the pair density calculated for the epilayer as shown in Fig. 11 (b:l . The ( IO) peak of the epilayer is found at approxiaccording to the lattice mismatch. At the mately 0.9 requ.OO, distance of ( 11) pairs for a relaxed lattice of the epilayer, no peak 1s found. Instead of that, owing to the changed structure, we obtain two ( 11) peaks: one for the ( 11) 11direction parallel to the interface at the position of the (11) peak of the unstrained substrate and one for the ( 11) I direction normal to the interface. With the values found for the pair distances in the ( 11) 11and ( 11) _L direction the angle for the rhombic distortion can be evaluated and gives a value of 5” which is also reflected in the lattice picture of Fig. 11 (a).
As a first test of the three-dimensional model using the LJ potential we investigated the influence of the binding energies between equal and unequal atoms on the growth mechanism. If we set Eo,Oo:Eo, X X:Eo,Ox = 2: 1: 1 and reqU,oo= rqU,x x we simulate a lattice-matched system where the binding energy between adatoms and the substrate are only half of the binding energy of the equal pairs. Fig. 12 shows the result of the simulation after the deposition of a number of particles, which corresponds to a full monolayer coverage. The upper layer of the substrate is represented by the ordered layer in Fig. 12 while the adatoms are shown with a different grey scale. It can be clearly seen that island growth occurs at these parameters. If we set the contrary relation between the binding energies E 0.00..E 0.x x *Eo,ox = 1: 1:2 and keep the other parameters . constant we simulate a lattice-matched system with a much stronger bond of the adatoms to the substrate than among each other. The obtained result is depicted in Fig. 13 in a top view to the crystal lattice. The regular array of light grey circles represents the top layer of the substrate, the dark grey and black circles depict the adatoms in the first and second monolayer above the substrate surface. In comparison with the previous situation we obtain layer growth. However, there are two different places on the substrate surface which can be occupied by the adatoms, the position of the nearest and the position of the next-nearest lattice site. That means that twins are formed in this growth direction in a layer-by-layer growth mode. It is interesting to observe that the adatoms of the second monolayer (black circles in Fig. 13) are deposited mainly at the twin boundaries where they find energetically preferable positions.
Fig. Il. (a) Lattice picture of a strained epilayer on a tilted substrate. (b) Pair density of the epilayer depicted in (a),
44
H. Sitter /Thin Solid Films 267 (1995) 3746
Fig, 12. Perspective lattice picture of island growth.
Fig. 13. Top view on a substrate
(light grey circles) covered with one monolayer
of adatoms (dark grey and black circles) showing twin formation
H. Sitter/Thin Solid Films 267 (1995) 37-46
Fig. 14 Top view of a substrate covered with one monolayer the directions of the strongest binding energies.
4s
of adatoms with four-fold symmetry in the directional LJ potential. The bars on the atoms represent
To investigate the stabilisation of the growth mechanism on vizinal surfaces we simulated the growth on substrates with atomic steps. When we chose the parameters to simulate &,OO = G. x x = &,o x, rqu.oo = r,,,, x x ) homoepitaxy we ohtain perfect layer-by-layer growth without any twins proving the experimentally well-known stabilisation effect of vicinal substrate surfaces. However, when we changed the radii of the adatoms keeping the binding energies constant we obtained different results. In the case requ,ogcqu. x x = 1: 1.03, which describes a lattice mismatch of 3%, the simulation results in stacking faults which are generated in the first monolayer and overgrown in the next monolayer when the terraces of the tilted structure are completed [ 17). If we simulate a lattice mismatch of 9% we obtain a destabilisation of the growth mode, because the energetictlly preferable lattice sites at the terraces cannot overcome the strain in the layers and additional lattice planes are incorporated to accommodate the lattice mismatch [ 171. Finally, we skipped to the three-dimensional directional LJ potential with the four-fold symmetry, as described in Section 2.3, to simulate a cubic primitive lattice. To show the directions of the strongest bonds among the atoms we added six white bars to any atom in the lattice pictures representing the directions of the bonds. The result of simulation with the parameters E,,oo:Eo, x x :Eo,o x = 1: 1:2 and rqu~oo = r,,,. x x is shown in a top view in Fig. 14. The well-ordered atoms characterised by the exactly aligned bonds represent the top layer of the substrate. Layer growth can be observed due to the chosen ratio in the binding energies but the adatoms order themselves only partly in the correct positions of the cubic
lattice as clearly seen in Fig. 14. This result is in agreement with the font that no crystals with covalent bonds exist in that symmetry class.
4. Conclusion In this paper we have demonstrated the use of a continuousspace Monte-Carlo method for studies of growth mechanisms in epitaxy. The restriction to a two-dimensional cross-section allows only a simplified picture of the reality, however, by using different binding energies for the deposited atoms island and layer growth was obtained as expected. For latticemismatched structures the incorporation of dislocations and strain could be simulated on flat substrates and also on vicinal surfaces, which can be achieved only with acontinuous space model. Also the ALE growth process could be simulated in the two-dimensional model.
Acknowledgements The author would like to thank W. Plotz and M. Tagwerker who wrote the simulation programs and implemented them on different computers. My sincere appreciation is extended to C. Leitner for typing the manuscript and to D. Stifter for the fabrication of the computer plots. This work was partly supported by the Fonds zur Fbrderung der wissenschaftlichen Forschung in esterreich, project number P9695.
46
H. Sitter/Thin Solid Films 267 (I 995) 37-46
References
[9] D.A. Faux, G. Gaynor, C.L. Carson, C.K. Hall and J. Bemholc, Phys. Rev. B, 42 ( 1990) 2914. [ 101 B.W. Dodson and P.A. Taylor, Phys. Rev. B, 34 ( 1986) 2112.
[ 1J M.A. Herman and H. Sitter, in M.B. Panish (ed.), Molecular beam epitaxy-fundamentals and current status, SpringerSeries in Materials Science, Vol. 7, Springer, Berlin, 1988. [2] W. Faschinger and H. Sitter, /. Cryst. Growth, 99 ( 1990) 566. [3] N. Metropolis, A.W. Rosenbluth, N.M. Rosenbluth, A.H. Teller and E. Teller, J. Chem. Phys., 21 (1953) 1087. [4] H. Mtiller-Kmmbhaar. in K. Binder (ed.), Monte Carlo Methodr in Statistical Physics, Springer, Berlin, 1979. [5] R.H. Swendsen, Phys. Rev. B, 1.5 (1976) 5421. [6] D. Kashchiev, J.P. van der Eerden and C. van Leeuwen, J. Cryst. Growth, 40 (1977) 47. [7] R.H. Swendsen, P.J. Kortman, D.P. Landauand H. Miiller-Krumbhaar. J. Cryst. Growth, 35 (1977) 73. [ 81 C.L. Carson and J. Bemholc, Appl. Phys. Lat.. 56 ( 1990) 1971.
[ll I B.W. Dodson and P.A. Taylor, Phys. Rev. B, 36 (1987) 1355. [ 121 B.W. Dodson, Phys. Rev. E, 30 (1984) 3545.
[ 131 M. Schneider, I.K. Schuller and A. Rahman, Phys. Rev. B, 36 ( 1987) 1340.
[ 141 W.M. Plotz, K. Hingerl and H. Sitter, J. Cryst. Growth, I15 (1991) 186. [ 151 W.M. Plotz, K. Hingerl and H. Sitter, Phys. Rev. B, 45 ( 1992) 12122. [16] W.M. Plotz, K. Hingerl and H. Sitter, Semicond. Sci. Technol., 9 ( 1994) 2224. [ 171 M. Tagwerker, W.M. Plotzand H. Sitter, J. Cryst. Growth, 146 ( 1995) 220. [ 181 F.A. Ponce, G.B. Anderson and J.M. Ballinger, Surf: Sci., 168 (1986) 564.
[ 191 H. Sitter, W. Faschinger, Thin Solid Films, 225 ( 1993) 250.
Thin Solid Films 267 ( 1995) 3746
MBE growth mechanisms; studies by Monte-Carlo simulation
’
H. Sitter Johannes
Kepler Universitiir Linz, Institute oj’Experimenta1 Physics, A-4040 Linz, Ausrria
-. Abstract A hlonte-Carlo-based method for simulation of epitaxial crystal growth is reported. The improvement of this concept in comparison with based on the solid-on-solid model is the possibility of simulating growth without the assumption of a discrete crystal lattice structure. The used continuous space model requires only the assessment of the Hamiltonian. In the first step the model was applied for studying growth in a two-dimensional cross-section through the substrate and the epilayer with a width of up to 80 atoms. Two different potentials can be chosen for the simulation: (a) Lennard-Jones and (b) directional Lennard-Jones to describe covalent systems. In addition to the selection of the potential, we can set the following parameters for different kinds of atoms: radii of the potential minimum, binding energy between the atoms for all possible pairings, growth temperature and substrate structure. Island and layer growth in the nucleation stage were investigated by selecting different binding energies of the epilayer atoms to the substrate atoms. The epitaxial growth of strained layers was simulated to study the incorporation of dislocations at the interface. The continuous space model also allowed calculation of the rhombic distortion of a strained layer grown on a tilted substrate with atomic steps. As a further extension of the two-dimensional model the re-evaporation of particles was allowed which is essential for the simulation of atomic-layer epitaxy. The obtained results for the growth rate in atomic layer epitaxy as a function of substrate temperature can be compared with experimental findings. In .I second attempt we extended the continuous space model to three dimensions based on the same assumption as tested in the twodimensional model. The three-dimensional model allowed the same simulations for islands and layer growth but also the selection of vicinal surfaces of the substrates with respect to the usually used ( 111 )-oriented surfaces. The simulation showed a clear stabilisation of the growth direct,on by the atomic steps on the substrate and even in lattice -matched structures the formation of twins was suppressed, the simulations
Keywc~rd~: Adsorption; Deposition process; Growth mechanism; Nucleation
1. Introduction Acvances
in crystal
growth
technologies,
in particular
in
beam epitaxy (MBE) [ 1 ] and atomic layer epitaxy (ALE) [ 21 have stimulated the search for realistic Irheoretical models of crystal growth. So far, phenomenological theories using thermodynamics were used to describe processes such as condensation and re-evaporation or misfit relaxation and incorporation of dislocations. With the ard of modern computer technology it became possible to realixe ab-initio calculations with realistic interaction potentials. These microscopic models increase the understanding of the physics of crystal growth and epitaxy. In principle, there are two possible simulation methods. First. one can solve the equation of motion of the molecules to track the system. This method is calledmoleculardynamics (MD). Second, it is possible to calculate the properties of a system at equilibrium from its distribution function in the the field of molecular
’ Presented as an invited talk at the Workshop on MBE-Growth and Technology, Warsaw, Poland, December 1994.
Physics
0040-6090/95/$09.50 0 1995 Elsevier Science S.A. All rights reserved .SSDIOO40-6000(95)06593-8
phase space. Because this is done with the use of random sampling, it is called the Monte-Carlo method (MC). The name Monte-Carlo method is used for all models requiring random numbers, but in this work we restrict ourselves to algorithms similar to the so-called important sampling method as described by Metropolis et al. [ 3 I. For both, MD and MC, one only needs to know or assume the Hamiltonian of the system. The main difference between MD and MC is the fact that MD can describe dynamic systems, while by MC no real time behaviour can be calculated because the algorithm relies on thermodynamic equilibrium. However, the MC calculations are performed in such a way that the system is disturbed by insertion of a new particle, e.g. an adatom which can be incorporated into the growing crystal, and then equilibrated until the change in the potential energy by more MC steps is negligible. This procedure can be repeated several times, and in this way dynamic processes with an artificial time scale can also be simulated by MC methods. On the other hand, the solution of the coupled equations of motion for any particle of the system in MD restricts the number of particles
38
H. Sitter / Thin Solid Films 267 (I 995) 3746
and also the range of real time simulation because of limitations in CPU time. Especially if the interaction potential is a complicated function, the MC method is favourable, because only one scalar number, the change of the potential energy, has to be calculated. In a huge number of computer simulations the growth of lattice-matched crystals is described by an Ising model [47]. All these models use a predefined lattice structure. Although these models are quite simple, they have been used successfully to study the surface roughening with increasing temperature [ 51, the transition between island and layer growth [ 61 and propagation of spiral dislocations [ 71. More sophisticated simulations based on Ising models allow the simulation of lattice-mismatched structures [ 8,9]. However, calculations using continuous space simulations are quite rare [ 10-17 1. Due to the complexity of continuous space models, their application has to be restricted to very small three-dimensional crystals or two-dimensional crosssections. Continuous space models using MD and MC methods provide some advantages in comparison with discrete space models, because the assumption of the Hamiltonian is sufficient to describe the physical behaviour of the system and the models are somewhat closer to reality. The most often used potential for the simulation of crystal growth is of the Lennard-Jones (LJ) type, which has been used for MC [ 10,15-171 and MD [ 111 simulations. For real semiconductor structures, a covalent bonding mechanism should be taken into account. This was achieved by a directional LJ (DLJ) potential, where the particles are defined by their locations and their orientation [ 12,15-171. By defining for each atom three bonds in an angle of 120” [ 121 or by using four bonds in an two-dimensional projection [ 161, this model was used to predict the stability range in strained layer superlattices. Another way to describe covalent bonds is to set up an interaction potential consisting of a LJ term together with a three-particle part, which favours a tetrahedral atomic configuration. With this model the epitaxial growth of silicon on silicon was simulated by MD calculations [ 131. At present, there is no general rule for the choice of the simulation model. All the models mentioned above have advantages and disadvantages. With the model presented in this paper it is principally possible to tackle the problem of lattice mismatch and miscut, substrate orientation, substratetemperature dependence, surface morphology, incorporation of dislocations and the flip from island to layer growth for heteroepitaxy, including the growth of binary and ternary compounds on binary substrates. In Section 2 we will discuss the Hamiltonian and certain restrictions for the simulation due to limitations in CPU time. In Section 3 we present some of the already achieved results. The most important results of our calculations are: the incorporation of misfit dislocations in a GaAs/CdTe heterostructure, and the change from island to layer growth as a function of binding energies between substrate and epilayer atoms in the two-and three-dimensional model. The prediction of the stability range in an ALE
process could be calculated in the two-dimensional and is comparable with experimental results.
2. Description
model
of the model
2. I. Two-dimensional
model
We use a continuous-space MC method to simulate epitaxial crystal growth in a two-dimensional cross-section through the substrate and the epilayer parallel to a (001) plane of L zinc-blende structure. In this cross-section the four tetrahedra: bonds of each atom are oriented in the crystal lattice as shown in Fig. 1 (a). It is evident that the projection of the tetragonal zinc-blende lattice into the plane gives a four-fold symmetry. This restriction to two dimensions was applied to the simulation to save CPU time during the first calculations in order to test the model. In our model each particle is represented by two twodimensional vectors in the plane of the cross-section, namely its location r and its orientation e with (e 1= 1. The pairwise interaction potential is characterised by three parameters: the binding energy E,, the binding length r,,, and the parameter p which describes the covalent character of the bond. The potential energy of two interacting particles is given by Eqn. (1).
‘\ ~100) (001)
Surface
(4
atom 2 Fig. 1. (a) (001) surface of a zinc-blende structure with the usual crystallographic directions and the corresponding bonds. (b) Schematic drawing of two particles at a distance of rj2, with their orientations e,, eZ and the definition of the angles cp,. cpz,
H. Sifter/ Thin Solid Films 267 (1995) 3746
39
r= Irl
r=r,-r2
r.e, cos
((0”)
=7
f(cp,)=cos2 o= 2$,6
Fig. Pairwise binding energy El2 (r. 9, 0) in units of EO, plotted in Cartesim coordinates for three different values of p.
I
/co~“-(1-P)~-“-Pf((p,>f((p2)m-6 r
a2
r
CUlOff
0 r 2 rculoff
(1) The pairwise interaction energy is a function of the variables r, 2, ‘p, and (p2which are shown in Fig. 1 (b) . The distance r12, the angular dependencef, the parameter (Y,and the cutoff radius rcut_Off are taken as
(2cp,)withnE(1,2}
(2)
rCuloff= 2Jr,,
A graphical representation of the potential energy for three different values of the parameter /3 is given in Fig. 2. The interaction potential can be tuned by the parameter /3 from the pure LJ potential with p = 0 (Fig. 2(a)) to the pure directional LJ potential with /3 = 1 (Fig. 2(c) ) . Since different species of atoms are allowed, E,, req,,, and j3 are specific for the atom species involved in the interaction. Therefore, we denote these parameters by E,,,j, requ,;, and PO to describe an interaction of an atom of kind i with an atom of kindj. To simulate the growth of compound epilayers on compound substrates our model allows one to select I?,,,,, rqu.,, and p,, for n different kinds of atoms. That means, in principle, 3n( n + 1) /2 parameters describing the interaction of atoms can be set at will. To save CPU time we set rq “,,, = ( rCq“,,, + r,.,,) /2 and pi,= pii = p,j= p. From now on the symbols 0, X and * denote the three kinds of atoms according to the representation of atoms in the lattice picture. By definition, 0 indicates a substrate atom while X and * indicate an atom of the epilayer. To simulate the growth of a crystal we use a substrate width of 80 atoms and periodic-boundary conditions in the horizontal direction. We start by placing some layers of substrate atoms in the distance of rqu OO. The substrate is equilibrated by means of a conventional’important sampling MC method [ 31; this means that one particle is selected at random and moved to a position within a rectangle with a horizontal length of 2.5 rqu,OO and a height of roqu,OO,The orientation e of the particle is chosen at random also. This move and rotation is accepted with the probability p = exp ( - AHIKJ). This procedure, commonly called the MC step, is repeated until an acceptable approach to the equilibrium is achieved. After the substrate is equilibrated, another particle is placed near the surface at a random position and the system is equilibrated again. After some particles are deposited, expectation values of the position of each particle are calculated. The difference between the actual and the average position is used as an indication of the stability of the lattice. This procedure is repeated until the required number of particles is deposited.
2.2. Model for atomic-layer epitaxy
Ultra-high vacuum (UHV) ALE of a binary compound on an elemental substrate is characterised by the following conditions.
H. Sitter/Thin
Solid Films 267 (1995) 3746
Fig. 3. Two-dimensional lattice used for the simulations of UHV-ALE.
(i> UHV ALE is only possible for binary systems where the bonds between the different atomic species are stronger than the bonds between the same atomic species, which is true for many II-VI compounds. (ii) In contrast to MBE growth, in an UHV ALE growth process a shutter sequence is applied to guarantee that only one species of atoms is present at the growing surface at any instant. As a consequence, in a certain temperature range only one monolayer of a single type of atom should be found on the surface after each type of evaporation pulse. For the simulation of the ALE process we use a ( 11) surface, corresponding to a cross-section along a ( 100) plane in a real crystal as seen in Fig. 3, in this case atoms 1, 3,5 and so on are atoms of the same species. Additional atoms of the same kind in position 2 and 4 are bound very weakly and can easily re-evaporate if the substrate temperature is high enough. If we deposit a different kind of atom in position 2 or 4 then two strong bonds to the underlying monolayer will give a new stable surface. If a monolayer is not complete and an atom of a wrong species would be deposited in this monolayer only two weak bonds would be available for this atom, so it is very improbable that the atom is incorporated in the wrong layer. Notice that an UHV ALE process consists of a series of different thermodynamic growth processes. For the simulation we split the UHV ALE process into these thermodynamic processes and use the result of one process as the initial setup of the following process. In addition to our two-dimensional model described in Section 2.1 we have to change the particle species to be deposited after each growth cycle. This is equivalent to the experimental situation in an UHV ALE process. Within each of these ALE growth cycles, atoms of one species sufficient for a coverage of four layers were deposited. By using different values of r,,, for substrate and epilayer atoms we simulate lattice-mismatched ALE. 2.3. Three-dimensional
model
Also in the three-dimensional case we use a LJ potential to describe the interaction between the adatoms and the substrate. So, we can choose the binding length and the binding energy between particles of species i and j as described in Section 2.1.
Fig. 4. Representation of the angles cpand 6 in the xy plane (a) and in the perpendicular plane which contains the vector r12 (h).
To simulate a cubic primitive structure we introduced a directional LJ potential with four-fold symmetry, which seemed to be the easiest case with respect to vectorisation in the computing. For a cubic symmetry we used the following analytical expression for the binding energy between atoms of species i and j.
r-‘2-f(~01,91).~(~*P2,~292).~i,.r-6 =Eo,!y.
(3)
ff;
with f((p,6)
= [cos2 (29) +cos2 &sin2
(2~)]
.cos2 (26)
(4)
The definition of cp and 8 follows the usual convention for polar coordinates. The meaning of the angles is depicted in Fig. 4 for the case where the distance between the atoms is described by the vector ti2 and the orientation of the atoms is represented by the vectors e, and e2, respectively. The parameter Q is defined as in the two-dimensional case in Eq. (2). In principle, it is possible to simulate with our three-dimensional model the growth of compound epilayers on compound substrates or to study the incorporation of impurities in a host lattice, because we can set the pairwise binding energies and length for n kinds of atoms at will. To simulate the growth of a crystal, we started by placing three layers of 30 X 30 substrate atoms according to their interaction potential, which means in a closed packed hexagonal structure for the spherical LJ potential and in a cubic primitive structure for the fourfold symmetric directional LJ potential. 2.4. Presentation of the results We present the results in three different ways. As long as we are working in a two-dimensional cross-section, it is possible to show a picture of the whole calculated crystal lattice.
H. Sitter/Thin Solid Films 267 (I 995) 37-46
However, for simplicity we show typical parts of the lattice picture, while the system taken into account is much larger. In the case of the three-dimensional model we can show the lattice picture in a perspective view or in a top view, where usually the atoms of the top layer of the substrate are shown as light grey circles. The adatoms are shaded with increasing darkness depending on their distance from the substrate. So, we usually have to restrict our presentations of the threedimensional model calculations to only some monolayers which are cut from the whole lattice and depicted in different projections. To demonstrate the lattice quality we plot the normalised pair density as a function of the pair distance r, which is defined as the number of pairs with a pair distance within r and r t A r. The resulting figure shows a series of peaks, each corresponding to a translation vector in the two-dimensional or three-dimensional lattice. The horizontal axis is always scaled in units of reqU,oo, which gives us the possibility of identifying the lattice misfit and the long range order in the lattice.
3, Applications
41
epilayer (E,, x 0 = 2.0E0, x x, Eo,oo = EO,x x ) and exact lattice match, the result is shown in Fig. 7. The lattice pictures Figs. 6(a) and 7(a) show clearly either island or layer growth. Figs. 6(b) and 7(b) show the pair density for island and layer growth. The insets define the unit cells with the base
+
t
4
1
Fig. 5. Typical part of the simulated lattice of an epilayer with a misfit of 14% simulating the situation present in CdTe on GaAs. The other parameters arep=O.S,E~~.,,= Ea.xo = E 0.30 and kaT= 0.04 Eo.,,. The resulting dislocations are marked by arrows.
and discussion
By setting the parameter /3 between 0 and 1, we could choose between LJ and DLJ potential. Whenever we set /3 = 0 we obtained a hexagonal lattice structure, while for /3 = 0.5 the results showed a quadratic lattice. It is obvious that there must exist a transition from a hexagonal to a quadratic lattice, but this transition was not investigated in detail. 3.1. Results
of the
two-dimensional
model
1’0 test our model we simulated the growth of a structure with a mismatch of 14% using DLJ with p= 0.5 and rc’,qu,>x = 1II 4 ~equ.00 The binding energies were chosen as E,, x < = EC, x 0 = EO,>:x and T= 0.04 E,,, x x kB ‘. The result is depicted in Fig. 5, which shows a typical section of the lattice but the simulation was calculated for a substrate width of 80 atoms. The resulting dislocations are marked by arrows. It is evident that the mean distance between two dislocations is seven substrate atoms. It is thus in excellent agreement with the TEM observations of CdTe/GaAs heterostructures [ 181. In the next set of simulations we investigated the transition from island to layer growth. The substrate width was kept constant during all simulations with 80 atoms and the temperature was set to T=O.O4 E,,, x k; I. The ratio of the bindmg energy between substrate and epilayer atoms E,,, 0 to the binding energy between epilayer atoms E,,. x was varied in seven steps from 0.5: 1 to 2: 1. For all the simulated systems we found a strong influence of the ratio of the binding energies on the nucleation stage. Out of 48 simulations we show two typical results. For a weak binding between substrate and epilayer E,,, o = OSE,, x x, E,,oO = E,_ x x and a lattice-mismatched (5%) system, we obtained the results shown in Fig. 6. For a strong bond between substrate and
Fig. 6. Result of a simulation to demonstrate island growth by selecting a weak bond between epilayer and substrate. The parameters are p=O, r,,.xx =0.95 reqU.oo, E0.x~=0.5 Eo.xx; Eo.~~=Eo,.,: (a) typical part of the lattice, (b) pair density vs. distance in units of T,.~,“.~~,
#x
XYXXXrYXXXXXXYXXX~XXxx~~~~xxxx~ XXX~XxrX 9000000000000000000000000000000000000 0000000000000000000000000000000000000~00000000
pooooooooooOOOOOoOOOoooooooooooooooooooooooooooo
o-
--
(a) :
Y
4
Fig. 7. Result of a simulation to demonstrate layer growth by selecting a strong bond between epilayer and substrate. The parameters are /3=0, f&O = I$. x x : (a) typical part of the r,,,. x x = re,.Oov &x0=2 Eo.xx; lattice, (b) pair density vs. distance in units of rcqU.oO.
42
H. Sitter/Thin
Solid Films 267 (1995) 3746
vectors for the identification of the translation vectors. Within an island one can find pairs for almost all translation vectors while in a single layer there exist only pairs for translation vectors parallel to the interface. Therefore, we find peaks for the translation vectors (11)) (21), (22)) and (31) in Fig. 6(b), whereas in Fig. 7(b) the (50), (60), and (70) peaks are clearly resolved. The lattice mismatch of 5% is evident by a shift of all peaks in Fig. 6(b) in comparison with the scale, which is given in lattice constants of the substrate.
The macroscopic parameter which can be compared with the experiment is the growth rate given in monolayers per growth cycle. For all the simulations, we calculated this
(a)
3.2. Results of ALE simulation To investigate the influence of the substrate temperature on the ALE process, we simulate a system with a width of 80 atoms consisting of an elemental substrate and a binary epilayer. The lattice pictures presented in the following figures show only typical parts of the whole calculated systems. The strong bonds between unlike epilayer atoms and the strong bonds to the substrate were obtained by setting
(b)
E 0,x x =0.2Eo.x * E O,** =0.2&,x
(d)
*
E o,x o = 0.2Eo.00 E 0,x * = Eo,oo
(o)
E o.*o--E 0.00 At first, we deposited an amount of atoms of one kind, which corresponds to 4 monolayers (MLs) . After equilibration of the system, we changed the atom species and repeated the same procedure to finish one digital cycle (1 ALE growth cycle). For each simulation we performed 2 ALE growth cycles. The temperature was chosen in the range from T= 0.01-0.035 Eo,Oo kI; ‘. The lattice pictures of the ALE simulation for different substrate temperatures are shown in Fig. 8. At the lowest temperature, no layer growth was obtained (see Fig. 8(a) ) . With increasing temperature, the self-regulatory process became increasingly pronounced (see Figs. 8 (b)-8 (d) ) . If the temperature was chosen too high, some parts of the completed monolayers began to evaporate (see Figs. 8 (e) and 8(f)). This temperature dependence of the ALE process is also evident in the pair densities as shown in Fig. 9; parts (a) and (b) correspond to the lattice pictures of Figs. 8(a) and 8 (d) . In Fig. 9(a) only peaks with low indices can be found together with a homogeneous density for larger distances which shows that only short-range order is present. Contrary to the inferior lattice quality represented in Figs. 8(a) and 9(a) the lattice shown in Fig. 8(d) has a pair density with very pronounced peaks, as seen in Fig. 9(b) . Because the lattice consists of only four monolayers on a ( 11) substrate surface are seen, high-indexed peaks as (40) and (50) and SO on cannot appear in the pair density. However, a peak at 5r eqU,ooarises due to the existence of pairs with a (43) lattice vector.
Fig. 8. Cross-sections of the resulting lattices from ALE growth simulations for different temperatures given in units of !I,., x k, ‘: (a) T= 0.01, (b) T=0.017, (c) T=0.020, (d) T=0.022, (e) T=0.027 and (f) T=0.035.
Fig. 9. (a) Pair density for the lattice plotted in Fig. 8(a). for the lattice plotted in Fig. 8(d),
(b) Pair density
H. Sitter/Thin
Solid Films 267 (1995) 37-46
43
3.3. Results of the three-dimensional model
00..
0.00
..!
0,Ol
-!
0.02 Growth
temperature
0.03
ON
LEg,,,okBm’l
Fig. IO. Calculated growth rate in monolayer (ML) per cycle in the ALE simulation as a function of temperature. parameter by dividing the number of deposited particles, given In monolayers, by the number of growth cycles. The calculated growth rate as a function of temperature is shown in Fig. 10. It is clear that there exists a temperature range from 0.02 to 0.03 E o,OOk; I, where a constant growth rate of 1 ML cycle- ’ can be obtained, reflecting the self-regulatory regime of ALE growth. The lower limit of the plateau is given by the incomplete re-evaporation of weakly bound atoms. The high temperature limit of the plateau is determined by the evaporation of the compound. The temperature dependence of the growth rate obtained by the MC simulation is in excellent agreement with the experimental results reported in Ref. [ 191. In another set of computer simulations we investigate the influence of a tilted substrate surface in connection with a large lattice mismatch between substrate and epilayer. This is achieved by setting r,,,, f . = rqU,x x = 0.9r,,,,, at T= 0.102E o,OOk; ‘, where we obtained layer-by-layer growth
for lattice-matched systems on flat substrates. No significant influence on the growth rate is found; however, a significant change of the lattice structure occurred. This can be seen in the lattice picture shown in Fig. 11 (a), where the atomic rows of the epilayer are tilted by about 5” with respect to the atomic rows Iof the substrate which stays quadratic. The epilayer is strained to meet the horizontal lattice constant of the substrate and all additional lattice plane is inserted close to the step in the substrate. As a consequence, the vertical lattice constant of the epilayer shrinks. As a result, the quadratic lattice is strained to a rhombic structure. This can be seen even better in the pair density calculated for the epilayer as shown in Fig. 11 (b:l . The ( IO) peak of the epilayer is found at approxiaccording to the lattice mismatch. At the mately 0.9 requ.OO, distance of ( 11) pairs for a relaxed lattice of the epilayer, no peak 1s found. Instead of that, owing to the changed structure, we obtain two ( 11) peaks: one for the ( 11) 11direction parallel to the interface at the position of the (11) peak of the unstrained substrate and one for the ( 11) I direction normal to the interface. With the values found for the pair distances in the ( 11) 11and ( 11) _L direction the angle for the rhombic distortion can be evaluated and gives a value of 5” which is also reflected in the lattice picture of Fig. 11 (a).
As a first test of the three-dimensional model using the LJ potential we investigated the influence of the binding energies between equal and unequal atoms on the growth mechanism. If we set Eo,Oo:Eo, X X:Eo,Ox = 2: 1: 1 and reqU,oo= rqU,x x we simulate a lattice-matched system where the binding energy between adatoms and the substrate are only half of the binding energy of the equal pairs. Fig. 12 shows the result of the simulation after the deposition of a number of particles, which corresponds to a full monolayer coverage. The upper layer of the substrate is represented by the ordered layer in Fig. 12 while the adatoms are shown with a different grey scale. It can be clearly seen that island growth occurs at these parameters. If we set the contrary relation between the binding energies E 0.00..E 0.x x *Eo,ox = 1: 1:2 and keep the other parameters . constant we simulate a lattice-matched system with a much stronger bond of the adatoms to the substrate than among each other. The obtained result is depicted in Fig. 13 in a top view to the crystal lattice. The regular array of light grey circles represents the top layer of the substrate, the dark grey and black circles depict the adatoms in the first and second monolayer above the substrate surface. In comparison with the previous situation we obtain layer growth. However, there are two different places on the substrate surface which can be occupied by the adatoms, the position of the nearest and the position of the next-nearest lattice site. That means that twins are formed in this growth direction in a layer-by-layer growth mode. It is interesting to observe that the adatoms of the second monolayer (black circles in Fig. 13) are deposited mainly at the twin boundaries where they find energetically preferable positions.
Fig. Il. (a) Lattice picture of a strained epilayer on a tilted substrate. (b) Pair density of the epilayer depicted in (a),
44
H. Sitter /Thin Solid Films 267 (1995) 3746
Fig, 12. Perspective lattice picture of island growth.
Fig. 13. Top view on a substrate
(light grey circles) covered with one monolayer
of adatoms (dark grey and black circles) showing twin formation
H. Sitter/Thin Solid Films 267 (1995) 37-46
Fig. 14 Top view of a substrate covered with one monolayer the directions of the strongest binding energies.
4s
of adatoms with four-fold symmetry in the directional LJ potential. The bars on the atoms represent
To investigate the stabilisation of the growth mechanism on vizinal surfaces we simulated the growth on substrates with atomic steps. When we chose the parameters to simulate &,OO = G. x x = &,o x, rqu.oo = r,,,, x x ) homoepitaxy we ohtain perfect layer-by-layer growth without any twins proving the experimentally well-known stabilisation effect of vicinal substrate surfaces. However, when we changed the radii of the adatoms keeping the binding energies constant we obtained different results. In the case requ,ogcqu. x x = 1: 1.03, which describes a lattice mismatch of 3%, the simulation results in stacking faults which are generated in the first monolayer and overgrown in the next monolayer when the terraces of the tilted structure are completed [ 17). If we simulate a lattice mismatch of 9% we obtain a destabilisation of the growth mode, because the energetictlly preferable lattice sites at the terraces cannot overcome the strain in the layers and additional lattice planes are incorporated to accommodate the lattice mismatch [ 171. Finally, we skipped to the three-dimensional directional LJ potential with the four-fold symmetry, as described in Section 2.3, to simulate a cubic primitive lattice. To show the directions of the strongest bonds among the atoms we added six white bars to any atom in the lattice pictures representing the directions of the bonds. The result of simulation with the parameters E,,oo:Eo, x x :Eo,o x = 1: 1:2 and rqu~oo = r,,,. x x is shown in a top view in Fig. 14. The well-ordered atoms characterised by the exactly aligned bonds represent the top layer of the substrate. Layer growth can be observed due to the chosen ratio in the binding energies but the adatoms order themselves only partly in the correct positions of the cubic
lattice as clearly seen in Fig. 14. This result is in agreement with the font that no crystals with covalent bonds exist in that symmetry class.
4. Conclusion In this paper we have demonstrated the use of a continuousspace Monte-Carlo method for studies of growth mechanisms in epitaxy. The restriction to a two-dimensional cross-section allows only a simplified picture of the reality, however, by using different binding energies for the deposited atoms island and layer growth was obtained as expected. For latticemismatched structures the incorporation of dislocations and strain could be simulated on flat substrates and also on vicinal surfaces, which can be achieved only with acontinuous space model. Also the ALE growth process could be simulated in the two-dimensional model.
Acknowledgements The author would like to thank W. Plotz and M. Tagwerker who wrote the simulation programs and implemented them on different computers. My sincere appreciation is extended to C. Leitner for typing the manuscript and to D. Stifter for the fabrication of the computer plots. This work was partly supported by the Fonds zur Fbrderung der wissenschaftlichen Forschung in esterreich, project number P9695.
46
H. Sitter/Thin Solid Films 267 (I 995) 37-46
References
[9] D.A. Faux, G. Gaynor, C.L. Carson, C.K. Hall and J. Bemholc, Phys. Rev. B, 42 ( 1990) 2914. [ 101 B.W. Dodson and P.A. Taylor, Phys. Rev. B, 34 ( 1986) 2112.
[ 1J M.A. Herman and H. Sitter, in M.B. Panish (ed.), Molecular beam epitaxy-fundamentals and current status, SpringerSeries in Materials Science, Vol. 7, Springer, Berlin, 1988. [2] W. Faschinger and H. Sitter, /. Cryst. Growth, 99 ( 1990) 566. [3] N. Metropolis, A.W. Rosenbluth, N.M. Rosenbluth, A.H. Teller and E. Teller, J. Chem. Phys., 21 (1953) 1087. [4] H. Mtiller-Kmmbhaar. in K. Binder (ed.), Monte Carlo Methodr in Statistical Physics, Springer, Berlin, 1979. [5] R.H. Swendsen, Phys. Rev. B, 1.5 (1976) 5421. [6] D. Kashchiev, J.P. van der Eerden and C. van Leeuwen, J. Cryst. Growth, 40 (1977) 47. [7] R.H. Swendsen, P.J. Kortman, D.P. Landauand H. Miiller-Krumbhaar. J. Cryst. Growth, 35 (1977) 73. [ 81 C.L. Carson and J. Bemholc, Appl. Phys. Lat.. 56 ( 1990) 1971.
[ll I B.W. Dodson and P.A. Taylor, Phys. Rev. B, 36 (1987) 1355. [ 121 B.W. Dodson, Phys. Rev. E, 30 (1984) 3545.
[ 131 M. Schneider, I.K. Schuller and A. Rahman, Phys. Rev. B, 36 ( 1987) 1340.
[ 141 W.M. Plotz, K. Hingerl and H. Sitter, J. Cryst. Growth, I15 (1991) 186. [ 151 W.M. Plotz, K. Hingerl and H. Sitter, Phys. Rev. B, 45 ( 1992) 12122. [16] W.M. Plotz, K. Hingerl and H. Sitter, Semicond. Sci. Technol., 9 ( 1994) 2224. [ 171 M. Tagwerker, W.M. Plotzand H. Sitter, J. Cryst. Growth, 146 ( 1995) 220. [ 181 F.A. Ponce, G.B. Anderson and J.M. Ballinger, Surf: Sci., 168 (1986) 564.
[ 191 H. Sitter, W. Faschinger, Thin Solid Films, 225 ( 1993) 250.