ABSTRACT Title of dissertation: A MACROSCALE PERSPECTIVE OF NEAR-EQUILIBRIUM RELAXATION OF STEPPED CRYSTAL SURFACES John Quah, Ph.D. 2009 Dissertation directed by: Professor Dionisios Margetis Department of Mathematics Institute for Physical Science and Technology Crystal surfaces serve a crucial function as building blocks in small electronic devices, especially for mobile communications technology and photovoltaics. In the history of computing, for example, a crucial innovation that hastened the demise of vacuum tube computers was the etching of patterns on surfaces of semiconduc- tor materials, which led to the integrated circuit. These early procedures typically worked with the materials at very high temperatures, where the thermally rough surface could be modeled from the perspective of continuum thermodynamics. More recently, with the drive towards smaller devices and the accompanying reduction of the lifetime of surface features, manufacturing conditions for the shaping of crys- tal surfaces have shifted to lower temperatures. At these lower temperatures the surface is no longer rough. In order to describe an evolving surface under typical experimental conditions today, we need to consider the processes that take place at the nanoscale. Nanoscale descriptions of surface evolution start with the motion of adsorbed atoms (adatoms). Because of their large numbers, the concentration of adatoms is a meaningful object to study. Restricted to certain bounded regions of the surface, the adatom concentration satisfies a diffusion equation. At the boundaries between these regions, the hopping of adatoms is governed by kinetic laws. Real-time obser- vation of these nanoscale processes is difficult to achieve, and experimentalists have had to devise creative methods for inferring the relevant energy barriers and kinetic rates. In contrast, the real-time observation of macroscale surface evolution can be achieved with simpler imaging techniques. Motivated by the possibility of experi- mental validation, we derive an equation for the macroscale surface height, which is consistent with the motion of adatoms. We hope to inspire future comparison with experiments by reporting the novel results of simulating the evolution of the macroscale surface height. Many competing models have been proposed for the diffusion and kinetics of adatoms. Due to the difficulty of observing adatom motion at the nanoscale, few of the competing models can be dismissed outright for failure to capture the observed behavior. This dissertation takes a few of the nanoscale models and sys- tematically derives the corresponding macroscopic evolution laws, of which some are implemented numerically to provide data sets for connection with experiments. For the modeling component of this thesis, I study the effect of anisotropic adatom diffusion at the nanoscale, the inclusion of an applied electric field, the desorption of adatoms, and the extension of linear kinetics in the presence of step permeability. Analytical conjectures based on the macroscale evolution equation are presented. For the numerical component of this thesis, I select a few representative simulations using the finite element method to illustrate the most salient features of the surface evolution. A MACROSCALE PERSPECTIVE OF NEAR-EQUILIBRIUM RELAXATION OF STEPPED CRYSTAL SURFACES by John Quah Dissertation submitted to the Faculty of the Graduate School of the University of Maryland, College Park in partial fulfillment of the requirements for the degree of Doctor of Philosophy 2009 Advisory Committee: Professor Dionisios Margetis, Chair/Advisor Professor Ricardo Nochetto Professor Manoussos Grillakis Professor Theodore Einstein Professor Janice Reutt-Robey Dedication To M. Z., who taught me to leave bridges intact and encouraged me to repair the ones I burned in moments of impulsiveness. ii Acknowledgments I owe my gratitude to all the people who have made this thesis possible and without whom my perseverence in the graduate program might not have endured six turbulent years. First and foremost I?d like to thank my advisor, Professor Dionisios Margetis (Dio), who presented me with challenging research problems and equipped me with the specialized techniques needed to solve them. If I had not been daring enough to enroll in his course on advanced analytical methods the first time he taught it at Maryland (Spring 2007), this collaboration and the resulting thesis might never have taken off. Apart from being an invaluable source of mathematical knowledge, Dio also made a special effort to involve me in the conferences and workshops he organized. These meetings brought me in contact with a larger network of research scientists and mathematicians, whose examples inspired me to pursue my own the- sis problem with more alacrity. The professional training I received under Dio?s supervision greatly exceeded any expectations I had upon entering the graduate program. Other faculty at the University of Maryland deserve credit for their role in the development of this thesis. Thanks to Professor Georg Dolzmann, whose course in numerical analysis imparted an appreciation for the beauty of the field. Thanks to Professors Eugenia Kalnay, Da-Lin Zhang, and Jim Carton for making accessible to applied mathematics students the world of atmospheric and oceanic science, along with the quantitative tools needed to understand the interactions among processes iii at different time and length scales. Professors Stuart Antman, Manoussos Grillakis, and Matei Machedon conveyed their knowledge of partial differential equations and nonlinear analysis with exceptional clarity. I am particularly indebted to several faculty in physics, chemistry, and ma- terials science. Professor Ted Einstein of the Physics Department co-organized a research interaction team (RIT) with Dio and Professor Tzavaras in the spring of 2008. This series featured many informative lectures and also gave me the opportu- nity to practice my own blackboard presentation skills. Professor Einstein contin- ued to provide guidance and physical intuition as my thesis research progressed into new areas. Similarly, conversations with Professor John Weeks of the Chemistry Department and IPST1, and with Professor Ray Phaneuf of the Materials Science and Engineering Department, helped point out some of the problems I would face in connecting my theoretical investigations with real materials. Professor Janice Reutt- Robey of the Chemistry Department has been gracious with her time by offering to serve as the Dean?s Representative on my thesis defense committee. I also appreciate the financial support of the Monroe Martin Graduate Fel- lowship in summer 2008, which would not have come to my attention were it not for Professor Jim Yorke of the Mathematics Department and IPST. Professors Ein- stein and Margetis wrote the letters to nominate me for this stipend. During that summer, collaboration with Dio and the REU student Jerrod Young led to my first publication in Physical Review E. My own undergraduate research experience, working out a problem in the 1Institute of Physical Science and Technology iv mathematicsofcompositematerialsundertheguidanceofProfessorYuryGrabovsky of Temple University, was decisive in encouraging me to pursue a graduate degree in applied mathematics. Little did I know at the time that a colleague of one of Professor Grabovsky?s collaborators would eventually supervise my thesis, after first introducing me to interesting mathematical problems motivated by a different branch of materials science. Other faculty who deserve mention for their helpful dis- cussions are Bob Kohn of New York University, Joachim Krug of the Institute for Theoretical Physics (Universit?at K?oln), John Venables of Arizona State University, Jonah Erlebacher of Johns Hopkins University, Jim Evans of Iowa State University, and Richard James of the University of Minnesota. I am grateful to Dr. Alex Prasertchoung of MRSEC2 for pushing me to be a mentor for the AIP3 Student Science Conference (a truly rewarding experience, not least because I got to see how much talent and fluid intelligence today?s middle school students possess), and for pushing me around on the tennis court. I also had the good fortune of receiving financial support from MRSEC through grant NSF- MRSEC DMR0520471 during the spring of 2009, for which I thank Professor Ellen Williams, director of the UMD MRSEC. By sheer coincidence, I happened to be finalizing my thesis precisely when the NSF4 was scheduled to conduct a site visit at the UMD MRSEC. Professor Williams invited me to present a poster at the site visit, which gave me the rare opportunity to have my work scrutinized in person by scientists from outside Maryland. 2Materials Research Science and Engineering Center 3American Institute of Physics 4National Science Foundation v Among all the faculty behind this thesis, Professor Ricardo Nochetto deserves recognition for sharing his extensive knowledge of the finite element method, both in private discussions and in a formal classroom context. For his patience with me in the latter setting, in a semester which challenged my management of time among such disparate tasks as thesis writing, teaching and grading, work at the food co-op, organizing numerical simulations, and finishing up a paper for publication, Professor Nochetto can hardly be thanked enough. Collaboration with him and his former postdoc, Andrea Bonito, always went smoothly, despite the differing perspectives we brought as mathematical modelers and numerical analysts. The staff of the math department also offered just the right degree of support, without drowning me in paperwork or other bureaucratic hurdles. Alverda McCoy of the AMSC5 program has been ever helpful with registration issues and completion of forms along the road to graduation. Mark Tilmes graciously applied his computer expertise to maintain all the software I needed for the numerical portion of this thesis. Bill Schildknecht accommodated several last-minute requests to change my teaching assignment, no small feat considering how many courses and other graduate teaching assistants were potentially affected. For their support in my personal life, the following people deserve a word of thanks. David Shoup, Carter Price, J. C. Flake, Charles Glover, Chris Manon, Jeff Heath, Avanti Athreya, Ryan Janicki, Toni Watson, I-Kun Chen, and other occu- pants of ?the Zoo? in 2003 helped with my acclimation to graduate life. Chris Zorn, David Bourne, Seba Pauletti, J. T. Halbert, and other more advanced graduate stu- 5Applied Mathematics, Statistics, and Scientific Computation vi dents showed by example how to balance work and play while still making adequate progress toward a finished thesis. Dr. Jonathan Kandell of the Counseling Cen- ter offered some much-needed relief when relations among peers became strained. Under the care of Dr. Stephen Feinberg, my astigmatic eyes received the proper corrective lenses and vision therapy to enable two years of intense near-focus work. Dr. Robert Delong and the neurosurgeons at Duke University Medical Center were instrumental in curing me of epilepsy. Dr. James Wainer of Raleigh, NC, and Dr. Eric Griffey of Lewiston, ME, kept me on psychotropic drugs long enough to delay until graduate school the broadening of social and intellectual horizons that I should have experienced in my undergraduate years. I also thank the badminton club and its dedicated faculty sponsor, Professor Peter Teuben, for providing access to a fun-filled team sport and a supportive en- vironment on the many occassions when I needed to release frustrations. On my very first visit to the badminton club, I made contact with one person in particular, whose friendship led to many learning experiences in disciplines ranging from the humanities to the fine arts. This friendship finally culminated in the healing of a rift between me and the Maryland Food Co-op, where I have worked or volunteered since summer 2007. The community of workers and volunteers at the food co-op has kept me anchored and well-nourished throughout the hectic two years of working on my thesis. Finally, I?d like to thank my family for providing a loving and nurturing envi- ronment throughout my formative years. Their continued support during graduate school has not gone unnoticed, although I have not always made the time to recip- vii rocate and keep in touch with them. When I emerge from the tunnel vision that will allow me to see the thesis through to the end, my parents deserve to be the first to hear that I?m still alive. viii Table of Contents List of Abbreviations xiii 1 Introduction 1 1.1 Chapter overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Physical description of surfaces at various length scales . . . . . . . . 1 1.3 Mathematical models of various length scales . . . . . . . . . . . . . . 2 1.4 Connecting step models to experiments . . . . . . . . . . . . . . . . . 10 1.5 Macroscale limit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.6 Overview of this thesis . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2 Macroscopic equations and numerical solutions 16 2.1 Macroscopic theories of crystal surfaces . . . . . . . . . . . . . . . . . 16 2.1.1 Effect of different mass transport mechanisms . . . . . . . . . 17 2.1.2 Macroscopic theory above the roughening temperature . . . . 18 2.1.3 Macroscopic theory below the roughening temperature . . . . 21 2.2 Numerical simulations . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3 One-dimensional nanoscale models of interacting steps 27 3.1 Geometry of straight steps . . . . . . . . . . . . . . . . . . . . . . . . 27 3.2 BCF model for interacting straight steps . . . . . . . . . . . . . . . . 28 3.3 1D interaction energy and the step chemical potential . . . . . . . . . 30 3.4 BCF model for interacting concentric circular steps . . . . . . . . . . 37 3.5 Interaction energy and step chemical potential for circular steps . . . 40 4 Elastic-dipole interactions between steps in 2 dimensions 43 4.1 From electrostatics to linear elasticity . . . . . . . . . . . . . . . . . . 43 4.1.1 Field and potential of an electric dipole . . . . . . . . . . . . . 44 4.1.2 Energy of a charge distribution in an electric field . . . . . . . 45 4.1.3 Interaction energy of two force dipoles . . . . . . . . . . . . . 47 4.2 Integral for the energy of interacting steps . . . . . . . . . . . . . . . 49 4.2.1 Mellin transform of the axisymmetric integral . . . . . . . . . 52 4.2.2 Asymptotic expansion using the inverse Mellin transform . . . 54 4.2.3 Corrections to the leading-order term under axisymmetry . . . 55 4.3 Elastic interaction energy between steps without axisymmetry: A brief discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 5 Review of equations for interacting steps in 2+1 dimensions, and macroscale limit 65 5.1 Geometry of a step train in 2+1 dimensions . . . . . . . . . . . . . . 66 5.2 BCF model with step interactions in 2+1 dimensions . . . . . . . . . 68 5.3 Approximations for slowly varying step train . . . . . . . . . . . . . . 73 5.4 Macroscale limit with isotropic diffusion in 2+1 dimensions . . . . . . 75 5.4.1 Adatom flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 ix 5.4.2 Macroscale step chemical potential . . . . . . . . . . . . . . . 78 5.4.3 Mass conservation for adatoms . . . . . . . . . . . . . . . . . . 81 5.4.4 Evolution equation for surface height . . . . . . . . . . . . . . 82 5.5 Scaling ansatz for the height profile . . . . . . . . . . . . . . . . . . . 83 5.6 Presence of facets in the variational formulation . . . . . . . . . . . . 84 6 Macroscale equation with terrace anisotropy and step edge diffusion 86 6.1 Approximations for fast and slow step variables . . . . . . . . . . . . 87 6.2 Macroscale adatom flux . . . . . . . . . . . . . . . . . . . . . . . . . 90 6.3 Alternative approach to macroscale limit: Taylor expansions . . . . . 92 6.4 Mass conservation law and total surface flux . . . . . . . . . . . . . . 93 6.5 Evolution equation for the surface height . . . . . . . . . . . . . . . . 95 6.6 New scaling laws with terrace anisotropy and edge diffusion . . . . . 96 6.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 7 Step permeability and macroscale limit 103 7.1 Step permeability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 7.2 Macroscale limit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 7.2.1 Straight steps . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 7.2.2 Circular steps with axisymmetry . . . . . . . . . . . . . . . . 110 7.2.3 Steps in 2+1 dimensions . . . . . . . . . . . . . . . . . . . . . 114 7.3 Discussion and interpretation . . . . . . . . . . . . . . . . . . . . . . 118 8 Electromigration and desorption 123 8.1 Macroscale limit and assumptions . . . . . . . . . . . . . . . . . . . . 125 8.2 One-dimensional step models . . . . . . . . . . . . . . . . . . . . . . 127 8.2.1 Straight steps . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 8.2.1.1 Electric field with no desorption . . . . . . . . . . . . 129 8.2.1.2 Desorption with no electric field . . . . . . . . . . . . 130 8.2.2 Concentric circular steps . . . . . . . . . . . . . . . . . . . . . 130 8.2.2.1 Electric field with no desorption . . . . . . . . . . . . 132 8.2.2.2 Desorption with no electric field . . . . . . . . . . . . 134 8.3 2-dimensional step geometry . . . . . . . . . . . . . . . . . . . . . . . 135 8.3.1 Step geometry and formulation . . . . . . . . . . . . . . . . . 136 8.3.2 Equations of motion . . . . . . . . . . . . . . . . . . . . . . . 136 8.3.3 Approximate solution of diffusion equation . . . . . . . . . . . 137 8.4 Evolution equations at the macroscale . . . . . . . . . . . . . . . . . 142 8.4.1 Macroscopic Fick?s law with drift . . . . . . . . . . . . . . . . 143 8.4.2 Evolution equation in Cartesian system . . . . . . . . . . . . . 146 8.4.3 Change of variables . . . . . . . . . . . . . . . . . . . . . . . . 147 8.5 Extensions of macroscopic limit . . . . . . . . . . . . . . . . . . . . . 149 8.5.1 Exponential law for step chemical potential . . . . . . . . . . 150 8.5.2 Adatom desorption . . . . . . . . . . . . . . . . . . . . . . . . 151 8.5.2.1 Solution for microscale diffusion . . . . . . . . . . . . 152 8.5.2.2 Limit a? 0 . . . . . . . . . . . . . . . . . . . . . . . 154 x 8.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 9 Facets and boundary layers 158 9.1 Faceting on crystal surfaces . . . . . . . . . . . . . . . . . . . . . . . 159 9.2 Formulation of faceted relaxation as a free-boundary problem . . . . . 159 9.3 Natural boundary conditions of the variational problem for periodic profiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 9.4 Scaling of the boundary layer width . . . . . . . . . . . . . . . . . . . 169 9.5 Effect of an electric field on the slope profile near a facet . . . . . . . 174 9.5.1 Straight-step morphology . . . . . . . . . . . . . . . . . . . . . 174 9.5.2 Axisymmetric structure . . . . . . . . . . . . . . . . . . . . . 176 9.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 177 10 Numerical results 179 10.1 Review of the finite element method . . . . . . . . . . . . . . . . . . . 182 10.2 Software and computers . . . . . . . . . . . . . . . . . . . . . . . . . 184 10.3 Zero line tension (?ideal? case) . . . . . . . . . . . . . . . . . . . . . 185 10.3.1 Simulations with zero E-field . . . . . . . . . . . . . . . . . . . 187 10.3.2 Simulations with nonzero E-field . . . . . . . . . . . . . . . . . 203 10.4 Continuous dependence on line tension of the 2D to 1D transition . . 207 10.5 Conclusions and open issues . . . . . . . . . . . . . . . . . . . . . . . 210 11 Epilogue: conclusions and open questions 213 11.1 Summary of contributions . . . . . . . . . . . . . . . . . . . . . . . . 215 11.2 Open questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 216 A Brief review of the Mellin transform 220 B Nondimensionalization of the evolution equations 224 Bibliography 227 xi List of Abbreviations symbol meaning dimensions ? aspect ratio of a biperiodic height profile ?x, ?y characteristic lengths at the macroscale length a atomic height length h macroscale surface height length r position vector in the (x,y) plane length ? transverse step coordinate ?i transverse coordinate of the ith step ? longitudinal step coordinate ?? metric coefficient along ? length ?? metric coefficient along ? length mi discrete step density E surface free energy energy ? free energy density energy/length2 ? step stiffness energy/length g1 line tension coefficient energy/length2 g3 step interaction coefficient energy/length2 Ci adatom density on the ith terrace 1/length2 C macroscale adatom density 1/length2 Ceqi equilibrium adatom density at the ith step 1/length2 Cs straight step equilibrium adatom den- sity 1/length2 D, Dad terrace adatom diffusivity (bold- face=tensor) length2/time Ded edge atom diffusivity length2/time Jad terrace adatom flux 1/(length?time) Jed edge atom flux 1/(length?time) J macroscale adatom flux 1/(length?time) k, ku, kd kinetic rates of attachment/detachment length/time M macroscale mobility (boldface=tensor) length2/(energy?time) ? dimensionless mobility tensor ?i chemical potential of the ith step energy ? macroscale chemical potential energy q material parameter comparing diffusion to attachment/detachment T temperature temperature v adatom drift velocity length/time u nondimensionalized drift velocity E electric field energy/(charge?length) Z?e effective charge of adatoms charge ? atomic volume length3 xii Chapter 1 Introduction 1.1 Chapter overview In this chapter, we will introduce the basic principles of the near-equilibrium morphological evolution of crystal surfaces. Depending on the length scale under consideration, these surfaces present different physical pictures to the viewer. To capture the relevant physics of these experimental pictures, mathematical mod- els make certain simplifying assumptions about the operative mechanisms at each length scale. The different mathematical models are typically analyzed indepen- dently of one another. In this chapter, we emphasize the connections between mod- els at different length scales. These connections constitute the primary focus of the next few chapters. Finally, we outline the main results of this thesis, indicating how this work contributes to the ongoing effort to understand and exploit the physical mechanisms that govern the evolution of crystal surfaces. 1.2 Physical description of surfaces at various length scales There are three distinct scales for crystal surfaces: atomistic scale, nanoscale, and macroscale. A typical picture on the atomistic scale has individually distin- guishable atoms arranged periodically in a crystal lattice. The atomistic picture 1 of a surface has atoms with missing bonds at the solid-vapor interface. The more missing bonds these surface atoms have, the rougher the resulting height profile looks. At the nanoscale, the viewer can identify not the individual atoms, but only the extended surface features they form in aggregate. Just as a viewer from an airplane can see the Great Wall of China but not its individual stones, the nanoscale viewer of crystal surfaces can see the meandering lines formed by a large number of surface atoms at one height, adjacent to a large number of surface atoms at a different height. These line defects are called steps, since they appear at the boundary between two regions whose heights differ by a lattice constant. From the perspective of the nanoscale viewer, the steps resemble smooth curves. For this viewer, the fact that steps can be counted individually is the only clue that the microscope is showing an essentially discrete object. At the macroscale, all clues of discreteness are taken away, and the viewer is left with what looks like a smooth topography. Until the mid-20th century, this resolution was the only one accessible to experimentalists. However, the mental picture of a discrete lattice structure was always available, and it helped establish a framework for thinking about crystal surfaces. 1.3 Mathematical models of various length scales One possible mathematical approach to the atomistic scale is to seek approx- imate solutions of the many-body Schroedinger equation for electrons and nuclei 2 [121]. These approximate solutions are valid over prohibitively short time scales, which impedes most attempts to reconcile the motion of individual atoms with coarser descriptions of the surface. By allowing for atomic motions governed by classical mechanics, rather than quantum mechanics, Vvedensky and Haselwandter [121] sought a connection between atomistic models and the nanoscale. This con- nection lies outside our present scope. In this thesis, we study analytical and numerical aspects of the connection between nanoscale and macroscale. The nanoscale model assumes that steps move in response to (i) adatom attachment/detachment at step edges, (ii) diffusion of adatoms on terraces, and (iii) energetic interactions between steps. This description of the microscopic physics, introduced in a simplified form by Burton, Cabrera, and Frank (BCF) [10], serves as the basis for our present discussion of surface morphological evolution. The BCF model and its generalizations support a rich variety of step dynamics. The patterns that emerge during surface evolution have motivated analytical and numerical attempts to understand what experimentalists observe based on the behavior of steps. The temperature below which the BCF model applies is called the rough- ening temperature, TR. Thermal wandering of steps is allowed below TR, but the steps themselves are stable objects, not momentary deviations from surface flatness. Above TR, creation of a step has zero free energy cost, and steps are created and destroyed spontaneously. Not only steps, but also point defects such as vacancies can be present on the surface. This roughening indicates that a description of the surface in terms of steps is not appropriate. A rough surface is more easily visualized 3 Figure 1.1: The stepped surface of Si(001) as seen through an STM, taken from [99]. at the atomistic scale, where bonds between atoms at the surface are broken and formed at random. The macroscopic theory of Herring and Mullins [34, 75], based on continuum thermodynamics and mass conservation, accurately predicts the observed behavior of surfaces above the roughening temperature. However, the fabrication of crys- tal surface morphologies is no longer restricted to temperatures above TR, which challenges the validity of the continuum theory of Herring and Mullins. Following the trend toward miniaturization of semiconductor devices, the length scales on which crystal surfaces are patterned have grown smaller, neces- sitating lower temperatures so that the desired pattern is more stable. These lower temperatures require a macroscale theory consistent with the flow of thermally sta- ble steps. Figure 1.1 shows the surface of Si(001) as seen through a Scanning Tunneling Microscope (STM). We can identify terraces where the surface is locally flat, except 4 Figure 1.2: Schematic of curved crystal steps on a vicinal surface. The arrows represent the different physical processes through which adatoms affect the flow of steps. for a few isolated bumps and vacancies of size about 3 ?A(one lattice constant). The line defects bounding the terraces are steps about one lattice constant high. The widths of terraces range from one to a thousand lattice constants. The presence of kinks and edge atoms, visible in Figure 1.3, reminds us that continuum steps are onlyan approximation of a fundamentallydiscrete phenomenon. To use an image with the same resolution as Figure 1.3 to generate initial conditions for a simulation of continuum steps, one would have to smooth out the kinks and edge atoms. Jeong and Williams [43] suggest the use of an integral mollifier for this purpose. The result is a representation of the surface by smooth steps as in Figure 1.2. In Figure 1.2, the relevant near-equilibrium physical processes are (i) adatom diffusion on terraces, (ii) diffusion of edge adatoms along steps, and (iii) attach- ment/detachment of adatoms at step edges. Diffusion processes such as (i) and (ii) arise when each terrace or edge adatom is assumed to perform a random walk; the 5 Figure 1.3: Image of a silicon (100) surface, taken from [99], showing kinks and edge atoms over many steps. Solid curves along the diagonal have been superimposed to illustrate the approximation of smooth steps. Note two kinds of steps on this surface. 6 resulting step-continuum equations represent the homogenized (averaged) effect of a large number of these adatoms moving stochastically throughout their domains. The overall evolution of nanoscale averages, such as the adatom density, is predictable according to the law of large numbers. The kinetic process (iii) is a manifestation of a near-equilibrium law, whose forcing term is similar to those that appear in dif- ferential equations for chemical reaction rates. In the case of surface steps, the rate at which adatoms attach or detach at step edges is proportional to a difference of concentrations. A step advances when adatoms attach to the step edge; it retreats when adatoms detach from the step edge. We now discuss a few of the atomic processes that change the shape of a step. A curved step, such as those depicted by the solid diagonal paths in Figure 1.3, actually features kinks that ?zig-zag? along the step edge. The expansion of kinks can result from nucleation of terrace adatoms, if the difference of concentrations is favorable to adatom attachment at the ends of a kink. The opposite effect, when concentration differences favor adatom detachment, consists of a step emitting adatoms. This change in the step shape reduces the free energy associated with unbonded edge atoms. Such energy-minimizing mass transport can be interpreted as an effect of step line tension. In general, a step can change its shape by adsorbing or emitting adatoms. It remains to explain how the equilibrium shape of a step is communicated to the adatoms. Absorption and emission of adatoms depend on the differences between con- centrations. One of these, the terrace adatom concentration, is in principle directly observable. The other concentration is an equilibrium quantity that depends on the 7 surface geometry via the step curvature and energetics of interacting steps. Intu- itively, steps with high curvature and small separations from their nearest neighbors have higher energy than straight steps widely separated from their nearest neigh- bors. We expect a highly curved step to require more energy because a greater number of broken bonds are associated with regions of high curvature, as one can see in Figures 1.1 and 1.3. The effect of step separation is more subtle and involves two different kinds of interaction, as we see shortly. The energy of a step train is conveniently analyzed using the concept of step chemical potential. The step chemical potential, reflecting the contribution of step line tension and step-step interactions, measures the propensity of a step to incor- porate adatoms. This physical interpretation follows from the monotone relation between adatom equilibrium density and step chemical potential; see the discussion around (3.10). Interactions between steps of the same sign (both either step-up or step-down) are of two types, entropic and elastic-dipole. Entropic repulsions enforce the non- crossing of steps, by restricting how far each step at nonzero temperature can deviate from its average position. This thermal wandering is an unavoidable consequence of nonzero temperature. However, the presence of other steps imposes a penalty for any wandering that brings a step into collision with one of its neighbors. This penalty takes the form of a loss in entropy, making the step?s free energy greater than it would be in the absence of neighboring steps. Elastic-dipole interactions are mediated by the strain fields that steps produce in the crystal [63, 76]. Any defect on the surface, including a step, induces a force distribution that in principle 8 has a multipole expansion of all orders. Comparison between predictions of this multipole expansion and data from classical atomistic simulations led Najafabadi and Srolovitz [76] to conclude that the dipole-dipole interaction is dominant over higher-order contributions for the typical range of terrace widths. The combination of elastic-dipole and entropic interactions yields the step chemical potential as the superposition of two terms, each term reflecting an inverse-square dependence on the distance between the step and one of its neighbors. This superposition is not linear, appearing in the equation for interaction energy as a renormalized parameter for the strength of repulsions [43]. The approximation of continuum steps, as noted above, deliberately glosses over the details of point defects on the steps themselves. Kinks and edge atoms, always present to some extent on steps at nonzero temperature, are not directly incorporated into the formula for step chemical potential. Some variations of the BCF model account for kinks and edge atoms by introducing auxiliary fields, such as kink density. Equations for these auxiliary fields are coupled with the other equations for step motion in a way that respects mass conservation [3, 11, 68]. We do not treat coupled field equations in this thesis, instead opting for the simplicity of a semi-empirical formula for the step chemical potential. By ?semi-empirical? we mean that the scaling of the step chemical potential with nearest-neighbor distance can be derived analytically, but in general the prefactor can only be approximated by comparisons with data, numerical or experimental. 9 1.4 Connecting step models to experiments The feedback loop between step models and experiments began with an obser- vation that puzzled the mid-20th century crystal growth community. According to the accepted models of crystal nucleation and growth in solution, the concentration of crystal atoms had to be sufficiently large in order for crystals to precipitate on the growing surface. When the threshold concentration was observed to be much lower than predicted, Burton, Cabrera, and Frank [10] hypothesized a new mech- anism for crystal growth in solution to explain the discrepancy. This explanation proposed the existence of steps, which serve as convenient nucleation sites for the atoms adsorbed from solution. In the BCF model, step motion is governed by the diffusion of adatoms across terraces and the attachment/detachment of adatoms at step edges. Their formulation of equilibrium density differs from the one we adopt following [37, 64, 113], in that the BCF model omits step interactions. The simple picture inspired by the BCF model enabled many successful predictions of surface phenomena, such as Spohn?s prediction of step density in the neighborhood of a facet [109], and the inverse linear decay of peak-to-valley distances found by Rettori and Villain [98] for certain biperiodic profiles. More recently, the BCF model has been expanded to account for step-step in- teractions and anisotropies in the material parameters. These extensions of the BCF model are often introduced to provide better interpretations of the observations of crystal surfaces. For example, the experimental determination of certain coefficients in the formula for surface free energy leads to interesting analogies between the ter- 10 race width distribution and other probabilistic phenomena. The successful interplay between step models and experiments confirms that step-based models still hold a valuable place among the various approaches to surface science. 1.5 Macroscale limit Advances in imaging techniques since the early investigations of epitaxial phe- nomena have expanded the range of time and length scales over which we can study crystal surfaces. Today with Atomic-Force Microscopy (AFM) and other imaging devices, it is common to observe samples with characteristic lengths of the order of 100 ?m or larger for durations of minutes, hours or more. Samples of this size comprise hundreds or thousands of steps, and the AFM probe will not resolve all the details fast enough to follow the individual steps as the surface relaxes. A meaning- ful connection between observations at this scale and the predictions of a step-based theory requires a systematic way to bridge the gap between the limited capabilities of an AFM probe and the difficulties of understanding large-scale surface morpho- logical evolution via the flow of steps. In applications it is also often desired to control the evolution of crystal surfaces, and this control is more readily achieved at larger scales. Because of the below-roughening temperatures under which modern small devices are made, the presence of steps at the nanoscale must be taken into account when attempting to control the surface evolution at the macroscale. It turns out that we can bypass a complete simulation of step motion in 2+1 dimensions by using a formal procedure known as coarse-graining. This method 11 regards the equations of step flow as a discrete scheme for some evolution law gov- erning macroscale quantities, such as the surface height. In the same way that decompilation of computer code reconstructs a plausible version of the higher-level source code, coarse-graining takes the low-level ?instructions? given by a nanoscale model and determines a consistent macroscale ?program? for which the nanoscale equations are one implementation. This analogy fails in several ways, most impor- tantly by suggesting that the output of a coarse-graining procedure is not unique (up to isomorphism), since different computer code decompilers can yield structurally different source code programs. In fact, coarse-graining is a well-defined operation, producing a unique partial differential equation (PDE) given a deterministic system for the finer-scale variables. The PDE generated by a coarse-graining operation describes the relation be- tween space and time derivatives of the surface height. This PDE is also known as the full continuum limit (or macroscale limit) of the nanoscale equations. This terminology reflects the facts that (i) a limit is taken when going to the coarser scale, and (ii) the resulting PDE is a fully continuum equation, making no explicit reference to the nanoscale variables from which it was derived. After verifying the hypotheses that justify using the macroscale limit, one might make direct compar- isons between the predictions of this PDE and the experimental results. Inthecaseofsurfacerelaxation, themacroscalelimitisafourth-orderparabolic nonlinear PDE for the height decay in the region of space where ?hnegationslash= 0 [37, 64, 96]. If the complementary region {x : ?h(x) = 0}?i.e., the set of macroscopically flat facets?has positive measure, then we need to determine the appropriate bound- 12 ary conditions at facet edges to obtain a well-posed problem from the macroscale limit [31, 109]. The choice of boundary conditions must weigh the sometimes con- flicting demands of easy numerical implementation on the one hand, and accurate microscale physics on the other hand. 1.6 Overview of this thesis The main results of this thesis are fully continuum equations for the surface height based on step models, and investigation of these equations from numerical and analytical perspectives. In Chapter 2, we introduce the basic ideas common to all macroscopic theories, and we review the numerical method we apply later to approximate the solutions of macroscale equations. In Chapter 3, we introduce the assumptions of near-equilibrium step flow in the context of one-dimensional step trains, to lay the groundwork for subsequent generalizations in 2+1 dimensions. In Chapter 4, we take a brief detour to derive the formula for elastic interaction energy between steps. The generalization of this formula is used in Chapter 5, where we review the (2+1)-dimensional macroscopic limit of [66]. The modeling component of this thesis extends the (2+1)-dimensional setting of [66] in several ways. First, the diffusion of adatoms on terraces is allowed to be anisotropic, which allows for a more realistic representation of surfaces where adatoms hop more readily perpendicular to steps than parallel to them, or vice versa. Second, we include the effect of step edge diffusion, where mass transport along the step edge is modeled as the result of longitudinal variations in the step 13 chemical potential. The macroscale limit with anisotropic terrace diffusion and step edge diffusion is presented in Chapter 6. Different scalings of approximate, separable solutions of the PDE are made possible by the additional parameters corresponding to step edge diffusion and anisotropic terrace diffusion. In Chapter 7, we generalize the boundary conditions at step edges to address the possibility of permeable steps. Adatoms approaching a permeable step are free to hop directly to the next terrace. In principle, step permeability can lead to nonlocal effects, as adatoms diffuse over many terraces before attaching to a step edge. The macroscale limit, however, reflects step permeability only through a renormalized kinetic rate. Our derivation in 2+1 dimensions generalizes the work of Pierre-Louis [88] for straight steps. Another extension of (2+1)-dimensional step models introduces a drift velocity in the terrace diffusion equation, which corresponds physically to an applied elec- tric field. The study of crystal surfaces evolving under the influence of an electric field began with an experiment by Latyshev et al. [59], who reported the onset of step bunching. This surface instability is characterized by wide terraces between clusters of closely-spaced steps. Electromigration-induced step bunching has been used as a practical method for fabricating quantum wires. More recently, the appli- cation of an electric field parallel to step edges has been shown to produce a step meandering instability. This phenomenon features the growth of harmonic pertur- bations of straight steps, as predicted by Dufay et al. [15] and observed by Degawa et al. [14]. Motivated by these possible applications, we introduce in Chapter 8 the physics of electromigration at the nanoscale. Then we use the BCF-type model with 14 electromigration to derive the corresponding macroscale equation. In Chapter 9, we discuss the inclusion of facets in the macroscopic theory. A re- laxing surface with a facet can be modeled as a free-boundary problem [64, 109]. The influence of step-step interactions on the slope profile depends on the distance from the facet edge. We apply boundary-layer ideas, similar to those employed in [64], to determine how the slope profile behaves near a facet edge in 2+1 dimensions. The boundary-layer solution is qualitatively different if we introduce an appreciably large electric field. We study this effect analytically in the case of straight steps, axisymmetry, and also more general step geometries with slowly varying curvature. In Chapter 10, we use the finite element method to simulate the macroscopic equations. Numerical simulations of the PDEs derived in Chapters 5?8 illustrate the rich variety of relaxation phenomena predicted by the model. The first noteworthy result, evident even in the case of isotropic terrace diffusion, is the morphological transition exhibited by an initially biperiodic profile, which relaxes to become almost one-dimensional [7]. We discuss the plausible mechanisms behind this transition. To isolate the effect of longitudinal fluxes, we simulate the macroscopic equation with electromigration for an initially biperiodic profile. In this case, the acceler- ated transition to an almost one-dimensional profile offers insight into the driving mechanism behind the observed morphological changes. 15 Chapter 2 Macroscopic equations and numerical solutions In this chapter, we provide the necessary background for the subsequent in- vestigation of macroscale crystal surface evolution. Our research relies on PDEs for the macroscale surface variables, as well as numerical simulations of these PDEs. These components of our research program build on previous work in modeling and numerics. In the following review of the relevant background material, we make special note of the limitations of previous macroscopic theories and simulations, in order to motivate the new topics of the following chapters. 2.1 Macroscopic theories of crystal surfaces In a fully continuum picture of crystal surfaces, we smooth out the kinks and approximate the discrete height by a continuous graph. The immediate advantage of this approach is the reduced number of dependent variables. The height profile in the macroscopic model is a single function of time and the two spatial variables. The PDE satisfied by the macroscale height function lends itself more readily to numerical simulations on modest hardware, whereas atomistic models require much more memory to store the positions of all the atoms in a sample. Fully continuum theories also have the advantage of enabling global predictions, such as scaling laws and similarity solutions, more easily than discrete models. These global predictions 16 are better suited for comparisons with experimental results, because the time scale of experimental observation is realized by macroscopic simulations but is typically unattainable by atomistic methods. Macroscopic theories are not without their limitations, however. One should be cautious in applying a fully continuum theory when detailed resolution of atomistic or nanoscale processes is desired. The configuration of kinks and edge atoms, for example, cannot be recovered from a macroscopic simulation. Moreover, kinetic rates manifest themselves in macroscopic theories only through combinations of material parameters. To determine any particular kinetic rate is outside the scope of our fully continuum models. 2.1.1 Effect of different mass transport mechanisms In total, four different physical mechanisms have been identified for mass trans- port in a crystal material. In [33], Herring considers separately the processes of (i) viscous flow, (ii) evaporation and condensation, (iii) volume diffusion, and (iv) surface diffusion. Each of these processes leads to a different scaling between the lifetime and size of a surface feature. For example, two geometrically similar sur- face features with characteristic length scales L1 and L2 will change their shape by evaporation/condensation with lifetimes ?1 and ?2, which satisfy ?2/?1 = (L2/L1)2. The exponent 2 appears because evaporation/condensation works at the solid-vapor interface, where total surface area (proportional to the square of the characteristic length) is relevant. In general, a surface feature with characteristic length ?L1 de- 17 caying via one of the mechanisms (i)?(iv) has a lifetime ? given by ? = ?1?p, where p ? {1,2,3,4} is associated with the operative mechanism of mass transport. The largest possible power, p = 4, is associated with surface diffusion. At a given temperature, any of the four possible mass transport mechanisms may be present in a crystal material. The lifetime of a structure varies most rapidly with characteristic size when surface diffusion is present. The trend toward fabricat- ing ever smaller surface features, as required for mobile electronic devices, would im- pose severe time constraints at the temperatures commonly used in previous decades for crystal surface patterning. For example, a miniaturized surface feature, half as large as one produced currently, might decay 16 times faster at the same tem- perature. To mitigate the difficulties of shorter lifetimes in structures subject to surface diffusion, experimentalists and manufacturers have brought their equipment to lower operating temperatures. At these lower temperatures, the continuum ther- modynamic description of crystal surface evolution is challenged. The appropriate theory, based on the motion of steps below roughening, is thus enjoying a resurgence in the research community. 2.1.2 Macroscopic theory above the roughening temperature Before describing the extension of macroscale theories to temperatures below roughening, we review the ingredients of surface evolution models that first appeared indescriptions above therougheningtemperature. Macroscale descriptionsofcrystal surfaces above roughening apply the ideas of continuum thermodynamics and mass 18 conservation, obtaining PDEs or variational principles for the surface height. A classic example is the fourth-order PDE, ?th? div?(?h), (2.1) derived by Mullins and Herring. This result holds if 1. the chemical potential ? is proportional to the mean curvature, ???1 +?2, 2. the surface flux is proportional to the gradient of chemical potential, J ???, and 3. mass conservation is enforced via ?th??divJ. In this derivation, Mullins introduces only one chemical potential,?, which is associ- ated with the surface geometry through the mean curvature. A second chemical po- tential for the vapor phase would be needed in the case of evaporation/condensation, but Mullins assumes only surface diffusion in his analysis. As noted in Chapter 1, it is reasonable to expect the mean curvature to determine the likelihood of adatom attachment/detachment on the surface. In the regime of small slope (|?h| lessmuch 1), the mean curvature definition of chemical potential yields ?? ?h, (2.2) where ? is the 2D Laplacian. An alternative formulation of the chemical potential is in terms of the surface energy E, via the variational formula ? = ?E?h, (2.3) 19 where E = integraldisplay R2 ?(|?h|)dx, ? = g0 + 12g2|?h|2 +o(|?h|2). (2.4) We remark that the above-roughening free energy density? is smooth for any surface slope |?h|. The relation between surface flux and the gradient of chemical potential is the macroscale analog of Fick?s law1. At the level of steps and terraces we have a flux Jad defined by Jad = ?Dad ??C, where C is the adatom concentration and Dad is the adatom diffusivity. At the macroscale we have the similar relation J ??M???, (2.5) where M is the mobility. Both Dad and M are in principle tensor-valued functions of the slope, restricted on physical grounds to be positive definite. The case of scalar Dad and M is perhaps easiest to understand intuitively, but we shall soon encounter regimes where this intuition is challenged. Finally, the law of mass conservation is written so that a net flux of adatoms across a level set for h is exactly balanced by the local time derivative of the height. In the absence of deposition, the material derivative (or total derivative) of the level set h = const. yields the mass conservation law, ?th??divJ. (2.6) When deposition is present, a forcing term must be added to the right hand side of (2.6). This case is not studied in this thesis. 1See the discussion around (5.3). 20 2.1.3 Macroscopic theory below the roughening temperature The three ingredients of macroscale models of surface morphological evolu- tion, first introduced in the study of materials above roughening, also carry over to temperatures below TR. Care must be taken in writing the analogs of (2.2)?(2.6) below roughening, in order to respect the presence of steps. For example, according to the BCF model, there are N different step chemical potentials for an array of N steps. It is not immediately obvious how one would obtain a fully continuum chemical potential ? in this setting, much less a relation analogous to (2.2). Sim- ilarly, the dependence of terrace flux Jad on adatom concentration is not readily connected with the macroscale relation (2.5). In Chapter 5 we review the details of this macroscale limit procedure for the simple case of isotropic terrace diffusion, where Dad is a scalar constant. The idea of a macroscale limit has a long history in the literature, producing PDEs based on BCF models alone [2, 51, 86, 94, 119], BCF models coupled with adatom density [16, 62, 104], or kinetic Monte Carlo models [120]. We follow the first approach and consider BCF models alone, allowing for some interesting gen- eralizations to account for different possible physics at the nanoscale. In Chapter 6 we include the effects of step edge diffusion and anisotropic terrace diffusion. In Chapter 7 we derive the macroscale limit for permeable steps, which allow adatoms to jump from one terrace to its neighbor without first attaching to the step between them. In Chapter 8 we include the effects of a drift velocity caused by an applied electric field, and desorption of adatoms from a terrace into the vapor phase. These 21 effects are easily accommodated in the below-roughening analogs of (2.2)?(2.6) with- out disrupting the structure of the macroscopic equations. In every case we study here, the macroscale limit can be written in the form ?th = div(M???), (2.7) ? = ?E?h (2.8) for some appropriate mobility M and surface energy E. Below roughening, the singular surface energy E can be expanded as a series in |?h|: E = g0 +g1|?h|+ g33 |?h|3. (2.9) By interpreting |?h| as the macroscale analog of the discrete step density, the co- efficients of the expansion for E reflect the relative contribution of energies (i) pro- portional to the number of steps (g1 term), (ii) independent of the number of steps (g0 term), or (iii) stemming from interactions that decay as the inverse square of step separation (g3 term). Because ? is not directly observable, but only serves as an auxiliary variable in the fully continuum PDE, it will sometimes be convenient to change variables so that the variational structure (2.7),(2.8) remains intact. 2.2 Numerical simulations In this section, we discuss the numerical treatment of evolution equations for the macroscale surface height. Our approach takes advantage of the variational structure (2.7),(2.8) to implement a weak form of the macroscale PDE using finite 22 elements. If the mobility M and surface energy E are chosen appropriately, we expect the solution of the macroscale PDE to approximate the near-equilibrium motion of steps. However, a full comparison between the predictions of nanoscale and macroscale models is outside the scope of this research. Lacking a numerical implementation of step flow equations in 2+1 dimensions, the only data available for comparison come from simulations of the coupled differential equations for the positions of straight or concentric circular steps [25, 45]. We now turn to the implementation of fully continuum models for the surface height. To simulate these macroscale PDEs, we must decide on a discretization, which introduces computable quantities to approximate the fully continuum data. The vector of step positions is, of course, one such discretization. However, the equations of step flow are in general too unwieldy to simulate directly. We take advantage of the flexibility afforded by the fact that any given differential equation has more than one possible discretization. Amongthepossiblediscretizationchoicesavailabletous, twoinparticularhave proven useful for the simulation of crystal surface relaxation. The first (Galerkin) method approximates the surface height by an expansion in trigonometric basis functions [108], which yields coupled differential equations for the coefficients after substitution in (2.7),(2.8). This system of ordinary differential equations is then integrated numerically. Due to its easy formulation and capacity to accommodate higher-order discretizations in time, the Galerkin method is a natural choice for simulating the evolution equation. We adopt the finite element method, which partitions the domain into non- 23 overlapping elements (e.g., triangles) and approximates at time n the surface height hn and the macroscale chemical potential ?n as a sum of basis functions ?i (some- times called hat functions), each of which is supported on a finite, contiguous cluster of elements, viz., hn = summationdisplay i Hni ?i, (2.10) ?n = summationdisplay i mni?i. (2.11) The approximations (2.10),(2.11) can only resolve the spatial dependence ofhand? over distances not less than the mesh size. By choosing the mesh size appropriately, we can resolve all the physically meaningful spatial variations in h, without intro- ducing unmanageably many degrees of freedom in the triangulation T . We denote by VT the finite element space spanned by the hat functions ?i for the triangula- tion T . The weak formulation of the evolution law, expressed in the space VT , will finally allow us to find approximate solutions of (2.7),(2.8). In order to write this weak formulation, we must first choose a discretization for the time derivative. The finite element method works by solving a linear system for the coefficients Hni and mni in (2.10),(2.11). The exact definition of ? in terms of h is, in general, nonlinear, so the first task of a time discretization is to find an approximate, linear version of the relation (2.8). We choose a semi-implicit Euler time step [32, 77], replacing (2.7) by hn ?hn?1 ?t = divM(h n?1)???n, (2.12) ?n = ?E?h vextendsinglevextendsingle vextendsinglevextendsingle hn?1 . (2.13) 24 This discretization allows us to find the unknown coefficients at any time by solving a matrix equation. To determine which linear system the coefficientsHni and ?ni satisfy, we take the inner product of (2.12) with an arbitrary test functiong ?VT , and also the inner product of (2.13) with an arbitrary test function? ?VT . Because VT is spanned by the basis functions ?j, it suffices to consider inner products with each ?j individually. The linear system that emerges is summationdisplay i (Hni ?Hn?1i ) integraldisplay B ?i?j = ?t summationdisplay i mni integraldisplay B M(hn?1)??i ???j, (2.14) summationdisplay i mni integraldisplay B ?i?j = integraldisplay B ?E ?h vextendsinglevextendsingle vextendsinglevextendsingle hn ?j, ??j, (2.15) where B is the domain over which the PDE (2.7),(2.8) is assumed to hold. The integral on the right-hand side of (2.15) is more commonly rewritten using integration by parts, as we shall see in Chapter 5 when considering the formula for macroscale chemical potential in more detail. Integration by parts has already been performed on the right-hand side of (2.14), with vanishing boundary terms2. We alert the reader that the surface energyE has a singularity at|?h| = 0. For the integration in (2.15) to make sense, we add a regularization parameter ? lessmuch 1 to the denominator |?h| that appears in the g1 term of ?E/?h. The qualitative behavior of the approximate solution is not significantly affected by ? [7]. Finally, the tensor mobility in the fixed coordinate system might also have a singularity at |?h| = 0. Different regularization schemes are possible to remove this singularity. For mathematical consistency we might wish that the regularized 2Namely, ?th and the flux J are both zero on the boundary of B; or, we apply periodic boundary conditions (assuming steps lie on a torus in R2). 25 mobility approaches a scalar as ?h ? 0. However, we find that the approximate solution is rather robust with respect to the mobility regularization, as discussed in Chapter 10. A number of linear solvers are available to compute the coefficients Hni , ?ni using (2.14),(2.15). Once a solution is obtained, we can quantify its usefulness by estimating the size of the truncation errors and interpolation errors (due to dis- cretization in time, the choice of mesh, and the choice of polynomial space for the basis functions). In the present work, we do not focus on a posteriori analy- sis in this vein. Our emphasis is on the global behavior of solutions to the PDE (2.7),(2.8), which we expect to be well captured by any reasonable numerical solution of (2.14),(2.15). The preceding illustration of the finite element method gives only a flavor of the mathematics behind this powerful tool. A much richer description appears in [9]. Evans [24] serves as a good reference for the necessary background in Sobolev spaces. 26 Chapter 3 One-dimensional nanoscale models of interacting steps In this chapter we review an extended version of the BCF step model. This classical picture of step motion, updated to include entropic and elastic-dipole step interactions, describes the relevant nanoscale physics of stepped surfaces. Here nanoscale reminds us that the motion of steps, which we use to derive a macroscale limit, is itself a coarse approximation of atomistic processes. 3.1 Geometry of straight steps We consider a descending train of steps with atomic heighta, separated by flat regions called terraces; see Figure 3.1. The step positions xi are treated as moving boundaries for the adatom diffusion of each terrace. To account for the fact that no overhangs or step crossings are allowed, we require that at each time t, the step position xi(t) increase monotonically with the step number i, xi+1(t) > xi(t). If this condition is satisfied for the initial step positions xi(0), then the subsequent evolution is expected to respect the monotonicity, due to the repulsive step-step interactions introduced in the equations below. We regard the monotonicity of step positions for t > 0 as a basic assumption, although an elementary proof is 27 Figure 3.1: Geometry of straight steps. possible for certain mass transport mechanisms.1 For each time t, the ith terrace {x : xi < x < xi+1} is a well-defined domain on which we seek solutions of the equation of motion for adatom diffusion. For the present case of one independent spatial variable, the elliptic problem corresponding to adatom diffusion reduces to an ordinary differential equation that admits an explicit solution. 3.2 BCF model for interacting straight steps A quantitative discussion of the BCF theory begins with the adatom density, Ci, on the ith terrace, xi < x < xi+1. This Ci satisfies the diffusion (or heat) equation [10], ?tCi = ?x(Dad ??xCi) , (3.1) where Dad is the diffusivity for adsorbed atoms on the terrace. Note that we have omitted from (3.1) terms that describe atom desorption, electromigration and ma- terial deposition from above, considering only the effect of surface relaxation. Equa- tion (3.1) can be simplified significantly by elimination of the time derivative, for the following reason. The time scale for step motion is much larger than the time scale 1For instance, by considering mass transport dominated by evaporation-condensation, we obtain a second-order PDE for the slope, to which the maximum principle applies. 28 for terrace diffusion. Thus, we simplify (3.1) via the ?quasi-steady approximation?, ?tCi ? 0; the time dependence in Ci enters through the boundary conditions at step edges. Solving the resulting differential equation for Ci by successive integrations yields Ci = Ai(x?xi) +Bi, xi 0 measures the strength of repulsive elastic interactions between steps. Another reason for steps to repel each other is the fact that steps are thermally rough. Just as steps are line defects on the surface, steps themselves can have point defects (edge atoms) or kinks, which slightly perturb their straight shape. The formation of an edge atom requires only the energyepsilon1? associated with broken bonds, while the creation of an entire step requires energy in proportion to its length. At any nonzero temperature T, it is more energetically favorable for a step to wander (due to formation of kinks and edge atoms) than to advance or retreat as a whole. This thermal wandering does not proceed unimpeded. Invoking the non-crossing condition for steps, Gruber and Mullins [30] observed that collisions with neighbors would restrict the thermal wandering of a given step in an array of parallel steps. This restriction leads to a smaller entropy for the step in question, which then has larger free energy. The overall effect is that of entropic repulsion between steps. To quantify the entropic interaction between steps, it is helpful to consider the in-plane correlation function ?|x(y)?x(yprime)|2? = kBT?? |y?yprime|, (3.21) where the angle brackets denote a thermal average with respect to the equilibrium distribution, the steps are aligned with the y-axis, and ?? is the step stiffness derived from the Ising lattice gas model [73], ?? = kBT a e epsilon1?/kBT. (3.22) Neighboring steps are assumed to have an average separation w (measured in the x-direction), and the average distance between collisions is Lc (measured in the 35 y-direction). Setting the left hand side of (3.21) equal to w2 and substituting Lc for |y?yprime|, we find Lc ?w2??/(kBT). If each collision decreases the entropy by a fixed amount kBC, then the step?s free energy per unit length is increased by the amount ?? ? CkBTL c = C(kBT) 2 w2?? . (3.23) To determine the coefficient C, it is useful to think of the steps as world lines of non-interacting fermions, with the y-axis playing the role of time [42, 44]. This representation yields C = pi2/6. The decay of entropic repulsion strength asw?2 is the same scaling we found for elastic-dipole interactions. The combination of both forces is easily accommodated in (3.20) by renormalizing the coefficient g. Here we consider ?i stemming from nearest-neighbor step interactions. Thus, ?i depends only on the distances |xi ?xi?1|. (In the macroscopic limit, adding up all the pairwise interactions has the same effect as renormalizing the parameter of interaction strength when considering only nearest-neighbor forces [66].) To deter- mine the chemical potential of the ith step, it suffices to compute the elastic-dipole forces exerted by neighboring steps only, not taking into account any interactions beyond nearest-neighbor. Suppose the ith step advances by an amount ?xi. This step advance requires the attachment of an additional ?N = a?xi/? atoms per unit length, where ? is the atomic volume. By analogy with the usual definition of chemical potential, we define the step chemical potential as the change in interaction energy upon addition 36 or removal of ?N atoms per unit length. Formally letting ?xi ? 0, we write ?i := ?a ddx i [Vi,i+1 +Vi,i?1]. (3.24) The formulas (3.20) and (3.24), together with (3.8) and (3.9), constitute a system of ordinary differential equations governing the flow of straight steps. A typical question of interest for this system is whether an equilibrium solution (i.e., that for equally-spaced steps or constant surface slope) is stable. If not, we try to characterize the types of instabilities that emerge as the surface relaxes. Stability analysis of the step flow equations in this geometry employs techniques from perturbation theory, such as the phase equation method [49], which studies the evolution equation satisfied by the phase ? of a given initial pattern. Research in this vein has led to general conclusions about the properties of steady state solutions for a large class of nonlinear PDEs [95]. The focus of our work is on step systems in 2+1 dimensions, where such a stability analysis is not practical. The reader is directed to [15, 95] for more details about the dynamics of (1+1)-dimensional step trains. 3.4 BCF model for interacting concentric circular steps The restriction that step positions be expressible in terms of only one spa- tial variable allows yet another geometry. Whereas the geometry of straight steps emerges from translational invariance, we obtain the geometry of concentric circu- lar steps (Figure 3.2) by requiring rotational invariance. Polar coordinates (r,?) are appropriate in this setting, where the dependence on the angular variable ? is 37 Figure 3.2: Geometry of concentric circular steps. eliminated due to the assumed axisymmetry. The point r = 0 is singled out as the center of the mound or valley. To establish some parallelism with the case of straight steps, we restrict our attention to axisymmetric structures where the height is a decreasing function of r (mounds). At any position inside the innermost step, in particular at r = 0, the height is a constant; h(r,t) = hf(t) for 0 0, r(s) = O(1), ?(s) = O(1), and epsilon1 lessmuch 1. Here epsilon1?1 denotes 57 Figure 4.1: Geometry of neighboring steps without axisymmetry, parametrized by r1(s) and r2(s). We define s := sin(?/2) to convert trigonometric functions of ? into algebraic functions of s. 58 the characteristic length scale, roughly analogous to the radius of circular steps. In addition, we assume ?(s),r(s) ?C1(?1,1), so that the tangent vector at each point on the two steps is well defined. A crucial quantity to calculate is the distance R between points on the two steps, which enters as R3 in the denominator of the elastic interaction energy inte- grand. We have R2 = r2(s)2 +r1(0)2 ?2r1(0)r2(s)cos? (4.53) = epsilon1?2{[?(s)??(0)]2 + 4s2?(0)?(s)}+ 2epsilon1?1?r(s)[?(s)??(0) + 2?(0)s2] +?2r(s)2. (4.54) Because ?,r are assumed to be C1, using Taylor series about s = 0 we can define the C0 functions ?,r by ?(s)??(0) =: s?(s), (4.55) r(s)?r(0) =: sr(s). (4.56) After adding and subtracting ?(0),r(0) as needed, we find R2 = epsilon1?24?(0)2s2 + 4epsilon1?1?r(0)?(0)s2 +?2r(0)2 +epsilon1?2[?(s)2 + 4s?(s)?(0)]s2 + 2epsilon1?1?s{r(s)[s?(s) + 2?(0)s2] +r(0)?(s)} +?2s[2r(0)r(s) +sr(s)2]. (4.57) In order to factor out the coefficient of s2, we define A(s) := {4[epsilon1?2?(0) +epsilon1?1?r(0)]?(0)}?1? parenleftbigepsilon1?2[?(s)2 + 4s?(s)?(0)]s+ 2epsilon1?1?{r(s)[s?(s) + 2?(0)s2] +r(0)?(s)} +?2[2r(0) +sr(s)]r(s)parenrightbig. (4.58) 59 Then R2 = 4[epsilon1?2?(0)+epsilon1?1?r(0)]?(0) braceleftbigg s2 +sA(s) + ? 2r(0)2 4?(0)[epsilon1?2?(0) +epsilon1?1?r(0)] bracerightbigg , (4.59) where A(s) ?O(1). To follow as closely as possible the notation of the axisymmetric case, we define the parameters ?2 := ? 2r(0)2 4?(0)[epsilon1?2?(0) +epsilon1?1?r(0)], (4.60) ?2 := 4[epsilon1?2?(0) +epsilon1?1?r(0)]?(0), (4.61) which feature in the asymptotic expansion of the integral for elastic interaction en- ergy. In terms of?and?, the squared distance between points with (s,r) coordinates (0,r1(0)), (s,r2(s)) is R2 = ?2{s2 +?2 +sA(s)}. (4.62) We note that the breaking of axisymmetry enters here through the termsA(s), which accounts for the curvature variation when steps deviate from circles. The geometry of Figure 4.1 is deceptive in its simplicity. We have chosen a coordinate system such that the vector normal to step 1 ats = 0 is aligned with thex-axis. By following the procedure of the axisymmetric case, we still obtain an elastic interaction energy per unit length, but now we cannot immediately generalize the resulting formula to other points along step 1. In principle, a different calculation with different definitions of A(s),?, and ? is required for each point on step 1. Despite this shortcoming, we expect the qualitative features of our interaction energy formula, in particular the inverse-square dependence on local step separation, to persist no matter which 60 point on the inner step is used to orient the Cartesian coordinate system. Below we outline the calculation for reasonably general 2D steps. A complete exposition of all the details would distract from the main point of this thesis (and chapter) and is hence omitted for the sake of clarity. For simplicity, we write the integral over step 2 (which in principle requires s to range over [?1,1]) as two separate integrals, dWint r1d?1 = 1??2 piY P 2 parenleftbiggintegraldisplay pi 0 + integraldisplay 0 ?pi parenrightbigg sin?dr 2/d?+r2 cos? (r21 +r22 ?2r1r2 cos?)3/2d?, (4.63) where P is the dipole moment per unit length of an atomic-height step. The first integral considers the energetic contributions given by (4.11) for ? ranging from 0 to pi. The evaluation of the second integral is essentially not different and is therefore omitted. We now face the task of calculating I2D = integraldisplay pi 0 sin?dr2/d?+r2 cos? (r21 +r22 ?2r1r2 cos?)3/2d?, (4.64) The numerator in (4.64) stems from the inner product of the two vectors normal to step 1 at s = 0 and normal to step 2 at s = sin(?/2). We expand the numerator in (4.64), substituting the definitions of ?2,?2, and A(s) wherever possible. After grouping terms according to the highest power of s they contain as a factor, we reduce our problem to calculating integrals of the form Ik(?2) := 2 integraldisplay 1 0 skBk(s) [A(s)s2 +?2r(s)2]3/2 ds? 1?s2, (4.65) where Bk(s? 0) negationslash= 0 and A(s? 0) negationslash= 0 and r(s? 0) negationslash= 0. We alert the reader that (4.65), unlike (4.26), in general cannot be evaluated in terms of elliptic integrals, especially since we don?t know much about the behavior 61 of A(s), r(s), or Bk(s). However, the Mellin transform with respect to ?2 does not need to know about these functions, since the main dependence on ?2 is already accounted for in (4.65). We apply (4.28) to (4.65), which yields ?Ik(?) = 2 integraldisplay ? 0 d(?2)(?2)?? integraldisplay 1 0 skBk(s) [A(s)s2 +?2r(s)2]3/2 ds? 1?s2. (4.66) We apply Fubini?s theorem to (4.66) in order to switch the order of integration. The integral with respect to?2 can be evaluated in terms of Gamma functions, which entails ?Ik(?) = 2?(1??)??(32 +?) ?(52) integraldisplay 1 0 ds? 1?s2s k?2??3Bk(s)A(s) ???3/2 [r(s)2]??+1 . (4.67) This integral converges for {? ?C : Re(?k+ 2? + 3) < 1} = {? ?C : Re? 0 is a constant quantifying the strength of step interactions, ?i corre- sponds to the radial distance in polar coordinates, mi is the discrete step density, and ? is a shape factor; note that ?(?i,?i) is a constant, independent of the distance ?i [66]. The important feature of (5.10) is the inverse-square dependence on step separation, which follows from the elastic interaction formulas of Chapter 4. As we saw in Chapter 4, the calculation of a shape factor even in the special case of axisymmetry requires a substantial amount of analysis. We alert the reader that 70 (5.10) is an approximation, omitting the ? dependence of the prefactor to avoid unwarranted complications in the subsequent coarse-graining. An alternative definition of the step chemical potential invokes the total energy, Esteps, of a step train in 2+1 dimensions. Each step contributes a line tension energy in proportion to its total arclength, because it takes energy to create the line defect of a step on a crystal surface. Moreover, the repulsive interaction between a step and its same-sign neighbors also adds to the total energy Esteps. Recall from Chapter 3 that the existence of a step introduces a mechanical strain field in the crystal. This strain field can be modeled as the result of a force dipole located at each point along the step. Each pair of neighboring step segments contributes an elastic interaction energy given by the energy of two interacting force dipoles, for example the dipoles dsie?|i, dsi+1e?|i+1. In addition, neighboring steps exert an entropic repulsion on each other. We consolidate these contributions into a single integral over the ith step, and then sum over the step train to obtain Esteps = summationdisplay i integraldisplay Li (? +Vi,i+1)ds. (5.11) To determine the step chemical potential ?i, we need to know the number of additional atoms required at the ith step to effect a given change in step train energy. The normal displacement vi(?)?t at each position ? along the ith step will require the addition of enough atoms to cover an area ofintegraltextLi vi(?)?tds. The typical cross-sectional area of an atom is ?/a. According to the definition of step chemical potential, these additional atoms at the ith step will change the total energy by an amount 1?/aintegraltextLi vi(?)?i(?)ds?t. Summing up the contributions from every step to 71 determine the change in total energy, ?E = ?Esteps?t, yields a weak formulation of the step chemical potential, ?i. We equate the sum over all steps with the direct calculation of ?Esteps using (5.11). The result is ?Esteps = 1 ?/a summationdisplay i integraldisplay Li vi?i(?)ds. (5.12) We will use this weak formulation of the step chemical potential when deriving the macroscale limit in Section 5.4.2. Before giving the final equations connecting the chemical potential to the step velocities, we take a moment to identify which of the preceding relations will carry through unchanged when we consider extensions of the terrace adatom kinetics. The definitions of local coordinates and metric coefficients, as noted above, will persist in their present form. In addition, the definitions of Ceqi and ?i are independent of the adatom kinetics on terraces, so their formulation here will carry through unaltered to the study of anisotropic terrace diffusion (Chapter 6) and electromigration (Chapter 8). Lastly, we present the step velocity law in 2+1 dimensions. By including diffusion of atoms along the step edge with constant edge diffusivityDed, the normal velocity of the ith step edge is [82] vi = e? ? dridt = ?a(Jadi?1,? ?Jadi,?) +a?s parenleftbigg Ded?s ?ik BT parenrightbigg , (5.13) where ?s is the space derivative along a step edge; ?s = ??1? ??. The first term in (5.13) is the contribution of terrace adatom fluxes. The second term is due to step edge diffusion and stems from the variation of the step chemical potential, ?i. We reason that the same ?i has to be used both in edge diffusion and in Ceqi by 72 virtue of the fact that ?i controls the equilibrium shape of a step. This equilibrium shape should be the same regardless of which kinetic pathway (edge diffusion or attachment-detachment) is operative in the deformation of the step shape. So, if mass exchange with the terrace is turned off and relaxation occurs via edge diffusion, the step attains the same shape as in the case where edge diffusion is turned off and relaxation is allowed only by attachment-detachment kinetics. This property implies that the thermodynamic driving force has to be the same chemical potential, ?i, in both cases [54]. Equations (5.2)?(5.13) in principle constitute a system of coupled differential equations for the step positions. As a discrete scheme of step flow, this system has been implemented numerically only in the case of straight or axisymmetric step trains. In this chapter we focus on (2+1)-dimensional settings where Ded = 0 (5.14) and Dad is a scalar constant. 5.3 Approximations for slowly varying step train We see from (5.13) that the adatom flux Jadi plays a pivotal role in connecting the step normal velocity vi to the step chemical potential. By solving approxi- mately the diffusion equation (5.2) following [66], we write an explicit formula for the adatom flux. We review the assumptions underlying the approximate solution of (5.2) sub- ject to (5.4),(5.5). We assume that the boundary data Ceqi varies more rapidly with 73 the step index i than the longitudinal coordinate ?. The slow variation in ? should also be expected of the solution for adatom density. We treat the derivative ?? as O(epsilon1) in comparison to the derivative ??, which is treated as O(1); epsilon1lessmuch 1. A possible geometric interpretation of epsilon1 is the ratio a/R where R = O(?) is a typical radius of curvature for steps and ? is a suitable macroscopic length [66]. Once the macroscale surface flux is derived, the assumptions for the ? and ? derivatives are relaxed: both derivatives are allowed to be O(1). By neglecting the? derivatives in (5.2), for constant Dad the diffusion equation for Ci reduces to ?? parenleftbigg? ? ????Ci parenrightbigg ? 0 , (5.15) which has the explicit solution Ci ?Ai(?,t) integraldisplay ? ?i ?? ?? d? prime +Bi(?,t) ?i t) depend on A(0) = A0 and H. The constant B3 is defined by B3 := DadCsg3?2k BT . ?2 lessmuch (1 +q?)?1 lessmuch 1 (1 +q?)?1 lessmuch?2 < 1 Step interaction A0e?CB3t A0(1 +cB3A0t)?1 Line tension A0radicalbig1?t/t? A0(1?t/t?) ing these two possibilities with the kinetic-geometric conditions on the mobility, Margetis [67] obtained the possible scaling laws listed in Table 5.1. We do not address analytically the spatial part of the approximate separable solution, H(r), which solves a nonlinear PDE. Numerical solutions of the PDE, such as those given in Chapter 10, remain our most viable option for visualizing the spatial dependence of a separable solution. 5.6 Presence of facets in the variational formulation The surface evolution law derived here conforms more closely to the behav- ior of steps than the PDEs obtained phenomenologically or by analogy with 1D macroscale models. However, the behavior of steps near a facet edge, where the sur- face energy has a singularity, is not respected by the PDE derived in this chapter. The variational formulation imposes its own natural boundary conditions at the facet edge, which in general are not identical to the boundary conditions stemming from 84 a microscale drop in height when the top step collapses [38]. We outline the natural boundary conditions in Chapter 9, in order to interpret the results of simulating the variational formulation in terms of strong PDE solutions informed by conditions at the free boundary. Only in the limiting case g1 ? 0 can we expect the variational formulation to respect the flow of steps near |?h| = 0. A hybrid approach, which informs the variational equations of the discrete height drop when the topmost step collapses, is beyond our scope in 2+1 dimensions, where an implementation of the step flow model is lacking. Comparison of the variational and hybrid approaches in the axisymmetric case [65] suggests that the discrepancy in boundary conditions is not likely to affect the observed scaling laws. To be cautious, however, we focus mainly on the case g1 = 0 in our subsequent numerical simulations. 85 Chapter 6 Macroscale equation with terrace anisotropy and step edge diffusion1 In this chapter we extend the theory of chapter 5 to cases with a tensor- valued terrace diffusivity Dad and a nonzero edge diffusivity Ded. These extensions offer a more realistic description of diffusion processes on terraces and steps, while remaining firmly based in the framework of continuum steps. Again we seek a PDE for the surface height. The main contribution of this chapter is the surface mobility M, which extends (5.24) by allowing for off-diagonal elements even in the local coordinate representation. For convenience we rewrite the quasi-steady equation (5.2) describing adatom diffusion on the ith terrace: div(Dad ??Ci) = 0 . (6.1) Equation (6.1) is complemented by the mixed boundary conditions dictated by linear kinetics: ?Jadi (?i,t) = ku[Ci(?i,t)?Ceqi (t)] , (6.2) Jadi (?i+1,t) = kd[Ci(?i+1,t)?Ceqi+1(t)] . (6.3) In this chapter the terrace diffusivity Dad is assumed to have the tensor form Dad = D11e?e? +D12e?e? +D21e?e? +D22e?e?. For the sake of some generality, 1Material in this chapter appeared previously in Quah and Margetis, J. Phys. A: Math. Theor. 41, 235004/1?18. 86 we do not impose the symmetry condition D12 = D21, although in most physical settings this equality holds. The surface flux Jadi depends on the adatom density Ci through the linear relation ? ?? ? Jadi,? Jadi,bardbl ? ?? ?= ? ? ?? ? D11 D12 D21 D22 ? ?? ?? ? ?? ? ??1? ??Ci ??1? ??Ci ? ?? ? ?i 0). The height profile H(r) 99 Table 6.1: Decay laws for height amplitude A(t) in ADL kinetics. Leftmost column indicates plausible conditions. Next two columns list decay laws for line tension and step interaction dominated ??. The time constant ? depends on A(0) and H. Line tension Step interaction |Dad|?D11D22 bgreatermuch max{(D22/D11)?2,D22/D11,(D22/D11)??2} A0radicalbig1?t/? A0 exp(?t/?) blessmuch min{D22/D11)?2,D22/D11,(D22/D11)??2} A0(1?t/?) A0/(1 +t/?) |Dad|/(ka) lessmuchaDed/(b?Cs) A0 exp(?t/?) A0/radicalbig1 +t/? solves a nonlinear PDE of the form CH ? div braceleftBigg (?xH)2 |?H|3 ? ?? ? mxx mxy myx myy ? ?? ??? bracketleftbigdivparenleftbig|?H|?Hparenrightbigbracketrightbig bracerightBigg . (6.43) The solution for A(t) is given in terms of the separation constant C and the initial amplitude A0: A(t) = A0e?Ct. Using a similar procedure, we derive other possible scaling laws for ADL kinetics under different restrictions. Our results are summa- rized in Table 6.1. We do not address the issue of solving (6.43) in this analysis. Particularly in- teresting is the case with facets. The macroscale limit breaks down near facet edges, and associated boundary conditions for H must take into account the discrete step flow equations [65]. A numerical scheme to implement these boundary conditions within a macroscopic simulation is a worthy subject of future research. A similar analysis can be carried out if terrace diffusion is the slowest process, i.e., q|?h| = |?h|D11/(ka) lessmuch 1. Then, b is approximately a constant, b ? 1. The 100 dominant terms in the mobility tensor scale as A0 or A1. Thus, we obtain ?A??Ap for p ? {0,1,2,3}, which yields four of the five decay laws already found for ADL kinetics. 6.7 Conclusion By interpreting a (2+1)-dimensional step flow model for a relaxing surface as a discretization of a macroscale evolution equation, we derived the relevant PDE for the surface height profile. The starting point is a step velocity law that accounts for anisotropic adatom diffusion on terraces, diffusion of atoms along step edges and atom attachment-detachment at steps. In the macroscale limit we obtained a relation between the surface flux and the step chemical potential. This relation involves a tensor surface mobility with nonzero off-diagonal elements even in the local coordinate representation. Combining the step velocity law with the constitutive relation between the sur- face flux and the step chemical potential resulted in a fourth-order nonlinear PDE for the surface height. Transforming the mobility tensor from local step coordinates to fixed coordinates introduced a dependence on the height partial derivatives. This dependence offers a plausible scenario of how an epitaxial surface can exhibit dif- ferent decay laws. We found separable solutions for the height that approximately satisfy the evolution equation under certain conditions. These separable solutions exhibit different decay and may be used as a guide in interpreting experimental observations from a macroscale perspective. 101 Our PDE only accounts for a part of the possible microscopic physics. In Chapters 7?8 we discuss the macroscale consequences of additional effects, such as step permeability, adatom drift under the influence of an electric field, and desorp- tion of adatoms. The calculations of this chapter serve as a prototypical example for later macroscale limits. 102 Chapter 7 Step permeability and macroscale limit 1 In this chapter, we extend the model of continuum steps in 2+1 dimensions to include the possibility of permeable steps. Also known as step transparency, this phenomenological mechanism allows adatoms to jump from one terrace to a neighboring terrace without first attaching and detaching at the intervening step. In principle, step permeability can lead to long-range interactions, since an adatom might cross many terraces before attaching to a step and influencing the step motion. If such an effect is realized in the equations of step flow, we might naturally wonder whether the macroscopic limit also shows evidence of long-distance interactions. The main achievement of this chapter is the derivation of a macroscale limit consistent with permeable steps in 2+1 dimensions. The macroscale PDE we derive here leaves unchanged the local character of step interactions, finding only a renormalized kinetic parameter as the result of including step permeability. This renormalized kinetic parameter generalizes to (2+1)-dimensions a formula in [88], which can be interpreted via an electric-circuit analog of the paths by which adatoms can jump between terraces. In this chapter we continue to assume that mass transport on terraces is due to diffusion only. We omit from the diffusion equation terms corresponding to de- 1Material in this chapter appeared previously in Quah, Young and Margetis, Phys. Rev. E 78, 042602/1?4. 103 position flux, electromigration current, and desorption of adatoms. The latter two effects will be addressed in Chapter 8. Here we discuss only the changed boundary conditions in the case of permeable steps. 7.1 Step permeability According to the original BCF framework for studying step motion, mass transport of adatoms between adjacent terraces requires the intermediate processes of attaching and detaching at the common step edge. By decomposing the mass transport mechanism in this way and allowing for an Ehrlich-Schwoebel barrier, we have the freedom to choose two kinetic rates ku and kd, which in principle can explain such phenomena as the appearance of different island sizes during multilayer growth [23]. Motivated more by intellectual curiosity than by discrepancies between theory and experiment, Ozdemir and Zangwill imagined the possibility of direct exchange of adatoms and vacancies from one terrace to another with kinetic rate p [81]. This exchange mechanism is called step permeability or step transparency. Pierre-Louis used this theoretical construct in his attempt to resolve the question of temperature-dependent conditions for electromigration-induced step bunching [88]. The effect of step permeability is to change the linear kinetic boundary con- ditions (3.4) and (3.5) to ?Ji,? = ku(Ci ?Ceqi ) +p(Ci ?Ci?1), ? = ?i, (7.1) Ji,? = kd(Ci ?Ceqi+1) +p(Ci ?Ci+1), ? = ?i+1, (7.2) where p is the permeability rate and Ji,? is the adatom flux component in the step- 104 normal direction; see Chapter 5 to review the geometry and notation. The same p appears in both boundary conditions by Onsager reciprocity relations [89]. 7.2 Macroscale limit In order to derive a macroscale limit from (7.1),(7.2), and the terrace diffusion equation, we impose the same geometric conditions as in chapter 6, along with the conditions on the kinetic rates given in the preceding sections. The terrace widths are assumed to be much smaller than [66]: (i) the length ? over which the step density varies; and (ii) the step radius of curvature, if applicable. The discrete step density a/??i with ??i = ??|?i??i is assumed to remain fixed while we formally let ??i ? 0. Recall that ? is a dimensionless coordinate with discrete isolines corresponding to steps. The restriction on the kinetic rates is that D/(?a) = O(1) for ? = ku,kd,p. We used a similar assumption in Chapter 5 to derive our first macroscale limit. Here the condition D/(?a) = O(1) plays the same role, allowing us to keep fixed the dimensionless parameters D/(?a) as various limits are taken. If we wanted to proceed in the same vein as our first macroscale limit in Chap- ter 5, we would begin by solving the diffusion equation ?2Ci = 0 on the ith terrace, subject to the boundary conditions (7.1),(7.2). Immediately we notice that step per- meability explicitly couples the resulting adatom concentrationCi with that of other terraces. For a step train with a finite number of terraces, the system of boundary conditions can be solved simultaneously by inverting a circulant matrix [43]. This 105 coupling at the step level suggests that the macroscopic limit of Ci must also be incorporated into the macroscale equation. A similar approach was taken in [16], although the authors provide only a brief discussion of the connection to the per- meability rate p. A direct solution of the terrace diffusion equation, in any case, would yield two integration constants that have to be found by applying the boundary conditions. We can bypass some tedious algebra by letting the values of adatom concentration and its normal derivative play the role of integration constants, for which the bound- ary conditions will serve as a linear system of equations. Accordingly, we assume the existence of the interpolating functions C and J for adatom concentration and transverse flux, respectively, which agree with Ci and Ji,? when evaluated at the ith step edge. We take C and J as our primary variables. Whether we evaluate the terrace concentration and its derivative at the upper step edge (?i) or the lower step edge (?i+1) is irrelevant in the macroscale limit, since ?i and ?i+1 approach the same value, corresponding to the transverse coordinate of a single level set. Our final goal is a system of equations for C and J?, which emerge as the macroscale limits of C and J, respectively. 106 7.2.1 Straight steps In the straight step setting with no Ehrlich-Schwoebel barrier, the boundary conditions (7.1),(7.2) reduce to ?Ji = k(Ci ?Ceqi ) +p(Ci ?Ci?1), x = xi (7.3) Ji = k(Ci ?Ceqi+1) +p(Ci ?Ci+1), x = xi+1. (7.4) We now identify Ci,Ci+1,Ci?1 with values of a continuous adatom concentra- tion function at the left step edge of the corresponding terrace. In particular, we define the interpolant C by Ci|xi =: C(i). As usual, the subscript xi to the right of the vertical bar denotes the point of evaluation. Similarly, the flux Ji|xi at the left step edge is identified with the value J(i) of an interpolating flux function. If a value of the terrace adatom concentration Ci at a point other than the left step edge is desired, then an approximation in terms of the concentration C(i) and flux J(i) can be found by applying a Taylor expansion along with Fick?s law: Ci|x = C(i)?D?1J(i)(x?xi) +O[(x?xi)2]. (7.5) For the purposes of the macroscale limit, only the first-order term in the expansion (7.5) needs to be retained. Using (7.5) atx = xi+1 and a similar expansion for Ci?1|xi, the boundary conditions now read ?J(i) = k[C(i)?Ceqi ] +p[C(i)?C(i?1) +D?1(xi ?xi?1)J(i?1)] (7.6) J(i) = k[C(i)?D?1(xi+1 ?xi)J(i)?Ceqi+1] +p[C(i)?D?1(xi+1 ?xi)J(i)?C(i+ 1)]. (7.7) 107 In the absence of an Ehrlich-Schwoebel barrier, direct subtraction of (7.6) and (7.7) leads immediately to ?2J(i) = k[Ceqi+1 ?Ceqi +D?1(xi+1 ?xi)J(i)] +p[C(i+ 1)?C(i?1) +D?1(xi+1 ?xi)J(i) +D?1(xi ?xi?1)J(i?1)]. (7.8) Bringing all the flux terms to the left side, we find J(i) bracketleftbigg 2 + kD(xi+1 ?xi) + pD(xi+1 ?xi?1) bracketrightbigg + pD(xi ?xi?1)J(i?1) = ?k[Ceqi+1 ?Ceqi ]?p[C(i+ 1)?C(i?1)]. (7.9) Assuming thatC and the macroscale equilibrium concentrationCeq both admit Taylor series expansions about xi, we can replace the differences Ceqi+1 ?Ceqi and C(i+ 1)?C(i?1) according to Ceqi+1 ?Ceqi = (xi+1 ?xi) ?C eq ?x vextendsinglevextendsingle vextendsinglevextendsingle xi +o(xi+1 ?xi), (7.10) C(i+ 1)?C(i?1) = (xi+1 ?xi?1) ?C?x vextendsinglevextendsingle vextendsinglevextendsingle xi +o(xi+1 ?xi?1). (7.11) Since the macroscale height profile is required to be differentiable, the mean- value theorem allows us to substitute 2a/|?xh|forxi+1?xi?1 anda/|?xh|forxi+1?xi, where the evaluation point of the x-derivative is left unspecified for now. The constitutive relation (7.9) then becomes J(i) bracketleftbigg 2 + kaD|? xh| + 2paD|? xh| bracketrightbigg + paD|? xh| J(i?1) = ? ka|? xh| ?xCeq? 2pa|? xh| ?xC. (7.12) In the limit as a/?? 0, the positions of steps i and i?1 are assumed to tend to a common point x, so the equation for the macroscale flux J is J(x) bracketleftbigg 2 + (k+ 2p)aD|? xh| bracketrightbigg = ka|? xh| ?xCeq + 2pa|? xh| ?xC. (7.13) 108 To eliminate the explicit dependence on C, we take the macroscale limit of the boundary condition (7.6) directly and then substitute for J(x) using (7.13). The macroscale limit of (7.6) reads ?J(x) = k(C?Ceq) +p bracketleftbigg ?xC a|? xh| + J(x)aD|? xh| bracketrightbigg . (7.14) Combined with the macroscale flux from (7.13), we find an equation for C ? Ceq. k(C?Ceq) + pa|? xh| ?xC = ka |?xh|?xC eq + 2pa |?xh|?xC 2 + (k+2p)aD|?xh| bracketleftbigg 1 + paD|? xh| bracketrightbigg parenleftbigg 1 + (k+ 2p)a2D|? xh| parenrightbigg (C?Ceq) = a2|? xh| parenleftbigg 1 + paD|? xh| parenrightbigg ?xCeq + p 2a2 kD|?xh|2?xC. (7.15) In principle, we regard (7.15) as a series in powers of a, which holds identically asa/?? 0. By equating the coefficients of like powers ofa, we obtain conditions on the coefficients. The dominant balance of terms in (7.15) indicates that we can take C = Ceq in the macroscopic limit. The factor multiplying C?Ceq on the left-hand side of (7.15) remains O(1) as a/? ? 0, while the right-hand side vanishes since p/k = O(1) does not scale with a. Replacing C by Ceq in (7.13), we find J(x) = ? D?xC1 + 2D (k+2p)a|?xh| = ?DCsk BT ?x? 1 +qeff|?xh|, (7.16) whereqeff := 2D/(keffa)andkeff := k+2p. Thus, the effective attachment-detachment rate keff appears as an additive renormalization of k to k+ 2p. We pause to mention a few previous derivations of this result in the literature. The unifying theme of these derivations is a replacement of the step-terrace system by a simpler mass transport mechanism with renormalized parameters. In this vein, 109 Figure 7.1: Electric circuit analogous to a terrace and a permeable step, adapted from [88]. Nozi`eres proceeds by analogy with an electric circuit [78], where the attachment and detachment processes are regarded as a series of resistors, and the effect of step transparency is to provide a parallel resistor through which an adatom current can flow across the step. The simplest illustration of this analogy appears in the subsequent paper by O. Pierre-Louis [88], adapted here as Figure 7.1. A more complicated electric circuit analogy, incorporating several steps and terraces, is given by Jeong and Williams [43]. 7.2.2 Circular steps with axisymmetry We use the familiar setting of axisymmetry to introduce the additional algebra needed to deal with an Ehrlich-Schwoebel barrier. This bias of adatom flux toward a preferred direction appears in the boundary conditions as an inequality in the kinetic rates: kd > ku if downhill drift is preferred, or ku > kd if uphill drift is 110 preferred. The boundary conditions (7.1),(7.2) read ?Ji = ku(Ci ?Ceqi ) +p(Ci ?Ci?1), r = ri (7.17) Ji = kd(Ci ?Ceqi+1) +p(Ci ?Ci+1), r = ri+1. (7.18) Again we identify Ci,Ci+1,Ci?1 with values of a continuous interpolating function C at the inner step edge of the corresponding terrace. This convention says C(i) := Ci|ri. Similarly, we define an interpolating flux function J(i) by J(i) := Ji|ri. An approximation of the adatom concentration at the outer step edge can be found in terms of the concentration C(i) and flux J(i) by applying Fick?s law: Ci|ri+1 = C(i)? J(i)D (ri+1 ?ri) +o(ri+1 ?ri). (7.19) Retaining only the terms of the expansion (7.19) that are linear in ri+1 ?ri, the boundary conditions (7.17),(7.18) become ?J(i) = ku[C(i)?Ceqi ] (7.20) +p[C(i)?C(i?1) + J(i?1)D (ri ?ri?1)], (7.21) J(i) +?rJi|ri(ri+1 ?ri) = kd[C(i)? J(i)D (ri+1 ?ri)?Ceqi+1] +p bracketleftbigg C(i)? J(i)D (ri+1 ?ri)?C(i+ 1) bracketrightbigg . (7.22) Eventually we plan to take a macroscale limit from these equations, which allows us to consider a/ri lessmuch 1. The nondimensional parameter kda/D, however, remains O(1) in the macroscale limit. Hence we have the inequality a ri lessmuch kda D . (7.23) 111 Multiplying (7.23) by |J(i)/a| and identifying the left-hand side as vextendsinglevextendsingle vextendsingle?rJi|ri vextendsinglevextendsingle vextendsingle by use of the diffusion equation, we find vextendsinglevextendsingle vextendsingle?rJi|ri vextendsinglevextendsingle vextendsinglelessmuch vextendsinglevextendsingle vextendsinglevextendsinglekdJ(i) D vextendsinglevextendsingle vextendsinglevextendsingle. (7.24) Therefore it is justified to neglect ?rJi|ri from (7.22), since our final goal is a macroscale limit relating the variables J and C. These macroscale variables sat- isfy C(i)?C(i?1) = (ri ?ri?1)?rC|ri +o(ri ?ri?1), (7.25) C(i+ 1)?C(i) = (ri+1 ?ri)?rC|ri +o(ri+1 ?ri), (7.26) J(i)?J(i?1) = o(1). (7.27) The discrete step density a/(ri+1 ?ri) is assumed to approach the macroscale surface slope |?h|. In fact, replacement ofri+1?ri (orri?ri?1) bya/|?rh|evaluated somewhere on [ri,ri+1] (or [ri?1,ri]) is justified by the mean value theorem and the assumed smoothness of the macroscale surface height h. The desired system of equations for the macroscale variables is ?J parenleftbigg 1 + pD a|? rh| parenrightbigg = ku(C?Ceq) + pa|? rh| ?rC, (7.28) J parenleftbigg 1 + p+kdD a|? rh| parenrightbigg = kd(C?Ceq)? kda|? rh| ?rCeq ? pa|? rh| ?rC, (7.29) where Ceq is the macroscale limit of Ceqi . To solve for the radial flux J, we divide (7.28) by ku and (7.29) by kd; the difference of the two resulting equations retains only the concentration gradients, not their pointwise values. We obtain the formula J = ?D(2p/k)?rC +?rC eq 1 + 2p/k+q|?rh| , (7.30) 112 where k := 2(k?1u +k?1d )?1 and q := 2D/(ka); see section 5.4.1. Unlike [16], we do not treat the macroscopic limit of adatom concentration as a variable in its own right. Indeed, we have seen above in the case of straight steps that the macroscale limit renders such a treatment unnecessary. A similar argument in the axisymmetric setting allows us to replace C by Ceq in (7.30). With this substitution, (7.30) becomes the desired relation between flux and chemical potential: J = ?D ?rC eq 1 +qeff|?rh| = ? DCs kBT ?r? 1 +qeff|?rh|, (7.31) where qeff = D/(keffa) and keff = k + 2p. Evidently, the effect of step permeability in the axisymmetric setting is precisely the same renormalization of the kinetic pa- rameter that we saw in the straight step case. The additional mass transport mech- anism of adatoms hopping directly between terraces is reflected in the macroscale equations by a modified attachment-detachment coefficient. Under the conventional interpretation of step permeability, where p > 0, the new attachment-detachment coefficient is greater than k. However, no physical principle is violated if p < 0. This generalization has been suggested as a means of introducing different kinetic regimes as the temperature or the supersaturation is varied [90]. We discuss below the implications of negative p in relation to a two-region model of step dynamics by Weeks et al. [125]. 113 7.2.3 Steps in 2+1 dimensions The geometry considered here is identical to that of Chapter 5, where? is used as the transverse coordinate and ? is used as the longitudinal coordinate. Starting with the boundary conditions (7.1),(7.2), we replace the adatom concentrations and the normal flux components with either the macroscale interpolants, e.g., C(i), J?(i), or series expansions about the upper step edge of the appropriate terrace. For example, the domain of the adatom concentration Ci?1 is ?i?1 D) leads to a negative coefficient of m0. Similarly, a permeability coefficient p sat- isfying p < ?k/2 leads to a negative q in the macroscale limit of this chapter. Evidently, the two-region model and the picture of permeable sharp steps both have enough flexibility in the choice of parameters to match a range of coefficients in the slope-dependent mobility. The agreement between these two nanoscale models at the level of macroscale equations is noteworthy, in view of the strikingly differ- ent pictures they present at the nanoscale. The persistence of a slope-dependent mobility with (possibly negative) material parameter q indicates just how robust 121 the macroscale limit is, emerging in the same form out of two conceptually distinct nanoscale models. Our interpretation of permeability, as an additive renormalization of the ki- netic rate in the macroscale equations, is simple enough to allow meaningful compar- isons between predictions of the macroscale PDE (chapter 10) and Monte-Carlo sim- ulations of stepped surfaces. An alternative interpretation by O. Pierre-Louis [89], in which the equilibrium concentration of the step flow equations is renormalized through a suitable combination of the boundary conditions, also admits direct com- parison with atomistic simulations. The question of whether step transparency itself is exhibited by realistic surfaces cannot be addressed without a further understand- ing of the physical mechanisms which combine with permeability to produce different morphological changes. Our discussion of electromigration in Chapter 8 helps com- plete the picture of possible nanoscale physics, in order to categorize the kinetic regimes under which a stepped surface can evolve. By studying those regimes where permeability would couple with electromigration to yield a qualitatively different evolution, we might in principle be able to detect the presence and strength of step transparency for realistic materials. 122 Chapter 8 Electromigration and desorption1 In this chapter, we consider another set of extensions to the BCF picture of terrace adatom diffusion. First we introduce an electric field, which biases the motion of adatoms in a prescribed direction. Next we allow desorption of adatoms into the vapor phase. Combined with our results for anisotropic terrace diffusion and step edge diffusion, the inclusion of electromigration and desorption provides a comprehensive model of the mass transport processes over terraces and steps. Our analysis aims to complement previous works in surface electromigration; see, e.g., [12, 15, 27, 28, 35, 52, 53, 55, 59, 60, 61, 72, 74, 88, 90, 92, 101, 102, 103, 111, 112, 124, 125]. These focus on the development of instabilities (e.g. step meandering and bunching) across scales, often from a dynamical system viewpoint. We, on the other hand, place electromigration in the context of connecting discrete schemes for interacting steps to global PDE laws and variational principles in 2+1 dimensions [16, 64, 66, 78, 80, 96, 98, 106, 109, 123]. The PDE we derive can be studied using local analysis, to see the possible control that an electric field might have on the spatial dependence of the slope profile near a facet. Investigating the possible global control of the surface evolution by an electric field is within the scope of our numerical method, as we discuss in Chapter 10. 1Parts of this chapter appeared previously in Quah and Margetis, Macroscopic Electromigration on Stepped Surfaces, submitted to SIAM Multiscale Modeling and Simulation. 123 The study of electromigration began with the observation by Latyshev et al. [59] that direct current heating of stepped Si(111) resulted in a step bunching instability. The appearance of step bunches depends both on the direction of the applied current and on the surface temperature. An early model by Stoyanov [111] offered a possible explanation for the link between electromigration current and step bunching. Stoyanov?s model works only up to 850?C, above which temperature the direction of electric current required for step bunching undergoes two distinct rever- sals. In particular, direct current in the step-down direction produces step bunching at temperatures around 935?C and 1275?C, but at 1190?C step bunching requires current in the step-up direction. More recent proposals to address the relation between electromigration and step bunching (although not directly relevant to our focus here) include but are not limited to the work by Pierre-Louis and M?etois [88, 90], and by Zhao, Weeks, and Kandel [124, 125]. Experimental studies found the bias in adatom motion?the drift velocity?to be always parallel to the imposed electric field E. This drift velocity v is related to E by [15, 28]: v = D(Z ?e)E T , (8.1) where D is the adatom diffusivity, Z?e is the effective adatom charge,2 and T is the Boltzmann energy3. At the level of steps, the electric field amounts to the addition of a convective (linear in the density gradient) term in the diffusion equation for the 2Here, we follow the notation of [15]. The symbol Z? (with an asterisk) should not be confused with the usual symbol for the conjugate of a complex number. 3In this chapter only we use units where kB = 1, thus replacing by T the kBT of other chapters. 124 adatom density; see (8.35). Regarding (8.1), |Z?| is larger than unity for metals but can be much smaller than unity for semiconductors [15]. We note in passing that E influences the temperature in an experimental setup, since it directly controls the electric current that flows through and eventually heats up the sample. In addition to electromigration, we will allow adatoms to desorb into the vapor with characteristic time ?. So, the term C/? is added in the adatom diffusion equa- tion where C is the adatom density. We study if and how this addition at the BCF level influences the macroscopic limit. Note that a phenomenological approach to desorption in crystal morphological evolution is offered by Villain [119]. We derive ?-dependent corrections of the large-scale flux, and describe how large ? needs to be so that desorption can be neglected in the macroscopic laws. 8.1 Macroscale limit and assumptions The macroscopic limit aims to describe a continuous surface at sufficiently large scales. Previous studies of electromigration and desorption apparently have not focused on global evolution laws for the surface height, and hence carry a perspective different from ours. These works (i) focus mainly on stability issues, (ii) make use primarily of nanoscale models where steps are everywhere parallel, and (iii) allow steps to interact mainly through terrace diffusion (with less emphasis on step-step entropic or elastic dipole interactions). Our derivation uses the multiscale expansion reviewed in Chapter 5, which handles a fully (2+1)-dimensional geometry just as easily as a straight-step geom- 125 etry. The expansion parameter is the step height, a. This formulation relies on the separation of local variables for interacting steps into fast and slow, and can encompass additional physical effects such as desorption and electromigration. As a result, we derive a fully continuum model in 2+1 dimensions that extends the model of Chapter 5 to a greater variety of experimentally relevant situations. The addi- tional effect of anisotropic terrace diffusion can also be incorporated subsequently, following the calculation of Chapter 6. We recall the geometric assumptions underlying the macroscopic limit. The (microscale) terrace width is assumed to be of the order of the step height, a, and small compared to: (i) the step radius of curvature; (ii) the length over which the step curvature varies; and (iii) the macroscopic length over which the step density varies. A step train satisfying these conditions is referred to as ?slowly varying?, cf. Sections 5.3 and 6.1. The macroscopic limit is reached formally as a/? ? 0 where ? is a macroscopic length. The step density is fixed and approaches the surface slope. We assume that the initial (at t = 0) ordering of steps is preserved by the flow (for t> 0),4 and that step trains are monotonic (say, with descending steps in the direction of increasing transverse coordinate). Furthermore, we require that certain kinetic parameter groups involving the step height, a, remain O(1) (as a/? ? 0). In particular, we take D/(ka) = O(1) where k is any relevant kinetic rate for attachment-detachment of adatoms at a step edge. 4This property should be derivable from microscale physics, such as the entropic repulsions that enforce the non-crossing condition. This aspect is not studied here. 126 Figure 8.1: Schematic of 1D steps (cross section) with an applied electric field, E. Steps are descending with increasing x; E is shown to be in the step-up direction (E < 0); and ku (kd) is the kinetic rate for atom attachment-detachment from a terrace to an up- (down-) step edge. 8.2 One-dimensional step models In this section, we review briefly ingredients of step motion, including repul- sive entropic and elastic dipole step interactions, for 1D geometries: straight and concentric circular steps. Related versions of the adatom diffusion equation with either an electric field or desorption are solved explicitly. The step configuration is shown in Figure 8.1. Both the diffusion of adatoms (density gradient) and electric field (drift velocity) contribute to the flux in each terrace. Each step advances or retreats in response to the net (normal) flux from the neighboring terraces. The surface evolves as the steps move, lowering its free energy. 127 8.2.1 Straight steps The position of the ith step is denoted by xi(t) (t is time), where xi+1 > xi for all t, and Ci(x,t) is the adatom density on the ith terrace, i.e., the region xi < x < xi+1. A constant electric field is applied externally in the x direction (perpendicular to steps), introducing a drift velocity v according to (8.1). The diffusion equation satisfied by Ci(x,t) is D?2xCi ?v?xCi ???1Ci = ?tCi ? 0 xi ri. The external electric field E is taken to be radial, viz., E = erE (er: radial unit vector) where E may vary with the polar coordinate, r. To simplify the exposition, we set E = E(r)/r.6 The electric field produces the radial drift velocity v = erv = er V/r according to (8.1), where V = DZ?eE/T. The radial adatom flux, Ji(r) = er ?Ji(r), is Ji(r) = ?D?rCi +vCi ri (1/2)|v|. Let Ci be a C2 (twice continuously differentiable) function on Ui and C1 (continuously differentiable) function on Ui, the closure of Ui. Then, for ? = O(1), Ci(?,?) = C(0)i (?,?) +o(1) where C(0)i (?,?) = Bi(?) +Ai(?) integraldisplay ? ?i dz ??|z? ?|z exp bracketleftbiggintegraldisplay z ?i (v???)|?prime D d? prime bracketrightbigg ; (8.40) o(1) ? 0 as epsilon1? 0. The integration constants Ai,Bi are given by Ai(?) = (1 +v?|?i/ku)C eq i+1 ?(1?v?|?i+1/kd)C eq i D ku??|?i parenleftbigg 1? v?|?i+1k d parenrightbigg + parenleftbigg 1 + v?|?ik u parenrightbiggbracketleftbigg D kd????fi|?i+1 + parenleftbigg 1? v?k d parenrightbigg fi|?i+1 bracketrightbigg , (8.41) Bi(?) = parenleftbigg 1 + v?|?ik u parenrightbigg?1bracketleftbigg Ceqi (?) + Dk u??|?i Ai(?) bracketrightbigg , (8.42) where fi(?,?) is defined by fi(?,?) := integraldisplay ? ?i ??|z ??|z exp bracketleftbiggintegraldisplay z ?i (v???)|?prime D d? prime bracketrightbigg dz (?,?) ?Ui . (8.43) In the above, Q|z denotes the value Q(? = z). The role of parameter epsilon1 is to express small variations of the step edge curvature, which in turn cause slow variations of Ci = Cepsilon1i(?,?) with respect to ?, permitting a perturbative treatment. In the limit epsilon1 ? 0, the step edges become 1D, approaching concentric circles or straight lines (depending on limits of ??, ??). We also assume that ku and kd are positive (as is typical in crystalline materials). For comments on the condition |v|/kl < 2 where l = u, d, see Remark 8.3.2. Proof. For convenience and notational economy, in this proof we set D = 1 (or, define the inverse length hatwidev := v/D and drop the hat) and suppress the terrace 11Alternatively, it can be assumed that v(r) varies slowly, e.g., v = v(epsilon11?,epsilon11?), epsilon11 lessmuch epsilon1. 138 index i unless noted otherwise. Assuming that a solution C(r) = Cepsilon1(r) exists and is unique, as implied below, we partly separate scales in PDE (8.39) for ? = O(1). So, we expand formally Cepsilon1(r) in an epsilon1-power series; each coefficient, C(j), is allowed to be a function of (?,?,?), which are treated as independent variables:12 Cepsilon1(?,?) = C(0)(?,?,?) + summationdisplay j?1 epsilon1j C(j)(?,?,?) = C0(?,?,?) +o(1) , (8.44) where the remainder approaches 0 as epsilon1? 0, assuming continuity of the solution with epsilon1.13 In addition, the operator ?? is replaced by the linear combination ?? ??? +epsilon1?? . By dominant balance of O(epsilon10) terms, (8.39) entails the (zeroth-order) PDE braceleftbigg 1 ???? bracketleftbigg ?? parenleftbigg? ? ???? parenrightbigg +?? parenleftbigg? ? ???? parenrightbiggbracketrightbigg ?v? ??? ? ?vbardbl ??? ? bracerightbigg C(0) = 0 (r ?U) , (8.45) which must be solved under conditions (8.37), (8.38) for C(0) (for epsilon1-independent ku, kd). Next, we show that, for given continuous Ceq(?), the above boundary value problem forC(0) has at most one solution. We apply a standard energy method [24]; the same argument carries through for proving uniqueness of Cepsilon1. First, by the transformation C(0) = ?Ce(1/2) R r ri v?dr prime, PDE (8.45) is converted to the Helmholtz equation ?r ?C ? 14(|v|2 ?2divv) ?C = 0 (r ? U: terrace), where ?C obeys the Robin 12The use of the extra slow variable ? = epsilon1?, although formally justifiable, is not deemed necessary. 13This physical property is invoked in conjunction with the assumption that eliminating the curvature variation (as epsilon1 ? 0) results in a well-defined adatom density, as suggested by Section 8.2. 139 boundary condition ?????C =: ??? ?C = K(r) ?C?k(r)Ceqe?(1/2) R r ri v?dr prime (r ??U); ? is the unit outward normal vector, K(r) := k(r)?(1/2)? ?v, and k(r) = ku for an up-step edge (? = ?i where ? = ?e?) while k(r) = kd for a down-step edge (? = ?i+1 where ? = e?). Consider sufficiently small |v| so that K(r) > 0, yet |v|2 ? 2divv, consistent with the neglect of the term (divv)Cepsilon1 in the diffusion equation. Now suppose there exist two solutions, say ?C1 and ?C2, of the boundary problem for C(0), with ? := ?C1 ? ?C2; this ? satisfies the given Helmholtz equation with boundary condition ???? = K?. Second, define the non-negative energy E[?] = 12 integraltextU |??|2 +pi12?2 where pi12 := (1/4)(|v|2 ? 2divv). By Green?s identity, E[?] =integraltext?U(????)? = ?integraltext?U K?2 ? 0; thus, ?? 0 which entails C1 ?C2. Based on this uniqueness assertion, we construct solution (8.40)?(8.42). The ?-dependence of Ceq in the boundary data suggests that we look for a C(0)(r) that depends on ? and ? (but not ?). Such a C(0)(r), if it can be constructed plausibly, is interpreted as the leading-order, unique solution of the boundary value problem. Hence, we solve (8.45) by dropping the ?-derivatives. Two successive integra- tions with respect to the variable ? immediately yield C(0)(?,?) = B(?)+ ?A(?) integraldisplay ? ?i exp braceleftbigg ? integraldisplay z ?i bracketleftbigg? ?prime ?? ??prime parenleftbigg? ? ??prime parenrightbigg ?(v?|?prime)??primeD bracketrightbigg d?prime bracerightbigg dz, (8.46) where ?A and B are integration constants. Equation (8.46) readily reduces to (8.40) by direct integration (in ?prime) of the first term in the exponent and the subsequent substitution A(?) := (??/??)|?i ?A(?). The coefficients A(?) and B(?) are determined through boundary conditions 140 (8.37), (8.38) with ? = ?prime. Accordingly, we obtain the system? ?? ?? ? Dk u??|?i 1 + v?|?ik u D kd????f|?i+1 + parenleftbig1?v ?/kd parenrightbigf| ?i+1 1? v?|?i+1 kd ? ?? ?? ? ?? ? A B ? ?? ?= ? ?? ? Ceqi Ceqi+1 ? ?? ? , where f(?,?) is defined by (8.43). Solving this system leads to (8.41) and (8.42). Henceforth, we drop the dependence of v? = v?e? (and vbardbl = e??v) on i, since v is considered macroscopic and e?, e? vary slowly with ?. It can be verified that (8.40) reduces to the solutions for 1D settings of sections 8.2.1.1 and 8.2.2.1. Remark 8.3.2. Thus far, we invoked the condition |v|/klscript < 2 (lscript = u, d). By definition (8.1) for v, this restriction amounts to imposing |v| 2klscript = D 2klscripta |Z?e||E|a T < 1 , lscript = u, d . In typical experimental situations, D/(2klscripta) is of the order of 102 or smaller [43], |eE|a/T is of the order of 10?5, and for semiconductors |Z?| ranges from 10?4 (or even smaller values) to about 10 [15, 28]. Hence, |v|/klscript would not exceed values of the order of 10?2. From the perspective of perturbation theory adopted here, D/(klscripta) is treated as an O(1) quantity whereas |v|/klscript will be considered as o(1) (see Proposition 8.4.1) when a ? 0. This view furnishes the distinct physical contribution of drift in the large-scale surface flux. Practically, setting |v|lessmuchklscript seems to more closely describe semiconductor surfaces, where |Z?| can be much smaller than unity [28]. Mathemat- ically, this view is a choice for obtaining the distinguished-limit contribution of v as a? 0. 141 In Section 8.4.1 we impose |v|/kl lessmuch 1; as a result, Fick?s law for the large-scale flux does not distinguish between up- and down-step edges, being symmetric in ku and kd. 8.4 Evolution equations at the macroscale In this section we find the evolution equations for the macroscale height in the presence of an electric field. The modification of Fick?s law (8.36) at the nanoscale leads to a convective term in the constitutive relation between the macroscale flux and chemical potential, as we show in Proposition 8.4.1. Unchanged from Chapters 5,6, and 7 are the mass conservation law and the formula for the macroscale chemical potential. The law of mass conservation, in the absence of edge atom diffusion and ma- terial deposition from above, reads ?th+ ?divJ = 0 ; (8.47) see (5.44). The macroscale chemical potential ? is regarded as a smooth interpolation, through a suitable Taylor expansion, of the microscale step chemical potential ?i. The presence of an electric field is not expected to change the surface energy due to steps, which experience the same line tensions and repulsive interactions as in the case of zero electric field. We adopt the formula from Chapter 5, which reads ? = ??div braceleftbigg ?m[m(g1 +g3m2)] ?h|?h| bracerightbigg , m := |?h| . (8.48) 142 8.4.1 Macroscopic Fick?s law with drift In this section, we derive the constitutive relation between macroscale flux and chemical potential in full 2D. The restrictions to 1D geometries then follow as special cases. We choose units where the macroscopic length ? is unity (? = 1) for convenience. Proposition 8.4.1. Suppose that in the macroscopic limit, a ? 0, the step density mi = a/(????i) and the kinetic parameters D/(kla) are O(1) while v/kl = o(1), where l = u, d and ??i := ?i+1 ??i. Then, the solution to diffusion equation (8.39), under boundary conditions (8.37), (8.38), gives rise to the macroscale constitutive relation J = ? ?? ? J? Jbardbl ? ?? ?= ?CsM? bracketleftbigg ??? TD v parenleftBig 1 + ?T parenrightBigbracketrightbigg parenleftbigg v = DZ ?eE T parenrightbigg , (8.49) where the (E-independent) mobility M (with units of length 2/energy/time) is a second-rank tensor. In the coordinate system (?,?), this M has the familiar repre- sentation of Chapter 5, namely, M = DT 11 +q|?h| ? ?? ? 1 0 0 1 +q|?h| ? ?? ? , q := 2D ka , k := 2 k?1u +k?1d . (8.50) Proof. The starting point is the solutionC(0)i derived in Proposition 8.3.1; see (8.40). We drop the superscript (denoting perturbation order) in C(0)i for ease of notation. The microscale adatom flux components are obtained by the formulas Ji,bardbl = ?D??1? ??Ci +vbardblCi , Ji,? = ?D??1? ??Ci +v?Ci . 143 The plan is to consider the restrictions of these components to ? = ?i, and view these restrictions as interpolations of (continuous) smooth functions, J?(r,t) and Jbardbl(r,t). First, we compute the requisite derivatives of Ci by (8.40): ??Ci ???Bi +??Ai integraldisplay ? ?i ??|z ??|z exp bracketleftbiggintegraldisplay z ?i ??prime v?|?prime D d? prime bracketrightbigg dz , ??Ci = Ai ??? ? exp bracketleftbiggintegraldisplay ? ?i ??prime(v?|?prime) D d? prime bracketrightbigg . It follows that Ji,bardbl = ?D??1? ??Bi +vbardblBi , (8.51) Ji,? = ?D??1? Ai +v?Bi ? = ?i . (8.52) Consider (8.41), (8.42) (in Proposition 8.3.1) for Ai, Bi. In the limit a? 0, or ??i ? 0 with????i = O(a), these formulas simplify via the expansionintegraltext?i+1?i F(?)d? = F(?i)??i +o(??i), where F(?) is any continuous function. After some algebra and neglect of o(??i) terms, (8.51) and (8.52) become Ji,bardbl|?i ? vbardbl bracketleftbigg D ????i parenleftbiggCeq i kd + Ceqi+1 ku parenrightbigg +Ceqi bracketrightbigg ?D bracketleftbigg ?bardblCeqi + D? ???i parenleftbigg? bardblC eq i kd + ?bardblCeqi+1 ku parenrightbiggbracketrightbigg D ????i parenleftbigg 1 ku + 1 kd parenrightbigg + parenleftbigg 1 + v?k u parenrightbigg , (8.53) Ji,?|?i ? D ????i(C eq i ?C eq i+1) +v?C eq i D ????i parenleftbigg 1 ku + 1 kd parenrightbigg + parenleftbigg 1 + v?k u parenrightbigg as ??i ? 0 ; ?bardbl := ??1? ?? . (8.54) In the above, all variables are evaluated at (the same) ? along the ith step edge. We seek further simplification of (8.53) and (8.54). In the macroscopic limit, we invoke the (assumed as well defined) C1 function Ceq(r) where Ceq(r)|?i ?Ceqi , 144 Ceqi+1 ? Ceq(r)|?i + (??Ceq)|?i??i +o(??i) and ?bardblCeq|?i ? ?bardblCeqi = ?bardblCeqi+1 +O(??i). We keep only those combinations of microscopic parameters that remain O(1). For example, the step density mi = a/(????i) approaches the positive surface slope, i.e., lima?0mi = |?h|. On the other hand, the ratio v?/ku, involved in the denominator of Ji,bardbl and Ji,?, is treated as negligibly small (compared to unity) by virtue of our hypothesis; see also Remark 8.3.2. Without further ado, we make the substitutions Ceqi+1 ?Ceqi = (??Ceq)????i + o(??i) and Ceq ? Cs(1 +?/T) (|?| lessmuch T), where ?? := ??1? ??. The resulting limits for the flux components read lim a?0 Ji,bardbl|?i =: Jbardbl(r)|?i = ?CsDT ?bardbl?(r) +Csvbardbl bracketleftbigg 1 + ?(r)T bracketrightbigg , (8.55) lim a?0 Ji,?|?i =: J?(r)|?i = ?CsDT ???? TDv? parenleftBig 1 + ?T parenrightBig 1 +q|?h| ? = ?i . (8.56) These relations are identified with (8.49) under definition (8.50) for the mobility M. It can be verified directly that the same limits emerge if the evaluation point of fluxes is at ? = ?i+1 [66]. Two remarks on the results of Proposition 8.4.1 are in order. Remark 8.4.2. In the absence of electromigration (v = 0), we recover the relation J = ?CsM ? ?? of Chapter 5. For nonzero drift (v negationslash= 0), the additional term derived above is affine with ?. This new relation is recast to the zero-drift form via an exponential transformation of ?; see Section 8.4.3. Remark 8.4.3. An alternate proof of Proposition 8.4.1 makes direct use of the 145 adatom concentration, Ci, and the normal flux component, Ji,?, avoiding entirely the use of integration constants Ai and Bi. This approach treats boundary conditions (8.37), (8.38) for ?prime = ? as a system of equations for Ci(?i) and Ji,?(?i). This argu- ment is also applicable to situations where the boundary conditions couple densities of different terraces, e.g., in the case with step permeability (see Chapter 7). By Proposition 8.4.1, we state the following corollary for cases of symmetry. Corollary 8.4.4. Consider 1D settings with v negationslash= 0, where translational or rotational symmetry causes all dependent variables to have zero?-derivatives. The macroscopic surface flux corresponding to (8.49) is J = J(?)e? , J(?) = ?CsDT ???1 +q|? ?h| + Csv1 +q|? ?h| parenleftBig 1 + ?T parenrightBig , (8.57) where ? = x for straight steps and ? = r for concentric circular steps; q = 2D/(ka). Recall that, in the microscale model underlying the limit of this section, the contribution (divv)Ci is left out from the terrace diffusion equation (see Section 8.2.2). For a discussion on a correction to the macroscopic limit, see Remark 8.5.4. 8.4.2 Evolution equation in Cartesian system In this section, we describe the PDE for the surface height by combining ingredients (8.47) and (8.48)?(8.50) and making use of Cartesian coordinates, (x,y). First, we recall the non-singular orthogonal matrix S, Eq. (5.26). S(?xh,?yh) = (e? e?) = 1|?h| ? ?? ? ??xh ?yh ??yh ??xh ? ?? ? (?hnegationslash= 0) . (8.58) 146 The mobility tensor M, which is defined by (8.49) in the (?,?) coordinate system, is now expressed in the Cartesian coordinate system (x,y) by M(x,y) = SM(?,?)ST ; (8.59) as usual, ST denotes the transpose of S (ST = S?1). The Cartesian components of the surface flux (8.49) read (elscript: orthonormal vectors, lscript = x, y): Jx =? Cs1 +q|?h| braceleftbiggD T bracketleftbiggparenleftbigg 1 +q (?yh) 2 |?h| parenrightbigg ?x??q (?xh)(?yh)|?h| ?y? bracketrightbigg ? parenleftbigg 1 + ?T parenrightbiggbracketleftbiggparenleftbigg 1 +q (?yh) 2 |?h| parenrightbigg vx ?q (?xh)(?yh)|?h| vy bracketrightbiggbracerightbigg , (8.60) Jy =? Cs1 +q|?h| braceleftbiggD T bracketleftbiggparenleftbigg 1 +q(?xh) 2 |?h| parenrightbigg ?y??q(?xh)(?yh)|?h| ?x? bracketrightbigg ? parenleftbigg 1 + ?T parenrightbiggbracketleftbiggparenleftbigg 1 +q (?xh) 2 |?h| parenrightbigg vy ?q (?yh)(?xh)|?h| vx bracketrightbiggbracerightbigg ; vlscript = elscript ?v (lscript = x, y) . (8.61) The PDE for the surface height follows from the mass conservation statement. Using the Cartesian representation of (8.47), we have ?th = ??(?xJx +?yJy) . (8.62) This relation leads to a nonlinear, fourth-order parabolic PDE for h after substitu- tion for Jx,Jy, and ? from (8.60), (8.61), and (8.48). 8.4.3 Change of variables We show that law (8.49) is recast to the form J = ?cM??? (c = const.) , (8.63) 147 leading to the PDE ??t?h = div?[M????] , (8.64) with suitable choices of the tensor M and variable ?[??], the transformed chemical potential; the dimensionless variables ?t, ?h, ?? and operators div?, ?? are defined in Appendix B. Constitutive relation (8.49) for the surface flux reads J = ?CsD? x tildewiderM?[?????u(1 + ??)] , (8.65) where the dimensionless mobilitytildewiderM ? in the (x,y) representation ? and drift velocity u are defined by tildewiderM(x,y) := tildewideS? ? ?? ?? 1 1 + ?q|???h| 0 0 1 ? ?? ???tildewideST , u := ?xD v ; ?q := qh0? x . (8.66) The corresponding change-of-basis matrix, tildewideS, is tildewideS ? S(??x?h,???y?h); cf. (8.58). Consequently, the PDE for the nondimensional height, ?h, reads ??t?h = div?{tildewiderM?[?????u(1 + ??)]} . (8.67) Next, we show (8.63) and (8.64). By inspection, we start with the transfor- mation ? = (1 + ??)f? , (8.68) where the (nonzero) f? = f?(?x,?y) is to be determined. By virtue of M???? = (Mf?)?[????+ (??f?/f?)(1 + ??)] and (8.65), we have the consistency relations Mf? = tildewiderM , (??f?)/f? = ?u , 148 the second one of which yields f?(?x,?y) = Ae?u?(?x,?y/?) ? M = A?1tildewiderM(???h)eu?(?x,?y/?) , A = const.negationslash= 0 . (8.69) Relations (8.63), (8.64) ensue. For definiteness, take A = 1. The transformed chemical potential ? and mobility M are ?(?x,?y) = braceleftbigg 1? ??0 div? bracketleftbigg ?g ?? ?h |???h| +|?? ?h|???h bracketrightbiggbracerightbigg e?(ux?x+uy?y/?) , (8.70) M(x,y) = eux?x+uy?y/?tildewideS?tildewiderM?tildewideST = eux?x+uy?y/?tildewideS? ? ?? ?? 1 1 + ?q|???h| 0 0 1 ? ?? ???tildewideST . (8.71) The finite element scheme introduced in Chapter 2 can easily accommodate the transformed chemical potential and mobility. The change-of-variables is not deemed likely to introduce numerical instability if the grid Peclet number Pe = |v|?x/D is much less than 1 [77].14 For easier readability of the code, however, we opted to include the convective term directly when simulating the macroscopic effect of an electric field. A thorough comparison of the two methods is still lacking. 8.5 Extensions of macroscopic limit In an attempt to formulate a reasonably general theory of macroscopic surface relaxation, we enrich the BCF model with additional microscale effects. These are: (i) the exponential law Ceqi = Cs exp(?i/T) in the place of its linearization, and (ii) atom desorption. Our broader goal with these extensions is to reconcile the macroscopic theory with realistic situations where an electric field is present. 14Pe compares the relative importance of convection to diffusion; ?x : mesh size. 149 We show that (i) can modify significantly the constitutive relation between surface flux and chemical potential. In contrast, effect (ii) arguably has vanishingly small influence on the macroscopic evolution law. 8.5.1 Exponential law for step chemical potential To derive the constitutive relation between large-scale flux and chemical po- tential, we approximated the difference of equilibrium concentrations Ceqi by use of the linearized form (8.6). The quantitative justification for this approximation is not clear in the literature, particularly since the chemical potential ?i is not measured directly. However, the PDE is easily modified to accommodate the complete, expo- nential law. For a similar modification in 1+1 dimensions withv = 0 and long-range step interactions, see [123]. By skipping details irrelevant to the modification at hand, we start from (8.54). By setting ?i = ?(r)|?i we expand the difference Ceqi+1 ?Ceqi as follows: Ceqi+1 ?Ceqi = Cs bracketleftbigg exp parenleftbigg? i+1 T parenrightbigg ?exp parenleftbigg? i T parenrightbiggbracketrightbigg = aCs|?h| exp parenleftbigg?(??) T parenrightbigg? ??(??) T , ?? ? [?i,?i+1] . Hence, the normal flux component Ji,? at (?i,?) becomes Ji,?|?i = ? D? ???i aCs |?h| exp parenleftbigg?(??) T parenrightbigg? ??(??) T +v?Cs exp parenleftbigg?(? i) T parenrightbigg D ????i parenleftbigg 1 ku + 1 kd parenrightbigg + parenleftbigg 1 + v?k u parenrightbigg . (8.72) Similarly, the ?-derivative of Ceqi appearing in (8.53) for the longitudinal flux com- 150 ponent, Ji,bardbl, now acquires exponential factors when expressed in terms of ?i : Ji,bardbl|i = bracketleftbigg ? D? ???i parenleftbigg 1 ku + 1 kd parenrightbigg ?1 bracketrightbigg?1braceleftbigg ?CsD? ???i vbardbl bracketleftbigg 1 kd exp parenleftBig?i T parenrightBig + 1k u exp parenleftBig?i+1 T parenrightBigbracketrightbigg ?vbardblCs exp parenleftBig?i T parenrightBig +CsDexp parenleftBig?i T parenrightBig?bardbl?i T + D 2Cs kd????i exp parenleftBig?i T parenrightBig?bardbl?i T + D2Cs ku????i exp parenleftBig?i+1 T parenrightBig?bardbl?i+1 T bracerightbigg . (8.73) The coarse-graining procedure carries through as in Proposition 8.4.1. The components of the flux J(r) are found to satisfy the equations J?(r) parenleftbigg2 k + a D|?h| parenrightbigg = exp bracketleftbigg?(r) T bracketrightbiggbraceleftbigg ?aCs|?h|???(r)T +v? aCsD|?h| bracerightbigg , (8.74) Cs exp bracketleftbigg?(r) T bracketrightbigg? bardbl?(r) T = ? 1 D braceleftbigg Jbardbl(r)?vbardblCs exp bracketleftbigg?(r) T bracketrightbiggbracerightbigg . (8.75) The last two equations result in the effective constitutive relation J = ?Cse?/TM? parenleftbigg ??? TD v parenrightbigg . (8.76) By inspection of (8.76), we have the following remark: Remark 8.5.1. The modified constitutive relation for the macroscopic surface flux results from the invariant under the law Ceq = Ceq[?] form (8.64) and the definition ?(r) := e?(r)/T?v?r/D , (8.77) which is a direct extension of the linearized version, (8.68) or (8.70). 8.5.2 Adatom desorption The effect of desorption is expected to affect only the relation between the large-scale surface flux and chemical potential. In this section, we show that, under 151 certain conditions, desorption does not appear in the macroscopic laws to leading order in a. 8.5.2.1 Solution for microscale diffusion Following the rationale of slow and fast step variables of section 8.3.3, we write the (quasi-steady) diffusion equation on the ith terrace with desorption time ? and a drift velocity v (v? = e? ?v) as ?? parenleftbigg? ? ????Ci parenrightbigg ? v???D ??Ci ? ?????DCi (?i 0, and yields a sufficiently smooth Ci. Since the kernel is proportional to ??1, the procedure yields a power series in 1/?, whose convergence rapidity is thus enhanced by increasing ?. The value of the Born-Neumann scheme becomes evident for ? ? ?i lessmuch 1. Since integraltext??i F(?prime)d?prime = F(?i)(???i)+o(???i) for any continuous F, the associated iterated integrals contribute respective ascending powers of ? ??i. So, evidently, in the limit maxi??i = maxi{?i+1 ??i} ? 0, constructing a solution to (8.79) via a Born-Neumann series corresponds to producing a Taylor series in ???i for Ci(?,?). This result is directly applicable to the macroscopic limit of the step system. 8.5.2.2 Limit a? 0 Next, we derive a macroscopic Fick?s law with desorption and external electric field. We follow the main procedure of Proposition 8.4.1 by use of formula (8.79) (in Proposition 8.5.2) for Ci. Suppressing the ? dependence, we define Cvi (?) := Bi +Aifi(?) , (8.83) which is the solution form without desorption.16 Thus, Ci satisfies Ci(?)?Cvi (?) = (?D)?1 integraldisplay ? ?i ?2?prime fi(?)?fi(? prime) ??primefi Ci(? prime) d?prime . (8.84) By differentiation of (8.84) with respect to ? we obtain the normal flux com- ponent, Ji,? = e? ?(?D?Ci +vCi): Ji,??Jvi,? = ???fi? integraldisplay ? ?i ?2?prime Ci(?prime) ??primefi d? prime + v? ?D integraldisplay ? ?i ?2?prime fi(?)?fi(? prime) ??primefi Ci(? prime) d?prime , (8.85) 16Ai and Bi entering Cvi in principle depend on ? via boundary conditions at step edges. 154 where Jvi,? := ?D??Cvi +v?Cvi ; recall that ?? = ??1? ??. Now consider the limit ?? := ???i ? 0. By Taylor expanding the right hand sides of (8.84) and (8.85), we readily obtain Ci(?)?Cvi (?) = (2D?)?1Cvi (?i)?2?i ??2 +O(??3) , (8.86) Ji,?(?)?Jvi,?(?) = ???1Cvi (?i)??i ??+O(??2) as ?? ? 0 . (8.87) Next, we apply conditions (8.37) and (8.38) for?prime = ?. With recourse toCvi (?i+1) = Cvi (?i) + (??Cvi )|?i??i +o(??i), and likewise for Ji,?(?i+1), we find ?k?1u Jvi,? = Cvi ?Ceqi , (8.88) k?1d bracketleftbiggparenleftbigg 1 + kdD?wi parenrightbigg Jvi,? + (??Jvi,?)?wi bracketrightbigg ? bracketleftbigg 1 + parenleftbiggv ? D + 1 kd? + ?wi 2D? parenrightbigg ?wi bracketrightbigg Cvi ?Ceqi+1 , (8.89) where ?wi := ??i ??i, and Cvi , Jvi,? and ??Jvi,? are evaluated at ?i. Note that (8.88) does not involve ? explicitly. In contrast, (8.89) manifests desorption terms in the right hand side. If these, ??1-scaled, terms can be neglected appropriately, the resulting system becomes identical to that without desorption; the corresponding terms in Ci and Ji,? can then be dropped. Remark 8.5.3. By inspection of the (attachment-detachment) boundary conditions (8.88) and (8.89) for the adatom flux and density, sufficient conditions for neglecting the desorption effect in the macroscopic limit [assuming k/klscript = O(1), lscript = u, d] are |v| D greatermuch 1 k? and a k? lessmuch 1 . (8.90) 155 Then, the large-scale adatom flux is not affected by ? uniformly in space coordinates. Both of these conditions are expected to be met in a wide range of physical situations. In particular, the first condition amounts to having |v|? greatermucha when D/(ka) = O(1). From the viewpoint of coarse graining and perturbation theory, formulas (8.90) sim- ply state that the effect of desorption is of higher order (in a). We conclude this section with the following observation. Remark 8.5.4. The (thus far neglected) term (divv)Ci in the terrace diffusion equation for the adatom density Ci can be treated on the same footing as the des- orption term Ci/?. Hence, in assessing the validity of dropping the former, it is reasonable to repeat the above procedure with ??1 being replaced by divv, in an appropriate sense. 8.6 Conclusions Starting with the BCF model of step flow, we derived macroscopic evolu- tion laws for the surface height in settings with an electric field. We considered microscopic processes of isotropic adatom diffusion on terraces, attachment and de- tachment of atoms at step edges with an Ehrlich-Schwoebel barrier, and desorption of atoms to the surrounding vapor. Energetic effects such as entropic and elastic dipole step interactions were included. Our central contribution is Fick?s law (8.49) relating surface flux J, chemical potential ?, and drift velocity v. The linear combination of ? and ?? of this law is consolidated into a single variable, ?, by an exponential transformation. Accord- 156 ingly, the PDE for the macroscale height is recast to form (8.64), which has the same structure as the evolution law ?th = ?Cs div{M???} derived in the absence of an electric field. We showed that desorption is a higher-order effect in a sense dictated by multiscale expansions in the step height, a. Many effects are absent from our analysis. For instance, time-dependent, coupled electromagnetic fields are not accounted for. In the same vein, the control of evolving surface morphologies by electric fields was barely touched upon. Our assumption of a slowly-varying step train restricts the validity of the PDE to regions where the steps do not experience drastic instabilities. We also leave out long-range step interactions, e.g., interactions mediated by bulk stress. We postpone until Chapter 9 a more thorough discussion of the appropriate boundary conditions for the derived PDE, and the effect of an electric field on the height profile near this boundary. 157 Chapter 9 Facets and boundary layers In this chapter, we discuss analytical aspects of faceted relaxation of crystal surfaces. As we see in the next chapter, the results of applying the finite element method to the previously derived PDEs have been validated only in the limiting case g1/g3 lessmuch 1. For larger values of g1/g3, our numerical scheme is challenged by the appearance of facets, which expand during the surface relaxation. The variational approach on which our numerics are based still offers a possible connection with the physically realistic case of nonzero line tension. To see this connection, we turn to analytical methods, focusing on the local behavior of the height profile near a facet edge. We study how the slope far from the facet approaches continuously its boundary value (|?h| = 0 at the facet edge), for different values of the material parameters and the electric field. This local analysis of this chapter complements the global analysis we used to derive decay rates for approximately separable PDE solutions, which are assumed to exist based on experimental observations and 1D step simulations. Our analysis is based on a conjecture (or ansatz) about the sep- arability of space and time dependence of the slope profile, which entails a locally factorizable solution of the PDE. 158 9.1 Faceting on crystal surfaces The presence of facets indicates clearly that a given surface orientation is below roughening, which is the context we assume when deriving a PDE for h from the motion of steps. However, strictly speaking the domain of this PDE does not include the facet itself, where the surface energy is singular due to the vanishing slope |?h|.1 We then face the problem of determining the appropriate boundary conditions at the facet edge. Analysis of the PDE and the boundary conditions near the facet edge provides insight into the local behavior of the slope profile. This analytically predicted behavior could be experimentally testable [115]. 9.2 Formulation of faceted relaxation as a free-boundary problem Spohn [109] first studied facet evolution as a free-boundary problem for nonlin- ear PDEs. Here, we summarize the main idea and extensions. The moving boundary is the facet edge, a curve we denote by rf(t) ? L0(t). More precisely, for a train of descending steps in the vicinity of a facet with height hf (Fig. 9.1), we identify the facet with the set {r = (x,y) : h(x,y,t) = hf(t)}, and rf(t) ?L0(t) is the boundary of this set. (Recall that Li was used in Chapter 5 to denote the ith step edge. Al- though, strictly speaking, the facet edge does not correspond to the topmost step, we use L0 for lack of a more suggestive notation.) 1The reason for this pathology is that step positions near the facet edge do not satisfy the usual assumptions of the macroscopic limit. Mathematically, the PDE can be extended on the facet by the ?subgradient (subdifferential) formulation? [24], as outlined below. 159 Figure 9.1: Orthogonal projection of a 2D step train in the vicinity of a facet. 160 The motion of the facet edge is found by imposing the requisite boundary conditions on the PDE solution. The slope profile at perpendicular distanced? from the facet edge has a (characteristic) d1/2? behavior related to a Pokrovsky-Talopov phase transition, as predicted by Jayaprakash et al. [41, 42] for equilibrium crystal shapes. Margetis, Aziz and Stone [64] applied the free-boundary approach in an axisymmetric setting (circular steps) to determine how the width of a boundary layer near the facet edge scales with the ratiog3/g1 of step interaction strength to the step line tension, when g3/g1 lessmuch 1. In the case of straight steps, Odisharia used the free- boundary approach to establish the time decay of self-similar profiles [79]. Shenoy et al. applied nonlinear Galerkin schemes to incorporate facets into a variational formulation of surface relaxation [108]. 9.3 Natural boundary conditions of the variational problem for pe- riodic profiles To investigate the effect of the free energy coefficients g1 and g3 on the ob- served slope profile near a facet, one must in principle write down the boundary conditions at the facet edge. A complete set of boundary conditions emerges ?nat- urally? by considering the variational approach in more detail. We alert the reader that not all of the natural boundary conditions play a role in our subsequent local analysis. For the sake of completeness, we present the full set of natural bound- ary conditions, so that the numerical results of Chapter 10 with g1 negationslash= 0 can be interpreted from the perspective of strong PDE solutions informed by conditions at 161 the free boundary. Not all of these boundary conditions are physically meaningful, but they possess the advantage of being enforced automatically by any numerical scheme based on the variational formulation. An ?unphysical? condition is that the chemical potential extends continuously across the facet edge. In fact, this latter condition can be replaced by a boundary condition respecting the discrete sequence of collapse times for the top step [37]. In this case, the resulting set of boundary conditions, while consistent with the microstructure of the crystal surface, unfor- tunately requires feedback from discrete step simulations to implement a numerical solution of the free-boundary problem. We now develop a brief argument for the natural boundary conditions arising from the variational formulation for evolution of a spatially periodic profile h. This approach follows closely the presentation of [79] in 1+1 dimensions. We restrict our attention to diffusion-limited kinetics, where the mobility M = CsD/(kBT) is a scalar constant. To address ADL kinetics we would need to invoke the semi- implicit time stepping of Chapter 2, which yields a weighted H?1 inner product in the extended gradient (subgradient) formulation. By using the height profile at a previous time step in the evaluation of the tensor mobility M, we determine explicitly the inner product through which the surface evolution is interpreted as a subgradient flow. This issue has been addressed in the doctoral dissertation by Odisharia [79]. The macroscale equations have the familiar structure of three key relations: (i) mass conservation, (ii) Fick?s law for surface diffusion, and (iii) the equation for chemical potential in terms of slope. We restate these three ingredients under the 162 assumption of diffusion-limited kinetics. ?h ?t = ??divJ, J = ?CsDk BT ??, ? = ?g3?E?h = ?g3? x div?N, where div? is the divergence operator with respect to nondimensional spatial vari- ables2, ?x is a characteristic length in the basal plane, and E = integraltext ?(?h)dx with ?(?h) = g1g3|?h|+13|?h|3. To emphasize the importance of the ratio g1/g3, we choose to treatE andNas nondimensional; the chemical potential acquires the correct units through the prefactor ?g3/?x. The auxiliary vector N is introduced to write the variation of E more explicitly: outside the facet the formula N = g1g3 ?h|?h| +|?h|?h holds. We take the height h to be a periodic function of the nondimensional spatial variables x,y, which entails [0,1]?[0,1]?[0,?) owner (x,y,t) mapsto?h(x,y,t) ? (?h0 ??,h0 +?); see Appendix B.3 Because the surface energy density ?(?h) has a corner singularity at ?h = 0, we need to make use of a generalized notion of differentiation when combining the three ingredients of the macroscopic theory to incorporate facets. This generaliza- tion of a gradient flow begins with the definition of subgradient for a real-valued 2Following the notation of Appendix B; in the calculations below we drop the subscript ? for the sake of readability. 3The range for h is chosen to accommodate the transient behavior in which the height slightly exceeds its initial amplitude maxh(?,?,0) = h0 at the beginning of a relaxation experiment. 163 function on R2. Here we follow the presentation of Odisharia [79] to introduce the subgradient in the context of surface relaxation. We use??,??to denote the Euclidean inner product in the following Definition 9.3.1. The subgradient of ? : R2 ?R, denoted ??, is the set ?? = {a?R2 : ?a,????(?h+?)??(?h) ?? ?R2}. To illustrate this definition, we consider first the singular case ?h = 0. Then ?? = {a?R2 : ?a,??? g1g 3 |?|+|?|3/3}, which is the ball of radius g1/g3 around 0. In the smooth case (?hnegationslash= 0), we have ?? = {a?R2 : ?a,??? g1g 3 |?h+?|? g1g 3 |?h|+ |?h+?| 3 3 ? |?h|3 3 }, which reduces to the single element a = g1g 3 ?h |?h| +|?h|?h, after expanding the norms in Taylor series about ? = 0. Evidently, the subgradient agrees with the usual notion of gradient when the free energy density is smooth enough. We note that the calculation of a subgradient depends on which inner product ??,?? is used in the definition. For the free energy density?, the natural choice for??,??is the Euclidean inner product. For a functional such as E[?h], we must use the inner product of the H?1 Sobolev space [24]. The essential property of the H?1 inner product is an identity that follows (formally) 164 from integration by parts: ?u,v?H?1 = ????1u,v?L2, (9.1) where??,??L2 denotestheusualinnerproductofsquare-integrablefunctions: ?f,g?L2 = integraltext [0,1]2 fgdA. The idea of the H ?1 inner product becomes more transparent by using the Fourier transform ?u of periodic u with zero mean. The variation of E must be expressed using the language of subgradients so that the chemical potential ? can be defined in the case of facets. In this vein, we compute the subgradient of the functional E[h] = integraltext ?(?h)dx with respect to the H?1 inner product. By definition, u??H?1E ?E[h+?]?E[h] ??u,??H?1 ??. (9.2) Having at our disposal the calculation of ??(?h) in both the smooth and the non- smooth cases, we let N(x,y) ???(?h) and consider u = ?divN. For such N we have ?(?h+??)??(?h) ??N(x,y),??? (9.3) pointwise, which can be integrated to obtain E[h+?]?E[h] ? integraldisplay ?N(x,y),???dx = ? integraldisplay ?divN,??dx (9.4) = ??divN,??L2 (9.5) = ??divN,??H?1. (9.6) Therefore u = ?divN ? ?H?1E. In fact, the subgradient of the functional E consists precisely of these elements [47]; i.e., ?H?1E = {?divN : N(x,y) ???(?h) for each (x,y)}. 165 Returning to the problem of surface evolution, with DL kinetics we have ?th = ?CsD? 2g3 kBT ?divN globally, (9.7) where N is the canonical restriction of the subgradient ?H?1E, i.e., the element ?nearest? to the origin with respect to the H?1 norm. As we observed above, for ?h negationslash= 0 the subgradient consists of a single element, and the canonical restriction does not come into play. Only in the singular case ?h = 0 do we need to invoke the canonical restriction to determine ht uniquely. The auxiliary function N ? ?H?1E can be found by solving the minimization problem min N??? bardbl?divNbardbl2H?1 ? min N??? bardbl?divNbardbl2L2. (9.8) The canonical restriction N which solves the minimization problem (9.8) must have a square integrable derivative ?divN. This condition implies that N and ? = ?g3 divN are continuous on [0,1]2. We are now in a position to state the natural boundary conditions arising from the variational formulation. The first boundary condition states that the slope |?h| approaches zero at the facet edge L0(t). Since N is restricted to the ball of radius g1/g3 on the facet, while |N| = g1g3(1 + |?h|2) outside the facet, continuity of N dictates that |?h| ? 0 as r ? rf ?L0, i.e., ?h = 0, (x,y) ?L0. (9.9) The second boundary condition in DL kinetics follows from enforcing continu- ity of surface flux (?divN in the discussion above) at the facet edge. Outside the 166 facet, we have a surface flux given by J = ?CsDk BT ?? = CsD?g3k BT ?div parenleftbiggg 1 g3 ?h |?h| +|?h|?h parenrightbigg . (9.10) The flux on the facet follows from mass conservation, ?hf +?divJf = 0, where hf(t) is the facet height. Solving for the divergence of Jf, we write divJf = ? ?hf ? on the facet. (9.11) Fix t, and integrate (9.11) over the interior of L0(t). We apply the divergence theorem to conclude contintegraldisplay L0(t) Jf ??ds = ???1?hfAf(t), (9.12) where Af(t) is the area of the facet at time t. Now we invoke continuity of J at L0(t) to replace Jf in (9.12) with its coun- terpart (9.10) restricted at L0(t). The result is ?hfAf(t) = ? contintegraldisplay L0(t) CsD?2g3 kBT ?div parenleftbiggg 1 g3 ?h |?h| +|?h|?h parenrightbigg ??ds. (9.13) The next boundary condition follows from continuity of the surface height h = h(x,y,t) at the facet edge. The height continuity dictates hf(t) = h(x,y,t)|r?L0(t). (9.14) We now differentiate with respect to t to obtain ?hf(t) = d dth(r,t) = ?h ?t vextendsinglevextendsingle vextendsinglevextendsingle r?L0(t)+ ??h|r ? ?r?t. (9.15) The value of ?th as (x,y) ? L0(t) from outside the facet is governed by the relaxation PDE for DL kinetics. ?h ?t = ? CsD?2g3 kBT ?div parenleftbiggg 1 g3 ?h |?h| +|?h|?h parenrightbigg . (9.16) 167 Using (9.16) together with?h|L0(t) = 0,we have the third boundary condition, ?hf = ?CsD?2g3 kBT ?div parenleftbiggg 1 g3 ?h |?h| +|?h|?h parenrightbiggvextendsinglevextendsingle vextendsinglevextendsingle L0(t) . (9.17) We now elaborate on the boundary conditions stemming from continuity of the chemical potential ? and the vector field N. These quantities have physical meaning outside the facet, where ?h negationslash= 0 and the step positions satisfy the usual assumptions conducive to a macroscale limit. Following Spohn [109], to determine how N and ? should be defined on the facet, we reverse the calculations that lead to surface flux outside the facet. In this vein, we start with the mass conservation law, ?hf = ??divJf, assuming that ?hf is a known function of t. Extending the constitutive relation (9.10) to the facet, we want to solve Jf = ?CsDk BT ??f for ?f. Due to the scalar proportionality constant in DL kinetics, the solution for ?f can be given as a line integral of Jf over a path contained in the facet: ?f(x,y,t) = ?0 ? kBTC sD integraldisplay (x,y) (x0,y0) Jf ?dx. (9.18) Then, continuity of ? implies ?f(rf(t),t) = ?(rf(t),t), or ?0 ? kBTC sD integraldisplay (x,y) (x0,y0) Jf ?dx = ??g3 div bracketleftbiggg 1 g3 ?h |?h| +|?h|?h bracketrightbigg , (x,y) ?L0(t). (9.19) Outside the facet we have ? = ?g3 divN, so the natural extension of N to the facet is a vector field Nf satisfying CsD?2g3 kBT ?divNf = ?hf. (9.20) 168 Continuity of N at the facet edge gives us Dirichlet boundary conditions to complement the third-order PDE (9.20): Nf = ? bracketleftbiggg 1 g3 ?h |?h| +|?h|?h bracketrightbigg , (x,y) ?L0(t). (9.21) Conditions (9.9), (9.13), (9.14), (9.19), (9.20), and (9.21) for J, ?, N, and h follow naturally from the treatment of macroscopic height evolution in terms of the variational framework. With this perspective, the smoothing of the height profile is described by a path of steepest descent for the energy with respect to the H?1 inner product. At the risk of redundancy, we emphasize that conditions (9.19), (9.20), and (9.21), for continuity of ? and N at the facet edge, rely on a nonphysical chemical potential defined on the facet.4 Implementation of the weak formulation for g1 negationslash= 0 must therefore be interpreted cautiously, since the boundary conditions it enforces are not entirely consistent with the discrete character of processes occurring on a facet [37]. 9.4 Scaling of the boundary layer width While studying numerical simulations of relaxing axisymmetric steps, Mar- getis, Aziz and Stone [64] noticed a rapid variation of the slope in the vicinity of the top facet (in contrast to the slowly varying slope farther away, in the bulk of the step train). This observation prompted them to investigate the possibility of a boundary layer near the facet edge. Intuitively, one might argue for the existence of a boundary layer by noting that a slowly varying slope profile outside the facet must 4This inconsistency has been reported previously by Israeli and Kandel [37]. 169 still connect continuously with the slope |?h| = 0 at the facet edge, for any solution that respects the boundary conditions given above. In the boundary layer, the step interaction (g3) term plays the role of a perturbation, influencing the slope profile through the balance among derivatives of different orders. This local analysis uses only the PDE for h and continuity of slope, not the other boundary conditions at the facet edge. As a result, the scaling we obtain is not altered if we replace ?-continuity by a more physical, microscale condition. To quantify the perturbation effect of the g3 term, Margetis et al. [64] introduced a small parameter epsilon1 = g3g 1 lessmuch 1. (9.22) The slope profile F := ?h? e? that solves the boundary value problem is then a function of space, time and epsilon1. The PDE for F is found by computing the directional derivative of (9.16) along the local normal: ?F ?t = CsD?2g1 kBT ? parenleftbig?div?bracketleftbige ? +epsilon1F2e? bracketrightbigparenrightbig?e ?. (9.23) Assuming a long-time similarity solution, which depends separately on t and the rescaled local coordinate, Margetis et al. found that the boundary layer width ? scales as epsilon11/3, independently of the axisymmetric initial conditions. In this section we expand the details of this calculation to accommodate ge- ometries in 2+1 dimensions. Allowing the facet to evolve in time, we assume that a boundary layer of (possibly time-dependent) width ? moves with the facet edge. Within this boundary layer, the slope profile is assumed to have an approximate factorization, which then leads to an ODE in the ?fast? spatial variable. The re- 170 sults we obtain for the boundary layer width take the form of inspired conjectures or speculations, awaiting verification by experiments and rigorous analysis. The fully 2-dimensional geometry is handled using local analysis, which sepa- rates the spatial coordinates into slow (?) and fast (?). The terminology reflects the separation of scales (short distances along the step normal direction, larger distances along the step tangential direction) over which comparable variations in the slope profile are observed. Suppose that the rapid variation of the slope F, from 0 at the facet edge ? = w(t) to its value where step interactions are less important, takes place within a boundary layer of width ?(t). We then define a rescaled transverse coordinate ? according to ? = ??w(t)?(t) , (9.24) so that?remainsO(1) inside the boundary layer. We express the spatial dependence of the slope profile in terms of this rescaled coordinate: F(?,?,t;epsilon1) = F(?,?,t;epsilon1) = F parenleftbigg??w(t) ?(t) ,?,t;epsilon1 parenrightbigg . (9.25) By substituting F for F in the PDE, we explicitly account for the rapid variation of the slope within the boundary layer. The time derivative ?tF is replaced by ?tF = ? parenleftBigg ??? ? + ?w ? parenrightBigg ??F +?tF, (9.26) while the normal derivative is replaced by ??F = ??F ? ???? = 1???F. (9.27) Motivated by analogy with observations of the step density from axisymmetric step simulations, we make the conjecture of a long-time solution that obeys the 171 following local ansatz to leading-order in epsilon1: 5 F(?,?,t;epsilon1) ?a0(?,t)?f0(?;epsilon1), |??w| 0; so the facet edge is at x = 0 as shown in Figure 9.2. We seek to characterize stationary solutions of the macroscale evolution equation. Setting ?th? 0 in the conservation law ?th+ ??xJ = 0, where 174 Figure 9.2: Cross-section of a facet and a descending train of straight steps. J = ?[?x??D?1vkBT(1 +?/(kBT))] entails the ODE ?x??l?1v ? = l?1v kBT?J0 ? ?(x) = (?0+kBT?lvJ0)ex/lv ?(kBT?lvJ0) x> 0 ; (9.34) wherelv := D/v,J0 := J(0)and?0 := ?(0). Bysubstituting? = ??g3?x(|?xh|?xh), g3 > 0, we obtain the relation ?g3(?xh)2 = (?kBT +lvJ0)x+lv(?0 +kBT ?lvJ0)parenleftbigex/lv ?1parenrightbig x> 0 ; ?g3 = ?g3 . (9.35) The above formula is simplified for 0 0 , radicalbig(?|l v|J0 ?kBT)x , xgreatermuch?lv > 0 . (9.36) 175 For xgreatermuch?lv > 0, a compatibility condition is ?DJ0 >kBT|v| . (9.37) By (9.36), changes of the electric field magnitude and direction can cause a drastic qualitative change in the slope behavior. In particular, a strong electric force Z?eE in the step-down direction causes an increase of the slope, as steps tend to bunch. (Ultimately, of course, m approaches O(?x) as x? 0.) This behavior is reversed by a strong electric force in the step-up direction, which restores the familiar O(?x) behavior. 9.5.2 Axisymmetric structure The radial case provides a model for the interplay of step edge curvature and electric field. Suppose the surface is axisymmetric, h = h(r), with the facet {(r,h) ? R+ ?R|h = const.} in 0 ? r ? rf while ?rh < 0 for r > rf; so, the facet edge is at r = rf. Consider a constant drift velocity v. The ODE for ?(r) is ?r??l?1v ? = l?1v kBT ?rfJ0/r, with solution ? = (?0 +kBT)e(r?rf)/lv ?kBT ?rfJ0 integraldisplay r rf e(r?rprime)/lv rprime dr prime r>rf , (9.38) where?0 := ?(rf),J0 := J(rf). Takingintoaccountthat??1? = g1/r+g3r?1?r[r(?rh)2], by direct integration we have ?g3(?rh)2 ? [lvrfJ0 (1?lv/rf)? ?g1 ?rfkBT](r?rf) + (kBT +?0)lvbracketleftbig(r?lv)e(r?rf)/lv ?rf +lvbracketrightbig?lv(rfJ0)(r?lv) integraldisplay r rf e(r?rprime)/lv rprime dr prime , 0 < r?rf rf lessmuch 1 ; (9.39) 176 ?gl := ?gl (l = 1, 3). Here, we consider distances r?rf from the facet boundary that are small compared to the facet radius of curvature. By analogy with Section 9.5.1, we simplify the formula for ?rh by imposing weak or strong drift: |?rh|? ?g?1/23 ? ??? ??? ? ??? ??? ? radicalBig (?0 ? ?g1/rf)(r?rf) , 0 0 , radicalbig(|l v||J0|?kBT ? ?g1/rf)(r?rf) , r?rf greatermuch?lv > 0 , (9.40) The approximation for r?rf greatermuch?lv > 0 is compatible with the condition ?DJ0 > (kBT + ?g1/rf)|v| . (9.41) Approximations analogous to the radial case can be worked out in 2D by invoking the separation of local variables into fast (?) and slow (?). Variations with respect to ? are dominant in the vicinity of the facet, which allows us to find stationary solutions by solving an ODE in ?. This approach provides formulas that directly generalize the radial case and is not further discussed here. 9.6 Conclusions In this chapter, we considered the inclusion of facets in a weak solution of the macroscale evolution equation. A weak solution with one or more facets can be interpreted as a strong solution of the PDE informed by natural boundary conditions at the moving facet edge. These natural boundary conditions include continuity of slope, continuity of flux, and continuity of chemical potential across the facet edge. To satisfy the boundary condition that the positive surface slope |?h| ap- proaches zero at the facet edge, while varying more slowly away from the facet, we 177 hypothesized the existence of a boundary layer near the edge of the facet. Under the assumption that the slope profile in this boundary layer satisfies a separability ansatz, we derived a scaling of the boundary layer width in terms of the ratio g3/g1 of step interactions to step line tension. Finally, we discussed the effect of an electric field on stationary solutions of the evolution equation in 1D step geometries. This approach yields qualitatively different slope profiles near the facet edge, depending on the strength and direction of the external electric field. 178 Chapter 10 Numerical results1 In this chapter, we illustrate numerically the rich variety of relaxation behav- iors corresponding to solutions of the macroscale PDEs derived above. Of particular interest are the simulations that demonstrate a morphological transition, from a 2- dimensional profile at the start to an essentially 1-dimensional profile at finite times. Factors affecting the onset time of this transition include (i) the material parameter q = D/ka, (ii) the aspect ratio ? = ?x/?y of the initial sinusoidal profile, where ?x and ?y are the wavelengths in x and y, respectively, and (iii) the presence of an external electric field. By varying these three factors independently, we tentatively conclude that transition phenomena can be explained as the result of longitudinal fluxes comparable in size to the transverse fluxes. We alert the reader, however, that we still lack a complete understanding of the reason for this behavior. Also, some of the results we obtain might not apply to more general periodic profiles. Initial data with additional harmonics superimposed on the fundamental sinusoidal profile might exhibit nonlinear mode coupling effects, which could destroy the separable behavior reported below. 1Material in this chapter appeared previously in Bonito, Nochetto, Quah and Margetis 2009, Phys. Rev. E 79, 050601 (R), and relies heavily on the finite element algorithms and codes gra- ciously contributed by Andrea Bonito and Ricardo Nochetto. 179 The height evolution is governed by the PDE ?h ?t = ?Cs div braceleftbigg M? bracketleftbigg ??? kBTD v parenleftbigg 1 + ?k BT parenrightbiggbracketrightbiggbracerightbigg , ? = ??div braceleftbigg g1 ?h|?h| +g3|?h|?h bracerightbigg , which we nondimensionalize prior to numerical implementation. This nondimension- alization (carried out in Appendix B) serves two main purposes: (i) to focus on the structure of the equations without the distraction of dimensional coefficients, and (ii) to isolate the relevant combinations of material parameters, in particular those that have significant temperature dependence. The second purpose furnishes five combinations of parameters to describe the experimental setting. These parameter combinations are: ? q = Dka, which compares diffusion and attachment/detachment rates; ? g1/g3, which compares the step line tension to step-step interactions; ? ? = ?x/?y, the aspect ratio of the initial data; ? h0/?x, the amplitude of the initial data; ? u = v?x/D, the nondimensional drift velocity (in the presence of an electric field). The last three parameters can in principle be adjusted independently for any given material and temperature. On physical grounds we expect some restrictions on the physically attainable profile dimensions and drift velocity magnitude. For example, a corrugated surface with aspect ratio ?x/?y = 1/10000 might be indistinguishable 180 from a one-dimensional profile, if imperfections in the patterning process ended up obscuring the long-wavelength period. Also, a sufficiently strong electric field will heat the sample, thereby changing the temperature unless the material is kept in contact with a suitably cold reservoir. The latter practical consideration is not accounted for in our simulations. The first two parameters depend strongly on the temperature and the choice of material. The parameterq is the ratio of two thermally-activated quantities with dif- ferent energetic barriers. If the energy barrier for adatom diffusion is larger (smaller) than that of attachment/detachment, then q is expected to increase (decrease) as temperature increases. The inequality between these energetic barriers need not be consistent from one material to another. In principle we can have one material where q increases with temperature, and another material where q decreases with temperature. Even two different orientations of the same material might exhibit contrasting temperature dependence of q. A good starting point for enumerating these possibilities is the table of energy barriers in [43] and the references cited therein. The ratio g1/g3 of line tension to step interactions is also sensitive to the choice of material and temperature. Roughly speaking, the line tension decreases as temperature increases, reaching g1 = 0 precisely at the roughening transition tem- perature. A more precise description of the temperature dependence of g1 is given by invoking the Ising model [73]. The strength of repulsive step interactions, quan- tified by g3, has two different contributions. Entropic repulsions become stronger as temperature increases, because steps have a greater tendency to wander and 181 collision distance is shorter. The strength of elastic dipole interactions, which are mediated through the crystal, varies with temperature mainly through the temper- ature dependence of the Young modulus for the material. Entropic repulsions tend to dominate at higher temperatures. As T approaches TR from below, we therefore expect g3 to increase. A formula for the step-step interaction coefficient g3 is given in [43]. We note that a meaningful connection with experiments requires physically attainable values for the material parameters. In particular, crystal surfaces well below the roughening transition temperature typically fail to satisfy g1/g3 lessmuch 1. To obtain a numerically reliable result, and to circumvent the difficulty of enforcing boundary conditions at the moving facet edge, we were forced to consider only g1/g3 = 0 or g1/g3 = O(10?8) in our simulations. Note that small (O(10?8)) values of g1/g3 were used in the numerical simulations in order to check the continuity of the observed phenomena withg1, as discussed below. Our preliminary investigations with g1/g3 lessmuch 1 merit careful scrutiny, in order to avoid erroneous conclusions about the relaxation of real materials. 10.1 Review of the finite element method With these limitations in mind, we now describe our numerical method, start- ing from the introductory material of Chapter 2. After nondimensionalization (see 182 Appendix B), the PDE we want to solve is ?th = divM???, (10.1) ? = ?E?h (10.2) for some appropriate mobility M and surface energy E. The spatial domain of this PDE is taken to be a rectangular region B = [0,1]?[0,1], and h,? at each time t are required to be spatially periodic. Given the initial data h(?,t) at t = 0, we want to determine the height profile at t> 0. We exploit the variational structure of (10.1),(10.2) to simulate surface relax- ation by using the FEM in space [9] and finite differences in time. The fourth-order PDE for h is conveniently written as two second-order equations: one for the height h, using mass conservation and the constitutive law for J, and another for the chem- ical potential ?, which is essentially the variation with respect to h of the energy functional E[h]. We apply a semi-implicit Euler scheme [77] to express (10.1),(10.2) as a system in the updated variables (hn+1,?n+1) where hn ?h(?,n?n) and ?n is the (adaptive) time step. The mobility M and the energy E[h] are evaluated by use of (hn,?n) to ensure linearity of the finite difference equations with (hn+1,?n+1). The equations for hn+1 and ?n+1 are recast to their ?weak form? via multi- plication by a periodic test function ? and integration by parts over B. Restricting hn+1,?n+1, and ? to the finite-dimensional space VT of continuous piecewise linear functions associated with the triangulation T , we obtain the FEM equations for (hn+1,?n+1), which are required to hold for all ??VT : 183 integraldisplay B parenleftbiggg 1 g3 ?hn+1 |?hn| +|?h n|?hn+1 parenrightbigg ??? = integraldisplay B ?n+1?, (10.3) integraldisplay B bracketleftbighn+1?+?nM?(?hn)??n+1 ???bracketrightbig=integraldisplay B hn?. (10.4) One major concern is the singularity of the g1 term and the similarity trans- formation Mx,y = SM?,?S?1, Eqs. (5.24)?(5.26), when |?h| = 0. Computing the fraction ?hn+1/|?hn| which multiplies g1/g3 is numerically unstable when |?hn| is close to zero. The elements of the change-of-basis matrix S given by (5.26) have |?h| in the denominator; consequently, the computation of S near ?h = 0 is also numerically unstable. We eliminate the singularity in these terms by adding a small regularization parameter so that the problematic denominators remain nonzero as ?h? 0; see Table 10.1 for possible regularization schemes. Mathematically it can be argued that the mobility becomes a scalar (identity tensor) wherever ?h = 0. Ideally, this limit should be respected by the regularization we choose, in order to avoid any spurious results arising from a mathematically inconsistent mobility. However, we do note that in a particular instance where one of our regularization schemes did not respect this limit, the numerical results were practically the same, as we discuss below. 10.2 Software and computers Our implementation of (10.3),(10.4) uses the FreeFem++ program, http://www.freefem.org/ff++/, developed by Olivier Pironneau, Fr?ed?eric Hecht, and Jacques Morice. The code specific to our PDE (10.3),(10.4) was written by Andrea 184 Bonito following the algorithm suggested by Ricardo Nochetto. Simulations were performed on a Dell desktop PC with 2.0GHz Intel Core 2 processor, 1GB of RAM, running Fedora Core 7 and Linux kernel 2.6.21. The following auxiliary libraries and programs were also used: ? The UMFPACK library, http://www.cise.ufl.edu/research/sparse/umfpack, was usedfordirectsolvers. UMFPACKimplementstheUnsymmetricMultiFrontal Method to solve unsymmetric sparse linear systems. ? The utility Gnuplot http://www.gnuplot.info/ was used to plot graphs. Gnu- plot offers a command-line interface for visualizing data and functions in 2D and 3D plots. ? The MATLAB software suite, http://www.mathworks.com/, was used for post- processing the FreeFem++ output files, when Gnuplot did not offer the needed capabilities for certain plots. 10.3 Zero line tension (?ideal? case) The macroscale equation we simulate for the height evolution has been shown in Chapter 5 to stem from the motion of steps, which exist as stable objects only below the surface roughening temperature. Temperatures below roughening dictate g1 > 0; otherwise there would be no free energy cost to create surface defects, and these would appear and disappear spontaneously. However, once we have found a PDE consistent with step flow, we are free to ask how its solutions behave in certain limiting cases. As a test case, we take g1 = 0 to circumvent the difficulty of treating 185 height evolution as a free boundary problem; see Chapter 9 for a discussion of the open problems in this direction. We also alert the reader that the case g1 negationslash= 0 with sufficiently high g1/g3 leads to an ill-conditioned linear system for our numerical method [6], which makes the results of a simulation unreliable. Thus far, we have not been able to overcome the numerical hurdles of larger values for g1/g3. The limit g1 ? 0 is a starting point, with the same capacity to illustrate the contributions of our model as the plausibly more physical choice with g1/g3 = O(1). In this limiting case we can see the effect of a tensor mobility and an electric field. One major concern addressed within this perspective is whether observables of the macroscopic theory (e.g., height, surface energy) are continuous solutions as functions of g1. In particular, we want to check whether the numerically observed transition for 2D to almost 1D profiles persists for nonzero (even small) g1/g3. In other words, we want to make sure that no ?spurious? singular behavior is developed by the macroscale theory as g1 ? 0. We begin this section by contrasting our tensor mobility predictions with the scalar mobility simulations of [108], both with g1 = 0. By analogy with 1D macroscale models, these simulations set M = Dk BT 1 1 +q|?h|, (10.5) i.e., a scalar mobility which does not distinguish between step-normal and step- parallel fluxes. Our purpose in conducting these tests is to see whether our finite element code generates height profiles consistent with the results of previous varia- tional methods, namely the Fourier series expansions within a Galerkin scheme used by Shenoy et al [107, 108]. 186 10.3.1 Simulations with zero E-field Our first test of the finite element implementation of the relaxation PDE used a scalar mobility, in order to validate our code against previously published results [107, 108]. In [107], Shenoy and co-authors outline their method of Galerkin expansions of the surface height profile. They start from a variational formulation of the height evolution equation and derive a system of ODEs for the coefficients in the Galerkin expansion of h. They simulate the relaxation with different values for the material parameters q and g1/g3, finding approximately separable solutions which decay either exponentially or algebraically in time when the initial height profile is a perfect sinusoid [107, 108]. We simulated a decaying bidirectional modulation, with different wavelengths in x and y. The height profiles from this simulation are plotted in Figure 10.1, with time increasing to the right. 187 Figure 10.1: Height profiles computed using scalar mobility M = Dk BT 1 1+q|?h|, for initial data h = h0 cos(k1x)cos(k2y), h0/?x = 0.03, k2/k1 = 11/24, and material parameters q = 104, g1/g3 = 0. The biperiodicity of the initial sinusoidal profile evidently persists throughout the relaxation, suggesting that the spatial dependence can be factored from the height in the approximate sense of (5.47). Assuming this factorization holds, a plot of the height peak versus time would reveal which decay law is exhibited by the solution. We plot in Figure 10.2 the evolution of energy and height peak using data from the same numerical experiment as Figure 10.1. With logarithmic axes for height and energy, the two decay curves in Fig- ure 10.2 are well-approximated by straight lines for most of the relaxation. To the eye, this trend is enough to classify the decay law as exponential. This finding is qualitatively consistent with [107, 108]. We are inclined to interpret their agreement with our plots as good evidence that Bonito?s code [7] is reliable for the parameters under consideration. However, we have not carried out comparisons of our finite element-based numerical results with Shenoy?s Fourier series-based numerics. Thus, our discussion here is restricted to qualitative (rather than precise quantitative) features of simulations. 188 Figure 10.2: Decay of height peak and energy computed using scalar mobility, for initial data h = h0 cos(k1x)cos(k2y), h0/?x = 0.03, k2/k1 = 11/24, and material parameters q = 104, g1/g3 = 0. t? and E? are measured in units where kBT?5xCsDg3?2h0 and g3h30?y?2 x are taken to be 1. 189 Table 10.1: Different regularization schemes for the tensor mobility. The first (?na??ve?) scheme was originally coded by Bonito. To respect the limit of scalar mobility as |?h|? 0, we discussed and implemented scheme 2. Nochetto suggested scheme 3, which takes a convex combination of the identity with the unregularized M. 1. M? = 1|?h|2 +?2 ? ?? ? h2x 1+q|?h| +h 2 y ? q|?h|hxhy 1+q|?h| ?q|?h|hxhy1+q|?h| h2y1+q|?h| +h2x ? ?? ? 2. M? = ? ?? ? (hx+?)2+(1+q|?h|)(hy+?)2 (1+q|?h|)(|?h|+?2?)2 ?q|?h|(hx+?)(hy+?) (1+q|?h|)(|?h|+?2?)2 ?q|?h|(hx+?)(hy+?) (1+q|?h|)(|?h|+?2?)2 (hy+?)2+(1+q|?h|)(hx+?)2 (1+q|?h|)(|?h|+?2?)2 ? ?? ? 3. M? = (1??2)I +?2M, ? = min(?,|?h|)/?. We now present the results of implementing a tensor mobility in the finite element code. Because the tensor mobility distinguishes between step-normal and step-parallel fluxes, its representation in a fixed coordinate system relies on the change-of-basis matrix S given by (5.26), which is singular at ?h = 0 as discussed above. We ran the simulations with different choices of regularization for the tensor mobility (see Table 10.1) , finding good agreement even between regularizations that gave different limits for M as ?h ? 0 [8]. This robustness with respect to regu- larization suggests that the novel phenomena reported below are insensitive to the exact behavior of fluxes around peaks and valleys. We expand on this interpretation after showing a few more plots. 190 Figure 10.3: Height profiles computed using tensor mobility, starting from the initial data h = h0 cos(k1x)cos(k2y), h0/?x = 0.03, k2/k1 = 11/24, and material parame- ters q = 104, g1/g3 = 0. Figure 10.3 shows the computed height profiles at three selected times during a relaxation with tensor mobility. The initial data is identical to that of Figure 10.1, biperiodic with different wavelengths in x and y. At later times (middle and right surface plots), the x-dependence becomes less pronounced. Even before the height peak decays to 1% of its initial amplitude, the computed height profile looks essen- tially like a 1-dimensional modulation. Because the spatial dependence of the height profile does not remain the same throughout the simulation, we suspect that a different separation ansatz h?H(x,y)A(t) is obeyed for different time intervals, e.g., h?H1(x,y)A(t), 0 0, where simulations are currently incomplete. In addition, the PDE simulated here is the limit of step motion. Thus, the limitg1 ? 0 of the macroscopic theory is essentially the limg1?0 lima?0 of the step flow equations. In principle, this limit is expected to be different from lima?0 limg1?0, in which the BCF theory is questionable or not valid. While the results forg1/g3 = 0 were still being collected and discussed, we also began preliminary trials with g1/g3 nonzero. Returning to the idea of partitioning 207 the parameter space into regions based on their eventual evolution, we wondered about the shape of the separatrix between transition-supporting and transition- suppressing regions. One cross-section that we considered trying to plot was the restriction of this separatrix by keeping fixed the aspect ratio ?x/?y, the initial am- plitude h0/?x, and the drift velocity u. The resulting cross-section would be a curve in the q?g1/g3 plane, which we hypothesized to look something like Figure 10.11. We conjectured a threshold q that increases with g1/g3, because the facet size grows more quickly for larger g1/g3, which would increase the contribution of transverse fluxes in the steeply-sloped regions between facets. The longitudinal fluxes can still be given time to effect a transition, if a larger value ofq is used to slow down the uphill and downhill mass transport. The threshold q has been determined only for the smallest values of g1/g3 in Figure 10.11. The trend of threshold q for largerg1/g3 is merely speculated. As it turned out, the calculation of these threshold q values for larger g1/g3 was never completed, due to the enormous computational time our algorithm seemed to require. After further investigation, we attributed the uncharacteristically slow convergence of the method to an ill-conditioned linear system [6]. Work to address this ill-conditioning and to extend confidently our finite element implementation in the case of large g1/g3 is still an open problem [6, 77]. The more modest goal of considering only a single nonzero value of g1/g3, small enough for the numerical method to yield reliable results, led us to take (g1/g3)(?x/h0)2 = 10?4. We simulated the relaxation using both (i) scalar and tensor mobilities, (ii) zero and nonzero drift velocities. Again the use of a tensor mobility produced a transition, which is marked by the abrupt change from pure 208 Figure 10.11: Speculation for the separatrix between transition-supporting and transition-suppressing regions, with fixed aspect ratio ? = 2/3, initial amplitude h0/?x = 0.03, and drift velocity u = 0. The upper (lower) curve would connect the computed (g1/g3,q) pairs for which a transition does (does not) occur. The uncertainty (vertical distance between the two curves) of the separatrix ordinates can in principle be decreased using bisection. In practice, a good preconditioner is needed to stabilize the numerics when g1/g3 is large [6]. 209 exponential behavior in the height peak and energy curves; see Figure 10.12. Note also that the presence of an electric field accelerates the transition. Evidently, the two major observations we made for zero g1 are also true for this particular choice of nonzero g1. It is tempting to conjecture that the transition and its acceleration by an electric field are not singular limits as g1/g3 ? 0. 10.5 Conclusions and open issues The two major contributions of our modeling in this thesis are (i) a tensor mobility, which distinguishes between step-parallel and step-normal components of the macroscale flux J, and (ii) an electric field, which leads to a convective term in the constitutive relation for J. These contributions inspired and led directly to the (hopefully) novel phenomena predicted here by finite element simulations of the macroscale evolution equation. A tensor mobility, which for ADL kinetics makes longitudinal fluxes predominant, leads to a transition from initially biperiodic data to an eventually one-dimensional profile. The inclusion of an electric field reinforces adatom fluxes along a specified direction, thereby causing the transition to be observed earlier. These two effects appear to persist continuously as g1/g3 is allowed to take nonzero values. The connection between our simulations and real materials is still elusive. The first reason is computational: at larger values of g1/g3 corresponding to real ma- terials, our finite element implementation tries to solve an ill-conditioned matrix problem. The second reason is physical: the boundary conditions enforced by ap- 210 Figure 10.12: Log plots for maximum height h?m = maxh/h0 (left axis) and sur- face energy E? (right axis), taken from [7]. These simulations used tensor mo- bility, material parameters (g1/g3)(?x/h0)2 = 10?4, q = 104, and initial data h = h0 cos(k1x)cos(k2y), k2/k1 = 11/24, h0/?x = 0.03. A transition still occurs in the presence of a facet. 211 plying the finite element method over the entire domain B are incompatible with the boundary conditions stemming from the flow of steps [65]. Even if we had trustwor- thy numerical results for larger values of g1/g3, the weak formulation itself would be questionable. We leave it for future numerical work to characterize how the predicted phe- nomena vary with increasing g1/g3. Pending issues include (i) finding a good pre- conditioner, so that the matrix problem given by our finite element equations is well-posed; (ii) enforcing more realistic boundary conditions, by restricting the do- main of the PDE and informing the solution of microscale effects at the facet edge; and (iii) developing a numerical scheme for step flow in 2+1 dimensions, in order to determine the discrete sequence of collapse times for the top step [37]. The most promising route for connecting our macroscale theory to experiments is analytical, at least until the ill-conditioning problem is resolved and a hybrid scheme is adopted to incorporate step collapse times into the weak formulation. We also hope that ongoing research in metamaterials might lead to an engineered material with the precise set of parameters that numerical constraints forced us to use in our simulations. Lacking a suitable experimental setting that can be simulated reliably by our numerics, we leave it for future research to make a more specific connection with real materials. For now we can only report the noteworthy numerical results for a hypothetical material, which appear precisely because of the novel contributions of our macroscale model; namely, a tensor mobility and the macroscale drift due to electromigration. 212 Chapter 11 Epilogue: conclusions and open questions In this thesis we have studied the macroscale consequences of different physical processes at the nanoscale, deriving PDEs to describe macroscopic surface relaxation in the absence of material deposition. Our main result is a macroscale analog of Fick?s law relating surface flux to the chemical potential. We revisit these results by discussing correspondences between the terms of this relation and the underlying nanoscale processes. For convenient reference, we give once again the macroscale analog of Fick?s law in the most general form found so far. J = ?CsM? bracketleftbigg ???kBTD?1v parenleftbigg 1 + ?k BT parenrightbiggbracketrightbigg (11.1) Energetic considerations yield the formula for the macroscale step chemical potential, which remains invariant under changes of the nanoscale kinetic processes considered in this thesis: ? = ??div braceleftbigg g1 ?h|?h| +g3|?h|?h bracerightbigg . (11.2) The inclusion or suppression of the different possible nanoscale processes is reflected in the macroscale equation (11.1) by a linear superposition of effects, or a renormalization of parameters. As an example, consider the effect of step edge diffusion. By allowing for adatoms to diffuse along the step edge, we saw in chapter 6 that the?? mobility element acquired an additional term proportional to the edge 213 Table 11.1: Macroscale consequences of several physical processes at the nanoscale. An asterisk next to ?none? alerts the reader that the macroscale effects are deemed negligible, being of higher order in the expansion parameter a. nanoscale process effect on M effect on v isotropic diffusion tensor character dis- tinguishes between transverse and longi- tudinal fluxes, with diagonal form in the local coordinate representation none anisotropic diffusion tensor form even in the local coordinate representation none step edge diffusion additional term inthe m ?? element none step transparency kinetic parame- ter q = 2D/(ka) renormalized via k mapsto?k+ 2p none electric field E none proportional to E adatom desorption none? none? diffusivity. Later, when studying the results of electromigration current, we saw that the macroscale equation acquired a convective term featuring the drift velocity v. These correspondences are summarized in Table 11.1. Although the description of step flow below roughening can vary appreciably depending on which nanoscale processes we include in the model, the macroscale limit is much more robust, taking the form of a standard conservation equation together with (11.1) and (11.2). Finding the same macroscale limit via homogeniza- tion theory [70] serves as further confirmation of the coarse-graining results shown here. 214 11.1 Summary of contributions One macroscale PDE, with appropriate modification of the parameters, en- compasses a wide range of step flow models within the BCF framework. We have studied representative versions of the macroscale evolution equation both numeri- cally and analytically. These complementary approaches address different regimes in the parameter space, and different questions about the evolving surface morphology. Numerically, we have simulated the macroscale PDE using a semi-implicit time scheme and the finite element method in space. We observe a 2D?1D transition under ADL kinetics, indicating that the tensor mobility distinguishes between step- normal and step-parallel fluxes. This transition appears to be accelerated by the application of an electric field in the direction of the longer initial wavelength. Our understanding of the transition phenomena is incomplete, but we hypothesize that a key contribution is the initially enhanced magnitude of longitudinal fluxes relative to transverse fluxes under ADL kinetics. The relevance of our numerical results to real materials is not immediate, owing to the restricted subset of the parameter space (g1/g3 lessmuch 1) that the finite element method could handle reliably. Preliminary computations with nonzero but small g1/g3 suggest that the 2D?1D transition persists continuously as g1/g3 is allowed to take on small positive values. This continuous dependence on g1/g3 offers the promise of an eventual connection with real materials, or even engineered metamaterials. Analytically, we have looked for approximately separable solutions of the macroscale PDE for surface height, predicting novel decay laws for some regions 215 of the parameter space. We have studied the possibility of a boundary layer near a facet edge, finding a boundary layer width that scales as a power of g1/g3 in the regime where g1/g3 is large. We have also characterized the effect of an electric field on the slope profile near a facet. These latter two predictions offer the most promis- ing connection with real materials, whose parameter values challenge the reliability of our numerical method. 11.2 Open questions An ambitious research program inevitably raises more questions than it an- swers. We list below some of the pending issues that this thesis could not address. ? Stabilizing the numerics for larger values of g1/g3. As indicated in Chapter 10, we lack a reliable simulation of the macroscale PDE for values of the parameter g1/g3 of the order of unity. The usual way to address this nu- merical instability is by means of a preconditioning matrix. An open problem is to develop an accurate and efficient preconditioner for the matrix equation given by the finite element method. ? Higher-order time stepping. Much of the discrepancy between different simulations could be attributed to the choice of adaptive time step [6]. In contrast, the calculated decay of height peak and energy were rather robust with respect to mesh size [6]. Perhaps a fruitful direction to pursue is the modification of our algorithm to reduce the time discretization errors. ? Comparison of macroscale simulations with step flow in 2D. Valida- 216 tion of the macroscale theory is commonly done by comparing its predictions with those of an atomistic simulation, as in the work of Shenoy et al. [107, 108]. Since our derivation of the macroscale PDE begins with the generalized BCF model rather than atomistic physics, it might be instructive to compare a step flow simulation with the results of this thesis. The step-flow model of Weeks et al. [122], which has been implemented numerically by Kan et al. [45], offers a promising means of comparison. A crucial question is whether the material parameters used in simulations of Weeks? model are within the scope of our finite element method. (Note that extensive comparisons of step flow with PDEs have been carried out for one-dimensional geometries [25, 79].) ? Boundary conditions at facet edges. As discussed in Chapter 9, the vari- ational formulation of our evolution equation does not respect the boundary conditions arising from the microstructure at a facet edge. An open problem is the development of a hybrid approach, which would couple the variational implementation of a fully continuum PDE with collapse time data from the flow of steps near facets. ? Further analysis of the macroscale PDE. All of our numerics assumed a spatially periodic solution whose period remained unchanged throughout the simulation. The analytical justification of this assumption is still lacking. Other qualitative features of the computed surface morphology, such as the non-monotonic decay of height peak, could likewise be subject to analysis. The question of backwards uniqueness in time is also of physical significance, 217 as indicated in the following. ? Inverse problem. To make our modeling more attractive to the materials engineering community, we might want to analyze whether a desired surface pattern can emerge from self-organization of an initial profile, which relaxes near equilibrium according to the flow of steps below roughening. This ques- tion is an obvious analog for our macroscale PDE of the backwards uniqueness property for the heat equation. We essentially want to characterize which ?fi- nal data? can successfully be evolved backward in time by our fourth-order nonlinear PDE. This characterization might involve the Fourier components of the final data (as in the case of the heat equation), or perhaps a more restric- tive set of conditions. If we manage by analytical means to address backwards uniqueness satisfactorily, we might then try to implement a numerical algo- rithm to approximate the solution of the inverse problem. ? Stochastic effects. In this thesis we performed coarse-graining for a com- pletely deterministic BCF-type model. Stochastic effects (especially the effect of noise at the nanoscale) were not studied. A more realistic picture of the nanoscale physics would allow for the fluctuation of step edges by adding noise terms to the step flow equations. An interesting problem is then to derive the macroscale limit from these stochastic differential equations. ? Far-from-equilibrium conditions. Throughout this thesis we assumed the surface to evolve near equilibrium. In this vein we have focused entirely on relaxation phenomena, omitting the effect of material deposition. The under- 218 standing we gain from studying relaxation is expected to carry through to far-from-equilibrium conditions, because relaxation by surface diffusion is al- ways present in evolving surfaces. However, the far-from-equilibrium setting is found to yield qualitatively different macroscale equations [11, 68, 69], which are of interest in their own right and for the connection they promise with modern experiments in epitaxial growth. 219 Appendix A Brief review of the Mellin transform In this Appendix, we introduce the Mellin transform starting from the more familiar Laplace integral transform. Motivated by the eventual use of the Mellin transform in asymptotic expansion of integrals, we include only those technical details needed to facilitate this application. The extension to iterated integrals and inversion formulas in higher dimensions are not needed for our purposes; the reader is directed to [100] for a more thorough discussion of these issues. Many problems in engineering and applied mathematics are cumbersome to solve when expressed in the relevant physical variables. Integral transforms were developed to address this difficulty. The role of an integral transform is to map an equation from the original domain (e.g., physical space or time) into another domain (e.g., wavenumber or frequency space) where the computations are simpler. The solution in the transformed domain is then pulled back to the physical domain by the inverse transform. The Laplace transform is well-known for its use in ordinary differential equa- tions (ODEs): by taking the Laplace transform of a linear ODE, we obtain an algebraic equation for the transform of the ODE solution. After solving the result- ing algebraic equation, we can apply the inverse Laplace transform to express the solution in the original variable. This idea is made rigorous by the following Definition A.0.1. For a function f : R ? R, the two-sided Laplace transform ?f(s) is given by: ?f(s) = integraldisplay ? ?? f(t)e?stdt. (A.1) 220 The original function f(t) is recovered from ?f(s) by the inverse Laplace transform: f(t) = 12pi? integraldisplay ?+?? ???? ?f(s)estds, (A.2) where ? is a real number so that the integration path lies in the region of convergence of ?f(s). For asymptotic expansion of integrals, it turns out to be more convenient to use an integral transform whose kernel is a power function rather than an exponential. In this vein, we let x = et, dx = xdt in (A.1), which yields ?f(s) = integraldisplay ? 0 f(lnx)x?sdxx . (A.3) We introduce F(x) to denote the composite function f(lnx), and ?? = ?s?1 to denote the power of x. The resulting integral, viewed as a function of ?, is called the Mellin transform and written ?F(?), i.e., ?F(?) = integraldisplay ? 0 F(x)x??dx. (A.4) We break up the integral (A.4) into two parts when studying its convergence properties: ?F(?) = parenleftbiggintegraldisplay 1 0 + integraldisplay ? 1 parenrightbigg F(x)x??dx. (A.5) The restrictions on ? for convergence of (A.5) emerge from enforcing boundedness of the two integrals separately. We have a bound for the first integral, vextendsinglevextendsingle vextendsinglevextendsingle integraldisplay 1 0 F(x)x??dx vextendsinglevextendsingle vextendsinglevextendsingle? integraldisplay 1 0 |F(x)|x?Re?dx, (A.6) where Re? denotes the real part of ?. Suppose that F(x) = O(xp) as x ? 0. This first integral converges if p?Re? >?1, or Re? < 1 +p. 221 On the other hand, we have for the second integral vextendsinglevextendsingle vextendsinglevextendsingle integraldisplay ? 1 F(x)x??dx vextendsinglevextendsingle vextendsinglevextendsingle? integraldisplay ? 1 |F(x)|x?Re?dx. (A.7) If F(x) = O(xq) as x ? ?, then the convergence of the second integral requires q?Re? < ?1, or 1 +q < Re?. Together with the first condition on ?, we find the fundamental strip 1+q< Re? < 1+p, which indicates the region in the complex ? plane where the integral defining the Mellin transform converges. The inverse Mellin transform is obtained from ?F(?) by the same contour integral as (A.2): F(x) = f(lnx) = 12pi? integraldisplay ?+?? ???? ?F(?)x??1d?. (A.8) For the asymptotic evaluation of integrals, whereF(x) is defined by an integral with a ?large? or ?small? variable x, this formulation is useful when the Mellin transform ?F(?) is meromorphic. Then the only singularities of ?F are poles in the complex ? plane, and an asymptotic expansion for F can be calculated by applying the Cauchy residue theorem. Care must be taken to ensure that the integration path is deformed to account for the sign of lnx and the strip of convergence. The Mellin transform is useful when F(x) is defined by a definite integral with ?large? or ?small? parameterx, and the main contribution to this integral cannot be attributed to isolated values of the integration variable. Indeed, classical asymptotic approaches, such as the steepest descent or the stationary phase methods, are usually sufficient in such cases. When the contribution to F(x) comes from an extended interval in the domain of integration, the Mellin transform helps by mapping this contribution to an isolated singularity in the transformed ? space. The required 222 calculation is often easier in the transformed domain, as we see in Chapter 4 when studying the integral for elastic interaction energy. 223 Appendix B Nondimensionalization of the evolution equations In this chapter we rewrite the evolution equations for the macroscale sur- face height in nondimensional form. These equivalent representations allow for a more systematic numerical treatment of relaxation experiments. We find that for relaxation of initially sinusoidal profiles, any collection of material and geometric parameters can be thought of as a point p in R6, with coordinates ?x/?y, h0/?x, ?xD?1v, q, g1/g3. If a certain crystal surface is specified, then q and g1/g3 are not independent; both quantities are related through the temperature T. If the crystal surface is not specified, then varying q and g1/g3 independently has the same effect as considering many different possible materials over a wide range of temperatures. The subset of R6 corresponding to physically realizable parameters can be fur- ther partitioned according to the distinct morphological changes that the macroscale PDE allows. One such partition, distinguishing between initial data that supports a 2D ? 1D transition and initial data that prohibits such a transition, has been illus- trated in [7] as cross-sections in the q,? plane, where ? = ?x/?y is the aspect ratio. To extract from these plots the initial conditions for an experimental comparison with real materials, we need to undo the nondimensionalization of the simulated equations. Undoing this nondimensionalization will determine all the terms of the 224 original macroscale PDE, rewritten here in the universal form ?h ?t = ?CsDdiv braceleftbigg ?? bracketleftbigg?? kBT ?D ?1v parenleftbigg 1 + ?k BT parenrightbiggbracketrightbiggbracerightbigg , (B.1) ? = ??g3 div braceleftbiggg 1 g3 ?h |?h| +|?h|?h bracerightbigg , (B.2) where ? is the dimensionless mobility tensor. The initial data has characteristic lengths ?x, ?y in the basal plane, and h0 in the vertical direction. Accordingly, we define nondimensional spatial variables by ?x = x/?x, (B.3) ?y = y/?y, (B.4) ?h = h/h0. (B.5) To compute spatial derivatives, we define the operators ?? and div? by ?? = (??x,???y), (B.6) div? = ??x +???y, (B.7) which are related to the usual operators ? and div in an obvious way: ?? = ?x? and div? = ?x div. (B.8) The nondimensional chemical potential ?? and drift velocity u are defined as ?? = ?k BT , (B.9) u = ?xD?1v. (B.10) Applying (B.3)?(B.8) to (B.2), we find ? = ??g3h 2 0 ?3x div? parenleftBigg g1?2x g3h20 ???h |???h| +|?? ?h|???h parenrightBigg . (B.11) 225 The prefactor in (B.11) requires us to find ?? by solving ?? = ? ?g3h 2 0 kBT?3x div? parenleftBigg g1?2x g3h20 ???h |???h| +|?? ?h|???h parenrightBigg . (B.12) Substituting (B.3)?(B.10) into the PDE (B.1), we find ??h ?t = ?CsD ?2xh0 div?{??[?????u(1 + ??)]}. (B.13) The prefactor in (B.13) suggests that we choose a time unit t0 given by t0 = h0? 2 x ?CsD, (B.14) so that ?t = t/t0 is nondimensional. In light of (B.12), (B.13), the numerical implementation of the macroscale height evolution with drift u negationslash= 0 requires that we rescale three coefficients in the code: the line tension coefficientg1/g3, the material parameterq appearing in ?, and the prefactor in (B.12). The nondimensional time can absorb only the prefactor in (B.13). In contrast, for u = 0 the nondimensional time can also absorb the prefactor of (B.12), thereby reducing the number of rescaled coefficients to two: g1/g3 and q. This procedure is adopted in Chapter 10, where our time unit is chosen as t0 = kBT?5x?2Csg3Dh0. 226 Bibliography [1] Abramowitz M and Stegun I 1964, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, (10th ed., National Bureau of Standards, 1972). [2] Atkinson C and Wilmott P 1994, Interactions of continuous distributions of ledges, Proc. Roy. Soc. London Ser. A 446, 277?287. [3] Balykov L and Voigt A 2006, A 2+1 -dimensional terrace-step-kink model for epitaxial growth far from equilibrium, Multiscale Model. Simul. 5, 45?61. [4] Blakely J, Umbach C and Tanaka S 1997, Atomic steps in the decay of 1- and 2-dimensional gratings, in Dynamics of Crystal Surfaces and Interfaces, ed. Duxbury P and Pence T, Plenum, New York, NY. [5] Boas M L 1984 Mathematical Methods in the Physical Sciences, Wiley, New York, NY. [6] Bonito A 2007?2009, private discussions. [7] Bonito A, Nochetto R, Quah J and Margetis D 2009, Self-organization of de- caying surface corrugations: A numerical study, Phys. Rev. E. 79, 050601/1?4. [8] Bonito A, Nochetto R, Margetis D, in preparation. [9] Braess D 2001, Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics, Cambridge University Press, Cambridge, UK. [10] Burton W K, Cabrera N, and Frank F C 1951, The growth of crystals and the equilibrium structure of their surfaces, Philos. Trans. R. Soc. London, Ser. A 243 299?358. [11] Caflisch R E, E W, Gyure M F, Merriman B and Ratch C 1999, Kinetic model for a step edge in epitaxial growth, Phys. Rev. E 3, 6879?6887. [12] Chang J, Pierre-Louis O and Misbah C 2006, Birth and morphological evolution of step bunches under electromigration, Phys. Rev. Lett., 96, 195901/1?4. [13] Danker G, Pierre-Louis O, Kassner K and Misbah C 2004, Peculiar effects of anisotropic diffusion on dynamics of vicinal surfaces, Phys. Rev. Lett. 93, 185504/1?4. [14] Degawa M, Minoda H, Tanishiro Y and Yagi K 2001, In-phase step wandering on Si(111) vicinal surfaces:Effect of direct current heating tilted from the step- down direction, Phys. Rev. B 63, 045309/1?8. 227 [15] Dufay M, Frisch T and Debierre J-M 2007, Role of step-flow advection during electromigration-induced step bunching, Phys. Rev. B 75, 241304/1?4. [16] E W and Yip N K 2001, Continuum theory of epitaxial crystal growth I, J. Statist. Phys. 104, 221?253. [17] Ehrlich G and Hudda F 1966, Atomic view of surface self diffusion: Tungsten on tungsten, J. Chem. Phys. 44, 1039?1099. [18] Erd?elyi A 1981, ed., Bateman Manuscript Project: Higher Transcendental Functions, Vol. I, Krieger, Malabar, FL. [19] Erd?elyi A 1981, ed., Bateman Manuscript Project: Higher Transcendental Functions, Vol. II, Krieger, Malabar, FL. [20] Erlebacher J, Aziz M J, Chason E, Sinclair M B and Floro J A 2000, Non- classical smoothening of nanoscale surface corrugations, Phys. Rev. Lett. 84, 5800?5803. [21] Eshelby J D 1956, The continuum theory of lattice defects, Solid State Phys. 3, 79?144. [22] Evans J W, Thiel P A and Bartelt M C 2006, Morphological evolution during epitaxial thin film growth: Formation of 2D islands and 3D mounds, Surf. Sci. Reports 62 1?128. [23] Evans J W 2008, Step dynamics modeling of multilayer growth, minisymposium 55, SIAM conference Philadelphia, PA. [24] Evans L C 2002, Partial Differential Equations, AMS, Providence, RI. [25] Fok P-W 2006, Simulation of Axisymmetric Stepped Surfaces with a Facet, MIT Ph.D. thesis. [26] Fok P-W, Rosales R R, and Margetis D 2007, Unification of step bunching phenomena on vicinal surfaces, Phys. Rev. B 76, 033408/1?4. [27] Frisch T and Verga A 2005, Kinetic step bunching instability during surface growth, Phys. Rev. Lett., 94, 226102/1?4. [28] Fu E S, Liu D-J, Johnson M D, Weeks J D, and Williams E D 1997, The effective charge in surface electromigration, Surf. Sci. 385, 259?269. [29] Ghez R and Iyer S S 1988, The kinetics of fast steps on crystal surfaces and its application to the molecular beam epitaxy of silicon, IBM J. Res. Develop. 32, 804?818. [30] Gruber E E and Mullins W W 1967, On the theory of anisotropy of crystalline surface tension, J. Phys. Chem. Solids 28, 875?887. 228 [31] Hager J and Spohn H 1995, Self-similar morphology and dynamics of periodic surface profiles below the roughening transition, Surf. Sci. 324, 365?372. [32] Hecht F, Pironneau O, Le Hyaric A and Ohtsuka K 2008, FreeFem++ (online documentation), accessed May 2008 at http://www.freefem.org/ff++/ftp/freefem++doc.pdf. [33] Herring C 1949, Effect of change of scale on sintering phenomena, J. Appl. Phys. 21, 301. [34] Herring C 1951, Surface tension as a motivation for sintering, in The Physics of Powder Metallurgy, ed W E Kingston, McGraw Hill, New York, NY, 143?179. [35] Ho P S and Kwok T 1989, Electromigration in metals, Rep. Prog. Phys. 52, 301?348. [36] Israeli N and Kandel D 1998, Profile scaling in decay of nanostructures, Phys. Rev. Lett. 80, 3300?3303. [37] Israeli N and Kandel D 1999, Profile of a decaying crystalline cone, Phys. Rev. B 60, 5946?5962. [38] Israeli N and Kandel D 2000, Decay of one-dimensional surface modulations, Phys. Rev. B 62, 13707?13717. [39] Israeli N and Kandel D 2002, Comment on ?Nonclassical smoothening of nanoscale surface corrugations?, Phys. Rev. Lett. 88, 169601. [40] Jackson J D 1999, Classical Electrodynamics, Wiley, New York, NY. [41] Jayaprakash C, Saam W F, and Teitel S 1983, Roughening and facet formation in crystals, Phys. Rev. Lett. 50, 2017?2020. [42] Jayaprakash C, Rottman C, and Saam W F 1984, Simple model for crystal shapes: Step-step interactions and facet edges, Phys. Rev. B 30, 6549?6554. [43] Jeong H-C and Williams E D 1999, Steps on surfaces: Experiments and theory, Surf. Sci. Reports 34, 171?294. [44] Jo?os B, Einstein T L, and Bartelt N C 1991, Distribution of terrace widths on a vicinal surface within the one-dimensional free-fermion model. Phys. Rev. B 43, 8153?8162. [45] Kan H-C, Kwon T and Phaneuf R J 2008, Effect of length scales in directing step bunch self-organization during annealing of patterned vicinal Si(111) surfaces: Comparison with a simple near-equilibrium model, Phys. Rev. B 77, 205401/1? 6. [46] Kandel D and Kaxiras E 1996, Microscopic theory of electromigration on semi- conductor surfaces, Phys. Rev. Lett. 76, 1114?1117. 229 [47] Kashima Y 2005, A subdifferential formulation of fourth order singular diffusion equations, Adv. Math Sci. Appl. 14, 49?74. [48] Keefe M E, Umbach C C and Blakely J M 1994, Surface self-diffusion on Si from the evolution of periodic atomic step arrays, J. Phys. Chem. Solids 55, 965?973. [49] Kevorkian J and Cole J D 1996, Multiple scale and singular perturbation meth- ods, Springer, New York, NY. [50] Kosevich A M 1972, Osnovy mechaniki kristallichesko?? (Fundamentals of Crys- tal Lattice Mechanics), Nauka, Moscow. [51] Krug J 1997, On the shape of wedding cakes, J. Statist. Phys. 87, 505?518. [52] Krug J 2005, Introduction to step dynamics and step instabilities, in Multiscale Modeling of Epitaxial Growth, ed A Voigt, International Series of Numerical Mathematics 149, Birkh?auser, Basel, Switzerland, 69?96. [53] Krug J, Tonchev V, Stoyanov S, and Pimpinelli A 2005, Scaling properties of step bunches induced by sublimation and related mechanisms, Phys. Rev. B, 71, 045412/1?15. [54] Krug J 2007, private communication. [55] Krug J 2009, Nonlinear dynamics of surface steps, preprint, cond-math.mtrl- sci/arXiv:0810.5749. [56] Kubo R and Nagamiya T 1969, Solid State Physics ed R S Knox, McGraw Hill, New York, NY. [57] Kukta R V, Peralta A, and Kouris D 2002, Elastic interaction of surface steps: Effect of atomic-scale roughness, Phys. Rev. Lett. 88, 186102/1?4. [58] Landau L D and Lifshitz E M 1965, Teoriya uprugosti, Nauka, Moscow. (Theory of Elasticity, 2nd ed., Pergamon Press, Oxford, 1970). [59] Latyshev A V, Aseev A L, Krasilnikov A B, and Stenin S I 1989, Transfor- mations on clean Si(111) stepped surface during sublimation, Surf. Sci. 213, 157?169. [60] Liu D-J and Weeks J D 1998, Quantitative theory of current-induced step bunch- ing on Si(111), Phys. Rev. B 57, 14891?14900. [61] Liu D-J, Weeks J D, and Kandel D 1998, Current-induced step bending insta- bility on vicinal surfaces, Phys. Rev. Lett. 81, 2743?2746. [62] Lo T S and Kohn R V 2002, A new approach to the continuum modeling of epitaxial growth: slope selection, coarsening, and the role of the uphill current, Physica D 161, 237?257. 230 [63] Marchenko V I and Parshin A Ya 1980, Elastic properties of crystal surfaces, Soviet Phys. JETP 52, pp 129?131. [64] Margetis D, Aziz M J and Stone H 2005, Continuum approach to self-similarity and scaling in morphological relaxation of a crystal with a facet, Phys. Rev. B 71, 165432/1?22. [65] Margetis D, Fok P-W, Aziz M J and Stone H A 2006, Continuum theory of nanostructure decay via a microscale condition, Phys. Rev. Lett. 97, 096102/1? 4. [66] Margetis D and Kohn R V 2006, Continuum relaxation of interacting steps on crystal surfaces in 2+1 dimensions, Multiscale Model. Simul. 5, 729?758. [67] Margetis D 2007, Unified continuum approach to crystal surface morphological relaxation, Phys. Rev. B 76, 193403/1?4. [68] Margetis D and Caflisch R E 2008, Anisotropic step stiffness from a kinetic model of epitaxial growth, Multiscale Model. Simul. 7, 242?273. [69] Margetis D and Kohn R V, Continuum theory of material deposition on stepped surfaces in 2+1 dimensions, in preparation. [70] Margetis D 2009, Homogenization of reconstructed crystal surfaces: Fick?s law of diffusion, Phys. Rev. E 79, 052601/1?4. [71] Margetis D 2009, Steps in 2D, lecture notes. [72] M?etois J J, Heyraud J C, and Stoyanov S 2001, Step flow growth of vicinal (111) Si surface at high temperatures: step kinetics or surface diffusion control, Surf. Sci. 486, 95?102. [73] Michely T and Krug J 2004, Islands, Mounds and Atoms, Patterns and Pro- cesses in Crystal Growth Far from Equilibrium, Springer, Heidelberg. [74] Misbah C and Pierre-Louis O 1996, Pulses and disorder in a continuum version of step-bunching dynamics, Phys. Rev. E, 53, R4318?R4321. [75] Mullins W W 1957, Theory of thermal grooving, J. Appl. Phys. 28 333?339. Mullins W W 1959, Flattening of a nearly plane solid surface due to capillarity, J. Appl. Phys. 30 77?83. [76] Najafabadi R and Srolovitz D J 1994, Elastic step interactions on vicinal sur- faces of fcc metals, Surf. Sci. 317, 221?234. [77] Nochetto R H 2007?2009, private discussions. [78] Nozi`eres P 1987, On the motion of steps on a vicinal surface, J. Phys. (France) 48, 1605?1608. 231 [79] Odisharia I V 2006, Simulation and Analysis of the Relaxation of a Crystalline Surface, NYU Ph.D. thesis. [80] Ozdemir M and Zangwill A 1990, Morphological equilibration of a corrugated crystalline surface, Phys. Rev. B 42, 5013?5024. [81] Ozdemir M and Zangwill A 1992, Morphological equilibration of a facetted crys- tal, Phys. Rev. B 45, 3718?3729. [82] Paulin S, Gillet F, Pierre-Louis O and Misbah C 2001, Unstable step meandering with elastic interactions, Phys. Rev. Lett. 86, 5538?5541. [83] Pedemonte L, Bracco G, Boragno C, Buatier de Mongeot F and Valbusa U 2003, Smoothing of nanoscale surface ripples studied by He atom scattering, Phys. Rev. B 68, 115431/1?6. [84] Pfeiffer R S and Mahan G D 1993, Mean-field theory of elastic dipoles on a face-centered-cubic lattice, Phys. Rev. B 48, 669?671. [85] Phaneuf R 2007?2009, private discussions. [86] Pierre-Louis O, Misbah C, Saito Y, Krug J and Politi P 1998, New nonlinear evolution equation for steps during molecular beam epitaxy on vicinal surfaces, Phys. Rev. Lett. 80, 4221?4224. [87] Pierre-Louis O 2001, Continuum model for low temperature relaxation of crystal steps, Phys. Rev. Lett. 87, 106104/1?4. [88] Pierre-Louis O 2003, Step bunching with general step kinetics: stability analysis and macroscopic models, Surf. Sci. 529, 114?134. [89] Pierre-Louis O 2003, Phase field models for step flow, Phys. Rev. E 68, 021604/1?19. [90] Pierre-Louis O and M?etois J 2004, Kinetic step pairing, Phys. Rev. Lett. 93, 165901/1?4. [91] Pierre-Louis O 2005, Dynamics of crystal steps, Comptes Rendus Physique, 6, 11?21. [92] Pierre-Louis O 2006, Local electromigration model for crystal surfaces, Phys. Rev. Lett. 96, 135901/1?4. [93] Pimpinelli A and Villain J 1998, Physics of Crystal Growth, Cambridge Uni- versity Press, Cambridge, UK. [94] Politi P and Villain J 1996, Ehrlich-Schwoebel instability in molecular-beam epitaxy: A minimal model, Phys. Rev. B 54, 5114?5129. 232 [95] Politi P and Misbah C 2006, Nonlinear dynamics in one dimension: A criterion for coarsening and its temporal law, Phys. Rev. E 73, 036133/1?15. [96] Quah J and Margetis D 2008, Anisotropic diffusion in continuum relaxation of stepped crystal surfaces, J. Phys. A: Math. Theor. 41, 235004/1?18. [97] Quah J and Margetis D 2009, Electromigration in macroscopic relaxation of stepped surfaces, submitted to Multiscale Model. Simul. [98] Rettori A and Villain J 1988, Flattening of grooves on a crystal surface: A method of investigation of surface roughness, J. Phys. (France) 49, 257?267. [99] http://www.sandia.gov/surface science/stm [100] Sasiela R J and Shelton J D 1993, Integral evaluation using Mellin transform methods, J. Math. Phys. 34, 2572?2617. [101] Sato M, Uwaha M, Saito Y, and Hirose Y 2002, Step wandering induced by the drift of adatoms in a conserved system, Phys. Rev. B, 65, 245427/1?6. [102] Sato M, Uwaha M, and Saito Y 2005, Evaporation and impingement effects on drift-induced step instabilities on a Si(001) vicinal surface, Phys. Rev. B 72, 045401/1?9. [103] Schimschak M and Krug J 1997, Surface electromigration as a moving bound- ary value problem, Phys. Rev. Lett. 78, 278?281. [104] Schulze T P and E W 2001, A continuum model for the growth of epitaxial films, J. Cryst. Growth 222, 414?425. [105] Schwoebel R L and Shipsey E J 1966, Step motion on crystal surfaces, J. Appl. Phys. 37, pp 3682?3686. [106] Shenoy V B and Freund L B 2002, A continuum description of the energetics and evolution of stepped surfaces in strained nanostructures, J. Mech. Phys. Solids 50, 1817?1841. [107] Shenoy V B, Ramasubramaniam A, Freund L B 2003, A variational approach to nonlinear dynamics of nanoscale surface modulations, Surf. Sci. 529, 365? 383. [108] Shenoy V B, Ramasubramaniam A, Ramanarayan H, Tambe D T, Chan W-L and Chason E 2004, Influence of step-edge barriers on the morphological relax- ation of nanoscale ripples on crystal surfaces, Phys. Rev. Lett. 92, 256101/1?4. [109] Spohn H 1993, Surface dynamics below the roughening temperature, J. Phys. I 3, 69?81. [110] Stangl J, Hol?y V, and Bauer G 2004, Structural properties of self-organized semiconductor nanostructures, Rev. Mod. Phys. 76, 725?783. 233 [111] Stoyanov S 1991, Electromigration induced step bunching on Si surfaces ? How does it depend on the temperature and heating current direction?, Jpn. J. Appl. Phys. 30, 1?6. [112] Stoyanov S 2000, Scaling properties of step bunches formed on vicinal crystal surfaces during MBE growth, Surf. Sci. 464, L715?L718. [113] Tanaka S, Bartelt N C, Umbach C C, Tromp R M and Blakely J M 1997, Step permeability and the relaxation of biperiodic gratings on Si(001), Phys. Rev. Lett. 78, 3342?3345. [114] Tersoff J, Johnson M D and Orr B G 1997, Adatom densities on GaAs: Evi- dence for near-equilibrium growth, Phys. Rev. Lett. 78, 282?285. [115] Th?urmer K, Reutt-Robey J E, Williams E D, Uwaha M, Emundts A, and Bonzel H P 2001, Step dynamics in 3D crystal shape relaxation, Phys. Rev. Lett. 87, 186102/1?4. [116] Tricomi F G 1985, Integral Equations, Dover, New York, NY. [117] Uwaha M 1988, Relaxation of crystal shapes caused by step motion, J. Phys. Soc. Japan 57, 1681?1686. [118] Uwaha M and Watanabe K 2000, Decay of an island on a facet via surface diffusion, J. Phys. Soc. Japan 69, 497?503. [119] Villain J 1991, Continuum models for crystal growth from atomic beams with and without desorption, J. Phys. I 1, 19?42. [120] Vvedensky D D, Zangwill A, Luse C N and Wilby M R 1993, Stochastic equations of motion for epitaxial growth, Phys. Rev. E 48, 852?862. [121] Vvedensky D D 2004, Multiscale modeling of nanostructures, J. Phys.: Con- dens. Matter 16, R1537?R1576. [122] Weeks J D, Liu D-J, and Jeong H-C 1997, Two-dimensional models for step dynamics in Dynamics of Crystal Surfaces and Interfaces, ed. Duxbury P and Spence T, Plenum, New York, NY. [123] Xiang Y 2002, Derivation of a continuum model for epitaxial growth with elasticity on vicinal surface, SIAM J. Appl. Math., 63, 241?258. [124] Zhao T, Weeks J D, and Kandel D 2004, Unified treatment of current-induced instabilities on Si surfaces, Phys. Rev. B, 70, 161303/1?4. [125] Zhao T, Weeks J D, and Kandel D 2005, From discrete hopping to continuum modeling on vicinal surfaces with applications to Si(001) electromigration, Phys. Rev. B 71, 155326/1?9. Zhao T and Weeks J D 2005, A two-region diffusion model for current-induced instabilities of step patterns on vicinal Si(111) surfaces, Surf. Sci. 580, 107?121. 234