Equations of motion for epitaxial growth

Surface Science Letters 274 (1992) L529-L534 North-Holland

surface science letters

Surface Science Letters

Equations of motion for epitaxial growth A . Z a n g w i l l a, C.N. L u s e a, D . D . V v e d e n s k y b a n d M . R . W i l b y a School of Physics, Georgia Institute of Technology, Atlanta, GA 30332, USA b The Blackett Laboratory, Imperial College, London SW7 2BZ, UK

Received 13 February 1992; accepted for publication 24 April 1992

A set of discrete, stochastic equations of motion which describe the epitaxial growth of a single crystal are derived beginning with a master equation description of the dynamics of a solid-on-solid model. The final set of coupled Langevin equations takes explicit account of the elementary microscopic processes of deposition, surface diffusion and desorption. The first of these contributes shot noise to the system while the last two produce configuration-dependent noise correlations. Direct numerical integration of these equations provides a formal alternative to Monte Carlo simulation of the growth process. Other applications and extensions are outlined in brief.

C o m p u t e r simulation has played a central role in the development of our understanding of epitaxial growth kinetics. Twenty years ago, M o n t e Carlo simulations of the so-called solid-on-solid (SOS) model [1] provided essential data for the test of then current analytic predictions for the macroscopic growth rate in both the nucleation and coalescence [2] and step flow [3] regimes of growth. This same methodology and model are being used today [4] with considerable success to provide a microscopic interpretation for the origin and nature of the distinct oscillations observed in reflection high energy electron diffraction ( R H E E D ) data collected in situ during molecular b e a m epitaxy (MBE). For example, these simulations provide detailed quantitative support for the notion [5] that the decay in the amplitude of R H E E D oscillations arises from a loss of correlation between spatially separated regions of the surface associated with a progressive roughening of the growth front as deposition proceeds. This type of result, combined with other developments in the intervening years such as M o n t e Carlo simulations of more realistic models [6] and molecular dynamics simulations [7] helps explain the growing popularity of computer methods in the study of epitaxiaI growth. O f course, there is still a place for analytic work in the theory of this problem. As an interesting example, we cite the recent suggestion by several authors [8-11] that the roughening p h e n o m e n a n o t e d above can be characterized at sufficiently large length and time scales by some form of nonlinear, stochastic, partial differential equation for the time evolution of the surface profile. This type of continuum (macroscopic) description has the great virtue that it lends itself to investigation by powerful analytic tools. On the other hand, it is not clear that such a description is appropriate for epitaxial growth where strong lattice anisotropy induces layer growth. Moreover, the terms omitted or retained in the final equations of motion generally have b e e n advanced on the basis of synunetry or phenomenological arguments alone. T h e purpose of the present L e t t e r is to describe an approach to epitaxial growth problems which retains the microscopic features of the computer simulation models yet assumes a form which permits comparison with the continuum equations of motion models. T o illustrate our method, we adopt a simple 0039-6028/92/$05.00 © 1992 - Elsevier Science Publishers B.V. All rights reserved

A. ZangwiUet aL / Equations of motion for epitaxialgrowth ppropriate to MBE growth. Therein, a column is associated with each site of the surface lattice. Deposition of an atom onto a site increases the corresponding column height by one unit. Desorption of an atom decreases the column height by one unit. A surface diffusion event decreases the height of one column and increases the height of a neighboring column. No overhangs or vacancies are permitted. Our aim is to derive a Langevin equation for the time evolution of each column height and thus for the surface morphology as a whole. These equations are constructed so that the time dependence of the associated column height joint probability distribution is identical to that found from the master equation for this quantity which takes exact account of the SOS elementary processes listed above. Although the study of crystal growth by master equation and related methods is not new [12-15], all previous investigations to our knowledge focus mainly on the statistical average of various spatially averaged quantities, e.g., the growth rate. The specifics of the model are as follows. A column height variable hij is assigned to every site of the (100) face of a simple cubic lattice of lattice constant a. The subscripts i and j refer to the x and y coordinates of the column. Since vacancies and overhangs are forbidden, every configuration of the system is specified completely by the set {hi1, hi2 , h21. . . . } - H . The joint probability that the system exhibits a particular configuration H at time t is denoted P(It; t). According to the general theory of stochastic processes [16], the evolution of P from a specified initial configuration is governed by the master equation

OP(H; t) ~t

~'~W(H', H ) P ( H ' ; t) - E W ( H , H ' ) P ( H ; t), n' H'

(1)

where the quantity W(H, H') denotes the transition rate from a configuration H to a configuration H ' . We choose the microscopic rate processes which govern the evolution of the column heights to be identical to those used in recent simulation studies [4]. In particular, growth is initiated by random deposition of atoms onto the substrate at a rate of r -1 per site. Any surface atom can migrate by hopping to one of its four nearest neighbors with a probability per unit time of k o exp(-E/kBT). The exponential prefactor k 0 is an attempt frequency which may be chosen equal to a typical adatom vibrational frequency. The migration barrier E is a configuration dependent quantity and is written as the sum of two factors: a term E s associated with the substrate and a contribution E N from each lateral nearest-neighbor. Thus, with/3 = 1/kBT, the total hopping from site (ij) is given by

k o exp(-~Es)z(hiy

, hi+l,j, h i _ l j , hi,j+l, hid_l) ,

(2)

where

z(hij, hi+l,j, hi-l,j, hij+l, hi,j-l) =

exp{-/~EN[0(hi+lj-hij) + O(hi_l,j - h i j ) + O ( h i . j + l - h i j ) + O ( h i , j _ l - h i j ) ] }

(3)

and the unit step function O(x) is defined by

O(x)=

1, 0,

if x_> O; if x < 0 .

(4)

The hopping rates to site (ij) from any of the four nearest neighbor sites are given by similar formulae. These expressions represent our ansatz for the transition probabilities W which appear in eq. (1). Thus, introducing the notation ro 1= k o e x p ( - ~ E s) for the hopping rate of an isolated adatom (no lateral

A. Zangwillet aL / Equations of motion for epitaxialgrowth neighbors), the master equation including only deposition and surface migration takes the

UP(H; t)

Ot

1

4"to ]~-'[P(hi'+a'hi+l'j-a'l~;t)z(hij+a'hi+l'y-a'hi-l"j'hi'j+l'hi'j-1) ij

+e(hij + a,

hi,j+ 1 - a , I~, l)z(hij ÷ a, hi,j+ 1 - a , hi_x,j, hi+l,j, hi,j_l)

+P(hij-a, hi+zj + a, a;

t)z(hi+l, j + a, h i j - a , hi+2,j, hi+l,j-1, hi+l,j+ 1)

+P(hij-a, hi,j+1 + a, 1~; t)z(hi,j+l + a, h i j - a , hi+l,j+ 1, hi_ld+l, hi,j+2) -e(hij, h,+ 1,j, 1~; t)z(hi+ 1,j,

hij, hi+2,j, hi+l,j-1, hi+x,j+l)

-P(h,j, h,,j+l, 1~; t)z(h,,j+l,

hi+l,j+ 1, hi_l,j+ 1, hi,j+ z, hij )

-2P(H;

t)z(hij , hi+l,j, hi_l,j, hi,j+l, hi,j-l)]

1

+-~ ~. [P(hij-a , I~; t ) - P ( H ; tj

t)].

(5)

The notation /-I signifies all of the hij not explicitly written in the argument of P. The terms in eq. (5) proportional to z o i describe site-to-site hopping and the terms proportional to -r-1 correspond to the deposition process. The eight hopping terms arise as follows. According to eq. (1), those terms with a positive sign correspond to processes which leave the (ij) column with a height exactly equal to hij while those terms with a negative sign correspond to processes which leave the (ij) column with a height different from hiy. In the first two terms, this column begins too high (by one unit) and an atom hops off to one of two neighboring sites [17]. In the second two terms, this column begins too low (by one unit) and an atom hops on from a neighboring site. In the next two terms, an atom hops onto the site (ij) so that the final column height is hij + a while the final two (equal) terms describe hopping-off processes where the final column is hij - a. The overall factor of 1/4 is present because the total probability to hop away from any site is equally divided amongst its nearest neighbors in an isotropic system. To take account of the possible evaporation of atoms back into the gas phase, we define a desorption barrier in analogy with the migration barrier which includes a term E~ from the substrate and a contribution E~ from each lateral nearest neighbour. Thus, using primes to denote quantities associated with desorption, we obtain the following terms to be added to eq. (5):

1 ~_~[P(hij+a,I~;t)z,(hij+ ..

q'!

q

a hi+lj, hi-l,j'hi,j-l'hi,j+l)

- P ( H; t)z'(hij, hi+l,j, hi_Lj, hid+l , hi,j_l) ] .

(6)

At this point, one could compute an equation of motion for, say, the statistical average of any particular column height directly from the master equation:

O
Ehijae(H; t) H

at

(7)

Using eqs. (1), (5) and (6), the right hand side of eq. (7) is seen to equal the statistical average of a quite complicated quantity. Formal evaluation of this object generates the average of yet more complex expressions and one is led rapidly to an infinite hierarchy of coupled equations of motion for successively

A. Zangwill et al. / Equations of motion for epitaxial growth correlation functions. To generate a closed set of equations, it is typical [12-13,18] to invoke some sort of factorization or truncation scheme. In practice, the problem rarely is pursued beyond the level of a mean field approximation. Alternatively, one can attempt to solve directly for the joint probability distribution P(H; t). In the limit of large system size, this can be done by solution of an appropriate Fokker-Planck equation [16,19]: O

1 82 -~-~P(H; t ) = - 8hi----~[g(i))P( H; t)] + 2 8hijOhlm [ gti2)/mP( H; t ) ] . 8

(8)

In this expression, summation over repeated indices is implied and the first and second moments of the master equation transition rates are defined by K,'(]) = E (h;y - hij)W(H, H')

(9)

H'

and

KC2)ij;lm= ~., ( h;j - hiy)( h;m - htm)W( H, H'). It'

(10)

The computation of these quantities is straightforward given the explicit form of the master equation in eq. (5). Noting that D s = a2/4% is the definition of the surface diffusion constant [19], one finds a

a

a

2

K~])= 4r0 (42 + a Y ) ' X i J - 7 ; t i J + - r

(11)

and

j 2 (l~i + ,~lm)(~jmA2 ~il + ~ilAy~jm)] 2 K(2)ij,tm=Ds[6it6/m( A2 + Ay)t~ij_

+ --TA'i/+ r -7 3il~Jm'

(12)

where we have introduced the condensed notation

Aij =z(hij, hi+ld,

hi_l,y,

hid+l, hi,j_l)

(13)

(and a similar expression for A'ij) and defined two discrete second derivatives as 2 __ Axfq = fi+ 1.j - 2fij +fi-l,j, A2 fij =-fid+ l -- 2f/j + fi,y-1,

(14)

for any double-indexed quantity fly. The symbols 3. and 8jm in eq. (12) denote Kronecker delta functions. It turns out that the numerical solution of eq. (8) is very demanding computationally. But as is well known [16], the resulting probability distribution P(H; t) can be associated uniquely with a Langevin equation for each component of H(t). Ignoring certain subtleties which arise when the statistical fluctuations are very large, these equations take the form [19]

d -~-~hii(t) =K~)+~ij(t),

(15)

where the Gaussian random variables nij have zero mean
(16)

A. Zangwill et al. / Equatior~ of motion for epitaxial growth

and are correlated as dictated by the covariances <'rlij( t )'tllm( t' ) >= K}2)mS( t - t ' ) .

(17)

The set of stochastic equations (15) constitutes our central result. They, the noise correlations defined in eqs. (16) and (17) and the explicit forms for the transition moments in eqs. (11) and (12) completely determine the growth dynamics of the SOS model as defined. Solutions for these equations can be obtained by well-tested numerical methods [21], and as such, provide an exact alternative to direct Monte Carlo simulation of the original model. Unfortunately, the computational labor involved is still considerable due to the complexity of the noise covariance matrix. Only the deposition process contributes simple, diagonal shot noise to the problem. Desorption contributes to the diagonal elements as well, but with an amplitude which depends upon the instantaneous disposition of the column heights. Surface diffusion is most complex of all-producing configuration-dependent contributions to both the diagonal and off-diagonal elements of the covariance matrix. Despite the foregoing, there are important virtues which accrue from our reformulation of the SOS model into a stochastic differential-difference equation. Not least among these is a fresh perspective on possible approximation schemes to reduce the computational complexity noted above. For example, most of the complications arise from the attractive lateral interactions (parameterized above by the quantities E N and E~) which induce the layer growth mechanism characteristic of epitaxy. One simple possibility is to replace the quantities Aiy and A'ij by constant average values A and A' in eq. (12) while retaining their exact values in the Langevin equation itself. In that case the noise correlations take the form ( ~ l i j ( t ) ' r / l m ( t ) ,> =

a 2 -"r'- + - -"i" t~il~jm

An even simpler version of this approximation was adopted by Weeks and Gilmer [1] in their discussion of a phenomenological stochastic equation of motion proposed for the SOS model. There, only the desorption contribution to eq. (12) was retained with A' set equal to twice the desorption rate from a kink (doubly-coordinated) site. Indeed, it may be sufficient to retain only the deposition shot noise in eq. (17). Numerical experiments designed to test these approximations will be reported elsewhere. Our set of discrete Langevin equations also forms the starting point for the search for a continuum equation of motion for a surface evolving according to the SOS model. That is, a partial differential equation for a quantity h(x, t) defined as the average of hij over a domain small on the macroscopic scale yet which contains many lattice points. Each such domain is associated with a single value of the two-dimensional continuous variable x. Upon performing such a spatial coarse-grain average on each term in eq. (15), it is not difficult to see that the resulting Langevin equation is ah 8 t = aDsV2A

a a --A'~.,+ --~.+ ~/,

(19)

where A, A' and ~/now denote the continuous functions of x and t obtained by averaging hij, h'ij and ~/~j over a macro-domain. The symbol V 2 is the usual Laplacian operator in two dimensions. Of course, this equation only attains the desired form if A(x, t) and h'(x, t) somehow can be expressed explicitly as functions of h(x, t). To date, we have not succeeded completely in this task and it is not obvious to us that such a program is even possible if eq. (19) is to be generally valid. On the other hand, if one is concerned exclusively with the properties of the surface at large length and time scales, we have shown elsewhere [22] that (apart from the precise form of the noise covariances) a one-dimensional version of eq. (19) neglecting desorption leads to a continuum equation of motion proposed for this problem by Villain [8] and Lai and Das Sarma [9]. To investigate growth phenomena in the non-asymptotic regime, one is forced to retain the full structure of eq. (19). Some insight into how to proceed may be gleaned from the observation that this

A. Zangwill et al. / Equations of motion for epitaxia[ growth

lentical in form to one proposed long ago by Mullins [23] for the time evolution of an non-facetted solid if one makes the identification h ( x ) = A'(x) =

exp[tz(x)/kBT],

(20)

where /x(x) is the local chemical potential of the surface. We exploit this observation in a companion paper [24] to derive a generalization of the classic analysis of Burton et al. [25] of growth in the step flow regime which takes account of lateral interactions between diffusing species. Comparison of our results with other approaches to this problem lends support to eq. (20) and suggests a search for a rigorous demonstration. Finally, we note that the basic methodology outlined in this paper can be applied to more general descriptions of epitaxial growth than the simple SOS model adopted here. Anisotropic surface diffusion and step attachment kinetics are particularly simple to incorporate. Other obvious avenues for further development include the incorporation of multiple surface species and epitaxial strain. Work on these issues currently is in progress. The authors acknowledge valuable conversations with Ronald Fox. Work performed at Georgia Tech is supported by the US Department of Energy under grant No. DEFG05-88ER45369. Work performed at Imperial College is supported by Imperial College and the Research Development Corporation of Japan under the auspices of the "Atomic Arrangement: Design and Control for New Materials" Joint Research Program. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25]

For a review, see: J.D. Weeks and G.H. Gilmer, Adv. Chem. Phys. 40 (1979) 157. G.H. Gilmer and P. Bennema, J. Appl. Phys. 43 (1972) 1347. C. van Leeuwen, R. van Rosmalen and P. Bennema, Surf. Sci. 44 (1974) 213. D.D. Vvedensky, S. Clarke, K.J. Hugill, A.tC Myers-Beaghton and M.R. Wilby, in: Kinetics of Ordering and Growth at Surfaces, Ed. M.G. Lagally (Plenum, New York, 1990) pp. 297-311. M.G. Lagally and R. Kariotis, Appl. Phys. Lett. 55 (1989) 960. A. Madhukar and S.V. Ghaisas, CRC Crit. Rev. Solid State Mater. Sci. 13 (1987) 1434. B.W. Dodson, CRC Crit. Rev. Solid State Mater. Sci. 16 (1990) 115. J. Villain, J. Phys. (Paris) 1 (1991) 19. Z.-W. Lai and S. Das Sarma, Phys. Rev. Lett. 66 (1991) 2348. L.-H. Tan and T. Nattermann, Phys. Rev. Lett. 66 (1991) 2899. D. Wolf, Phys. Rev. Lett. 67 (1991) 1783. Y. Saito and H. Miiller-Krumbhaar, J. Chem. Phys. 70 (1979) 1078. W. Haubenreisser and H. Pfeiffer, in: Modern Theory of Crystal Growth I, Eds. A.A. Chernov and H. Miiller-Krumbhaar (Springer, Berlin, 1983) pp. 43-73. K. Wada, M. Kaburagi, T. Uchida and R. Kikuchi, J. Stat. Phys. 53 (1988) 1081. C. Garrod, A.C. Levi and M. Touzani, Solid State Commun. 75 (1990) 375. C.W. Gardiner, Handbook of Stochastic Methods (Springer, Berlin, 1985); N.G. van Kampen, Stochastic Processes in Physics and Chemistry (North-Holland, Amsterdam, 1981). Here, and below, the other two sites available for the process in question are counted when the summations over i and j in eq. (4) are performed. A particularly clear exposition of this approach in the context of a two-dimensional kinetic lattice gas can be found in: H.J. Kreuzer and J. Zhang, Appl. Phys. A 51 (1990) 183. R.F. Fox and J. Keizer, Phys. Rev. A 43 (1991) 1709. See, e.g., A. Zangwill, Physics at Surfaces (Cambridge University Press, Cambridge, 1988) Ch. 14. R.F. Fox, J. Stat. Phys. 54 (1989) 1353. A. Zangwill, C.N. Luse, D.D. Vvedensky and M.R. Wilby, in: Interface Dynamics and Growth, Eds. K.S. Liang, M.P. Anderson, R.F. Bruinsma and G. Scoles (Materials Research Society, Pittsburgh, 1992) pp. 189-198. W.W. Mullins, J. Appl. Phys. 28 (1957) 333. C.N. Luse, A. Zangwill, D.D. Vvedensky and M.R. Wilby, Surf. Sci. Lett. 274 (1992) 1_535. W.K. Burton, N. Cabrera and F.C. Frank, Philos. Trans. R. Soc. A 243 (1951) 299.