UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Epitaxial growth dynamics in gallium arsenide Ballestad, Anders 2005

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


831-ubc_2005-103787.pdf [ 14.67MB ]
JSON: 831-1.0092343.json
JSON-LD: 831-1.0092343-ld.json
RDF/XML (Pretty): 831-1.0092343-rdf.xml
RDF/JSON: 831-1.0092343-rdf.json
Turtle: 831-1.0092343-turtle.txt
N-Triples: 831-1.0092343-rdf-ntriples.txt
Original Record: 831-1.0092343-source.json
Full Text

Full Text

Epitaxial Growth Dynamics in Gallium Arsenide by Anders Ballestad B .A.Sc , The University of British Columbia, 1996. M . A . S c , The University of British Columbia, 1998. A THESIS S U B M I T T E D IN P A R T I A L F U L F I L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R O F P H I L O S O P H Y in The Faculty of Graduate Studies (In Physics) T H E U N I V E R S I T Y O F BRITISH C O L U M B I A March 31, 2005 © Anders Ballestad, 2005 Abstract n A b s t r a c t The problem of a complete theory describing the far-from-equilibrium sta-tistical mechanics of epitaxial crystal growth remains unsolved. Besides its academic importance, this problem is also interesting from the point of view of device manufacturing. In order to improve on the quality and per-formance of lateral nanostructures at the lengthscales required by today's technology, a better understanding of the physical mechanisms at play dur-ing epitaxial growth and their influence on the evolution of the large-scale morphology is required. In this thesis, we present a study of the morphologi-cal evolution of GaAs (001) during molecular beam epitaxy by experimental investigation, theoretical considerations and computational modeling. Ex-perimental observations show that initially rough substrates smooth during growth and annealing towards a steady-state interface roughness, as dic-tated by kinetic roughening theory. This smoothing indicates that there is no need for a destabilizing step-edge barrier in this material system. In fact, generic surface growth models display a much better agreement with exper-iments when a weak, negative barrier is used. We also observe that surface features grow laterally, as well as vertically during epitaxy. A growth equa-tion that models smoothing combined with lateral growth is the nonlinear, stochastic Kardar-Parisi-Zhang (KPZ) equation. Simulation fits match the experimentally observed surface morphologies quite well, but we argue that this agreement is coincidental and possibly a result of limited dynamic range in our experimental measurements. In light of these findings, we proceed by developing a coupled growth equations (CGE) model that describes the full morphological evolution of both flat and patterned starting surfaces. The resulting fundamental model consists of two coupled, spatially dependent rate equations that describe the interaction between diffusing adatoms and the surface through physical processes such as adatom diffusion, deposition, and incorporation and detachment at step edges. In the low slope, small amplitude limit, the C G E model reduces to a nonlinear growth equation similar to the K P Z equation. From this, the apparent applicability of the K P Z equation to surface shape evolution is explained. The C G E model is based on fundamental physical processes, and can therefore explain the un-Abstract iii derlying physics, as well as describe macroscopic pattern evolution during growth. Contents iv C o n t e n t s Abstract ii Contents . iv List of Tables vii List of Figures viii List of Notation . xi Acknowledgements xiii 1 Introduction 1 2 The crystal growth problem 6 2.1 Adatom random walk 7 2.2 The solid-on-solid model 9 2.3 The Burton-Cabrera-Frank model 12 2.4 Inter layer transport 14 2.4.1 The Ehrlich-Schwobel step edge barrier 15 2.4.2 Step-traversal by insertion 16 2.4.3 Downhill funneling / knock-out process 18 2.4.4 Step traps 19 2.5 Interlayer transport in the Burton-Cabrera-Frank model . . . 21 2.5.1 The "Magic Slope" 24 2.6 Adatom desorption 25 2.7 A surface growth model derived in the Burton-Cabrera-Frank framework 27 2.8 Summary 29 3 Conventional continuum modeling of homoepitaxial growth 30 3.1 Continuum growth equations 30 3.2 Linear growth equations 31 Contents v 3.2.1 The Edwards-Wilkinson equation 32 3.2.2 The Ehrlich-Schwobel barrier in the continuum limit . 33 3.2.3 Mullins' diffusion 36 3.3 Nonlinear growth equations 37 3.3.1 The Kardar-Parisi-Zhang equation 38 3.3.2 The molecular beam epitaxy equation 39 3.3.3 Conservative Kardar-Parisi-Zhang equation 39 3.3.4 The Das Sarma equation . . 40 3.4 Effect of growth equation terms on surface shape evolution . 41 3.5 Universality classes 41 3.6 Power counting and relevance of terms 44 4 Epitaxial growth experiments on GaAs (001) 47 4.1 Simulations 48 4.2 Experimental method . 49 4.3 Results and analysis 51 4.3.1 Growth on thermally desorbed substrates 51 4.3.2 Growth on hydrogen etched substrates 62 4.4 Discussion 67 4.5 Summary 68 5 Modeling of growth on patterned surfaces 71 5.1 Introduction 72 5.2 Conventional modeling of weak surface texture 72 5.3 Sub-monolayer rate equations 75 5.4 Multilayer mean-field models 77 5.4.1 "Birth-death" models and diffusive growth 77 5.4.2 Evolution of the step density 79 5.4.3 Reflection high energy electron diffraction oscillations 80 5.5 Coupled growth equations model 82 5.5.1 Notes on the solution of the coupled growth equation model 84 5.6 Coupled equations in the asymptotic limit 85 5.7 Evolution of patterned surface gratings 92 5.7.1 Film thickness evolution 92 5.7.2 Temperature evolution 94 5.7.3 Transient amplitude overshoot 96 5.8 Summary 99 6 Conclusion 101 Contents v i A The kinetic Monte Carlo algorithm in solid-on-solid models 105 B Details of the Burton-Cabrera-Frank calculations 109 B . l Boundaries with Ehrlich-Schwobel barriers 109 B.2 Downhill funneling, Ehrlich-Schwobel barriers and detachment 110 B.3 Downhill funneling, Ehrlich-Schwobel barriers, detachment and desorption I l l C Finite difference scheme for nonlinear derivatives 112 D Scaling of coupled growth equations 115 D . l Lateral scaling: x —* bx 115 D.2 Vertical scaling: h —> ch 117 D.3 Time scaling: t -* dt 118 D.4 Composite scaling: (x,h,t) —> (bx,ch,dt) 119 D.5 Spatial frequency dependence 119 Bibliography 121 List of Tables vii L i s t o f T a b l e s 2.1 Standard Parameters used in kinetic Monte Carlo solid-on-solid simulations 10 3.1 Critical exponents associated with various growth equations in (2+1) dimensions 43 4.1 Summary of growth parameters for GaAs homoepitaxy ex-periments 50 5.1 Parameter table for coupled growth equations calculations . . 94 List of Figures viii L i s t o f F i g u r e s 2.1 Sample of simulated solid-on-solid simulation on flat starting surface 9 2.2 Unit cell of zincblende GaAs 11 2.3 Solution to diffusion equation on ideal terrace 13 2.4 Illustration of the Ehrlich Schwobel potential barrier 15 '2.5 Solid-on-solid kinetic Monte Carlo simulation showing effect on the large scale morphology of growth on flat starting sur-faces on the sign of the Ehrlich-Schwobel barrier 17 2.6 Diagrammatic illustration of "downhill funneling" 19 2.7 Alternate step edge potential diagrams: the step trap . . . . 20 2.8 The Bales-Zangwill instability 21 2.9 Adatom concentration on terrace with incorporation radius, Ehrlich-Schwobel barrier and detachment . 23 2.10 Effect of Ehrlich-Schwobel barrier and incorporation radius on adatom current 25 2.11 Homoepitaxy of Pt on P t ( l l l ) ; example of "magic slope" . . 26 3.1 Effect of positive Ehrlich-Schwobel barriers on vicinal surfaces (ID) 34 3.2 Solid-on-Solid simulation showing effect of the Ehrlich-Schwobel barrier on the morphology of a vicinal surface (2D) 35 3.3 Effect of Mullins' term for surface diffusion: surface peak-to-valley amplitude overshoot 37 3.4 Effect of the Das Sarma term on a square starting surface. . . 41 3.5 Effect of linear and nonlinear spatial derivatives in the growth equation on the surface evolution 42 4.1 Atomic force microscopy images of a GaAs wafer (a) after thermal removal of the surface oxide and (b) after 75 minutes growth at 550° 51 List of Figures i x 4.2 10 x 10 fxm2 atomic force microscopy images of a GaAs wafer with growth times (a) 10 minutes; (b) 37.5 minutes and (c) 150 minutes at 550° 52 4.3 10 x 10 nm2 simulations generated using the conservative form of the Kardar-Parisi-Zhang equation with growth times of (a) 10 minutes, and (b) 37.5 minutes • . 53 4.4 10 x 10 fxm2 simulations generated using the fourth order molecular beam epitaxy equation with growth times of (a) 10 minutes, and (b) 37.5 minutes 54 4.5 Comparison of power spectral density after different growth times on thermally cleaned substrates, measured along or-thogonal crystal directions. Both experimental results and simulations based on the Kardar-Parisi-Zhang equation are shown 56 4.6 Comparison of power spectral density. after different growth times on thermally cleaned substrates, measured along or-thogonal crystal directions. Both experimental results and simulations based on the molecular beam epitaxy equation are shown 57 4.7 Diffuse light scattering signal at 16 and 41 / im" 1 measured during growth of sample prepared by thermal oxide desorption 58 4.8 Scan lines from (a) measured atomic force microscopy im-ages, (b) the corresponding Kardar-Parisi-Zhang simulations, (c) molecular beam epitaxy equation simulations, and (d) Edwards-Wilkinson equation simulations 60 4.9 10 x 10 p 2 (a) atomic force microscopy image of thermally desorbed substrate, grown at 600° C along with (b) Kardar-Parisi-Zhang simulation 62 4.10 5 x 5 ^ m 2 atomic force microscopy images of (a) GaAs wafer which has had the surface oxide removed by hydrogen etching; (b) 75 minutes growth at 595°C and (c) Kardar-Parisi-Zhang simulation 63 4.11 Power spectral density of samples grown on hydrogen etched substrates measured along orthogonal crystal directions. Mea-surements and Kardar-Parisi-Zhang simulations are compared 64 4.12 2 x 2/ im 2 atomic force microscopy images illustrating the ef-fect of reducing the As2 overpressure, (a) V d l l ratio — 6.5, (b) V d l l ratio = 3.0, and (c) V:III ratio = 1.0 66 List of Figures x 5.1 Light scattering during growth corresponding to surface power spectral density at 41 /xm - 1 showing the effect of atom depo-sition on the smoothing rate of the surface 73 5.2 Atomic force microscopy images of (a) a sample quenched (fast cooled) after 69 minutes of growth at 600° C and (b) a sample annealed for 15 minutes at growth temperature 595° C after 40 minutes of growth . . 74 5.3 Evolution of sub-monolayer rate equations 76 5.4 Simple layer filling for the birth-death model 78 5.5 Layer filling for SOS simulation at 450°C and 625°C 79 5.6 Simple model of reflection high energy electron diffraction: SOS simulation with F = l ML/s 81 5.7 Comparison of growth on thermally desorbed surface with simulation based on the reduced coupled growth equation model 89 5.8 Solid-on-solid simulations showing strong coupling transition for increasing aspect ratio starting surfaces 91 5.9 Comparison of atomic force microscopy scanlines with cou-pled growth equation calculations and solid-on-solid simula-tions: film thickness dependence 93 5.10 Comparison of atomic force microscopy scanlines with cou-pled growth equation calculations: temperature dependence . 95 5.11 Temporal peak-to-valley amplitude evolution according to cou-pled growth equation model. Experimental values are also shown 97 5.12 Peak-to-valley amplitude evolution for corrugated surfaces: (a) ID calculations based on coupled growth equation model, and the inset shows actual surface evolution; (b) experimental data from regrowth on pitted starting surface 98 C . l Illustration of normal growth algorithm on discretized sur-faces:- (a) ID scheme and (b) 2D*scheme . . . . . . . . . . . . 112 List of Notation xi L i s t o f N o t a t i o n a . . . . roughness exponent: Wsat ~ La; also adatom-to-adatom incorporation coefficient ay . . . . in-plane inter-atomic spacing aj_ . . . monolayer height A F M . atomic force microscope (3 growth exponent: W ~ t@, t small; also adatom-to-step-edge incorporation coefficient BCF . Burton-Cabrera-Frank CGE . coupled growth equation d spatial dimension index D, Do diffusion constant, diffusion constant prefactor Ds ... single particle diffusion constant T](x,t) stochastic, random Gaussian deposition noise Eact • • surface activation energy Ees .. ES potential barrier Eiat .. lateral activation energy Esub . substrate binding energy ES . . . Ehrlich-Schwobel EW .. Edwards-Wilkinson F flux term or growth rate of surface r n , T c (Non-) conservative noise prefactor in continuum growth equations Q . . . . perturbation operator introduced in Chapter 5 h Planck's constant, surface height HE . . . hydrogen etched List of Notation xii 7i spatial Fourier transform of the surface height, h j net current of surface adatoms K coefficient of linear fourth-order term in Mullins' equation ks Boltzmann's constant K step edge release rate K M C . . kinetic Monte Carlo K P Z . . . Kardar-Parisi-Zhang KS Kuramoto-Sivashinsky A coefficient of non-linear term in K P Z equation A i 3 . . . . coefficient of non-linear term in Das Sarma equation A2,2 • • • • coefficient of non-linear term in Villain equation I terrace width lea Ehrlich-Schwobel length I* terrace width for which the adatom current is zero u, chemical potential M B E . . molecular beam epitaxy M L . . . . monolayer v coefficient for surface relaxation term in E W PSD . . . power spectral density q spatial frequency p, P, n . adsorbate particle density; also lattice site occupation probability R H E E D reflection high energy electron diffraction Rinc • • • incorporation radius SOS . . . solid-on-solid t, At . . . time, time increment T D thermally desorbed Tdes average time between adatom desorption events 6j filling of jth layer W interface width given by rms fluctuation (x, y) . . surface coordinates z dynamic exponent: z = a/(3 Acknowledgements xiii A c k n o w l e d g e m e n t s Thanks to Tom Tiedje, my supervisor, who is so famously easy to get along with and whose criticism always comes in the form of constructive encourage-ment. In spite of our fundamental differences regarding views on grammar and sentence structure, and commas, we seem to have been able to turn out a couple of papers where we agree in full on the scientific issues in question. Tom loves physics, and his knack for finding new and interesting projects is the primary reason I'm graduating as early as I am. He also loves the mountains of British Columbia, and there has many occasions where I have thought that I might have to finish my degree under a new supervisor: Tom and the boulder; Tom in the cravass; other stories come to mind. Tom is a great supervisor, and I consider myself very lucky to have worked with him. Thanks to Birger Bergersen, who took over the supervisorly role after Tom fled to Japan just prior to my graduation. Birger is living proof that physics is indeed a lot of fun. Jeff Young has been extremely supportive over a long time, and actually introduced me to Tom back a long time ago. His interest in motorcycles has probably been a long time motivator for me to want to be like Jeff; he's had them all: the Norton Commando, the V T R , the RC51... but not the Ducati 900SS. Thanks for not coming to my final defense! Walter Hardy came in at the last minute to help out as a committee member. He managed to dissect my thesis within two weeks and ask the most gruelling questions on my departmental defense. How do you do it? Walter deserves additional thanks for the many parties that he has organized over the years, and of course for purchasing the espresso machine in 242. All of my committee members have been very generous with their time and wisdom, and I am very grateful. In addition to my committee members, there are many people that de-serve thanks. These are people that I have harassed either by email, or by meeting them at conferences and other scientific occasions and that have helped me with comments and suggestions along the way. They are: Acknowledgements xiv Drs. Maria Bartelt 1, Matt Choptuik, Garry Clarke, Phil Cohen, Martin Grant, Mark Gyure, Makoto Itoh, Bruce Joyce, Hung-Chi Kan, Karen Ka-vanagh, Itzhak Kelson, Joachim Krug, Thomas Michely, Joanna Mirecki-Millunchick, Anthony Pierce, Mike Plischke, Christian Ratsch, Lou Ter-minello, Tony van Buuren, Raffaelle Vardavas, John Venables, Dmitri Vve-densky, Brian Wetton and Mao Hai Xie. All my love to Marianne: you have believed in me all this time. I couldn't have done this without you! There are so many people that have helped and influenced me in the course of my university life, too many to mention. Some of these people in-clude the current members of the MBElab: Mike Whitwick, Erin Young, Eric Nodwell, Nikolaj Zangenberg, Shawn Penson, Dan Beaton (from Mabou!), Kevin Mitchell, Scott Webster, Lilian Fan, Raveen Kumaran, Arvin Yazdi, and of course, Jim McKenzie; graduates of the lab: Jens Schmid, Martin Adamcyk, Tom Pinnington, Bayo Lau, Richard Mar, Mya Warren, Robin Coope, Christine Nicoll and others; Post-Docs: Ben Ruck, Mario Beaudoin, Yuval Levy, Sebastien Tixier; staff, professors and students in the depart-ment of Physics and Astronomy at U B C ; Tony Walters, our grad secretary and Maria Mannella in Finance. Thanks to Dr. Balzarini for organizing Friday hockey for some 34 years. Thanks to friends in Vancouver and elsewhere: Sune Hoifeldt, Arne Hoogensen, Dmitri Grandmaison and Mehrdad Max Mafmezam; thanks to Linux, Lego, the Internet and my iPod. Thanks to all those that I have forgotten in the moment of writing this. Finally, I am grateful to my family: mamma and Oluf, pappa and Gunn, mormor, Svanborg and Bjarne, Hanna, Aina, Helene and Jan-Oluf; Mr. and Mrs. Chang; Liam, Cloe, Joanne and Warren. 1 Deceased. Chapter 1. Introduction 1 Chapter 1 In t roduc t ion The shape of the growth surface is a basic property of epitaxial films that contains information about the atomic scale processes which take place dur-ing growth by molecular beam epitaxy (MBE). Theoretical investigations have shown that in the long wavelength and small amplitude limit, surface shape evolution during growth can be described by continuum equations whose terms reflect the underlying microscopic processes. The asymptotic scaling behaviour of various possible growth equations has been studied extensively theoretically, with both analytical and numerical methods. Sim-ilarly, experimentalists have attempted to verify the theoretically predicted growth exponents and scaling behaviour in real systems. However, experi-mental studies of scaling are inevitably limited by practical constraints, such as rather small ranges of time and distance, so that accurate measurements of growth exponents are difficult to achieve. As a result, there are few, if any, experimental studies that have been successful in using scaling analysis to show a quantitative match between experimentally grown surfaces and simulated surfaces based on theoretical growth models. We have investigated the morphology of single crystal semiconductor surfaces during epitaxial growth and attempted a qualitative understanding along with a quantitative model of its evolution. Having such a model is im-portant for many reasons. From a technological point of view, it allows fab-rication of lateral nanostructures by epitaxy. It also enables one to predict and control surface structures in devices, and thus to make better materials through better understanding of epitaxial growth. From a theoretical point of view, since there is still no underlying theory for this far-from-equilibrium problem, one can use the large-scale topography in order to understand mi-croscopic atomic processes. Our work aims to relate theoretical efforts in surface morphology modeling to experimental work on the technologically important GaAs (001) surface during homoepitaxial crystal growth on both nominally flat and patterned starting surfaces. There has been a large amount of research on surface morphology over the last 15 years, and we will use the greater parts of Chapters 2 and 3 to de-scribe some of it. Villain [Vil91] developed a framework for the use of contin-Chapter 1. Introduction 2 uum growth models by the means of partial differential equations, where the interface velocity was related to free energy gradients. Venables [VSH84], Bartelt [BE92], and others developed sub-monolayer rate equations that have seen great success for the early stages of growth and nucleation. Co-hen [CPP + 89] extended the rate equations to account for coalescence of islands on terraces and, thus, enabled the description of multilayer growth. However, morphological information is lacking in a rate equation approach. Many groups have investigated a variety of spatially dependent con-tinuum and discrete growth models and analysed their scaling behaviour in the hydrodynamic limit. A complete picture of the work done is best obtained by reading the review articles and books by for example Halpin-Healy [HHZ95], Krug [Kru97], Barabasi and Stanley [BS95], and Pimpinelli and Villain [PV98]. Recently, there has been a renewed interest in surface morphology studies due to experimental studies of overgrowth on textured surfaces for quantum wires and distributed feedback laser design [ B K T + 0 1 , TWY+03]. The nature of the adatom dynamics near step edges has been a theoreti-cal favourite in many studies. It is believed that the descent of adatoms from one terrace to the next ("interlayer transport") is inhibited by an asymmetric potential barrier at the step edges. Villain [Vil91] showed that such a poten-tial barrier (also referred to as the Ehrlich-Schwobel barrier [EH66, SS66]) causes preferential incorporation of adafoms to the ascending step on a terrace, along with reflection off the descending step, and that large-scale mounds would emerge as a result. Pioneering work by Johnson [JOH +94] and Orme [OJS +94, OJL+95] showed that the evolution of GaAs surfaces during growth exhibited the formation of such mounds and they concluded that this destabilizing effect could be attributed to the Ehrlich-Schwobel barrier and that GaAs growth was unstable. From an experimental point of view, this conclusion came as somewhat of a surprise. Long before theorists started modeling GaAs homoepitaxy, crystal growers had been growing thick buffer layers (on the order of microns) prior to device manufacturing in order to smooth the starting surface. This experimental knowledge indicated the exact opposite of what the potential barrier studies had found: namely that homoepitaxial growth on GaAs (001) was stable and that continued growth would make rough starting surfaces smoother, rather than filled with mounds. The disparity between these two views on the stability of GaAs epitaxy remained for some time. A stable, smoothing morphological evolution means that the net di-rection of adatom drift is downhill. This can be achieved in two ways: (1) by letting the step edge barrier be negative, which leads directly to Chapter 1. Introduction 3 downhill interlayer transport, or (2) by having a positive Ehrlich-Schwobel barrier along with other mechanisms that oppose the effects of this desta-bilizing mechanism. Some such mechanisms have been described in the literature, for example "downhill funneling" [WV90] and "the knock-out effect" [Eva91]. In fact, Smilauer and Vvedensky [vV93] showed that in order to accurately describe reflective high energy electron diffraction os-cillations of GaAs homoepitaxy, a downhill funneling mechanism was re-quired in addition to the positive Ehrlich-Schwobel barrier [OJS+94]. Krug et al. [KPS93] showed that this combination of opposing mechanisms at the step edge could in some cases lead to a net downhill drift of adatoms. The combined results of Refs. [vV93] and [KPS93] lends additional support to the notion that the evolution of GaAs homoepitaxially grown surfaces is stable. New studies have recently emerged that suggest that the Ehrlich-Schwobel barrier in GaAs might be small [CCM+98, CCOO, IO00] or even negative [BTS04a, BTS+04b]. In this thesis, we approach the problem of surface shape evolution by direct comparison of the morphology of epitaxially grown surfaces with sur-faces generated by computer simulation of various growth models. Experi-mental surface morphologies were measured with atomic force microscopy, which gives atomic resolution in the surface height over a wide range of sub-strate sizes. We have also tried to match the power spectral densities as a function of time in growth experiments on nominally flat surfaces. We dis-covered [ABP+00, BRA+01, BRS+02] along with other groups [CCM+98, CCOO] that the observed mounding can be explained as being a transient smoothing phenomenon of a rough starting surface, rather than being due to a growth instability inherent in the GaAs surface caused by a positive Ehrlich-Schwobel barrier. We therefore agree with the view that the mor-phological evolution of this material system is stable during epitaxial growth. This realization enabled us to model the surface shape evolution with growth equations that contained stable (non-diverging) terms. This modeling, along with accompanying experimental studies is discussed in Chapter 4. Continuum models are intended to describe surface morphology in the hydrodynamic limit, i.e. for low slopes and small surface roughness ampli-tudes. The question arises whether our modeling approach will work for lithographically patterned starting surfaces. The answer to this question was made particularly evident in recent work by Kan et al. [KSTEP04b], in which modeling of regrowth on 50 nm deep lithographically defined pits was attempted using conventional continuum growth equations, like the Kardar-Parisi-Zhang (KPZ) equation [KPZ86] and the so-called molecular beam Chapter 1. Introduction 4 epitaxy equation1 [SGG89]. Kan's calculations showed that several impor-tant and very noticeable characteristics of the surface shape evolution were completely missed by attempting to use these equations as a description of real crystal growth. Our analysis of this extended problem led us to formu-late a growth model that describes the interaction between adatoms and step edges in the form of a coupled growth equation model [BTS04a, BTS + 04b]. This work can be seen in Chapter 5. The progression of the Chapters is as follows. Due to the extensive body of work related to modeling of epitaxy, we feel that a thorough intro-duction to the field is appropriate. This appears in the form of a crit-ical review in Chapters 2 and 3. The material in Chapters 4 and 5 is new, and has resulted in a number of publications on growth modeling (see Refs. [ABP+00, BRA+01, BRS+02, BTS04a, BTS + 04b , BLST05]) and also the reverse process of etching (see Refs. [SBR+02, STMB04]). In Chapter 2, we review some basic aspects of adatom diffusion on crystal surfaces and introduce key elements in an approach to a microscopic sur-face growth model called the "solid-on-solid" model. We also introduce the Burton-Cabrera-Frank (BCF) [BCF51] theory along with the accompany-ing mathematical framework. We discuss in some detail the potential-energy landscape of the GaAs surface, expecially close to the steps, and discuss its effect on macroscopic surface evolution. A conventional approach to large-scale surface morphology modeling is through the use of continuum growth equations. We describe some common examples of such in Chapter 3, and we also outline the concept of scaling, which is a typical approach used to classify growth models. We believe that the ultimate test of a model is to try and apply it to the experimental system that it seeks to describe2. This is attempted in Chap-ter 4, where two surface growth models for M B E are applied and compared in an exhaustive experimental and simulational study: We investigate the dependence of surface evolution on film thickness and temperature. We show that the K P Z equation is successful in describing surface shape evolution of nominally flat starting surfaces in the low amplitude and long wavelength limit, whereas the M B E equation, ironically, is not. We also study the effect 1We use the acronym "MBE" to describe "molecular beam epitaxy". We will also follow literature in referring to this equation as the "MBE equation". We hope that it is clear from the context whether we are referring to the growth technique (MBE) or the growth model (MBE equation). 2This might seem obvious, however, the reader will be surprised to learn that this is not typically done in morphology studies. Quite often, there is not even an image of a surface! Chapter 1. Introduction 5 of the starting surface on the surface evolution and show that the As to Ga ratio dictates the anisotropy of the evolving surface features. This study proves to be valuable, in that we learn many things about our particular growth system, and also come to the conclusion that GaAs homoepitaxy is in fact a stable growth system. We conclude Chapter 4 by noting that conventional asymptotic growth equations leave many questions unanswered. We progress to Chapter 5, where we develop a complete physical model based on the underlying phys-ical processes and mechanisms that take place on the surface during growth and annealing. We show that this new model can describe regrowth exper-iments on lithographically patterned surfaces, including the complex shape evolution associated with large surface slopes. It. can also describe the evo-lution of the nominally flat surfaces that we studied in Chapter 4. The new model reduces to a nonlinear evolution equation in the low amplitude and long wavelength limit reminescent of the K P Z equation. The four appendices describe details of the various mathematical tech-niques used in this work. Chapter 2. The crystal growth problem 6 C h a p t e r 2 The crystal growth problem Before we start the task of modeling surface-shape evolution, we need to describe the phenomena believed to be present at the crystal surface. We will discuss important physical mechanisms that influence the evolution of the surface morphology, both during epitaxial growth and annealing. This Chapter will be historical in nature, and aims to discuss some major ideas that have affected the way we understand the crystal growth problem. Each square micron of GaAs surface contains approximately 6 million exposed atoms, most of them locked into the 3D crystal lattice and some diffusing freely on the surface, constantly rearranging and working hard to develop macroscopic surface shapes that keep on fascinating surface scien-tists. Due to the large number of particles involved, the development of complete, atomistic models for surface evolution during growth inevitably results in very calculationally intensive models. Practical models should be computationally tractable, while not jeopardizing their applicability to real surface growth evolution. We start this Chapter by discussing two common approaches to growth modeling. An old and often quoted theory is the one developed by Bur-ton, Cabrera and Frank (BCF) [BCF51]. It describes the motion of loosely bonded adatoms on a crystal surface consisting of flat terraces separated by atomic step edges. This model was developed over half a century ago, long before any scanning tunneling or atomic force microscopes revealed to us the secrets of the crystal surface. Since then, it has formed the basis for many theoretical models of adatom motion and nucleation on crystal surfaces, such as sub-monolayer rate equations. In this Chapter, we will incorporate many realistic physical processes into the BCF framework, and from that, develop a continuum growth model, thereby showing the strengths and weaknesses of this simple, yet remarkably useful formalism. The other approach that we will cover in this Chapter is the solid-on-solid (SOS) model, in which the stochastic nature of crystal growth is explored. In the development of both these models, we will pay special attention to the physics of adatoms near step-edges. The variations in the potential landscape near steps is found to have a profound impact on the evolution of Chapter 2. The crystal growth problem 7 large-scale surface morphology during growth. 2.1 A d a t o m r a n d o m w a l k The success of molecular beam epitaxy (MBE) as a method for crystal growth is in part due to the clean, high vacuum environment, but also in great part due to the fast diffusion of adatoms on the surface. While growth is done at only about half the melting temperature, which is at 1237° C in the case of GaAs, this, is still well above room temperature, and it allows sufficient adatom migration to achieve relatively flat films. An "adatom" is defined as a single atom, whose only bonds are to the crystal lattice beneath it. It can easily overcome the potential barrier re-quired to jump diffusively from one lattice position to a neighbouring site without desorbing back into the gas phase, since it is still in intimate con-tact with the surface. This diffusion is Arrhenius activated and allows the adatoms to seek out low-energy sites on the surface. During early stages of growth, when there are fewer steps on the surface, adatoms nucleate dimer-islands along with other adatoms. As the growth progresses, or on surfaces with some initial roughness, the most common way that adatoms minimize their potential energy is by attaching to step edges. Adatoms can also desorb back into the vapour. The diffusion constant D = a|ajn exp (—E a c t/kBT), where Eact is the energy that must be overcome in order for a diffusion event to occur, kg is Boltzmann's constant, T is the temperature of the system, h is Planck's constant and ay is the in-plane lattice constant. The attempt frequency LUQ has a weak temperature dependence through the equipartition theorem as u0 = 2kBT/h. The simplest models of GaAs growth usually do not make the distinction between atomic and molecular constituents, and an adatom is considered as a representation of a Ga-As complex. A "dimer" (two Ga-As complexes) is considered as part of the surface, although it can diffuse as a unit with a much larger activation energy than that of a single adatom. We will consider dimers as immobile in this thesis. Mobile islands have been treated for example in Refs. [VSH84, LGOO, KEOO]. We define P(x,y,t) to be the probability of an adatom occupying a lattice site at position (x, y) at time t. After a diffusional step, we can write Chapter 2. The crystal growth problem 8 P(x,y,t + At) as [Zan88]: P(x,y,t + At) = i ( P(x + ahy,t)+P(x-ahy,t) +P(x,y + a\\,t)+P(x,y-a\\,t) ) (2.1) which can be rewritten as P(x,y,t + At)-P{x,y,t) = ^l ( P(x + all,y,t)-2P(x,y,t) + P(x-ahy,t) ) + P(,x,y + a\\,t)-2P(x,y,t)+P(x,y-a\\,t) ) (2.2) Dividing this equation by At, we recognise it as a second order finite differ-ence scheme for the heat equation: dtP(x,y,t) = DsV2P(x,y,t) (2.3) where Ds — a2/AAt is the diffusion constant of a single particle and dt = d/dt. In the case of non-interacting adatoms on a surface, Fick's law states that the flux of diffusing adatoms is proportional to the density gradient j = -DVp (2.4) where p = p(x, y, t) is the adsorbate particle density. Adding the continuity equation d 4 p = - V - j (2.5) we get the diffusion equation (Eq. 2.3) for D = Ds and p = P, which is true when there is no interaction between the adsorbate particles. A further discussion on this can be found in Ref. [Zan88] (p. 375), Ref. [PV98] (p. 112), as well as an in-depth analysis in Ref. [BMT90] (p. 529). We will use P and p interchangeably throughout this Chapter 1. When these are not interchangable, we will note this explicitly. 1 In Chapter 5 , we will refer to the adatom concentration by the variable name n, so that it should be clear from the context what formalism, we are using. Chapter 2. The crystal growth problem 9 .4 , # , i ( * < i> i (a) (b) Figure 2.1: Sample of simulated SOS growth on flat starting surface, (a) Perspective view after 1.9 M L growth; (b) top view after 65 ML growth on a 300x300 atom surface; higher layers in lighter shading. Growth parameters were F = l M L / s , T=500°C, E S U6=1.2oeV, E i a t =0.35eV and E e s =-0.05eV. 2.2 T h e s o l i d - o n - s o l i d m o d e l A widely used method in simulating epitaxial growth is the solid-on-solid model. Its name originates from the fact that it assumes that the smallest particles used in the model land and diffuse on top of other particles, i . e. no overhangs or vacancies are allowed. This means that all the surface atoms are sitting on top of a solid column of atoms, and that the surface can be fully described by a 2D matrix, containing only the height information of the surface atoms. A small area of a typical SOS simulation is depicted in Fig. 2.1. In this Section, we outline the specifications of our particular SOS model. In the simplest model, the surface atoms are defined on a square lat-tice with periodic boundary conditions, and individual atoms are deposited randomly around the surface at a rate F. Each surface atom is bonded to the substrate directly below it with a bonding potential Esub- Atoms with m lateral neighbours feels an additional bonding of mEiat, where m ranges from 1 to 4 in a cubic lattice [MG88]. The total activation energy for a surface atom is then Eact = Esub + mEiat. The surface atoms diffuse at rates dictated by their bond-strength: a weakly bonded atom, such as an adatom, is more likely to move than an atom with many lateral neighbours. An atom with coordination number m = 4 is part of a flat surface, and is also allowed to diffuse, although at a much lower rate. We also include an Ehrlich-Schwobel (ES) barrier, which describes the additional energy bar-rier that an atom encounters as it approaches a step edge. This barrier is characterized by a potential Ees. We discuss this barrier in more detail in Chapter 2. The crystal growth problem 10 Table 2.1: Standard Parameters used in K M C SOS simulations in this work F Elat Fes l . O M L / s 1.25 eV 0.35 eV -0.05 eV Section 2.4. In this work, unless otherwise noted, we will use a fixed set of parameters in all our SOS simulations, which we will refer to as the stan-dard parameters. These parameters are listed in Table 2.1 and were obtained through optimization and comparison between a wide variety of simulation scenarios, and related experiments, whenever available. The kinetic Monte Carlo ( K M C ) algorithm is well suited for simulating the inherent randomness of the SOS dynamics. More details on the K M C approach are outlined in Appendix A . Once the computational method has been established, we can focus on the growth model itself, or the "rule-base" that makes this simulation imitate epitaxial growth as closely as possible. We use a one-component model, which means' that the As-atoms are not distinguished from the Ga-atoms. The smallest unit of surface particle is thus a Ga-As complex. This is reasonable, as GaAs epitaxial growth is usually done with excess As present in the growth chamber, which makes the growth Ga-limited. However, the As-overpressure is known to affect the anisotropy in the evolving surface morphology, where surface features are stretched in the [110] direction when the As:Ga ratio decreases [ABP+00]. In our growth experiments, we usually set this .ratio to be > 3, for which the surface feature anisotropy is minimal. The surface, anisotropy could be modeled indirectly by assuming an anisotropic diffusion constant for the Ga-As complex or possibly also an anisotropy in the incorporation dynamics on the surface [ItoOl]. The use of a one-component model effectively reduces the crystal struc-ture from a zincblende structure (two face-centered cubic (FCC) crystals intertwined, offset by a quarter lattice spacing in the [111] direction) to a single F C C crystal with a basis representing the Ga-As complex. In our model, this is further simplified to a simple cubic crystal, which can be justified by interpreting each layer in the simulation as a bi-layer, and let-ting the simulation (x, y) directions be the crystal [110] and [110] directions, respectively, see Fig. 2.2. We have implemented the restricted SOS model [KK89, RSLP91], in which there can be no double-height steps on the surface. This is a reflection of the fact that a Ga-As complex right at the edge of a step represents a high-energy, and therefore improbable configuration. B y restricting the Chapter 2. The crystal growth problem 11 Figure 2.2: Unit cell of zincblende GaAs, with As-atoms shown in black and Ga-atoms in grey. By letting x in the K M C simulations symbolize the [110] direction, and y the [110] direction, the cubic crystal used in K M C can be justified. height at the step edge, the FCC lattice of GaAs is more closely simulated. The macroscopic effect is a restriction of the slope to below 45°. We do not consider this to be a limitation of the model's accuracy, as the next low index planes, (110) and (111), are 45° and 54.7° from the surface normal of the (001) plane. The surfaces investigated experimentally in this thesis have slopes less than ~ 5° for the randomly rough surfaces described in Chapter 4 and less than ~ 30° for the patterned surfaces in Chapter 5. Elkinani et al. [EV94] and Politi et al [PV96] showed that an SOS model without the double-height step restriction leads to infinitely steep sections, due to step-bunching induced by diffusion of adatoms over steps. The double-height restriction is implemented in the following way. Any surface atom with less than 4 lateral neighbours is considered to be a step. Adatoms that land on top of such a step site, either by diffusing or being deposited there, is forced to make another diffusional step until it lands on top of an atom that has 4 lateral neighbours. Surface reconstruction, where the exposed atoms settle into reduced-potential configurations, is not included in our models. Surface reconstruc-tion is dependent on the substrate temperature, and also on the As2:Ga ratio. The inclusion of this effect could be done at a large computational expense, but would be better handled by a truly atomistic model imple-mented for example with molecular dynamics, in which surface atoms can be moved off their bulk lattice positions. In this work, we keep the min-imum number of mechanisms necessary to be consistent with experiments and remain consistent with generally accepted models of crystal growth. In spite of the many simplifications that go into our SOS model, the re-Chapter 2. The crystal growth problem 12 suits that emerge from this computational tool have proven to be predictive and accurate in describing a wide range of surface dynamics phenomena. If one views the SOS model as a box of Lego-pieces, i.e. not as an accurate description on the small scale, but great at resembling large-scale surface shapes, then one could call this method a computational playground for the surface theorist. In order to explore phenomena that are not readily available through experiment, we will resort to the SOS model, and con-sider systems with up to about ~1000xl000 atoms and growth times on the order of hours. 2.3 T h e B u r t o n - C a b r e r a - F r a n k m o d e l The fundamental assumptions in the theory of adatom diffusion according to B C F [BCF51] are that the surface consists of flat terraces separated by steps, and that the concentration of adatoms on a terrace is determined fully by the position of the step edges. This means that the re-ordering of the adatom concentration on a terrace is fast compared to the motion of the steps. A surface that is not perfectly flat can be considered as a collection of steps separated by top and bottom terraces. If the surface average normal is slightly tilted off a low-index crystal orientation, such as (001) or (111), it is called vicinal or stepped. A flat surface that is not vicinal is called singular. Figure 2.3 illustrates a vicinal surface along with some definitions that we will use in the following. With adatoms impinging on the surface at an average rate F, the equation that governs their steady-state configuration on a terrace of width Z is 0 = dtp(x, t) = DV2p(x, t) + F/af • (BCF equation) (2.6) where a | is the d-dimensional lattice site area2. In the simplest ID (d = 1) model, we will assume that the ascending step is perfectly absorbing, which gives us the boundary condition p(x = 0) = 0, (Absorbing boundary) (2.7) On the descending step, we allow drift of adatoms over the step without any additional barrier, which gives us the following condition: —DVp(x = I) = —p{x = I) (Zero barrier boundary) (2.8) 2The units of Eq. 2.6 in d dimensions is: [p] = 1/L d , [D] = L 2 / T , [F] = 1/T, where L is length and T is time. Chapter 2. The crystal growth problem 13 P(x) x=0 x=[ Figure 2.3: Solution to diffusion equation (Eq. 2.6) of vicinal terrace with perfectly ab-sorbing boundary at x = 0 (Eq. 2.7) and a zero energy barrier boundary at x = I (Eq. 2.8). where the left hand side represents the particle current at the step edge, and at the right hand side we have the number of jump attempts at the step edge. The breach of symmetry in the boundary conditions is a result of this sys-tem being in a nonequilibrium steady-state situation. A simple experimental proof of this fact is that the morphology evolves: without the asymmetric boundary conditions, the net current of adatoms off a terrace would be zero, and the surface shape would remain fixed3. The solution to Eq. 2.6 with boundary conditions Eqs. 2.7 and 2.8 is a simple quadratic: and is plotted on the diagram in Fig. 2.3. Fick's law (Eq. 2.4) will now give us the current of adatoms at the as-cending step, j u p , and the descending step, jdown- The destabilizing current j u p drives adatoms towards the ascending step, and hence make the terrace shorter. On the descending step, jdown will represent a measure of the num-ber of adatoms that fall over the edge and hence lengthen the terrace in question; such a stabilizing term acts to flatten the overall surface. For the adatom density defined by Eq. 2.9, we get the following expressions for the 3 Another example of a nonequilibrium steady-state system is that of a gas contained in a cylinder, with the two ends of the cylinder being held at different temperatures. One could here easily obtain a situation where the temperature reached a steady-state. Yet, the pressure of the gas near the two ends of the cylinder will be different, even in the steady-state, due to the fact that the system is in non-equilibrium. (2.9) Chapter 2. The crystal growth problem 14 two contributing currents: j U p ( 0 = -DVp\x=0 Fl ( 2 a n + Q 2a|| (ay + Z) j d o w n ( 0 . = - D V p | x = z F Z 2 ~~ 2ay (ay + Z) . The total current is the sum of these: j ( 0 = j u p ( 0 + jdcmm(Z) F Z ~ . ay + Z This expression is negative for all values of Z > 0. On the surface depicted in Fig. 2.3, where the gradient is negative, this means that the net adatom drift is uphill. Thus, the surface will get rougher as the growth progresses, by building up 3D structures. A net downhill adatom current, on the other hand, will cause an initially rough surface to smooth. In the next Section, we will allow for permeable boundaries at both ends of the terrace, and we will describe in more detail the physics of adatoms traversing step edges, which in the context of surface growth modeling is referred to as "interlayer trans-port" . The physics determining the permeability of the step edges is crucial in order to understand the.surface shape evolution on all interesting length scales. In a recent Nature article, aptly named "Thin-Film Cliffhanger", La-gally and Zhang [LZ02] write that: "A proper accounting for these barriers in film growth will determine the ultimate quality, reliability and stability of films grown for a great variety of purposes". We continue by discussing physical mechanisms that favour adatom drift in either of the two directions, uphill or downhill, and the consequenses of incorporating such mechanisms into our surface model. 2.4 I n te r l aye r t r a n s p o r t In order to create smoother surfaces from a rough starting surface, the adatoms have to flow downhill, making terraces longer, thereby depleting the top terraces and filling the troughs. There must be physical processes on the surface that allow for the adatoms to jump to terraces above or below the one they are currently on. In this Section, we discuss some important (2.10) (2.11) (2.12) Chapter 2. The crystal growth problem 15 Surface coordinate Figure 2.4: Illustration of the Ehrlich-Schwobel potential barrier. An adatom feels weaker bonding to the surface when traversing the step in a perpendicular direction. Immediately beneath the step, this situation is reversed, and an atom will feel additional bonding to the surface. The ES barrier, Ees, is indicated. mechanisms that would influence the nature of interlayer transport, includ-ing potential barriers at step edges (Section 2.4.1 and 2.4.2), local relaxation of newly deposited atoms (Section 2.4.3), as well as some of the consequenses of having competing mechanisms for interlayer transport (Section 2.4.4). F i -nally, we discuss some alternative theories (Section 2.5.1). 2.4.1 The Ehrlich-Schwobel step edge barrier A naive rendering of the potential landscape near the step-edge of a ID cubic lattice is illustrated in Fig. 2.4. On its journey across the surface, a diffusing adatom continously encounters potential energy barriers. Away from steps, this activation energy is caused mostly by the interaction be-tween the adatom and the atoms in the substrate immediately beneath it, characterized by Esui, in Fig. 2.4. Near to an ascending step, the adatom feels added bonding from the lateral bond to the step-atom. It is based on this argument that crystal growth is often described as being a dynamic and reversible process of attachment and detachment of step edge atoms (BCF theory). In order to descend over a step, the adatom has to overcome a poten-tial barrier that is somewhat larger than that which characterizes a plain diffusion event on a flat surface. The added potential barrier in this naive picture is the ES barrier of magnitude Ees, named after the experimental Chapter 2. The crystal growth problem 16 observations by Ehrlich and Hudda [EH66]. Some of the consequences of such a barrier were analyzed by Schwobel and Shipsey [SS66, Sch69]. We also note the lowered activation energy associated with diffusing to the as-cending step from the lower terrace. This activation energy towards the step edge is referred to as the "incorporation barrier" [XLT02], and is indicated as Einc in the figure. 2.4.2 Step-traversal by insertion Another possibility for step-traversal was suggested by Hansen et al. [HSJN91]. The diffusing adatom approaching the step edge is inserted be-hind the edge atom, and hence pushes the current edge atom over to the next lattice site. This process is commonly referred to as "insertion" or "the exchange process". In this case, two atoms are rearranged, instead of just one, as was the case for over-the-edge step-traversal. Little is known about the microscopies of the ES barrier on a polar semiconductor surface, such as GaAs, including whether step-traversal is over-the-edge or via insertion. Pimpinelli and Villain (p. 94 in [PV98]) in-dicate that in the case of (100) face of metals, Au favours the insertion mechanism before over-the-edge traversal. Effective medium theory calcu-lations by Stoltze [Sto94] and direct field-ion microscopy observations by Wang and Ehrlich [WE91]) indicate that the opposite is true for Cu, Ag, Ni, Pd and Pt, as well as for (111) facets. Ab initio calculations by Stumpf and Scheffler [SS94] indicate that for aluminium, the insertion mechanism is important, but also predict that for certain orientations of the steps, there is no ES barrier, "in the sense that the potential is maximum at the center of the terrace rather than on the terrace edge'''' (p. 95 in [PV98]). Itoh et al. [ItoOl, IO00] argue a similar point, that polar semiconductors like GaAs might not have an ES barrier at all. In this thesis, we will not attempt to derive the exact magnitude of this barrier. Instead, we will look at the effects of incorporating such a barrier into analytical models (this Chapter) and computational models (Chapter 4) and compare simulation results of large-scale morphology to experimental data in the form of atomic force microscopy (AFM) images. However, it is important to realize the immediate and very dramatic effect that this barrier has on the morphological evolution of a flat starting surface, even for very small barriers. As an example, in Fig. 2.5 we show the effect on the large-scale morphol-ogy of simulating the surface evolution with an SOS model that incorporates the ES barrier. Part (a) shows that even after a few monolayers of growth, Chapter 2. The crystal growth problem 17 (a) (b) Figure 2.5: Solid-on-solid simulations of growth on a flat starting surface, showing the large-scale morphological dependence on the sign of the ES barrier, (a) E e s =100 meV, total growth is 5 M L , z-range is 16 ML; (b) Ees = —100 meV, total growth is 53 M L , z-range is 3 M L . Parameters used were Esub =1.3 eV, Elat =0.35 eV, T=550°C and F = l M L / s . Substrate size was 300x300 atoms, and red means higher in both topographs. protrusions ( ~ 1 6 M L high) form when a positive ES barrier of 100 meV is employed; these are the initiators of mound-growth. For comparison, part (b) of the Figure shows the effect of growing over 50 M L on a surface where a negative ES barrier of —100 meV is employed. Only a slight deviation from the ideal layer-by-layer growth mode is observed. It is clear that the ES barrier in the common sense as implemented in part (a) destabilizes the evolution of the surface. On the other hand, a negative ES barrier, as seen implemented in part (b), will result in a roughness evolution that is even less pronounced than in the case of kinetic roughening with a neutral ES barrier. The physical meaning of a negative ES barrier is that on GaAs (001), step edges are more attractive to adatoms from the top than they are from the bottom. Due to the complex geometry of reconstructed step edges on GaAs (001) it would not be surprising that intuitive arguments for ES barriers based on SOS models and double height steps might give the wrong sign. In the epitaxial growth experiments done in our laboratory, instabilities like those in Fig. 2.5 (a) have never been observed in GaAs. This important observation could mean one of two things: (1) that the ES barrier is negative, Chapter 2. The crystal growth problem 18 or (2) that there are other mechanisms on the step edge that counteract and stabilize the destabilizing effects of a positive ES barrier. We discuss two such stabilizing mechanisms in the next Section. 2.4.3 Downhill funneling / knock-out process The ES barrier is commonly implemented as a positive barrier that impedes downhill drift, and for many growth systems, especially epitaxy of metal films, it is quite clear that this is a correct assumption (see for example work by Michely et al. [MKC+02] and references therein). It was shown in the previous Section that such a mechanism causes an instability. How-ever, as we will show in Chapter 4, the immediate consequenses of the ES barrier, namely mounding, is not observed in GaAs growth experiments. Crystal growers had been growing smooth GaAs films years before the ES barrier became a theoretical favourite, and additional mechanisms for inter-layer transport were introduced in order to reach agreement between surface growth models and experiments. The two mechanisms for surface relaxation that we will describe in this Section only occur during growth, i.e. to recently deposited atoms. The first mechanism is that of "downhill funneling", introduced by Wolf and Vi l -lain [WV90]. This mechanism describes the local relaxation that a freshly deposited atom can undergo while dissipating its kinetic energy. It is typi-cally implemented as a search for a highest coordination number site within some distance L of the deposition site. Smilauer et al. [vV93, vV95] im-plemented the downhill funneling effect.in combination with a positive ES barrier as part of an SOS model in order to describe RHEED oscillations in (001) GaAs homoepitaxy. In their study, the deposited atoms were allowed to relax within a 7x7-atom area of the drop points. The downhill funneling effect is a smoothing mechanism, just like the inverse ES barrier (i.e. Ees < 0, as in Fig. 2.5 (b)), except that it only applies to atoms that have just dropped from the vapour. Including this relaxation mechanism can cancel the effect of a positive ES barrier during growth, and the net current in such a model can be downhill [KPS93]. It should therefore be noted that the combined effect of a positive ES barrier and a relaxation mechanism could equally well be interpreted using a single mechanism, namely a negative ES barrier. We will show this analytically in the framework of BCF in Section 2.5, where a surface relaxation mechanism is explicitly implemented. Another mechanism for growth induced surface relaxation, often referred to as the "knock-out" [Eva91] process, has an effect on the surface morphol-Chapter 2. The crystal growth problem 19 A Figure 2.6: "Downhill funneling" according to [vV93]. An incoming atom (A) on the upper terrace relaxes to a lower energy site (B) on the lower terrace. t B ogy similar to that of the downhill funneling mechanism described above: an energetic incoming atom that lands close to a step edge can insert be-hind one or more of the atoms closest to the step edge, or "knock them out" as the name indicates, much like the "insertion" mechanism described in Section 2.4.1. The overall effect on the surface morphology of the two mechanisms is similar. 2.4.4 Step traps The potential landscape surrounding a step edge is more complex than the naive diagram in Fig. 2.4 indicates. Some particularily interesting stud-ies were reported by Wang and Ehrlich [WE91, WE93] and Kyuno and Ehrlich [KEOO]. Diffusion of iridium atoms on small islands on the (001) plane was imaged by a field-ion microscope, and the adatoms were found to get temporarily trapped on the upper side of the descending step. This lead to the speculation that for certain surfaces, there might be a poten-tial minimum there. Furthermore, the temporarily trapped atoms at the step were more likely to descend over the step than return to the upper terrace. Based on these observations, the potential diagram in Fig. 2.7 was suggested [WE91, WE93, KEOO]. The reduced Eg in Fig. 2.7 serves to make traversal across the descending step more attractive, as a return to the upper terrace is through the larger barrier ER [Mar94]. Markov [Mar96] suggested that an adatom that got stuck in the trap on the upper side of the step could more easily insert and hence effectively incorporate to the lower step that way, even if the barrier for downhill diffusion might be larger than Esub (note that the barrier labelled by EB here equals Esub + Ees). The potential diagram in Fig. 2.7 might have merit for other surfaces, in addition to the metal surfaces studied in [WE91]. Figure 2.8 (b) shows a top-view A F M image of homoepitaxial growth on GaAs. This surface is prepared with a hydrogen etch and was grown at 590° C for 75 minutes at a Chapter 2. The crystal growth problem 20 Step Trap Surface coordinate Figure 2.7: Potential diagram near a step edge, showing the "step trap", as suggested by Wang and Ehrlich [WE91]. growth rate of l M L / s with an As:Ga ratio of 3.5. The rich morphology of this image requires further study, but it is tempting to point out that the step edges appear white in several places, an indication of additional height. This might be a real effect and hence indicate that even GaAs steps could have potential traps. The white ridges in the A F M image could also simply be an underdamping effect from the oscillating A F M tip. Another interesting feature of this surface is the elongated finger-shapes pointing towards the upper right hand corner of the image. The surface is slightly vicinal with an approximate slope of about 0.05° in the [110] direc-tion, the same direction the fingers are pointing. Figure 2.8 (a) illustrates a possible explanation for these fingers-shapes. Using the terminology from Fig. 2.7, we note that if the potential barrier for. descent, EB, is larger than the incorporation barrier Einc to the step from the lower terrace, the step front will receive more adatoms from the lower terrace than from the up-per terrace. In addition, as the adatom isoconcentration lines (dashed) in Fig. 2.8 (a) indicate, the protruding points on the lower side of the step have a higher chance of being visited by these adatoms and will grow at an average rate that will be larger than the growth rate of the receding parts of the step front, v\. As a result, the already protruding points of the step front will get bigger, eventually causing the fingering that we ob-serve in Fig. 2.8 (b). This effect has been discussed in the literature (the "Bales-Zangwill instability" [BZ90]), and is qualitatively similar to that of dendritic growth. A nice overview of the Bales-Zangwill instability can be Chapter 2. The crystal growth problem 21 Figure 2.8: The Bales-Zangwill instability: part (a) illustrates in top-view the mechanism where protruding points on the step front receives more adatoms from the lower terrace, and therefore grows at a faster rate V2 than the receding points that grow at a rate vi; part (b) shows 2x2/im A F M image of hydrogen-etched GaAs grown for 75 minutes at 590°C. found on p. 170 in [PV98]. Interestingly, the Bales-Zangwill instability is typically associated with a positive ES barrier that allows more adatoms to approach the step from the lower terrace. However, as we pointed out, this effect is also consistent with a negative ES barrier, as long as the incorporation barrier Einc is more easily overcome than EB- The fingering that we see in Fig. 2.8 (b) might also be explained by an anisotropy in the adatom diffusion constant. We will return to this point in Chapter 4 (p. 65), where we show how the As:Ga ratio affects the degree of anisotropy of the growing surface features. The finger-shaped formations seen in Fig. 2.8 (b) may be explained by the potential-energy landscape of the GaAs surface near the step edges, by the anisotropy caused by the As:Ga ratio, or by a combination of the two effects. 2.5 I n te r l aye r t r a n s p o r t i n t he B u r t o n - C a b r e r a - F r a n k m o d e l The interlayer transport mechanisms can be introduced into the BCF frame-work by altering the boundary conditions. The ES barrier can be introduced Chapter 2. The crystal growth problem 22 by noting that it changes the probability of diffusion over the step edge: -DVp(x = 0) = ~^-p(x = ° ) (Ascending step with ES barrier) (2.13) —DVp(x = I) = y—p(x — I) (Descending step with ES barrier) (2.14) where the right-hand side of Eq. 2.14 denotes the number of jump attempts at the descending step edge (~Dp) times a relative rate for crossing the ES barrier, expressed as [PV96] les = a\\ekBT. (2.15) The ES length les is a .measure of the bias that the ES barrier introduces to an adatom placed near the step edge: at a length ay away from the ascending step edge, the adatom behaves as if it were placed at les. When the ES barrier is negative, then les < ay, and vice versa when the ES barrier is positive. In the event that there is no ES barrier associated with the step edge, then the diffusion probability over the step equals that of a diffusion event on the flat terrace away from the edges (les = ay). Solving Eq. 2.6 with the boundary conditions of Eqs. 2.13 and 2.14 gives an adatom current of the form: Fl (an — hs) Ul) = — / 11 e s J (2.16) J W ay (ay + ks + 0 The details of this calculation can be found in Appendix B . l . By allowing the boundaries to be permeable, we see that the ES barrier fully determines the sign of the adatom current: a positive barrier (les > ay) creates an uphill, unstable current and a negative barrier (les < ay) creates a downhill, stable current. Next, we discuss the local relaxation effects from Section 2.4.3. Schinzer et al. [SKR00] incorporated the local relaxation mechanisms described above into the ID B C F formalism by assuming that "there exists an incorporation radius [Rinc] such that all particles arriving close to a downward step within this radius immediately jump down the step edge." This adds to the downhill current, and stabilizes the growth process. Mathematically, little difference can be seen between the effects of the downhill funneling and knock-out mechanisms, and Schinzer's approach can be seen as an implementation of either effect. In order to incorporate Rinc into the B C F model, we split the terrace into two regions: in region 1 (x < I — Rinc) the diffusion equation holds Chapter 2. The crystal growth problem 23 A = U A — I Figure 2.9: Solution to diffusion equation on vicinal terrace with permeable step edges, step edge detachment, and an incorporation radius Rinc, within which incoming atoms fall over the edge in the form of Eq. 2.6, and in region 2 (I — R{nc < x < I) we use Eq. 2.6 with F = 0. This is indicated in Fig. 2.9. The atoms that drop within Rinc of the descending step edge are added directly to the downhill current expression. We refer to the adatom densities in these two regions as p\ and P2, respectively. We also include the detachment at rate K of atoms from the step edges , so the new boundary conditions become: x = 0 : -DVpi = K - —pi (2.17a) a I x = l - R i n c : | v £ l £ L (2.17b) Vp '2 x = I : -DVp2 = —fit - Kp- (2.17c) The solution to the adatom concentration on a step where the incorporation radius is included, is depicted in Fig. 2.9. A full derivation can be found in Appendix B.2. The adatom current off the terrace can be written as: = {-DVPlmup+[-DVp2{l) + ^f) down F [l(a\\ — lea) + Rinc(2les + Rinc)] — * i ? [2.181 0|| 0| | + ks + I and we will discuss this expression in the next Section. Notice that the step edge detachment does not enter into the expression for the surface current. This indicates that the steps release adatoms back onto the terrace, independent of the size of the terraces. Chapter 2. The crystal growth problem 24 A n approximation can be made by noticing that Rinc is usually small (on the order of ay), so by assuming that the solution is pi{x) ^ pi(x) in region 2 (x > I — Rinc)) the current expression including the incorporation radius is equal to the current expression without Rinc, plus a term FRinc/a\\. The effect of including an "incorporation radius" is therefore simply an increase in the surface current, which in turn acts to stabilize the surface growth. This can lead to a "magic slope", that we will describe in the next Section. 2.5.1 The "Magic Slope" We have plotted the current expression from Eq. 2.18 in Fig. 2.10 for various combinations of Ees and Rinc- For a zero incorporation radius and a positive ES barrier, we see that the current is negative for all values of the terrace width. However, when a positive incorporation radius is included, the curve crosses zero. This happens at a terrace width, I*, that is given by ^* Rjncj^es ~f" Rinc) ^ ^Q^ les ^|| A positive "magic-length" can only happen when the ES barrier is positive, so that les > This form for the adatom current is interesting in that it causes the surface to stabilize at a fixed ("magic") slope: terraces of width I < I* will get longer by the downhill current and terraces wider than I* will get shorter through the uphill current. Eventually, all terraces will have a width I*. Note that without the incorporation radius, there is no stable slope. Thus, two opposing mechanisms are necessary for the current to have a zero-point; in this case, in the presence of a stabilizing incorporation radius, the ES barrier must be positive. "Magic slopes" have not been observed on GaAs, however, they can easily be identified in homoepitaxy of Pt on P t ( l l l ) [MKC+02, KPVC00] , where the emergence of pyramid-like mounds occurs from growth on initially flat surfaces, see Fig. 2.11. Similar work on many other metal systems shows similar results, for example Cu(100) [EFFL94, ZW97], Ge(001) [NCH+95], Fe(001) [SPSZ95, TKWR95] and Rh(001) [TWUC96]. The "mounds" in Fig. 2.11 have straight, linear slopes, and for that reason we will refer to them as "pyramids". These pyramids therefore appear to be inversion-symmetric under h —• —h. However, closer inspection of the images in Fig. 2.11 (c and d) reveals that the top layers of the pyramids are flat and become relatively large before further nucleation of new islands occur on them. On the other hand, the bottoms appear to be sharp or cuspy, a pinning effect that we will discuss further in Section 3.2.2 (p. 31). Chapter 2. The crystal growth problem 25 1 • . _ E >0, R. =0 -es inc E >0, R =a„ es inc E <0, R. =0 — — es inc E <0, R. =a„ es inc • ^ i 1 • • • 4 6 8 Terrace width, I (nm) 10 12 Figure 2.10: Adatom current as function of terrace width. The ES barrier and the incor-poration radius is varied. The current crosses zero only in the case when Ees > 0 and Rinc > 0, and the zero-current point is indicated: terraces shorter than Z* will grow longer and terraces longer than I* will shorten. The inversion symmetry is therefore not complete, even for metal epitaxy and positive ES barriers. These pyramid-shapes differ markedly from the typical mounds that we observe in GaAs homoepitaxy both in shape and temporal evolution [BRA+01, NCC+96]. The "magic slope" has been investigated by, among others, Schinzer et al. [SKR00]. Johnson et al. [JLB+96] did similar work, where a stable-current expression was derived by customizing the surface adatom current. By stabilizing the current, the uncontrollable growth due to the instabilities induced by the ES barrier are avoided. The magic slope therefore serves two purposes: (1) on one hand, it allows pyramid-growth on initially flat surfaces, and (2) it prevents such perturbations in the surface height from diverging. 2.6 A d a t o m d e s o r p t i o n It is typically assumed that desorption of adatoms is negligable in the case of GaAs homoepitaxy at standard growth temperatures (~600°C), as the rate of incoming atoms at typical growth fluxes is high enough to bury an adatom before it has a chance to escape back into the vapour. However, when the flux decreases and we anneal the surface, this situation changes. We will Chapter 2. The crystal growth problem 26 >>>^> N ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ 4 i f Figure 2.11: Growth of Pt on P t ( l l l ) at T = 440 K. The scanning tunneling microscopy topographs show surfaces after (a) 0.3 M L , (b) 3 ML, (c) 12 M L , (d) 90 ML. The scan size in (a)-(d) is 2590x3450 A 2 and the growth rate is F = 7 x 1 0 " 3 M L s _ 1 . From T. Michely [MKC+02]. implement desorption in this Section with annealing or low-flux growth in mind and see how it affects the adatom current. We incorporate desorption into the BCF model by DV2p{x,t) = - - - r ^ 9 - (2.20) where Tdes is the average time between desorption events. We include step edge detachment as before and solve Eq. 2.20 with the boundary conditions defined by Eq. 2.17. However, in order to reduce the complexity of the solution somewhat, we make use of the fact that Rinc ~ that is, Rinc is small compared to the typical terrace width in the weak surface texture limit. This means that we can approximate pi by p\ in region 2 of Fig. 2.9, but still add the downhill funneling to the current term. The expression for Chapter 2. The crystal growth problem 27 the adatom current is j (0 = • I H « \\)\r^ ajj_ where A = e l / V D T d e s (2.22) The details of the calculation can be found in Appendix B.3. The term ( A 2 — 1) in Eq. 2.21 is always positive, and hence does not contribute to the overall net direction of the surface adatom current. The denominator in Eq. 2.21 is also always positive. The term (les — ay) changes sign with the ES barrier, as before. The last factor in the numerator of Eq. 2.21 is negative for Taes > a^K/DF, which it is in the case of GaAs growth, so the expression for the current that we have derived is positive (stable) during growth. However, when the flux goes to zero, and we anneal the surface, desorption could have a destabilizing effect on the surface morphology. 2.7 A sur face g r o w t h m o d e l d e r i v e d i n t h e B u r t o n - C a b r e r a - F r a n k f r a m e w o r k Earlier in this Chapter, we derived expressions for the adatom current in the B C F framework (see for example Eqs. 2.12, 2.16, 2.18 and 2.21) of the form J(0 = § | (2.23) where P and Q are polynomials in I of order a and b, respectively, with prefactors that are determined by the growth parameters F, T, Ees, and Rinc, but not the diffusion constant D or the step edge release rate K. Noting that the local surface derivative m can be written as m± = ±a±/l (for V / i sj; 0), we can rewrite j(Z) as j ( m ± ) , and obtain a surface growth equation in terms of the surface derivatives. We Taylor expand j ( m ± ) in order to get the leading terms in m±. Using the expression Eq. 2.18 for the current, where all the physics of B C F was involved except the adatom desorption, we expand j(m+) for positive slopes and get . , = ^ = F ^ - U + ^ ( a i M 1 - ^ ( 2 . 2 4 ) • I  i = i ^ I  ' Chapter 2. The crystal growth problem 28 where V = 4 ( 4 - 4 + Rinc(2les + Rinc) ) (2.25) QJ I Replacing m+ with V / i , we can derive the growth equation directly from this expression as dth = Vx - j(V/i) = Vzh x 2u„dj(Vh) d(Vh) i=0 1 + (a.i\+les) \Vh\)' (2.26) (2.27) where in Eq. 2.27 we let Vh —> | V / i | in order to account for the negative slope expansion m_ = —a±/l. This closed form equation is valid as long as |V7i| < — ^ y - (2.28) a\\ + les Inserting the expression for les from Eq. 2.15, we note that the slope re-striction for a zero ES barrier is V / i < 1/2, or surfaces less steep than about 26°. A similar form for the growth equation was proposed by John-son et al. [JOH +94]. In their case, the ES barrier was assumed positive, and the incorporation radius was not considered, so the prefactor u turned out negative. Their growth equation was therefore unstable, and a stabilizing "diffusional" term — K V 4 / I was added, without which the equation would be ill-behaved (see discussion on p. 231 in Ref. [BS95]). By truncating Eq. 2.26 at i = 2, a more common form for the continuum growth equations is revealed: dth ~ uV2h + i / i | V / i | V 2 / i + z / 2 (V/ i ) 2 V 2 / i , (2.29) a = 1,2,3,... (2.30) Remarkably, Eq. 2.29 shows that to the lowest order in V / i , we have re-derived the Edwards-Wilkinson (EW) equation. A close inspection of the prefactor u (Eq. 2.25) shows that for a zero incorporation radius (Rinc — 0), Chapter 2. The crystal growth problem 29 the sign of v is completely determined by the ES barrier: a positive barrier makes les > oy and hence leads to a negative u, and vice-versa for a negative barrier. It also shows that when we set the ES barrier to zero, i.e. when interlayer diffusion is unbiased, the presence of the incorporation radius indeed stabilizes the surface, as u > 0 whenever Rinc > 0. The third term of Eq. 2.29 can be recognized from Ref. [SK94], and is equivalent in form to V(V7i) 3 . 2.8 S u m m a r y In this Chapter, we reviewed some basic concepts in crystal growth dy-namics, including mechanisms that have a strong influence on the evolution of the surface morphology, such as adatom diffusion, the ES barrier and downhill funneling. These mechanisms were introduced into two growth models, specifically the analytical B C F theory and a stochastic SOS model. Through analysis of these models, it was found that a positive ES barrier causes mounding on initially flat surfaces. Conversely, a negative barrier was found to smoothen initially rough surfaces. It was also argued that the com-bined effect of a positive, destabilizing ES barrier and a surface relaxation mechanism could be equally well explained by a negative ES barrier. Some consequences of various interlayer transport mechanisms were dis-cussed. As an example, the "magic slope" arose when competing uphill and downhill adatom drift mechanisms were included in the growth model. This was found to give rise to slope-limited pyramid-shapes in the evolving sur-face morphology. We also discovered that the potential-energy landscape of GaAs near step edges can give rise to finger-shaped protrusions in the direction perpendicular to the step-direction. This effect was also consistent with the ES barrier being negative. The B C F formalism assumes that the motion of step edges is much slower than diffusion of adatoms between the steps. This way, the adatom concen-tration on the terraces can be assumed to be in steady-state. Furthermore, it does not include the nucleation of islands on terraces, a limitation that we will investigate further in Chapter 3. Although simple, the B C F theory proves to be a remarkably useful formalism. Based on surface adatom cur-rent expressions calculated in this Chapter, we derived a continuum growth equation that to the lowest order coincides with the E W model. This deriva-tion required that the net surface adatom current was downhill. Chapter 3. Conventional continuum modeling of homoepitaxial growth 30 Chapter 3 Conventional continuum modeling of homoepitaxial growth A particularily useful approach for modeling the surface morphology during epitaxial growth is to use continuum growth equations of the type that we referred to in Chapter 2 (see for example Eqs. 2.27 and 2.29 of Section 2.7). These equations use a coarse-grained version of the surface, and hence ex-hibit a lower spatial resolution than for example SOS models. However, their relative calculational simplicity allow them to complete in minutes, rather than hours or days which is the case for microscopic models like SOS. Also, they allow for an understanding of the effects of microscopic mechanisms on the macroscopic morphology. Many books and review articles are available for a more complete introduction to the topic of continuum models, see for example Refs. [Vil91, BS95, HHZ95, PV98] and references therein. 3.1 C o n t i n u u m g r o w t h equa t i ons In a continuum approximation, the discrete nature of the atoms is averaged, and the vertical steps are coarse-grained into a smooth surface, whose slope determines the density of incorporation sites (or steps) as S oc Vh. In this approach, the time rate of change in the surface height h(x, t) is expressed in terms of derivatives of h [Vil91, BS95, HHZ95, Kru97, ZL99], of the form dth oc V2h + V 4 / i + (V/i) 2 + V 2 ( V / i ) 2 + V ( V / i ) 3 + ... + F + r?(x,t) (3.1) The various terms on the right hand side are preceeded by constants (left out of Eq. 3.1) that further reflect the underlying physical processes dur-ing growth. The source term ?7(x, t) simulates the random arrival of atoms ("noise") at an average rate F. Generally, the noise is assumed to be uncor-related, with (?7(x,i)77(x',t')) = D<5(x — x')<5(£ — t'), where (...) represents an average over possible noise configurations. It is common to write the growth Chapter 3. Conventional continuum modeling of homoepitaxial growth 31 equations in a moving reference frame, by the transformation h —> h — Ft, which removes F from the growth equation. We will adopt this convention here. In this general growth equation, symmetry arguments account for the absence of certain terms on the right hand side, like the direct dependence on h, x, and t, as well as V n / i and (V/Y) n for n odd. For a more in-depth discussion on symmetry, see for example Ref. [BS95]. Certain symmetries will be broken in some growth systems. For example, sloped parts of the surface grow faster laterally than vertically, causing a mounding effect, and this asymmetry in the growth direction is a sign of the presence of nonlinear terms in the growth equations that inherently break the up-down symmetry h-+-h. Continuum equations in the context of crystal growth have received much attention in the literature, especially in the last decade. The emphasis has been on the understanding of universality classes of various growth models when applied to perfectly flat starting surfaces. Less effort has been put towards applying such, models to actual surface growth experiments, with realistic starting surface shapes. In fact, investigations into various growth equations typically do not even have a figure showing a surface. Obviously, stronger communication between crystal growers and surface theorists will lead to better growth models. 3.2 L i n e a r g r o w t h equa t i ons The linear part of a growth equation consists of even derivative terms of the form dth = uV2h - KV% + C6V6h + ... + T} . (3.2) where the equation is written in such a way that positive prefactors u, K and Ce lead to a stable growth equation. For small surface slopes (long wavelength surface height perturbations), terms of order 6 and higher are very weak and their impact on the surface evolution harder to verify, and therefore they are often not included. The stability of a linear growth equa-tion can easily be investigated by taking the spatial Fourier transform of the deterministic part of Eq. 3.2: = -(vq2 + Kq4)H (3.3) Chapter 3. Conventional continuum modeling of homoepitaxial growth 32 where 7Y=7i(q, t) and q is the spatial frequency. The straightforward solu-tion is dtH(q, t) = e'^+^niq, 0) (3.4) As long as the prefactors are all positive, the evolution of each Fourier component is stable. If the highest order prefactor (in this case K) is nega-tive, then the q4 term dominates, and the higher-order Fourier components of the surface diverge. A special case occurs when the lower order term is unstable, but the higher order term stabilizes it, as is the case in the Kuramoto-Sivashinsky (KS) equation [Kur94, Siv77, SM80] dth = -uV2h - K V 4 / I + V (3-5) In this case, the lowest order term — uV2h grows when u > 0, but is stabilized by the higher order term — K V 4 / I . This results in the emergence of a fixed wavelength Lc = 2ir/qc with spatial frequency qc = y/v/2n [PT00]. Higher order analogies can be made. For example, a stable sixth order term can stabilize an unstable fourth order term with a critical qc = ^2K/3CQ, or even an unstable second order term —vV2h, with a critical qc = {y/3Ce)1/4. The KS equation was used in models of turbulence in chemical reactions and in fluid mechanics. It has been used to describe for instance flame fronts and unstable flow of a viscous fluid on a vertical plane [FBSK04]. Barabasi and Stanley [BS95] showed that the KS equation could explain the pattern formation that occur during sputtering. However, it does not apply to the case of GaAs homoepitaxy, since the surface amplitude decays with time and there is no fixed wavelength. 3.2.1 The Edwards-Wilkinson equation The lowest order linear equation that could be expected to describe surface growth is the E W equation [EW82]: b\h(x, t) = i/V2fc(x, t) + ??(x, t) (3.6) The linear term can arise from evaporation/condensation dynamics, as sug-gested by Mullins [Mul59] who described the interaction between a vapour with a constant chemical potential \iv and a surface with an average chemi-cal potential \XQ < u.v. In the continuum approximation, (XQ is determined by the surface tension os (free energy per unit surface area). To lowest order, variations about are induced by the local surface curvature [Kru97] as /z(x) = Mo - cr s V 2 / i (x ) (3.7) Chapter 3. Conventional continuum modeling of homoepitaxial growth 33 The incorporation rate can be determined using the ansatz: dth = — T/j, for the relaxation, where T > 0 is the interface mobility. Inserting this expres-sion for u. (Eq. 3.7) and subtracting the constant growth rate T/J,Q yields the E W equation (Eq. 3.6), with v = asT, where v is always positive. The E W equation is stable, and any perturbations on the starting surface will smooth out towards a long-time steady-state level dictated by kinetic roughening, i.e. a balance between the smoothing effects of V 2 / i and the fluctuations due to the growth noise 77. The deterministic version of the E W equation is easy to solve, and gives a monotonic decay of each Fourier coefficient of the surface according to exp(—vq2t), where q = 2-K/L for surface features of lateral size L. 3.2.2 The Ehrlich-Schwobel barrier in the continuum limit In an influential paper, Villain [Vil91] argues that the ES barrier can lead to the E W equation. This argument contradicts the results of the simulations that we showed in Figure. 2.5, where the conventional ES barrier caused mounding on a flat starting surface. His argument considers the case of a slightly miscut (vicinal) surface. Intuitively, in the presence of "diffusion bias" (interlayer transport in the uphill direction) caused by the ES barrier at the step edge, a stepped surface with terraces of different lengths would even out, and successive approximations result in the linear smoothing term vV2h with v > 0. The argument in Ref. [Vil91] can be summarized as follows: assume an infinite ES barrier at the descending step, and assume a stepped surface where the terraces alternate between being wide and narrow, and where no island nucleation is allowed. The wide terraces will receive more atoms due to their larger areas than the narrow terraces, and all of these atoms will eventually attach to the ascending step. The ascending steps on a wide terrace will therefore move faster than the ascending step on a narrow terrace, thus the terrace widths equalize. This Gedanken-experiment is verified by a ID calculation in Fig. 3.1, where the artificially created rough and vicinal starting surface straightens out during growth. The growth-mode reached when the terrace widths equalize is called step-flow. Villain continues his discussion of growth equations by deriving the E W term for growth on a singular surface, but this time the conclusion is that the prefactor v is negative (see Appendix C in Ref. [Vil91]). This contradiction is not resolved in the paper by Villain, although a consensus seems to be that v must be negative due to the commonly accepted destabilizing effect of the ES barrier. Villain estimated that the size of the ES barrier can be no more than a room temperature value of fcgT (~26 meV), and at the Chapter 3. Conventional continuum modeling of homoepitaxial growth 34 i , , , Terrace length —h—h—'l ' l 'l 1 1 1 l i l l — i i i — i - » I I 1 ¥ 1 1 i 1 i 1 i 1 i ' i 1 i 1 i 1 — H — i 1 1 ' • 1 ' i 1 ' i 1 ' • • ' / — / — Substrate ^—h—<—,—1 H H Figure 3.1: Effect of positive ES barriers on ID vicinal surfaces according to Villain [Vil91]. Each terrace receives atoms at a rate proportional to its width. An infinitely strong ES barrier is considered, so all adatoms stay on the terrace where they land. The equivalent of 10 M L was deposited in this calculation at a growth rate F = l M L / s . The boxes represent the material added each 0.5 s. Hence, the wider terraces have larger boxes, as they receive more material per unit time. time of the paper (1991), "no such instability has been reported, presumably because experimentalists just try to avoid it" so that at common epitaxial growth temperatures (~500-600°C), the effect of the ES barrier is "probably negligible". In the following, we will show where it is that Villain's argument for a positive ES barrier's stabilizing effect on a vicinal surface breaks down. We already showed in Fig. 2.5 that introducing a positive ES barrier into an SOS simulation triggered instabilities in the evolution of the morphology of a singular, flat starting surface. Fig. 3.2 shows a similar SOS calculation, but this time using a vicinal starting surface. The results indicate that the ID calculation in Fig. 3.1 must be ignoring important physics. Fig. 3.2 (a) shows a starting surfaces that we evolve in Fig. 3.2 (b) and (c) for ES barriers with opposite sign. The starting surface was created so as to include both large flat terraces where nucleation is likely to occur, as well as steep sections where step-flow is the preferred growth mode. These simulations show that it is a negative ES barrier that stabilizes vicinal surfaces, and not a positive one, in contrast to Villain's argumentation. Figure 3.2 (b) shows the evolution of this starting surface using a negative ES barrier. The preferred growth mode is by step-flow, and there is only a minor indication of 2D growth on the wider terraces on top and bottom of this image. The negative ES barrier is stabilizing, and the initially very wide terraces at the top and bottom of the surface have shortened and the shorter terraces in the middle have grown wider. Chapter 3. Conventional continuum modeling of homoepitaxial growth 35 (a) (c) (d) Surface height Figure 3.2: SOS simulations showing the effect of the ES barrier on the surface morphol-ogy of a vicinal surface. Figures (a) through (c) show 400x200 atom topographs of (a) the starting surface; (b) 28 ML growth with E e s=—50meV and (c) 20 M L growth with E e s =50meV (z-range 48 ML). Figure (d) shows, from left to right, crossections down the middle of the images (a)-(c). Apart from the ES barrier, the standard parameters were used (Table 2.1). The growth temperature was 600°C. Figure 3.2 (c) shows the corresponding simulation for a positive ES bar-rier. In this image, there is enhanced nucleation on the wider terraces when compared with the simulation in part (b), as one might expect. What might be surprising here is that large mounds form on the steep parts of the sur-face. This is in direct disagreement with Villain's ID thought experiment (Fig. 3.1). The adatom concentration, n, on the narrow terraces in the steep regions will increase due to the blockage at the descending steps. This in turn increases the probability of dimer nucleation, which goes as ~ n 2 . It is evident that even a few occurrences of nucleation on the narrow terraces is enough to destabilize this surface and lead to 3D growth. While it is the 2D effect of nucleation of small dimer-islands that triggers this effect, it is the formation of "advacancies" between the small islands and the step edges of the macroscopic morphology that causes the dramatic 2D to 3D transition. Figure 3.2 (d) shows, from left to right, scanlines down the middle of the im-ages (a)-(c). Again, it is quite clear that a negative ES barrier (middle line) preserves the step-flow growth mode, while the positive ES barrier (right line) causes mounding in the steeper section of the surface. As the terraces are already narrow, small advacancies (or trenches) will quickly develop when dimers nucleate. A conventional ES barrier inhibits Chapter 3. Conventional continuum modeling of homoepitaxial growth 36 downhill flow of adatoms, and will therefore prevent adatoms from entering the advacancies. In fact, with a large ES barrier, a single advacancy will fill in at an average rate of F. This means that it obeys the statistics of, and evolve similarily to, a surface subject to a random deposition model without surface diffusion. Surface sites away from advacancies are bombarded by the same amount of atoms from the vapour, but also get visited by diffusing adatoms. The added visits from diffusing adatoms allow regular surface sites to fill in at a faster rate than the advacancies, or vice-versa, the advacancies can be seen as pinning centres for the surface height. Recent theoretical work recognizes the instability in step-flow growth on vicinal surfaces due to large, positive ES barrier [RvK96, KK04]. It has also been suggested that the instability on vicinal surfaces is a realization of the Bales-Zangwill mechanism1, which causes the steps to "meander" [KKK02]. We summarize this Section as follows: (1) a positive ES barrier is desta-bilizing both to singular (Fig. 2.5) and vicinal surfaces (Fig. 3.2) and cause an obvious emergence of 3D growth. Conversely, a negative ES barrier is stabilizing in both cases, and 3D growth is not observed; (2) Villain's simple thought experiments [Vil91] in ID ignore nucleation of dimer islands, and the resulting morphology from his model is in stark contradiction with the 2D analysis done in this Section. Based on the arguments made in this Sec-tion, the unstable E W equation (pt = — |^|V 2 / i ) is now consistently linked to a positive ES barrier, whether it is applied in growth on vicinal or singular surfaces. 3.2.3 Mull ins ' diffusion The next possible linear growth term is the diffusion term attributed to Mullins [Mul59]: dth = -KV2(V2h) (3.8) or simply —K,V4h. The formulation is based on the mass conservation argu-ment, where it must be possible to write the growth equation as a continuity equation b\h = -V-i (3.9) The current density j is a vector parallel to the average surface direction (and not the local surface). Landau & Lifshitz [LL67] showed that j = Const, x Vfi(r, t) (3.10) 1 J . Krug, private communication. Chapter 3. Conventional continuum modeling of homoepitaxial growth 37 o1.2 Q. E CO t CO L2 (b) 1*0 Surf, coordinate (arb. units) Time (arb. units) Figure 3.3: Mullins' term for surface diffusion causes surface peak-to-valley amplitude overshoot, (a) shows the square starting surface (at to) along with linescans of the surface at two later times (ti and t^)- The average growth rate F is subtracted, (b) shows the peak-to-peak amplitude of the surface. Equation (3.9) and (3.10), together with the expression in Eq. 3.7 for the chemical potential, lead to Mullins' diffusion term, Eq. 3.8. A n interesting consequence of Mullins' diffusion term is that for certain starting surfaces, namely those with steep sections and hence significant higher-order Fourier components, it causes the initial peak-to-peak ampli-tude to temporarily overshoot, as indicated by Fig. 3.3. This might seem surprising, as Mullins' equation (3.8) should have its Fourier coefficients decay in a fashion similar to the E W equation, namely as exp(—KqH). How-ever, Fig. 3.3 (a) shows the evolution of a square grating surface (at time to) during growth, according to Mullins' equation. Soon after growth com-mences, holes are formed close to the sidewalls in the valleys, as well as bumps close to the top of the grating. This causes the early increase in the surface amplitude, labelled t\ in Fig. 3.3 (b). As these valley-holes move slowly towards the centre of the valley, they temporarily start to level out, but then superimpose to form one deeper cusp (at time t2), which causes the second rise in the surface amplitude. The bumps on the top of the grating do the same in this symmetric problem. Once the bumps come together, the valley fills monotonically. Interestingly, this form of amplitude overshoot (bimodal, non-monotonic) is observed in a later Chapter (Ch. 5), where the starting surfaces are pre-patterned lithographically with shapes similar to those depicted in Fig. 3.3 (a). 3.3 N o n l i n e a r g r o w t h e q u a t i o n s Linear growth equations like Eq. 3.2 are analytically solvable in the deter-ministic case, but the inclusion of noise or nonlinear terms requires the use of Chapter 3. Conventional continuum modeling of homoepitaxial growth 38 numerical solution techniques. The complex properties of even the simplest nonlinear growth model, the Kardar-Parisi-Zhang (KPZ) equation [KPZ86], continue to fascinate researchers, and new publications on this equation it-self as well as its applicability to crystal growth occur frequently. In this Section, we describe some of the simplest nonlinear growth terms. The ques-tions we seek to answer include (a) their effects on the surface morphology, and (b) their physical origins. 3.3.1 The Kardar-Parisi-Zhang equation One of the most studied continuum equations is the K P Z equation [KPZ86, Yak81]: dth = uV2h + | (V/ i ) 2 + 77 (3.11) The nonlinear term is an approximation which is often used to model the sit-uation where growth proceeds outward from the local surface normal [BS95, Kru97], although, other interpretations are possible. Villain [Vil91] showed that growth on a stepped surface where desorption was' allowed could lead to the K P Z nonlinearity. It has also been shown that the nonlinear term comes out naturally when overhangs and vacancies are allowed in the growth model [Bal98], for example in models like the ballistic deposition model. Overhangs and vacancies are not" observed in M B E , and in GaAs homoepi-taxy, neither is desorption. A common objection to the K P Z nonlinearity in the context of M B E is that it does not conserve mass. This is a valid concern in the case of GaAs homoepitaxy, although for the relatively flat surfaces that we treat in this Chapter, the added material due to the K P Z nonlinearity has a minimal influence on the average growth rate. The nonlinear term also breaks inversion symmetry in the growth di-rection: positive values of A cause round mounds on the tops with sharper cusps in the valleys, and vice versa for negative values. Intuitively, one would assume that the mounded shapes are related to the added growth on the sloped parts of the surface, where the incorporation of adatoms to steps takes place, as the steps are in abundance there. In the case of homoepitaxy on GaAs (001), the presence of a second order linear term means that there is a downhill current of adatoms, caused for example by a negative ES barrier, as we discussed in Section 3.2.2. Even if the ES barrier is positive, the nonlinear term is capable of balancing the instability introduced by a negative E W term (u < 0) [KPS93], but it is not known whether weaker nonlinearities that conserve mass (i.e. V 2 ( V / i ) 2 ) have this ability [SKJ+92]. Chapter 3. Conventional continuum modeling of homoepitaxial growth 39 3.3.2 The molecular beam epitaxy equation Under typical GaAs M B E growth conditions, atoms arriving at the surface undergo significant diffusion before incorporating at favourable sites. If the sticking coefficient is unity, then the growth rate is determined solely by the flux and not by the surface shape. Therefore, nonconservative terms should not be present in the growth equation. To this end, the M B E equation has been formulated [SGG89, WV90, LS91, Sar93, SK94, BS95, Kru97] where r\ may be either nonconservative flux noise like we described in Sec-tion 3.1, or conservative diffusive noise with correlation function [Kru97] (rj(x,t)rj(x!,t/)) = D c V 2 £ ( x - x.')8(t - t'). The fourth order linear term is the same as Mullins' diffusion term [Kru97, Mul57]. The fourth-order non-linear term is generally believed to be generated by a nonequilibrium surface chemical potential induced by the deposition flux [Kru97]. In this case, one expects the chemical potential to be an even expansion of V / i , of the type //(x) oc (V/T.)2 + (V/I) 4+. •., instead of the linear form in Eq. 3.7. A discussion on this expansion can be found in Ref. [Kru97], pp.245-252. This term can also be attributed to a slope dependent adatom concentration [Vil91]. For a highly symmetric orientation, Villain proposed an adatom concentration dependency of the form p = po + pi(V/Y) 2 , from which the nonlinear term emerges through the application of Fick's Law, Eq. 2.4, and a continuity equation, like Eq. 3.9 [Vil91]. We will refer to the fourth-order nonlinearity as the "Villain-term". 3.3.3 Conservative Kardar-Parisi-Zhang equation A conservative version of the K P Z equation can be generated by modifying the nonlinear term. To achieve this, we first note that when the growth proceeds outward from the surface normal at a constant rate A, the growth rate can be projected onto the vertical direction, leading to demonstrating the link between normal growth and the K P Z term. The approximation is valid for | V / i | -C 1. For a rough surface of side length L, Eq. (3.13) leads to the addition of a total volume: dth = - f i V 4 / i - f V 2 ( V / i ) 2 + 77 (3.12) (3.13) 8V (3.14) Chapter 3. Conventional continuum modeling of homoepitaxial growth 40 in an infinitesimal time 8t, where (...) represents an averaging over the sur-face. On a perfectly flat surface with the same dimensions, the total volume added is simply 5VQ = XStL2. A conservative growth equation is obtained by normalizing the volume added, then subtracting the growth rate. This is achieved through the substitution ^ ( V r * ) 2 - A (3.15) where we have used the approximation in Eq. (3.13). Physically, in the case of GaAs growth discussed here, the correction describes the situation where the density of mobile atoms on the surface adjusts itself such that the net incorporation rate balances the arrival rate of atoms from the growth flux. The resulting equation is non-local and is difficult to investigate using stan-dard analytical techniques, so we have simulated the behavior numerically. One could also analyze the K P Z equation in its original form (Eq. 3.11) by viewing it as a lowest order approximation to an actual, conservative term. Our particular treatment results in a non-local equation, which might very well belong to a different universality class than the K P Z equation itself (see Section 3.5). 3.3.4 The Das Sarma equation The last term that we will discuss is the conservative, nonlinear term A i i 3 V ( V / i ) 3 = 3Ai i 3 (Vfr) 2 V 2 rt . Das Sarma and Kotlyar [SK94] discovered that in the absence of the E W term V2h, this term will dominate the asymp-totic, long-time scaling behaviour of a general growth equation2 of the type dth = -nVAh + ^V2(Vh)2 + A 1 ; 3 V ( V / i ) 3 + T? (3.16) We have already covered the first two terms on the right hand side, so in Fig. 3.4 we show the effect of the third term A i ) 3 V ( V r i ) 3 on a square starting surface, with A i ) 3 > 0. We see that this term ensures smooth surface evolution: in the case that the E W term vV2h is zero, the surface still smooths. No physical mechanism is.known to generate this term in the absence of the V2h term [LS91]. 2See Section 3.6. Chapter 3. Conventional continuum modeling of homoepitaxial growth 41 Surf, coordinate (arb. units) 3.4 E f fec t o f g r o w t h e q u a t i o n t e r m s o n sur face s h a p e e v o l u t i o n A n educational way of looking at the various linear and nonlinear growth terms is to evolve them one time step and see what effect they have on the surface evolution. In order to achieve this, we have included the often-used illustration [BS95, SK94] depicted in Fig. 3.5. In part (a), a typical surface segment is shown. Parts (b) to (f) show the various growth terms discussed in this Section when applied to the surface in (a). The nonlinear K P Z term in (b) adds to the steepest parts of the surface, but it is always positive, and hence nonconservative. The E W term in (c) smooths the surface by removing material from points and adding it to kinks. The Das Sarma term in (f) is seen to act similarly to the E W term, but at a slightly shorter length scale, or higher spatial frequency, as expected from the higher-order of derivatives present in this term. Mullins' term in (d), with a negative pre-factor, adds material to kinks by drawing from nearby regions. The Villain term in (e) conservatively draws adatoms from kinks and points and adds it to the steeper parts of the surface. 3.5 U n i v e r s a l i t y c lasses The different kinds of growth processes, as represented by the different growth equations, can be characterized by their universality class. Asso-ciated with each universality class is a set of exponents a (roughness expo-nent), j3 (growth exponent), and z = a/(3 (dynamic exponent), describing the scaling properties of the surface. These scaling parameters were intro-duced in the previous Section. We give a quick overview of the theory of scaling in this Section [Bal98]. c t=0 Figure 3.4: Das Sarma's term causes no surface amplitude overshoot. The square starting surface evolves into the succes-sively smoother surfaces. The average growth rate F is subtracted. At is an arbitrary time-increment. Chapter 3. Conventional continuum modeling of homoepitaxial growth 42 (d) (e) (f) Figure 3.5: Effect of linear and nonlinear growth equation terms on surface evolution, (a) Typical surface segment; (b)-(f) various derivatives of (a) as indicated. Al l x-axes span equal ranges; y-axes are normalized individually. The interface width, W, of a continuous surface is defined by , W2 = y f dx (h(x) - hf (3.17) L JL where h is the average'height of the surface, and equals Ft. This quantity gives us a statistical measurement of the surface roughness. Many studies of thin film growth are based on the dynamic scaling behaviour of the surface statistics that depends on the various surface relaxation and roughening processes [FV85, FV91]. By plotting the interface width as a function of time, two distinct regions are observed: at times less than some 'crossover' time tx, W increases as a power of time [BS95]: WfL,t)~tP, [t4Ztx] (3.18) where (3 characterizes the time-dependence of the growth dynamics. This interface width continues to increase until a certain value is reached: W saturates according to Wsat(L)~La, (3.19) The crossover time depends on the system size as tx ~ Lz, [t.~tx\ (3.20) Chapter 3. Conventional continuum modeling of homoepitaxial growth 43 Table 3.1: Critical exponents associated with various growth equations in (2+1) dimen-sions. The second column indicates whether the noise is conservative (C), or nonconser-vative (N). Equation Noise ot. P z E W N 0 0 2 K P Z N 0.385 0.24 1.58 M B E N 0.667 0.2 3.333 M B E C 0 0 4 The three parameters, a, (3 and z are dependent on one another, and they can be extracted from a given system by normalizing W and t [FV85]. We obtain the relation W ™ • / ( * • ) (3.21) Wsat{L) Substituting from above we find that for an initially smooth surface, the in-terface width scales according to the so called Family- Vicsek scaling relation [BS95, FV85, FA85]: W(L, t) = Laf (t/L"'^ (3.22a) / ( x < l ) ~ ^ (3.22b) f(x > 1) -> constant (3.22c) By plotting this relationship, we see that systems of different sizes L will collapse onto one curve. Solving for z around t = tx we get This scaling law links the three exponents, and is valid for any growth system that follows the Family-Vicsek scaling relation. Table 3.1 summarizes the known exponents for some of the growth equations [BS95]. Experimental measurements of the surface shape can be used to extract the growth exponents, thereby giving an indication of the particular equation governing the surface evolution. One useful measure is the power spectral density (PSD), defined by PSD(q,t) = [ft(g,i)]2 (3.24) Chapter 3. Conventional continuum modeling of homoepitaxial growth 44 where H(q, t) is the Fourier transform of /i(x, t) at spatial frequency q. After a long period of growth the PSD is expected to approach an asymptotic form determined by the critical exponents of the growth equation. In (2+1) dimensions the asymptotic form is PSD(t ^ oo) ~ q - 2 ( 1 + a ) (3.25) where q = |q|, and the overall amplitude depends on both the strength of the noise and the coefficients in the growth equation. For finite growth times, the asymptote will be observed only for q larger than some cutoff value qc ~t~l/z. Neither the PSD nor the interface width W contain any phase informa-tion. This means that W does not distinguish between a surface that has mounds in the up-direction from one that has them in the down-direction. In Fig. 3.5, we showed that certain terms in a general growth equation expan-sion attribute to a breach in the inversion symmetry for the surface height, h. The prefactor of such terms determines whether the mounds point up or down. This means that the sign of the prefactors for the nonlinear terms cannot be determined by analyzing the interface width or PSD alone. 3.6 P o w e r c o u n t i n g a n d re levance o f t e r m s By scaling the growth equation in question, we can determine the relevance of the various terms in the long-time, long-wavelength limit ("the hydro-dynamic limit"). We scale the relevant variables according to the following transformation: h = bah(bx,bzt) (3.26) In the new coordinates, the interface velocity can be written as: b\h = ba~zdth (3.27) The E W term scales as: V2h = ba-2V2h (3.28) We can therefore write the scaled E W equation (3.6) as: ba~zdth = vba-2V2h + b-(d+z)/2r) (3.29) 8th = vbz~2S72h + b^-d)l2-ar) (3.30) Chapter 3. Conventional continuum modeling of homoepitaxial growth 45 where a derivation of the scaling of the noise-term in Eq. 3.29 can be found in Refs. [BS95, Kru97]. Dynamical scaling invariance requires that Eq. 3.30 can be written as the original equation we started with, the E W equation. This is obtained by "counting the powers" in Eq. 3.30. In order to regain the E W equation, the prefactor for the V2h term must equal v, which can be achieved when z = 2. Substituting this value into the prefactor for the noise-term gives the general scaling relation for the E W equation: ( z = 2 E W scaling : { a = ^ (3.31) In order to see the relevance of the nonlinear term with respect to the linear term in the K P Z equation, we scale it: (V/*) 2 = b2{-a~1\Vh)2 (3.32) Inserting (3.32) into (3.29), and using the results for the scaling parameters (a, /?, z) from the linear case (Eq. (3.31)), we arrive at > A = baX = b{2'd)/2X (3.33) This means that the K P Z nonlinearity is relevant as long as the prefactor amplifies during scaling, i.e. as long as in Eq. 3.33 we have b^-W2 > 1. This creates the condition that d < 2, and we refer to this as the "critical dimension" for the K P Z equation. • By applying power counting to the K P Z equation, we have learned that the nonlinear term is irrelevant in the hydro-dynamic limit for the practical case d = 2. There is, however, a situation where power counting does not give a fully reliable answer, namely in the strong coupling regime. This occurs when the prefactor, A, exceeds some critical value. It has been shown that for the K P Z equation, such a situation occurs whenever Xho/v > 1, where ho is proportional to the amplitude of the largest perturbation in the surface height [NB96, NK96]. It is not clear at present whether strong coupling applies to other nonlinear interface models [Kru97]. The various growth terms considered in this Chapter scale as follows: E W Mullins K P Z Villain Das Sarma V2h -• ba~2V2h V4h 6 a ~ 4 V 4 / i (\7h)2 b2a-2(Vh)2 V 2 ( V r i ) 2 b2a~4V2(Vh)2 V ( V / i ) 3 - » 6 3 a " 4 V ( V / i ) 3 (3.34) Chapter 3. Conventional continuum modeling of homoepitaxial growth 46 From these relationships, it can be seen that the Mullins term is irrelevant when compared to the E W term, because V4h : V2h ~ b~2 (3.35) This means that a completely linear growth equation scales according to the parameters determined in (3.31), i.e. the E W scaling parameters. We also see that the K P Z term dominates the Villain term (~ 62), and the Das Sarma term dominates both the Mullins term (~ b2a) and the Villain term (~ ba), as indicated in Section 3.3.4. In the absence of the E W term then, the Das Sarma equation (3.16) scales according to the term V(V/z) 3 in the hydro-dynamic limit. It might appear as if the order of the linear term must be equal or smaller than the order of the nonlinear term in a growth equation in order to agree with scaling arguments. This is not always the case, as was shown by Golubovic et al. [GB91], who investigated the following equation: dth= -KV4h + ±(Vh)2 + r) (3.36) It was found that the nonlinear term was relevant ford < 8. However, as the nonlinear K P Z term is nonconservative, it must always be accompanied by'an evaporation/condensation term proportional to V 2 / i [Kru97]. Chapter 4. Epitaxial growth experiments on GaAs (001) 47 Chapter 4 Epitaxial growth experiments on GaAs (001) Theoretical analysis of different growth mechanisms yields different contin-uum equations [Vil91, BS95, HHZ95, Kru97], so determining the equation that best fits the experimental data can provide valuable insight into the dominant physical processes involved in the growth. As an example, the presence of an ES barrier [EH66, SS66], which inhibits downward diffu-sion of mobile atoms at step edges, leads to a lowest order growth equa-tion dth oc —V 2 / i . This equation is unstable, and a small initial pertur-bation in h will grow with time. This seems like a suitable mathematical description of mounding, and it has been interpreted that way by several groups [OJS+94, JOH+94, HOW+94, AHD+00]. In these references, coun-tering mechanisms were introduced into the growth models in order to limit the growth and extent of the mounds. We discussed some of these mecha-nisms in Chapter 2. In this Chapter, we will investigate mounding and large-scale morphol-ogy of GaAs by looking into the the low-amplitude, long-wavelength limit of homoepitaxial growth. We put two continuum models to the ultimate test by comparing their predictions to the real world of epitaxial growth. This is done by growing thin films on top of GaAs substrates with various starting surfaces, at different temperatures and for different thickness films. We then compare the surfaces and power spectra obtained from A F M images of the experimental surfaces and the corresponding images from the growth sim-ulations based on the K P Z equation and the M B E equation1. The epitaxy experiments that we report in this Chapter were performed in great part by 1Versions of this Chapter have been published. A. Ballestad, B. J. Ruck, M . Adamcyk, T. Pinnington and T. Tiedje (2001) Evidence from the surface morphology for nonlinear growth of epitaxial GaAs films, Phys. Rev. Lett. 86: 2377-2380 [BRA+01]; A. Ballestad, B. J. Ruck, J. H . Schmid, M . Adamcyk, E. Nodwell, C. Nicoll and T. Tiedje (2002) Surface morphology during GaAs molecular beam epitaxy growth: Comparison of exper-imental data with simulations based on continuum growth equations, Phys. Rev. B 65: 205302 [BRS+02]. Chapter 4. Epitaxial growth experiments on GaAs (001) 48 Martin Adamcyk, with some aid from Ben Ruck, Tom Tiedje and myself. Further experimental details beyond those outlined here can be found in Martin Adamcyk's thesis [Ada02]. 4.1 S i m u l a t i o n s The critical exponents given in Table 3.1 for the different growth equations are difficult to determine experimentally,. due to a limited dynamic range accessible during growth. Furthermore, for real systems the morphology of the starting surface may also affect the scaling properties of the system at finite growth times. During this transient regime, measurements of the scaling exponents may not reveal the universality class. Instead, it is desir-able to compare the measured surfaces to simulations of the different growth equations. We have discretized the continuum growth Eq. 3.11 (KPZ equation) including the modification described by Eq. 3.15, as well as Eq. 3.12 ( M B E equation). Numerical instabilities restrict the range of parameters which can be used in typical finite difference implementations of the nonlinear terms [NB96]. Therefore, we have used a novel, more stable implementation based on the normal growth approximation to the K P Z term (Eq. 3.13). The algorithm translates all points on the surface outwards along the normal by a constant amount, thereby providing an excellent way to approximate the (V/i) 2 term. Details of the implementation of this term and the other terms in the growth equation are provided in Appendix C. Nonconservative noise is included in the simulations by adding an amount rncT rj\/l2AtC/(t) to each point on the surface at each time step [BS95], where a 2 = 2D/(Ax)2, Ax is the spacing between the lattice points in the simu-lation, At is the simulation time step, U(t) is a random number uniformly distributed between —0.5 and 0.5, and Tn is a dimensionless fitting parame-ter. For flux noise, one expects Fn = 1, and D = Fa±a2 in our units [BS95], where F is the flux in nm/s, ay is the in-plane atomic spacing of'0.4nm and a± is the GaAs monolayer height of 0.28 nm. Conservative noise can also arise during the growth process due to ther-mal fluctuations in the binding site of surface atoms. We have implemented an algorithm to approximate conservative noise in the following way. At each time step, we add/remove an amount Tcac^/V2U(t) from each point on the surface and transfer exactly this amount to one of the nearest neighbor sites. This is repeated for each of the four nearest neighbors, using peri-odic boundary conditions. The coefficient o~c = y/At/(Ax)4 ensures that Chapter 4. Epitaxial growth experiments on GaAs (001) 49 the amplitude of the PSD generated by the simulations does not change if we vary the lattice spacing Aa; or time step At. Fc is a fitting parameter with dimensions n m 3 / s 1 / 2 that adjusts the overall noise level to optimize the match between the simulations and experimental data. Generally, Fc is expected to increase as either the number of mobile atoms or the indi-vidual adatom diffusion constant increases, by increasing the growth rate F or the growth temperature T, respectively. Although other, more detailed implementations of conservative noise are certainly possible, our algorithm is adequate for testing the continuum growth equations. We have simulated the K P Z and the M B E equations numerically in order to generate the PSD at a range of different times, using a flat surface as the initial condition. The fourth order M B E equation was simulated using both conservative and nonconservative noise. In all cases the critical exponents extracted from the PSDs agree within uncertainty with the values given in Table 3.1, confirming the accuracy of our numerical schemes. 4.2 E x p e r i m e n t a l m e t h o d All samples in this study were prepared on [001] oriented GaAs substrates in a VG-V80H M B E chamber equipped with solid source effusion cells for both group III and group V elements. The growth rate, determined by the Ga flux, was kept constant at around 1 fim/hr during each run. A valved cracker was used as the As source, under conditions which give almost entirely A s 2 as the group V flux. The group V to group III beam equivalent pressure (BEP) flux ratios were estimated from measurements made with an ion gauge placed in front of the samples. The substrate temperature was monitored by optical bandgap thermometry [Joh95, JLT93, JT97] throughout the growth, with an absolute accuracy of about ± 5 ° C . Prior to loading into the M B E chamber, the substrates were exposed to ultraviolet light and ozone to remove carbon contaminants. The resulting surface oxide was removed in-situ, either by thermal evaporation at 600° C under an A s 2 overpressure, or by an atomic hydrogen etch. The hydrogen used in the etch was cracked with a W filament placed in front of the sample. During the etch the substrate was nominally set to 150°C, but radiation from the W filament caused the temperature to rise somewhat above this value. Details of the parameters associated with each sample can be found in Table 4.1. During growth, the surface roughness was monitored by elastic light scattering. Light from a chopped Hg arc lamp was incident through a quartz viewport onto the growing sample. The surface roughness diffusely Chapter 4. Epitaxial growth experiments on GaAs (001) 50 Table 4.1: Summary of growth parameters for the samples described in the text. Tsut>. is the substrate temperature and V:III is the Ga to As flux ratio (BEP) during growth. The oxide removal method is either thermal desorption (TD) or hydrogen etch (HE). Sample Ts ub. Growth Time V:III Surface Prep. (°C) (mins.) 730 TO - 0 - T D 744 HI 595 75 6.5 H E 755 T4 550 75 2.9 T D 769 HO - 0 - H E 776 H2 552 10 5.5 H E 780 H3 553 37.5 6.5 H E 836 T6 600 69 8.3 T D 838 T2 550 37.5 7.8 T D 839 T l 550 10 8.3 T D 841 T3 550 150 8.4 T D 852 T5 550 3 8.5 T D 899 H5 550 30 1.0 H E 912 H4 550 30 3.0 H E scatters some of the incident light out of the chamber through two addi-tional viewports, where it is detected by photomultiplier tubes fitted with wavelength selective filters. The chamber geometry, and the wavelength of light monitored, allows us to measure the roughness at spatial frequen-cies of 16 and 41/xm _ 1 , corresponding to length scales of 393 and 153 nm, respectively. It can be shown that the measured light scattering signal is proportional to the PSD of the surface at the spatial frequency moni-tored [ABP+00, CJZ77, CJZ79]. The scattering vectors for the two signals are 22.5° apart in the plane of the wafer, so the roughness is probed along slightly different crystallographic directions on the film surface. A F M images were obtained from the samples soon after removal from U H V , using a Digital Instruments Multimode Scanning Probe Microscope. Scans ranging from l x l / /m 2 to 100 x 100 pm2 were taken in tapping mode using Si tips with approximate tip radius 30 nm [Pin99]. Therefore, we can obtain reliable data up to spatial frequencies of at least 100 / i m - 1 before tip convolution effects become significant. Chapter 4. Epitaxial growth experiments on GaAs (001) 51 Figure 4.1: A F M images of a GaAs wafer (a) after thermal removal of the surface oxide (sample TO; z-range 50 nm) and (b) after 75 minutes growth at 550° (sample T4; z-range 8nm, brighter is higher). The arrow points along the [110] direction, and the image size was 10 x 10 jim2. 4.3 R e s u l t s a n d ana l ys i s 4.3.1 Growth on thermally desorbed substrates Time dependence Figure 4.1 shows a GaAs surface immediately prior to growth, after thermal cleaning of the surface oxide. The thermal desorption process produces a surface covered with pits, with the largest pits separated by a characteristic distance of around 1 fim, and having a maximum depth of around 40 nm. Figure 4.2 shows a set of A F M scans from samples grown for different times at 550°C on thermally cleaned substrates, using a similar A s 2 flux for each growth. Each of the grown surfaces is covered with mounds elongated along the [110] crystal axis, and separated by sharp V-grooves. The mounds have a characteristic separation similar to that of the pits on the starting surface, strongly suggesting that they are a remnant of the initial condition. The root-mean-square (RMS) roughness of the surfaces shown in Figs. 4.1 and 4.2 decreases steadily with time during growth, progressing through 4.9 nm, 1.0 nm, 0.7 nm, and 0.5 nm after 0, 10, 37.5, and 150 minutes of deposition, respectively. This smoothing is inconsistent with the idea of unstable growth [OJS +94, JOH + 94] , in which the mounds would be ex-Chapter 4. Epitaxial growth experiments on GaAs (001) 52 Figure 4.2: A F M images from thermally desorbed samples grown under nominally identical conditions, but for different times, (a) 10 minutes growth (sample T l , z-range 11 nm). (b) 37.5 minutes growth (sample T2, z-range 7nm). (c) 150 minutes growth (sample T3, z-range 5nm). The arrows point along the [110] direction, and the image size was 10 x 10/xm2. Chapter 4. Epitaxial growth experiments on GaAs (001) 53 Figure 4.3: 10 x 10 //m 2 simulations generated using the conservative form of the K P Z equation with growth times of (a) 10 minutes, and (b) 37.5 minutes. The z-range is 7.5 nm for both images, and the arrows point along the [110] direction. The simulations compare well with the real surfaces shown in Fig. 4.2 (a) and (b), respectively. pected to increase in amplitude as time progresses. A careful comparison of Figs. 4.2 (a)-(c) shows also that the number of mounds decreases as the surface smooths, while the lateral length scale of the remaining mounds in-creases. This is consistent with the behavior expected from the second order nonlinear term in the K P Z equation (Eqs. 3.11 and 3.15); the mounds grow outward in time such that the larger mounds absorb the smaller ones. If a single continuum equation is to be used to describe the surface evolution, then the equation must be capable of reproducing the morphology after a range of different growth times. The persistence of V-grooves on the surface suggests that a nonlinear equation, such as the K P Z or M B E equations, is the correct choice. Figure 4.3 shows two surfaces simulated using the K P Z equation with parameters vx = 10nm 2/s, vy = l n m 2 / s , and A x = Xy = 12nm/s, where the subscripts x and y correspond to the [110] and [110] directions, respectively. Nonconservative noise was included in each of the simulations with strength Tn = 10, considerably larger than the value T n = 1 expected from the random arrival of atoms from the flux (see Section 4.1). Figure 4.4 shows two surfaces simulated using the M B E equation with parameters 10 5 nm 4 /s , Ky — 10 5 nm 4 /s , Ax = 10 6 nm 3 /s , and Ay = 10 5 nm 3 /s . Conservative noise was added to these simulations with strength Fc = 100nm 3 /s 1 / 2 . In both cases the simulations Chapter 4. Epitaxial growth experiments on GaAs (001) 54 (a) (b) Figure 4.4: 10 x 10 pm 2 simulations generated using the fourth order molecular beam epitaxy equation with growth times of (a) 10 minutes (z-range 14 nm), and (b) 37.5 minutes (z-range l l nm) . The arrows point along the [110] direction. The fourth order equation also produces a mounded surface, although the depth of the cusps between the mounds is greater than on the real surfaces. A comparison between these images and the real surfaces shown in Fig. 4.2 (a) and (b), are not as favourable as the K P Z simulations shown in Fig. 4.3 were. are performed on a 256 x 256 grid, using an A F M scan of the desorbed wafer as the initial condition. For both sets of simulations the times match the shortest two growth times in Fig. 4.2 (10 and 37.5 minutes), and the coefficients have been ad-justed to provide the best possible agreement with the experimental data, as determined by visual inspection. Note that the coefficients used for the K P Z simulations are slightly different to those used in Ref. [BRA +01] to model a 75 minute growth. This is due to slight variations in the surface morphology caused by the higher A s 2 overpressure used during the growth of the present samples. Without the nonlinear terms both growth equa-tions generate inversion symmetric surfaces [Vil91, BS95, HHZ95, Kru97], whereas with our inclusion of the nonlinear terms both equations cause the etch pits on the starting surface to develop into mounds separated by V-grooves, in agreement with the experimental surface morphology. However, while the surfaces generated by the K P Z equation have similar roughness to the real surfaces at all length scales, the M B E equation generates less smoothing of the largest features than is seen in the experiment. Modifying the parameters of the M B E equation to enhance smoothing of the deepest Chapter 4. Epitaxial growth experiments on GaAs (001) 55 pits causes the smaller mounds to be less prominent than those on the real surface. Figure 4.5 (a) shows the PSDs of a thermally desorbed substrate (sample TO), a sample grown for 10 minutes (sample T l ) , and a sample grown for 150 minutes (sample T3), measured along the [110] direction. Figure 4.5 (b) shows the PSD from the same surfaces measured along the [110] direction. PSDs calculated from A F M scans ranging in size from l x l pirn2 to 100 x 100 /am2 were combined to generate these Figures. The P S D shrinks rapidly as a function of time during growth, until it reaches a saturated level for spatial frequencies greater than a crossover frequency qc (indicated by the vertical dashed lines). The crossover frequency decreases monotonically as a function of growth time, going from around 36 u,m~l to around 9 / i m - 1 ([110] direction), or from around 60 jim~l to around 17/_im_ 1 ([110] direction) after 10 and 150 minutes growth respectively. In the saturated region (q > qc) the PSD is well described by a power law, with slope close to —2, and a magnitude which is the same for all samples grown under similar conditions. The peak in the PSD of the thermally desorbed surface, located at around 4 /jm" 1 , moves gradually towards lower spatial frequencies during growth as the average size and spacing of the mounds increases. The solid lines on Fig. 4.5 show the PSDs obtained from simulations using the conservative form of the K P Z equation (Eq. 3.11 with the correc-tion in Eq. 3.15). Figure 4.6 shows the PSD of the same samples shown in Fig. 4.5, but in this case the solid lines are PSDs obtained from simulations using the M B E equation. The simulation times match the experimental data, and the parameters are the same as those used to generate the images shown in Figs. 4.3 and 4.4. As above, the PSDs are generated from a combination of simulated surfaces of different sizes. In the saturated region q > qc the magnitude of the PSD is determined by both the coefficients in the growth equation and the strength of the noise included in the simulations. The values of the parameters Fn (KPZ equation) and Tc ( M B E equation) quoted above were determined by matching the simulated and measured PSDs in the saturated region at high spatial frequencies. The K P Z simulations in Fig. 4.5 are in excellent agreement with the experimental data over the entire q range. For q > qc the M B E equation with conservative noise predicts the correct slope on the log-log plot (Fig. 4.6). However, unlike the K P Z equation, it does not predict the time dependence of the cutoff frequency, and for q < qc it tends to predict a slope steeper than the experimental data. Very little smoothing takes place at the lowest spatial frequencies in the M B E equation simulations, whereas the experimental data Chapter 4. Epitaxial growth experiments on GaAs (001) 56 10' Q CO CL , 10' 10 A (a). • — ' ° * =s » v \ * \^ X. \ ^ \ \ ^ -A 0 minutes ! ^^ _^>Sv * 10 minutes | o 150 minutes \ ^ ^ ^ . ^ ^ B ^ [110] S l o p e - 2 | \ ^ | j | ^ 10 10 q ( urn-1 ) 10' Q CO CL 10 A 0 minutes * 10 minutes o 150 minutes 10 10' q ( urn-1 ) 10 Figure 4.5: PSD after different growth times on thermally cleaned substrates, measured along (a) the [110] direction, and (b) the [110] direction. The symbols represent the experimental data. The solid lines, representing 10 and 150 minute simulations generated with the K P Z equation, are in excellent agreement with the experimental data. The vertical dashed lines indicate the cutoff frequencies qc at 10 and 150 minutes (right to left). Chapter 4. Epitaxial growth experiments on GaAs (001) 57 10' i10< Q CO CL , 10' 10 (a). * H*. \ ^ °° \ >A X . A 0 minutes * 10 minutes o 150 minutes [110] Slope - 2 10 10 q (u.rrf1 ) 10 10" Q 00 CL , 10' 10' 10 A 0 minutes * 10 minutes o 150 minutes 1 0 _ q (urn - 1 ) Figure 4.6: PSD of the same samples shown in Fig. 4.5, measured along (a) the [110] direction, and (b) the [110] direction. The solid lines are from 10 and 150 minute simula-tions generated with the M B E equation. The simulated PSD has a stronger q dependence than the PSD of the grown films. The vertical dashed lines indicate the cutoff frequencies qc at 10 and 150. minutes (right to left). Chapter 4. Epitaxial growth experiments on GaAs (001) 58 0 10 20 30 40 50 60 70 Time (mins.) Figure 4.7: Diffuse light scattering signal at 16 and 41 /j,m _ 1 measured during growth of sample T4. For clarity, the 41/xm _ 1 signal has been offset. The insets show semi-logarithmic plots of portions of the same data during the growth, after subtraction of the constant background level. The solid lines on both plots show the predicted light scattering from simulations using the K P Z equation. clearly show a decrease corresponding to filling in of the thermal desorb pits. The time dependence of the surface roughness can also be monitored by elastic light scattering. Figure 4.7 shows the measured light scattering signal at spatial frequencies of 16 and 41 / j m _ 1 during the thermal cleaning and subsequent growth of sample T4. The thermal desorption of the surface oxide at 600° C, which takes place at around 10 minutes in the Figure, is accompanied by a rapid increase in the scattered light intensity at both spatial frequencies, caused by the appearance of the desorption pits seen in Fig. 4.1. The sample is maintained at 600°C for several minutes to complete the oxide removal, during which time the surface smooths considerably at 41 / jm _ 1 but stays approximately constant at 16 fim" 1 , corresponding to the annealing of some short scale features. The sample is brought down to growth temperature (550°C) after 23 minutes, and growth begins at a time of 31 minutes. The surface immediately begins to smooth at both length Chapter 4. Epitaxial growth experiments on GaAs (001) 59 scales, until a background level is reached after about 10 minutes of growth. The light scattering data indicate that, during growth, the surface smooths at a rate which is approximately independent of the probed length scale. In fact, the data are well described by an exponential decay with time constants Tie and T41 both equal to around 3 minutes. A n exponen-tial decay would be expected from purely linear equations such as the E W equation or the fourth order linear term in the M B E equation. However, the E W equation predicts the characteristic smoothing times to be in the ratio T41/T16 = (16/41)2, and the linear M B E equation predicts the ratio (16/41)4, both inconsistent with the observed ratio. The predicted smoothing rates of nonlinear equations, such as the K P Z equation, are harder to determine analytically, but for a given starting sur-face they can be determined from simulations. The inset in Fig. 4.7 shows the measured light scattering signal during growth and the scattering sig-nal calculated from simulations using the K P Z equation. The agreement between the calculated and measured data is quite reasonable, as the only parameter in the calculation is the initial amplitude. Al l other parame-ters were taken to be identical to those employed in Ref. [BRA +01] to fit the A F M images from this sample, namely vx = 10nm 2/s, uy = l n m 2 / s , A x = \ = 5nm/s, and Yn = 10. Most notably, the calculated time con-stant is only weakly q dependent (much less so than for simulations based on the linear growth equations), and is within a factor of two of both of the measured values. During the growth of this sample, the plane of incidence of the Hg lamp was at 45° to the [110] and [110] axes. The 16 / j m - 1 signal is measured in the plane of incidence, while the 41 / j m - 1 signal is measured at an angle 22.5° out of the plane of incidence, rotated towards the [110] direction. Therefore, the 16 / im" 1 signal measures the roughness along a direction closer to the [110] axis, where the value of v is largest, than the 41/am - 1 signal, which measures along a direction with a smaller value for v. This partially off-sets the expected faster smoothing rate at 41 /xm _ 1 . More importantly, the nonlinear K P Z term has a large effect on the smoothing rate early in the growth when the surface is at its roughest. The nonlinear term tends to fa-vor a weaker q dependence in the smoothing rate than the linear term in the growth equation. Further simulations have shown that this weak sensitivity to length scale cannot be reproduced by the M B E equation. Finally, in Figure 4.8 (a) we show a series of scan lines from A F M im-ages of samples grown at 550° C. From bottom to top the cross sections come from the thermally desorbed starting surface, and samples grown for 3, 10, 37.5, and 150 minutes (samples TO, T5, T l , T2, and T3, respectively). The Chapter 4. Epitaxial growth experiments on GaAs (001) 60 0 2 4 . 6 8 Position (u.m j Figure 4.8: Scan lines along the [110] direction from (a) measured A F M images, (b) the corresponding K P Z simulations, (c) M B E equation simulations, and (d) EW equation simulations. From bottom to top the curves correspond to growth times of 0, 3, 10, 37.5, and 150 minutes. The lines in (a) are from five separate samples, whereas the scan lines in (b), (c), and (d) are taken at the same position in the evolving simulations. Scans are offset for clarity. Chapter 4. Epitaxial growth experiments on GaAs (001)' 61 scan lines are taken along the [110] direction and are offset by 15 nm for clarity. Figure 4.8 (b) shows scan lines from the K P Z simulations, where the different lines are extracted at the same times as those in Fig. 4.8 (a). Figure 4.8 (c) shows a similar set of scan lines from the M B E equation sim-ulations. The initial condition for the simulations is the thermally desorbed surface (sample TO). These lines are all taken from the same position on the surface as the simulation progressed. It is clear that while the K P Z sim-ulations are in excellent agreement with the real surfaces at all times and length scales, the M B E equation does not correctly predict the rate at which the deepest pits fill in. This, again, indicates that the rate of change of the surface morphology during growth of this material system is only weakly q dependent. In Fig. 4.8 (d) we also include a series of scan lines extracted from a simulation using the E W equation with vx = 20nm 2/s, uy = 5nm 2 /s , and r n = 10 (i.e., the K P Z equation with A^ = Xy — Onm/s). As before, the parameters have been selected to optimise the similarity between the simu-lated and real surfaces. Although we obtain reasonable agreement with the real data after long growth times, the large features associated with the des-orption pits clearly smooth too slowly in the E W simulation. Furthermore, the cusp-like features reproduced in the K P Z simulations become rounded in the purely linear simulation, demonstrating the importance of the nonlinear K P Z term. Temperature dependence Figure 4.9 (a) shows an A F M scan from sample T6, grown for 69 minutes at 600° C on a thermally desorbed substrate. The surface still shows large scale mound-like features related to the initial roughness. However, the V -grooves between the mounds are not apparent, and the surface is much more inversion symmetric than the samples grown at 550° C. To simulate the 600°C growth, we have scaled the linear coefficients in the K P Z equation by a factor of three relative to their values at 550°C. This is the ratio that would be expected if the linear term represents a thermally activated process with activation energy of around 0.7 eV. Using valuesfor the nonlinear coefficients of \ x = \ y = 5nm/s, and Fn = 10 results in the simulation shown in Fig. 4.9 (b). As above, the K P Z simulation generates an excellent likeness of the experimental data. The increase of the linear relative to the nonlinear coefficients enhances the anisotropy of the surface structure, and at the same time reduces the inversion asymmetry. Overall, the grown surface is smoother at 600° C than at 550° C. Chapter 4. Epitaxial growth experiments on GaAs (001) 62 (a) (b) Figure 4.9: High temperature growth: (a) A F M image from sample T6, grown on a thermally desorbed substrate (z-range 4 nm); (b) simulation using the KPZ equation with vx = 30nm 2/s, vv = 3nm 2/s, and A* = Aj, = 5nm/s (z-range 5nm). The morphology of the 600°C grown surface is reproduced. The arrows point along the [110] direction, and the image size was 10 x lO/um2. 4.3.2 Growth on hydrogen etched substrates Kinetic roughening Figure 4.10 (a) shows an A F M image of a substrate from which the surface oxide has been removed by hydrogen etching. Unlike the thermal desorp-tion process, the hydrogen etch leaves a relatively smooth surface, with an RMS roughness of less than 0.2 nm in this sample. Figure 4.10 (b) shows an A F M image from sample HI, grown for 75 minutes at 595°C on a hydrogen etched substrate. The surface remains relatively smooth, with RMS rough-ness equal to 0.2 nm. The large amplitude mounds seen on the thermally desorbed samples are absent in this case, but the roughness does appear to be correlated over larger length scales than on the initial substrate. A 75 minute growth simulation generated using the K P Z equation with parame-ters vx = 30nm 2/s, vy = 3nm 2 /s , Xx = Xy = 5nm/s, and Tn = 10, is shown in Fig. 4.10 (c), where the hydrogen etched surface was used as the initial condition. These parameters, which are the same as those used above to successfully model the growth at 600° C on thermally desorbed substrates, provide an excellent likeness of the grown surface. Because the starting surface is almost flat, we can investigate the kinetic Chapter 4. Epitaxial growth experiments on GaAs (001) 63 Figure 4.10: A F M images of GaAs surfaces which have had the surface oxide removed by hydrogen etching: (a) Starting surface (sample HO, z-range 1.8nm); (b) 75 minutes growth at 595°C (sample HI, z-range 1.2nm); and (c) K P Z simulation using sample in (a) as initial condition (z-range 1.1 nm). The parameters are the same as those used to simulate the growth on thermally desorbed substrates at 600° C. The arrows point along the [110] direction, and the image size was 5 x 5 / W 2 . Chapter 4. Epitaxial growth experiments on GaAs (001) 64 10' 10" V°2 Q CO CL 10 (a) q ( urn 1 ) „ o o * \ '(b) _ ° o °° ;\°° ** 1 ^sSx$»> — i s^ro i * 10 minutes i o 37.5 minutes ] [110] Slope - 2 N. N 10 q ( urn 1 ) 10' Figure 4.11: PSD of samples grown on hydrogen etched substrates, measured along (a) the [110] direction, and (b) the [110] direction. The growth times are 10 minutes and 37.5 min-utes for samples H2 and H3, respectively. The solid lines are from simulations generated with the K P Z equation, using the same parameters used to model the growth on thermally desorbed substrates. The vertical dashed lines in (b) indicate the cutoff frequencies qc at 10 and 37.5 minutes (right to left). Chapter 4. Epitaxial growth experiments on GaAs (001) 65 roughening- of the' films grown on hydrogen etched substrates. The PSD of two such films, grown at 550°C for 10 and 37.5 minutes, are shown in Figures 4.11 (a) and (b) (samples H2 and H3 respectively). The PSD of the starting surface, not shown in the figure, lies slightly below the PSDs of the grown films, and has a sharper roll-off at high spatial frequencies. In the [110] direction, the PSD of the grown films follows a power law with exponent approximately —2 over the entire q range. The amplitude of the PSD in this direction is independent of time, indicating that a saturated roughness has been reached. In the [110] direction the PSD also follows a power law with exponent —2, but only for spatial frequencies greater than a cutoff frequency qc of about 20 after 10 minutes growth, and about 10/xm _ 1 after 37.5 minutes growth (indicated by the vertical dashed lines on the Figure). Below the cutoff frequency the PSD is approximately independent of q, and has an amplitude which increases with time. This behavior is exactly what is expected from kinetic roughening; correlations develop at longer and longer length scales as time progresses, leading to a saturation in the PSD which shifts to smaller and smaller spatial frequencies. The PSD saturates more quickly along the [110] axis due to the more rapid smoothing rate in this direction. The solid lines on Figs. 4.11 (a) and (b) show simulated PSDs gener-ated with the K P Z equation using the same parameters used to generate the simulated PSDs in Figs. 4.5 (a) and (b), but with a hydrogen etched starting surface. The simulations are in excellent agreement with the mea-sured data. Most notably, the slope and amplitude of the saturated PSDs are reproduced, and the position of the cutoff frequency qc is in good agree-ment with the data. Therefore, we conclude that the modified K P Z equation describes the growth on both smooth and rough initial surfaces. The exponent of the power law describing both the measured and the simulated PSDs is close to —2 in both directions. This is in agreement with the exponents measured on the thermally desorbed surfaces, although it is somewhat surprising, as it differs from the K P Z prediction of —2.8. We return to this point in Section 4.4, below. Dependence on As overpressure The anisotropy of the surface structure is sensitive to the A s 2 overpressure during growth. This is demonstrated in Figure 4.12 which shows three high resolution A F M scans from samples grown at 550° C on hydrogen etched sub-strates under varying A s 2 flux, (samples H3 (a), H4 (b), and H5 (c)). Well defined atomic steps are visible in all three images. At a group V to group III Chapter 4. Epitaxial growth experiments on GaAs (001) 6 6 Figure 4.12: A F M images illustrating the effect of reducing the As2 overpressure, (a) V:III ratio = 6.5 (sample H3), (b) V:III ratio = 3.0 (sample H4), and (c) V:III ratio = 1.0 (sample H5). The gray scale in all images are 2.4 nm, and the surface features are elongated along the [110] direction (along the arrow), and the image sizes were 2 x 2/j.m2. The 2D PSD of each image is included as insets, with spatial frequencies ranging from —300 to 300 / i m " 1 in each direction. Chapter 4. Epitaxial growth experiments on GaAs (001) 67 flux ratio of 6.5 [Fig. 4.12 (a)] the surface is covered with islands which show a moderate amount of elongation along the [110] direction, as demonstrated by the relatively isotropic two dimensional PSD shown in the inset. As the V:III ratio is decreased to 3.0 [Fig. 4.12 (b)] the anisotropy increases, and the PSD becomes much more elongated. In terms of the continuum growth equations this corresponds to an enhancement in the anisotropy of the coeffi-cients. Physically, the elongation results from an anisotropic mobility of the adatoms, and possibly also an anisotropy in the incorporation dynamics on the surface [ItoOl]. Once the V:III ratio approaches unity [Fig. 4.12 (c)] the surface structure changes considerably, possibly due to an imminent change in the reconstruction from 2 x 4 to 2 x 2. The surface is now covered by terraces that are almost continuous along the [110] direction, each with a characteristic width of around 80 nm. This length scale is reflected in the PSD, which now shows two peaks on either side of the origin perpendicular to the direction of elongation. The GaAs surfaces described in Refs. [JOH+94, OJS+94, HOW+94], on which mounds were attributed to unstable growth, also show a large degree of anisotropy. This indicates that they were grown under an effectively lower As overpressure (AS4 was used as the source) than most of the samples displayed here. It is therefore interesting to consider whether an appreciable Ehrlich-Schwobel barrier may emerge in the [110] direction in conditions of low As flux. This could explain the length scale that we observe emerging in Fig. 4.12 (c), and is consistent with the flat tops that we observed in Pt growth in Fig. 2.11. The extensions of surface features in the [110] direction indicates that the ES barrier 'in this direction remains negative, or at least relatively unaffected by the As overpressure. 4.4 D i s c u s s i o n Our simulations show that both the K P Z and the M B E equations can quali-tatively reproduce the mounds on the surface, and clearly prove that mounds can emerge as a transient shape during smoothing, rather than as a mani-festation of unstable growth. Quantitatively, only the K P Z equation is able to reproduce the PSD of the surfaces. This is because the measured PSD smooths at a rate which is fairly insensitive to the length scale being probed. The stronger spatial frequency dependence inherent in the higher order M B E equation makes it difficult to reproduce the relatively q independent decrease in the measured PSD (see Figs. 4.5 and 4.6). While the K P Z simulations provide the better fits to the experimental surface evolution, it should be Chapter 4. Epitaxial growth experiments on GaAs (001) 68 mentioned that the values of A required in the simulations are considerably larger than the growth rate of around 0.3nm/s. Thus, the nonlinear K P Z term cannot be accounted for by simply assuming the surface grows outward along the surface normal. Our interpretation of the terms in the K P Z equation also provides a natural explanation for the excess noise in the K P Z simulations. The noise represents fluctuations in the density of adatoms at a given point on the surface. The fluctuations associated with adatom incorporation and dis-sociation from step edges can be much larger than the fluctuations in the arrival rate of adatoms from the flux. Furthermore, due to the existence of the 2x4 surface reconstruction, it may not be favorable for Ga adatoms to incorporate individually [ItoOl]. Instead, the basic unit which must be added to the surface may consist of several Ga atoms, enough to recreate a unit cell of the 2x4 reconstruction, for example. In this case the coeffi-cient of the noise correlation function will be modified from D = Fa±a2 to D = nFa±a2, where n is the number of Ga atoms forming the basic incor-poration unit. Therefore the value Tn = 10 used in the K P Z simulations is physically possible. We now comment on the power law observed in the PSD at high spatial frequencies. Despite the importance of the nonlinear terms for reproducing the surface morphology, the PSD tends to display the power law exponent ex-pected for the purely linear equations (i.e. —2). This should not be taken as indicating a disagreement with the theoretical predictions. Instead, it simply indicates that we are in a transient regime in which the relatively large lin-ear terms are dominating the nonlinear terms at short length scales [NB96]. The small surface roughness at short length scales means that the nonlinear term can be neglected in this regime. This conclusion is supported by nu-merical simulations of the K P Z equation. For example, using the isotropic parameters v = 5nm 2 /s , A = 5nm/s, and Fn = 10 the simulations exhibit a power law exponent of —2 in the PSD even after 1000 minutes of growth on a flat 10 x 10 /jm 2 substrate. However, using the parameters v = 0.01 nm 2 /s , A = 5nm/s, and Tn = 10, an exponent of —2.8 is observed for times greater than about 60 minutes, consistent with the expected value 2(1 — a)=2.76 as obtained from Eq. 3.25. 4.5 S u m m a r y Based on our comparison of experimental data with simulations of two pop-ular continuum growth equations, we conclude that the evolution of the sur-Chapter 4. Epitaxial growth experiments on GaAs (001) 69 face morphology of M B E grown GaAs is well described by the K P Z equation with a stable linear term. By contrast, the M B E equation fails to reproduce the measured surface morphology of the grown films. These results lead us to the following interpretation of the dominant smoothing mechanism during film growth. The surface is covered by mo-bile adatoms which diffuse randomly on the surface, forming an effective vapor phase. Incorporation of these adatoms takes place preferentially at sites on the surface with positive curvature, as described by a second order, linear growth equation. Nonlinear corrections to the incorporation rate are observed to be important. We attribute these corrections to the preferred in-corporation of adatoms to the sloped parts of the surface, where the steps are in abundance. The breach of inversion symmetry in the growth direction, as indicated by the transient mounds separated by cusps, is well approximated by the K P Z nonlinearity, although we speculate that this second-order term is a lowest order approximation to a conservative nonlinear term of higher order. We will return to this important point in Chapter 5. In M B E growth, the incorporation process must necessarily be of a con-servative nature, i.e. the growth rate should not depend on the surface shape. To account for this fact we have made a correction to the K P Z equa-tion. It is important to point out that, due to the rather small slopes on even the thermally desorbed substrates, this modification has only a minor effect on the morphology of the simulated surfaces. Therefore, the simulations although conservative are essentially indistinguishable from simulations of the unmodified K P Z equation. For a given growth condition, a single set of parameters successfully de-scribes growth on surfaces that are initially rough or smooth. Regardless of initial condition, after long growth times the surface tends towards a steady-state roughness level determined by the interplay between the smoothing rate and random noise in the system, as predicted by kinetic roughening theory. We have found that increasing the temperature or As overpressure leads to smoother surfaces, implying increased values for the coefficients of the smoothing terms in the K P Z equation. Increasing the As overpressure also reduces the anisotropy of the surface morphology. Our experimental studies lead us to draw new conclusions as to the stability of the GaAs (001) surface during epitaxial growth. While com-monly believed to be unstable, as first reported by Johnson [JOH+94] and Orme [OJS+94], we find that the long-time behaviour shows the opposite: GaAs (001) homoepitaxy is stable and the observed mounds are simply tran-sient remnants of the starting surface roughness. Since the mounds are no longer a result of an instability, there is no need for a strong, positive ES Chapter 4. Epitaxial growth experiments on GaAs (001) 70 barrier to explain the long-time surface morphology evolution. In fact, mod-eling the surface evolution with the K P Z model showed that the net adatom drift in this growth system must be downhill. It remains an unresolved issue how one can relate the K P Z coefficients v and A to real physical parameters, like F, T and the As overpressure, as well as physical processes, such as step edge incorporation and detachment, and nucleation. While it is clear that the E W parameter v is dependent on growth chamber temperature, growth rate and the ES barrier, it is unclear exactly how it is related to the underlying physics. We also have concerns about the origin of the growth terms themselves. For example, the linear E W term V 2 / i has been associated with either an evaporation/condensation mechanism, or an inverse ES effect. Since these processes may occur simultanously and independently, can their effects be expressed mathematically as multiplicative effects, or additive effects? And while we are fairly confident about the existence and origin of the linear E W term, what is the underlying mechanisms causing the nonlinear K P Z term, and is it at all valid in the context of M B E growth? If not valid, then why does it do such a good job at simulating the evolution of GaAs surface morphology? In Chapter 2, we discussed surface growth mainly from the point of view of the diffusing adatoms on flat terraces separated by steps. In Chapter 3, we discussed surface morphology directly in terms of the surface height itself and through terms that depend on it through spatial derivatives of various linear and nonlinear orders. We showed in this Chapter how to apply continuum equations in order to model surface growth. In the next Chapter, we will try and motivate a formulation that describes the interplay between the adatoms and the surface, and hence combines the approaches from these previous Chapters. We arrive at a surface growth model that includes physically based, continuum, morphological terms, and we put this model to test by applying it to the evolution of flat and patterned surfaces.. , Chapter 5. Modeling of growth on patterned surfaces 71 Chapter 5 Modeling patterned of growth on surfaces We showed in Chapter 4 that the K P Z equation adequately describes the evolution of GaAs (001) surfaces during epitaxial growth for nominally flat starting surfaces, where the roughness induced by common surface prepara-tion techniques is relatively small, and the local surface slopes do not exceed a few degrees. We also showed that in spite of its applicability to GaAs ho-moepitaxy, the K P Z equation faces some theoretical problems surrounding the interpretation and validity of the nonlinear term. In this Chapter we investigate patterned surfaces, with lithographically textured large-scale features like gratings and pits 1. Such surfaces are in-teresting from the point of view of device manufacturing, but conventional continuum equations derived in the low slope approximation fail to describe their evolution. We describe experiments where large scale gratings are pre-pared, epitaxy is performed and the evolution of the shape of the gratings is closely investigated with A F M . The epitaxy*regrowth experiments, as well as the grating preparation, were done by Jens Schmid and is described in detail in his Ph .D. thesis [Sch04]. We show that the K P Z equation breaks down for the large slope sur-faces, and develop a coupled equation continuum model that explains the complex surface shapes observed in these regrowth experiments. This model describes the dependence of the surface morphology on film thickness and growth temperature in terms of a few simple atomic scale processes including adatom diffusion, step edge attachment and detachment, and a net downhill migration of surface adatoms. The continuum model reduces to a nonlinear 1Versions of this Chapter have been published. A. Ballestad, T. Tiedje, J. H. Schmid, B. J. Ruck and M . Adamcyk (2004) Predicting GaAs Surface Shapes During M B E Re-growth on Patterned Substrates, J. Cryst: Growth 271: 13-21 [BTS+04b]; A. Ballestad, T. Tiedje and J. H. Schmid (2004) Comment on "Transient Evolution of Surface Rough-ness on Patterned GaAs (001) during Homo-Epitaxial Growth", Phys. Rev. Lett. 93: 159601 [BTS04a]; A. Ballestad, Bayo Lau, J. H. Schmid and T. Tiedje (2005) Nonlin-ear growth in GaAs molecular beam epitaxy, Mater. Res. Soc. Symp. Proc. 859E, JJ9.6 [BLST05]. Chapter 5. Modeling of growth on patterned surfaces 72 equation reminiscent of the K P Z equation in the long wavelength limit, with a smoothing rate that is dependent on the growth rate, the magnitude of the ES barrier, and temperature. 5.1 I n t r o d u c t i o n The problem of the time evolution of the shape of crystal surfaces has a long history dating back to Mullins and Herring who considered relaxation during annealing above the roughening temperature [PV98]. More recently, shape relaxation below the roughening temperature has been studied ex-tensively [DP97, EAC+00, SRF03]. Below the roughening temperature the problem is complicated by the need to keep track of the dynamics of atomic steps and the fact that the surface free energy of crystal facets is singular. Biasiol et al. [BGLK02] have extended the theory of shape relaxation below the roughening temperature to include the effects of atom deposition, and use this theory to explain the self limiting V-grooves observed in organo-metallic chemical vapour deposition growth on corrugated GaAs substrates. In this Chapter, we present a new continuum model which we use to inter-pret measurements of the shape of corrugated (001) GaAs surfaces under growth conditions which do not produce faceting. Facets are not present in our experiments due to atomic scale roughness associated with atom deposi-tion in the island growth mode, and the fact that the surface topography is sufficiently weak that the surface slope does not reach the low energy [111] facets. We show that this model reproduces the surface morphology that develops during M B E regrowth on ID surface gratings. 5.2 C o n v e n t i o n a l m o d e l i n g o f weak sur face t e x t u r e The nonlinear term in K P Z is associated with growth along the outward normal, as in chemical vapour deposition. In this case, which is not obvi-ously applicable in M B E , A should be equal to the growth rate F. Also, the value for A needed to simulate the experimental results is almost two orders of magnitude larger than F (see Ref. [BRS+02] and Chapter 4). This raises the question as to the physical origin of the non-linear term in the context of M B E growth. The K P Z nonlinearity is non-conservative, whereas M B E growth is conservative with a growth rate that is independent of the surface shape. In practice, for the low surface slopes where K P Z is applica-ble, the change in growth rate associated with the non-linear term is very Chapter 5. Modeling of growth on patterned surfaces 73 i 1 1 1 1 1 1 1 r 0 10 20 30 40 50 60 70 80 90 100 Time (minutes) Figure 5.1: Light scattering during growth corresponding to surface power spectral density at 41 ixnT1 showing the effect of atom deposition on the smoothing rate. The sample roughens during a temperature ramp to remove the surface oxide at about 5 minutes in the figure, which is followed by relatively fast smoothing during a high temperature (620° C) anneal for about 7 minutes, and then slower smoothing during annealing at growth temperature (550° C). small. Although the K P Z nonlinearity gives a good approximation to the surface shape, its physical origin is obscure in the case of M B E growth. One would prefer to have a model that can be derived from underlying physical phenomena. In addition, the K P Z description with constant coefficients is not con-sistent with experiments which show that the smoothing rate depends on the growth rate. For example, in Fig. 5.1 we show the scattered light inten-sity from a GaAs surface during an interruption in growth on a randomly textured substrate. The intensity of scattered light is proportional to the power spectral density of the surface topography at a spatial frequency q determined by the optical wavelength and geometrical factors [BRS+02]. In this case g=41 /xm _ 1 , corresponding to a lateral surface length-scale of about 150 nm. For low amplitude surface textures, in the K P Z model the surface should smooth exponentially with a characteristic rate given by uq2 where q is the spatial frequency of the surface roughness [PV98]. As shown in the Chapter 5. Modeling of growth on patterned surfaces 7 4 Q.5n m (a) (b) Figure 5.2: A F M images of (a) a sample quenched (fast cooled) after 69 minutes of growth at 600° C and (b) a sample annealed for 15 minutes at growth temperature 595° C after 40 minutes of growth. inset of Fig. 5.1, the smoothing rate responds immediately to changes in the growth flux; it is faster during deposition and slower during annealing, sug-gesting that v is flux dependent. This continued smoothing of the surface in the absence of an atom flux indicates that the physical mechanisms at play on the surface favor a net downhill migration of surface adatoms, even when the flux of atoms is turned off. This suggests that if there are two competing mechanisms including a positive ES barrier that the stabilizing mechanism that creates the downhill flow is not associated with energetic adatoms deposited from the vapor such as step edge knockout or downhill funneling. To summarize, the K P Z equation provides an accurate description of the surface morphology under certain restricted conditions (constant growth rate, low surface slope and long wavelength surface structures). In addition, the physical origin of the non-linearity in K P Z is unclear in the case of M B E growth. Therefore, we would like to develop a new continuum model based on a few simple physical processes that describes the surface morphology over a broader range of conditions than the K P Z equation, and that can be compared with experimental data. To illustrate the physical processes involved in the epitaxial growth pro-cess and to motivate the new model, we compare atomic force microscope Chapter 5. Modeling of growth on patterned surfaces 75 (AFM) images of a growth surface which is fast cooled after growth is termi-nated (quenched) with one that has been annealed (see Fig. 5.2 (a) and (b)). In this Figure and in Fig. 5.1, the As:Ga flux ratio was 10:1. The quenched sample (a) is covered with small islands, whereas the annealed sample (b) has broad terraces with few islands. The small islands must coalesce into the step edges during annealing. The kinetic barrier to the adatom coales-cence into the step edges, causes the growth process to be non-local. This means that the growth rate at one location will be affected by the morphol-ogy at another location in contrast to K P Z , for which the growth rate only depends on the local surface slope and curvature. This is a problem not only for K P Z but also for alternative growth equations such as the M B E equation [PV98]. Physically, a high density of steps at one location that absorb adatoms will affect the adatom density and hence the growth rate at another nearby location. We hope to have motivated the necessity for a model that includes the diffusion of adatoms. We proceed by discussing an approach that describes the concentration of islands of all sizes and that has been a popular description of nucleation during the early stages of growth. 5.3 S u b - m o n o l a y e r ra te equa t i ons Sub-monolayer rate equations (SMREs) describe the nucleation and growth of 2D islands on an initially flat surface though mean-field concentrations ns of islands of size s. The source of adatoms is as usual through a flux term F. The total number of islands with s atoms in them is Ns = ns*A, where A = 1? is the substrate area. If only, adatoms are mobile on the surface, then islands with s atoms can form either by the attachment of an adatom to an island of size (s — 1), or by the detachment of an adatom from an island of size (s + 1). The reverse processes reduce ns, and the total nucleation rate of islands of size 5 can conveniently be described by dn —j^- = Dm (n s _i - n s ) - K (ns - ns+1) (5.1) where D is the diffusion constant and K is the detachment rate of an atom from an island. The concentrations ns for various values of s are coupled: ns-i decreases when hit by an adatom that increases ns. The coupling leads to the following set of self-consistent SMREs [BE92, Kru97, KKvOO, vS16, Chapter 5. Modeling of growth on patterned surfaces 76 x 10" 5^ •5 2 £= 0) • D • D c « 1 0 3 1 1 1 0.05s — i 1 1 1 y ^ 0.10s X \ 0.15s 0.20s 0.25s v ~ \ \ \ \ i i i i i i i 0 10 20 30 40 50 60 70 # atoms in island, s Figure 5.3: Evolution of sub-monolayer rate equations (Eqs. 5.2) for Esub =1.3eV, mEiat =0.9eV, T = 600°C, F = 1 ML/s and as = -ys = 1 for all s. Growth with (solid lines) and without (dashed lines) detachment is shown, and the growth times are indicated. vS17, VveOO, GRG+01, Ven87, BC94, RZv94]: ^jT = F - Dni (2crim - T,s>2asns) at dns ~dt = Dni(as-ins-i - asns) -Kirfsns - 7 s + i n s + i ) , s = 2,3, (5.2a) (5.2b) where as is the "capture coefficient" [VSH84] for islands of size s, and are related to the size, stability and spatial distribution of islands. "fs is a corresponding " release coefficient" from a step edge. Al l concentrations ns have the initial value of 0. An example of the evolution of theses SMREs are shown in Fig. 5.3. We see that the average size of islands grows as more material is deposited on the surface, as well as their total count. The success and applicability of SMREs is well documented for a variety of growth systems, as evidenced by the large amount of literature on the subject (see earlier references in this Section). We have introduced these equations here for completeness, but we do not intend any further investi-gation into their current format, as we are mainly concerned with growth of surfaces that span several layers. The SMREs are limited to low coverages Chapter 5. Modeling of growth on patterned surfaces 77 (less than ~ 20%) on flat starting surfaces, as coalescence of islands is unac-counted for in its formalism. We now address the question of how to extend this formalism to account for coalescence within the layer and multilayer growth. 5.4 M u l t i l a y e r m e a n - f i e l d m o d e l s In order to extend the rate equation approach to a multilayer description, we choose to express the filling of each layer rather than the densities of each island size in each layer. We denote the filling of the jth layer above the substrate by 9j G [0,1]. The index j = 0 denotes the substrate, for which $o(t = 0) = 1, and a flat starting surface is described by {Vj, j > 0,6j(t — 0) = 0}. As with the SMREs in the previous Section, this formulation does not take into account the morphology of the growing surface. Also, the incoming atoms are considered to incorporate directly into the growing film, so there is no explicit variable for the adatom concentration. A multilayer description in 6 couples the growth of all active (that is non-empty and non-full) layers on the surface. 5.4.1 "Birth-death" models and diffusive growth In one extreme case, layers fill in at a rate proportional to their exposed area, and several layers will be active at any given time. This is called nondiffusive or Poisson growth. In this limit, all atoms that land in a layer incorporate to the ascending step edge. The filling of successive layers according to this model can be described by This is a birth-death (BD) model [WGJ76] "since the growth of an unfilled layer is rapid, while the growth of a nearly completed layer is slow" [CPP +89]. Nondiffusive growth refers to an experimental scenario where the adatom diffusion is low. This could be caused by either a strong binding energy between the adatom and the surface, or by the presence of a high growth flux that prevents the adatoms from diffusing before they are hit by another incoming atom from the vapour. In a simple model, such as the one de-scribed by Eq. 5.3 where no morphology is specified, it can also be used to describe a strong ES barrier at the step edge. (5.3) Chapter 5. Modeling of growth on patterned surfaces 78 Figure 5.4: Simple layer filling for the birth-death model, Eq. 5.3, on an initially flat starting surface. The growth rate used was F = l M L / s . 2 3 4 Time (s) The solutions for the first two layers are = 1 - e " -Ft -Ft ..-Ft l - e ~ r i - F t e -and the general solution for the jth layer can be written as: = l - e -Ft n=0 (Ft)n n\ • (5.4) (5.5) (5.6) where the sum, incidentally, is the Taylor series expansion of eFt, truncated after the jth term. The evolution of the first few layers according to this expression is shown in Fig. 5.4. The leftmost line in that plot is the filling of the first layer, the second line from the left is the second layer, and so forth. In the other extreme, each layer fills in completely before the next layer is started, a model often referred to as layer-by-layer growth. The filling of the active layer is then simply 9 = Ft. This extreme corresponds to the case where the temperature is high, the growth flux low or the ES barrier is strongly negative. This makes the descending step not only transparent, but attractive, since all atoms that land on a higher terrace funnel down to the lowest terrace on the surface. Nondiffusive and layer-by-layer growth represent two extremes for layer-filling during epitaxial growth. In Figure 5.5 we show the' layer filling ac-cording to SOS simulations, using the standard .parameters. In part (a), low-temperature growth is simulated, and in part (b), the temperature is set to 625°C, close to that typically used in GaAs growth. The insets show the actual surface after some 5 seconds of growth. Both simulations show a layer-filling behaviour that is very close to that found in the layer-by-layer growth model. Chapter 5. Modeling of growth on patterned surfaces 79 0 " HE /////// /////// 0 1 2 3 4 5 6 0 1 2 3 4 5 6 Time (s) Time (s) Figure 5.5: Layer filling for SOS simulation on an initially flat starting surface for growth at (a) 450°C and (b) 625°C. Standard parameters were used. The insets in each part show 75x150 atom sections from the simulations at t = 5.23 s. 5.4.2 Evolution of the step density In Section 5.3, the early stages of growth within one monolayer was described as an evolution of islands of all sizes. Earlier in this section, all the islands in a given layer were summed into a single parameter as the filling of that layer, 9j. If we for the time being assume that the surface grows like the layers of a circular wedding-cake, then we can define the number of steps within one layer as being equal to the length of the perimeter of the island, Sj = 2nrj, where Tj = y/Oj/n is the radius of the (single) island that has incorporated on the jth layer. We then get that Sj oc y/Q]. As the islands grow, coalescence of islands become a factor, and we can eventually expect that the step density will decrease as the filling exceeds ~0.5 and proceeds towards a full layer: Sj oc (1 — y/dj). A formula that incorporates both the initial and later filling of a layer was proposed by Cohen [CPP+89]: Sj oc PjiBj-x - 6j)q (5.7) where (p, q) € [0.5,1]. This expression was then used in a model that en-capsulates the basics of a multi-layer mean-field description [CPP +89]: h = F - W - « i - i ) + (% - 0J+I)<*;] (5-8) where — 6j) can be recognized as the area of the jth terrace, and aj Chapter 5. Modeling of growth on patterned surfaces 80 describes a perimeter effect that we define by: "  ai = irrt- (5-9) where A is a constant that determines the degree of interlayer transport present in the model. 5.4.3 Reflection high energy electron diffraction oscillations Reflection high-energy electron diffraction (RHEED) is a popular technique used in M B E for measuring surface roughness as well as reconstructions on the growing surface. The occurrence of intensity oscillations of specular R H E E D has been linked to the layer-by-layer filling of monolayers during growth [SJS+92, SJE + 93]. Experimentally, it was found that as initial sur-face roughness smoothed out, and the surface entered a layer-by-layer growth mode, the specular intensity would oscillate with a frequency that appeared similar to the growth rate, F. It has been shown that the intensity variations of the specular amplitude of the R H E E D signal follows [CPP+89]: m / o c ^ j - ^ + i ) 2 . (5-10) i=o where m represents the index for the uppermost level with a non-zero fill-ing. This expression is based on the notion that the step-edges scatters the R H E E D signal and hence reduces its specular intensity. Several groups have used the relationship in Eq. 5.10 in order to study R H E E D oscillations, and also the transition to step-flow growth which occurs at high temperatures and low Ga fluxes [NDJZ85, SVW+92, VC90]. We illustrate the relationship between the layer-filling and the intensity oscillations calculated from Eq. 5.10 in Figure 5.6. In part (a), we show the layer-filling from an SOS simulation using the standard parameters. The intensity oscillations in (b) are very similar to typical experimental findings for R H E E D oscillations. As the growth is stopped at 19.4 s in the figure, the intensity slowly increases, an effect typically referred to in the literature as "step-recovery". We saw this effect in Figure 5.2. Closer inspection of the signal Fig. 5.6 (b) indicates that at 15 seconds, less than 15 full oscillations are completed. The frequency of these intensity oscillations is often associated with the completion time of a layer, and hence is used to calibrate the growth rate: one oscillation per monolayer deposited. Chapter 5. Modeling of growth on patterned surfaces 81 If this were the case, then the signal should oscillate with a frequency of 1 Hz, which it clearly does not. For smooth growth, as indicated here, the oscillation frequency of the diffracted intensity is close, but not identical to the growth rate F. This is because as the surface is subjected to kinetic roughening, several layers are active at any given time. Therefore, a phase-lag is introduced into the layer-completion time, with the effect of decreasing the oscillation frequency of the specular R H E E D intensity. Recently, an alternate explanation for R H E E D oscillations has been sug-gested [BDP98]. In this model, the oscillations are believed to be caused by an interference effect in elastic multiple scattering of electrons between the top and bottom of the reconstructed adlayers. The decreased R H E E D oscillation frequency is consistent with this idea. In fact, there is no in-consistency between the two models: in one case, the R H E E D intensity is linked to the step-density, and in the other it is linked to the layer-filling. We showed in Eq. 5.7 that the step-density S is linked to the layer-filling 8, and we speculate that the two explanations for R H E E D oscillations are ultimately different formulations of the same idea. In these last two sections, we have introduced some mean-field mod-els that were shown to be useful in simple comparisons with experimental observations, such as the specular R H E E D oscillations. We showed how attachment and detachment could be incorporated in a simple fashion. The obvious weakness in the mean-field models is their lack of spatial informa-tion. This makes it difficult to use them for modeling of growth on vicinal Chapter 5. Modeling of growth on patterned surfaces 82 substrates, or on surfaces with some pre-defined roughness. This will be the topic of the next Section, where we will derive a growth equation from first principles of statistical motion of adatoms on a surface filled with available incorporation sites. 5.5 C o u p l e d g r o w t h e q u a t i o n s m o d e l The growth phenomena discussed in the previous sections can be explained in a natural way if we extend the growth model to include the adatom dynamics explicitly with two ID coupled growth equations ( C G E ) 2 : C G E : cVi + V - j = F-a^dtioT^h) (5.11a) dt{a-]}h) = aDn2 + (3DnS -al{KS (5.11b) Eq. 5.11a is a continuity equation for the adatom density n(x,t) (in units of n m - 1 ) with source and sink terms, while Eq. 5.11b describes the time dependence of the surface height h(x,t) (in nm), which depends on the dimer nucleation rate and the net adatom attachment rate at steps. The constants are defined as follows: F - deposition rate from the vapor (in 1/nm-s), D - adatom diffusion coefficient (in nm 2 /s), S(x,t) - density of steps (in 1/nm), K - rate of thermal detachment of atoms from step edges into the adatom phase (in 1/s), a and /5 are the incorporation coefficients for an adatom to another adatom or a step edge site, respectively. The lattice parameters are ay (in-plane) and a± in the growth direction. A n adatom is defined as a diffusing unit on the surface, which could be a Ga atom or a Ga-As complex. Equations 5.11a and 5.11b are continuum equations in the sense that'the variables, namely the surface height, adatom density and step density are macroscopic quantities averaged over a number of atomic units. In Eq. 5.11b, any adatom that attaches to a step edge is assumed to have incorporated into the film. We also define j = -D{a]}C,nVh + Vn) (5.12a) S = ^S2 + ( a ^ V / i ) 2 . (5.12b) where in Eq. 5.12a, j is the surface current of adatoms and £ is a propor-tionality constant which describes the net.directional drift of adatoms. A 2 Related coupled equations have been used earlier to describe growth on faceted sur-faces [KNT+94, BKT+01, HIWK90]. In this work the adatom incorporation rate was dependent on the crystal facet orientation. Chapter 5. Modeling of growth on patterned surfaces 83 positive value for £ gives a net downhill adatom current, consistent with the surface smoothing that is observed experimentally for GaAs (001) [CCOO, CCM+98, BRS+02]. At this point the form for the surface current in Eq. 5.12a can be regarded as a hypothesis, motivated by the success of the second order linear term in the K P Z equation in fitting the GaAs sur-face morphology in the limit of low surface slopes [BRA +01]. The second term in Eq. 5.12a represents adatom diffusion. Equation 5.12b is a physically plausible expression for the dependence of the rms step density on the surface slope. In this equation So is the random step density on a flat surface due to growth related phenomena such as island nucleation and thermally induced disorder. This step density will depend on temperature, growth rate, arsenic flux and on time, if the growth rate is not constant (see Fig. 5.2 (a) and (b)) [BJK +03]. If the random step density is not too large, one can define a random local slope associated with the local configuration of the steps. In the presence of a macroscopic topography the random local slope associated with So adds to the macroscopic slope V / i . If the two slopes are uncorrelated then the rms step density is given by Eq. 5.12b. In general there may be correlations between So and V / i . Nevertheless, one might expect this expression for the step density to be relatively insensitive to correlations, as it has the correct limiting behavior for large and small surface slopes. Therefore one could also regard Eq. 5.12b as a convenient interpolation formula. In the limit that Vh < So there will be numerous up and down steps (i.e. small islands) between successive net upward (or net downward) steps associated with the macroscopic surface slope. The net downward flux asso-ciated with the step edge ES barriers or other mechanisms which drive the flow of adatoms downhill will depend on the macroscopic average surface slope; For example, a monolayer island located on a terrace on a vicinal surface, will hot cause a net macroscopic downhill flux of adatoms, even in the presence of ES barriers as the down flow on one side of the island will balance the down flow in the opposite direction on the other side. For this reason the adatom current in Eq. 5.12a is proportional to the macroscopic surface slope Vh. The simple picture of a surface consisting of flat terraces separated by atomic steps, can be expected to provide a good description of the surface as long as the surface slope does not reach the next low index crystal planes, namely (110) and (111). These planes are 45° and 54.7° from the surface normal, and are beyond the range of surface slopes that we have explored experimentally (less than ~ 30°) . Chapter 5. Modeling of growth on patterned surfaces 84 5.5.1 Notes on the solution of the coupled growth equation model Only a few partial differential equations (PDEs) can be solved analytically, and only the simplest ones at that, and it is common to solve more complex PDEs with numerical approximation techniques, for example by the finite difference approach. The C G E model consists of two coupled, nonlinear, stiff partial differential equations, to which no analytical solution exists. In this Section, we show how we would approach the solution of a simple P D E using numerical methods, and outline the extension to the solution of the C G E model with a i r its complexities. We recall the Taylor formula for a general function f(x): f(x ± Ax) = f(x) ± Axf'(x) + ^rf"(x) + ...' (5.13) We can easily verify the following approximations: /'(x) = /(* + ^ ) - / ( * - A x ) + 0 ( A x ) { 5 1 4 ) /»(x) = /(* + A * ) - 2 J M + / ( x - A x ) + 0 ( A i 2 ) ( M 5 ) where O(-) indicates the order of the truncation error of the approxima-tion scheme. Applying these techniques, we can write the space-discretized version of the heat equation as dtf = V2f(x,t) dtfi ^ 3 ^ {fi+i - 2/i + / U ) (5.16) where / , = f(iAx). Al l of the quantities on the right hand side are now known, and we are left with a set of coupled ordinary differential equations (ODEs), one for each spatial grid point, x = iAx. A further complication inherent in the C G E model and its underlying physics is the great spread in timescales of interest: while the diffusion of adatoms operates on the order of microseconds, the evolution of the large scale morphology of the surface is on the order of minutes to hours. This translates directly to the mathematical nature of the C G E model in that it is a stiff model, and more sophisticated techniques are required in order to solve their time evolution in a reasonable amount of calculation time. To shed some light on this issue, consider the adatom diffusion term that occurs by combining Eqs. 5.11a and 5.12a: <9jn = D V 2 n , commonly referred Chapter 5. Modeling of growth on patterned surfaces 85 to as the heat equation. The maximum time step that the conventional forward Euler time stepping scheme will allow for this equation is Atmax = Ax2/2D. The diffusion parameter D can easily be on the order of 105 nm 2 / s or more. With a grid spacing of typically 10 nm, that gives a max time step of less than a millisecond. Comparing this result to for example the E W equation, which has the same functional form as the heat equation, but a prefactor no larger than v ~10nm 2 /s , we can easily get a time step on the order of tens of seconds. The presence of the adatom diffusion term in the C G E model therefore requires more advanced time-stepping algorithms. We chose to use Matlab's implementation [SR97, SRK99] of the "Liv-ermore Solver for Ordinary Differential and Algebraic equations" (LSODA) [RH93], one of nine solvers in O D E P A C K [Hin83], which is freely available from the Lawrence Livermore National Laboratories website3. LSODA's strength lies in its adaptive time-stepping abilities, and hence its abilities to solve stiff ODEs. By combining the use of L S O D A with the decoupling of the spatial grid points that we did in Eq. 5.16, the solution of the C G E model was straight forward, and would typically be completed in a few seconds using Matlab on a typical desktop computer. 5.6 Coupled equations in the asymptotic limit Before we go on to calculating the evolution of large scale surface shapes according to our C G E model, we will try and reduce it to a single continuum growth equation. The first and simplest case assumes no adatom nucleation, n 2 ~ 0. This approximation should not affect our calculations very strongly, as we aim to describe the morphology of epitaxial growth on surfaces with large slopes, such that the incorporation of adatoms into the film is mainly mediated through the large amount of steps on the surface, rather than the nucleation of dimer islands. Also, one should keep in mind that although the C G E model is fully deterministic in its current format, there will still be incorporation of adatoms to small scale surface structures through the thermally generated background steps, represented by the variable SQ. The typical choice of initial condition for the adatom concentration is n(x, 0) = 0. As soon as the growth begins, the adatom concentration quickly builds to a constant level that only varies slightly in time and space, n(x,t > 0) ~ n + 5n. For a sufficiently small variation 5n, we can make the approximation c^n ~ dtn = 0. For low amplitudes and long wavelengths, 3The Center for Applied Scientific Computing provides free of charge Alan Hindmarsh's ODEPACK at http://www.llnl.gov/CASC/odepack/ Chapter 5. Modeling of growth on patterned surfaces 86 the large-scale surface morphology is weak, or V/t < So, so we also assume that the spatial variation of n is small. In this case, V n ~ V 2 n ~ 0, and of course S ~ So. Furthermore, we assume for now that the product n V 2 / i is small. We incorporate these simplifications into Eqs. 5.11, 5.12 and get 0 = F-a^dtia^h) = F-a^tfDnSo-awKSo) t A = ^ F + KSo) (5.17) In this case, we can reduce Eqs. 5.11a, 5.11b to a^dtia^h) = F + Da7}tnV2h b\h = a[{a±F + uV2h (5.18) where a 2C / F \ P \S0 This reproduces the E W equation, and we see the explicit dependency of the linear smoothing coefficient v on the deposition rate and the downhill drift parameter £. In addition, it shows that in the absence of growth (F = 0) the linear smoothing term is independent of the background step density So- This agrees with the light scattering data in Fig. 5.1, which shows that the smoothing rate is relatively constant during a growth interruption even though the A F M images in Fig. 5.2 indicate that the step density drops dramatically during annealing. We showed in Chapter 4 that the presence of the E W term was well justified experimentally. We would now like to extend our analysis to include the next higher order terms in the surface gradients. We use a perturbation method devised by Bayo Lau 4 , similar to that used in quantum mechanics [Sch68]. Writing the adatom concentration as n = no + en\ + e2ri2 + • • • (5.20) we look for a growth equation in h from the perturbation 3th = a t (/ i 0 + e k i + . . . ) ( 5 - 2 1 ) 4This work was performed as part of a summer term project investigating analytical and numerical aspects of the CGE equations. His report " Brief Report on Coupled Growth Equation Investigation" is available on request. Chapter 5. Modeling of growth on patterned surfaces 87 We insert the expression for n from Eq. 5.20 into Eq. 5.11b and get: dtia^h) = (aDnl +(3Dn0S-a\\KS) +e{Dm(2an0 + (3S)) + . . . (5.22) We now insert the formulation for the surface current (Eq. 5.12a) into Eq. 5.11a, while letting dtn —> 0, and find that: oo dt[aj}h) = ay (F + a^D^hn) + ^ e^^ i - i (5.23) t=i where the operator Q is defined as g = : a | | L > ( a J 1 C V W + V 2 ) (5.24) We get n-i from equating the ith order terms in e from Eqs. 5.22 and 5.23. To the zeroth order, this perturbation method gives the following expression for n o : 0 = aDnl + D((3S-a\\a^(y2h)n0-a\\(F + KS) t ^ • SoP ( L , ^ ( F + KSQ) \ A S 0 a N ° " PDS$ -l)=w{T~1) (5-25) When P2DSQ » 2aa\\(F + KSQ), which is the limit in which nucleation of dimer islands is not important, n o reduces to the expression for n as defined in Eq. 5.17. We proceed to the first order in the perturbation approach. From Eqs. 5.21 and 5.23 we see that dthi can be obtained from a±Qno. After a little bit of algebra, we arrive at the following expansion ^ = - ^ p V 2 ( V / i ) 2 - nV4h + A i ; 3 V ( V n ) 3 + O (V*(V/iy) (5.26) where i and j are integers and their sum exceeds 4. The prefactor for the Villain term is defined by A 2 , 2 = J D ( I - l ) ( a | | C 2 - ( l + i ) | ) - 2 ( l + I ) ^ (5.27) where T was defined in Eq. 5.25. We can also here use the approximation that (3DS2 » 2ay(F + KSQ) (ignoring dimer nucleation), and thereby find a simpler expression for the Villain prefactor: ^ ( J r f ( £ - ) ) Chapter 5. Modeling of growth on patterned surfaces 88 The sign of this prefactor changes when the growth flux goes to zero. The prefactor for the Mullins term V 4 / i is which is negative for typical values of D, ( and T (see Table 5.1 in Sec-tion 5.7.1). This means that its contribution in the growth equation is a destabilizing one. Therefore, we cannot attribute this term to the diffusion of surface adatoms, like Mullins did. We speculate that by further expand-ing the C G E model, higher order terms will emerge that stabilizes this term. We discussed such a scenario in Section 3.2 with respect to the KS equation. We see that the Das Sarma term V(Vfa) 3 also emerges from the pertur-bation method, and we showed in Section 3.6 that this term was dominant in the absence of the E W term V 2 / i in the hydrodynamic limit. We have not attempted to work out the prefactor A i ^ , and further analytical in-vestigation is required into the nature of this term, as well as higher order terms. The lowest order nonlinear growth equation that emerges from the per-turbation of the C G E model is therefore where all the prefactors are positive under typical growth conditions, except at high T and low F, where A2}2 is found to change sign. It should be pointed out explicitly that in Eq. 5.30, we have dropped the terms of order e2, as well as terms from hi that contained higher order derivatives, like V 4 ( V / i ) 2 and the like. We notice immediately that the adatom current proportionality factor £ is multiplicative in the linear term, so that in the case of a zero ES barrier, there is no linear smoothing during growth. The Villain term, however, still contributes in the case of a net zero adatom current. In Fig. 5.7 we compare a simulation of Eq. 5.30 to a homoepitaxy ex-periment. In the experiment, the surface was prepared with a 10 minute O 2 U V treatment, and followed by a thermal desorption at 620°C. A buffer layer was then grown at T = 550°C for 90 minutes, yielding a 1.2 u,m thick film. The simulation only includes the first two terms in Eq. 5.30, and the. parameters used were vx = 10nm2/s, vy = \nm?/s and A2 ] 2 ,x = &-2,2,y = 106 nm 3 / s . Like the simulations performed in Chap-ter 4, the added growth noise used in these simulations was non-conservative. Fig. 5.7 (a) shows the surface after growth, with the typical weakly elongated (5.29) b\h = a\\axF + uV2h - %?-V 2 (V/ i ) 2 (5.30) Chapter 5. Modeling of growth on patterned surfaces 89 Wavevector q Gun" ) Figure 5.7: Comparison of simulation based on reduced C G E equation (Eq. 5.30) and regrowth experiment [Sch]. (a) shows actual surface after 90 minutes growth; (b) shows the corresponding simulation and (c) shows the PSD of both, along with that of the starting surface. mounds along the [110] direction. Part (b) shows the corresponding simu-lated surface for an initial condition similar to that shown in Fig. 4.2 (a). The length scale of the mounds is in good agreement with the surface in (a). Part (c) shows the PSD along the [110] direction, and again excellent agreement is found between the experimental results and the simulation. Our work in Chapter 4 and Ref. [BRS+02], showed that surface shape comparisons, as well as comparisons of the power spectra from the K P Z simulations and experimental data showed excellent agreement. However, in Chapter 5. Modeling of growth on patterned surfaces 90 Fig. 5.7 we see that by replacing the K P Z nonlinearity with the Villain non-linearity, the resulting simulated surface shapes and PSDs also show good agreement with experimental results. We conclude that the K P Z and the Villain nonlinearities have similar effects on the surface shape. Fig. 3.5 (b) and (e) showed that the effects of the terms (V/i) 2 and — V 2 ( V / i ) 2 on the surface evolution are very similar. We conclude that the spatial frequency range accessible to us experi-mentally is not large enough to allow us to distinguish the two different nonlinear terms in the growth equation. If we compare the magnitude for these two nonlinear terms at the dominant spatial frequency in Fig. 5.7 (c) (approximately 4 yum), we find that they are approximately equal. This can be done by considering the two nonlinear terms in Fourier space. The fourth order nonlinearity becomes: A 2 ,2 {-iqfF [{Vh)2] (5.31) Comparing this result to the Fourier transform of the K P Z nonlinearity, ^[(Xxpz/^iVh)2}, we arrive at: <72A2,2 = XKPZ (5.32) Substituting in values for the dominant spatial frequency along with the parameters used for the simulations in Figure 5.7, we get: QmaxA2,2 = (4/^m)2 x (106 nm3/s) = 16 nm/s (5.33) which compares well to the value used for XKPZ under similar growth con-ditions in Chapter 4 (see Figure 4.3), namely Ylnrajs. This shows that the magnitude of the K P Z nonlinearity obtained from the K P Z fits to experimental data is similar to the magnitude of the Villain nonlinearity inferred from the parameters used in the coupled growth equa-tions. Although the g-dependence of the two non-linearities is different, the range of spatial frequencies accessible experimentally is too small to distin-guish them. This result provides a natural explanation for why A ^ F in the K P Z fits, as discussed earlier. The nonlinearity is not due to growth along the surface normal but rather to the more complex growth process discussed by Eqs. 5.11a and 5.11b. While the randomly roughened substrates that we investigated in Chap-ter 4 show clear sign of nonlinear shape evolution, they are not rough enough for the nonlinear term to dominate the growth equation. When the non-linear term is dominant, the surface is said to be in the strong coupling regime [NB96]. Chapter 5. Modeling of growth on patterned surfaces 91 30 60 100 200 300 Spatial frequency, q (1/um) Figure 5.8: SOS simulations showing strong coupling transition for increasing aspect ratio starting surfaces (a-c). Part (d) shows the smoothing rate dependence on spatial frequency in the linear regime depicted in (a). The standard parameters were used, and the growth temperature was 575° C. In order to explore the relative importance of the nonlinear term, we have performed 2D K M C simulations using the SOS model with standard param-eters on a ID sinusoidal starting surfaces with varying pitch and amplitude (see Fig. 5.8). In part (a), a pitch of 600 atoms was used, and the initial peak-to-peak amplitude was 8 atoms. This corresponds to a steepest surface slope of 2.4°. This surface is found to evolve linearly, with each successive growth scanline being very similar to the starting shape on the bottom, but with a slowly decaying amplitude. Part (b) shows the evolution of a similar Chapter 5. Modeling of growth on patterned surfaces 92 surface with half the pitch, and hence a steepest surface slope of 4.8°. We see that this is steep enough to warrant nonlinear surface shape evolution, and that the successive growth scanlines evolve into rounded mounds sepa-rated by cusps, a growth behaviour typical of the K P Z equation. The fact that a shorter pitch at the same amplitude leads to more nonlinear growth suggests that the Villain nonlinearity is operative. In the case of the K P Z equation, the linear and nonlinear terms have the same g-dependence. Fig. 5.8 (c) shows a taller grating, with a 16 atom peak-to-peak ampli-tude, and a steepest slope of 9.5°. The surface shape evolution is definitely nonlinear, but does not exhibit K P Z behaviour until the uppermost scanline after 90 M L growth, where the amplitude has been sufficiently reduced. The transient surface shapes after 30 and 60 MLs have mounds and cusps, like those we see in (b), but in addition, the sidewalls have positive curvatures, something the K P Z equation cannot produce. We speculate that the Vil-lain nonlinearity is required in order to model such behaviour. In part (d), we. show the dependence of the smoothing rate on the spatial frequency of the grating for low aspect ratio surfaces, like the one displayed in (a). It is seen to vary as q2-2, close to the q2 dependency expected for the linear E W equation. In the next Section, we will perform regrowth on patterned substrates in the form of ID gratings with a characteristic spatial frequency of q ~ 4 / i m _ 1 , corresponding to pitches in the low micron range. These surfaces are found to evolve unmistakeably in a strong nonlinear fashion. 5.7 E v o l u t i o n o f p a t t e r n e d su r face g ra t i ngs 5.7.1 F i lm thickness evolution Growth on substrates with larger amplitude surface slopes, up to ~ 30°, show complex surface shapes before evolving into parabolic mounds, as shown in Fig. 5.9 (a). At intermediate times, the valleys are V-shaped with concave rather than convex sidewalls and distinct shoulders near the top of the sidewalls. Note the absence of (001) facets which are predicted theo-retically for annealing below the roughening temperature in the absence of deposition [SRF03]. Equations 5.11 and 5.12 can be solved in seconds with a finite difference scheme and a coupled differential-algebraic system solver, and a ID solution is shown in Fig. 5.9 (b) with parameters adjusted to match the experimental data in Fig. 5.9 (a) (see Table 5.1 for parameters). The agreement with the experimental surface shapes is striking. In particular, the model reproduces the inverted " Gothic window" shape of the valley for Chapter 5. Modeling of growth on patterned surfaces 93 <D . C <D O t CO 1 1 1 ( a ) 2600nm (b) 2000nm 600nm 200nm regrowth Starting surface n n i i i -2 x (urn) Figure 5.9: Film thickness dependence: (a) A F M scan lines for regrowth on 100 nm deep gratings oriented perpendicular to the [110] direction (substrate temperature was 580°C, F = l M L / s and the As:Ga pressure was 3:1); (b) Scan lines from C G E calculation; (c) Scan lines from 2D K M C simulation of 10 nm high grating structure, where one At equals 5.6 M L of growth. A l l offsets arbitrary. the 600 nm growth and the KPZ-like cusps in the 2600 nm growth where the grating amplitude has reduced significantly. A continuum model cannot include the microscopic details of the atomic scale phenomena, such as the geometrical configuration of the step edges. We therefore compare the continuum model described by Eqs. 5.11 and 5.12 with K M C simulations of the SOS model outlined in Ch. 2, using the standard parameters and a negative ES barrier. K M C simulations produce a random step distribution automatically due to the statistical nature of the model. In K M C , the binding energy for an atom at a step edge depends on how many neighbors it has (~ mEiat), whereas in the C G E continuum model a single average value is used for the step edge binding energy. SOS simulations of M B E growth by K M C are limited by available com-Chapter 5. Modeling of growth on patterned surfaces 94 Table 5.1: Parameter table for C G E calculations. The gratings in Fig. 5.9 (b) are _ L [110] and in Fig. 5.10 (b) _!_ [110]. The numerical values for a match those of 0 in our simulations. Figure F T D /10 5 K So /3/10- 4 C / i o - 3 (ML /s) (°C) (nm 2/s) (Vs) (1/nm) 5.9 (B) 1.0 580 2.5 4.3 0.07 3.3 1.4 5.10 (B) 0.8 420 0.04 0.004 0.49 6.6 1.5 500 0.4 0.2 0.38 1.6 1.4 550 1.3 1.4 0.19 0.4 1.2 610 4.6 12 0.03 13 1.1 puting power to small scale structures, and become intractable for realistic, high temperature growth scenarios where 2D systems have sides up to mi-crons and growth times on the order of hours. In Fig. 5.9 (c), we show a K M C simulation for a surface, grating that is a fraction of the size of the experimen-tal structures. The simulated grating profiles in Fig. 5.9 (c) were obtained by projecting 2D K M C simulations onto a line at each time point by taking the average elevation perpendicular to the scan line. The agreement with the experimental shapes is excellent, reproducing all of the main features, except they are on a smaller size scale. The substrate and lateral binding energies are similar to values reported earlier in the interpretation of reflection high energy electron diffraction (RHEED) data [Joy03, SVW+92] and compati-ble with the fitting parameters found in the continuum model. Kratzer et al. [KMS98] found from density-functional theory calculations that the diffu-sional barriers are anisotropic, with £>[no].= (0-1 cm2/s)exp [—1.5eV/ksT] and -D[iio] = (0.07 cm2/s) exp [—1.2 ey/fcgT], i.e. faster diffusion along the [110] direction. It is plausible that similar shapes could be obtained for the larger size scales relevant to the experiments by scaling the parameters appropriately. In the case of the C G E model (Eqs. 5.11, 5.12) we find that the parameters can indeed be scaled to yield similar surface shapes at different length scales (see Appendix D). 5.7.2 Temperature evolution In Fig. 5.10 (a), we show the dependence of the surface topography on growth temperature, for a fixed layer thickness together with (b) the simulated surface topography using Eqs. 5.11 and 5.12 and parameters as outlined Chapter 5. Modeling of growth on patterned surfaces 95 \ / \ / 5 0 0 ° c \ / \ / \ / 550° c \ / Start ing Su r f a c e i i i Start ing S u r f a c e 0 5 10 Su r f a c e Coo rd i n a t e x (u.m) Figure 5.10: Growth temperature dependence: (A) A F M scan lines for regrowth on 100 nm deep gratings oriented perpendicular to the [110] direction; (B) C G E calculation. The grating pitch is 5 fim. A l l growths are 1 hour at 0.8 M L / s and the As:Ga pressure was 3:1. Al l offsets arbitrary. in Table 5.1. The experimental data is obtained from growths on 100 nm deep gratings oriented perpendicular to [110]. This is the faster diffusion direction in this material system [BRS+02]. The anisotropy in the diffusion depends on the A s 2 / G a ratio during growth, which was equal to three for the data shown in Figs. 5.9 and 5.10. This observation is consistent with the values used for the downhill drift parameter in our calculations, where the best fits were obtained using a larger £ when the gratings were oriented perpendicular to [110] (Fig. 5.10 (b)) than for gratings perpendicular to the [110] direction (Fig. 5.9.(b)). There is some uncontrolled variation in the pitch and depth of the gratings in the experimental data in Fig. 5.10 (a). With that in mind, we see that the C G E model reproduces the main features in the temperature dependent data, namely the small secondary mound in the valley at 420°C and 500°C, the KPZ-like cusp at 550°C and the inverse Gothic window shape for the valleys at 610°C. The shoulders at the edges of the ridges at 610°C are also reproduced by the model, although they are Chapter 5. Modeling of growth on patterned surfaces 96 not as sharp as in the experimental data. The parameters used in these calculations are based on the same energies used in the K M C simulations in Fig. 5.9 (c). The diffusion constant is related to the substrate binding energy by D = aj(2kT/h) exp (-E s u b /kT). The step edge detachment rate is calculated from K = D x exp (—rhEiat/kT), where m is an average number of neighbors for atoms at a step edge which we set equal to 2.5. The declining value of SQ with temperature is reason-able under growth conditions; one might also expect £ to decrease weakly with temperature. Satisfactory fits to the data were also obtained with a larger activation energy for D (1.8 eV rather than 1.25 eV) together with a smaller prefactor and somewhat different (yet still physically reasonable) temperature dependences for the other parameters. Experimental and the-oretical work suggests that the activation energy for D is in the 1.2-2.0 eV range [SVW+92, Joy03, YLB+99, KPS02]. 5.7.3 Transient amplitude overshoot In a pair of recent publications, Kan et al. [KSTEP04b] and Shah et al. [SGL+03] (see also Refs. [BTS04a, KSTEP04a]) investigated the shape evo-lution of patterned GaAs surfaces during epitaxial growth. The starting shapes in their studies consisted of circular holes or pits ~50nm deep, or-dered in a square array across the surface. Using the diameter of the holes as a control parameter, homoepitaxial regrowth experiments were then per-formed on the patterned surfaces and the grown surfaces were imaged at various growth times with A F M . The diameter of the holes was equal to one half their centre-to-centre spacing. Kan and Shah observed that the peak-to-valley amplitude of their pat-terned surfaces underwent a non-monotonic decay: the amplitude of the surfaces with pits larger than about 1 pirn diameter would first increase be-fore declining, whereas the amplitude of the smaller diameter holes would decrease immediately upon commencement of growth. This effect is also present in our regrowth experiments on ID gratings. For example in Fig. 5.9 (a) the data shows that the surface peak-to-valley amplitude has in-creased after 200 and 600 nm film growth, but has decreased again at longer growth times. Figure 5.10 (a) shows the same effect, where the data shows an increase in the surface peak-to-valley amplitude, except possibly in the case of growth at 620°C, where the amplitude is about the same as the starting surface. In Fig. 5.11 we show the temporal evolution of the peak-to-valley am-plitude of a ID grating according to the C G E model at T = 5 8 0 ° C , with Chapter 5. Modeling of growth on patterned surfaces 97 l.1-2h « 1.1 y Q . i 0.9^  °-0.7r .§ 0.6 -o 0 0.5 1 1.5 2 2.5 Time (hr) Figure 5.11: Temporal peak-to-valley amplitude evolution according to C G E model. Ex-perimental amplitudes from regrowth experiments shown in Fig. 5.9 (a) are superimposed. parameters as outlined in Table 5.1. Cross-sectional linescans from this cal-culation were shown in Fig. 5.9 (b). The.sidewalls of the starting surface are set to 30° in order to match the experimental starting surfaces. The mea-sured values for the amplitudes of the experimental gratings from Fig. 5.9 (a) are also plotted in Fig. 5.11, and found to follow the trend of the calculation very well. It is interesting to note that a.stable.model, with a net-downhill drift of adatoms can create a characteristic of unstable growth, namely the peak-to-valley amplitude overshoot. In Figure 5.12, we investigate the full amplitude overshoot dependence on growth time and pit diameter. Part (b) shows the systematic experimen-tal investigation done in Refs. [KSTEP04b, SGL+03], where the pit diam-eter, and hence the pit spacing, was varied in the range 700 nm to 8 fj,m. The different lines reflect the peak-to-valley amplitudes for different pit di-ameters at four different growth times, as indicated in the legend. These data values were digitized from the original Fig. 2 (a) of Ref. [KSTEP04b]. The experimental uncertainties in these measurements were not reported. As the growth time increases, we find that the maximum amplitude over-shoot occurs for larger diameter holes. Figure 5.12 (a) shows the corre-sponding curves calculated from the C G E model, as originally published in Ref. [BTS04a]. The starting surfaces are indicated in the lowest curve in the inset in Fig. 5.12 (a). We tried to match the starting shapes of Kan's Chapter 5. Modeling of growth on patterned surfaces 98 Pit Diameter (u.m) Figure 5.12: Peak-to-valley amplitude evolution for corrugated surfaces: (a) ID calcula-tions based on C G E model, and the inset shows actual surface evolution; (b) experimental data from regrowth on pitted starting surface (from Kan et al. [KSTEP04b]). pits, with a fixed sidewall slope of 30°. The remaining curves in the inset show the evolution of a typical grating according to our calculations. The parameters used in the C G E calculation were again the same as those used for the calculation of Fig. 5.9 (b). The agreement between Kan's data and our calculations are very good. The surface profiles shown in the inset of Fig. 5.12 (a) suggest an expla-nation for the non-monotonic decrease in surface amplitude. The sidewalls expand laterally, which means that they grow faster vertically than the ridge tops and valleys. This is due to adatom attachment at the high density of steps on the sidewalls. As the width of the valley shrinks, the adatoms in the valley have a higher probability of migrating to the sidewalls than the adatoms that are deposited on top of the ridges which are now wider. This means the ridge tops grow faster than the valley bottoms, leading to a net increase in peak amplitude, as observed experimentally by Kan et al. A n attempt was made in Ref. [KSTEP04b] at modelling the evolution Chapter 5. Modeling of growth on patterned surfaces 99 of these patterned surfaces using conventional continuum growth equations, including the K P Z and the M B E equations, without success: neither the surface shape evolution nor the amplitude evolution could be reproduced. In fact, according to numerical calculations in Ref. [KSTEP04b] using the K P Z and the M B E equations, the amplitude evolution of their simulated surfaces was monotonically decreasing without the initial overshoot. Kan et al. concluded that none of the growth equations investigated could success-fully describe the surface evolution of GaAs during epitaxial growth. The continuum growth equations considered by Kan are asymptotic models, ap-plicable for long times and low surface slopes (V/i <C 1), and therefore it is not surprising that these models do not fit the relatively large amplitude surface topography studied experimentally by Kan and Shah. We showed in Section 5.6 that the nonlinear term that applies to our growth system is in fact V 2 ( V / i ) 2 . Both the K P Z equation and our reduced growth equation (5.30) produce very similar surface shape evolutions in the low amplitude limit. This is evident when comparing simulations based on the K P Z model (Figs. 4.3 and 4.5) and the C G E model (Fig. 5.7 (b) and (c)) to experimental A F M images (Figs. 4.2 and 5.7 (a), respectively). Our growth experiments and accompanying K P Z simulations in Chapter 4 and Ref. [BRS+02] provide sound experimental and computational evidence that K P Z provides an accurate description of the low amplitude surface morphological evolution for GaAs. As is seen in Figs. 5.9, 5.10 and 5.11, the full C G E model does a very good job at reproducing the surface shape evolution of large-scale, patterned surfaces, like Kan's pits and our gratings. It is also evident from the scanlines of the longer growth time gratings in Fig. 5.9 (a,b,c) that all of the experimental data, and calculations based on the C G E model and K M C simulations produce surface shapes that resemble typical K P Z mounds in the low amplitude limit. 5.8 Summary Our analysis in Chapter 4 lead us to investigate the smoothing behaviour of randomly patterned surfaces closer. Light scattering experiments along with A F M images showed that the smoothing rate of the surface roughness must be dependent on the deposition rate F, and hence also the adatom concentration, n. A coupled growth equation model was devised, incorporating the in-terplay between the adatom concentration and .the surface height. In this model, the adatoms were deposited, and were allowed to diffuse, attach and Chapter 5. Modeling of growth on patterned surfaces 100 detach at steps, as well as nucleate with other adatoms. The step density was modeled by a simple spatially varying interpolation formula. A net adatom current factor, £, was also included, and calculations showed that this the net direction of the adatom current was downhill. Using a perturbation method, we reduced the coupled growth equations to a simple, nonlinear growth equation: dth = alia±F + vV2h - % ^ V 2 ( V / i ) 2 (5.34) When applying this formula to the randomly patterned surfaces that we investigated in Chapter 4, they were found to evolve in a fashion very similar to the simulations based on the K P Z model. The prefactors in Eq. 5.34, given by Eqs. 5.19 and 5.28, were found to depend on the growth flux, the ES parameter £, the thermally generated steps, as well as the step edge release rate, K. SOS simulations supported the findings that it is the Villain nonlinearity that is operative during growth, and not the K P Z nonlinearity. The coupled growth equation model explains the replication of litho-graphically patterned surfaces very well, both as a function of film thickness and growth temperature. Unexpected support for our model came in the form of an excellent match with published, yet unexplained data describing the peak-to-valley amplitude overshoot for a pitted surface during epitaxial growth. In conclusion, we have shown that the complex surface shapes which de-velop during epitaxial regrowth on patterned (001) GaAs substrates can be explained by the dynamics of the deposited adatoms, including step edge attachment and detachment, adatom diffusion, and a stable, down-hill adatom drift. Although we attribute the downhill drift to a negative Ehrlich-Schwobel barrier we cannot rule out the possibility that this effect is caused by some other mechanism, or a combination of other mechanisms, as discussed in Chapter 2. Our analysis has been specific to GaAs, and all experiments are done on GaAs, however, the concepts are generic and should be applicable to other material systems. Chapter 6. Conclusion 101 Chapter 6 Conclusion GaAs is a technologically important material system used to make devices such as microwave integrated circuits, light-emitting diodes and laser diodes. It is a model system for epitaxial studies, since its physical and electronic properties are generally well understood. Something that is not well un-derstood, however, is the evolving surface morphology of epitaxially grown films, and this has been the concern of this thesis. We expect that an understanding of this basic property of growing films and the underlying microscopic mechanisms that govern it will become increasingly important as device dimensions continue to shrink. Examples of this are found in the design and fabrication of distributed feedback lasers and photonic crystals, where lateral shape control is crucial. In Chapter 2, we described some of the basic concepts related to adatom dynamics on crystal surfaces, and we introduced the Burton-Cabrera-Frank theory. We showed that interlayer transport of adatoms between the differ-ent terraces was an important factor in determining the overall evolution of the surface morphology during growth and annealing. If the net direction of adatoms is uphill, then the growth is unstable, and mounds will form due to nucleation of new dimer-islands on top of currently existing islands. How-ever, if the adatoms trickle downhill, then the growth is stable, and rough surfaces will smooth during epitaxial growth. In this case, we were able to derive the stable Edwards-Wilkinson (EW) equation. A commonly discussed mechanism that governs interlayer transport is the Ehrlich-Schwobel (ES) step edge barrier. This barrier inhibits adatom diffusion over steps, and is a destabilizing mechanism that effectively intro-duces an uphill current of adatoms. A n example of a stabilizing mechanism is downhill funneling, where newly deposited atoms relax into a low energy position upon impact. We showed that it is the net direction of the adatom current that deter-mines the long-time outcome of a growth. If, for example, the ES barrier is the stronger step edge mechanism, then large-scale roughness in the form of mounds will grow unabatedly. No such tendencies in the evolving sur-face morphology were observed experimentally. The next possibility is that Chapter 6. Conclusion 102 the destabilizing effects are balanced by the smoothing effects and that at a given "magic slope", the net adatom current is zero. This results in pyramid-shapes on the growth surface, a result also not observed in experiments of GaAs homoepitaxy. The final possibility is that the smoothing mechanisms are stronger than the destabilizing mechanism for terraces of any width. In this case, rough surfaces will smooth until they are balanced by an rms in-terface width determined by' kinetic roughening. This option is consistent with all our experimental observations. Several groups have modeled GaAs homoepitaxy using a positive ES bar-rier. In these studies a stabilizing mechanism such as downhill funneling was required in order to get agreement between simulation results and reflective high energy electron diffraction intensity oscillations. We found that all of our experimental results could be reproduced to high accuracy using only one mechanism, namely an inverse (negative) ES barrier. Although similar results could have been obtained by a combination of a positive ES barrier and a stronger downhill funneling, we believe that the continued smoothing of the surface after growth is ended indicates that the ES barrier is indeed negative. This is because the local relaxation mechanisms, such as downhill funneling, are associated with growth only. Based on these findings, we challenged the notion that the experimen-tally observed mounds in GaAs form due to a positive ES barrier, and that epitaxy on this growth system is unstable. In Chapter 4, we reported ho-moepitaxial growth experiments on GaAs (001) that showed that it was the surface preparation technique that caused the observed mounding. A com-mon procedure prior to growth on GaAs is to thermally desorb the native oxide from the surface. This leaves a rough, pitted surface, with pits as deep as 40 nm. These pits, in conjunction with the nonlinear growth that is ob-served in this material system, causes the surface to go through a mounded stage as it smooths. This is therefore a transient smoothing phenomenon, and not an instability caused by a physical mechansism on an initially flat starting surface. When preparing the samples with a hydrogen-etch tech-nique, a much smoother starting surface was obtained, and in this case, no mounding was observed during growth. We implemented a solid-on-solid model and simulated growth on singu-lar and vicinal surfaces using both a positive and a negative value for the ES barrier. As expected, a simulation of growth on the singular surface showed an increasingly rough surface as growth progressed when a positive ES barrier was used. On the vicinal starting surface, however, our findings contradicted those by Villain [Vil91]. We found that a positive barrier had a destabilizing effect on the surface morphology also here and that large Chapter 6. Conclusion 103 mounds formed on the steep sections of the random starting surface. Again, a negative ES barrier equalized the terrace widths. These results are also consistent with the stable E W equation. We also showed that the surface morphological evolution of both the thermally desorbed and the hydrogen-etched starting surfaces were well modelled by the Kardar-Parisi-Zhang (KPZ) equation, but not by the molec-ular beam epitaxy (MBE) equation. The second order linear term in the K P Z equation describes the net downhill diffusion of adatoms across step edges. A problem with the K P Z nonlinearity is that it does not conserve mass, which is a nonphysical scenario in M B E growth. Another issue with the nonlinear term is that it cannot be a manifestation of growth along the surface normal, because the prefactor A should equal the growth rate F and in our simulations the value of A that matches the data was on the order of 50F. We also questioned how one could relate the K P Z coefficients v and A to real physical parameters, like F, T and the As overpressure, as well as microscopic physical mechanisms on the surface. We have presented a new formulation for surface.growth based on two coupled partial differential equations (GGE). In this model, we took into account the spatial and temporal evolution of both the adatom concentra-tion and the surface height. We expect this model to be applicable to other material systems, as well. The new model was found to be a good descrip-tion of the evolution of randomly textured, nearly flat surfaces as well as lithographically patterned surfaces. For low surface slopes, the new model reduces to a KPZ-like equation, except with a fourth order, conservative nonlinear term. The similarity between the results obtained with the C G E model, and the K P Z model confirmed that the two conservative and noncon-servative nonlinear terms behave very similarly for the relevant experimental parameters. The C G E model was also able to predict the complex shape evolution of lithographically patterned gratings with sidewalls as steep as 30°. This is well beyond the range which asymptotic continuum growth equations such as the K P Z are applicable. The model was also able to re-produce the temporal evolution of the peak-to-valley amplitude reported recently by another group in regrowth on lithographically patterned GaAs surfaces. It is clear that the model accounts well for relatively disparate sce-narios designed independently by different groups. This said, there is still much to be learned about the fundamental physical mechanisms that are at play during epitaxial growth, not only in GaAs, but also in other material systems. We observed experimentally several interesting features that were not in-cluded in our modeling. Specifically, these were the surface feature anisotropy Chapter 6. Conclusion 104 associated with the As:Ga ratio used during growth and the fingering per-pendicular to the step direction, believed to be a manifestation of the Bales-Zangwill instability. We have developed a model that reproduces the surface shape evolution of GaAs (001) during epitaxial growth. One of the reasons for our success is that we have continuously attempted to match experimental results to the findings obtained by simulations. Continued investigation into this problem will have positive consequences for our fundamental understanding of the semiconductor surface. It will also benefit the device manufacturer, as an increased refinement is required in design on short length-scales. Appendix A. The kinetic Monte Carlo algorithm in solid-on-solid models 105 Appendix A The kinetic Monte Carlo algorithm in solid-on-solid models The kinetic Monte Carlo (KMC) algorithm is particularly well suited for the simulation of epitaxial growth. This algorithm is applicable to externally driven systems and is based on the work by Metropolis et al. [MRRT53] on equilibrium systems. Events on the surface, like diffusion and step edge de-tachment, are picked randomly from a distribution weighted by their rate of occurrence. This approach is best illustrated by a simple example. Imagine a system where only two events are possible (events 1 and 2) and one of these events are twice as likely to occur as the other: P( l ) = 2P(2). In other words, P(l) = 2/3 and P(2) = 1/3, since the sum of the probabilities must equal 1. According to the K M C algorithm, one would now pick a ran-dom number, r, from a uniform distribution between 0 and 1. If r < 2/3, event 1 is chosen and the system is updated. Otherwise, one proceeds with event 2. From the central limit theorem, we know that as the total number of events N increase, the number of type 1 events N± approaches P(1)N. The use of the K M C algorithm for implementation of an SOS model is well documented [JOH+94, vV93, Var02, NB99]. Our system is a 2D surface specified in an L x L matrix H with pe-riodic boundaries, where the matrix entries describe the surface height at the (i,j)th lattice site. During every step of the growth, we keep track of the number of sites J\fm that have a surface atom surrounded by m lateral neighbouring atoms. In the SOS model that we introduced in Chapter 2, we use a cubic lattice, and only the nearest neighbours are accounted for (m € {0,4}). The microscopic process in our system include deposition and diffusion, as well as detachment from any surface-site. This makes for a total of 6 different types of events that can take place: deposition at rate F (in ML/s ) , as well as diffusion of atoms with m neighbours at rates km. These rates follow Arrhenius behaviour, and are calculated accordingly to Appendix A. The kinetic Monte Carlo algorithm in solid-on-solid models 106 km = J^oexp [-fiEoct], with an activation energy E a c t = E s u b + mEiaU where Esub is the binding energy to the substrate and Eiat is the lateral binding energy. The prefactor VQ used in this work is derived from the equipartition theorem: VQ = 2kT/h, where h is Planck's constant. Deposition is assumed uniform across the surface. The K M C algorithm repeats and time is incre-mented [MG88, Wil89] until the desired growth time is reached. Due to the temperature dependence of the rates km, the high-temperature simulations will require the calculation of more events than a lower-temperature growth in order to get to a certain film-thickness. We outline the K M C approach according to Vardavas [Var02]. Step I. The first step is to calculate the total rate of activity on the surface based on the possible physical microscopic processes that can take place: 4 Ktot = FL2 + ^ 2 kmNm (A.l). 771=0 This allows us to associate a probability for a deposition event P(dep) and the diffusion events P(m): . P(dep) = FL2/Ktot (A.2a) P{m) = kmNm/ntot (A.2b) We continue by calculating the cumulative probabilites: Co - 0, Cm+1 = Cm + P(m), m e {0,4} (A.3) A random number is now generated from a uniform distribution r i [0,1) and an m-event is chosen if: C m < r i [ 0 , l ) < C m + i A deposition event is chosen if r i [0,1) > C5. Step II. Picking the type of event to execute is the easy part in K M C ; finding a site to which we can apply that given transition is harder. In the case that a deposition event is chosen, a site (i, j) is chosen randomly, and the surface height is incremented by 1 atomic unit at that site: (A.4) H ( i , j ) - > H ( i , j ) + l (A.S) Appendix A. The kinetic Monte Carlo algorithm in solid-on-solid models 107 When the chosen event is not deposition, but rather the diffusion of a surface atom with m neighbours, we need to find an atom that is in such a configuration. Searching randomly around the surface matrix for such an atom can be done, but at a high computational cost. The reason for this is twofold: first, the atoms with the fewest neighbours are the ones that are most likely to diffuse, and secondly, the fast-diffusing atoms are scarse, as they most likely diffuse step edges in a relatively short time, where they are less likely to detach and continue diffusing. Our approach to finding these rare atoms involves updating 4 binary trees1 Tm (m G {0, 3}) containing the coordinates of all atoms that have m neighbours. Given an event m, an appropri-ate set of coordinates is then randomly chosen out of the binary tree Tm. For the atoms that have 4 neighbours, a corresponding matrix N was used, where the (i, j)th entry contains a "1" if that site has 4 neighbours, and a "0" otherwise. This proved to be the most cost-effective approach for our particular growth sys-tem, where even the very rough or textured surfaces were found to have an abundance of 4-neighbour atoms. This way, an aver-age of less than two searches were required in order to find the appropriate coordinates. Now that we know what atom to hop, we randomly choose a site (i',f) which is within one jump of in the x- or y-direction only, and update H as follows: After the hop or deposition event has been performed, and the surface matrix has been updated, the quantities A/", Tm and N must be updated. Step III. The last step in the K M C algorithm is to increment the time, Atevent. In K M C , the time-step is not constant, as the microscopic events are independent of time, and follow a negative 1 Actually, AVL trees were used, where the computational cost of adding, deleting and searching the tree is less than that of a perfectly balanced binary tree. An AVL tree is a balanced binary search tree where the height of the two child subtrees differ by at most one, otherwise known as height-balanced. Look-up, insertion and deletion are all 0(log(n)) in both the average and worst cases. Additions and deletions may require the tree to be rebalanced by one or more tree rotations. The AVL tree is named after its inventors, Adelson-Velskii and Landis [AVL62]. ... H ( i , j ) ^ H ( v j ) - l H(i',j')->H(i',j') + l (A.6a) (A.6b) Appendix A. The kinetic Monte Carlo algorithm in solid-on-solid models 108 exponential distribution. It can therefore be shown that the time interval between successive events over the whole system also follows a negative exponential distribution with parameter TZtot-By producing another random variable r"2'[0,1), we can produce a time increment as follows: Atevent = — (A. 7) ><-tot Time can also be calculated by dividing the total coverage 8 on the surface by the growth rate, t — d/F, however, this approach will only increment the time whenever there is a deposition event. Appendix B. Details of the Burton-Cabrera-Frank calculations 109 Appendix B Details of the Burton-Cabrera-Prank calculations In Chapter 2, we introduced the B C F model for adatom density on vicinal surfaces. In this Appendix, we show some of the details of the calculations that were left out of that Chapter. B . l B o u n d a r i e s w i t h E h r l i c h - S c h w o b e l b a r r i e r s When there is no desorption of adatoms back into the vapour, the adatom density in ID obeys DV2p = -F/a\\ (B.l) whose solution regardless of boundary conditions is Fx2 • pW = -^p+c^x + C2 . (B-2) where c\ and c2 are constants to be determined from the boundary condi-tions. Including the ES barrier can be done by the following boundary condi-tions (Eqns. 2.13 and 2.14): x — 0 : —DVp = ——p, (Ascending step with ES barrier) (B.3) a x = I : —DVp = -—p. (Descending step with ES barrier) (B.4) IES • • Solving for (c\,c2) gives C l = _ZL (* + ***) (B.5) c 2 = anci (B.6) Appendix B. Details of the Burton-Cabrera-Frank calculations 110 so that ( \ - 4- F l & + a\\)(2lES + I) ( , P { X ) ~ 2a\\D + 2a\\D (l + a\\+lES) 1 j and the adatom current is j = -D (Vp(0) + Vp(0) (B.8) Fl - 2 D c x + — aj| F Z ( a n - lEs) = - x  — (B.9) an ay (ay + lEs + I) which is the result of Eqn. 2.16. (B.10) B . 2 D o w n h i l l f u n n e l i n g , E h r l i c h - S c h w o b e l b a r r i e r s a n d d e t a c h m e n t Eqn. B.2 still holds in zone 1 of Fig. 2.9 for p = pi, but in zone 2 we let F —> 0 in the B C F equation, with the result that, the adatom density in that region becomes P2(x) = C 3 X + C4 • (B- l l ) The boundary conditions (Eqns. 2.17) are repeated here x = 0 : -DVpi = K — —pi (B.12a) a l l x = l-Rinc: { v ; ; : ^ 2 (B.12b) x = l: -DVp=^-p-K^- (B.12c) l>ES <>ES which gives the following values for the c\ and C 4 : F (I - Rinc)(2lES + I + Rinc) M , O N C l - 2 ^ D (a,, + lES + Z) ( B ' 1 3 ) a\\K F {l-Rinc)(lES + l){l + ^ - R i n c ) C 4 = ^ + 2 a ^ ( a „ + ^ - M ) ~ ( R 1 4 ) and C2 and C3 can be expressed as C2 = a||Ci + ^ - - (B.15) C 3 = Cl - ^-{l - Rinc) (B.16) anL> Appendix B. Details of the Burton-Cabrera-Frank calculations 111 The adatom current is then FR-j = - £ ( V M 0 ) + V p 2 ( 0 ) + - ^ p (B.17) = -D(cl + c3) + ^ ^ (B.18) (B.19) F [l(a\\ - IES) +Rinc(^ES + Rinc)] ay (ay + IES + 1) which is Eqn. '2.18. B . 3 D o w n h i l l f u n n e l i n g , E h r l i c h - S c h w o b e l b a r r i e r s , d e t a c h m e n t a n d d e s o r p t i o n When desorption is introduced we have to apply Eq. 2.20. The solution for the adatom density p(x) is p(x) = C l exp (-^=) + C2 exp (—^=) + ^ (B.20) where ci = y(a\K-DTdesF){a\\L\ + l E S - y + L\y)/Y (B.21) c 2 = Ay{a\K - DrdesF){a\\ + AlES - y + Ay)/Y (B.22) and r = a\\D [ (A 2 + l)y(lES + ay) + ( A 2 - l)(y 2 + a^s)] (B.23) with y = \fD~Tfes and A as defined in Eq. 2.22. This leads to a very man-agable expression for the surface adatom current: ( A 2 - l ) ( / ^ - a y ) ( ^ - ^ f ) j(Z) = , ^4 r (B.24) which is Eq. 2.21. Appendix C. Finite difference scheme for nonlinear derivatives 112 Appendix C Finite difference scheme for nonlinear derivatives Figure C . l : (a) Application of the normal growth algorithm to a discrete representation of a ID surface. The surface is translated outwards by an amount ds = Adt, leading to a growth rate in the vertical direction of dh. (b) The 2D stencil [ABCDP] used to calculate the (V/i ) 2 term at point P. Conventional finite difference schemes for the nonlinear term (V/i) 2 in the K P Z equation (Eq. 3.11) are usually based on a centered difference approximation, such as the second order accurate V h = ^ ± ^ ± + 0{Ax2) (C.l) where hi is the height at the ith point on the ID surface, Ax is the spac-ing between the points and 0 is the usual "big-O" notation that signifies "approximation error on the order of {j" 1. As noted in Ref. [NB96], this : A version of this Chapter has been published. A. Ballestad, B. J. Ruck, J. H. Schmid, M . Adamcyk, E. Nodwell, C. Nicoll and T. Tiedje (2002) Surface morphology during GaAs Appendix C. Finite difference scheme for nonlinear derivatives 113 implementation fails to include the grid point hi, and can be highly unsta-ble. This restricts the range of parameters which can be used in simulations, more specifically it sets an upper limit to the ratio of the K P Z parameters X/v to a value that is too small for what we require in our morphology sim-ulations. Therefore, we have used an alternate implementation based on the normal growth interpretation of the K P Z term (Eq. 3.13), where the surface is translated outwards from the local surface normal by a constant amount. Consider the ID discrete representation of a surface shown in Figure C. la . The dashed lines show the positions of each surface element after translation outwards by a uniform amount. The thick solid line shows the new surface generated from the dashed lines, by choosing the maximum at any point where there is ambiguity in the choice of the new height. This procedure can be generalized to 2D. Figure C . lb shows the stencil used for the 2D calculation, where the relevant surface elements have been shaded. At most only one of the four shaded surface elements will actually decide the final height increment at point P. To determine which this is, we find the largest of HA, hs, he, and ho, and call this hi. Of the two remaining points of A , B, C, and D which are closest to this point, we find the one with the next largest h, and call this h2. Then, assuming Arc = Ay, The final step is to subtract Xdt from the matrix of height increments, thus leaving a matrix of values closely approximating A(V/i) 2 /2 . This algorithm is stable for any ratio X/v, as long as the time step is not excessive. Fur-thermore, the simulation code can be fully vectorized, leading to a 16-fold decrease in calculation time over an implementation using nested for loops to access the matrix elements. The algorithm has been fully tested on artificial surfaces for which (V/i) 2 can be calculated exactly. The conservative form of the K P Z equation (Eq. 3.15) is implemented via an extension to the nonconservative K P Z simulation. After calculating Ah on an N x N matrix for a given time step using Eq. C.2, the total volume represented by Ah is determined: dV — ^ Ax2Ah. This is done molecular beam epitaxy growth: Comparison of experimental data with simulations based on continuum growth equations, Phys. Rev. B 65: 205302 [BRS+02]. f Xdt, if [hi <hP kh2< hp] if [hi > hp & h2 < hp] (C.2) Appendix C. Finite difference scheme for nonlinear derivatives 114 prior to subtracting Xdt from Ah. The new matrix Ahc is then calculated: Ahc = Ah(dV/dV0) - Xdt, where dVb = X(NAx)2 is the volume that would be added to a flat surface during the same time step using the growth rule Eq. C.2. Note that this is not the same as ensuring a conservative growth term by simply subtracting the average surface height at each time step. The second order linear derivative term is implemented using a five point stencil on the 2D lattice, i.e. V h « vr Ax2 Ay2 (C.3) where Ax = Ay in our simulations. The fourth order linear term is imple-mented with a nine point stencil as V 4 / i hi+2,j — 4/ij+ij + 6/iij — 4/ij_ij + hi-2j ~ ~ ~ ~ ~ 4/it,j+i + 6hjj — Ahjj-i + hjj-2 ~~~ AT , + (C.4) The different coefficients (ux, vy, K X , K V ) allow for the inclusion of anisotropy in the simulations. Linear transformations of the terms in the numerical calculation allow the anisotropy axis to be rotated such that they match the elongation axis of the A F M images. It is obvious from Eq. 3.12 that this algorithm can also be used to simulate the fourth order nonlin-ear term, by simply applying the V 2 scheme to the matrix of (V/ i ) 2 values generated in the manner described above. Appendix D. Scaling of coupled growth equations 115 Appendix D Scaling of coupled growth equations The importance of scaling equations should be obvious to epitaxial growth modelers: we are trying to unite the physics of the atomic scale with the observable evolution of the large scale surfaces on the meso- to macroscale. Much of the work in this thesis is dedicated to exactly this problem: we utilize SOS models in order to simulate the underlying physics at the mi-croscopic or atomistic scale, while we use coarse grained partial differential equations in order to model what we observe with A F M images of surfaces grown (and etched) by M B E . In this Appendix, we outline the simple scal-ing behaviour of the coupled growth equations model that we developed in Chapter 5. D . l L a t e r a l s ca l i ng : x —>• bx We start by reminding the reader of the full form of the C G E model 1 t^n + V - j = F — ajj"19 t(aj1/i) (D.la) 3 t ( a J » = aDn2 + (J3Dn-a\\K)S (D.lb) j = -D{all(nVh + Vn) (D.lc) S = xjs2. + (aj /V/z) 2 (D.ld) While it is a well known fact that a large pitch grating will smooth at a slower rate than a small pitch grating, it is not equally clear whether Eqs. D . l will reproduce the complex shape evolutions that we witnessed in Ch. 5 if the lateral pitch of a grating is changed. This is the question that we wish to address in this Section: will changing the pitch of the surface pattern simply result in a different decay rate during growth, or will 1 The complete units are listed here: [n]=l/L, [j]=l/T, [F]=1/LT, [aj.]=[a||]=L, [h]=L, [a}=[0}=[(}=l, [D]=L2/T, .[K]=l/T, [5]=[50]=1/L, where L is length and T is time. Appendix D. Scaling of coupled growth equations 116 the shape evolution also change? We address this question by defining the desired variable transformation x -> x = bx (D.2) where b is the scaling factor with which we wish to scale the initial surface, and x, is the spatial variable in the new coordinate system. The remaining variables h, n and t all stay the same under this transformation. The variable transformation in Eq. D.2 requires the derivatives in the C G E to be changed according to: V V = 6"1 V, V 2 -> V 2 = b~2V2 (D.3) The candidates for a transformation are V —> V , F —> F, a —> a, D -» D, 0 p, K -> K, C -> C and S0 -> 50. We rewrite Eq. D . l a in the transformed variable set (denoted with the tilde): dtn = F - V -j-a^dtia^h) = F - V • (.-a±DCnVh - DVn) - a^dtia^h) = F + a±D(V • (nVh) + DV2n - a^d^a^h) b\n = F + a±D(b-2V • (nVh) + Db~2V2n - a^dtia^h) (D.4) The following transformation for our free parameters will do: F = F, D = b2D, ( = ( (D.5) The surface height itself (Eq. D.lb) undergoes the following transformation .dt{a1lh) = -aDn2 + ppriS - a[\KS dtia^h) = a(b2D)n2+~p{b2D)nS-a\\KS (D.6) We can immediately pick out the transformation a = b~2a (D.7) however, for the determination of P and So, we need to analyze the step density S2 = Srf + ia^Vh)2 = S02 + b-2{allVh)2 (D.8) Appendix D. Scaling of coupled growth equations 117 so the following transformation will do S0 = b^So, S = 6 _ 1 5 (D.9) which gives us (3 = b^p, K = bK (D.10) All our findings will be tabulated at the end of this chapter. D . 2 V e r t i c a l s c a l i n g : h —» ch Our second transformation deals with the increase in surface height of the surface structures: h^h = ch (D.l l ) where c is the vertical scaling factor. It is important to point out now that the lateral variable x remains the same, as does time t. However, from Eq. D . l a we see that whatever transformation is done to h must also be done to n, so vertical scaling implies that n^h = cn (D.12) All the temporal and spatial derivates remain the same under this trans-formation. Again, we rewrite the first part of the C G E in the transformed variable set: dth = F + bV2h + a-A}DC,V • (hVh) - a'1 dt{a-lh) cdtn = F + cbV2n + c2a^bQV-(nVh)-ca^ldt{a^h) (D.13) This time, the following transformations are necessary: F = cF, D — D, C = c - x C (D.H) The step density itself sees the following change: S2 = S02 + (a^Vh)2 = So2+ c2(allVh) S0 = cSo, S = cS (D.15) 2 Appendix D. Scaling of coupled growth equations 118 We then get dt{aji) = &Dn2 +J3DnS-a\\KS cdt(a±h) = aD{cn)2 + J3D{cn){cS) - a\{k(cS) (D.16) so that the complete transformation is: a = c~la, (3 = c-l(3, K = K (D.17) D . 3 T i m e s c a l i n g : t —• dt Our final transformation deals with stretching of time: t-^i=dt (D.18) where d is the scaling factor. In this transformation, the surface itself, and hence also the adatom concentration n, undergo no transformation. The spatial derivatives remain unchanged, but the time derivative change according to b\^di = dr1dt (D.19) We proceed as usual: dpi = F + DV2n + a J ^ C V • (riVh) - a^d^a^h) d^dtn = F + Z ) V 2 n + a J ^ C V • (nV/i) - cT1 d^dt^h) (D.20) from which we find that F = d~1F, D = d~1D, C = C (D.21) Since all the spatial derivatives are unchanged, the step-density remains unchanged, and this makes the second part of our analysis simple: d-^a-^h) = &Dh2 + pDhS-a\\kS d-ldt{a^-h) = a(dr1D)n2 + $(dr1D)nS-ai\kS (D.22) so we do not need to change a or J3, and the only remaining change required is: K = drxK (D.23) Appendix D. Scaling of coupled growth equations 119 D . 4 C o m p o s i t e s c a l i n g : ( x , h , t ) —>• (bx,ch,dt) The various transformations can be applied successively, so that in the same calculation, we can let all of x, h and t scale simultaneously. In that case, the total transformation on our system parameters are as follows: x h (n t This summary can be rewritten in the conventional scaling form (See for example Barabasi [BS95]), where x —» x' = bx, h —> h! = bah (where a is the roughness exponent that we discussed in Chapter 3) and t —> t' — bzt, simply by replacing c by ba and d by bz. This way, D scales as D —» D' = b2~zD, F F' = ba-zF, and so forth. ' D tfd^D K bd~lK F -> cd-xF c - c~\ a —> 6 2 c 1 a P , So 6 - ^ 5 0 D . 5 S p a t i a l f r e q u e n c y d e p e n d e n c e In the low amplitude limit, the morphology of a growing surface should be described well by the (deterministic) E W equation. By letting the surface undergo the composite scaling of section D.4, we obtain the following: dh ~dt = DV2h (D.24) dlTt = "H^h (D-25) so the required scaling of v is b2 v=-jV. (D.26) Appendix D. Scaling of coupled growth equations 120 We compare with our expression for v, derived in Ch. 5: v = 4(£ + K ) which is the same result that we got for E W above. In this linear approx-imation, the spatial frequency dependence goes as q2, as we have come to expect from the E W equation. (D.27) (D.28) (D.29) Bibliography 121 Bibliography [ABP+00] M . Adamcyk, A . Ballestad, T . Pinnington, T . Tiedje, iVI. Davies, and Y . Feng. Smoothing of textured GaAs surfaces during molecular beam epitaxy growth. J. Vac. Sci. Technol. B, 18:1488-1492, 2000. [Ada02] Martin B. Adamcyk. Epitaxial growth of dilute nitride-arsenide compound semiconductors by molecular beam epitaxy. PhD thesis, University of British Columbia, 2002. [AHD+00] G . Apostolopoulos, J . Herfort, L . Daweritz, K . H . Ploog, and M . Luysley. Reentrant mound formation in GaAs (001) ho-moepitaxy observed by ex situ atomic force microscopy. Phys. Rev. Lett, 84:3358-3361, 2000. [AVL62] G. M . Adelson-Velskii and Y . M . Landis. A n algorithm for the organization of information. Soviet Math. Dokl., 3:1259-1262, 1962. [Bal98] Anders Ballestad. Smoothing of patterned gallium arsenide surfaces during epitaxial growth. Master's thesis, University of British Columbia, 1998. [BC94] G. S. Bales and D. C. Chrzan. Dynamics of irreversible island growth during submonolayer epitaxy. Phys. Rev. B, 50:6057, 1994. [B.CF51] W. K. Burton, N. Cabrera, and F. Frank. The growth of crystals and the equilibrium structure of their surfaces. Phil. Trans. Roy. Soc, 243(866):299-358, 1951. [BDP98] W. Braun, L . Daweritz, and K . H . Ploog. Origin of electron diffraction oscillations during crystal growth. Phys. Rev. Lett., 80:4935-4938, 1998. Bibliography 122 [BE92] M . C. Bartelt and J. W. Evans. Scaling analysis of diffusion-mediated island growth in surface adsorption processes. Phys. Rev. B, 46:12675-87, 1992. [BGLK02] G. Biasiol, A. Gustafson, K . Leifer, and E . Kapon. Mecha-nisms of self-ordering in nonplanar epitaxy of semiconductor nanostructures. Phys. Rev. B, 65:205306-205320, 2002. [BJK+03] W. Braun, B. Jenichen, V . M . Kaganer, A . S. Shtukenberg, L . Daweritz, and K . H . Ploog. Island and pit kinetics on the growing GaAs (001) surface studied by synchrotron X-ray diffraction. J. Cryst. Growth, 251:56, 2003. [BKT+01] W. Braun, V . M . Kaganer, A . Trampert, H.-P. Schonherr, Q. Gong, R. Notzel, L . Daweritz, and K. H . Ploog. Diffu-sion and incorporation: shape evolution during overgrowth on structured substrates. J. Cryst. Growth, 227-228:51-55, 2001. [BLST05] A . Ballestad, Bayo Lau, J . H . Schmid, and T . Tiedje. Non-linear growth in GaAs molecular beam epitaxy, volume 859E. Mater. Res. Soc. Symp. Proc, 2005. [BMT90] V . Bortolani, N. H . March, and M . P. Tosi. Interaction of atoms and molecules with solid surfaces. Plenum Press, New York, 1990. [BRA+01] A . Ballestad, B. J . Ruck, M . Adamcyk, T . Pinnington, and T . Tiedje. Evidence from the surface morphology for nonlinear growth of epitaxial GaAs films. Phys. Rev. Lett., 86(11):2377-2380, 2001. [BRS+02] A. Ballestad, B. J . Ruck, J . H . Schmid, M . Adamcyk, E . Nod-well, C. Nicoll, and T . Tiedje. Surface morphology during GaAs molecular beam epitaxy growth: Comparison of ex-perimental data with simulations based on continuum growth equations. Phys. Rev. B, 65(20):205302, 2002. [BS95] A . - L . Barabasi and H. E . Stanley. Fractal concepts in sur-face growth. Cambridge University Press, Cambridge, United Kingdom, 1995. [BTS04a] A . Ballestad, T . Tiedje, and J . H . Schmid. Comment on "Transient evolution of surface roughness on patterned Bibliography 123 GaAs (001) during homo-epitaxial growth". Phys. Rev. Lett., 93(15):159601, 2004. [BTS+04b] A . Ballestad, T . Tiedje, J . H . Schmid, B. J . Ruck, and M . Adamcyk. Predicting GaAs surface shapes during M B E regrowth on patterned substrates. J. Cryst. Growth, 271:13-21, 2004. [BZ90] G . S. Bales and A. Zangwill. Morphological instability of a ter-race edge during step-flow growth. Phys. Rev. B, 41(9):5500-5508, 1990. [CC00] V . R. Coluci and M . A . Cotta. Influence of rough substrates on the morphology evolution of epitaxial films. Phys. Rev. B, 61:13703, 2000. [CCM+98] V . R. Coluci, M . A. Cotta, C. A . C. Mendonca, K . M . L -Landers, and M . M . G. de Carvalho. Surface morphologies in GaAs homoepitaxy: Mound formation and evolution. Phys. Rev. B, 58:1947, 1998. [CJZ77] E . L . Church, H . A . Jenkinson, and J . M . Zavada. Measure-ment of the finish of diamond-turned metal-surfaces by differ-ential light-scattering. Opt. Eng., 16(4):360-374, 1977. [CJZ79] E . L . Church, H . A . Jenkinson, and J . M . Zavada. Rela-tionship between surface scattering and microtopographic fea-tures. Opt Eng., 18(2):125-136, 1979. [CPP+89] P. I. Cohen, G . S. Petrich, P. R. Pukite, G . J . Whaley, and A. S. Arrott. Birth-death models of epitaxy : I. Diffraction oscillations from low index surfaces. Surf. Sci., 216:222-248, 1989. [DP97] Ed. P. M . Duxbury and T . J . Pence. Dynamics of Crystal Surfaces and Interfaces. Plenum Press, New York, USA, 1997. [EAC+00] J . Erlebacher, M . J . Aziz, E . Chason, M . B. Sinclair, and J. A. Floro. Nonclassical smoothening of nanoscale surface corrugations. Phys. Rev. Lett, 84(25):5800-5803, 2000. [EFFL94] H. -J . Ernst, F . Fabre, R. Folkerts, and J . Lapujoulade. Obser-vation of a growth instability during low temperature molec-ular beam epitaxy. Phys. Rev. Lett., 72:112, 1994. Bibliography 124 [EH66] G. Ehrlich and F. G. Hudda. Atomic view of surface diffusion: tungsten on tungsten. J. Chem. Phys., 44:1039, 1966. [EV94] I. Elkinani and J . Villain. Growth roughness and instability due to the Schwoebel effect: a one-dimensional model. J. de Physique I, 4:949, 1994. [Eva91] J . W. Evans. Factors mediating smoothness in epitaxial film growth. Phys. Rev. B, 43:3897, 1991. [EW82] S. F . Edwards and D. R. Wilkinson. The surface statistics of a granular aggregate. Proc. R. Soc. Lond. A, 381:17-31, 1982. [FA85] F . Family and J . Amar. The morphology and evolution of the surface in epitaxial and thin film growth: a continuum model with surface diffusion. World Scientific, Singapore, 1985. [FBSK04] S. Facsko, T . Bobek, A . Stahl, and H. Kurz. Dissipative con-tinuum model for self-organized pattern formation during ion-beam erosion. Phys. Rev. B, 69:153412, 2004. [FV85] F. Family and T . Viscek. Scaling of the active zone in the Eden process on percolation networks and the ballistic deposition model. J. Phys. A, 18:L75-L81, 1985. [FV91] F . Family and T . Viscek. Dynamics of Fractal Surfaces. World Scientific, Singapore, 1991. [GB91] L . Golubovic and R. Bruinsma. Surface diffusion and fluc-tuations of growing interfaces. Phys. Rev. Lett., 66:321-324, 1991. [GRG+01] F . G. Gibou, C. Ratsch, M . F . Gyure, S. Chen, and R . E . Caflish. Rate equations and capture numbers with implicit islands correlations. Phys. Rev. B, 63:115401-4, 2001. [HHZ95] T . Halpin-Healy and Y . - C . Zhang. Kinetic roughening phe-nomena, stochastic growth, directed polymers and all that. Aspects of multidisciplinary statistical mechanics. Phys. Rep., 254(4-6) :215-400, 1995. [Hin83] A. C. Hindmarsh. ODEPACK, A Systematized Collection of ODE Solvers, vol. 1 of IMACS Transactions on Scientific Com-puting. North-Holland, Amsterdam, 1983. Bibliography 125 [HIWK90] M . Hata, T . Isu, A; Watanabe, and Y . Katayama. Distribu-tions of growth rates on patterned surfaces measured by scan-ning microprobe reflection high-energy electron diffraction. J . Vac. Sci. Technol. B, 8(4):692-696, 1990. [HOW+94] A. W. Hunt, C. Orme, D. R. M . Williams, B. G . Orr, and L. M . Sander. Instabilities in M B E growth. Europhys. Lett., 27(8):611-616, 1994. [HSJN91] L . Hansen, P. Stoltze, K . W . Jacobsen, and J . K. N0rskov. Self-diffusion on copper surfaces. Phys. Rev. B, 44:6523, 1991. [IO00] M . Itoh and T . Ohno. Absence of a step-edge barrier on a polar semiconductor surface with reconstruction. Phys. Rev. B, 62:1889, 2000. [ItoOl] M . Itoh. Kinetic Monte Carlo study of step asymmetry and stable step orientations on GaAs (001). Phys. Rev. B, 64:045301, 2001. [JLB+96] M . D. Johnson, K. T . Leung, A. Birch, B. G . Orr, and J . Ter-soff. Adatom concentration on GaAs (001) during M B E an-nealing. Surf. Sci., 350:254-258, 1996. [JLT93] S. R. Johnson, C. Lavoie, and T . Tiedje. Semiconductor sub-strate temperature measurement by diffuse reflectance spec-troscopy in molecular beam epitaxy. J. Vac. Sci. Techol, 11 (3): 1007-1010, 1993. [JOH+94] M . D. Johnson, C. Orme, A . W. Hunt, D. Graff, J . Sudijono, L. M . Sander, and B. G . Orr. Stable and unstable growth in molecular beam epitaxy.. Phys. Rev. Lett., 72:116-119, 1994. [Joh95] S. R. Johnson. Optical bandgap thermometry in molecular beam epitaxy. PhD thesis, University of British Columbia, 1995. [Joy03] B. A . Joyce. In situ studies of the M B E growth of III-V sys-tems using R H E E D and S T M . J. Mat. Sci. - Mat. in Elec-tronics, 14:591-598, 2003. [JT97] S. R. Johnson and T . Tiedje. Effect of substrate thick-ness, back surface texture, reflectivity, and thin film inter-Bibliography 126 ference on optical band-gap thermometry. J. Cryst. Growth, 175/176:273-280, 1997. [KEOO] K . Kyuno and G. Ehrlich. Cluster diffusion and dissociation in the kinetics of layer growth: A n atomic view. Phys. Rev. Lett, 84:2658-2661, 2000. [KK89] J . M . Kim and J . M . Kosterlitz. Growth in a restricted solid-on-solid model. Phys. Rev. Lett, 62:2289-2292, 1989. [KK04] J . Kallunki and J . Krug. Breakdown of step flow growth in unstable homoepitaxy. Europhys. Lett, 66:749-755, 2004. [KKK02] J . Kallunki, J . Krug, and M . Kotrla. Competing mecha-nisms for step meandering in unstable growth. Phys. Rev. B, 65:205411, 2002. [KKvOO] M . Kotrla, J . Krug, and P. Smilauer. Submonolayer epitaxy with impurities: Kinetic monte carlo simulations and rate-equation analysis. Phys. Rev. B, 62:2889, 2000. [KMS98] P. Kratzer, C. G. Morgan, and M . SchefHer. Density-functional theory studies on microscopic processes of GaAs growth. Prog. Surf. Sci., 59(1-4): 135-147, 1998. [KNT+94] S. Koshiba, Y . Nakamura, M . Tsuchiya, H . Noge, H . Kano, Y . Nagamune, T . Noda, and H. Sakaki. Surface diffusion pro-cesses in molecular beam epitaxial growth of GaAs and AlAs as studied on GaAs (001)-(111)B facet structures. J. Appl. Phys., 76:4138, 1994. [KPS93] j . Krug, M . Plischke, and M . Siegert. Surface diffusion cur-rents and the universality classes of growth. Phys. Rev. Lett., 70(21) :32.71-3274, 1993. [KPS02] P. Kratzer, E . Penev, and M . SchefHer. First-principles studies of kinetics in epitaxial growth of III-V semiconductors. Appl. Phys. i4,'75(l):79-88, 2002. [KPVC00] R. Kunkel, B. Poelsema, L . K. Verheij, and G. Comsa. Reen-trant layer-by-layer growth during molecular-beam epitaxy of metal-on-metal substrates. Phys. Rev. Lett., 65:733-736, 2000. Bibliography 127 [KPZ86] M . Kardar, G. Parisi, and Y . - C . Zhang. Dynamic scaling of growing interfaces. Phys. Rev. Lett, 56(9):889-892, 1986. [Kru97] J . Krug. Origins of scale invariance in growth processes. Adv. Phys., 46(2):139-282, 1997. [KSTEP04a] H. -C . Kan, S. Shah, T . Tadyyon-Eslami, and R. J . Phaneuf. Reply to Phys. Rev. Lett. 93 (2004) 159601. Phys. Rev, Lett, 93:159602,2004. [KSTEP04b] H . -C . Kan, S. Shah, T . Tadyyon-Eslami, and R. J . Pha-neuf. Transient evolution of surface roughness on patterned GaAs (001) during homoepitaxial growth. Phys. Rev. Lett., 92:146101-4, 2004. [Kur94] Y . Kuramoto. Chemical Oscillations, Waves, and Turbulence. Springer-Verlag, Berlin, 1994. [LG00] S. B. Lee and B. C . Gupta. Influence of small-cluster mobility on the island formation in molecular beam epitaxy. Phys. Rev. B, 62:7545-7552, 2000. [LL67] L . Landau and E . Lifshitz. Physique statistique. Edition Mir, Moscow, 1967. [LS91] Z.-W. Lai and S. Das Sarma. Kinetic growth with surface relaxation: Continuum versus atomistic models. Phys. Rev. Lett, 66:2348-2351, 1991. [LZ02] M . G . Lagally and Z. Zhang. Thin-film cliffhanger. Nature, 417:907-910, 2002. [Mar94] Ivan Markov. Kinetics of surfactant-mediated epitaxial growth. Phys. Rev. B. Rapid Comm., 50(15):11271-11274, 1994. [Mar96] Ivan Markov. Method for evaluation of the Ehrlich-Schwoebel barrier to interlayer transport in metal homoepitaxy. Phys. Rev. B., 54(24):17930-17937, 1996. [MG88] A . Madhukar and S. V . Ghaisas. The nature of molecular beam epitaxial growth examined via computer simulations. CRC Crit. Rev. Sol. State and Mater. Sci., 14:1, 1988. Bibliography 128 [MKC+02] T . Michely, M . Kalff, G . Comsa, M . Strobel, and K . - H . Heinig. Coarsening mechanisms in surface morphological evolution. J. Phys.: Condens. Matter, 14:4177-4185, 2002. [MRRT53] N. Metropolis, A . W. Rosenbluth, M . N. Rosenbluth, and A. H . Teller. Equation of state calculations by fast computing machines. J. Chem. Phys., 21:1087, 1953. [Mul57] W. W. Mullins. Theory of thermal grooving. J. Appl. Phys., 28(3):333-339, 1957. [Mul59] W. W. Mullins. Flattening of a nearly planar solid surface due to capillarity. J. Appl. Phys., 30:77, 1959. [NB96] T . J . Newman and A. J . Bray. Strong-coupling behavior in dis-crete Kardar-Parisi-Zhang equations. J. Phys. A, 24(29):7917, 1996. [NB99] M . E . J. Newman and G . T . Barkema. Monte Carlo Methods in Statistical Physics. Oxford University Press, Oxford, United Kingdom, 1999. [NCC+96] J . E . Van Nostrand, S. J . Chey, D. G. Cahill, A . E . Botchkarev, and H . Morkog. Surface morphology of GaAs (001) grown by solid- and gas-source molecular beam epitaxy. Surf. Sci., 346:136-144, 1996. [NCH+95] J . E . Van Nostrand, S. J . Chey, M . A. Hasan, D. G . Cahill, and J. E . Greene. Surface morphology during multilayer epitaxial growth of Ge (001). Phys. Rev. Lett, 74:1127-1130, 1995. [NDJZ85] J . H . Neave, P. J. Dobson, B. A. Joyce, and J . Zhang. Reflec-tion high-energy electron diffraction oscillations from vicinal surfaces: a new approach to surface diffusion measurements. Appl. Phys. Lett, 47(2):100-102, 1985. [NK96] T . J . Newman and H. Kallabis. Strong-coupling behavior for the Kardar-Parisi-Zhang equation. J. Phys. I France, 6:373-383, 1996. [OJL+95] C. Orme, M . D. Johnson, K . T . Leung, B. G . Orr, P. Smilauer, and D. D. Vvedensky. Studies of large scale unstable growth formed during GaAs (001) homoepitaxy. J. Cryst. Growth, 150:128-135, 1995. Bibliography 129 [OJS+94] C . Orme, M . D. Johnson, J . L . Sudijono, K . T . Leung, and B. G . Orr. Large-scale surface-structure formed during GaAs (001) homoepitaxy. Appl. Phys. Lett, 64(7):860, 1994. [Pin99] Thomas H. Pinnington. Surface morphology dynamics in strained-layer epitaxy. PhD thesis, University of British Columbia, 1999. [PT00] P. Politi and A. Torcini. Coarsening in surface growth models without slope selection. J. Phys. A: Math. Gen., 33:L77-L82, 2000. [PV96] P. Politi and J . Villain. Ehrlich-Schwoebel instability in Molecular-Beam-Epitaxy: a minimal model. Phys. Rev. B, 54:5114, 1996. [PV98] A. Pimpinelli and J . Villain. Physics of crystal growth. Cam-bridge University Press, Cambridge, United Kingdom, 1998. [RH93] K . Radhakrishnan and A. C. Hindmarsh. Description and use of L S O D E , the Livermore Solver for Ordinary Differential Equations. Lawrence Livermore National Laboratories report, (UCRL-ID-113855), 1993. [RSLP91] Z. Racz, M . Siegert, D. Liu, and M . Plischke. Scaling prop-erties of driven interfaces: Symmetries, conservation laws and the role of constraints. Phys. Rev. A, 43(10):5275-5283, 1991. [RvK96] M . Rost, P. Smilauer, and J . Krug. Unstable epitaxy on vicinal surfaces. Surf. Sci., 369:393-402, 1996. [RZv94] C. Ratsch, A . Zangwill, and P. Smilauer. Saturation and scal-ing of epitaxial island densities. Phys. Rev. Lett., 72:3194, 1994. [Sar93] S. Das Sarma. Kinetic surface roughening and molecular beam epitaxy. Fractals, 1:784, 1993. [SBR+02] J . H . Schmid, A . Ballestad, B. J . Ruck, M . Adamcyk, and T . Tiedje. Kinetic roughening of GaAs(OOl) during thermal C l 2 etching. Phys. Rev. B, 65:155315-9, 2002. [Sch] Jens H . Schmid. Unpublished. Bibliography- 130 [Sch68] L . I. Schiff. Quantum Mechanics, 3rd ed. McGraw-Hill, New York, USA, 1968. [Sch69] R. L . Schwoebel. Step motion on crystal surface II. J. Appl. Phys., 40:614-618, 1969. [Sch04] Jens H . Schmid. Evolution of surface texture in thermal chlo-rine etching and molecular, beam epitaxy of gallium arsenide. PhD thesis, University of British Columbia, 2004. [SGG89] T . Sun, H . Guo, and M . Grant. Dynamics of driven interfaces with a conservation law. Phys. Rev. A, 40:6763, 1989. [SGL+03] S. Shah, T. J . Garrett, K . Limpaphayom, T . Tadyyon-Eslami, H . -C . Kan, and R. J . Phaneuf. Pattern-based investigation of the length-scale dependence of the surface evolution during multilayer epitaxial growth. Appl. Phys. Lett., 83:4330-4332, 2003. [Siv77] G. I. Sivashinsky. Non-linear analysis of hydrodynamic insta-bility in laminar flames -I. derivation of basic equations. Acta Astronautica, 4:1177-1206, 1977. [SJE+93] J . Sudijono, M . D. Johnson, M . B. Elowitz, C. W. Snyder, and B. G. Orr. A n S T M study of molecular-beam epitaxy growth of GaAs. Surf Sci., 280:247, 1993. [SJS+92] J. Sudijono, M . D. Johnson, C. W. Snyder, M . B. Elowitz, and B. G . Orr. Surface evolution during molecular-beam epitaxy deposition of GaAs. Phys. Rev. Lett, 69:2811-2814, 1992. [SK94] S. Das Sarma and R. Kotlyar. Dynamical renormalization group analysis of fourth-order conserved growth nonlinearities. Phys. Rev. E, 50:R4275-R4278, 1994. [SKJ+92] K . Sneppen, J . Krug, M . H . Jensen, C. Jayaprakash, and T . Bohr. Dynamic scaling and crossover analysis for the Kuramoto-Sivashinsky equation. Phys. Rev. A, 46:7351, 1992. [SKR00] S. Schinzer, S. Koehler, and G. Reents. Ehrlich-Schwoebel barrier controlled slope selection in epitaxial growth. Eur. Phys. J. B, 15:161-168, 2000. Bibliography 131 [SM80] G . I. Sivashinsky and D. M . Michelson. On irregular wavy flow of a liquid film down a vertical plane. Prog. Theor. Phys., 63:2112, 1980. [SPSZ95] J . A . Stroscio, D. T . Pierce, M . D. Stiles, and A. Zangwill. Coarsening of unstable surface features during Fe (001) ho-moepitaxy. Phys. Rev. Lett, 75:4246-4249, 1995. [SR97] L . F . Shampine and M . W. Reichelt. The M A T L A B O D E Suite. SIAM J. of Scient Comp., 18:1-22, 1997. [SRF03] V . B. Shenoy, A . Ramasubramaniam, and L . B. Freund. A variational approach to nonlinear dynamics of nanoscale sur-face modulations. Surf. Sci., 529:365-383, 2003. [SRK99] L . F . Shampine, M . W. Reichelt, and J . A . Kierzenka. Solv-ing index-1 DAEs in M A T L A B and Simulink. SIAM Review, 41:538-552, 1999. [SS66] R. L . Schwoebel and E . J. Shipsey. Step motion on crystal surface. J. Appl. Phys., 37:3682-3686, 1966. [SS94] R. Stumpf and M . Scheffler. Theory of self-diffusion and growth of A l ( l l l ) . Phys. Rev. Lett, 72:254, 1994. [STMB04] J . H . Schmid, T . Tiedje, R. Mar, and A. Ballestad. Surface pattern transfer in GaAs with molecular beams of CI2. Phys. Rev. B, 70:045315, 2004. [Sto94] P. Stoltze. Simulation of surface defects. J. Phys.: Condens. Matter, 6:9495, 1994. [SVW+92] T . Shitara, D. D. Vvedensky, M . R. Wilby, J . Zhang, J . H . Neave, and B. A . Joyce. Misorientation dependence of epi-taxial growth on vicinal GaAs (001). Phys. Rev: B, 46:6825, 1992. [TKWR95] K. Thiirmer, R. Koch, M . Weber, and K . H . Rieder. Dynamic evolution of pyramid structures during growth of epitaxial Fe (001) films. Phys. Rev. Lett, 75:1767-1770, 1995. [TWUC96] F. Tsui, J . Wellman, C. Uher, and R. Clarke. Morphology transition and layer-by-layer growth of Rh (111). Phys. Rev. Lett, 76:3164-3167, 1996. Bibliography 132 [TWY+03] Yasushi Takahashi, Shinichi Watanabe, Masahiro Yoshita, Hi-rotake Itoh, Yuhei Hayamizu, Hidefumi Akiyama, Loren N. Pfeiffer, and Ken W. West. Imaging of emission patterns in a t-shaped quantum wire laser. Appl. Phys. Lett., 83:4089, 2003. [Var02] Raffaele Vardavas. Fluctuations and scaling in ID irreversible film growth models. PhD thesis, Imperial College, University of London, 2002. [VC90] D. D. Vvedensky and S. Clarke. Recovery kinetics during interrupted epitaxial growth. Surf. Sci., 225:373, 1990. [Ven87] J. A . Venables. Nucleation calculations in a pair-binding model. Phys. Rev. B, 36:4153, 1987. [Vil91] J . Villain. Continuum models of crystal growth from atomic beams with and without desorption. J. de Physique I, 1:19-42, 1991. [vS16] M . von Smoluchowski. Drei Vortrage iiber Diffusion, Brownsche Molekularbewegung und Koagulation von Kol-loidteilchen. Z. Phys. Chem., 17:557, 1916. [vS17] M . von Smoluchowski. Versuch einer mathematischen Theorie der Koagulationskinetik kolloider Losungen. Z. Phys. Chem., 92:129,1917. [VSH84] J . A . Venables, G. D. T . Spiller, and M . Hanbiicken. Nucle-ation and growth of thin films. Reports on Progress in Physics, 47:399-459, 1984. [vV93] P. Smilauer and D. D. Vvedensky. Step-edge barriers on GaAs (001). Phys. Rev. B, 48:17603, 1993. [vV95] P. Smilauer and D. D. Vvedensky. Coarsening and slope evo-lution during epitaxial growth. Phys. Rev. B, 52(19):14263-14272, 1995. [VveOO] D. D. Vvedensky. Scaling functions for island-size distribu-tions. Phys. Rev. B, 62:15435-38, 2000. [WE91] S. C. Wang and G. Ehrlich. Atom incorporation at surface 'clusters: an atomic view. Phys. Rev. Lett., 67(18):2509-2512, 1991. Bibliography 133 [WE93] S. C. Wang and G. Ehrlich. Atom condensation at lattice steps and clusters. Phys. Rev. Lett, 71:4174-4177, 1993. [WGJ76] J . D. Weeks, G. H . Gilmer, and K. A . Jackson. Analytical theory of crystal-growth. J. Chem. Phys., 65(2):712-720, 1976. [Wil89] M . R. Wilby. Kinetics of Systems far from Equlibrium. PhD thesis, Imperial College, University of London, 1989. [WV90] D. E . Wolf and J . Villain. Growth with surface-diffusion. Eu-rophys. Lett, 13:389-394, 1990. [XLT02] M . H . Xie, S. Y . Leung, and S. Y . Tong. What causes step bunching - negative Ehrlich-Schwoebel barrier versus positive incorporation barrier. Sur. Sci. Lett, 515:L459-L463, 2002. [Yak81] V . Yakhot. Large-scale properties of unstable systems gov-erned by the Kuramoto-Sivashinski equation. Phys. Rev. A, 24:642, 1981. [YLB+99] H . Yang, V . P. LaBella, D; W. Bullock, Z, Ding, J . B. Smath-ers, and P. M . Thibado. Activation energy for Ga diffusion on the GaAs (001)-(2 x 4) surface: an M B E - S T M study. J. Cryst. Growth, 201/202:88, 1999. [Zan88] Andrew Zangwill. Physics at Surfaces. Cambridge University Press, Cambridge, United Kingdom, 1988. [ZL99] Z. Zhang and M . G. Lagally. Morphological Organization in Epitaxial Growth and Removal. World Scientific Publishing, New York, 1999. [ZW97] J . -K. Zuo and J . F . Wendelken. Evolution of mound morphol-ogy in reversible homoepitaxy on Cu (100). Phys. Rev. Lett., 78:2791-2794, 1997. 


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items