Disassembly of Faceted Macrosteps in the Step Droplet Zone in Non-Equilibrium Steady State

A Wulff figure—the polar graph of the surface tension of a crystal—with a discontinuity was calculated by applying the density matrix renormalization group method to the p-RSOS model, a restricted solid-on-solid model with a point-contact-type step–step attraction. In the step droplet zone in this model, the surface tension is discontinuous around the (111) surface and continuous around the (001) surface. The vicinal surface of 4H-SiC crystal in a Si–Cr–C solution is thought to be in the step droplet zone. The dependence of the vicinal surface growth rate and the macrostep size 〈n〉 on the driving force ∆μ for a typical state in the step droplet zone in non-equilibrium steady state was calculated using the Monte Carlo method. In contrast to the known step bunching phenomenon, the size of the macrostep was found to decrease with increasing driving force. The detachment of elementary steps from a macrostep was investigated, and it was found that 〈n〉 satisfies a scaling function. Moreover, kinetic roughening was observed for |∆μ| > ∆μR, where ∆μR is the crossover driving force above which the macrostep disappears.


Introduction
The faceted macrosteps on a crystal surface are known to degrade the grown crystal [1].Although studies have investigated methods of dispersing faceted macrosteps, an effective method has not yet been established.For example, in the case of solution growth for 4H-SiC, the faceted macrosteps remain near equilibrium.This formation of faceted macrosteps near equilibrium has been considered to be due to the effects of anomalous surface tension.To control the dynamics of macrosteps, the fundamentals of the phenomena with macrosteps must be clarified.
The connection between the surface tension and the instability with respect to macrostep formation is also scientifically interesting.In 1963, Cabrera and Coleman [2] phenomenologically studied anomalous surface tension and studied its effect on a vicinal surface near equilibrium.They assumed several anomalous surface tensions, and then discussed the possible equilibrium crystal shapes (ECSs).They also declared the instability with respect to the macrostep formation to be the result of anomalous surface tension [3].However, at that time, the microscopic model used to determine the anomalous surface tension was not provided.Jayaprakash et al. [4] studied the faceting transition of the ECS using a terrace-step-kink (TSK) model with the long-range step-step interaction expressed by the equation ∑ i =j ∑ ỹ g 0 /[ xi (ỹ) − xj (ỹ)] 2 [5][6][7], where g 0 is the coupling constant of the elastic repulsion, xi is the location of the ith step on a vicinal surface normal to the mean running direction of steps, and ỹ is the location along the steps.They showed that the step-step interaction affects the coefficient of the O(ρ 3 ) term in the surface free energy, which is given by f (ρ) = f (0) + γρ + Bρ 3 + O(ρ 4 ), (1) where ρ is the step density, γ is the step tension, and B is step interaction coefficient.This ρ-expanded form of the free energy excluding the quadratic term with respect to ρ is called the Gruber-Mullins-Pokrovsky-Talapov (GMPT) [8][9][10][11] universal form.It is well known that the ECS can be obtained using the Landau-Andreev method [12,13] or from the polar graph of the surface tension using the Wulff theorem [14][15][16][17].When the long-range step-step interaction is attractive, B becomes negative at low temperatures, and the slope of the surface on the ECS has a jump at the facet edge.This jump in the surface slope is referred to as a first-order shape transition [4,18] after the fashion of a phase transition.
The step bunching for Si near equilibrium in ultrahigh vacuum has been studied for Si(111) [19][20][21][22][23][24] and Si(113) [25][26][27][28].Williams et al. [19,20] have experimentally demonstrated that the step bunching of Si(111) is caused by the competition between the polymorphic surface free energies for the Si(111) (1 × 1) and (7 × 7) structures based on the GMPT surface free energy (Equation ( 1)).Additionally, Song et al. obtained the temperature-slope phase diagram for the step bunching on Si(113) by analyzing the results of X-ray diffraction, and have reported that an anomalous surface free energy caused by a step-step attraction may explain the step bunching phenomenon.Lössig [29] and Bhattacharjee [30] have stated that the TSK model with a short-range step-step attraction and with the long-range step-step repulsion can represent the step bunching phenomenon.By analyzing the scanning tunneling microscopy results on a vicinal surface of Si(113) tilted toward the [1 10] direction, van Dijken et al. [27] demonstrated that there is a step-step attraction to condense steps and a large long-range step-step repulsion.Shenoy et al. [31,32] have shown that the TSK model with a short-range step-step attraction and with the long-range step-step repulsion causes the periodic array of the n-merged steps using a renormalization group method.Einstein et al. [33] introduced the idea of the random matrix to the terrace width distribution to assist in the determination of the strength of the long-range step-step repulsion.
Step bunching or step faceting near equilibrium occurs without long-range step-step repulsion [18,[34][35][36][37][38][39][40][41][42][43][44][45].Rottman and Wortis [18] studied an ECS using a three-dimensional ferromagnetic Ising model with both nearest-neighbor (NN) and next-nearest-neighbor (NNN) interactions.They calculated ECS using mean-field theory.When the NNN interaction is negative, they showed that the ECS has the first-order shape transition at the (001) facet edge at low temperatures.Using a lattice model with a point-contact-type step-step attraction [38][39][40][41][42][43][44][45], Akutsu showed that a faceted macrostep self-assembles at equilibrium in association with the morphological change resulting from the discontinuous surface tension.The lattice model was a restricted solid-on-solid (RSOS) model with point-contact-type step-step attraction (p-RSOS model, Figure 1).The term "restricted" means that the height difference between NN sites is restricted to {0, ±1}.It was considered that the origin of the point-contact-type step-step attraction is the orbital overlap of the dangling bonds at the meeting point of the neighboring steps.The energy gained by forming the bonding state is regarded as the attractive energy between steps.The surface tension of the model was numerically calculated using the density matrix renormalization group (DMRG) method [46][47][48][49][50][51], and it was demonstrated that the surface tension of the vicinal surface tilted toward the 111 direction is discontinuous at low temperatures because of the point-contact-type step-step attraction.
This connectivity of the surface tension for the p-RSOS model is directly linked to the faceting diagram (Figure 2).There are two transition temperatures in the p-RSOS model-T f ,1 and T f ,2 .At temperatures of T < T f ,1 , the surface tension around the (111) surface is discontinuous; for T < T f ,2 (< T f ,1 ), the surface tension around the (001) surface is discontinuous.Based on the connectivity of the surface tension, the temperature region T < T f ,2 is called the step-faceting zone, the region T f ,2 ≤ T < T f ,1 is called the step droplet zone, and the region T f ,1 ≤ T is called the GMPT [8,9] zone.
Moreover, the plot of the roughening transition temperature of the (001) surface divides the step droplet and GMPT zones into step droplet zones I and II and GMPT zones I and II, respectively [43,44].Step droplet II  Triangles: calculated values of T f ,2 .Open circles: calculated roughening transition temperatures of the (001) surface.Solid line: zone boundary line calculated using the two-dimensional Ising model.For the definitions of and details about the QI Bose solid, liquid, and gas, please refer to Akutsu [43].This figure was reproduced from Akutsu [43].GMPT: Gruber-Mullins-Pokrovsky-Talapov.
In our previous study [45], the height profile of a faceted macrostep at equilibrium was investigated on the p-RSOS model, and it was demonstrated that the characteristics of the height profile of a macrostep can be classified by the connectivity of the surface tension.The characteristics of the height profile of a macrostep are irrelevant to the details of the crystal structure.Hence, the height profile of a macrostep can be used to determine in which zone the surface exists.For example, the height profile of the macrostep in the case of 4H-SiC [1] is similar to the profile in the step droplet zone [45].This suggests that the surface tension of 4H-SiC around the faceted side surface is discontinuous.
This paper is organized as follows.In §2, the model and the discontinuous surface tension are briefly explained.In §3, the results obtained using the Monte Carlo method are presented.The process of detaching the steps from a macrostep is discussed in §4 through the analysis of the results in the case of small |∆µ|.In §5, the case of large |∆µ| is considered; for this case, the size of a macrostep is modeled using a scaling function, and a crossover point ∆µ R associated with kinetic roughening is introduced.The implications of the results are discussed in §6, and §7 concludes the paper.

Restricted Solid-on-Solid Model with Point-Contact-Type Step-Step Attraction
The microscopic model considered in this study is the p-RSOS model (Figure 1).The Hamiltonian 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 step energy, δ(a, b) is the Kronecker delta, and int is the microscopic step-step interaction energy.The summation with respect to (n, m) is taken over all sites on the square lattice.The RSOS condition is required implicitly.When int is negative, the step-step interaction becomes attractive (sticky steps).
It should be noted that the p-RSOS model (Equation ( 2)) automatically includes the "entropic step-step repulsion".Since the p-RSOS model is an RSOS model, the overhang structures with respect to the height of the surface are inhibited by the geometrical restriction.This is the microscopic origin of the entropic step-step repulsion.

Discontinuous Surface Tension
The polar graphs of the surface tension and the surface free energy were calculated using the DMRG method [46][47][48][49][50][51], and are shown in Figure 4.In low-dimensional cases, more precision is required than can be provided by a mean field calculation of the partition function [56].Hence, to obtain reliable results, the DMRG method-which was developed for one-dimensional (1D) quantum spin systems [46]-was adopted.The transfer matrix version of the DMRG method-known as the product wave function renormalization group (PWFRG) method [49][50][51]-was used in this study.Details of the calculation method for the surface tension and the surface free energy are given in Appendix A. In Figure 4, the angle θ is the tilt angle from the (001) surface toward the 111 direction.The surface gradient p is related to θ as |p| = ± tan θ.The surface tension was calculated from the surface free energy f (p) as The calculated surface tension and the surface free energy at k B T/ = 0.63 and int / = −0.9 are shown in Figure 4.The surface tension contains discontinuities near the (001) and (111) surfaces in the case of int / = −0.9,k B T f ,2 / = 0.613 ± 0.02, and k B T f ,1 / = 0.709 ± 0.02.At k B T/ = 0.63, the surface is in the step droplet zone (Figure 2), where the surface tension is continuous around the (001) surface but discontinuous around the (111) surface.
The ECS calculated using the DMRG method (Appendix A) is plotted in red in Figure 4a.This ECS result agrees well with the ECS obtained using the Wulff theorem [14][15][16][17].In step droplet zone I, the (001) surface meets the curved area without a discontinuity in its slope point P in Figure A1) (The Gauss curvature, which is the determinant of the curvature tensor, jumps at the facet edge of the (001) surface), whereas the (111) surfaces meet the curved areas with a discontinuity in their slopes (point Q in Figure A1).
At the zone boundary lines in Figure 2, the following conditions are satisfied [40,43]: where q(T) is the surface gradient on the vicinal surface near the (111) surface, f (111) (q) is the surface free energy of the vicinal surface near the (111) surface, γ q,n (T) is the step tension of an n-merged (negative) step [43], and B q,eff (T) and C q,eff (T) are coefficients.In GMPT zones I and II, since B q,eff (T) > 0, the surface free energy f (111) (q) has a form similar to Equation (1).As the temperature decreases, B q,eff (T) decreases and C q,eff (T) increases.For T < T f ,1 , where the vicinal surface exists in the step droplet zone, B q,eff (T) becomes negative and the first-order transition occurs.Hence, the upper zone boundary line of T = T f ,1 is a critical curve.
The key points to obtain Equations ( 4) and ( 5) are the meeting of neighboring steps and the inhomogeneity of the vicinal surface [40,43].For the vicinal surface tilting toward the 111 direction, two neighboring steps can occupy one site at the same time, and no more than three steps can occupy one site at a time because of the geometrical restrictions of the RSOS model.Hence, the surface cannot be mapped to a 1D fermion model [57][58][59].The double occupancy of a site gives rise to the point-contact-type step-step interaction.When the interaction is repulsive, the term C q,eff (T)|q| 4  is present [43,60,61], whereas in the case of attractive interaction, the vicinal surface becomes inhomogeneous and can be expressed as a mixture of the various n-merged steps (macrosteps) [40].Since the population of the n-merged steps depends on the surface slope, B q,eff (T)|q| 3 is affected by the point-contact-type step-step attraction through the slope dependence of the size of the macrosteps.

Monte Carlo Method
To study the non-equilibrium steady state with macrosteps, the vicinal surface of the following Hamiltonian with a fixed number N step of steps was investigated using the Monte Carlo method with the Metropolis algorithm: where ∆µ = µ ambient − µ crystal is the driving force, µ crystal is the chemical potential of the bulk crystal, and µ ambient is the chemical potential of the ambient phase.The explicit form of ∆µ is given in Markov [62].When ∆µ > 0, the crystal grows because its chemical potential is lower than that of the ambient phase, whereas when ∆µ < 0, the crystal recedes (evaporate, dissociates, or melts).
The explicit procedure of the application of the Monte Carlo method in this study is as follows.At the initial time, the steps are positioned at equal distances.Then, the lattice site to be updated is randomly selected.The surface structure is updated non-conservatively using the Monte Carlo method with the Metropolis algorithm.With the RSOS restriction taken into consideration, the structure is updated with probability 1 when ∆E ≤ 0 and with probability exp(−∆E/k B T) when ∆E > 0, where ∆E = E f − E i , E i is the energy of the present configuration and E f is the energy of the updated configuration.The energy is calculated from the Hamiltonian (Equation ( 6)).A periodic boundary condition was imposed in the direction parallel to the steps.In the direction normal to the steps, the lowest side of the structure was connected to the uppermost side by adding a height with a number N step of steps.
Snapshots of the simulated surfaces after 2 × 10 8 Monte Carlo steps per site (MCS/site) are shown in Figure 5.The height profiles of the surfaces for the cross section along the bottom of the top-down views are also shown as side views.

Macrostep Size and Surface Growth Rate
At equilibrium (∆µ = 0), the vicinal surface showed a homogeneous stepped surface for a small mean surface slope p (| p| < p 1,eq , Figure 4b) because the surface tension around the (001) surface is continuous in the step droplet zone.In contrast, when the mean surface slope satisfies p 1,eq < p < √ 2, homogeneous stepped surfaces are thermodynamically unstable [40].Then, the surface is realized through two-surface coexistence; the two surfaces are the surface with a slope equal to p 1,eq and the (111) surface [40][41][42][43][44][45].This is illustrated in Figure 5a.Because the (111) surface is smooth (which means a small number of kinks exist on it), it hardly moves.From the time-dependent Ginzburg-Landau equation of the surface [63,64], the smooth surface does not move because the surface stiffness is divergent [65].The reason the faceted macrosteps move at equilibrium is the finiteness of the system size [45].
For ∆µ > 0, the size of a single macrostep decreases as ∆µ increases, as shown in Figure 5. Furthermore, for ∆µ/ 0.006, the macrostep disassembles to form a homogeneous rough surface.Because the temperature was the same for all cases shown in Figure 5, the surfaces with ∆µ/ 0.006 roughen kinetically.Interestingly, the change in the size of macrostep is symmetric in the case of ∆µ < 0, where the steps are receding.The patterns obtained in the case of ∆µ < 0 (e.g., ∆µ/ = −0.003) is quite similar to the patterns obtained in the case of |∆µ| (e.g., ∆µ/ = 0.003).
To study the characteristics of the vicinal surface on a mesoscopic scale (approximately 10 nm to 1 µm) in detail, the size of the macrosteps and the growth rate of the surface were measured during the Monte Carlo simulation.To evaluate the size of a macrostep, the number n of elementary steps in a locally merged step [40] was introduced.The average size of the locally merged steps is obtained as where N n is the number of n-merged steps on the vicinal surface.The time evolutions of n at different values of ∆µ/ are shown in Figure 6a.As shown in Figure 6a, for |∆µ| ≥ 0.003, n is constant near 2 × 10 8 MCS/site.Hence, surfaces with |∆µ| ≥ 0.003 are in non-equilibrium steady state.The time evolutions of n in non-equilibrium steady state were also obtained for |∆µ| = 0.001 and 0.002, the results of which are not shown in Figure 5.The ∆µ-dependence of n at 2 × 10 8 MCS/site is shown in Figure 3a.For small |∆µ|, n decreases linearly as |∆µ| increases.
To estimate the growth rate v, the average surface height h(t) was calculated as h(t) = (1/N ) ∑ n,m h(n, m).The time evolutions of h(t) at different values of ∆µ/ are shown in Figure 6b.As shown in the figure, h(t) increases or decreases linearly with increasing t.Hence, v is defined as The ∆µ-dependence of v is shown in Figure 3b.|v| increases as |∆µ| increases.

Size of a Macrostep
As shown in Figure 5b,c, for sufficiently small |∆µ|, n can be approximated as ) where n m is number of elementary steps contained in the most dominant size of the macrosteps, N m is the number of n m -merged steps, N 1 is the number of elementary steps outside of the macrostep, z is the ratio N 1 /N step , and p 1 is the slope of the "terrace" surface that is in contact with the (111) surface.From the snapshots, it was assumed that N m = 2.
From the definitions of z and Equation ( 9), z can be calculated from n as Thus, the curve of best fit was obtained as z = 0.307 + 17.1|∆µ|/ + 6.86 × 10 3 (|∆µ|/ ) 2 by applying the method of least squares to the values of z, which was estimated from n using Equation (11) on a system of size 240 The values of n reproduced by the best fit equations are plotted as blue lines in Figure 3a.The lines agree well with the values of n for small |∆µ|.The best fit equation reveals the following.In the case of p = 0.707, 69.0% of all N step steps self-assemble to form a macrostep in the limit of |∆µ| → 0. For small |∆µ|, n decreases linearly as |∆µ| increases.This indicates that n has a cusp singularity at |∆µ| = 0, because the slope of the line in Figure 3 in the limit of ∆µ → +0 (from the right-hand side) is different from the slope of the line in the limit of ∆µ → −0 (from the left-hand side).

Growth Rate
Next, the growth rate was investigated.Here, the following model for the time evolution of n was considered: where n + is the rate at which elementary steps join the macrostep, and n − is the rate at which elementary steps detach from the macrostep.In the growth condition (∆µ > 0), when n + < n − , the macrostep dissociates.In this case, n + limits the surface growth rate.In contrast, when n + > n − , n increases up to N step , and n − limits the surface growth rate.When the surface is in steady state, n + = n − .Because a surface with a slope of p 1 shows step flow growth, n + is considered to be p 1 v 1 for small |∆µ|, where v 1 is the growth rate of a single step.Hence, the growth rate in steady state is expressed as p 1 is obtained from Equation (10) as At equilibrium, the terrace surface with a slope of p 1,eq = 0.3532, which was calculated using the DMRG method, coexists with the (111) surface, which forms the side surface of the macrostep.For |∆µ| > 0, p 1 can be obtained from Equation ( 14) using the value of z.In the same manner as for z, the curve of best fit was obtained for p 1 on a system of size 240 √ 2 × 240 √ 2 as p 1 = 0.332 + 15.6|∆µ|/ + 4.43 × 10 3 (|∆µ|/ ) 2 .In this manner, it was found that p 1 increases linearly with |∆µ| for small |∆µ|.Equation ( 13) with v 1 = 0.132∆µ is plotted in Figure 3b as blue lines.
Since n + = n − in steady state, the equation for n − may be used.To model the mechanism of the detachment of an elementary step from a macrostep, the two-dimensional (2D) nucleation mode was considered.In this case, v should be proportional to exp(−const./|∆µ|) [66].However, the Monte Carlo data could not be fit by this equation.The 2D multi-nucleation mode was also considered.In this case, v should be proportional to |∆µ| 2/3 exp(−const./|∆µ|)[62,66].However, the Monte Carlo data could not be fit by this equation either.From this, the detachment of an elementary step from a macrostep was considered to be caused by the "noise" of the attachment and detachment of atoms ("Atoms" in the model correspond to unit cells) from the ambient phase in association with thermal fluctuations.
Analysis of the results obtained by the Monte Carlo calculations for n yielded the scaling function Y(x) for n , and the power law behavior for v (Figure 7).The scaling function where x = ln |∆µ|/ .x R is determined as the point of intersection of the blue and the pink lines in Figure 7b.The value of x R corresponds to ∆µ R / = 0.005.This yields This means that n − n ∞ shows power law behavior, with the power depending on the system size.The lines based on the power law equation for a system size of 240 √ 2 × 240 √ 2 are shown as pink lines in Figure 3a.The values obtained using the Monte Carlo method agree well with the power law functions for large |∆µ|.
v also shows power law behavior.ln |v| (Figure 7b) increases linearly with increasing ln |∆µ|/ for |∆µ| > ∆µ R .Thus, v can be expressed as Here, the choice of β as the symbol for the exponent is in accordance with Reference [70].The power law equation given by Equation ( 17) is plotted as pink lines in Figure 3b.The values obtained using the Monte Carlo method agree well with the lines for large |∆µ|.As shown in Figure 7, ∆µ R is not a roughening transition point in the Kosterlitz-Thouless universality class [71][72][73][74], but a crossover point from the two-surface coexistence state to a homogeneous rough surface.To clearly observe this crossover, the case with int / = −1.3 and k B T/ = 0.83-which is also in step droplet zone I-was calculated (Figure 8).In this case, n and v showed behaviors quite similar to the case of int / = −0.9 and k B T/ = 0.63.In the case of int / = −1.3,∆µ R = 0.0045 was obtained.The shift from the power law behavior for v at small |∆µ| in the case of int / = −1.3 was more evident than in the case of int / = −0.9.
It should be noted that the surface structure contains locally merged steps to some extent, though the surface is homogeneous at the thermodynamic limit (L → ∞).For |∆µ| > ∆µ R , the power of v is similar to that in the case of the original RSOS model ( int = 0, Figure 7).In spite of this, the growth rate remained lower than that in the original RSOS model at large |∆µ| (Figure 3b).To see the structure at a more microscopic scale, the surface was simulated with a small system size, as shown in Figure 5f.Though the surface appeared homogeneous for a system size of 240 √ 2 × 240 √ 2, several locally merged steps were observed at a system size of 40 √ 2 × 40 √ 2. These locally merged steps cause the kink density of the surface to decrease.Hence, the growth rate in this case was lower than that for the original RSOS model.

Discussion
At equilibrium, planar surfaces such as the (001) surface are smooth at temperatures less than the roughening transition temperature T R .The height-height correlation function G(r) = (h(r) − h ) 2 of the planar surface is constant for any r on a smooth surface.In contrast, for T ≥ T R , planar surfaces are rough, and G(r) is logarithmically divergent with respect to |r| [11,[71][72][73][74].For the vicinal surface, G(r) is logarithmically divergent with respect to |r|, though the terrace is smooth [71,72].Hence, in the present case, the height-height correlation function of the surface with a slope of p 1,eq is logarithmically divergent, whereas that of the (111) surface-which is the side surface of the macrostep-is constant (non-divergent) at equilibrium.For 0 < |∆µ|, the height-height correlation function for a surface with a slope of p 1 is considered to be logarithmically divergent with respect to |r|.This is also supported by the exponent β introduced in §5.The growth rate of the rough surface is known to increase linearly as ∆µ increases [66,75].The exponent β being close to 1 is consistent with the results obtained through calculations using the discrete Gaussian model [75].
The kinetically roughened state for |∆µ|/ > ∆µ R is somewhat different from the state in step droplet zone II [43,44].In step droplet zone II, the (001) surface is rough because the temperature is higher than its roughening transition temperature.Hence, the correlation length ξ for G(r) on the (001) surface is divergent.Furthermore, n depends on the surface slope p as n = n 0 + Cp 2 + O(p 3 ) in the limit p → 0, where n 0 and C are constants.In contrast, the correlation length ξ for the (001) surface in step droplet zone I is considered to be finite.n varies with the surface slope p as n = n 0 + Cp + O(p 2 ) in the limit p → 0. The locally merged steps remain to some extent (Figure 5f).
The relaxation time required to reach steady state increases as |∆µ| decreases (Figure 6a).The time to relax to steady state is associated with the power law behavior of the surface [55,69,[76][77][78].
The focus of the present study was the dynamics affected by the surface tension.As indicated in the cases of 4H-SiC [1] and Si(113) [21,[25][26][27][28][29][30][31][32][33], elastic step-step repulsion expressed by the formula ∑ i =j ∑ ỹ g 0 /[ xi ( ỹ) − xj ( ỹ)] 2 [5][6][7] is important in real systems near equilibrium.The elastic effects are considered to weaken the step-step attraction in the present study.In the step droplet zone, the elastic step-step repulsion produces two effects: the shift of the zone boundary lines [43,45] and the formation of a regular array of n-merged macrosteps [31,32].Since the elastic step-step repulsion causes B q,eff (T) (Equation ( 4)) and γ q,n (T) (Equation ( 5)) to increase, the zone boundary lines shift to lower temperatures.In the short range, the step-step attraction dominates and sticks steps together, whereas in the long range, the step-step repulsion dominates and separates n * -merged macrosteps, where n * is the size of the macrostep with the lowest surface free energy.Thus, the following guidelines were obtained to disperse macrosteps: (1) the temperature should be raised; (2) the elastic step-step repulsion should be increased; and (3) the absolute value of the driving force should be increased.
In real systems, other effects, such as surface diffusion [52][53][54][55] and polymorphic effects [19,20], should also be taken into consideration.The combination of these effects and the effect of the discontinuous surface tension will be considered in future work.

Conclusions
The effect of the driving force ∆µ on the size of a faceted macrostep and the growth rate of the vicinal surface in non-equilibrium steady state were investigated using the Monte Carlo model.
Step droplet zone I for the p-RSOS model was the focus of this study.

•
As |∆µ| increases, the size of the macrostep n decreases, whereas the growth rate |v| increases.

•
At small |∆µ|, the ∆µ-dependence of n and v can be explained by the attachment and detachment of elementary steps to and from the macrostep.

•
When |∆µ| ≥ ∆µ R , the macrostep disassembles, and the surface roughens kinetically.∆µ R / = 0.005 is the crossover point from the two-surface coexistent state to the rough surface state.

Figure 1 .
Figure 1.(a) Perspective view of the restricted solid-on-solid (RSOS) model tilted toward the 110 direction; (b) Top-down view of the RSOS model.Thick blue lines represent surface steps.This figure was reproduced from Akutsu [45].

Figure 2 .
Figure 2.Faceting diagram of the p-RSOS model (RSOS model with point-contact-type step-step attraction) for a vicinal surface obtained using the density matrix renormalization group (DMRG) method.Squares: calculated values of T f ,1 .Triangles: calculated values of T f ,2 .Open circles: calculated roughening transition temperatures of the (001) surface.Solid line: zone boundary line calculated using the two-dimensional Ising model.For the definitions of and details about the QI Bose solid, liquid, and gas, please refer to Akutsu[43].This figure was reproduced from Akutsu[43].GMPT: Gruber-Mullins-Pokrovsky-Talapov.

Figure 4 .
Figure 4. (a) Polar graph of the surface tension (Wulff figure) and Andreev free energy (equilibrium crystal shape, ECS) calculated using the DMRG method.This figure was reproduced from Akutsu [45]; (b) Surface free energy.k B T/ = 0.63, int / = −0.9(step droplet zone I).Red lines: Z(R) calculated using the DMRG method.Blue lines and squares: polar plots of (a) the surface tension γ surf (θ)/ (−54.74 • < θ < 54.74 • ), where θ is the tilt angle of the vicinal surface from the (001) surface toward the 111 direction; and (b) the surface free energy.Pale blue lines: values for the metastable surfaces for (a) the surface tension; and (b) the surface free energy.End points of the pale blue lines (p * and −p * ): approximate spinodal points.Point O: Wulff point.surf was assumed to equal .An enlarged figure of the ECS near the facet edges is shown in Appendix A.