Surface Roughness Changes Induced by Stoichiometric Deviation in Ambient Phase for Two-Component Semiconductor Crystals

The effects of a deviation in the fraction of the components in the ambient phase of a stoichiometric AB compound, such as GaN or SiC crystals, on the surface roughness and step self-assembly and disassembly on a vicinal surface are studied using the Monte Carlo method based on a staggered restricted solid-on-solid (st-RSOS) model at equilibrium. The (001) and (111) surfaces are typical examples of non-polar and polar surfaces, respectively. Although a stoichiometric deviation of the ambient phase does not affect the surface energy of a non-polar surface, it affects that of polar surfaces such as the 
 
 
 
 
 ( 
 111 
 ) 
 
 A 
 
 
 
 and 
 
 
 
 
 ( 
 111 
 ) 
 
 B 
 
 
 
 surfaces. We found that the vicinal surface of an AB compound is atomically smooth and globally rough. Globally, the vicinal surface is not affected by a stoichiometric deviation in the ambient phase. In contrast, in a small area, the structure of the vicinal surface is affected by a stoichiometric deviation in the ambient phase. The vicinal surface consists of local double and quadruple steps. The characteristic length 
 
 
 
 L 
 
 M 
 F 
 L 
 
 
 
 
 , which separates the global length scale region and the local length scale region, has a maximum value of 
 
 
 
 156 
 a 
 
 
 
 in the present study, where a is the lattice constant. When temperature decreases, 
 
 
 
 L 
 
 M 
 F 
 L 
 
 
 
 
 can become large.


Introduction
Stoichiometric two-component semiconductors such as SiC, GaN, and GaAs (called AB compounds) have potential applications in power devices. In the fabrication of bulk semiconductor crystals, the formation of macrosteps on the crystal surface hinders the preparation of good-quality crystals, such as 4H-SiC, which are required for electrical devices [1].
In our previous studies, we showed that discontinuous surface tension, where surface tension is the surface free energy per unit area, causes faceted macrosteps at equilibrium [2][3][4][5][6][7][8][9]. We studied the restricted solid-on-solid (RSOS) model with point-contact-type step-step attraction (p-RSOS model). This attraction, which is the microscopic origin of the discontinuous surface tension, is considered to result from quantum interactions between steps at step meeting points [10]. Our non-equilibrium steady-state study [11][12][13] showed that the height of the faceted macrosteps decreases as the degree of non-equilibrium increases. At the same time, the vicinal surface becomes rougher due to the detachment of elementary steps from the faceted macrosteps.
Macrosteps can also form on a vicinal surface when a (7 × 7) Si(111) surface coexists with the (1 × 1) surface [14]. When the slope of the vicinal surface is large, the surface free energy of the (1 × 1) Si(111) surface has lower free energy than that of the (7 × 7) Si(111) surface. The (1 × 1) side surface of a macrostep coexists with the (7 × 7) terrace surface.
Here, the surface roughness in the statistical mechanics is defined in terms of a surface width W, given by where h(x, y) is the height of the surface at a site (x, y) and · is the thermal average. In the limit of L → ∞, where L is the linear size of the system, the global surface roughness is defined as follows [16,33]: In addition, based on the scaling concept of growing surfaces [34], the roughness exponent α is defined by For α = 0, the surface width diverges logarithmically. Macrostep formation occurs on the vicinal surface. For the original RSOS model, the vicinal surface at temperatures sufficiently less than T R (001) is well described by the TSK model; that is, the vicinal surface consists of smooth terrace surfaces and an approximately parallel train of steps with entropic repulsion. The entropic repulsion is caused by a zig-zag structure on the edge of single steps, which is expressed by the 1D free fermion's random walk [29][30][31][32]. These characteristics of 1D free fermions of elementary steps lead to GMPT behavior on the vicinal surface; that is, the surface free energy per x − y area f (p) has the following GMPT form: where f (0) is the free energy of the terrace surface, γ(φ) is the step tension (free energy per unit length), a is the unit step height, p is the surface gradient, and B(φ) is the step interaction coefficient. Here, φ is the tilt angle in the mean step-running direction relative to the y-axis in the terrace plane. The vicinal surface corresponds to the tangential plane at a point on the curved area of an equilibrium crystal shape (ECS), which is the shape with the lowest total surface free energy [22,26,35], around the (001) or (111) facet. The squared surface width of the tangential plane is logarithmically divergent with respect to the linear system size L. This logarithmic divergence of the squared surface width also indicates the roughness exponent α = 0. In the limit of L → ∞, the divergence results from surface slope fluctuations with small wave numbers caused by step wandering on the vicinal surface [33,35]. W 2 diverges, such as W 2 = (1/2π 2 ) ln L, in the limit of the surface slope p → 0 [30,31,36,37]. Hence, the vicinal surface of the original RSOS model is rough according to the definition in Equation (2) at all temperatures. Although the vicinal surface is rough, the ratio W/L ∼ ln L/L converges to zero in the limit of L → ∞. Therefore, the vicinal surface can be considered to be a thin flexible film on the mesoscopic scale.
For an AB compound, the polar A surface has lower energy than that of the B surface if the ambient phase is A-rich [38][39][40]. In this circumstance, macrostep formation may occur due to component deviation in the ambient phase. The deviation of the components in the ambient phase has been experimentally studied. It is known that a small deviation in the components in the ambient phase leads to good-quality crystals in crystal growth from the vapor phase [41]. However, for crystal growth of SiC from Si solution, the concentration of C is extremely low because of the poor solubility of carbon in silicon [1]. The deviation of the components in the ambient phase may be responsible for macrostep formation. Surface roughness measurement can be used to determine whether macrosteps are faceted.
The aim of this paper is to study how component deviation affects surface step assembly and disassembly and surface roughness. We adopt a staggered restricted solid-on-solid (st-RSOS) model [42][43][44] for the Monte Carlo method.
The basic quantity of the st-RSOS model is the microscopic ledge energy and the bulk chemical potentials of the ambient and crystal phases, namely µ ambient and µ crys , respectively. The st-RSOS model is a slightly more coarse-grained model compared to the atomic model for the first-principles density functional theory (DFT) method [45][46][47][48][49][50]. An "atom" in the st-RSOS model is about a half-unit cell in the DFT model. The and µ crys of the st-RSOS model should be obtained from the surface free energy per half-unit cell of the DFT model. Hence, and µ crys depend on temperature. Throughout this paper, is assumed to be constant because the main focus is the effects caused by stoichiometric deviation in the ambient phase. We will return to this in section 2.1.
Because surface phenomena are complex, we adopt the strategy of divide-and-conquer to study them step-by-step. To clearly observe the effects of the change due to component deviation in the ambient phase at equilibrium, we do not consider the second-nearest-neighbor interaction between atoms [51] in crystals, the Ehrlich-Schwoebel effect [52,53], elastic interactions [54], surface reconstruction [14], and point-contact-type step-step attraction [2].
It should be noted that there are other concepts associated with surface roughness, such as surface diffuseness and sharpness. The 4 He solid-super fluid interface [55][56][57][58][59] is an example of a diffuse interface. The facets on the ECS of a 4 He crystal droplet are an example of an atomically rough [60] but globally smooth surface. Surface diffuseness or sharpness is considered to be related to the crystal-melt interface or surface melting phenomena [61,62].
The rest of this paper is organized as follows. In Section 2, the surface model and chemical potentials for AB compounds are briefly explained. The Monte Carlo method is also explained. In Section 3, the Monte Carlo results concerning the behavior of the vicinal surface in a small area are shown and discussed. In Section 4, the Monte Carlo results concerning the global behavior of the vicinal surface are shown and discussed. The border between local behavior and global behavior with respect to vicinal surface roughness is also discussed.

Surface Hamiltonian
An AB compound is made of A and B atoms in the ambient phase as follows: Let us define µ A,ambient (µ B,ambient ) to be the bulk chemical potential of A atoms (B atoms) in the ambient phase. Let us consider an NaCl-type crystal structure ( Figure 1) of an AB compound, where µ A,crys (µ B,crys ) is the bulk chemical potential of A atoms (B atoms) in the AB compound (Appendix A). The {111} and {001} surfaces are the polar and non-polar surfaces, respectively. Here, the chemical potential of the AB compound µ AB,crys is defined by [38,40,47] µ AB,crys = µ A,crys + µ B,crys . To study the surface roughness of a vicinal surface, the st-RSOS model [42], which is a staggered 19-vertex model [43,44], is adopted for the microscopic model in the Monte Carlo method. Here, "restricted" means that the height difference between nearest-neighbor sites is restricted to {0, ±1}. The surface Hamiltonian is given by the following equation: where h(m, n, σ) is the height of the surface at a site (n, m) with atomic species σ = A or B, is half the lateral effective bond energy between nearest-neighbor atoms (microscopic ledge energy), N is the total number of unit cells on the (001) surface, and E A,surf (E B,surf ) is the surface energy at a site for A atoms (B atoms) per half-unit cell. The RSOS condition is implicitly required. Here, ∆µ A , ∆µ B , and ∆µ are introduced as follows: Here, the number of atoms in the crystal or ambient phase is not fixed. The atoms in the ambient phase are captured by the vicinal surface; the atoms on the outermost layer of the crystal can escape to the ambient phase. Hence, in the statistical system under consideration, a grand canonical ensemble is adopted. The third and fourth (fifth) terms on the right-hand side of Equation (7) indicate the energy gain that an A-(B-) atom transfers from the ambient phase to the crystal. When h(m, n) + m + n is even (odd), the transferred atom is A (B). The partition function in Equation (7) Z (T, ∆µ A , ∆µ B ) is the grand partition function, and the thermodynamic potential Φ(T, ∆µ A , ∆µ Because an AB compound is studied, the equilibrium condition at two-phase coexistence is given by [38][39][40]63,64] µ A,crys + µ B,crys = µ A,ambient + µ B,ambient . From Equations (8) and (9), we obtain the following Gibbs-Duhem equation: Here, the steps covered by A atoms (B atoms) are referred to as A steps (B steps). For ∆µ A > 0, an A step has lower energy than a B step because of the fourth and fifth terms on the right-hand side of Equation (7). Because the (111) surface is a close-packed arrangement of steps, the (111) A surface has lower energy than that of the (111) B surface. This has been experimentally confirmed by Nishizawa and Pons [41]. Based on the energy Hamiltonian in Equation (7), the surface roughness for a vicinal surface tilted towards the 111 direction is studied.
When the ambient phase is a gas, ∆µ A and ∆µ B are expressed as follows assuming an ideal gas (Appendix A): where P A (P B ) is the partial pressure of A atoms (B atoms) in the bulk gas phase. Instead of using µ A,crys or µ B,crys , µ A,gas and µ B,gas are used so that P A0 (P B0 ) is the partial pressure determined by the condition in Equation (12) (Appendix A). Hence, we obtain Equation (11).
It should be noted that the dominant contribution to ∆µ A or ∆µ B is the mixed entropy in the ambient phase. Irrespective of the atomic species, the thermodynamic effect of the concentration deviation in the ambient phase can be determined from the change in ∆µ A or ∆µ B . As an example, when the ambient phase is a solution, ∆µ A and ∆µ B are expressed as follows assuming an ideal solution: where c A (c B ) is the concentration of A atoms (B atoms) and c A0 (c B0 ) is their concentration under the condition in Equation (12). The bulk chemical potentials and the surface energies for the flat (001) surface should be calculated using the first-principles DFT method [45][46][47][48][49][50]. However, the surface model expressed by the Hamiltonian in Equation (7) is slightly more coarse-grained compared to the surface structure considered in the first-principles DFT calculations. In the model of Equation (7), or E A,surf (E B,surf ) is considered to correspond to the surface free energy per half-unit cell in the DFT calculations. Hence, µ A,crys (µ B,crys ) ≈ z with z = 6, where z is the number of nearest-neighbor sites in the bulk crystal, E A,surf ≈ E B,surf ≈ . According to a recent study by Kempisty et al. [50], the surface free energy decreases as the temperature increases due to the entropy caused by distortions or vibrations. Hence, , P A0 , P B0 , c A0 , and c B0 depend on temperature.

Monte Carlo Method
To study the effects of a stoichiometric deviation in the ambient phase, the Monte Carlo method for non-conserved systems is used. Here, "non-conserved" indicates that the number of crystals is not conserved. Atoms on the crystal surface can escape to the ambient phase or be captured from the ambient phase. Hence, surface and volume diffusion are not taken into consideration.
The vicinal surface around the (001) surface tilted towards the (111) direction is considered using the Monte Carlo method with the Metropolis algorithm. The external parameters are temperature T, ∆µ A , number of steps N step , and the linear size of the system L.
Initially, N step steps run in the mean directionỹ = 1 10 , for which the periodic boundary condition is required. Thex direction is assigned to the 110 direction. The configurations on the lower side (right side of the top-down view of the surface, e.g., Figure 2a) are connected to the upper side (left side of the top-down view of the surface) by adding N step a. We consider two types of initial surface configuration: a surface with a macrostep where all steps assemble, and a surface with a train of steps.
The lattice sites on the surface at which an event occurs are selected randomly. "Capture" or "escape" for the lattice site is selected randomly with a probability of 1/2. The energy change between the states before and after an event is calculated using Equation (7). If the energy change is negative, the surface configuration is updated with a probability of 1. If the energy change is positive, the surface configuration is updated with a probability of exp[−(E after − E before )/k B T]. Here, E after (E before ) is the surface energy calculated using Equation (7) after (before) the surface configuration update.
The capture or escape of atoms on the (001) terrace surface occurs at equilibrium. Because surface diffusion of atoms is not taken into consideration, the advance and recession of an elementary step are, respectively, caused by the capture and escape of atoms at kinks on step edges. Snapshots of the top-down and side views of the simulated surface at 4 × 10 8 Monte Carlo steps per site (MCS/site) are shown in Figure 2. It should be noted that the elementary steps have few kinks on step edges for ∆µ A / = 2. We will return to this point in the following sections.
The mean surface profile normal to the mean running direction of elementary steps is calculated as follows: where theỹ andx directions are, respectively, parallel and normal to the mean running direction of elementary steps. In the present study,ỹ andx are, respectively, in the 1 10 and 110 directions .
In Figure 3a, an example of a mean surface profile is shown with a size of 320 √ 2a (a = 1) at 4 × 10 8 MCS/site. The homogeneous slope of the surface can be seen.  (14)). t = 4 × 10 8 MCS/site. (b) Dark full lines:x dependence of gW 2 (x, L) (Equation (16)). Light thick line: gW 2 (L) obtained using the least squares method. ∆µ A = 0. (c) The relationship between W and √ gW.
h (x, t) and gW 2 (L) are averaged over 2 × 10 8 MCS/site after initial 2 × 10 8 MCS/site. The mean surface heighth is calculated as follows: In Figure 4, time evolutions of the mean surface height are shown for several ∆µ A values. The fluctuations of the location of the surface are relatively large for ∆µ A / = 0 and 4. As shown, there are long-time fluctuations around a mean value. We will return to this point in Section 3.3. However, there is no systematic increase or decrease in surface height. Hence, the systems are at equilibrium as expected.
The square of the surface width W 2 (x, L) is defined by the variance of the surface height, which is calculated as where g is a metric for a tilted surface [35], θ is the angle of the tilt towards the 111 direction, τ 0 = 2 × 10 8 and τ max = 4 × 10 8 MCS/site are the averaged times, and p x and p y are the surface gradients.
In the present Monte Carlo simulations,p An example of the calculated gW 2 (x, L) is shown in Figure 3b. Although the surface is tilted, the square of the surface width is constant with respect tox. We determine gW 2 (L) as the averaged gW 2 (x, L), fitted using the least squares method.

Mean Height of Local Merged Steps
To study whether a macrostep is self-organized, we calculate the mean height of merged steps n as follows: where n step (ỹ) is the number of macrosteps atỹ along thex direction. n is averaged over from τ 0 to τ max . The obtained n is shown in Figure 5a. The values of n with respective ∆µ A agree with the values of n with −∆µ A within the errors. Near ∆µ A ≈ 0, n is about 1.1, which means that the elementary steps are well separated. This can be seen in Figure 2a.
Around ∆µ A / ≈ 2, however, there is a broad peak that diminishes when the size of the system increases. Double steps covered by A atoms appear in Figure 2b. There are few kinks on the edges of the elementary steps. Because the peak disappears in large systems, it is considered to be caused by a finite size effect. The zig-zag structure of steps causes entropic repulsion between steps. The small number of kinks in the small system weakens this repulsive interaction between steps. In addition, once neighboring steps in the small system arrive at a close-packed location, kink formation becomes difficult due to the RSOS geometrical restriction.  (17)). Averaged over 1 × 10 8 MCS/site. (b) ∆µ A dependence of the double and quadruple step ratios C 2 and C 4 , respectively (Equations (18) and (19)). k B T/ = 0.2.p = 0.530.
For 1 ≤ ∆µ A / ≤ 3, n is about 1.25 for a large system. The fraction of double steps c is estimated using the following equation: We have C 2 = 0.25 for n = 1.25 (Figure 5b). For a small system, n has a maximum value of 1.7 near ∆µ A / ≈ 2. We have C 2 = 0.7. This means that 70% of elementary steps form double steps for a small system.
For ∆µ A / ≥ 4, n increases to 2.5 as ∆µ A increases. A snapshot of the surface for ∆µ A / = 5.0 is shown in Figure 2d. Almost the whole vicinal surface is covered by A atoms. Assuming that the vicinal surface consists of double steps and quadruple steps, we have where C 4 is the fraction of quadruple steps (Figure 5b).

No Step Faceting
As described in the previous subsection, double and quadruple merged steps are observed locally. However, macrosteps that run from edge to edge are not observed. For a macrostep to form on a vicinal surface at equilibrium, a first-order transition on the surface free energy with respect to the surface orientation is required [3,65]. This can be seen in the Wulff figure, which is the polar graph of surface tension [66][67][68].
At finite temperature for the original RSOS model, the surface tension γ surf decreases owing to surface entropy per area s surf as the temperature increases, because γ surf = E surf − Ts surf [5]. The (001) surface is thermally roughened at the roughening transition temperature T R (001) [15,20]. The order of the roughening transition is the Kosterlitz-Thouless-type continuous transition [20]. The T R of the original RSOS model was calculated by den Nijs to be approximately k B T R (001) / = 1.580 [69].
The step tension on the (001) surface decreases to zero at T R (001) as the temperature increases [20]. The (001) facets on the ECS shrink with increasing temperature, vanishing at T R (001) . This shape transition is called the faceting transition [22,23,26]. The temperature dependence of the (001) facet-shapes on the ECS has been studied by several researchers [22,70,71]. For the st-RSOS model, the (001) surface and the ECS also cause roughening and faceting transitions at T R (001) , respectively. In addition, T R (001) depends on the component deviation in the ambient phase ∆µ A . We first study the stoichiometric change in the Wulff figures at T = 0K. Cross-sections of the Wulff figures obtained at T = 0 K for the st-RSOS model are shown in Figure 6. The figures show how the change in relative chemical potential in the ambient phase affects the surface tension for ∆µ A ≥ 0. The ECS is also shown. At ∆µ = 0, the surface tension of the A surface equals that of the B surface and the ECS for the (111)A surface agrees with that for the (111)B surface. For 0 < ∆µ A / < 4, the surface tension of the A surface (dark blue lines) is lower than that of the B surface (light blue lines). Because the surface tension of the A surface is smaller than that of the B surface, the ECS with the (111)A surface (red full lines) is realized. For ∆µ A / ≥ 4, the (001) surface does not appear because it is roughened stoichiometrically. At ∆µ A / = 6, the surface energy of the (111) surface becomes zero. For ∆µ A / > 6, the surface energy of the (111) surface is negative, and thus, the system becomes unphysical. Because surfaces and steps are low-dimensional objects, thermal fluctuations are often sufficient to destroy the ordered phase at finite temperature in a system based on short-range forces [72]. The calculation of the surface free energy using a mean field or quasi-chemical approximation for low-dimensional objects often leads to incorrect results. Hence, calculations more precise than those obtained from mean field approximation are required. However, such calculations are quite difficult for the st-RSOS model.
At temperatures sufficiently lower than T R (001) , the (001) surface is well approximated by a 2D lattice gas model [73,74]. The lattice gas model causes a second-order 2D Ising-type phase transition at T c . Although T R (001) is known to be higher than T c , we can calculate the ∆µ A dependence of T c using the imaginary-path-weight random walk (IPW) method [75][76][77][78][79] based on a 2D lattice gas model (Appendix B). The ∆µ A dependence of T c is shown by the broken line in Fig. 13(a) in Ref. [40] (H/J = ∆µ A / ) as a critical curve under the assumption that is constant. The temperature on the critical curve at ∆µ A = 0 equals k B T c / = 1/ ln(1 + √ 2) ≈ 1.13459 (= k B T c /(2J)). ∆µ Ac / on the critical curve at T = 0 and k B T/ = 0.2 equals 4 and 3.7227, respectively [40].
The ∆µ A dependence of facet shapes at k B T/ = 0.2 is shown in Figure 7 [40,74]. As seen from the figures, the anisotropy of the surface tension and the facet shape depends on ∆µ A . The area of the facet shape shrinks as ∆µ A increases, becoming zero at a temperature k B T/ = 0.2 with ∆µ Ac / = 3.7227. The morphology of the surface obeys not only the symmetry of the energy Hamiltonian but also the geometrical symmetry of the crystal structure. This geometrical situation causes non-monotonic changes of gW 2 (L), n , and step quantities with respect to ∆µ A . The formation of macrosteps at equilibrium requires a first-order shape transition [23], which leads to discontinuous surface tension [3,8]. However, the Wulff figures do not show such a discontinuity in the present case. To obtain discontinuous surface tension at low temperature, the point-contact-type step-step attraction that originates from the quantum mechanical interactions between steps [10] is required. The present Monte Carlo study does not produce macrosteps. Therefore, a stoichiometric deviation in the ambient phase does not cause global macrosteps.

∆µ A Dependence of Surface Roughness
The obtained ∆µ A dependence of gW 2 (L) for the st-RSOS model is shown in Figure 8. The values obtained for different initial configurations agree well with each other. In contrast to n , gW 2 (L) has a minimum at around ∆µ A / = 2. Regarding the roughness of the vicinal surface, this behavior of gW 2 (L) is consistent with long-term height fluctuation up to ∆µ A / ∼ 4 ( Figure 4). Because the number of atoms in the surface area is an internal parameter, it fluctuates. Although the (001) surface is smooth, the capture and escape of atoms (units) on the (001) surface occur at equilibrium. Wandering of the parallel train of steps [33,35] on the vicinal surface is caused by the atomic capture and escape processes at kinks on step edges. The change in the (001) surface roughness is thus not the dominant factor in the change in gW 2 (L). More importantly, the change in kink density on step edges and the change in the fraction of double and quadruple steps are mainly responsible for the gW 2 (L) change for |∆µ A / | < ∼ 4.  (Figure 2d), which are typical for a thermally roughened (001) surface. The intrinsic width of a step [80,81] is infinitely large for a thermally roughened surface.
However, the mean height fluctuation shown in Figure 4d is small. This is because the concentration of B atoms in the ambient phase is insufficient to grow or recede the surface.

Effect of Component Deviation on Roughness Exponent
The size dependence of gW 2 (L) for |∆µ A / | < 4 is shown in Figure 9. From Figure 9b, it can be clearly seen that gW 2 (L) increases logarithmically as L increases for large L. The logarithmic divergence of the surface width for high |∆µ|, such that |∆µ/ | = {4, 5}, is shown in Figure 10. Therefore, α = 0 for all ∆µ A . This behavior is consistent with previously reported results for thermally roughened surfaces (interfaces) [26,37]. These results indicate that a stoichiometric deviation in the ambient phase does not affect the roughness exponent for a large system.

Mean Free Length
According to the linear theory of molecular beam epitaxy (MBE), the time evolution of the surface can be expressed by the following equation [34]: where h is the surface height, ν is the surface tension, K is a constant related to surface diffusion, and η is a random force. The competition between the terms ∇ 2 h and ∇ 4 h generates a characteristic length L 1,MBE , which is expressed by [34] For L < L 1,MBE , ∇ 4 h is dominant. Hence, the roughness exponent α is (4 − d)/2, where d is the spatial dimension. For L > L 1,MBE , ∇ 2 h is dominant, and thus α is Because Equation (20) is derived from the symmetry principle [34], the linear theory with K = 0 is applicable to a crystal and solution system, such as SiC in Si fluid. In the present study, because d = 2 and K = 0, we have L 1,MBE = 0 from Equation (21). Hence, gW 2 (L) is expected to behave logarithmically for all L. However, as shown in Figure 9, there exists another characteristic length L 1 . Figure 9b, the calculated gW 2 (L) deviates from the logarithmic behavior on the short length scale. First, the Monte Carlo data of L ≥ 160 √ 2 in Figure 9b are fitted to the following equation:

As shown in
where A is the amplitude and L MFL is the mean free length defined in Equation (22). The obtained A and L MFL are listed in Table 1. L MFL has a maximum value at around ∆µ A / = 2.0. The ratio between L MFL and L MFL | ∆µ A =0 is shown in the fourth column in Table 1. As shown in Figure 2b, there are few kinks on a given step. To see what happens on the surface at ∆µ A / = 2.0, we consider a 2D lattice gas model for the two-component system. Table 1. Characteristic lengths, step tension, and step stiffness.

Amplitude
Step

Analysis Based on Two-Dimensional Lattice Gas Model
Step tension γ, which is the step free energy per unit length, and step stiffnessγ = γ + ∂ 2 γ/∂φ 2 , where φ is the tilt angle between the mean step running direction and the y direction ([010] direction), can be calculated based on the surface configuration energy in Equation (A11). Because a step is a 1D object, the entropy is significant in low dimensions [72]. Hence, entropy calculation that is more precise than mean field or quasi-chemical approximation is required.
The step tension and step stiffness are calculated precisely using the IPW method [40,[75][76][77][78][79] (Figures 7 and 11). The calculated step tension and step stiffness are shown in Table 1 in the fifth and sixth columns, respectively.
Due to the mixing of single and double steps for 1 ≤ ∆µ A / < 4, γ m andγ m are considered to be The ratios γ m /γ| ∆µ A =0 andγ m /γ| ∆µ A =0 are shown in the seventh and eighth columns of Table 1, respectively.
According to Equation (20), surface tension may explain the ∆µ A dependence of L MFL . The surface tension γ surf , which is the surface free energy per normal area, is expressed by γ surf ≈ γ m ρ, where ρ is the step density, for the lowest order with respect to ρ. However, the ∆µ A dependence of the surface tension or the step tension (seventh column in Table 1  The ∆µ A dependence of the ratioγ m /γ| ∆µ A =0 (eighth column in Table 1) is similar to that of L MFL /L MFL | ∆µ A =0 . In the following subsection, the relationship between step stiffness and the mean free length is examined.

Relationship between Surface Width and Step Width
For single steps, the step width w s,single (L) is estimated using a 1D free random walk such that w 2 s,single (L) = Lσ, where L is the mean running distance of a step along the step and σ is the scaled step width. In the 1D free random walk, L corresponds to time and σ corresponds to the diffusion constant. σ is related to the stiffness of a single step [82] as follows: where λ is related to the volume of a 2D island with an equilibrium shape and κ is the curvature of the edge of the island. For a vicinal surface, steps collide with neighboring steps. This step collision changes the squared fluctuation width of a single step from linear L to log L [36,37]. Recall the relationship between surface height-height correlation G(x,ỹ) and the step width ∆ 2 s,vicinal (ỹ)) [36,37] (Figure 12a) on a vicinal surface: where G(0,ỹ) is whereÃ is the amplitude. G(0,ỹ) = 2gW 2 in a large L limit. Then, let us introduce L MFL,g as a length, with which an individual step makes contact with neighboring steps (Figure 12b): where is the mean step distance. Substituting Equation (24) into Equation (27), we have Forỹ < L MFL,g , we have the following equations: Forỹ > L MFL,g , G(0,ỹ) increases logarithmically with respect toỹ, as shown in Equation (26).
we have Here, s is a parameter related to the mixture of merged steps. Equations (28), (30), and (31) well explain the ∆µ A dependence of L MFL and gW 2 (L) for small L.

Small Systems
When is assumed to be a (=1), L MFL,g expresses the mean distance between kinks on the step edge. We denote this quantity as L kink . From Equations (24) and (27), we have Around ∆µ A = 2, L kink estimated from the value of the step stiffness is 26a. This result is consistent with Figure 2b. When the system size is sufficiently small, there are no kinks on the step edge. In such a situation, the surface is similar to a smooth surface. Hence, gW 2 (L) is almost constant, rather than proportional to L, as shown by Equation (29). The deviation from the logarithmic behavior around ∆µ A / = 2 ( Figure 9) confirms this.
L kink can also be observed using the method proposed by Swartzentruber et al. [83][84][85]. For L < L kink , phenomena similar to 1D nucleation [86] may occur at an elementary step edge. Due to the entropy effect in low dimensions [72], 1D nucleation does not typically occur [87]. However, when the system size is sufficiently small (L < L kink ), 1D-nucleation-like phenomena may occur in a practical case of crystal growth.
When the temperature is higher than that in the present study, step stiffness decreases, becoming more isotropic. L MFL and L kink decrease according to Equations (28) and (32), respectively. When the temperature is lower, the anisotropy of step stiffness becomes stronger. Hence, L MFL around ∆µ/ ∼ 2 becomes substantially larger. This is consistent with the models used in the first-principles quantum mechanical calculation of the surface. In the quantum mechanical calculation, the model is assumed to contain several unit cells in a flat area. The surface is assumed to be smooth. However, on a mesoscopic scale, the vicinal surface is rough. Our results resolve this conflict and connect the quantities calculated using quantum mechanics to the hydrodynamic model.

Conclusions
For a stoichiometric two-component crystal, this study investigated the effects of the component deviation in the ambient phase on the assembly and disassembly of elementary steps on a vicinal surface using the Monte Carlo method based on the st-RSOS model.
The global behavior of the vicinal surface can be summarized as follows: • A component deviation in the ambient phase ∆µ A does not cause macrostep self-organization. • ∆µ A does not affect the roughness exponent; i. e., the roughness exponent α = 0. For L > L MFL , the squared surface width W 2 (L) logarithmically diverges with respect to L → ∞.
The behavior in a small area can be summarized as follows: • ∆µ A ( = 0) stabilizes polar surfaces by lowering the surface energy. This changes the morphology of the crystal through the anisotropy of the surface tension and step tension. This also causes non-monotonic changes of gW 2 , n , L MFL , σ, and step stiffness with respect to ∆µ A . • Double steps and quadruple steps form locally for |∆µ A / | 1 and |∆µ A / | 4, respectively. • For L kink < L < L MFL , W 2 (L) is proportional to L; for L < L kink , W 2 (L) is constant, making the surface similar to a smooth surface.

•
L MFL is about 160a, where a is the lattice constant, around |∆µ A / | = 2. At lower temperatures, L MFL becomes substantially longer.
The st-RSOS model is a slightly more coarse-grained model compared to that used for the first-principles DFT calculations. This study links the values obtained by first-principles DFT calculations to fundamental quantities required by coarse-grained continuous models on a mesoscopic scale.
The chemical potential of an ideal gas is expressed as follows [63,64]: where Z trans , Z rot , and Z vib are the partition functions for the molecular translation, molecular rotation, and molecular vibration, respectively, g e is the degree of degeneracy of the electron energy level, and P is the gas pressure. We can rewrite Equation (A1) as follows: where f gas (T) is the terms in Equation (A1) except for k B T ln P.
The chemical potential of the flat (001) crystal surface can be expressed as follows [49]: where T • , P • , µ • crys (P • , andT • ) are the standard temperature, standard pressure, and standard chemical potential, respectively. s and v are the crystal entropy and crystal volume per unit cell, respectively. The term µ crys (T, P) primarily depends on temperature and weakly depends on pressure [49]. We thus adopt µ crys (T) .
For an AB mixture, assuming an ideal gas, the free energy of the gas mixture is obtained using the following equations [64], which are similar to Equation (A1): Hence, we have 1 N A + N B G mixture (T, P A , P B ) = x A µ A,gas (P A , T) + x B µ B,gas (P B , T), where µ A,gas (P A , T) = −k B T ln (Z Atrans Z Arot Z Avib ) = k B T ln P A + f A (T), µ B,gas (P B , T) = −k B T ln (Z Btrans Z Brot Z Bvib ) = k B T ln P B + f B (T).
From the phase coexistence condition in Equation (12), we have the following equations, which are similar to Equation (A4): µ A,crys (T) = k B T ln P A0 + f A,gas (T), µ B,crys (T) = k B T ln P B0 + f B,gas (T).

Appendix B. Two-Dimensional Lattice Gas Model
The energy of the 2D lattice gas model of the configuration of the (001) surface is expressed as follows [40]: where is half of the chemical bond energy between A and B particles, A and B are the surface chemical potentials for a lattice gas particle on the A and B sub-lattices, respectively, and C A (x) = {0, 1} is the lattice gas variable at a lattice site x on the sub-lattice A . When C A (x) = 1, the sub-lattice site A at the lattice point x is occupied; when C A (x) = 0, it is empty. Similarly, C B (x) is assigned for a particle on the B sub-lattice site. The surface chemical potential is also calculated using simple broken bond counting. For the (001) surface, the surface chemical potentials are expressed by where z = 4 is the number of nearest-neighbor sites in the lateral plane. At equilibrium with two-phase coexistence, the energy with all sites occupied should be equal to that with all sites empty. Hence, we have the same equation as Equation (10): For the (001) surface at equilibrium, let us introduce anti-ferromagnetic Ising spin variables Substituting Equation (A14) into the lattice gas Hamiltonian in Equation (A11), at equilibrium, we have the spin Hamiltonian where J is the anti-ferromagnetic coupling constant between nearest-neighbor spins and H is the uniform magnetic field. Polar graphs of the interface stiffness calculated using the IPW method for the 2D lattice gas model are shown in Figure 6 [40].