Kardar–Parisi–Zhang roughening associated with nucleation-limited steady crystal growth

The roughness of crystal surfaces and the shape of crystals play important roles in multiscale phenomena. For example, the roughness of the crystal surface affects the frictional and optical properties of materials such as ice or silica. Theoretical studies on crystal surfaces based on the symmetry principle proposed that the growing surfaces of crystal growth could be classified in the universal class of Kardar–Parisi–Zhang (KPZ), but experiments rarely observe KPZ properties. To fill this the gap, extensive numerical calculations of the crystal growth rates and the surface roughness (surface width) have been performed for a nanoscale lattice model using the Monte Carlo method. The results indicate that a (001) surface is smooth within the single nucleation growth region. In contrast, the same surface is atomically smooth but thermodynamically rough in the poly-nucleation growth region in conjunction with a KPZ roughness exponent. Inclined surfaces are known to become Berezinskii–Kosterlitz–Thouless (BKT) rough surfaces both at and near equilibrium. The two types of steps associated with the (001) and (111) terraces were found to induce KPZ surface roughness, while the interplay between steps and multilayered islands promoted BKT roughness.

of kinetic roughening, and Cuppen et al. 23 , considering a non-equilibrium definition of kinetic roughening that also includes MC simulations for the nearest neighbor (nn) Kossel solid-on-solid (SOS) model.For so-called classical criteria 22 , �µ c is calculated based on the relationship πγ 2 s /(3k B T�µ c ) ∼ 1 25 , where γ s is the step ten- sion, which is the step free energy per unit length.
Meanwhile, in the field of statistical mechanics, the term kinetic roughening is used to refer to Kardar-Parisi-Zhang (KPZ) roughening 28 .For fluctuating surfaces (or interfaces), the Family-Vicsek scaling relationship [29][30][31][32][33][34][35] has been widely used to describe the self-affine surface.The Family-Vicsek scaling relationship for a surface can be expressed as where α , β and z are the roughness, growth and dynamic exponents, respectively.In the non-equilibrium steady state, the surface width is characterized by the roughness exponent α .Based on the symmetry principle, a surface growth equation including a nonlinear term obtained from the KPZ model 28 can be derived as where h t is the surface height at time t, v 0 is the constant surface velocity, ν > 0 is a coefficient related to surface tension, is a coefficient for the nonlinear term and η t is white noise in space and time.In the case of a 2D sur- face in a 3D system, the values of these exponents are predicted to be α = 0.3869 , β = 0.2398 and z = 1.6131 34,35 (KPZ-rough surface).The experimentally determined values of various systems such as directed polymers are known to agree with these exponents, indicating that these systems belong to the KPZ universality class.However, in the case of crystal growth, the observed roughness exponents tend to differ from those predicted by the KPZ model 30,35,36 , with the exception of several special surface systems 37,38 .
In our previous work 39 , we found crossover phenomena between BKT rough and KPZ rough surfaces at �µ cr for inclined surfaces using the Monte Carlo method based on the nn restricted solid-on-solid (RSOS) model on a square lattice.Here, the term "restricted" indicates that the surface height difference between nearest neighbor sites is limited to {0, ±1} .In the work, thermally excited structures such as adatoms, adholes, islands or negative islands (clusters of adholes) on terraces are found to cause the crossover phenomena.These thermally excited structures are thought to be irrelevant and they only renormalize the step tension for the thermal roughening transition at equilibrium.However, several Monte Carlo results on the surface width W suggested a complex surface slope dependence of W.
The aim of the present work is to clarify what makes the growing crystal surface KPZ rough using the RSOS model with and without surface steps.In addition, we obtain the detailed surface slope dependence of W explicitly.For this purpose, extensive numerical data on growth rate and surface roughness of planar or inclined surfaces were collected using the Monte Carlo method for a non-equilibrium steady state.
The RSOS model is more restricted than the Kossel SOS model studied by Van Veenendaal et al. 22 or the absolute SOS (ASOS) model 4,40 .However, the model has a special characteristic in that an almost ideal terracestep-kink (TSK) model 41,42 is realized 39,43 around a (111) surface due to the RSOS restriction.In the ideal TSK model, thermally excited structures are prohibited on the terrace surface 44,45 .
The RSOS model is a more microscopic model than the models used in phase field calculations 46 , whereas it is a coarse grained model for the purposes of first principles quantum mechanical calculations 47 .
It should also be noted that the RSOS model is equivalent to a 19-vertex model 48,49 and, because the latter represents a non-integrable system, the RSOS model cannot be solved exactly using the Bethe ansatz approach 50 .This is one of the reasons why the present work chose to study the RSOS model numerically.
At equilibrium, the RSOS model employed in the present work is equivalent to that previously used to determine roughness exponents by Amar and Family 51,52 .It is also important to note that a model used in the field of nonlinear dynamics with the restriction of the height difference being an integer is also sometimes referred to as the RSOS model but is known as the ASOS model in the field of roughening transition studies 4,40 .
During the present work, surface diffusion, volume diffusion and elastic effects were not taken into consideration.

RSOS model
The surface energy of a surface with an orientation close to (001) exhibiting (001) terrace roughness can be expressed by the discrete Hamiltonian 39 where h(n, m) is the surface height at site (n, m) on a square lattice, N is the total number of lattice points, ε surf is the surface energy per unit cell on the planar (001) surface and ε is the microscopic ledge energy associated with nearest neighbor (nn) interactions.The summation with respect to (n, m) is over all sites on the square lattice.The RSOS condition, meaning that the height difference between nearest neighbor sites is restricted to {0, ±1} , is required implicitly.In this equation, �µ is the driving force for crystal growth.To exclude diffusion effects, the ambient phase is assumed to be uniform.In the nanometer length scale near equilibrium, this assumption can be realized.In the case that the ambient phase is an ideal solution, �µ = k B T ln C/C eq 53 , where k B is the Boltzmann constant, T is temperature, C is the concentration of the solute and C eq is the concentration of the solute at saturation.If the ambient phase is an ideal gas, �µ = k B T ln P/P eq 54 , where P is the gas pressure and P eq is the gas pressure at equilibrium.ε and �µ in the present model correspond to φ and �µ in Van Veenendaal et al. 's work 22 .
Since the RSOS model is a coarse grained model used for the purpose of first principles quantum mechanical calculations, ε surf and ε relate to the surface free energy in the atomic model include the entropy for lattice vibra- tions and distortions 47 .Thus, these variables are affected by temperature.However, the present work assumes constant values for ε surf and ε in all calculations.

Monte Carlo calculations
In this work, the surface configuration was updated using the Metropolis algorithm and the energy difference, E , was calculated based on Eq. ( 4).The first 2 × 10 8 Monte Carlo steps per site (MCS/site) were ignored and each quantity was averaged over the subsequent 2 × 10 8 MCS/site.The surface slope p and the surface growth velocity V were calculated as macroscopic variables, such as where N step is the number of steps, which was fixed during the simulation, a = 1 is a lattice constant and τ is set to 2 × 10 8 MCS/site.
When considering an inclined surface, the squared surface width was calculated as where W is a surface width normal to the inclined surface, g is a geometrical factor defined as 1 + p 2 x + p 2 y with p x = ∂�h�/∂x and p y = ∂�h�/∂y 55 , x and ỹ are the [110] and [ 110] directions, respectively, and �•� ỹ and �•� x are the averages over the ỹ and x directions.
Periodic boundary conditions were adopted in the vertical ( [ 110] ) direction.In the horizontal ([110]) direc- tion, periodic boundary conditions were adopted while also adding the number of steps, N step .
Crystal growth proceeds by the attachment/detachment of specific units.As such, the number of units in the crystal does not have to be conserved during the process, making this a non-conserved system.The present work also did not include unit exchange on the surface, meaning that surface diffusion was neglected.At equilibrium, the unit attachment rate will equal the detachment rate.The attachment rate automatically increases whereas the detachment rate decreases as �µ increases.

Monte Carlo results
Figure 1 presents the Monte Carlo results for the surface growth velocity, V, and the scaled surface width, W, with regard to the (001) surface.Figure 1a indicates that the surface grows exponentially with respect to �µ during the 2D nucleation process for �µ/ε < 2.0 .In contrast, for 2.0 ≤ �µ , the surface grows linearly via an adhesive growth process.Because the temperature value of k B T/ε = 0.4 is far less than the thermal roughening temperature of k B T R /ε = 1.578 48,56 , the (001) surface is atomically and thermodynamically smooth at equilib- rium.It should be noted that the 2D critical nucleus sizes on the (001) surface were determined to be 2a and a for �µ/ε = 1 and 2, respectively, assuming that each nucleus was a square.In addition, the 2D critical nucleus sizes at the edges of the straight (01) steps were less than a for 1 < �µ/ε .In these processes, an atom (that is, the growth unit) attached at the edges of the steps associated with an island will increase the island's size on average in the case of 1 < �µ/ε .The attachment of an adatom to the (001) surface, which is also regarded as an island, will increase its size on average for 2 < �µ/ε.
Figure 1b,c provide the scaled surface width data.Near equilibrium and for �µ/ε < 0.55 , W = 0 and the surface is atomically and thermodynamically smooth.However, in the region defined by 0.55 ≤ �µ/ε < 2.0 , the surface width, W, increases as the system size, L, increases, meaning that the surface is thermodynamically rough.
To our surprise, in the present work the roughened surface was found to have a KPZ roughness exponent.In the non-equilibrium steady state, the roughness exponent α determines the universality class for a 2D growing surface (Eq.( 2)).In Fig. 1c the surface widths are scaled by L α and exhibit good agreement with the KPZ rough- ness exponent ( α ≈ 0.385 ).Herein, the KPZ roughening point is designated as �µ (001) KPZ = 0.55ε , while the end of the KPZ rough surface is designated as �µ (001) KtoB = 2.0ε.In Fig. 2a-c, images acquired at 4 × 10 8 MCS/site are shown for �µ/ε = 0.8 , 1.4 and 2.6, respectively.From Fig. 2a,b, it is evident that poly-nucleated multilayer islands appeared on the (001) surface and the perimeter of each island is also apparent.At �µ/ε = 0.8 (Fig. 2a), island-on-island structures can be seen, otherwise referred to as multilayer islands 57 , having distorted square morphologies.The side view presented in Fig. 2b also clearly indicates an island-on-island structure associated with �µ/ε = 1.4 .These islands are known to coalesce to complete the growth layer [24][25][26] , and so a self-affine surface exhibiting KPZ roughness was formed based on the island-on-island structure.
It should be noted that this KPZ rough surface represents a type of faceted rough surface that is atomically smooth but thermodynamically rough.The authors previously proposed the new concept of a so-called faceted rough surface and provided numerical evidence related to kinetic roughening 1 based on evaluating the surface roughness of surface systems between the atomic and mesoscopic length scales.A faceted rough surface is defined as being atomically smooth but thermodynamically rough, even though such surfaces grow via a 2D poly-nucleation process.Our prior work 1 using the Monte Carlo method demonstrated that an inclined surface meeting the criteria will be thermodynamically rough with a roughness exponent of α = 0.60 in the non-equilibrium steady state.This result provided evidence for the possibility of a large roughness exponent.�µ (001) KPZ is also a crossover point between the single and poly-nucleation growth processes.For small �µ values, V will be proportional to the single nucleation rate per area, I n ∝ exp(−G * /k B T) , according to the relationship V = L 2 I n .Here, G * is the free energy change for the formation of a critical nucleus.Based on the thermodynamics of a 2D island (see the "Analysis of 2D single nucleation" subsection in the Methods section), G * /k B T can be expressed as G * /k B T = γ 2 s,total /(4Sk B T�µ) ≡ g * /�µ (Eqs. ( 12 and ( 13)), where γ s,total is the total step free energy at the perimeter of the 2D equilibrium crystal shape with the Lagrange multiplier being 1 To better indicate the shapes of the steps on the crystal surfaces, the surface height is represented by 10 degrees of brightness, with a brighter color corresponding to a greater height.www.nature.com/scientificreports/and S is the corresponding area 58 .In this work, γ s,total and S were calculated based on the 2D Ising model using the imaginary path-weight random walk method 59,60 .Then, γ s,total a/ε = 6.441 and S/a 2 = 3.001 were determined at k B T/ε = 0.4 .On this basis, a value of g * /ε = 8.640 was obtained.Figure 3 plots ln(V ) as a function of ε/�µ and ln(V /v 2/3 s ) as a function of ε/�µ based on the Monte Carlo results in Fig. 1a, where v s is the step velocity (Eq.( 15) in the Methods section).For �µ < �µ (001) KPZ , the slopes of the lines obtained using different system sizes are in good agreement and are close to the value of 8.640 calculated using the Ising model.For larger �µ , the logarithm of V /v 2/3 s was plotted against 1/�µ (Fig. 3b) because for 2D poly-nucleation 25,26 .
The Monte Carlo results obtained for 0.55 < �µ/ε < 1.2 form a suitably straight line for L = 320 √ 2 , 240 √ 2 and 160 √ 2 .The least squares fit to these values gave a slope of −3.67 , the absolute value of which is larger than the expected value of g * /3 = 2.88 based on 2D poly-nucleation theory 25,26 Eq. (7).In contrast, in the case of 1.2 < �µ/ε , the Monte Carlo data deviate from a straight line.Here it is important to recall that �µ (001) kr = 1.15ε 39 for a relatively high growth velocity, V (Fig. 1a).
As is evident from Fig. 1b, for 2.0 ≤ �µ/ε , the surface approaches BKT roughness as �µ is increased.In addition, at �µ/ε = 2.6 (Fig. 2c), island-on-island structures are still seen but the size of the islands decreases and the side view of the surface indicates a greater number of fine irregularities.Furthermore, the surface velocity in this region increases linearly as �µ increases (Fig. 1a).Based on the results concerning this structure and the surface growth velocity data, the surface can be said to be kinetically, atomically and thermodynamically rough.
For large �µ , various types of kinetic roughening are known to occur in the crystal growth field [22][23][24][25][26]57 . Onekinetic roughening point with the classical criteria 22 , �µ c , is defined as πγ 2 s /(3k B T�µ c ) ∼ g * /(3�µ c ) ∼ 1 25 or g * /�µ c ∼ 1 26 .By substituting g * of the Monte Carlo result in the equation, we have �µ c /ε ∼ 2.88 or 8.64, which is a large value compared with the Monte Carlo result of �µ (001) KtoB = 2.0 at which adhesive growth starts for �µ (001) KtoB < �µ.Another criterion for �µ c is obtained from the zeros of the non-equilibrium step free energy 23 . Leving aside whether non-equilibrium step free energy is well defined, the zeros of the step free energy explain why the surface is atomically and thermodynamically rough, and why the surface steps are not clear.We estimated �µ c /ε for the Cuppen et al. results for the present case using Fig. 2 in Ref. 23 , giving a value of 1.3, where φ/k B T is 2.5, and It should be noted that the size of the critical nucleus is 1 at �µ/ε = 2.0 if we assume that the island shape is square (Eq.( 14) in the Methods section).As stated by Saito 25 , another criterion for �µ c is that the size of the critical nucleus becomes sufficiently small.We may therefore adopt "the size of the critical nucleus is one" as the criterion for �µ (001) KtoB .According to the 2D single nucleation theory (see the "Analysis of 2D single nucleation" subsection in the Methods section), the quantities at equilibrium are �µ and γs,total , where γs,total is the total step free energy around the 2D island.The step free energy at equilibrium should be calculated in the long step length limit to correctly take into consideration the entropy for the zig-zag structure of a step.However, when the size of the critical nucleus becomes 2 or 1, the entropy for the zig-zag structure is almost excluded because the step is too short to take a zig-zag structure.Since the lattice structure is simple cubic, the shape of the nucleus can be assumed to be a square with straight sides for a small nucleus, resulting from the finite size effect.Then, from Eq. ( 14) in the Methods section, the size of the critical nucleus c is given by c = 2ε/�µ .For 1.0 < �µ/ε < 2.0 , the critical size of a square on a (001) terrace is less than 2 but larger than 1.This means that on average a single atom on the (001) terrace detaches.For 2.0 < �µ/ε , the critical size of a square on a (001) terrace is less than 1.Hence, on average less than a single atom detaches from the (001) terrace.Therefore, we adopt "the size of the critical nucleus is one" as the criterion for �µ (001) KtoB .We will return to this point in the Discussion section.

AKPZ criteria
Here it is helpful to examine the extent of agreement between the Monte Carlo results and the KPZ equation.The crossover from BKT roughness to KPZ roughness on a surface can be discussed using the arguments proposed in Refs. 52,61.The relationship between the surface velocity and the fluctuation width was discussed by Wolf 61 using the renormalization group method based on the anisotropic KPZ (AKPZ) model.The values of x and ỹ are given by where p is the surface slope in the x direction for an inclined surface.Wolf 61 reported the criteria for the clas- sification of surface width as follows: For the limit p → 0 , the surface growth velocity on an inclined surface can be expressed as V ≈ v s p for �µ < �µ (001) KPZ (Fig. 4a,b).In the case of this �µ range, the surface grows based on a TSK process.In contrast, for �µ (001) KtoB , islands are frequently formed on the (001) terraces.Using a magnified version of x ỹ > 0 W ∝ L α (algebraic rough) x ỹ ≤ 0 α = 0, W 2 ∝ ln L. www.nature.com/scientificreports/Fig. 4a, we confirmed that V = V 0 + c p p 2 + O (p 3 ) for p → 0 , where c p is a positive coefficient at �µ/ε = 0.8 with L = 320 √ 2 .Note that this effect of the slope on V in the vicinity of p = 0 is revisited in the subsection titled "Kinetic shape changes of a crystallite." For �µ (001) KtoB , ∂ 2 V /∂p 2 > 0 and (∂V /∂p)/p > 0 near p = 0 , giving a coefficient value of x y > 0 .Hence, the surface should be algebraically rough and the Monte Carlo data also show KPZ roughen- ing on the surface.
In the case of �µ (001) KtoB < �µ , we have x y > 0 for the same reason as in the case of �µ (001) KtoB .Consequently, the surface should be algebraically rough but the Monte Carlo results suggest BKT roughening.The surface growth velocity scaled by V 1.061 (the surface growth velocity at p = 1.061 ) as a function of the slope value.Line: Eq. ( 10).Symbols: L = 240 √ 2a , 160 √ 2a and 80 √ 2a with a = 1 .Note that V is independent of system size.(b) A polar graph of surface velocity normal to the inclined surface, V n = V / √ g , where g = 1 + p 2 x + p 2 y .Taking angle θ = 0 as the 001 direction, V n are plotted from the origin to the normal direction of the surface between the � 11 1� ( −54.74 • ) and 111 (54.74 • ) directions.Dark shaded area: surface orientations less than −54.74 • and larger than 54.74 • .Light shaded area: terrace-step-kink (TSK) regions with 0.9 < |p| .L = 240 √ 2 . (c) Surface growth velocity as a function of �µ for several slopes.L = 160 √ 2a with a = 1 .k B T/ε = 0.4 .�µ cr /ε = 0.3.www.nature.com/scientificreports/Therefore, we conclude that the AKPZ criteria (Eq.( 9)) are partly consistent with the Monte Carlo results for a (001) surface.

Inclined surfaces
This section clarifies the difference in surface roughness between a TSK 41,42 model surface and a surface with terrace roughness for inclined surfaces.In the ideal TSK model, an inclined surface consists of a train of elementary steps having unit heights with no islands or negative islands (that is, clusters of adholes) on the terraces.Figure 5 summarizes the effect of �µ on the scaled surface width for several surface slopes and system sizes at a temperature of k B T/ε = 0.4 .Here, the upper subfigures show the squared surface widths scaled by ln L , while the lower subfigures present the surface widths scaled by L α .In Fig. 6, the calculated scaled surface widths for several �µ and system sizes at a temperature k B T/ε = 0.4 are plotted as functions of slope.The upper subfigures show the squared surface widths scaled by ln L , while the lower subfigures show the surface width scaled by L α with α = 0.385 (that is, the value of the KPZ roughness exponent).
The logarithmic divergence of gW 2 for an inclined surface at equilibrium is well established 4,5,18,19,39 .Here, g = 1 + p 2 x + p 2 y is the first fundamental quantity in the differential geometry 10,55 .For all slopes near equilibrium, the surface with �µ/ε ≤ 0.3 exhibits BKT roughening (see Figs. 5 and 6a).Hence, the present work confirms that the crossover point between BKT and KPZ roughening of an inclined surfaces is �µ cr /ε = 0.3.

Surfaces with almost ideal TSK structures
Hereafter, we consider an inclined surface for which �µ cr < �µ .As is evident from Fig. 4a, in the case of large surface slopes ( 0.9 < p ) near the (111) surface, the surface growth velocities are in good agreement with one another because the surfaces have almost the ideal TSK structure due to the RSOS restriction 39,43,62 .Figure 4c indicates the effect of �µ on V and demonstrates that V for p = 1.061 plateaus at 0.8 ≤ �µ/ε .Physically, surface growth occurs via the attachment and detachment of the atoms (here represented as cubes or growth units) at the step edges.
The lower subfigure in Fig. 5a demonstrates that, at p = 1.061 , a stepped surface without terrace islands becomes algebraically rough at �µ cr < �µ .Here, the roughness exponent is α = 0.33 and so is slightly smaller than the expected value for a KPZ roughened surface.From Fig. 6b-d, it is apparent from the surface widths for surfaces with 0.9 < p < 1.25 that the surfaces are algebraically rough.In addition, the roughness exponent appears to gradually decrease as p increases.
At this point, it is helpful to ascertain agreement with the AKPZ criteria in Eq. ( 9).Since ∂(V /V 1.601 )/∂p < 0 and ∂ 2 V /V 1.601 /∂p 2 < 0 , x ỹ > 0 based on Eq. ( 9).Hence, the surfaces should be algebraically rough, which is consistent with the Monte Carlo results for 0.9 < p .For the limit p → √ 2 , the present numerical results con- firm that , where v neg s is the step velocity for a negative step, meaning a step associated with a (111) terrace.These results demonstrate that the contribution of the nonlinear terms in the KPZ or AKPZ equation are reduced as p approaches √ 2. On this basis, we conclude that the effects of the slope and �µ on the surface width for 0.9 < p as obtained using the Monte Carlo method are consistent with the KPZ or AKPZ criteria.www.nature.com/scientificreports/Kinetic shape changes of a crystallite Figure 4c presents results for surfaces with p = 0.247 and 0.530 and shows that V increases steeply for �µ (001) kr < �µ as �µ increases, except for the surface for which p = 1.061 .This steep increase in V (other than the above exception) provides evidence that island formation on (001) terraces resulting from the 2D nucleation process causes a steep increase in V, as is also the case for p = 0.
When assessing crystal growth, it is interesting that the anisotropy with respect to V for �µ < �µ (001) KPZ is large compared with that in the corresponding Wulff figure 63 showing the polar graph of the surface tension.This can be seen from Fig. 4b.The significant anisotropy of V indicates that the crystallite grows such that it has a wider (001) facet compared with the equilibrium shape for �µ < �µ (001) KPZ .This effect occurs because there are two kinds of steps around p ∼ 0.7 associated with two kinds of terraces: the (001) and (111) terraces.As noted in the previous section, steps with (111) terraces and (001) side surfaces can be considered as negative steps.In this scenario, steps with small p values will grow to the right (e.g.see Fig. 7), whereas negative steps with large p values will grow to the left.For p values of 0.5-0.8, a surface having a mixture of steps including negative steps will grow in both directions (see the Supplementary Movie 1 S1).In the case of �µ < �µ (001) KPZ , the surface velocity, V, can be written as where w (001) and w (111) are the statistical weights for the number of steps determined by the surface slope, p.Because v s is approximately equivalent to v neg s (see Eq. ( 15)) for �µ < �µ (001) KPZ , the slope dependence of V is almost symmetrical along with p = 1/ √ 2 (Fig. 4a).Figure 4a plots the line obtained from Eq. ( 10) with and this line is seen to be in good agreement around both limits p → 0 and p → √ 2 .It should also be noted that the values obtained using the Monte Carlo method for p ∼ 0.7 are higher than those produced by Eq. (10).
For �µ (001) KPZ < �µ , because ∂V /∂p = 0 at p = 0 , the planar shape on the (001) surface becomes unstable in conjunction with infinitesimally small fluctuations in the slope.Hence, a kinetic version of the faceting transition is expected to occur at �µ kr , V continues to exhibit a high degree of anisotropy.Hence, the (001) surface should not be planar but rather should exhibit some small degree of curvature.For �µ (001) kr < �µ , the anisotropy in V is drastically reduced in the vicinity of the (001) surface.In particular, in the region �µ (001) KtoB < �µ , the concentration of adatoms on the (001) terraces is so large that the elementary step as the edge of the (001) terrace cannot be well defined (see Fig. 7c,c'), similar to the behavior of a surface with T R < T. 111) , www.nature.com/scientificreports/Nonlinear effect For p < 0.4 , the roughness change is complex due to the interplay between multilayered islands and steps.Near p = 0 , the surfaces are algebraically rough and √ gW decreases as p increases, as shown in Fig. 6b-d.With increases in p, a transition to a BKT roughened surfaces occurs near p ∼ 0.05 that is dependent on �µ .Here, the crossover point for p is denoted as p (p<0.4) KtoB .For �µ/ε = 0.8 , 1.4 and 2.2, p (p<0.4) KtoB ≈ 0.017 , 0.051 and 0.073, respectively.In the vicinity of p = 0.3 , there is another transition to an algebraically roughened surface and this crossover point is denoted as p (p<0.4) BtoK .For �µ/ε = 0.8 , 1.4 and 2.2, p (p<0.4) BtoK = 0.11 , 0.28 and 0.30, respectively.Both p (p<0.4) KtoB and p (p<0.4) BtoK increase as �µ increases.Here, it is helpful to assess the level of consistency between Monte Carlo results and the AKPZ criteria.Using the data in Fig. 4a, the inflection point with respect to p was calculated based on the Monte Carlo data.This process gave 0.056, 0.17 and 0.14 for �µ/ε = 0.8 , 1.4 and 2.2, respectively and these inflection points are close to those for p (p<0.4) KtoB .That is, when p < p (p<0.4) KtoB , ∂ 2 V /∂p 2 > 0 .On the basis of the criteria given by Eq. ( 9), the sur- face should therefore be algebraically rough.In contrast, if p (p<0.4) KtoB < p < 0.4 , ∂ 2 V /∂p 2 < 0 and Eq. ( 9) suggests that the surface should exhibit BKT roughening.Figure 7a-c demonstrate the coalescence of terrace islands to steps and this process enhances the step fluctuations for small values of p ( N step = 8 ).Therefore, we conclude that the nonlinear effect associated with surface growth makes the inclined surface near p = 0 algebraically rough.
Physically, the inflection point can be explained by an effect in which the inclined steps hinder the formation and free growth of multilayered islands.Figures 7a′-c′ present images of the surfaces for p = 0.247 ( N step = 28 ) and it is apparent that these surfaces appear different from those in (a), (b) and (c).In Fig. 7 a′-c′, fewer multilayered islands appear than in Fig. 2 and so it is evident that the KPZ structure was changed to a BKT structure.
In the case of p (p<0.4) BtoK < p , the transition from BKT roughened to algebraically roughened cannot be explained by the KPZ criteria.The following subsection discusses the origin of this crossover point.

Competition between step growth velocity and nucleation rate on the terrace
From Fig. 5c, it is apparent that, for �µ/ε ∼ 1.5 at p = 0.247 , the surface width becomes BKT rough.A peak around �µ/ε ∼ 0.8 is also apparent.In our previous work 39 , the peak was thought to be related to a kinetic roughening.However, from the extensive calculations in this study, that conclusion was found to be incorrect.In Fig. 7a′, the surface structure at �µ/ε = 0.8 is shown and this surface appears to be a TSK-type stepped surface.Few adatoms or adholes are seen on terraces in this area.Rather, as shown in the case of p = 1.061 , the TSK like surface crossovers from BKT rough to a power law rough for �µ cr < �µ.
While, for the surface structure at �µ/ε = 1.4 , which is presented in Fig. 7b′, the elementary steps show numerous overhang structures and there are small numbers of islands or negative islands on the terraces.The poly-nucleated clusters on the terraces merge with growing steps, and then the step edges have generated overhang structures.Because islands having different heights or negative islands cannot merge completely with steps, the higher islands or lower negative islands act as obstacles to the growing steps.In this manner, fluctuations of the steps are reduced by the multi-height islands or negative islands.For �µ (001) kr < �µ , there is a non-negligible reduction in step fluctuations that produces a BKT roughened surface.As a result, the surface width decreases to form a peak as seen in Fig. 5c.We denote this crossover point for �µ as �µ

Discussion
In the case that of an inclined surface for which the temperature is higher than k B T/ε = 0.4 , the crossover point �µ cr becomes larger 39 .Hence, for k B T/ε = 0.63 and 1.7, �µ cr /ε becomes 0.5 and 1.2, respectively 39 .The population of adatoms on the (001) terraces will therefore be larger at higher temperatures and these adatoms will partially block the step fluctuations that generate surface BKT roughening.
In contrast, �µ KPZ decreases rapidly as temperature increases.Since the step free energy at equilibrium decreases as the temperature increases, 2D nuclei are formed more frequently, which leads to a crossover to poly nucleation easily.�µ (001) KPZ becomes almost zero at k B T/ε ∼ 0.9 in our preliminary simulations.Figures 4a,b reproduce qualitatively the characteristics of the angle dependence of the growth rate shown in Fig. 4 in the paper of Van Veenendaal et al. 22 , although the models are different from each other.At small �µ (Fig. 4 a in Ref. [22]), the angle dependence of the surface growth rate has a cusp singularity at p = 0 for φ/k B T = 2 and �µ/k B T = 1 ( k B T/ε = 0.5 and �µ/ε = 0.5 ).For large �µ (Fig. 4 c in Ref. [22]), where φ/k B T = 2 and �µ/k B T = 3 , the angle dependence of the surface growth rate is parabolic at p = 0 .These agreements between the present results and Van Veenendaal et al. 's results around p = 0 indicate that the growth shape change of a crystalline near the (001) facet is a universal phenomenon associated with KPZ roughening.
On the other hand, for �µ c for the (001) surface, our preliminary results were complex.If we adopt "the size of the critical nucleus is one" as the criterion for �µ KtoB does not depend on temperature if ε does not change.Referring to Fig. 2 in Ref. 23 , �µ c decreases almost linearly from the value for k B T/ε = 0.5 to zero at about the T c of the 2D nn Ising model.From our preliminary Monte Carlo results, there is a possibility that �µ (001) KtoB may be different from �µ c .For the surface growth velocity, ∂V /∂�µ seems to change discontinuously around �µ/ε = 2.0 for several temperatures up to k B T/ε = 1.4 , which is lower than T (001) R but higher than T c of the 2D nn Ising model.In addition, another characteristic point �µ where Ŝ is the area of the island and γs,total is the total step free energy at the perimeter of the island.In the case that the island is the critical nucleus, ∂G/∂ = 0 and we have Since G * /k B T ≡ g * /�µ , we can write When the shape is a square with side length a , we have Ŝ = ( a) 2 /a 2 and γs,total = 4( a)(ε/a) , which leads to S = 1 and γ s,total = 4ε .Then, we have

Estimation of v s
To obtain the explicit form for v s , this work used parameters that provided the best least squares fit to the Monte Carlo results in Fig. 4 (c) for a negative step velocity v neg s at p = 1.061 .Here, a negative step is defined as a step with a (111) terrace and (001) side surface, such that The v s value at p = 0.09 was confirmed to equal v neg s at p = 1.061 within a difference of 5%.The results obtained with p = 1.061 were employed to determine v s because there was a lack of nucleation on the (111) terraces for negative steps due to the RSOS restriction.

Figure 1 .
Figure 1.Surface growth velocity and scaled surface widths at the (001) surface ( p = 0 ) as functions of �µ .(a) Surface growth velocity with unit of a/τ , where a ( = 1 ) is the unit height and τ is the time interval for 1 MCS/site.Line: V = 0.0643�µ/ε − 0.0412 . (b) The squared surface width, W 2 = �[h(� x) − �h�] 2 � , scaled by the logarithm of the system size, L. (c) The surface width scaled by L α with the roughness exponent α ≈ 0.385 determined from the KPZ model.k B T/ε = 0.4.

Figure 2 .
Figure 2. Images of simulated p = 0 surfaces.Upper figures: overhead views.Lower figures: side views showing the height along the lower perimeters of the upper figures for�µ/ε = (a) 0.8, (b) 1.4 and (c) 2.6.k B T/ε = 0.4 .L = 320 × √ 2 .�µ cr /ε = 0.3 39 .To better indicate the shapes of the steps on the crystal surfaces, the surface height is represented by 10 degrees of brightness, with a brighter color corresponding to a greater height.