Height of a faceted macrostep for sticky steps in a step-faceting zone

The driving force dependence of the surface velocity and the average height of faceted merged steps, the terrace-surface-slope, and the elementary step velocity in the non-equilibrium steady-state are studied using the Monte Carlo method. The Monte Carlo study is based on a lattice model, the restricted solid-on-solid model with point-contact type step--step attraction (p-RSOS model). The temperature is selected to be in the step-faceting zone where the surface is surrounded by the (001) terrace and the (111) faceted step at equilibrium. Long time simulations are performed at this temperature to obtain steady-states for the different driving forces that influence the growth/recession of the surface. A Wulff figure of the p-RSOS model is produced through the anomalous surface tension calculated using the density-matrix renormalization group method. Although the p-RSOS model is a simplified model, the model shows a wide variety of dynamics in the step-faceting zone. There are four characteristic driving forces, $\Delta \mu_y$, $\Delta \mu_f$, $\Delta \mu_{co}$, and $\Delta \mu_R$. For $\Delta \mu_{co}<|\Delta \mu|<\Delta \mu_R$, the surface grows/recedes with the successive attachment-detachment of steps to/from a macrostep. When $|\Delta \mu|$ exceeds $\Delta \mu_R $, the macrostep vanishes and the surface roughens kinetically. Classical 2D heterogeneous multi-nucleation was determined to be valid with slight modifications based on the Monte Carlo results of the step velocity and the change in the surface slope of the"terrace". The finite size effects were also determined to be distinctive near equilibrium.


I. INTRODUCTION
Faceted macrosteps have sometimes been considered to degrade the quality of grown crystals. In the case of solution growth for 4H-SiC [1] as an example, which is expected to be used for future power devices, the faceted macrosteps near equilibrium hinder the preparation of good quality crystals that satisfy the requirements for electrical devices. Therefore, to control the dynamics of macrosteps, the fundamentals regarding the formation of macrosteps should be clarified.
For smooth surfaces, the approaches based on the timedependent Ginzburg-Landau (TDGL) equation for surface motion [2][3][4][5][6], which is used to study rough surfaces, are not valid. The nucleation model is instead known to be more effective. In Saito's solution for the TDGL equation of a smooth surface with a modified discrete Gaussian (MDG) model [5], the two-dimensional (2D) nucleus is not included. Based on the TDGL equation of the surface, the surface cannot grow until ∆µ exceeds a "spinodal" value, ∆µ c . Here, ∆µ = µ ambient − µ crystal is the driving force, where µ crystal is the chemical potential of the bulk crystal and µ ambient is the chemical potential of the ambient phase. This suggests that the excitation * nori3@phys.osakac.ac.jp, nori@phys.osakac.ac.jp of islands on the surface, which can be 2D nuclei, corresponds to a higher order response to the driving force with respect to the crystal growth.
The phase field method [7] is also known as a powerful tool to study the solidifications or other non-equilibrium phenomena accompanied by phase changes. In the phase field method, the phase boundary is assumed to be analytic and differentiable. The faceted structure on the surface can be simulated with the phase field method by the introduction of a strong anisotropy to the interface tension. The planes can be imitated by curved surfaces with small curvature. However, this slight difference in the curvature causes a significant difference in the long time behavior, i.e., the non-equilibrium steady-state behavior. For example, the singularity of the flat smooth interface inhibits the growth/recession of the crystal without 2D nucleation processes or screw dislocations [8].
With respect to nucleation, the nucleation model interprets how the large clusters (domains) are formed; therefore, the nucleation model is widely accepted for study of the dynamics around the first-order phase transitions [6,[8][9][10]. However, for quantitative study, the classical nucleation theory [8,11,12] often disagrees with the experimental observations or the results of large-scale molecular dynamics by an order of 10 3 [13]. Therefore, many improvements in the nucleation theory have been reported [9,10]. In recent years, the notion of modern "multi-nucleation" [13][14][15][16], where the growing clus-  [26], with the permission of Hindawi Publishing Corporation.
ters change the crystal structure at a certain size during growth, has attracted attention and opened a new research area.
A vicinal surface with faceted macrosteps is surrounded by smooth surfaces; therefore, the surface is considered to grow/recede by way of 2D nucleation [8,11,12]. However, macrosteps are considered to be unstable at equilibrium [17,18] without impurities or adatoms. Therefore, the dynamics of a vicinal surface with faceted macrosteps in the non-equilibrium steadystate have not been studied sufficiently.
With respect to the stable faceted macrostep, we reported a study on a restricted solid-on-solid (RSOS) model with a point-contact-type step-step attraction (p-RSOS model, Fig. 1) [19][20][21][22][23][24][25][26]. Here, "restricted" means that the height difference between nearest neighbor sites is restricted to {0, ±1}. The origin of the point-contacttype step-step attraction is considered to be the orbital overlap of the dangling bonds at the meeting point of neighboring steps. The energy gained by the formation of a bonding state is regarded as the attractive energy between steps.
The characteristic of the p-RSOS model is the discontinuous surface tension at low temperatures [21][22][23]. Macrosteps are stabilized when the surface free energy has anomalous anisotropy [27,28]. Therefore, it seems a simple task to calculate the surface free energy explicitly with a standard method of the statistical mechanics. However, due to the large contribution of thermal fluctuations in a low-dimensional substance [29], it is difficult to obtain reliable results theoretically using the mean field approximation. Therefore, to obtain reliable results, the density-matrix renormalization group (DMRG) method [30][31][32][33][34][35] is used for calculation of the surface tension (the surface free energy per normal unit area).
A faceting diagram (Fig. 2) that corresponds to the connectivity of the surface tension was obtained [24]. The faceted macrostep is stabilized in the step faceting zone and in the step-droplet zone. In the Gruber-Mullins-Pokrovsky-Talapov (GMPT) zone, there is no faceted macrostep; the vicinal surface obeys the GMPT universal behavior [36][37][38]. In the step faceting zone, the vicinal 0 -0.5 Step-faceting (QI Bose solid +Vacuum) Step droplet II ( [24]. Reproduced from [24], with the permission of AIP Publishing. surface consists of (001) terraces and a single (111) surface which forms the side surface of a faceted macrostep. In contrast, the vicinal surface in the step droplet zone consists of a single (111) surface and surfaces with slope p 1 . The characteristic of the height profile of the faceted macrostep is also classified by the connectivity of the surface tension at equilibrium [26].
In this article, we study the driving force dependence of the surface velocity V , the average height of faceted merged steps n , the terrace-surface-slope p 1 , and the elementary step velocity v step in the non-equilibrium steady-state using the Monte Carlo method for a vicinal surface with a faceted macrostep in the p-RSOS model. The temperature is selected to be in the step-faceting zone where the surface is surrounded by a (001) terrace and (111) faceted step at equilibrium. The Wulff figure of the p-RSOS model is produced from the anomalous surface tension calculated with the DMRG method.
This article is organized as follows. In §II, the model Hamiltonian and the surface tension calculated with the DMRG method are shown. In §III, we present the results obtained from long time Monte Carlo simulation for an average-sized faceted macrostep n and the growth rate of the surface, V . In §IV, the Monte Carlo results are analyzed based on the classical 2D heterogeneous nucleation and the classical 2D heterogeneous multi-nucleation. The driving force dependence of the slope of the "terrace" and the step velocity are also presented. The crossover to the kinetic roughened surface is discussed in §IV D. Further discussions are given in §V, and conclusions are presented in §VI.

A. The p-RSOS model
The microscopic model considered in this study is the p-RSOS model ( Fig. 1) [19][20][21][22][23][24][25]. In this model, "an atom" corresponds to a unit cube. 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).

B. Discontinuous surface tension
The surface tension is the surface free energy per unit normal area. The surface tension γ surf (p) was calculated from the surface free energy f (p) per projected x-y area for the vicinal surface as where p = (p x , p y ) is the surface gradient of the vicinal surface [60]. The surface free energy f (p) is calculated from the Andreev free energy [61], which is the grand potential on the grand partition function with respect to the number of steps (Appendix A). The DMRG method was applied for calculation of the grand partition function. The transfer matrix version of the DMRG method, which is known as the product wave function renormalization group (PWFRG) method [33][34][35], was used in Surface tension γ(p)/ǫ. Pale lines: Andreev's free energy calculated using the DMRG method [30][31][32][33][34][35]. ǫ surf is assumed to equal ǫ.
this study. Details of the method for calculation of the surface tension and the surface free energy are given in Appendix A.
The polar graphs of the surface tension are shown in Fig. 3. The surface gradient p is related to the tilt angle θ as p = ± tan θ. It should be noted that the Andreev's free energy is similar to the equilibrium crystal shape (ECS), which is the shape with the least total surface free energy. The ECS is obtained by the Landau-Andreev method [61,62]. Alternatively, the ECS is obtained from the surface tension with the Wulff construction [8,[63][64][65] based on the Wulff theorem.
In Fig. 3, there are values for the (001) and (111) surfaces [21]. The vicinal surfaces between the (001) surface and the (111) surface did not appear, because these vicinal surfaces are thermodynamically unstable. This is the characteristic of the profile of the faceted macrostep in the step-faceting zone at equilibrium [24].

A. 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 t is the time measured by the Monte Carlo steps per site (MCS/site). When ∆µ > 0, the crystal grows, whereas when ∆µ < 0, the crystal recedes (evaporates, dissociates, or melts).
The explicit procedure for application of the Monte Carlo method in this study is as follows. At the initial time, the vicinal surface is set with an initial configuration. The lattice site to be updated is then 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 using the Hamiltonian (Eq. (3)).
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 the addition of a height with a number N step of steps.
Two types of initial configuration for the steps were prepared: a train of elementary steps with equal distance (the TS configuration) and one macrostep with the (111) side surface (the MS configuration). In both configurations, the mean surface slopep = N step /L is kept constant, where L is the linear size of the system. Figure 4 shows snapshots of the vicinal surfaces at 4 × 10 8 Monte Carlo MCS/site, where the initial configuration is TS.

B. Time evolution of surface height
To study the characteristics of the vicinal surface at the mesoscopic scale (20 nm to 500 nm), the average height of merged steps [21] and the growth rate of the surface were calculated using the Monte Carlo method. To evaluate the size of a macrostep in detail, we consider the number of elementary steps in a locally merged step n. The average height of the merged steps is obtained as follows: wherex is selected as the 110 direction (normal to the mean step-running direction),ỹ as the 1 10 direction (along the mean step-running direction), N step is the total number of elementary steps, and n step is the number of merged steps. Time evolutions of n are shown in Fig. 5(a). After 2 × 10 8 MCS/site, n is almost constant for ∆µ/ǫ > 0.05. Therefore, the values obtained within 0-2 × 10 8 MCS/site were discarded, and n was averaged over successive 2 × 10 8 MCS/site. Figure 6 shows the ∆µ dependence of n .
To estimate the growth rate of the surface V , the average surface heighth(t), was calculated, wherē Time evolutions ofh(t) are also shown in Fig. 5  The growth rate of the surface V , is defined as:

C. Characteristic driving forces
There are several characteristic driving forces in the non-equilibrium steady-state, which are listed in Table  I. Each value will be explained in detail in the following sections. Here, we will explain them briefly.
First, ∆µ R (L) is defined as the minimum |∆µ| so that the macrostep disappears ( Fig. 4)(f), §IV B 2). L is the linear size of the system. In this region, the surface roughens kinetically. The velocity of the surface V , and the average height of the merged step n , exhibits power law behavior with respect to |∆µ| − ∆µ co (L) if we introduce ∆µ co (L) ( §IV D). ∆µ co (L) is determined so that the log-log plot with respect to V vs. |∆µ| − ∆µ co (L) becomes linear. The Monte Carlo results show that ∆µ R (L) − ∆µ co (L) is not significantly dependent on L (Table I).
Next, we define ∆µ f (τ, L), where τ is the observation time. ∆µ f (τ, L) is defined as the maximum |∆µ| so that growth/recession is inhibited during the observation time. In the region |∆µ| < ∆µ f (τ, L), the   growth/recession of the vicinal surface freezes due to the finite size and time effect ( §IV A). The surface does not reach the true non-equilibrium steady-state; therefore, the surface morphology is dependent on the initial configuration of the surface. n in Figs. 6(a) and (b) for |∆µ| < ∆µ f (τ, L) is strongly dependent on the initial configuration.
For ∆µ f (τ, L) < |∆µ| < ∆µ co (L), the surface grows/recedes intermittently in the manner of 2D heterogeneous nucleation ( §IV A, Figs. 4(b) and (c)). n of the TS initial condition is smaller than n of the MS initial configuration. On the other hand, the growth rates with the TS initial configuration agree well with those for the MS initial configuration. An island on the (001) surface or an island on the (111) surface is formed at the edge of the faceted macrostep, as shown in Figs. 4(b) and (c).

Growth/recession rate of the surface
For ∆µ f (τ, L) < |∆µ| < ∆µ co (L) (Table I), the surface grows/recedes by step-detachment through 2D "heterogeneous" nucleation (Fig. 8). This means that the growth/recession of the surface occurs intermittently. The nuclei were created at the lower/upper side of the macrostep-edge (Figs. 4(b) and (c)) in the growth/recession mode, respectively. We describe the side length of the critical nucleus on the (001) surface as l c , and that on the (111) surface as l ′ c (Fig. 8). Figure  8 shows l c ∼ l ′ c of the 2D critical nucleus for the stepdetachment with the growth of the surface (∆µ > 0). For ∆µ < 0, an elementary step detaches from the upper side of the (111) surface by forming a critical negative nucleus. The shape of the critical negative nucleus is similar but reversed to the shape shown in Fig. 8(c) (also see Fig. 4(c)).
The Gibbs free energy G(l) (or G(l ′ )) of the island attached to the faceted macrostep is expressed as: where S(l) is the area of the island, Γ(l) is the total step free energy at the edge of the island of the elementary step, and γ (110) 2 is the step free energy of the doubly merged step. For the critical nucleus, the Gibbs free energy is a minimum with respect to the shape, but is a maximum with respect to the size [11]. The shape of the critical nucleus is similar to the equilibrium island shape (Fig. 8(c)). The size of the critical nucleus is determined as explained in Appendix B, i.e., l c = l c,0 ǫ/|∆µ|, G(l c ) = G(l c,0 )ǫ/∆µ, From the classical nucleation theory [6,11,12], the nucleation frequency I n is expressed as follows:  where Z is the Zeldovich factor, N 0 = 1/( √ 2a) is a lattice-point density, and C is a coefficient relating to geometry. The waiting time for a single nucleation t n is where g * /∆µ = G(l c )/(k B T ). In the limit of L → ∞, t n reduces to zero because 2D nuclei are formed somewhere at the step edge of the macrostep. The growth rate V is expressed using t n as ln(|V |) is shown with the horizontal axis being 1/(|∆µ|/ǫ) in Fig. 9(b). The Monte Carlo results gave by fitting a line in Fig. 9(b) to the values of ∆µ y (L) < |∆µ| < ∆µ co (L) with L = 240 √ 2. Therefore, we have g * MC = 0.509ǫ. This line is shown by the dotted line in Fig. 9(b). Z/( √ 2C) = 4.5 × 10 −3 was also obtained from the prefactor on the right-hand side of Eq. (12). Using this value with Eq. (11), we show the explicit values for t n with L = 240 √ 2 for several driving forces (Table II). To numerically obtain l c , l c,0 and G(l c,0 ) were calculated using the 2D Ising model [66][67][68] g * obtained from Eq. (14) is Table II shows the size of the l c for several ∆µ. g * Ising is slightly smaller than g * MC obtained by the Monte Carlo method. The reason for this seems to be a reduction of the entropy due to finite size of the critical island. The values of l c in Table II indicate the size of the critical island to be less than approximately 100 in the present simulations. With such a short length, the entropy term in the step tension becomes smaller than the entropy term in the step tension with infinite length. Therefore, g * MC becomes larger than g * Ising calculated for the infinite length of the domain boundary line.
Due to the long waiting time for the 2D heterogeneous nucleation at the edge of the faceted macrostep, the surface growth occurs intermittently. The intermittent growth/recession can be observed explicitly in the case of |∆µ/ǫ| = 0.03, as shown in Fig. 5(c).

Average height of merged steps
To understand the ∆µ dependence of the vicinal surface morphology we consider a step-attachmentdetachment model for the time evolution of n [69]: where n + is the rate when the elementary steps catch up to a macrostep, and n − is the rate when the elementary steps detach from a macrostep. When n + < n − , a macrostep dissociates, whereas when n + > n − , n increases up to N step , where N step is the total number of elementary steps on the surface. In this case, n − limits the growth/recession rate of the surface. At steady-state, n + = n − = V /a, where a is the height of the elementary step.
In the region of ∆µ f (τ, L) < |∆µ| < ∆µ co (L), n + is considered to be n + ≈ ρ 1 v 1 , where ρ 1 is the density of elementary steps on the "terrace", and v 1 is the step velocity of an elementary step perpendicular to the mean running direction of the step. On the other hand, n − is proportional to the 2D heterogeneous nucleation rate. Therefore, n + > n − is expected because the growth rate of an elementary step v 1 is relatively large (for example, §IV B 3, Fig. 12). Therefore, after a sufficiently long time, elementary steps merge to form a single macrostep (Figs. 4(b) and (c)).
It is noted that the morphology of the surface also freezes for |∆µ| < ∆µ f (τ, L). n is strongly dependent on the initial configuration of the surface (Fig. 4(a), Figs. 6(a) and (b)).

Multi-nucleations
From Fig. 9(b) of the growth/recession rate of the surface for ∆µ co < |∆µ| < ∆µ R (L), the slope of the line changes at a crossover point, ∆µ co . In addition, n is no longer constant but decreases as 1/|∆µ| decreases. In this subsection, we will apply the classical multiple heterogeneous nucleation theory to the Monte Carlo results.
We assume some critical nuclei arise at the edge of a macrostep with a mean equal distance l d (Fig. 10). The step-detachment time t d , is then approximated as where v t is the step unzipping velocity. From Eq. (17), l d is expressed using v t and I n as By substituting Eq. (18) into Eq. (17), we obtain The explicit values of l d and t d are calculated using Eqs. (18) and (19) with g * being g * Ising and Z/( √ 2C) being 4.5 × 10 −3 , where v t is assumed to equal v step,RSOS (|∆µ|) ( §IV B 3, Eq. (26)), which are given in Table II.
The growth/recession rate of the surface is then obtained by From v t ∝ ∆µ, we expect approximately for ∆µ co (L) < |∆µ| < ∆µ R (L). In Fig. 9(b), the slope of the Monte Carlo results in this region seems to be smaller than the slope for |∆µ| < ∆µ co (L). However, the Monte Carlo results for |V | bend around ∆µ co < |∆µ| < ∆µ R (L). It is interesting that t d < ∼ t n for ∆µ co (L) < |∆µ|. The 2D nuclei are randomly formed relatively often at the edge of the macrostep. The elementary step formed by 2D nucleation advances/recedes by detachment from the macrostep. After a certain time, the elementary step may be pulled back to the facet edge with some probability. Therefore, for ∆µ co (L) < |∆µ| < ∆µ R (L), the detachment of an elementary step from the edge of the macrostep limits the growth/recession rate of the vicinal surface.
In this manner, classical 2D heterogeneous multinucleation is found to explain the phenomena roughly. However, with a slight modification, agreements between the Monte Carlo results and the expressions based on classical 2D heterogeneous multi-nucleation with respect to |V | and n are significantly improved, as demonstrated in the following sub-subsections.

Surface slope of the "terrace"
Let us call the surface that contacts the faceted macrostep the "terrace". At equilibrium in the stepfaceting zone, the "terrace" is the (001) surface, which ∆µ/ε characterizes the profile of the faceted macrostep [26]. In the case of a non-equilibrium steady-state, the change of the "terrace" slope changes the dynamics of the vicinal surface. In this sub-subsection, we explain how the surface slope p 1 , is connected to the size of the merged step n , the yielding point ∆µ y (L), and the crossover point to a kinetically roughened surface, ∆µ R (L).
For every t d , an elementary step is detached from the macrostep at the edge of the macrostep (Fig. 10). The detached elementary steps form a vicinal surface with the slope p 1 , which contacts the macrostep. After several calculations (Appendix C), the surface slope is described by n , as follows: [69]: where N m is the number of macrosteps in the simulated system. Using Eq. (22) with the assumption N m = 1.75, p 1 is calculated from n obtained from the Monte Carlo calculation ( Fig. 11(a)). Keeping Eq. (21) in mind, we present ln[p 1 |∆µ/ǫ|] vs. 1/|∆µ/ǫ| (Fig. 11(b)). Here, the Monte Carlo results are not straight lines. In addition, for small |∆µ|, the Monte Carlo results reveal the size dependence.
Here, let us introduce ∆µ y (L) so that the Monte Carlo results are well reproduced by a straight line (Fig. 11(c)). For µ co (L) < |∆µ| < µ R (L), the best fitted line is obtained as ∆µ y (L) for respective L are shown in Table I. It should be noted that g * p is very close to g * Ising . The lines for Eq. (23) with L being 240 √ 2a are shown by the dark solid lines in Fig. 11(b). The lines for Eq. (23) with L being 160 √ 2a and 80 √ 2a are shown as pale solid lines in Fig. 11(b). Although the only modification is the introduction of ∆µ y (L), the lines for Eq. (23) reproduce the p 1 based on the Monte Carlo results quite well for ∆µ co < |∆µ| < ∆µ R (L).
For large |∆µ|, p 1 based on the Monte Carlo results depart from Eq. (23). The departing point |∆µ|/ǫ = 0.12, agrees well with ∆µ R (240 √ 2a) in Table I. For ∆µ y (L) < |∆µ|, the "terrace" has a slope p 1 , so that the profile of the faceted macrostep becomes similar to that in the step droplet zone at equilibrium. That is, the characteristic profile of the faceted macrostep in the step-faceting zone changes to that in the step droplet zone. This ∆µ y (L) indicates a yielding point with respect to the self-detachment of steps from the macrostep.
It is interesting that p 1 is singular at the point ∆µ y (L). ∆µ y (L) is a candidate for the non-equilibrium phase transition point. However, ∆µ y (L) < ∆µ f (τ, L) in the present study, so that the vicinal surface freezes around the yielding point. Therefore, phenomena regarding ∆µ y (L) were not realized.

Step velocity
In the steady-state, n + = n − = V /a and n + = p 1 v step /a, where a is the height of the elementary step. Therefore, we obtain another key quantity v step calculated with v step = V /p 1 (24) using the Monte Carlo results for V and n . The step velocity v step is approximately proportional to ∆µ for ∆µ co (L) < |∆µ|. Figure 12 shows the ∆µ dependence of v step /|∆µ/ǫ|. v step based on the Monte Carlo results can be fit by the following equations: These lines are shown in Fig. 12. The reason for the steep decrease in the step velocity as |∆µ| increases is the meeting of steps, which inhibits the growth/recession of steps substantially [22]. Figures 4(d) and (e) show that the detached steps meet at several sites on the surface due to thermal fluctuations. The sticky character of steps merges these steps locally. The step velocity of the merged steps is substantially small [22,23], so that the merged steps pin the growth/recession of the steps. The density of steps becomes larger as |∆µ| increases because t d becomes shorter. The steps then meet and are pinned more frequently for large |∆µ|. Therefore, the step velocity becomes smaller as |∆µ| increases. It is interesting that the phases of the waves on the detached steps (meandering) are often coherent, although the surface does not contain dislocations or impurities. Besides, the present model does not take the surface diffusion into account. Nevertheless, the phases of the waves on the detached steps appear coherent. This is because the advance/recession of the embryo formed at the edge of the macrostep is blocked more often by the preceding step when the location of the embryo is nearer to the unzipping point of the preceding step. The embryo formed at almost the center between the two unzipping points of the preceding step is more likely to survive.
For |∆µ| < ∆µ co , l d exceeds the linear size of the system, L (Table II). Therefore, the successive nucleation at the step edge breaks. The surface moves intermittently through 2D heterogeneous nucleation at the step edge, which is consistent with the interpretation in the previous subsection ( §IV A 1). In this region, V /p 1 does not indicate v step , because n − < v step p 1 . v step in this region should be the lesser of the two values obtained by Eqs. (25) and (26).

C. Consistency
Using the equations for p 1 and v step (Eqs. (23)-(25)), we can reproduce n and V for ∆µ co (L) < |∆µ| < ∆µ R (L). V is then expressed as follows: Equation (27) In the limit of p 1 → 0, n converges to N step /N m . The N step /N m values also reproduce the constant values of n for ∆µ f (τ, L) < |∆µ| < ∆µ co (L). In the case of p 1 = 0, z is expressed by using Eqs. (23) and (28). The lines of n are shown in Figs. 6 and 9(a). These lines also reproduce the Monte Carlo results well.

D. Kinetic roughening
For ∆µ R (L) < |∆µ|, the vicinal surface is kinetically roughened and the faceted macrostep disappears (Fig. 4(f)). Although there is no large-scale macrostep, the inhomogeneous bumpy structure remains on the surface. This bumpy structure is formed by thermal noise (Fig. 4(f)).
Here, the choice of β as the symbol for the exponent is in accordance with Ref. [70]. The Monte Carlo results for all sizes agree with the two lines of Eq. (27) or Eq. (31), as shown in Fig. 13(b). From the cross point of the two lines, ∆µ R (L)/ǫ − ∆µ co (L) is determined as 0.071 ± 0.005. For L = 240 √ 2a, we have ∆µ R (240 √ 2a)/ǫ = 0.121 ± 0.012. The step velocity is obtained from the Monte Carlo results in Fig. 12. The results are fitted to the following equation: as represented by the pale solid line in Fig. 12. Using Eqs. (32) and (31), we obtain p 1 from p 1 = V /v step : which is shown by pale solid lines in Fig. 11. The line reproduces the Monte Carlo results for p 1 well. The crossover from the vicinal surface with the faceted macrostep to the kinetically roughened surface is essentially caused by the change of the |∆µ| dependence of p 1 (Fig. 11(a)). By definition (Appendix C), p 1 indicates the density of the elementary steps. The meeting of steps occurs more frequently when the step density is larger. In Eq. (23), the merging of steps is not taken into consideration; Hence, the increase of the locally merged steps contributes to a decrease in p 1 ; i.e., it changes the |∆µ| dependence of p 1 .

V. DISCUSSION
Near equilibrium, the finite size effect is prominent. The question then arises, what happens with the infinite system size? The underline in the Table II shows the border value that exceeds the system size or the observation time in the present study. From Eq. (11), the waiting time t n for 2D heterogeneous nucleation converges to zero as L → ∞. We then have lim L→∞ ∆µ f (τ, L) = 0. However, t n increases so rapidly as |∆µ| decreases that the system size should be approximately 10 14 or more for t n < 10 8 with |∆µ|/ǫ = 0.01. Therefore, in an actual system with a length of ca. 1 mm, the non-negligible frozen region with respect to the surface growth/recession remains near equilibrium.
For l d and t d , there is no explicit size dependence. The waiting time t d , for the step detachment seems to be linked to ∆µ y (L). The Monte Carlo results in this study show a slight L dependence on ∆µ y (L). If ∆µ y (L) converges to a non-zero value ∆µ y (∞) in the limit of infinite system size, then the point may be a nonequilibrium phase transition point. However, to clarify whether ∆µ y (∞) is finite or zero is a future problem.
In our previous work [69], we studied the ∆µ dependence of the size for faceted merged steps n and the growth rate of the vicinal surface V in the step droplet zone I for the non-equilibrium steady-state. Some results are similar to the results of this study. n decreases as |∆µ| increases. In addition, for |∆µ| > ∆µ R (L), the vicinal surface with faceted macrosteps crosses over to the kinetically roughened surface without a macrostep. In the results of the previous study, the freezing region lacks a near equilibrium state, and the elementary step self-detaches by thermal noise without 2D nucleation processes. The morphology of the kinetically roughened surface is somehow different from that of the present study. The bumpy structures on the vicinal surface that were obtained in the previous study were so small that they cannot be discerned in the images of the simulated surface without magnification of the images [69]. Therefore, the scaling behavior in the kinetically roughened surface is slightly different from that in the present study.
It is interesting that the figure for V (Fig. 7) is analogous to the Fig. 4(b) in Ref. [70], which shows the velocity of the particles where plastic depinning occurs. The system has a depinning threshold, V ∝ (F D − F c ) β with β = 1.5, where V is the average velocity of the particles, F D is the driving force, and F c is the depinning threshold. The plasticity is said to be relevant to the charge-density wave systems [70,71]. As noted in §IV B 3 and §IV C, the elementary steps that detach from the faceted macrostep meet neighboring steps due to thermal noise. The steps are sticky, so that they merge at the meeting point. The locally merged steps with substantially low velocity then pin the motion of the elementary steps. Therefore, the step attachment-detachment motion is analogous to the motion of particles with plastic depinning. This begs the question, is there is a common mathematical framework? However, this question has yet to be answered.

VI. CONCLUSIONS
• Steps on the vicinal surface self-assemble to form faceted macrosteps in the steady-state for |∆µ| < ∆µ R (L).
• For ∆µ co (L) < |∆µ| < ∆µ R (L), the vicinal surface grows/recedes in the manner of attachment-detachment of steps at the macrostep edge (Figs. 4(d) and (e), Fig. 10). The attachmentdetachment of steps is understood based on successive classical 2D heterogeneous multi-nucleation at the edge of the macrostep. The absolute value of the step velocity |v step |, the absolute value of the surface velocity |V |, and the slope of the "terrace" p 1 increase with |∆µ|, according to Eqs. (25), (27), and (23), respectively. The characteristic length l d , the step-detachment time t d , and the average height of the merged step n decrease as |∆µ| increases, according to Eqs. (18), (19), and (28), respectively.
• The finite size and the finite time effects are distinctive for |∆µ| < ∼ ∆µ co (L).
Therefore, we have λ 2 = ∆µ 1 /∆µ 2 . ∆µ can be selected arbitrarily, so that ∆µ can be one, i.e., where N 1 is the number of single steps on the surface in contact with the (111) surface, and N m is the number of macrosteps. At the temperature k B T /ǫ = 0.4, we assume N m ≈ 1.75. Next, let us introduce z and x so that where N macro is the number of elementary steps that compose macrosteps, L 1 is the linear length of the "terrace", and L macro is the linear length of the macrosteps. N 1 = zpL/a = xp 1 L/d, so that and N macro = (1 − z)pL/a = √ 2(1 − x)L/a, so that From Eqs. (C1)-(C4), we obtain Eq. (22).