Scaling and Dynamic Stability of Model Vicinal Surfaces - Crystal

Dec 12, 2018 - Based on a novel model we analyze the step bunching instability on ... an integrated modeling approach to the fundamental problem of vi...
0 downloads 0 Views 1MB Size
Subscriber access provided by Columbia University Libraries

Article

Scaling and dynamic stability of model vicinal surfaces Filip Krzy#ewski, Magdalena Za#uska-Kotur, Anna Krasteva, Hristina Popova, and Vesselin Tonchev Cryst. Growth Des., Just Accepted Manuscript • DOI: 10.1021/acs.cgd.8b01379 • Publication Date (Web): 12 Dec 2018 Downloaded from http://pubs.acs.org on December 14, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Crystal Growth & Design

Scaling and dynamic stability of model vicinal surfaces Filip Krzyżewski1, Magdalena Załuska-Kotur1, Anna Krasteva2, Hristina Popova3 and Vesselin Tonchev 4, 1 * 1Institute

of Physics, Polish Academy of Sciences, al. Lotników 32/46, 02-668 Warsaw, Poland. of Electronics, Bulgarian Academy of Sciences, 72 Tzarigradsko chaussee blvd, 1784 Sofia, Bulgaria 3R. Kaischew Institute of Physical Chemistry, Bulgarian Academy of Sciences, Acad. G. Bonchev Str., block 11, 1113 Sofia, Bulgaria 4 Faculty of Physics, Sofia University, 5 James Bourchier Blvd., 1164 Sofia, Bulgaria *Corresponding author: [email protected] 2Institute

We propose an integrated modelling approach to the fundamental problem of vicinal crystal surfaces destabilized by step-down (SD) and step-up (SU) currents with focus on both the initial and the intermediate stages of the process. We reproduce and analyze quantitatively the step bunching (SB) instability, caused by the two opposite drift directions in the two situations of irreversible step motion mediating sublimation and growth. For this reason we develop further our atomistic scale model (vicCA) of vicinal crystal growth (Gr) destabilized by SD drift of the adatoms in order to account for also the vicinal crystal sublimation (Sbl) and the SU drift of the adatoms as an alternative mode of destabilization. For each of the four possible cases - Gr+SD, Gr+SU, Sbl+SD, Sbl+SU, we find a self-similar solution - the time-scaling of the number of steps in the bunch N, N  2 T 3 , where T is the time, rescaled with a combination of model parameters. In order to study systematically the emergence of the instability, we use N further as a measure and probe the model’s stability against SB on a dense grid of points in the parameter space. Stability diagrams are obtained, based on simulations running to fixed moderate re-scaled times and with small-size systems. We confirm the value of the numerical pre-factor in the time scaling of N, 2

3 by results obtained from systems of ordinary differential equations for the step velocity that

contain, in contrast to vicCA, step-step repulsions. This last part of our study provides also the possibility to distinguish between diffusion-limited and kinetic-limited versions of the step bunching phenomena.

1 Introduction Step bunching (SB) is a phenomenon archetypical for a class of surface instabilities - SB (of straight steps), step meandering (SM), simultaneous SB and SM, step bending, surface faceting. These and, ACS Paragon Plus Environment

Crystal Growth & Design 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 23

with increasing share, their practical aspects and applications, were intensively studied in the years after 1990, in parallel with the dramatic development and sophistication of the experimental techniques for surface analysis, the available computational power and the ever-increasing demand from technological side. Although the first report on surface patterning under electric fields was reported as early as 1938 by Johnson1 who studied tungsten filaments used in the incandescent lamps, the term step bunching itself appears in 50’s2-4 (and the references therein). Later the interest concentrated on almost flat on atomic scale tungsten surfaces5 and the first theory of electromigration induced surface instability was advanced6. In a rather different context, the interest in SB is closely connected with its, sometimes dramatic, impact on the processes of nanostructure layer-by-layer epitaxial growth as realized in various deposition techniques having their industrial realizations. Nowadays, the keyword step bunching bridges studies on material systems as diverse as CH3NH3PbI3 7, GaN 8, AlN 9, AlGaN

10,

SiC

11-14,

tetracarboxylic acid dianhydride)/Ag17, pyronaridine/heme ferritin

24,

graphene

18,

15, 16,

PTCDA (perylene

19 SrRuO3/(001) SrTiO3 , KDP

20-23

,

etc. The phenomenon however is not limited to the epitaxial growth only. It has been

shown that SB development in crystal growth from solution SB development is controlled by fluid flow25,26. Models of step bunching are now part of novel approaches to exotic localizations such as the sidewalls of the nanowires27. The onset of the modern stage of SB studies could be dated back to 1989 when Latyshev et. al28 report the observation of the phenomenon on sublimating Si(111) surfaces destabilized by both step-up and step-down direct currents used to heat the surface but by only one of these for any fixed temperature. It was hypothesized soon after based on a simple Burton-Cabrera-Frank29 (BCF) type of model that it is the directional asymmetry (bias) in the diffusion of the charged surface atoms (adatoms) that makes the motion of the steps from the vicinal surface unstable30-32. As a result, the regular step distribution of the vicinal surface is broken and a surface that is sequence of groups of steps (bunches) separated by large terraces almost free of steps, hills and valleys structure, results. In the next years the phenomenon of SB was investigated actively by experimental techniques33-37 analyzed theoretically38-49, and by numerical simulations50-58 It becomes clear only quite recently that SB can be caused at a given temperature by both current directions across the steps but on a metal surface – that of tungsten, W(110)59. Besides semiconductor and metal crystal surfaces also the insulator ones exhibit the same, adatom electromigration induced, type of behaviour60. A special case of SB is the macrostep formation – when the distances between the steps do not exist anymore due to a strong destabilization effect and/or weak step-step repulsion, this phenomenon was studied both experimentally61 and numerically57,58,62. The mechanisms of all enlisted above phenomena are far from being unified based on clear understanding of the interplay of the surface elementary processes. In the narrower domain of SB ACS Paragon Plus Environment

Page 3 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Crystal Growth & Design

induced by adatom electromigration, the one caused by step-down flux of adatoms is easy to understand based on the concepts within the BCF paradigm - the step velocity is obtained from the gradients at the steps of the diffusion field on a single terrace. Thus, when the contributions of the two terraces encircling the step are uneven the step motion could be unstable. As a result bunches of steps moving close to each other (SB) start to generate if the contribution of the terrace that is behind the step is larger than from the one ahead. This mechanism was used also to explain the step-up drift induced SB, for which change of the sign of adatom charge, the nature of the bias, etc., were evoked. The present paradigm was formed with the important input from S.Stoyanov63 – the destabilization due to a step-up adatom drift is stemming from a specific step property – the so called step transparency63-65 (or permeability66,67) – the adatoms, during their diffusion along the surface, almost do not feel the transparent steps as special points on the surface, and this prevents the steps from being sinks/sources of adatoms. This property allows for the adatom fluxes to carry information of the local ordering along the surface over many initial terrace distances across the steps 65. None of the theories was able to comprise the mechanism of uneven terrace contribution to the step motion and the step transparency from a unified viewpoint. Some other deficiencies of the existing paradigm are the unfinished program of finding the full set of scaling exponents/relations quantifying the intermediate stages of the instability, initiated by A. Pimpinelli et al.40 and partially continued by J. Krug et al.68 based on model equations of hydrodynamic type for the local height of the crystal surface that are assumed applicable for different sources of the instability – SD electromigration and presence of Ehrlich-Schwoebel effects (normal and inversed)40,68, and the incorporation of the simultaneous step bunching and step meandering69-71. Developing further the recently introduced 1D atomistic scale model of growing vicinal surface destabilized by SD drift of adatoms57,58, a realization of fine-tuning of step transparency bound to the adatom diffusion and step kinetics, we are now able to reproduce the instability caused by the two opposite drift directions in the two fundamental situations of step motion – vicinal sublimation and vicinal growth. 2 The Model Our model (vicCA)57,58 is a conceptual realization of a vicinal crystal surface with adatoms on it and combines Monte Carlo (MC) with Cellular Automaton (CA) modules. It allows for quantitative analysis of the pattern formation dynamics by applying a modified monitoring protocol42. While the automaton rule acts in a parallel fashion providing simultaneous realization of all growth/sublimation events, the rule of the adatom diffusion realized in between is simulated by the MC module which acts in a serial mode, dealing with adatom after adatom. The adatoms are stored in a separate table with surface concentration c0 and diffuse on top of the vicinal surface without “feeling” the steps. Their diffusion is influenced by a jump asymmetry (bias), thus the jump ACS Paragon Plus Environment

Crystal Growth & Design 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 23

probability pJ in each direction (left, right) is dependent on δ – pJ  1 2  ,1 2   , with   1 2 . Negative (positive) values of δ stand for bias pointing in the step-up (step-down) direction. In a crystal table only the top atoms from the crystal lattice are stored. Their corresponding heights are controlled, forming the vicinal stairway descending to the right with initial, vicinal distance between sequential steps of l0. When crystal height changes from site to site by one unit cell it is understand as a single step, and when it changes by more than one unit cell as a macrostep. Crystal growth events are realized with the CA module which prescribes that an adatom right to a step sticks to it with a probability pG, the step moves one lattice site to the right and the adatom is deleted from the table of adatoms. Only during the execution of the growth rule adatoms enter into contact with the steps. In the present version of vicCA another part of the composite CA module rules the sublimation - a process opposite to the growth and causing step motion to the left. Technically, the sublimation of the crystal surfaces is realized by checking whether there is an adatom on top of a site where single or macrostep is present. If not, adatom detachment is performed from the step with probability pS. When detachment is accepted the lattice height at that site is reduced by one and new adatom is added to the table with the adatoms at that position. In such a way we can simulate sublimation and growth as a reversible processes. Between two growth modules all adatoms try to make nDS diffusional hops but only those that point at a neighboring unoccupied lattice site are performed. When increasing nDS one departs from the diffusion-limited (DL) growth mode towards the kinetics-limited (KL) one and, simultaneously, increases the step transparency. It is important to note that setting the attachment probability pG  1 and detachment probability pS  1 at the steps is another way of increasing step transparency even when nDS=1. Each time step of simulation is finished by setting the adparticle concentration to its initial value c0. Thus, in the sublimation case the density of adatoms at the surface will exceed c0 and the excess is deleted from the table with adatoms at random. In the growth case, when concentration of the adatoms is expired below c0, we add necessary adatoms at random on the unoccupied sites thus maintaining the adatom concentration at c0. More details of the above routine are available in57,58. Summarizing, а single time step of a vicCA model of sublimation(growth) process starts with sublimation (growth) update and then nDS diffusional steps are made after which growth (sublimation) update is performed. The attachment (sticking) and detachment probabilities vary from 0 to 1. Finally, after all growth, diffusional and sublimation processes, adatom concentration is recovered to c0. In order to study the stability of this model against step bunching (SB) we scan the parameter space consisting of c0, nDS,  and probabilities of adatom attachment (detachment) to (from) the steps, while fixing the initial vicinal distance l0=15. Such distance has been chosen as quite small, so our simulations could be performed in finite time and at the same time not very small so system shows full spectrum of possible behavior. In all simulations below two different ACS Paragon Plus Environment

Page 5 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Crystal Growth & Design

system sizes were used: 10000 sites for time scaling plots and 1000 sites for phase diagrams calculations. Note that in the model we do not incorporate repulsion between steps and as a result bunches formed consist of densely located single steps and also of macrosteps. 3 Results 3.1. Time scaling of the bunch size We quantify the similarity of SD and SU biased bunching process by studying the time evolution of the average number of steps in the bunch N and showing that it tends to the same, universal scaling law, can be written in the general form as N  t   a  t 



and we find both the

parameter combinations (characteristic timescale)  and the time-scaling exponent β, as well as the numerical pre-factor (amplitude) a. The dependencies found for each of the four cases studied SD+Gr, SD+Sbl, SU+Gr, and SU+Sbl are shown in Fig 1. The curves were obtained after averaging over three runs for each set of parameters chosen. At first place, it is seen that in every case studied bunch size always achieves, sooner or later, a power law dependence on the time with   1 2 . We could speculate on this observation further – in a growth model of single crystal SiC growing further in Si-Ni flux conditions, Onuma et al.14 found step bunching to occur in diffusion-limited conditions but not in the kinetics limited ones. Based on our results, we could suggest that in the KL case the process may be was not carried to long enough times that would permit the observation of SB. The time-scaling exponent of 1/2 was obtained in many different studies, both theoretical and experimental, and may even seem trivial. We go beyond by showing that the dependencies obtained in all four cases collapse onto the single master curve:

N 2 T 3

(1)

where T  t  , and the parameter combinations  used to rescale the time depend on the bias direction but do not depend on the vicinal process - growth or evaporation. The following two timescales  are used to re-scale the time: SD bias  4l0 nDS ,    2l c 1 c  | | n , SU bias 0 DS  0 0

(2)

The universality and source of the number prefactor 2 3 will be discussed in Section 3.5. Note that the timescale (2) we use for step bunching in the SD case is slightly different than the one we used recently for SD+Gr case58 - in the present study a numerical factor of 4 appears instead of c0. For densities around 0.2 that were studied before the difference is negligible. After investigating wider span of surface densities we have found that the time scaling for the SD cases does not ACS Paragon Plus Environment

Crystal Growth & Design 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 23

depend on the density.

Figure 1. Time dependent bunch size N for the four studied cases and different simulation parameters. Timescales (parameters combinations) used for growth and sublimation differ for SD and SU cases, see eq.(2). Note that here we go well beyond the value of Tgr=500 used to plot most of the stability diagrams in the next Section. At the same time too large values of T, above 5 ∙ 104 lead to the visible fluctuations caused by too small statistics of bunch numbers in the system. In all examples it is set pG=1-co+0.1 and p S  c0 in the case of growing systems, while for sublimating systems pS=co+0.1 and pG=1-co. Other choices of attachment/detachment probabilities give the same scaling, if only allow for step bunching development.

Figure 2 Stability diagrams plotted for three consecutive rescaled times Tgr = 100, 250, 1000, for two different situations in the case of SD bias – irreversible growth/sublimation (upper row) and reversible ones (bottom). It is seen that in general the SB in the KL part (high nDS) emerges at later times than in DL part (low nDS) of the diagrams. It can be also seen that both choices of timescales in eq.(2), are symmetric for c0 and (1-c0)

ACS Paragon Plus Environment

Page 7 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Crystal Growth & Design

interchange what has its consequences in the reflexion symmetry of the diagrams presented in the next Section. There we will use equation (2) to estimate the simulation time necessary to obtain credible stability diagrams – on the one side, the evolution time should be long enough to verify which set of parameters leads to emergence of SB and which does not, and, on the other side, it should not be too long, thus permitting to perform quite more runs with different parameters. Assuming rescaled evolution time Tgr =500 and thus maximal bunch size to be expected is N  25 , we calculate the corresponding number of not rescaled time steps tgr for each parameter combination on the grid. 3.2. Stability diagrams The stability diagram of the vicCA model is studied in the four possible combinations of vicinal process, growth and sublimation, and destabilizations – bias (drift) oriented step-down or step-up. In each point we start with a vicinal surface consisting of equally distributed steps separated by terraces of width l0=15. We start with fixed values of δ=0.2 for the SD cases and δ=0.2 for the SU one, and perform runs of each of them for 500 rescaled time units monitoring the average bunch size N. Thus our stability diagrams contain information about the average bunch size in comparable stages of systems evolution which not always comes after the same time of simulation. In Figure 2 we show time evolution for the stability diagram in the SD case (δ=0.2) in two different situations. Upper panels illustrate the evolution of diagram for the simplest, irreversible dynamics with pG=1 and pS=0 for growth and pS=1 and pG=0 for sublimation. We have checked that SB is obtained only in such conditions for which rates of growth and sublimation are not too high. Those systems that grow or sublimate at high velocity do not form bunches but roughen. In order to ensure slow step motion and observe step bunching for the case of SU drift, one needs to make the process reversible – to realize both attachment and detachment of adatoms to/from the steps at each time step of the run. Thus we assumed that pG=1-co+0.1 and p S  c0 in the case of growing systems. When the system is sublimating we fix step detachment probability as pS=co+0.1 and set pG=1-co. The lower panels of Figure 2 show stability diagrams in the described above reversible dynamics for SD bias. It is clearly seen that the shape of upper and lower diagrams is different but its character does not change much. On increasing time size of bunches become larger and the areas of stable and unstable part of diagram become more pronounced. We decided that Tgr=500 units is a good choice. Contrary to the situation presented in Figure 2 the system under SU bias in irreversible case does not bunch at all, while the stability diagram for the reversible situation evolves in time similarly to the case of SD bias. Below we concentrate our study on the reversible situation only. The results for the stability of the reversible model for Tgr=500 are summarized in Figure 3. ACS Paragon Plus Environment

Crystal Growth & Design 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 23

Panel (a), upper left, presents diagram of the mean bunch size N as a function of nDS and c0 for the case of growth with SD bias. Note that bunching happens for densities lower than 0.5 and number of diffusion steps nDS less than 100. The most effective bunching process is concentrated rather in the corner, for lower nDS and c0. Similar behaviour with (1-c0) playing role of c0 is observed at the sublimated system with SD bias. Stability diagram for that case is presented in Figure 3b. It is a mirror image (with respect to the c0-axis) of the diagram for the growth+SD case, Fig 3a. The stability diagram for this sublimation process shows that SB is strongest at low nDS values and concentrations greater than 0.5. It can be also summarized that low values of nDS (below 40) are the optimal regime for step bunching in both growth and sublimation processes destabilized by SD drift. The results for SU particle drift are shown in Figure 3c and 3d. The areas where SB develops during growth and sublimation are separated. When crystal grows, Figure 3c, the most intensive SB is present at higher adatom concentrations varying between 0.75 and 1 and for broad range of nDS from 2 to values higher than 100. For the sublimating system, Figure 3d, bunches emerge at lower concentrations ranging from 0 to 0.25 and the diagram is again as in the case of SD drift, a mirror image of Figure 3c. It should be stressed that when switching from SD to SU bias bunching instability “moves” to the opposite corner of the diagram: from DL (small nDS) to KL (large nDS) and from c0 to 1-c0. Such symmetry is to be explained by mechanisms that are responsible for step bunching processes and will be discussed below. One can also notice that the shape of parameter areas on the plot where SB is most intensive do not match the ones obtained for the opposite direction of the current. Hence in the present model, the best way to obtain SB at SU and SD bias is to change the other parameters of the simulation. The possible overlapping of SD and SU biased bunching can be examined by investigating the stability diagrams as a function of the drift direction. In Figure 3 such diagrams as a function of the concentration c0 and value of the bias δ are presented. Positive (negative) values of δ cause step down (step up) drift of the adatoms. For switching from growth to sublimation the model parametrization (values of pG and pS) is changed as described above. It is seen that step bunching occurs in the same system, growing or sublimating, for both bias directions, however the ranges of favorable surface densities (reflecting external flux of adatoms) are almost separated. There is only small range of c0 values close to 0.4 where one finds SB with both drift directions.

ACS Paragon Plus Environment

Page 9 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Crystal Growth & Design

a

b

c

d

Figure 3 Stability diagrams of bunch size N in the space of nDS and c0 for the four investigated cases: growth destabilized by SD bias (a), sublimation - by SD bias (b), growth - by SU bias (c) and sublimation - by SU bias (d). Tgr = 500. Note the mirror symmetry between growth and sublimation case for given destabilization factor. As we can see in Figure 3 and Figure 3 the step bunching process is quite general phenomenon and happens for all combinations of growth/sublimation and step-up/step-down drifts. However the parameters for which the instability is observed are different. In general step-down bias induces bunching in DL/step impermeability regime, i.e. for low nDS values, whereas step-up bias should be combined with high nDS values, corresponding to KL/step permeability, in order to produce step bunching. Moreover depending on the type of process and bias low density or high density conditions are preferred. Let us discuss in more details the instability development case by case. The first and the simplest one is growth (sublimation) at a low (high) density destabilized by bias in the SD direction. The step motion there is controlled by the adatoms that stick (detach) to (from) the steps. Thus, only few diffusional hops are to be executed before each sticking (detachment) of an adatom to (from) a step and it is highly probable that each adatom (vacancy) ACS Paragon Plus Environment

Crystal Growth & Design 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 23

arriving at the step will be incorporated into the crystal.

Figure 3 Stability diagram as a function of c0 and δ, l0=15, nDS=10, Tgr = 500. The positive δ values cause SD drift of adatoms. We can assume further that all adatoms (vacancies) that are present on a terrace eventually stick at one of the step edges. With bias  directed SD, >0, the result is that the upper (lower) terrace contributes with more adatoms (vacancies) to the velocity of the i-th step vi than the lower (upper) one and one can write for the case of nDS = 1 when no contributions from the further located terraces are present:

vi 

dxi  c  b1xi  b2 xi1  dt

(3)

where is xi  xi  xi1 the length of the upper terrace, comprised by the steps at xi and xi-1, b1 and b2 are (still unknown) functions of  and

c  c0 ( c  1 c0 ) for growth (sublimation). A simplest

linear stability analysis in which xi  xi  x (and thus xi1  xi1  x ) shows that any positive (negative) perturbation x will grow (decrease) further the length of xi provided b1>b2, thus resulting into the classical instability scenario – “these that have, will have more” (and vice versa), called in the general scientific context “Matthew effect”

71after

Matthew’s 25:29. Since in

the model step-step repulsions are not incorporated, the development of the instability in our vicCA model leads to the formation of macrosteps. Beyond the initial destabilization caused by the uneven adatom fluxes from the two neighboring terraces that attach to (detach from) there is another mechanism of destabilization that comes into the play - when the lengths of the free from steps terraces between the dense step formations are increasing this leads to further destabilization ACS Paragon Plus Environment

Page 11 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Crystal Growth & Design

because the hop asymmetry is amplified by the cumulation of many such hops. Thus, the machinery of the step bunching in the case of SD bias for crystal growth (sublimation) process is well understood in the case when the adatom concentration c0 is low (high). Then the question arises: why high density of adatoms for growth and high vacancy density for sublimation destroys the tendency for step bunching. The reasoning here is as follows: high (low) c0 for growth (sublimation) means that steps move faster and this, combined with the low diffusion rate (nDS close to 1) causes that length of the terraces does not count and no step-step correlation is possible. The mechanism of step bunching caused by SU bias is more complicated. It happens at high (low) adatom concentration for crystal growth (sublimation) and for high diffusion rate, nDS >> 1, which causes the steps to become transparent, because adatoms bypass them during long diffusional pathways. But even then the condition for large nDS is not enough to induce the bunching process but only surface roughening results. Our study shows that it is necessary to slow down the step motion by making it reversible - increasing the probability of sublimation during the growth process and the probability of growth during the sublimation process, Figure 3 and Figure 3 are obtained under such conditions. The scenario of the bunching instability in such a case can be understood as a communication between steps by induced by step motion excess of holes (adatoms) produced at one step, and then transported toward the next one. When velocities of steps and adatoms are correlated speeding up of one step slows down the next one and we observe step bunching process. The development of the instability in the SU case makes it more and more similar to a SD case at the same stage – larger terraces amplify the instability by cumulating the diffusional asymmetry over many jumps. In the same time, at both ends of an increasing terrace there are step formations with increasing collective impermeability – a step bunch of many permeable steps as a coarse grained object is itself increasingly impermeable proportionally to the number of steps in the bunch. Thus the similarity between SD and SU cases increases making them indistinguishable in the intermediate stages of the instability. At the same time the mechanism of bunching for SD and SU bias is completely different. The first originates from the step movement instability caused by density fluctuations, while the second case can be explained by particle and step velocity correlations and in general is more subtle phenomenon. As a result we don’t expect full symmetry between SD and SU cases in stability diagrams in Figures 3, 4. 3.3 Growth and sublimation rate In Figure 4 we show how the growth rate vgr defined as: n vgr  gr t gr

ACS Paragon Plus Environment

(4)

Crystal Growth & Design 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 23

depends on the adatom concentration c0 in all four investigated cases SD and SU for growth and sublimation. In this definition ngr is the number of mono-atomic layers grown for time tgr. In all cases growth rate is low. It increases slowly with concentration of particles for crystal growth or holes for sublimation. In the region of bunching process the growth rate drops down to almost zero in DL regime, while in KL regime it fluctuates around slightly larger value. Above c0=0.5 the number of sites occupied by adatoms is larger than number of empty sites. Diffusion of holes (i.e. empty lattice sites) rather than that of adatoms is realized effectively. At the same time more adatoms are consumed at steps and more new adatoms are added during the compensation stage. Positions of the adatoms become not so well defined. In the KL regime positions of adatoms, due to the fast diffusion, are “spread” also at low concentrations and growth rate increases as the probability of adatom encountered at steps increases. It can be seen however that in the region where step bunching occurs the rate of crystal growth decreases. The negative values of the growth rate in Figure 4b mean sublimation. The system evolution is slower at high adatom concentration. It means slow detachment which is blocked by the presence of adatoms on top of the steps. The decrease of c0 (increase of (1-c0)) leads to a faster sublimation because less adatoms are present at the surface and adatom detachment from steps happens more often. In Figure 4, c and d, we can see how growth and sublimation rate behave as a function of the density in the presence of SU drift. Again two different values of nDS are shown. The rate of both processes is smaller in the bunching region and grows slightly when system evolves keeping the steps regularly distributed. Both regimes can be clearly seen in the velocity plots.

ACS Paragon Plus Environment

Page 13 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Crystal Growth & Design

Figure 4 Concentration dependent growth rate in KL and DL regime at four investigated cases: growth at SD bias (a), sublimation at SD bias (b), growth at SU bias (c) and sublimation at SU bias (d). 3.4 Step trajectories In order to compare both scenarios of step bunching we can analyse also the time evolution of step position. Figure 5 shows step (red lines) and macrostep (green lines) trajectories obtained from the four investigated cases of model evolution. In all plots macrosteps dominate the system’s behavior. Single steps are very rare. They are faster than macrosteps and wander between them mediating a step exchange mechanism between neighbouring bunches. In Figure 5a, b (SD bias), single steps detach from bunches at random moments of time and wander fast to the next ones, while in the system grown with SU current (Figure 5c, d) single steps are visible mostly at the moments preceding bunch connections. This observation is interrelated to the fact that the bunching mechanism for SD bias is based on absorbing all adatoms or holes from the terrace, hence terraces are important, while for SU bias steps communicate between themselves. Hence in this last case connection of two neighbouring bunches together happens with the intensive exchange of individual steps.

ACS Paragon Plus Environment

Crystal Growth & Design 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 23

Figure 5 Step and macrostep trajectories observed at four investigated cases: growth at SD bias (l0=15, c0=0.15, nDS=10, δ=0.2) (a), sublimation at SD bias (l0=10, c0=0.9, nDS=5, δ=0.5) (b), growth at SU bias (l0=15, c0=0.85, nDS=10, δ=−0.2) (c) and sublimation at SU bias (l0=15, c0=0.2, nDS=15, δ=−0.5) (d). 3.5 Numerical analysis of ODE based models In what follows we provide a perspective on just obtained results for the time-scaling of the bunch size from vicCA in order to justify the numerical pre-factor. For this purpose, two models are studied, both based on systems of ordinary differential equations (ODE) for the velocity of each step in the step train. The first one, that of Liu and Weeks (LW) 72 is aimed at the KL regime of the instability and is defined through its equations for the step position xi as: dxi  1  b  1 b    xi    xi1  Ri  2  dt  2 

(5)

where the interplay of the first two terms containing neighboring terrace widths xi  xi  xi1 destabilizes the surface when the parameter b>0. The destabilizing impact of the surface kinetics is opposed by the omnipresent step-step interactions that tend to equilibrate the lengths of the two terraces adjacent to the step. In (a non-dimensionalized version of) LW the term that accounts realistically the role of the step-step repulsions is

Ri  2Fi(n1)  Fi 1(n1)  Fi 1(n1) , with

n1 Fi  n1  xi n1  xi1 and the canonical value of n, the power of the step-step distance d in the

ACS Paragon Plus Environment

Page 15 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Crystal Growth & Design

repulsions law U ~ 1 d n is 2 73, stemming mainly from elastic, at low temperatures, or entropic, at high temperature, sources. We confront the results from LW with an ad hoc, minimal model (MM)49 in which the first two terms (the destabilizing or positive feedback) are the same as in LW, plus a simpler negative, opposing the destabilization, feedback R  F ( n1) . Here we show results i

i

only for the one-sided versions of the two models, i.e. b=1, then the destabilizing part simplifies to xi and no parameter remains to enter eventually the time-scaling. The time-scaling of the bunch size N in LW and MM is presented in the Figure 7. It is the same as in vicCA, eq. (1) and also independent of n (see the supporting information for further discussion on this independence):

N

 1/ T 3

(6)

where a further conjecture is made -  is the order of the curve “surface slope vs. spatial coordinate” in the bunches, in LW =2 while in MM it is 1. Note that MM provides also a counterexample to the eventual statements of a triviality of the value   1 2 . We compare these results with data obtained when integrating the equations of step motion as deduced first by Sato and Uwaha51 for the case of growth destabilized by the presence of inverse (negative) EhrlichSchwoebel effect (Gr+iSE). The same monitoring scheme42 is used for all numerical studies.

Figure 6 Time-scaling of the bunch size in the one-sided version of LW (KL version of the phenomenon) and MM models, left panel. L0 is a rescaled value propoportional to the mean terrace width l0, defined in the Supporting Information. Larger L0 means a larger destabilization. For MM n = 2. In the right panel are shown results for the Sato and Uwaha model (Gr+iSE)53 for two regimes of the instability – kinetics-lmited (KL) and diffusion-limited (DL), time t is rescaled properly, as shown in the Supporting Information.

ACS Paragon Plus Environment

Crystal Growth & Design 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 23

Figure 8 Parallel results on the average bunch width as function of the rescaled time. It is the bunch width that distinguishes in between the values of the step-step repulsions exponent n, B  1 b  1 b  . The time-scaling of N is only one part of the whole scaling picture as drawn for first time in Ref. 40 and modified in Ref. 67 in the three ODE-system based models the SB phenomenon results in selfaffine surfaces and two length-scales are necessary to describe them thoroughly 43. Thus, while the bunch size N corresponds to the vertical length-scale (height), an additional, horizontal length-scale enters the problem – the bunch width W40,43. For LW we show in Figure 8 that: W ~T

n1 2( n1)

with the important exception that for n=2 the time scaling exponent of W is 1/5 instead of 1/6. The exponent in eq. (6), together with the exponent of N, β=1/2, corresponds to the prediction of Pimpinelli-Tonchev-Videcoq-Vladimirova (PTVV) 40 for the ρ=-1 universality class but with a shift k=2 in the values of n. With such a shift the exponents in the time-scaling of the bunch width predicted by PTVV for values of n’ are obtained from the LW model with values of n = n’+2 (note here again the exception of bunch width time-scaling exponent of 1/5 instead of 1/6 as predicted formally by PTVV for n’=0). Thus the program started by PTVV

40

and modified by Krug et al.

67

by introducing the correction k in order to distinguish between the DL and KL modes of the instability and transforming n’ into n’’=n’+k, with values of k being 0 (DL) and 1 (KL), should be modified further – our study shows that for the LW-model (KL) k=2 in general. For the Gr+iSE model the shifts from the PTVV prediction depend on whether the KL or DL version of the model is studied. For the KL version the shift is the same as for LW while for the DL version the shift is n = n’+1 hence k = 1 in the DL regime of the instability. 4 Conclusion It has been shown that our model, a combination of Cellular Automaton and Monte Carlo modules, significantly contributes towards building a unifying viewpoint on a class of phenomena,

ACS Paragon Plus Environment

Page 17 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Crystal Growth & Design

the surface instabilities, and, more specifically, on an archetypical instability – the step bunching (SB). By use of the single conceptual atomistic scale vicCA model with the decisive support of results from models on the next scale - based on systems of ODE for the velocity of each step from the vicinal crystal surface we address two deficiencies of the theoretical description of surface instabilities. First of them is the lack of unified treatment of the SB instability that would comprise both directions of the direct current or fluid flows above the surface, step-up and step-down. The second is the unfinished program of finding the full set of scaling exponents/relations quantifying the intermediate stages of the instability, started by A. Pimpinelli et al. 40 and partially continued by J. Krug et al.

46.

For each concrete case it is of special importance to decide what is the relevant

number of length-scales, one or two, towards quantifying the intermediate regime of the instability revealing the surface self-similarity/self-affinity. We provide an innovative version of stability analysis in a situation where analytical expressions are missing – based on the bunch size N as a measure and using small-size systems we probe the parameter space on a dense grid up to moderate times that are rescaled using the time scaling of N. Out of this analysis we obtain information organized into color charts for the four possible cases combining one of the vicinal processes, growth or sublimation, with one of the destabilizations, step-up or step-down adatom drift. These plots, possessing mirror symmetry for given destabilization, show where the system is most unstable. The value of the numerical pre-factor in the time scaling of N we find from a “classical”, Burton-Cabrera-Frank type of model – the most studied in the field model of Liu and Weeks (LW)72, taking the special case when no parameter remains in the model equation(s) – the so called one-sided model, then we integrate the ODE-system numerically. Important role plays further the parallel study of an ad hoc model with different feed-back mechanism to oppose to the step bunching/merger but with the same destabilization mechanism. We confirm this value with studies on the Sato and Uwaha model of vicinal growth destabilized by negative Ehrlich-Shwoebel effect (Gr+iSE), it permits studies of the two versions of the phenomenon – diffusion-limited and kineticslimited. In conclusion, we have shown that our atomistic-scale model vicCA is able to reproduce the instability caused by the two opposite drift directions in the two fundamental situations of step motion – vicinal sublimation and vicinal growth. By performing a cross-validation between the vicCA model as extended here, on the one side, and, on the other side, LW44,45,67,72 and Gr+iSE53 models, based on the time-scaling of N, the universality of intermediate stages of bunch dynamics is revealed. We show unambiguously that it is independent of the presence or absence of step-step repulsion and, further – independent also of the concrete form as determined by the power of the distance in the step-step repulsions law. The second length-scale, macrostep size for vicCA58 or bunch width for LW/Gr+iSE, is the one that brings the specific information on the presence and ACS Paragon Plus Environment

Crystal Growth & Design 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 23

type of step-step repulsions. Our study opens up a pathway for studying of a further extension of vicCA model in 2+1 dimensions, in order to quantify from a unified viewpoint the step bunching + step meandering instability. Acknowledgements The present research is supported by the National Science Centre (NCN) of Poland (Grant No. 2013/11/D/ST3/02700). VT, HP and AK acknowledge partial financial support from T02-8/121214 with the Bulgarian National Science Fund. Most of the calculations were done on HPC facility Nestum (BG161PO003-1.2.05). VT acknowledges the useful discussions and collaborations on the ODE based models with S. Stoyanov, C. Misbah, J. Krug, A. Pimpinelli, N. Akutsu, H. Omi, Y. Homma, Y. Saito, B. Ranguelov, D. Staneva and others through the years. Supporting Information. Supporting Information file is available. It contains - (i) dimensional analysis of the continuum equation of Pimpinelli-Tonchev-Videcoq-Vladimirova40 for the time evolution of height of the crystal surface resulting in a new prediction – that not only β but also the pre-factor in the time-scaling of the bunch size N (in the ρ=-1 universality class) is independent of the step-step repulsion exponent n; (ii) dimensional and further numerical analysis of the LW and Gr+iSE models illustrating in particular the role of the second length-scale – the bunch width and (iii) brief discussion on the universality classes in step bunching.

ACS Paragon Plus Environment

Page 19 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Crystal Growth & Design

References (1) Johnson, R. P., Construction of filament surfaces. Phys. Rev. 1938, 54, 459. (2) Chernov, A., The spiral growth of crystals. Phys.-Usp. 1961, 4, 116-148. (3) Schwoebel, R. L.; Shipsey, E. J., Step motion on crystal surfaces. J. Appl. Phys. 1966, 37, 3682-3686. (4) Mischgofsky, F. H.; Delhez, R., Optical and X-ray characterization of the layer type perovskites (CnH2n+1NH3)2CuCl4 in relation to their growth. J. Cryst. Growth 1978, 44, 145-156. (5) Zakurdaev, I. V., Change in the geometry and emission parameters of crystal surfaces at directed processes of mass transfer. Izv. Akad. Nauk SSSR Ser. Fiz. 1976, 40, 7. (6) Geguzin, Y. E.; Kaganovski, Y. S., Diffusionnie Processi na Poverkhnosti Kristallov, Moscow, Energoizdat: 1984. (7) Ding, J.; Zhao, Y.; Sun, Y.; Du, S.; Cui, H.; Jing, L.; Cheng, X.; Zuo, Z.; Zhan, X., Atomic force microscopy investigation of a step generation and bunching on the (100) facet of a CH3NH3PbI3 crystal, grown from γ‐Butyrolactone. Cryst. Res. Technol. 2017, 52. (8) Matsuoka, K.; Yagi, S.; Yaguchi, H., Growth of InN/GaN dots on 4H-SiC (0001) 4° off vicinal substrates by molecular beam epitaxy. J. Cryst. Growth 2017, 477, 201-206. (9) Bellmann, K.; Pohl, U. W.; Kuhn, C.; Wernicke, T.; Kneissl, M., Controlling the morphology transition between step-flow growth and step-bunching growth. J. Cryst. Growth 2017, 478, 187-192. (10) Hou, M.; Qin, Z.; Zhang, L.; Han, T.; Wang, M.; Xu, F.; Wang, X.; Yu, T.; Fang, Z.; Shen, B., Excitonic localization at macrostep edges in AlGaN/AlGaN multiple quantum wells. Superlattices Microstruct. 2017, 104, 397-401. (11) Syväjärvi, M.; Yakimova, R.; Janzén, E., Step-bunching in 6H-SiC growth by sublimation epitaxy. J. Phys.: Condens. Matter 1999, 11, 10019-10024. (12) Syväjärvi, M.; Yakimova, R.; Janzén, E., Step-bunching in SiC epitaxy: Anisotropy and influence of growth temperature. J. Cryst. Growth 2002, 236, 297-304. (13) Tabuchi, Y.; Ashida, K.; Sonoda, M.; Kaneko, T.; Ohtani, N.; Katsuno, M.; Sato, S.; Tsuge, H.; Fujimoto, T., Wide (0001) terrace formation due to step bunching on a vicinal 4H-SiC (0001) epitaxial layer surface. J. Appl. Phys. 2017, 122, 075702. (14) Onuma, A.; Maruyama, S.; Komatsu, N.; Mitani, T.; Kato, T.; Okumura, H.; Matsumoto, Y., Quantitative Analysis of Nanoscale Step Dynamics in High-Temperature Solution-Grown Single Crystal 4H-SiC via In Situ Confocal Laser Scanning Microscope. Cryst. Growth Des. 2017, 17, 2844-2851. (15) Bao, J.; Yasui, O.; Norimatsu, W.; Matsuda, K.; Kusunoki, M., Sequential control of stepbunching during graphene growth on SiC (0001). Appl. Phys. Lett. 2016, 109, 081602. (16) Yi, D.; Luo, D.; Wang, Z.-J.; Dong, J.; Zhang, X.; Willinger, M.-G.; Ruoff, R. S.; Ding, F., What drives metal-surface step-bunching in graphene chemical vapor deposition? Phys. Rev. Lett. 2018, 120, 246101. (17) Pollinger, F.; Schmitt, S.; Sander, D.; Tian, Z.; Kirschner, J.; Vrdoljak, P.; Stadler, C.; Maier, F.; Marchetto, H.; Schmidt, T., Nanoscale patterning, macroscopic reconstruction, and enhanced surface stress by organic adsorption on vicinal surfaces. New J. Phys. 2017, 19, 013019. (18) Sullivan, D. J., Quinolines block every step of malaria heme crystal growth. Proc. Natl. Acad. Sci. USA. 2017, 114, 7483-7485. (19) Gura, A.; Bertino, G.; Bein, B.; Dawber, M., Transition regime from step-flow to stepbunching in the growth of epitaxial SrRuO3 on (001) SrTiO3. Appl. Phys. Lett. 2018, 112, 182902. (20) Booth, N. A.; Chernov, A. A.; Vekilov, P. G., Step bunching in potassium dihydrogen phosphate crystal growth: Phenomenology. J. Mater. Res. 2002, 17, 2059-2065. (21) Booth, N. A.; Chernov, A. A.; Vekilov, P. G., Characteristic lengthscales of step bunching in KDP crystal growth: In situ differential phase-shifting interferometry study. J. Cryst. Growth 2002, 237-239, 1818-1824. (22) Bredikhin, V.; Galushkina, G.; Kulagin, A.; Kuznetsov, S.; Malshakova, O., Competing ACS Paragon Plus Environment

Crystal Growth & Design 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 23

growth centers and step bunching in KDP crystal growth from solutions. J. Cryst. Growth 2003, 259, 309-320. (23) Li, W. D.; Wang, D. L.; Gu, Q. T.; Wang, S. L. In Surface Micromorphology and Growth Rate Related to Supersaturation of KDP Crystals Grown at Different Temperature, Mater. Sci. Forum, 2017; Trans Tech Publ: 2017; pp 111-119. (24) Gliko, O.; Booth, N. A.; Vekilov, P. G., Step bunching in a diffusion-controlled system: Phase-shifting interferometry investigation of ferritin. Acta Crystallogr., Sect. D: Biol. Crystallogr. 2002, 58, 1622-1627. (25) Can Zhu, Shunta Harada, Kazuaki Seki, Huayu Zhang, Hiromasa Niinomi, Miho Tagawa, and Toru Ujihara, Influence of Solution Flow on Step Bunching in Solution Growth of SiC Crystals Cryst. Growth Des. 2013, 13, 3691−3696 (26) Ariyawong, K.; Shin, Y. J.; Dedulle, J.-M.; Chaussende, D., Analysis of macrostep formation during top seeded solution growth of 4H-SiC.Cryst. Growth Des. 2016, 16, 3231-3236. (27) Filimonov, S. N.; Hervieu, Y. Y., Model of step propagation and step bunching at the sidewalls of nanowires. J. Cryst. Growth 2015, 427, 60-66. (28) Latyshev, A. V.; Aseev, A. L.; Krasilnikov, A. B.; Stenin, S. I., Transformations on clean Si(111) stepped surface during sublimation. Surf. Sci. 1989, 213, 157-169. (29) Burton, W.-K.; Cabrera, N.; Frank, F., The growth of crystals and the equilibrium structure of their surfaces. Philos. Trans. R. Soc., A 1951, 243, 299-358. (31) Stoyanov, S., Heating Current Induced Conversion between 2× 1 and 1× 2 Domains at Vicinal (001) Si Surfaces–Can it be Explained by Electromigration of Si Adatoms? Jpn. J. Appl. Phys. 1990, 29, L659. (31) Stoyanov, S., Electromigration induced step bunching on si surfaces-how does it depend on the temperature and heating current direction?. Jpn. J. Appl. Phys. 1991, 30, 1-6. (32) Stoyanov, S. S.; Nakahara, H.; Ichikawa, M., Dynamics of step bunching induced by dc resistive heating of Si wafer. Jpn. J. Appl. Phys. 1994, 33, 254. (33) Zheng, L. X.; Xie, M. H.; Xu, S. J.; Cheung, S. H.; Tong, S. Y., Current-induced migration of surface adatoms during GaN growth by molecular beam epitaxy. J. Cryst. Growth 2001, 227228, 376-380. (34) De Theije, F. K.; Schermer, J. J.; Van Enckevort, W. J. P., Effects of nitrogen impurities on the CVD growth of diamond: Step bunching in theory and experiment. Diamond Relat. Mater. 2000, 9, 1439-1449. (35) Golovin, A. A.; Davis, S. H.; Voorhees, P. W., Step bunching in the absence of an EhrlichSchwoebel barrier during nanowire growth. J. Appl. Phys. 2010, 107, 024308. (36) Fujita, K.; Ichikawa, M.; Stoyanov, S. S., Size-scaling exponents of current-induced step bunching on silicon surfaces. Phys. Rev. B 1999, 60, 16006-16012. (37) Xie, M. H.; Cheung, S. H.; Zheng, L. X.; Ng, Y. F.; Wu, H.; Ohtani, N.; Tong, S. Y., Step bunching of vicinal GaN(0001) surfaces during molecular beam epitaxy. Phys. Rev. B 2000, 61, 9983-9985. (38) Uwaha, M.; Sato, M., Theoretical study of step bunching and step wandering induced by drift of adatoms. J. Vac. Soc. Jpn. 2002, 45, 557-563. (39) Sato, M.; Uwaha, M.; Saito, Y., Step bunching induced by drift of adatoms with anisotropic surface diffusion. J. Cryst. Growth 2002, 237-239, 43-46. (40) Pimpinelli, A.; Tonchev, V.; Videcoq, A.; Vladimirova, M., Scaling and universality of selforganized patterns on unstable vicinal surfaces. Phys. Rev. Lett. 2002, 88, 2061031-2061034. (41) Stoyanov, S.; Métois, J. J.; Tonchev, V., Current induced bunches of steps on the Si(111) surface - a key to measuring the temperature dependence of the step interaction coefficient. Surf. Sci. 2000, 465, 227-242. (42) Tonchev, V.; Ranguelov, B.; Omi, H.; Pimpinelli, A., Scaling and universality in models of step bunching: The "c +-C -" model. Eur. Phys. J. B 2010, 73, 539-546. (43) Tonchev, V., Classification of step bunching phenomena. Bulg. Chem. Commun. 2012, 44, 124-130. ACS Paragon Plus Environment

Page 21 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Crystal Growth & Design

(44) Popkov, V.; Krug, J., Shape and scaling of moving step bunches. Europhys. Lett. 2005, 72, 1025-1031. (45) Popkov, V.; Krug, J., Dynamic phase transitions in electromigration-induced step bunching. Phys. Rev. B 2006, 73, 235430. (46) Krug, J., Nonlinear Dynamics of Surface Steps. In Nonlinear Dynamics of Nanosystems, 2010; pp 143-164. (47) Pierre-Louis, O.; Danker, G.; Chang, J.; Kassner, K.; Misbah, C., Nonlinear dynamics of vicinal surfaces. J. Cryst. Growth 2005, 275, 56-64. (48) Chang, J.; Pierre-Louis, O.; Misbah, C., Birth and morphological evolution of step bunches under electromigration. Phys. Rev. Lett. 2006, 96, 195901. (49) Ranguelov, B.; Tonchev, V.; Misbah, C., Step Bunching in Conserved Systems: Scaling and Universality. In Nanoscience and Nanotechnology, Balabanova, E.; Dragieva, I., Eds. Heron Press: Sofia, 2006; Vol. 6. (50) Sato, M.; Uwaha, M., Pattern formation in the instability of a vicinal surface by the drift of adatoms. Phys. Rev. E 1999, 60, 7120-7125. (51) Sato, M.; Uwaha, M., Growth of step bunches formed by the drift of adatoms. Surf. Sci. 1999, 442, 318-328. (52) Sato, M.; Uwaha, M.; Saito, Y., Instabilities of steps induced by the drift of adatoms and effect of the step permeability. Phys. Rev. B 2000, 62, 8452-8472. (53) Sato, M.; Uwaha, M., Growth law of step bunches induced by the Ehrlich-Schwoebel effect in growth. Surf. Sci. 2001, 493, 494-498. (54) Załuska-Kotur, M.; Krzyzewski, F.; Krukowski, S., Emergence of regular meandered step structure in simulated growth of GaN(0001) surface. J. Cryst. Growth 2012, 343, 138-144. (55) Załuska‐Kotur, M. A.; Krzyżewski, F., Step bunching process induced by the flow of steps at the sublimated crystal surface. J. Appl. Phys. 2012, 111, 114311. (54) Krzyzewski, F.; Załuska-Kotur, M. A., Coexistence of bunching and meandering instability in simulated growth of 4H-SiC(0001) surface. J. Appl. Phys. 2014, 115, 213517. (57) Krasteva, A.; Popova, H.; Krzyżewski, F.; Załuska-Kotur, M.; Tonchev, V. In Unstable vicinal crystal growth from cellular automata, AIP Conf. Proc. 2016, 220014. (58) Krzyżewski, F.; Załuska-Kotur, M.; Krasteva, A.; Popova, H.; Tonchev, V., Step bunching and macrostep formation in 1D atomistic scale model of unstable vicinal crystal growth. J. Cryst. Growth 2017, 474, 135-139. (59) Toktarbaiuly, O.; Usov, V.; Ó Coileáin, C.; Siewierska, K.; Krasnikov, S.; Norton, E.; Bozhko, S. I.; Semenov, V. N.; Chaika, A. N.; Murphy, B. E.; Lübben, O.; Krzyżewski, F.; Załuska-Kotur, M. A.; Krasteva, A.; Popova, H.; Tonchev, V.; Shvets, I. V., Step bunching with both directions of the current: Vicinal W(110) surfaces versus atomistic-scale model. Phys. Rev. B 2018, 97, 035436. (60) Toktarbaiuly, O. Electromigration of atoms on vicinal surfaces at high temperatures, PhD thesis, University of Dublin, Trinity College Dublin, Dublin, 2015. (61) Lutsko, J. F.; Van Driessche, A. E. S.; Durán-Olivencia, M. A.; Maes, D.; Sleutel, M., Step Crowding Effects Dampen the Stochasticity of Crystal Growth Kinetics. Phys. Rev. Lett. 2016, 116, 015501. (62) Stoyanov, S., New type of step bunching instability at vicinal surfaces in crystal evaporation affected by electromigration. Surf. Sci. 1998, 416, 200-213. (63) Stoyanov, S.; Tonchev, V., Properties and dynamic interaction of step density waves at a crystal surface during electromigration affected sublimation. Phys. Rev. B 1998, 58, 1590-1600. (64) Ranguelov, B.; Stoyanov, S., Instability of vicinal crystal surfaces with transparent steps: Transient kinetics and non-local electromigration. Surf. Sci. 2009, 603, 2907-2911. (65) Ranguelov, B.; Altman, M. S.; Markov, I., Critical terrace width for step flow growth: Effect of attachment-detachment asymmetry and step permeability. Phys. Rev. B 2007, 75, 245419. (66) Ranguelov, B. S.; Markov, I. V., Adatom diffusion on vicinal surfaces with permeable steps. Cent. Eur. J. Phys. 2009, 7, 350-355. ACS Paragon Plus Environment

Crystal Growth & Design 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 23

(67) Krug, J.; Tonchev, V.; Stoyanov, S.; Pimpinelli, A., Scaling properties of step bunches induced by sublimation and related mechanisms. Phys. Rev. B 2005, 71, 045412. (68) Ohtani, N.; Katsuno, M.; Aigo, T.; Fujimoto, T.; Tsuge, H.; Yashiro, H.; Kanaya, M., Step bunching behaviour on the {0001} surface of hexagonal SiC. J. Cryst. Growth 2000, 210, 613-622. (69) Néel, N.; Maroutian, T.; Douillard, L.; Ernst, H. J., From meandering to faceting, is step flow growth ever stable? Phys. Rev. Lett. 2003, 91, 226103. (70) Omi, H.; Homma, Y.; Tonchev, V.; Pimpinelli, A., New types of unstable step-flow growth on Si(111)-(7×7) during molecular beam epitaxy: Scaling and universality. Phys. Rev. Lett. 2005, 95, 216101. (71) Merton, R. K., The Matthew effect in science: The reward and communication systems of science are considered. Science 1968, 159, 56-63. (72) Liu, D. J.; Weeks, J. D., Quantitative theory of current-induced step bunching on Si(111). Phys. Rev. B 1998, 57, 14891-14900. (73) Nozières, P., Shape and growth of crystals. In Solids far from equilibrium, Godrèche, C., Ed. Cambridge University Press: 1992; Vol. 56, pp 1-154.

ACS Paragon Plus Environment

Page 23 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Crystal Growth & Design

For Table of Contents Use Only

Scaling and dynamic stability of model vicinal surfaces Filip Krzyżewski, Magdalena Załuska-Kotur, Anna Krasteva, Hristina Popova and Vesselin Tonchev

Synopsis Based on a novel model we analyze the step bunching instability on vicinal crystal surfaces destabilized by step-down and step-up currents in sublimation and growth. After finding a self-similar solution - the time-scaling of the bunch size N, N  2 T 3 , we use N further as a measure to probe the model’s stability on a grid of points in the parameter space.

ACS Paragon Plus Environment