Multilayer Scaling and Universal Behavior of Molecular Beam Epitaxy

width and reflection high energy electron diffraction (RHEED) intensity have been determined as ... characterization techniques, like reflection high-...
0 downloads 0 Views 598KB Size
Download PDF
Langmuir 1996, 12, 29-35

29

Multilayer Scaling and Universal Behavior of Molecular Beam Epitaxy Grown Surfaces† D. P. Landau* and S. Pal Center for Simulational Physics, The University of Georgia, Athens, Georgia 30602 Received September 1, 1994. In Final Form: December 12, 1994X Extensive kinetic Monte Carlo simulations have been used to model the growth of epitaxial films by molecular beam epitaxy (MBE). We allow both deposition and diffusion (single particle hopping events) to occur simultaneously and have investigated several different rules for hopping. Simple solid-on-solid models have been used with L × L substrate sizes and with nearest neighbor interactions only. Interfacial width and reflection high energy electron diffraction (RHEED) intensity have been determined as measures of the roughness of the surface. We find that in the long time limit the interfacial width has a logarithmic dependence on time and substrate size, although there is a pronounced crossover region for moderate times which depends on both temperature and flux. These results suggest that these models fall in the universality class of the Edwards-Wilkinson model. Comparison with experimental results shows similar behavior up to intermediate times and suggests that experimental data are needed over longer time scales.

1. Introduction Attempts to understand crystal growth, and the attendant nucleation processes which play a role, have been actively pursued for most of this century. Much of our understanding of microscopic crystal growth modes has come about in recent years as powerful experimental and analytical tools have become available due to rapid progress in technology. Additional progress has been made by the recent appearance (by historical standards) of a now widely used method that can grow films with desired specifications on the atomic level; molecular beam epitaxy (MBE). Through accurate control of the growth environment MBE enhances our ability to test theories of crystal growth with the help of modern surface characterization techniques, like reflection high-energy electron diffraction (RHEED) and scanning tunneling microscopy (STM). The precise control over growth parameters and the actual determination of observables may still be rather difficult in experiments, and often the results may be skewed by irrelevant parameters over which the experimenter may have limited control. The growth of a surface is a “nonequilibrium” process in which particle deposition competes with particle diffusion. The incident flux and the temperature together control the deviation of the surface from the equilibrium structure. The surface dynamics at equilibrium have been studied using theoretical techniques and Monte Carlo simulations, and the various issues are well understood. (For reviews of equilibrium and nonequilibrium studies see refs 1 and 2.) It is now well-known that under equilibrium conditions the surface fluctuations are finite † Presented at the symposium on Advances in the Measurement and Modeling of Surface Phenomena, San Luis, Argentina, Aug 24-30, 1994. X Abstract published in Advance ACS Abstracts, January 1, 1996.

(1) Kawasaki, K. In Phase Transitions and Critical Phenomena; Domb, C., Green, M. S., Eds.; Academic Press: New York, 1983; Vol. 2, p 443. Binder, K. In Phase Transitions and Critical Phenomena; Domb, C., Lebowitz, J. L., Eds.; Academic Press: New York, 1983; Vol. 8, p 1. Gunton, J. D.; San Miguel, S.; Sahni P. S. In Phase Transitions and Critical Phenomena; Domb, C., Lebowitz, J. L., Eds.; Academic Press: New York, 1983; Vol. 8, p 269. Landau, D. P. In Applications of the Monte Carlo method in Statistical Physics; Binder, K., Ed.; Springer Verlag: New York, 1987; p 93. Structure and Dynamics of Surfaces; Schommers, W., Von Blanckenhagen, P., Eds.; Springer Verlag: New York, 1986; Vol. 1. (2) Solids Far From Equilibrium; Godre`che, C., Ed.; Cambridge University Press: New York, 1992. Dynamics of Fractal Surfaces; Family, F., Vicsek, T., Eds.; World Scientific: Singapore, 1991. Family, F. Physica 1990, A168, 561.

0743-7463/96/2412-0029$12.00/0

below a temperature TR, known as the roughening transition temperature, above which the fluctuations diverge. The study of surface growth provides an opportunity to understand how the surface fluctuations are modified under nonequilibrium conditions. Various theoretical approaches have been pursued in an attempt to provide a framework for the understanding of MBE growth. These have included a variety of simple models as well as stochastic Arrhenius models, i.e., which include an activation barrier for diffusion, which employ more realistic surface dynamics. Extensive references to earlier simulational work on these models are given in refs 3 and 4. The main hurdle to the complete solution of stochastic Arrhenius models is the requirement of enormous computing resources, so most of the earlier simulations for the study of the macroscopic behavior were done on one-dimensional interfaces.5 Very few results are available for these models in 2 + 1 dimensions for the long time behavior.3 The problem of determining the asymptotic behavior has also been tackled from a different perspective using numerical solution of continuum growth equations.6,7 Note that the growth equations themselves have to be derived from purely theoretical considerations. The comparison of the simulational results of computer models with those of the growth equations may serve to shed light on what the essential terms are in any theoretically derived growth equation. In our work, we have simulated simple 2 + 1 dimensional models, with coarse grained “effective” interactions, depositing a large number of particles to produce behavior which occurs on mesoscopic length scales. This approach allows us essentially complete control over the relevant parameters of the problem, thus permitting a straightforward interpretation of the results. In the absence of first-principle calculations (which are by far too difficult to complete with current computational resources), one can carry out classical simulations using parametrized potentials in the hope that the essential physics will still be captured. In principle, one can then perform molecular dynamics studies to follow the time development of these (3) Pal, S.; Landau, D. P. Phys. Rev. 1994, B49, 10597. (4) Metiu, H.; Lu, Y. T.; Zhang, Z. Science 1992, 255, 1088. (5) Tamborenea, P. I.; Das Sarma, S. Phys. Rev. 1993, E48, 2575. (6) Amar, J. G.; Family, F. Phys. Rev. 1990, A41, 3399. (7) Grossmann, B.; Guo, H.; Grant, M. Phys. Rev. 1991, A43, 1727.

© 1996 American Chemical Society

30

Langmuir, Vol. 12, No. 1, 1996

systems, but here too, we are many orders of magnitude away from being able to achieve macroscopic time scales. In this article we are interested in pursuing the development of structures for relatively long times and large length scales using kinetic Monte Carlo methods. Such simulations consider lattice models and describe diffusion in terms of simple particle hopping from site to site. This approach also uses very simple, near neighbor interactions which govern the hopping. Our expectation is that the resultant structures will not depend upon many of the details of the models and that the dynamic behavior will fall into one of a small number of universality classes. The models may be also used to study a coarse grained microstructure on the surface, which in turn may help in understanding mechanisms that lead to such a structure. The next section briefly reviews that relevant theory, prior simulational work, and the solution of diverse growth equations. Our simulational models and methods are described in section 3 while section 4 presents the results and analysis. We draw comparisons to experimental results with our simulations and discuss the factors affecting the roughness of a growing surface and then add a few concluding words in section 5. 2. Background Growth by MBE is carried out in an ultrahigh vacuum chamber to prevent contamination and to preserve a beam of mass flow toward the substrate. A substrate of a particular material with a given orientation is maintained at some temperature Ts, while a gas of impinging atoms falls on the surface. When the gas atoms are the same as that of the substrate, the growth is homoepitaxial, and when the species differ it is called heteroepitaxial. The competition between the incident flux and the surface diffusion controls the surface structure, which can be tuned to achieve a particular growth morphology. In experiments the roughness of the growing surface is measured by using reflection high-energy electron diffraction (RHEED), where the intensity is proportional to the smoothness of the surface. It has been observed that for a system subject to constant flux and held at fixed, low temperature, the RHEED intensity smoothly decreases with time, indicative of rough growth. As the temperature is elevated, the intensity starts to oscillate with time, which is indicative of layer-by-layer growth. When deposition is carried out on metal substrates, the surface may show oscillations at low temperatures (depending on the flux), which are found to decay when the temperature is raised,8 only to return when the temperature is raised still further; this anomalous behavior has been attributed to the existence of barriers to atoms hopping down a step edge.9 It may be sometimes quite difficult to grow a crystal in the layer-by-layer growth mode as the growing surface naturally tends to become rougher as more layers are grown. Apart from tuning the temperature and the flux, one can add surfactants10 or bombard the surface with high-energy ions to improve epitaxy of the growing surface.11 The adoption of such methods alters the true interactions and structure of the surface, which in turn complicates the understanding of the mechanisms of film growth. In the initial stages of the growth, small clusters (8) Kunkel, R.; Poelsema, B.; Verheij, L. K.; Comsa, G. Phys. Rev. Lett. 1990, 65, 733. (9) Ehrlich, G.; Hudda, F. J. Chem. Phys. 1966, 44, 1039. Schwo¨bel, R. L. J. Appl. Phys. 1969, 44, 614. (10) Li, D.; Freitag, M.; Pearson, J.; Qiu, Z. Q.; Bader, S. D. To appear in J. Appl. Phys. van der Vegt, H. A.; van Pinxteren, H. M.; Lohmeier, M.; Vlieg, E.; Thornton, J. M. C. Phys. Rev. Lett. 1992, 68, 3335. Copel, M.; Reuter, M. C.; Horn von Hoegen, M.; Tromp, R. M. Phys. Rev. 1990, B42, 11682. (11) Stroud, P. T. Thin Solid Films 1972b, 11, 1.

Landau and Pal

are formed from the individual adatoms on the surface. As time progresses, these clusters may grow, nucleate and finally coalesce to form a continuous film. At any instant of time, the flux increases the number of nucleation centers and surface diffusion reduces them by increasing the size and merging of clusters. In the layer-by-layer growth regime, nucleation and growth occur simultaneously and most often nucleation processes already start occurring on top of a layer before it is complete. Thus nucleation at higher layers is in competition with the growth of the bottom incomplete layers, and this competition plays an important role in determining the resultant surface structure. There is also a theoretical framework for the understanding of MBE growth. The deposition of atoms and their movement on the surface, along with defect creation, result in fluctuations of the surface structure. The root mean square of these fluctuations at any instant of time is a measure of the interfacial width W(L,t). The fluctuations reach a temporally and spatially scale invariant state when growth is allowed to occur for very long times and satisfy a scaling relation12

W(L,t) ) LRf

( ) t LR/β

(1)

where L is the lateral extent of the substrate and the growth exponents R and β define the surface fluctuations of the growth model and characterize the universal behavior. The ratio z ) R/β is called the dynamic exponent, and the scaling function f(x) obeys the relation f(x) ∼ xβ, for x , 1, and f(x) f constant, for x . 1. The exponents R and β lie in the range of 0.08 < R < 0.4 and 0.04 < β < 0.25 for the Eden, ballistic deposition13 and the restricted solid-on-solid model.14 In order to understand the physical processes responsible for the values of the exponents, one has to compare them with those obtained from a growth equation. Instead of writing multiple growth equations separately, we write a comprehensive equation15

∂H/∂t ) ν∇2H + λ(∇H)2 + K∇2(∇2H) + σ∇2(∇H)2 + R + ξ (2) where R is the renormalized deposition flux, ξ is the zero mean random fluctuation in the flux, and H is the height, measured relative to the average height of the surface. When the coefficients λ ) K ) σ ) 0, the equation reduces to the Edwards-Wilkinson (EW) equation.16 Normally ν is regarded as the surface tension, but in the EW formulation the coefficient ν is proportional to the flux F, which in turn is proportional to the interface velocity. The EW equation appears to be symmetric about the interface, thus giving a notion that it corresponds to equilibrium dynamics. Note that the equilibrium derivation of surface current assumes that the condition of detailed balance is satisfied, which, however, does not hold in the presence of flux (R + ξ).17 Thus any growth equation cannot be related to an equilibrium formulation. When the coefficients K ) σ ) 0, we get the equation proposed by Kardar-Parisi-Zhang (KPZ),19 where ν is the surface tension, while λ is proportional to the incident flux. Mapping the KPZ equation to the directed polymer problem allows it to be solved in the large coupling limit (12) Family, F.; Viscek, T. J. Phys. 1985, A18, L75. (13) Krug, J.; Spohn, H. In Solids Far From Equilibrium; Godre`che, C., Ed.; Cambridge University Press: New York, 1992; p 532. (14) Kim, J. M.; Kosterlitz, J. M. Phys. Rev. Lett. 1989, 62, 2289. (15) Villain, J. J. Phys. (Paris) 1991, I 19, 1. (16) Edwards, S. F.; Wilkinson, D. R. Proc. R. Soc. London 1982, A381, 17. (17) Hohenberg, P. C.; Halperin, B. J. Rev. Mod. Phys. 1977, 49, 435.

Growth of Epitaxial Films

Langmuir, Vol. 12, No. 1, 1996 31

(i.e., large λ). There have been a number of numerical studies of the KPZ equation and models which are believed to be equivalent, and there is a diversity of opinion about what the actual exponent values are.13,18 (The original estimates were that R ) 0.53 ( 0.08, β ) 0.33 ( 0.03, and z ) 1.6 ( 0.119 in 2 + 1 dimensions; more recent studies yield different but not mutually consistent values.13) The universality class of Wolf-Villain and Das SarmaTamborenea (WV-DT)20 is obtained when ν ) λ ) σ ) 0. In this case R ) (d - 5)/2, β ) (d - 5)/8, thus z ) 4.21 Nozie`res and Gallet (NG)22 have analyzed a theory of crystal growth with a pinning potential (implying that growth is carried out below the roughening transition temperature) using

1 ∂H ) γ˜ ∇2H + φ + R + ξ Γ ∂t

(3)

where Γ is the mobility, γ˜ is the stiffness tensor, and φ is the pinning force. For large fluxes the theory concludes that the asymptotic behavior of their equation is the same as that of the equation of Edwards and Wilkinson, which predicts a logarithmic time and substrate size dependence of the surface fluctuations. It is worth comparing the different times in the EW, NG, and the KPZ equation. The forms of the NG and the EW equations are essentially the same, except a pinning force is included in the NG equation. In the absence of flux, the NG equation describes the equilibrium situation, while the EW equation becomes invalid (as flux is the coefficient of the ∇2h term). Above the equilibrium roughening transition temperature, the NG equation has a dynamic exponent z ) 2. The dynamic exponent can be obtained using the standard “equilibrium” methods of the dynamic critical phenomenon.23 In the presence of flux, the EW and NG equations describe nonequilibrium growth equations and thus it is doubtful whether any equilibrium methods will be successful to obtain the dynamic exponent. To overcome the problems in the current theoretical understanding of non-equilibrium systems, Kardar, Parisi, and Zhang,19 theorized that the addition of nonlinear terms to the linearized equilibrium equation can account for the nonequilibrium effects. In this way the growth equation looks like an expansion around the equilibrium fluctuations,13,24 where the KPZ term is the first-order correction to the nonequilibrium effects. This important deduction can be viewed as a way to solve nonequilibrium equations by expanding around the equilibrium fluctuations, such that the current “equilibrium” dynamical renormalization group methods can be used to solve the problem. These methods do not work well for the KPZ equation in 2 + 1 dimensions, and thus the equation has only been solved numerically. Large scale simulations by Pal and Landau3 found EW universality in a stochastic Arrhenius model of crystal growth with z ) 1.61 ( 0.02, which, perhaps surprisingly, is similar to several of the estimates for the KPZ equation.13,19 (It has been suggested, however, that z should equal 2 for the EW equation in 2 + 1 dimensions.13) The similarity of the dynamic exponent of a nonequilibrium model of crystal growth with that of a nonlinear stochastic equation would appear to confirm the validity of the expansions of the growth equation around equi(18) See for example, Meakin, P. Phys. Rep. 1993, 235, 189. (19) Kardar, M.; Parisi, G.; Zhang, Y. Phys. Rev. Lett. 1986, 56, 889. (20) Wolf, D. E.; Villain, J. Europhys. Lett. 1990, 13, 389. Das Sarma, S.; Tamborenea, P. Phys. Rev. Lett. 1991, 66, 325. (21) Lai, Z. W.; Das Sarma, S. Phys. Rev. Lett. 1991, 66, 2348. (22) Nozie`res, P.; Gallet, F. J. Phys. (Paris) 1987, 48, 353. (23) Forster, D.; Nelson, D. R.; Stephen, M. J. Phys. Rev. 1977, A16, 732. (24) Rost, M.; Spohn, H. Phys. Rev. 1994, E49, 3709.

librium to take into account nonequilibrium effects. However, inconsistencies still exist in eq 2 regarding whether the asymptotic behavior is dominated by the KPZ term, i.e., (∇H)2, or the EW term (∇2H). Although it has been hypothesized that the long-time behavior of a growing surface should be controlled by the KPZ term,13,19 our simulations indicate that in 2 + 1 dimensions the asymptotic behavior is dominated by the EW term which converts the algebraic dependence on time and substrate sizes of the interface width to a logarithmic one. It also seems possible, then, that previous numerical integrations6,7 of the KPZ equation may not actually have reached the asymptotic limit (we will demonstrate the difficulties which one may encounter in reaching the asymptotic limit using our method in the results section), and thus it is not clear whether the exponents R and β were really reliable. It is thus possible that the EW and KPZ universality are identical, with the ∇2H term controlling the asymptotic behavior. This is not surprising, as on a macroscopic scale the ∇2H term arises due to the chemical potential difference between the solid and the vapor phase created by the incident flux. 3. Model and Method We consider a solid-on-solid model growing on an L × L square lattice substrate with periodic boundary conditions. The height of the surface above the two-dimensional substrate at position r and time t is denoted by h(r,t). Deposition is carried out randomly on the surface and each deposition event increases the height at a site by unity. We assume that the deposition of L2 particles per Monte Carlo time unit is equivalent to the growth of one layer per second for a system with a lattice constant of a ) 10 Å and a flux of 1014/cm2 s. This choice allows us to compare the time scales in the experiments with those of the simulations. A small number of particles are deposited randomly on the surface, in a small time step, and then a fraction of particles on the surface is allowed to attempt to hop to nearest neighbor sites. The size of the time step δt depends on the flux (R), and varies as δt ) 0.02/Ra2, while the number of particles deposited at a time step is R(La)2δt. All particles are treated identically when it comes to hopping, and no attempt is made to dissipate the energy of incident particles. Hopping attempts are really two-step procedures. First, a particle attempts to overcome an activation barrier, and then, if it succeeds, it must choose a neighboring site to which to move. The thermally activated hopping probability is given by the equation Phop ) exp(-EA/kbT), where EA ) nJ is a site dependent activation energy, J is the “effective” bond energy, and n is the number of nearest neighbors. Here J is an effective bond energy that results dynamically from correlated motions of the neighbors. In this study we confine ourselves to the case of homoepitaxy, so all nearest-neighbors are equivalent. We studied several different models with different hopping rules for atoms which have to overcome the activation barrier. In the DP (down-preferential) model a particle which undergoes diffusion moves may hop to a nearest neighbor column subject to the restriction that the final position is not higher than the initial one and with a probability which depends on the number of bonds that will be formed in the new site. In the DR model (down-random) the particle moves to a randomly chosen nearest neighbor site once it has overcome the activation energy. In the UDP (up-downpreferential) and UDR (up-down-random) models there is no restriction on the final height of the adatom and the motion is preferential or random, respectively. The question of the nature of the hopping, particularly at step edges has been the topic of some discussion for considerable

32

Langmuir, Vol. 12, No. 1, 1996

Landau and Pal

time. From the purely theoretical perspective, it is also interesting to see if all four models belong to the same universality class. In the case of preferential hopping, the unnormalized probability of a particle hopping to the k th site is proportional to

Pk ) exp(-Ek/kBT)

(4)

where Ek is the site energy available at the kth site. The total probability of moving to any of the nearest neighbors is denoted ∑P. A random number “rnd” is generated such that, 0 < rnd < ∑P, and the hop is carried out depending on which “interval” in the range of probabilities the random number lies. This algorithm ensures that the adatom prefers to move to the site where it will find the highest coordination. We emphasize again that we do not expect that these simple models will describe the behavior of physical systems in a quantitative manner but hope rather that they will nonetheless yield the correct functional form for the behavior of MBE grown films as well as the same growth exponents. Simulations were carried out for L ranging from 10 to 1000 for a wide range of temperatures and fluxes which varied from 1010/cm2 s to 1015/cm2 s using IBM RS/6000 workstations. As many as 108 particles were deposited in a single run and runs were repeated many times with different random number sequences with the data then averaged over multiple runs. To determine the roughness of the growing surface, we measured the interface width W(L,t) from the configurations produced by the simulation. The width is defined as

W(L,t) ) [〈h2(r,t)〉 - 〈h(r,t)〉2]1/2

(5)

where the thermal averages 〈hn〉 are calculated using 〈hn(r,t)〉 ) L-2∑rhn(r,t), and the overbar denotes the average over the statistically independent runs. The RHEED intensity is measured in the out-of-phase condition using the first top five layers above the completely filled layer, and is given by 5

4

3

2

c2i - 2∑cici+1 + 2∑cici+2 - 2∑cici+3 + ∑ i)1 i)1 i)1 i)1

I ) G0(

4

2

cici+1 + 4G3∑cici+3 ∑ i)1 i)1

2c1c5) + 4G1

(6)

where ci is the coverage of the ith layer above the highest fully filled layer. G0, G1, and G3 are functions that need to be adjusted in order to improve the sensitivity of the RHEED intensity due to step edges. Following ref 25 for our measurements, we chose G1 ) G3 ) 0.025G0. Visual images of the surfaces are generated using a Silicon Graphics (Indigo) workstation to aid intuitive understanding of the growth characteristics. 4. Results With the deposition of just a few layers of material, the simulation data are qualitatively quite close to what is observed experimentally. At quite low temperature, diffusion is suppressed and growth is essentially by ballistic deposition; hence the width increases smoothly with time and the RHEED intensity shows only subtle structure. As the temperature increases, oscillations appear, in a correlated way, in both the width and RHEED intensity, and those continue to grow as the temperature (25) Lent, C. S.; Cohen, P. I. Surf. Sci. 1984, 39, 121.

Figure 1. Evolution of interfacial width and the RHEED intensity shown for the (DP) model for relatively early times. Note that the number of layers deposited is proportional to the deposition time. The arrows in the figure indicate the time when the flux is turned off. The RHEED intensities have been shifted vertically to clarify portions of the graph. A lattice size of L ) 800, and flux ) 1014/cm2 s is chosen for a range of kBT/J.

is raised. Sample results are shown in Figure 1 for the DP model on quite large substrates with a flux of 1014/cm2 s. Since the RHEED data curves all lie quite close to one another, we have shifted each curve upward in order to separate them. The oscillations suggest layer-by-layer growth, but the slight decay in the amplitudes suggests that successive layers are becoming slightly less complete. We can see that in the temperature region where the RHEED oscillations are pronounced, the interfacial width does not grow dramatically and, in fact, seems to oscillate with a period which is commensurate with the time it takes to complete a layer. The arrows in Figure 1 show where the flux is turned off; beyond this point the diffusion attempts to make the surface smooth. Other simulational work by Vvedensky et al.26 has shown that in some cases this recovery of RHEED intensity can be made quantitative by the inclusion of appropriate step edge barriers in the hopping rules. The data for a flux of 1012/cm2 s, shown in Figure 2, indicate that there are quite noticeable differences in the behavior as compared with results at the higher flux. Here, as the temperature increases the interfacial width (at constant temperature) first decreases and then begins to increase again. In addition, the oscillations which appear show a pronounced tendency to increase in value with increasing time. The RHEED curves are now well separated so there is no shifting of curves for successive temperatures. The oscillations in I become quite large, but there is also a clear tendency for (26) Vvedensky, D. D.; Clarke, S. Surf. Sci. 1990, 225, 373.

Growth of Epitaxial Films

Langmuir, Vol. 12, No. 1, 1996 33

Figure 4. Evolution of the surface width plotted for the (UDP) model, up to intermediate times. Two different temperatures and fluxes are studied for L ) 160. The solid lines are shown for comparison, and the slopes are indicated by the β values.

Figure 2. Time evolution of the surface width and the RHEED intensity shown for the (UDP) model. Lattice size of L ) 160, and flux ) 1012/cm2 s is chosen for different values of kBT/J. Note that in this case the RHEED intensities have not been displaced vertically, to emphasize the crossing of the intensities as the temperature is raised.

Figure 3. Interfacial width plotted with the number of deposited layers, for the (UDP) and the (DP) model. The flux is kept fixed at 1012/cm2 s, and a lattice of size L ) 160 is chosen for the simulation. Two different temperatures are used to show the differences in the time evolution of the models.

their amplitude to decrease with increasing layer number. This strongly suggests that successive layers are clearly increasingly imperfect. In Figure 3 we compare the behavior of the interfacial width for the two models with preferential hopping, UDP and DP, at two different temperatures. At the lower temperature, kT/J ) 0.3, W oscillates but grows fairly

rapidly and is essentially identical for both models. At the higher temperature, however, after less than a monolayer has been deposited, the width increases substantially faster for the UDP model than for the DP model. We interpret this behavior in the following way. At low temperatures it is unlikely that particles at a step edge will want to hop up to the next layer, where they are unlikely to have lateral nearest neighbors to attach to, even if they are permitted to. At higher temperatures diffusion is enhanced so more particles migrate to the step edges where they can drop down to help complete the lower layer; the result is better layer-by-layer-like growth. Because of the higher temperature, however, upward hops become more likely, so the UDP model develops rougher surfaces. The long-time behavior of these models is even more interesting. We believe that the UDP model is probably the most physical so we show only data for this case in Figure 4; the surface width at “long times” seems to grow with a power law whose exponent depends on both flux and temperature. The nature of the surface also changes dramatically as time progresses. The initial impression is that there is no universality and that different exponents result from different growth conditions. In contrast, our earlier data for the DP model showed rather convincingly that the surface width grew logarithmically with time. The evolution of the surface structure for the UDP model growing at a low flux is shown in Figure 5 where we see “snapshots” of the surface for four different times. For short times the morphology of the surface is one with extensive fine detail. The characteristic size of the structure continues to increase with time, and after 600 layers have been grown, the surface shows rather coarse bumps. (The height of these bumps is also much greater than that which is observed for the DP model.) When the long time behavior is carefully examined, however, we see that in those cases where we can grow a sufficient number of layers, the power associated with the increase in the width slowly decreases. In fact, when plotted on a semilogarithmic scale, the asymptotic behavior becomes linear signifying a logarithmic growth law. Figure 6a shows the behavior for the DP model which very quickly becomes logarithmic. Data have been obtained for a number of different lattice sizes and temperatures for a flux of 1014/cm2 s and finite size scaling of the saturated width has been tested. We find that scaling is obtained 2 ) W20 + m exp(3.8J/ with the width obeying the form Wsat

34

Langmuir, Vol. 12, No. 1, 1996

Landau and Pal

a

b

Figure 5. Evolution of the surface structure shown for the (UDP) model after deposition of (top to bottom) 0.25, 1.0, 20.0, and 600.0 monolayers. The parameters used are L ) 160, flux ) 1012/cm2 s, and kBT/J ) 0.6. The lighter shades indicate an area of greater height above the substrate.

kBT) ln(L), where W0 ) 0.50 ( 0.01. The resultant scaling plot, for a wide range of temperatures, is shown in Figure 6b. In Figure 7a we show data for the UDP model for a very high flux of 1015/cm2 s. Also plotted in this figure are some recent experimental results27 for the vapor deposition of Ag films on quartz surfaces at room temperature and numerical integration data6 for the KPZ equation. Pronounced finite size effects are obvious for long times in the simulations, but for the largest lattice, over a decade and a half, the agreement in the slope with the experiment and numerical solution of KPZ equation with (λ ) 316.2) is quite close. With increasing time, however, the results from the simulation slowly bend over by an amount which cannot be explained by any conventional finite size rounding. If we examine the L ) 40 data on a semilogarithmic scale, shown in Figure 7b, we can see that the asymptotic behavior is indeed logarithmic but that there is a long crossover region in time where the data are consistent with an “effective” power law. It thus appears likely that the exponent which can be deduced from existing experimental data or the numerical solution of the KPZ may not really describe the asymptotic regime. Possible phase transitions and crossover effects for the KPZ model have been discussed by Nattermann and Tang.28 There is, of course, no reason to expect that our model describes the size of the crossover region in the physical system, so it is unclear just how thick the films will have to be grown in the laboratory before the true long-time regime is encountered. Similar effects were seen with a smaller flux, and there the asymptotic regime occurred sooner. This was not altogether helpful, however, since with a lower flux it takes much longer to deposit a given number of atoms. A full dynamic scaling analysis can best be made by first defining a “scaled interfacial width” W ˜ as W ˜ ) W(L,kBT/J,t)/W(L,kBT/J,tf∞), where the time t is proportional to the number of deposited layers. (27) Palasantzas, G.; Krim, J. Mater. Res. Soc. Symp. Vol. 317, p 111. (28) Nattermann, T.; Tang, L.-H. Phys. Rev. 1992, A45, 7156.

Figure 6. (a) Long-time data of the interfacial width shown for the (DP) model. A flux of 5 × 1012/cm2 s, and kBT/J ) 1.2 is used for a range of lattice sizes. The approximate location of the saturation time τc is shown by the arrow for L ) 80. For the smallest lattice size L ) 20, a fit is made through the data points to show the oscillations in the interfacial width. (b) Square of the saturated interface width plotted for a range of kBT/J, using the flux of 1014/cm2 s. The x-axis is a variable that depends on both the lattice size and kBT/J values. The straight line fit through the data points has a slope of m ) (5.64 ( 0.05) × 10-3. Thus the saturated interface width satisfies a relation 2 Wsat ) W20 + m exp(3.8J/kBT) ln(L), where W0 ) 0.50 ( 0.01.

The scaled time ˜t is defined as ˜t ) 2t(kBT/J)φ/Lz, where φ ) 3.85, and the choice of z is used to obtain scaling. In Figure 8 we test full dynamic finite size scaling for the DP model for which extensive data exist. The scaled interfacial width has been plotted against the scaled time and temperature. When plotted in this way, the data for a wide range of temperatures, times, and lattice sizes collapse onto a single curve, thus verifying the scaling ansatz. In Figure 9 dynamic finite size scaling is also tested for a flux of 1015/cm2 s, on the UDP model. These data scale over a substantial range of scaled time only for the largest lattices, indicating that corrections to scaling are quite important for the two smaller lattice sizes. Given the small values of L used, it is not altogether surprising that this is the case. 5. Conclusions We have shown that kinetic Monte Carlo methods are now able to simulate the growth of simple MBE model films which begin to achieve mesoscopic size. With the use of several simple solid-on-solid models we have

Growth of Epitaxial Films

Langmuir, Vol. 12, No. 1, 1996 35

a

b Figure 8. Dynamic scaling plot for flux ) 1014/cm2 s shown for different lattice sizes and temperatures for the (DP) model. The scaled width W ˜ is defined as W ˜ ) W(L,kBT/J,t)/W(L,kBT/ J,tf∞), where the time t is proportional to the number of deposited layers. The scaled time ˜t is defined as ˜t ) 2t(kBT/ J)φ/Lz, where φ ) 3.85, and z ) 1.61 is used to obtain the scaling. The data are consistent with a logarithmic scaling function.

Figure 7. (a) Long-time data for the (UDP) model shown for a range of lattice sizes in a logarithmic scale. A flux of 1015/cm2 s, and kBT/J ) 0.6 is used for the simulations. The filled circles indicate experimental results.27 The data points with the curve fits show numerical simulation results6 for the KPZ equation in the same time scale as our simulations. Note that in 1 s 10 layers are deposited when flux ) 1015/cm2 s. (ν ) 1) is used in eq 2, with λ ) 316.2 (top) and λ ) 14.14 (bottom) for the numerical solutions on L ) 1024 grid size. The integration time steps used are ∆t ) 0.000 25 and 0.005 s for each of the two cases, respectively. Note that the widths are shifted vertically by a factor of 120 and 25 respectively for the upper and lower fits. (b) Time evolution of the interfacial width in a semilogarithmic scale for L ) 40 plotted using the same simulation parameters as the top picture. A straight line through the data points indicates the logarithmic character of the time evolution of the width.

determined the growth of the interfacial width, the RHEED intensity, and have generated “snapshots” of the surfaces. In all cases we find logarithmic asymptotic behavior which is consistent with the Edwards-Wilkinson model, although in some instances there is slow crossover from effective power law behavior for short times. This result strongly suggests that any experimental determination of the nature of the growth of the interfacial width must cover a number of decades (in the number of layers grown) before any firm conclusion can be drawn. This observation also raises the possibility that existing numerical integration studies of the 2 + 1 dimensional KPZ equation may not have reached the long time regime. There are clearly additional complexities which must be added

Figure 9. Dynamic scaling results shown for the (UDP) model in a semilogarithmic scale. The parameters used for the simulation are flux ) 1015/cm2 s and kBT/J ) 0.6. The scaled width W ˜ is plotted against the scaled time ˜t ) number of layers added/Lz. The deposition time is proportional to the number of deposited layers. z ) 1.61 is used to obtain the scaling of data, which shows good scaling for the larger lattice sizes and obeys the logarithmic scaling function.

to these models to make them material specific, but we believe that the capabilities of the approach have been amply demonstrated. Acknowledgment. We wish to thank J. Krug for enlightening discussions. This research was supported in part by NSF Grant No. DMR-9405018. LA940705C