A Discrete Solid-on-Solid Model for Nonequilibrium Growth Under Surface Diffusion Bias S. Das Sarma and P. Punyindu Department of Physics, University of Maryland, College Park, MD 20742-4111 8 (February 1, 2008) 9 A limited mobility nonequilibrium solid-on-solid dynamical model for kinetic surface growth is 9 1 introduced as a simple description for the morphological evolution of a growing interface under random vapor deposition and surface diffusion bias conditions. Large scale stochastic Monte Carlo n simulations using a local coordination dependent instantaneous relaxation of the deposited atoms a produce complex surface morphologies whose dynamical evolution is not consistent with any of J the existing continuum dynamical surface growth equations. Critical exponents for coarsening and 9 roughening dynamics of the morphological evolution are quantitatively calculated using large scale 2 simulations. ] h PACS: 05.40.+j, 81.10.Aj, 81.15.Hi, 68.55.-a c e m An atom moving on a free surface is known to en- - counter an additional potential barrier, often called a PL PU t surfacediffusionbias[1](oranEhrlich[2]-Schwoebel[3] P P a L L st barrier), as it approaches a step from the upper terrace PL PU . – there is no such extra barrier for an atom approaching t a the step from the lower terrace (the surface step sepa- m rates the upper and the lower terraces). Since this diffu- sion bias makes it preferentially more likely for an atom - d to attach itself to the upper terrace than the lower one, n it leads to mound (or pyramid) - type structures on the o surface under growth conditions as deposited atoms are c [ probabilistically less able to come down from upper to lower terraces. This dynamical growth behavior under a 1 surfacediffusionbiasissometimescalledan“instability” v because a flat two dimensional surface growing under a 3 strong surface diffusion bias is unstable toward three di- 1 3 mensional mound/pyramid formation. There has been a 1 great deal of recent interest [1,4–20] in the morphologi- 0 cal evolution of growing interfaces under nonequilibrium 8 growth conditions in the presence of such surface diffu- 9 sionbias. InthisLetterweintroduce,whatwebelieveto / t be, the minimal nonequilibrium atomistic growth model a for ideal molecular beam epitaxial - type random vapor m deposition growth under a surface diffusion bias. Exten- - d sive stochastic simulation results presented in this paper n establish the morphological evolution of a surface grow- o ing under diffusion biasconditions to be rathercomplex. c Various critical growth exponents [21,22], which asymp- : v totically describe the large-scale dynamical evolution of i the kinetically growing surface in our minimal discrete X growthmodel,areinconsistentwithalltheproposedcon- r tinuumtheoriesfornonequilibriumsurfacegrowthunder a diffusion bias conditions. Our results lead to the con- clusion that a continuum description for nonequilibrium FIG.1. (a)Schematicconfigurationdefininggrowth rules growth under a surface diffusion bias does not exist and in 1+1 dimensions. (b) Morphologies (PU =1;PL =0.5) for may require a theoretical formulation which is substan- L=100×100 at 103 ML and (c) 106 ML. tially different from the ones currently existing in the literature. In Fig. 1(a) we schematically describe our solid-on- solid (SOS) nonequilibrium growth model : (1) Atoms 1 are deposited randomly (with an average rate of 1 layer/unit time, which defines the unit of time in the growth problem – the length unit is the lattice spacing taken to be the same along the substrate plane and the growth direction) and sequentially on the surface start- ingwithaflatsubstrate;(2)adepositedatomisincorpo- ratedinstantaneouslyifithasatleastonelateralnearest - neighbor atom (i.e. if it has a coordination of 2 or more since there must always be an atom underneath to satisfy the SOS constraint); (3) singly coordinated de- posited atoms (i.e. the ones without any lateral neigh- bors) could instantaneously relax to a neighboring site within a diffusion length of l provided the neighboring site of incorporation has a higher coordination than the original deposition site; (4) the instantaneous relaxation process is constrained by two probabilities P and P L U (with 0≤P ,P ≤1)whereP is the probabilityfor L U L(U) theatomtoattachitselftothelower(upper)terraceafter relaxation (note that a “terrace” here could be just one other atom). The surface diffusion bias is implemented in our model by simply taking P >P , making it more U L likely for atoms to attach to the upper terrace. Under the surface diffusion bias, therefore, an atom deposited at the top of a step edge feels a barrier (whose strength iscontrolledbyP −P )incomingdowncomparedwith U L an atom at the lower terrace attaching itself to the step. Our surface diffusion bias model is well-defined for any value of the diffusion length l including the most com- monly studied situation of nearest - neighbor relaxation (l = 1). We have carried out extensive simulations both in1+1and2+1dimensions(d)varyingP ,P aswellas L U l, also including in our simulations the inverse situation (the so-called‘negative’biascondition)withP >P so L U that deposited atoms preferentially come down attach- ing themselves to lower steps producing in the process a FIG. 2. (a) The mound evolution in 1+1 dimensions smooth growth morphology. Because of lack of space we showing a section of 500 middle lattice sites from a sub- do not present our negative bias results here except to strate size of L = 10000 at 104 ML, 105 ML, and 106 ML note that the smoothdynamical growthmorphologyun- (PU = 1;PL = 0.5). (b,c) The d=2+1 growth morphologies der our negativebias modelobeys the generic Edwards- on a 100×100 substrate (b) 103 ML, and (c) 106 ML. The darker(lighter) shadesrepresentlower(higher)pointsonthe Wilkinson universality [1,21,22]. surface. Before presenting our numerical results we point out two important features of our growth model : (1) For PL =PU =1ourmodelreducestotheoneintroducedin In Fig.1 we show representative d=2+1 simulated ref. 23(andstudiedextensively[21–27]inthe literature) growth morphologies for our positive bias model. The as a minimal model for molecular beam epitaxy in the diffusionbiasproducesmoundedstructureswhicharevi- absence of any diffusion bias; (2) we find, in complete sually statistically scale invariant only on length scales agreement with earlier findings [23,24] in the absence of much larger (or much smaller) than the typical mound diffusion bias, that the diffusion length l is an irrelevant size. In Fig.2 we show morphological evolutions for variable which does not affect any of our calculated crit- d=2+1 and 1+1 in order to visualize the coarsen- ical growth exponents (but does affect finite size correc- ing/steepening dynamics of the mound structures. We tions–increasinglrequiresaconcomitantincreaseinthe mention that in producing our final results we utilize a system size to reduce finite size effects). We, therefore, noise reduction [28]technique which accepts only a frac- mostly present our l =1 simulation results here empha- tion of the attempted kinetic events, and in the process sizing that our critical exponents are independent of l produces smoother results (reducing noise effects) with- provided finite size effects are appropriately accounted out in any way affecting the critical exponents. The cor- for. Our calculated exponents are also independent of responding results without noise reduction are visually the precise values of PU and PL as found in ref. 23 and more noisy with identical growth exponents. 24 for the P =P case. U L 2 To proceed quantitatively we now introduce the dy- namic scaling ansatz which seems to describe well all our simulated results. We have studied the root mean square surface width or surface roughness (W), the av- erage mound size (R), the average mound height (H), andtheaveragemoundslope(M)asfunctionsofgrowth time. We have also studied the various moments of dy- namicalheight-heightcorrelationfunction,andthesecor- relation function results (to be reported elsewhere) are consistentwiththeonesobtainedfromourstudyofW(t), R(t), H(t), and M(t). The dynamical scaling ansatz in the context of the evolving mound morphologies can be writtenaspowerlawsingrowthtime(whichisequivalent to power laws in the averagefilm thickness) : W(t)∼tβ ; R(t) ∼ tn ; H(t) ∼ tκ ; M(t) ∼ tλ ; ξ(t) ∼ t1/z, where ξ(t) is the lateral correlation length (with z as the dy- namical exponent) and β, n, κ, λ, z are various growth exponents which are not necessarily independent. In the > steady state, when ξ(t) ∼ L where L is the lateral sub- strate size, effective β vanishes as the surface roughness saturatesto avalue W (L)≡W(L,t→∞)∼Lα, where s α = βz is the roughness exponent. We find in all our simulations n ≈ z−1, and thus the coarsening exponent n, which describes how the individual mound sizes in- crease in time, is the same as the inverse dynamical ex- ponent in our model. We also find β = κ in all our results, which is understandable in a mound-dominated morphology. In addition, all our results satisfy the ex- pected exponent identity β = κ = n + λ because the mound slope M ∼ H/R. The evolving growth morphol- ogyisthuscompletelydefinedbytwoindependentcritical exponentsβ(thegrowthexponent)andn(thecoarsening exponent), which is similar to the standard(i.e. without anydiffusionbias)dynamicscalingsituationwhereβ and z (= n−1 in the presence of diffusion bias) completely define the scaling properties. We note also that our neg- ative bias results (not shown here for lack of space) are completelyconsistentwiththeexpected[1,22]Edwards- Wilkinson universalityclass [1,21,22]with our numerical findings being β =0.25 and z =2 in d=1+1, and β =0 (i.e. W ∼ lnt) and z = 2 in d=2+1. This is expected because our negative bias model explicitly introduces an inclinationdependentdownhillsurfacecurrent[29]inthe FIG.3. (a) The surface roughness W in 1+1 dimensions problem. as a function of deposition time t in the L = 10000 system. In Fig.3 we show our representative scaling results for Left inset: the growth exponent β calculated from the local nonequilibrium growth under surface diffusion bias con- derivativeoflog W withrespecttolog t. Rightinset: Av- ditions. It is clear that we consistently find β ≃ 0.5 in 10 10 erage mound size as a function of deposition time. (b) The bothd=1+1and2+1forgrowthunderasurfacediffusion surface roughness in the 100×100 system. Left inset: the bias. This β ≃ 0.5 is, however, very different from the localgrowthexponentβ. Rightinset: Averagemoundheight usualPoissongrowthunderpurerandomdepositionwith vstime. (c)Theaveragemoundsizeinthe100×100system. no relaxation where there are no lateral growth correla- Left inset: the local coarsening exponent n calculated from tions. Our calculated asymptotic coarsening exponent n the local derivative of log R with respect to log t. Right 10 10 inbothd=1+1and2+1isessentiallyzero(<0.1)atlong inset: theaverage mound slope vstime. times. In all our results we find the effective coarsening exponentshowingacrossoverfromn≈0.2atearlytimes toarathersmallvalue(<0.1)atlongtimes –webelieve the asymptotic n to be zero in our model. 3 In comparing with the existing continuum growth publications on this topic, is premature at this stage be- equation results (of which there are quite a few) we cause the real Ehrlich - Schwoebel barrier in experimen- findthatnonecanquantitativelyexplainallourfindings. tal systems [33]is expected [34]to be quite complicated, Golubovic [18]predicts β =0.5, which is consistent with and simple growth models used by us and others [4–20] our simulations, but his finding of n = λ = 1/4 in both most likely do not apply. What we have established in d=1+1and2+1isinconsistentwithourresults(n<0.1, this paper is that even a very simple limited mobility λ≃0.4−0.5)exceptinthelimitedtransientregime. The nonequilibriumgrowthmodelleadstoextremelycomplex analyticresultsofRostandKrug[19]alsocannotexplain dynamicalinterfacemorphologiesundersurfacediffusion ourresults,becausetheypredict,inagreementwithGol- bias conditions. The fact that even a deceptively simple ubovic, that if β = 1/2, then n = λ = 1/4, which is in limited mobility growth model such as the one studied disagreementwith ourfindings. We alsofindourβ to be in this paper seems to defy a theoretical continuum de- essentially0.5independentoftheactualvalueofP −P , scriptionisastrongindicationthatourunderstandingof U L whichdisagreeswithref. 19. Ourmodelobviouslyhasno nonequilibrium growth under a surface diffusion bias is slope selection (λ≡0), and therefore theories predicting far from complete. slope selections do not apply. This work is supported by the US-ONR and the NSF- Before concluding we point out some reasons for why MRSEC. we believe our nonequilibrium limited mobility growth model to be a good zeroth order description for kinetic growth under a surface diffusion bias. The main rea- son for this is the success of the corresponding mini- mal growth model, introduced in ref. 23, in provid- ing a good zeroth order description of molecular beam epitaxial growth in the absence of any surface diffusion [1] J. Villain, J. Phys. I 1, 19 (1991). [2] G. Ehrlich and F. G. Hudda, J. Chem. Phys. 44, 1039 bias. The d=2+1 critical exponents [27] in the unbi- ased model [23] are β = 0.25−0.2 and α ≃ 0.6−0.7, (1966). [3] R. L. Schwoebel and E. J. Shipsey, J. Appl. Phys. 37, which are in quantitative agreement with a number of 3682 (1966). experimental measurements [21,22] where surface diffu- [4] A. W. Huntet al., Europhys.Lett. 27, 611 (1994). sion bias is thought to be dynamically unimportant. An [5] Z.Zhang,J.Detch,andH.Metiu,Phys.Rev.B48,4972 equally significant theoretical reason is that the corre- (1993). spondingunbiasedgrowthmodel[23–27]istheonlyexist- [6] I. Elkinani and J. Villain, J. Phys. I 4, 949 (1994). ing nonequilibriumgrowthmodel whichis known[22,27] [7] P.SmilauerandD.D.Vvedensky,Phys.Rev.B52,14263 not to have the linear Edwards - Wilkinson ‘∇2h’ term (1995). [1] in its coarse-grained long wavelength continuum de- [8] M.SiegertandM.Plischke,Phys.Rev.E53,307(1996). scription (an equivalent statement [29] is that this is the [9] M. D.Johnson et al., Phys. Rev.Lett. 72, 116 (1994). onlylimitedmobilitymodelwhichhasavanishingsurface [10] P. Politi and J. Villain, Phys.Rev. B 54, 5114 (1996). currenton a tilted substrate)– in fact, our negative bias [11] J. A.Stroscio, et al. Phys.Rev.Lett. 75, 4246 (1995). model introduces precisely this ‘∇2h’ term by produc- [12] C. J. Lanczycki and S. Das Sarma, Phys. Rev.Lett. 76, ing a downhill current on a tilted substrate. (The other 780 (1996). limited mobility nonequilibrium growth models [30,31] [13] M. Siegert and M. Plischke, Phys. Rev. Lett. 73, 1517 introduced in the literature are known to belong to the (1994). Edwards-Wilkinsonuniversalityclass,andaretherefore [14] J. Krugand M. Schimschak,J. Phys.I 5, 1065 (1995). unsuitable for diffusion bias studies [32].) Therefore, the [15] M. Rost, P. Smilauer, and J. Krug, Surf. Sci. 369, 393 minimally biased version of this model which we study (1996). in this paper should be the appropriate zeroth order de- [16] J. Amarand F. Family, Phys.Rev.B 54, 14742 (1996). [17] J. Krug, J. Stat. Phys.87, 505 (1997). scription for nonequilibrium growth under surface diffu- [18] L. Golubovic, Phys. Rev.Lett. 78, 90 (1997). sion bias conditions. Since an approximate continuum [19] M. Rost and J. Krug, Phys.Rev. E 55, 3952 (1997). descriptionforthe unbiasedmodel[23]hasrecentlybeen [20] M. Siegert, M. Plischke, and R. K. P. Zia, Phys. Rev. developed [26], one could use that as the starting point Lett. 78, 3705 (1997). to construct a continuum growth model for the biased [21] A. L. Barabasi and H. E. Stanley, Fractal Concepts in growth situation. Such a continuum description is, how- Surface Growth(CambridgeUniversityPress,NewYork, ever, extremely complex [26] as it requires the existence 1995). ofaninfinite numberoftermsinthe growthequation. It [22] J. Krug, Adv.Phys. 46, 139 (1997). remainsunclearthatameaningfulcontinuumdescription [23] S.DasSarmaandP.I.Tamborenea,Phys.Rev.Lett.66, for nonequilibrium growth under a surface diffusion bias 325 (1991). is indeed possible even for our simple limited mobility [24] P. I. Tamborenea and S. Das Sarma, Phys. Rev. E 48, growth model. 2575 (1993). Weconcludebymentioningthatanycomparisontoex- [25] J.Krug,Phys.Rev.Lett.72,2907(1994); S.DasSarma periments, as has been done in severalrecent theoretical 4 et al., Phys. Rev.E 53, 359 (1996). [26] C. Dasgupta et al., Phys.Rev.E 55, 2235 (1997). [27] S. Das Sarma and P. Punyindu, Phys. Rev. E 55, 5361 (1997). [28] P.Punyinduand S.Das Sarma, unpublished(1997). [29] J. Krug, M. Plischke, and M. Siegert, Phys. Rev. Lett. 70, 3271 (1993). [30] F. Family, J. Phys.A 19, L441 (1986). [31] D.E.WolfandJ.Villian,Europhys.Lett.13,389(1990). [32] Ourmeasured surface current on tilted substrates is up- hill(downhill) for positive(negative) bias conditions be- comingzeroonlyintheunbiasedsituationofPU =PL – themeasuredcurrent,however,isanontrivialnonmono- tonic function of PU −PL. In the other limited mobility growth models, as in ref. 30 and 31, the surface current ontiltedsubstratesisalwaysdownhill,andtheytherefore must belong asymptotically to the Edwards - Wilkinson universality class (as in our negative bias model) inde- pendentof theirup/down incorporation probabilities. [33] K.Kyunoand G. Ehrlich, Surf.Sci. 383, L766 (1997). [34] S.Kodiyalam,K.E.Khor,andS.DasSarma,Phys.Rev. B 53, 9913 (1996). 5