Relationship Between Macrostep Height and Surface Velocity for a Reaction-Limited Crystal Growth Process

This work examined the effect of macrostep height on the growth velocity of a vicinal surface during reaction- (interface-) limited crystal growth under non-equilibrium steady state conditions. The Monte Carlo method was employed, based on a restricted solid-on-solid (RSOS) model with point-contact-type step-step attraction (termed the p-RSOS model). Although this is a simple lattice model, the model surface shows a variety of distinctive configurations depending on the temperature and the driving force for crystal growth. The results demonstrate that the surface velocity decreases as the height of the faceted macrostep increases. In addition, the significant variation in surface velocity recently reported by Onuma {\it et al.} in a study based on 4H-SiC was reproduced. This work also shows that the terrace slope, elementary step velocity and elementary step kinetic coefficient are all affected by the faceted macrostep height.


I. INTRODUCTION
The production of high-quality SiC crystals is an important prerequisite for the production of advanced power devices with low power consumption rates. Although self-organized faceted macrosteps are known to lower the quality of crystalline SiC [1], dislocations penetrating the crystal have been shown to end at the side surfaces of macrosteps [2][3][4]. Hence, the intentional introduction of macrosteps can effectively decrease the dislocation density in a SiC crystal. In addition, controlling the self-assembly and disassembly of faceted macrosteps is an important aspect of fabricating semiconductor crystals.
As an example, Chernov [10] and co-workers [11,12] performed detailed studies of the diffusion field around macrosteps or trains of macrosteps during steady-state growth in which volume diffusion was the rate-determining process. This same prior work also generated mathematical expressions for the velocity of macrostep advancement V m under these conditions, i. e. V m ∝ 1/h m where h m is the height of the macrostep [10]. The essential point to derive this equation is the mass conservation during the diffusion limited crystal growth. Recently, Onuma et al. [13] determined the rela- * nori3@phys.osakac.ac.jp, nori@phys.osakac.ac.jp tionship between the macrostep velocities and heights in 4H-SiC using in situ confocal laser scanning microscopy and found that, when employing a Si-Ni flux, the experimental data could be explained by Chernov's equation [10]. Thus, they concluded that the rate-determining process in 4H-SiC crystal growth is volume diffusion. In contrast, when Al was added to the Si-Ni flux, the data were too scattered to be explained by the equation, and so the rate was thought to have been determined by the interfacial reactions.
Despite this prior work, reaction-(that is, interface-) limited growth has not been studied sufficiently. In addition, such interfacial (or surface) reactions are always accompanied by volume diffusion or surface diffusion, and so it is difficult to separate the effects of interfacial reactions from those of diffusion in an actual system. For this reason, the present work employed computational experiments involving extreme conditions, using the Monte Carlo method. We have been studying the self-assembly/disassembly of macrosteps on a theoretical basis, neglecting the effects of the surface and volume diffusion. The main driving force for the assembly of steps in such cases is anomalous surface tension at low temperatures. The anomaly in surface tension is caused microscopically by step-step attraction. This prior work employed the restricted solidon-solid (RSOS) model in conjunction with point contact step-step attraction, termed the p-RSOS model [14][15][16][17][18][19]. Point contact step-step attraction is considered to result from the energy gain associated with the point at which neighboring steps meet. At such meeting points, the dan- gling bonds of the neighboring steps overlap to merge to create actual bonds between atoms. The resulting energy gain leads to the step-step attraction. Although the p-RSOS model involves a simple lattice, the vicinal surface of the model exhibits a variety of surface configurations with respect to the self-assembling/disassembling of macrosteps depending on the temperature and the driving force for crystal growth [20][21][22][23].
The most important aspect of the p-RSOS model is that it allows reliable calculations of a polar graph of the surface tension (the Wulff figure) and the equilibrium crystal shape (ECS). The surface tension is surface free energy per unit normal area. Since surfaces and steps are low-dimensional objects, thermal fluctuations are so severe that they destroy the ordered phase for the system with the short-range force, generally [24]. The calculation of the surface free energy with the mean-field or the quasi-chemical approximation often leads to wrong results. Hence, the calculations more precise than the mean field approximation are required. In the present work, we applied the transfer matrix version of the density matrix renormalization group (DMRG) method [25][26][27][28][29][30][31][32] to calculate the surface tension. Therefore, the morphological phenomena resulting from the anisotropy of the surface tension can be analyzed by referring to the calculated Wulff figures.
The aim of this work was to demonstrate the effect of the macrostep height on the surface velocity, the terrace slope, the elementary step velocity, and the elementary step kinetic coefficient in the case of reaction (interface)-limited crystal growth in the non-equilibrium steady state. This involved simulating the vicinal surface based on the p-RSOS model and using the Monte Carlo method in the non-conserved system of mass.

II. MICROSCOPIC MODEL AND CALCULATION METHOD
A. The p-RSOS model In this study, the p-RSOS model was employed as the microscopic model ( Fig. 1) [14]- [23]. In this approach, the total energy of the (001) surface can be written as where N is the total number of lattice points, ǫ surf is the surface energy per unit cell on the planar (001) surface, ǫ is the microscopic ledge energy, δ(a, b) is the Kronecker delta and ǫ int is the microscopic step-step interaction energy. ǫ int contributes to the surface energy at the meeting point of neighboring steps when the height difference in the diagonal direction is ±2, and is assumed to originate from the energy gain when overlapping dangling bonds at step edges form bonding states. The summation with respect to (n, m) is performed over all sites on the square lattice and when ǫ int is negative, the step-step interaction becomes attractive.
It should also be noted that the RSOS model implicitly requires that the height difference between the nearest neighbor sites be restricted to {0, ±1}. The surface free energy density, surface tension and ECS can be calculated using the transfer matrix version of the DMRG method [25][26][27], termed the product-wavefunction renormalization group (PWFRG, or tensor network) method [28][29][30][31]35].
After detailed calculations of the surface free energy density and the surface tension, we found out that the surface undergoes the first order transition with respect to the self-assembling/disassembling of the elementary steps in the vicinal surface ( Fig. 2). At high temperatures, the vicinal surface exhibits the Gruber-Mullins-Pokrovsky-Talapov (GMPT) universal behavior [31,[33][34][35] (Fig. 2 (c)). As the temperature is decreased, the surface tension around the (111) surface becomes discontinuous ( Fig. 2 (b)) below T f,1 , while the surface tension around the (001) surface becomes discontinuous ( Fig. 2  (a)) below T f,2 . The faceting diagram shown in Fig. 3 summarizes the surface tension discontinuities.
The profile of the faceted macrostep on the vicinal surface at equilibrium is also determined by the surface tension continuity [21]. At T > T f,1 (GMPT zone), a faceted macrostep does not appear on the vicinal surface. In contrast, at T f,1 ≥ T > T f,2 (that is, the step-droplet zone), a faceted macrostep with the side surface being the (111) surface coexists with vicinal surfaces having the slope p 1 . Note that vicinal surfaces with the slope p 1 < p(111) do not appear, because they are thermodynamically unstable. At T < T f,2 (that is, the step-faceting zone), the vicinal surface consists of only the (001) and the (111) surfaces.
Step-faceting (QI Bose solid +Vacuum) Step droplet II ( [20]. Red triangles indicate T f,2 values. Below T f,2 , the polar graph of the surface tension is discontinuous around the (001) surface [16,17,20]. Purple squares indicate T f,1 values. Below T f,1 , the polar graph of the surface tension is discontinuous around the (111) surface [16,17,20]. Pink circles indicate the roughening transition temperatures of the (001) surface. All symbols show values calculated using the PWFRG method. As for the "QI Bose solid", "QI Bose liquid", and "AI Bose gas", please refer to the Ref. [20]. This figure is taken from [20].

B. Monte Carlo method
The Monte Carlo method was adopted together with the Metropolis algorithm for the non-conserved systems to study the non-equilibrium steady state. The microscopic surface energy for a fixed number of steps, N step , was determined using the equation where t is the time in units of Monte Carlo steps per site (MCS/site) and ∆µ is the driving force for the crystal growth (that is, the chemical potential difference between the bulk crystal and the ambient phase). When ∆µ > 0, the crystal grows, whereas when ∆µ < 0, the crystal recedes (evaporates, dissociates or melts). The site on the surface for an Monte Carlo event is chosen randomly. Addition or removal of an growth unit is determined with the probability 0.5. The energies of the surface before or after of the event E i or E f are calculated based on Eq. (2).
Periodic boundary conditions were required in the direction of the mean step-running. In the direction normal to the mean step-running, the surface heights of the lowest side line were connected to the surface heights at the topmost line by adding the number of elementary steps.
The initial surface configuration involved the preparation of a parallel train of steps and a single macrostep.
It should be noted that the external parameters are the microscopic ledge energy ǫ, microscopic step-step attraction ǫ int (< 0), temperature T , total number of elementary steps N step , linear system size L, and driving force for the crystal growth ∆µ.

C. Macrostep height
Images of the simulated surfaces at 4 × 10 8 MCS/site are presented in Fig. 4. The average macrostep height was obtained using the equation wherex is the 110 direction (normal to the mean steprunning direction),ỹ is the 1 10 direction (along the mean step-running direction) and n step is the number of merged steps. When determining the average value, the data over the first 2 × 10 8 MCS/site were discarded and the values over the subsequent 2 × 10 8 MCS/site were averaged. The surface velocity V was estimated by finding the average surface heighth(t), using the formulā Here, V is defined as where t 0 and t max are 2 × 10 8 MCS/site and 4 × 10 8 MCS/site, respectively.

A. Macrostep height dependence of the surface velocity
The effects of height n on the surface velocity are shown in Fig. 5. These data demonstrate that surface velocity V decreases as the mean size of the macrostep, n , increases. It should be noted that both V and n are functions of driving force for the crystal growth ∆µ in the non-equilibrium steady state [23]. According to our previous work [23], there are several characteristic driving forces, which are summarized in Table I. In the step-droplet zone, for |∆µ| < ∆µ R , the vicinal surface has the structure shown in Fig. 4 (a), i. e. a faceted macrostep with a (111) side surface in conjunction "terrace" surfaces having the slope p 1 . The characteristics of this structure are provided in Fig. 6 (c). In contrast, in the case of |∆µ| > ∆µ R , the vicinal surface is kinetically roughened and a macrostep is not formed. The characteristics of this structure are shown in Fig.  6 (d). Although there is no macrostep, the elementary steps merge locally to form the faceted local macrosteps. In addition, n and V exhibit power law behaviors [22].
In the step-faceting zone, the change in dynamics is more complex than in the step-droplet zone, as there are characteristic driving forces near equilibrium in addition to ∆µ R . At equilibrium, only the (001) and (111) surfaces are thermodynamically stable, and the (001) surface forms the terrace, while the (111) surface forms the side surface of the macrostep. In the event that |∆µ| < ∆µ f , the surface does not grow, because there is a lack of nucleation growth, meaning that the mean waiting time for the formation of a single nucleus exceeds the observation time (4 × 10 8 MCS/site in this work). In the case that a 2D island is formed at the edge of the macrostep as a result of thermal fluctuations, the island shrinks because the elementary step that surrounds the island recedes back to the macrostep due to the step tension ( Fig. 6 (a)). In addition, since the vicinal surface freezes for |∆µ| < ∆µ f , the height of the faceted macrostep is largely determined by the initial configuration. As a result, the relationship between n and V also depends on the initial configuration, and so the n -V plot in Fig. 7 exhibits significant scatter in the region between the solid lines.
In the case of ∆µ f < |∆µ| < ∆µ co , the crystal growth proceeds via 2D single nucleation (Fig. 6 (b)). If ∆µ > 0, then the crystal grows and nucleation occurs at the lower edge of the macrostep while growth occurs on both the (001) and (111) surfaces ( Fig. 8 (a)). Conversely, if ∆µ < 0, then the crystal evaporates or dissolves. A "negative" 2D island, which is the cluster of holes on a surface, (Fig. 4 (c)) is formed as a nucleus at the upper edge of the macrostep, and this negative island increases in size (Fig. 8 (b)). If ∆µ co < |∆µ| < ∆µ R , then the crystal grows via 2D multi nucleation (Fig. 6 (c)). If ∆µ R ≤ |∆µ|, then the surface is kinetically roughened, and both V and n demonstrate power law behavior with respect to |∆µ|. Although the relationship between V and n in the stepfaceting zone is similar to that in the step-droplet zone, the details of the relationships are different. To understand this phenomenon on the microscopic scale, these Monte Carlo results are further analyzed in the following subsections.

B.
Step attachment/detachment model In the case that two surfaces coexist, elementary steps will attach to or detach from the edge of the faceted side surface. The effect of ∆µ on the vicinal surface morphology for |∆µ| < ∆µ R was assessed by considering a step-attachment-detachment model as n evolved over time [22]. This was calculated as where n + is the rate at which the elementary steps catch up to a macrostep and n − is the rate at which the elementary steps detach from a macrostep. When n + < n − , the macrostep dissociates (as is the case for ∆µ R < |∆µ|), whereas when n + > n − (the case for ∆µ f < |∆µ| < ∆µ co ), n increases up to N step . In this case, n − limits the growth/recession rate of the surface. Under steady-state conditions, n + = n − = V /a, where a is the height of the elementary step. Thus, if |∆µ| < ∆µ R , then n + can be expressed as n + = ρ 1 v 1 , where ρ 1 and v 1 are the step density and the step velocity of the elementary steps, respectively. Based on the step-attachment-detachment model, surface velocity V can be expressed using more microscopic quantities, as [22,23] where p 1 is the slope of the terrace surface and v 1 is the mean step velocity of the elementary steps. Since total number of elementary steps N step is conserved, the slope p 1 can be expressed using n as wherep is the mean slope of the vicinal surface, N m is the number of macrosteps and a is the height of the elementary steps. It should be noted that p 1 = 0 for n = N step /N m = n * .
In the case of n ∼ n * , p 1 is approximately { √ 2p/( √ 2 − p)}(1/ n − 1/n * ). In addition, if n << n * , then 1/n * becomes negligible, such that Mean step velocity v 1 can be obtained on the basis of V and n together with Eqs. (7) and (8) as The kinetic coefficient k = v 1 /∆µ is also calculated using V and n , as Because n depends on ∆µ in the non-equilibrium steady state, a large n value is associated with |∆µ| < ∆µ y ≈ ∆µ f and p 1 ≈ 0. This result indicates that the v 1 and k terms in Eqs. (10) and (11) are divergent near equilibrium (Eq. (9)). However, V is also affected by ∆µ and decreases to zero near equilibrium. Whether v 1 and k are divergent or not near equilibrium depends on the rate of decreases of V and p 1 .

C. Step-droplet zone
In the step-droplet zone, the vicinal surface for |∆µ| < ∆µ R has a configuration similar to that in Fig. 6 (c). It is interesting to observe that the macrostep moves in the opposite direction compared to that of the elementary steps ( Fig. 9 (a)). This occurs because the side of the faceted macrostep is a (111) smooth surface, such that the number of kinks on this surface is extremely low. In addition, the vicinal surface grows via the step detachment/attachment mode, and so the effect of n on the surface velocity V is the same as the effect on the macrostep velocity V m .
In contrast, in the case of step bunching caused by diffusion-limited growth, the side surface of the bunched step is rough and there are many kinks on the side surface. Consequently, the center of the bunched step side surface moves in the same direction as the elementary steps ( Fig. 9 (b)). These results show that the ratelimiting process can be determined if the growth directions of the macrostep and the elementary steps are identified. Specifically, when the growth direction of the macrostep is the same as that of the elementary steps, the crystal growth is diffusion-limited, while if the opposite is true, then the crystal growth is reaction (interface)limited. The slope value, p 1 , can be calculated using Eq. (8) based on the n value together with N m = 2 andp = √ 2/2. The resulting p 1 values are shown in Fig. 10 (a), which demonstrates that p 1 decreases as n increases. From Fig. 10 (b), we can see that p 1 can also be expressed as  Fig. 10 (c). These data demonstrate that the absolute value of v 1 decreases as n increases. It is interesting to note that the kinetic coefficient of the elementary steps, k, calculated using Eq. (11) (Fig. 10 (d)) is almost constant (k ≈ .0.132) with respect to n , i. e. v 1 = 0.132∆µ. Since p 1 is also correlated with |∆µ| according to the relationship [22] where a = 0.332, b = 15.6, and c = 4.43 × 10 3 , we can solve Eq. (13) with respect to |∆µ|/ǫ to obtain Substituting Eq. (12) into Eq. (14), we obtain The values of v 1 calculated from Eq. (15) are shown by the solid lines in Fig. 10 (c), and these lines accurately reproduce the Monte Carlo data for |∆µ| < ∆µ R . Therefore, the surface velocity can be expressed as |V | = (kǫ)p 1 |∆µ|/ǫ = 0.132 × (0.674 − 0.299 ln n ) −1.76 × 10 −3 + 8.03 × 10 −5 − 6.75 × 10 −5 ln n The lines obtained from Eq. (16) are indicated in Fig.  5 (a) by solid lines, and also accurately reproduce the Monte Carlo data. In addition, for |∆µ| > ∆µ R (in which case only local macrosteps are present (Fig. 6 (d))), and using L = 240 √ 2, we obtain n − n ∞ = 9.16 × 10 −6 (|∆µ|/ǫ) −ζ , ζ = 1.89 ± 0.07, (17) |V | = 0.0855(∆µ)(|∆µ|/ǫ) β , β = 1.06 ± 0.06.
Therefore, we have The lines generated using Eqs. (19) and (20) are shown by the broken lines in Fig. 5 (a) and Fig. 10 (c), respectively. Although a large macrostep does not form, step velocity v 1 is lower than that produced by the original RSOS model because the locally merged steps pin the motion of the elementary steps. As shown in Fig. 6 (d), the elementary steps form a network of steps and, at the meeting points of neighboring steps, locally merged steps form local faceted macrosteps with finite lifetimes. Because the side surface of the local faceted macrostep has few kinks, the growth of this region is quite slow [17,18]. In this manner, the local faceted macrostep pins the elementary steps.

D. Step-faceting zone
In the case of the step-faceting zone, the vicinal surface consists of (001) terrace surfaces and the (111) side surface (Fig. 8)  In this subsection, we analyze the data shown in Fig. 6 (c). If ∆µ co < |∆µ| < ∆µ R , growth occurs via step-detachment/attachment, as in (Fig. 6  (c)). Since multi-nucleation occurs on both the (001) and (111) surfaces at the edge of the macrostep, the velocity of the macrostep is almost zero, and surface velocity V does not coincide with that of the macrostep.
In Fig. 11 (a), the p 1 values calculated from n data using Eq. (8) in conjunction with the Monte Carlo method are indicated by symbols, along with the p 1 values calculated using Eq. (8). These results demonstrates that p 1 falls to zero at n * . The line p 1 ∼ 1/ n is shown as a light green solid line in Fig. 11 (b). This line reproduces the symbols accurately except for the data close to n ∼ n * . In Figs. 11 (c), the v 1 values obtained by the Monte Carlo method in conjunction with Eq. (10) by using n and V values are indicated by symbols. In Figs. 11 (d), the k values obtained by the Monte Carlo method in conjunction with Eq. (11) by using n and V values are indicated by symbols.
The solid line in Fig. 11 (d) [23] represents where a ′ = 0.094, b ′ = 3.2 × 10 −3 , and c ′ = 0.18. |∆µ| is related to n [23] in association with p 1 according to and Here, ∆µ y (L) is one of the characteristic driving forces (Table I), g * P is the total step free energy of the 2D critical island at the macrostep edge (calculated using the 2D Ising model) times |∆µ/ǫ|, and c p is a coefficient related to the Zeldovich factor for the 2D nucleation. The solid line obtained from Eq. (21) reproduces the Monte Carlo results well, except in the case that n is close to n * . Here, k increases as n increases due to the decreased step density on the terrace, i. e. p 1 → 0. In the case of a high step density on the terrace, the number of locally merged steps is also high. Since each locally merged step is faceted and has few kinks, the velocity of the locally faceted macrostep is substantially lower than that of an elementary step. As a result, the motions of steps are pinned randomly by the locally faceted macrosteps (Fig. 6 (c)). Conversely, when n is large, |∆µ| is small, the step density on the terrace is nearly zero, and the elementary steps meet neighboring steps less frequently. Thus, the elementary steps travel like free steps and are not pinned as often.
It is interesting to observe that v 1 is almost constant except near n /N step ≈ 0 and n /N step ≈ n * /N step . The solid line in Fig. 11 (c) was calculated using the equation v 1 = k∆µ/ǫ, where k is that in Eq. (21) and ∆µ is related to the value of n by Eqs. (22) and (23). As n increases, k increases but |∆µ| decreases, and so v 1 is almost unchanged over the approximate range of ∆µ co < |∆µ| < ∆µ R . Because p 1 ∝ 1/ n , surface velocity V becomes over the approximate range of ∆µ co < |∆µ| < ∆µ R , except for n ∼ n * . At n n * , where |∆µ| ∆µ f , there is one faceted macrostep and several elementary steps traveling on the terrace surface. In this case, because V approaches zero more rapidly than does p 1 (due to the finite size effect), the k values obtained from Eq. (11) are reduced to zero.
Similar to the results obtained for the step-droplet zone, the local macrosteps pin the step motion of the elementary steps. Hence, the k and v 1 values are both substantially smaller than those in the original RSOS model.

IV. CONCLUSIONS
• The surface velocity V decreases as the height of a faceted macrostep increases during reaction-(interface-) limited crystal growth, and this re-duced velocity is related to a decrease in surface roughness or kink density.
• In the step-droplet zone, the faceted macrostep grows (recedes) in the opposite direction to the growth (recession) direction of the elementary steps. The decrease in the velocity of the faceted macrostep V m is also proportional to ln n (Eq. (16), |V m | ∝ |V |).
• In the step-droplet zone, the kinetic coefficient of the elementary steps is almost constant. This coefficient has a value approximately one-third that of the kinetic coefficient of the original restricted solid-on-solid (RSOS) model because of the local faceted macrosteps on the terrace.
• In the step-faceting zone, the absolute value of the growth velocity of the surface, |V |, varies widely within the range defined by |V | = constant ×1/ n , because the height of the macrostep depends on the initial configuration of the vicinal surface.
• In the step-faceting zone, the surface grows (re-cedes) via 2D multi-nucleation on both the (001) and (111) surfaces at the lower (upper) edge of the faceted macrostep. This nucleation on both surfaces inhibits the advancement of the faceted macrostep (|V m | ∼ 0).

V. ACKNOWLEDGMENTS
The author wishes to acknowledge helpful discussions with Prof. Y. Matsumoto, Prof. T. Nishinaga, Prof. T. Ujihara, Prof. K. Fujiwara, Prof. S. Uda and Prof. K. Kakimoto. This work is supported by Kakenhi Grants-in-Aid (Nos. JP25400413 and JP17K05503) from the Japan Society for the Promotion of Science (JSPS).