Driving force dependence of the height of a faceted macrostep in non-equilibrium steady-state crystal growth

In order to understand the dynamics of the self-organized macrosteps, the vicinal surface with faceted macrosteps is studied by the Monte Carlo method based on a microscopic lattice model, the restricted solid-on-solid model with point-contact-type step-step attraction (p-RSOS model). We focus on the dynamical effects caused by the change of the surface roughness or the change of the kink density which are masked by the effect of the surface and volume diffusion of the crystal atoms in the ambient phase. Contrast to the step-bunching in the diffusion-limited cases, the height of the faceted macrostep decreases as the driving force for the crystal growth increases.


Introduction
Making good quality crystals of SiC is an urgent problem for the production of future power devices with low power consumption. Self-organized macrosteps are known to hinder making good quality of crystal of SiC [1]. We studied the largeness of the macrosteps theoretically based on a lattice model, the restricted solid-on-solid model with point-contact-type step-step attraction (p-RSOS model) [2,3]. We focus on the surface dynamics with the change of kink density or surface roughness, which has long been overlooked in studies of surface and volume diffusion. In our lattice model, a faceted macrostep exists on a vicinal surface stably, even at equilibrium at low temperatures. The characteristics of the profile of the macrostep are classified by the connectivity of the surface free energy (or the surface tension) which was calculated by the density matrix renormalization group method [4]- [9]. Driving force dependence of the profile of vicinal surface with faceted macrosteps is analyzed by the Monte Carlo method [14,15]. Contrast to the step-bunching in the diffusion-limited cases [10]- [13], the height of the faceted macrostep decreases as the driving force for the crystal growth increases.

The model and statistical mechanical calculations 2.1. The p-RSOS model
The microscopic model considered in this study is the p-RSOS model ( Fig. 1) [2,14]- [20]. The total energy of the (001) surface can be written as (a) (b) Figure 1. Perspective view of the surface. This figure is taken from [16].
where 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. The RSOS condition, which is the height difference between the NN sites restricted to {0, ±1}, is required implicitly. Here, δ(a, b) is the Kronecker delta, and ϵ int is the microscopic step-step interaction energy. ϵ int contributes to the surface energy only at the meeting point of neighboring steps because the height difference can be ±2 only at the meeting points of neighboring steps. The summation with respect to (n, m) is taken over all sites on the square lattice. When ϵ int is negative, the step-step interaction becomes attractive (sticky steps).

Calculations of the surface tension: DMRG (PWFRG) calculations
The surface energy of a vicinal surface is calculated by the following equation [21]: where η = (η x , η y ) is the the Andreev field [22]. The Andreev field is a chemical potential with respect to a single step. From statistical mechanics, the grand partition function Z of the vicinal surface is calculated as follows: The summation with respect to {h(m, n)} is taken over all possible values of h(m, n). The thermodynamic grand potential calculated from the grand partition function Z is the Andreev free energyf (η) [22], which is obtained as follows: In low-dimensional cases, more precision is required than that provided by a mean field calculation of the partition function [23]. Hence, to obtain reliable results, we adopt the DMRG method [4]- [9], which was developed for one-dimensional quantum spin systems. We use the transfer matrix version of the DMRG method, which is known as the product wave-function renormalization group (PWFRG, tensor network) method [7,8,9].
The surface free energy per projected x-y area f (p) is obtained fromf (η) using the following equation: where p = (p x , p y ) is the surface gradient of the vicinal surface. The surface tension is the surface free energy per unit normal area. Hence, the surface tension γ surf (p) was calculated from the surface free energy f (p) as [24] γ surf (p) = f (p)

Faceting diagram obtained by the anomalous surface tension
In this manner, we obtained the surface tension of the p-RSOS model numerically. The Fig. 2 shows the Wulff figures. The blue lines show the polar graph of surface tensions. The pink lines show the Andreev free energy which is similar to the equilibrium crystal shapes. These light blue lines show the surface tension for metastable states.
As seen from Fig. 2, there are two transition temperatures. At temperatures lower than T f,1 , the surface tension becomes discontinuous around (111) facet. At temperatures lower than T f,2 , the surface tension becomes discontinuous around (001) facet. So, we call the zone at the temperatures less than T f,2 the step-faceting zone; the zone at the temperatures T f,2 < T < T f,1 , the step droplet zone; and the zone at the temperatures T f,1 < T , we call the Gruber-Mullins-Pokrovsky-Talapov (GMPT) zone.
According to the connectivity of the surface tension, we can obtain the faceting diagram (Fig.  3). In the step-faceting zone, only the (001) surface and the (111) surface are thermodynamically stable. Hence, the vicinal surface is formed by the giant steps of the smooth (001) and the (111) surfaces. Whereas, in the step droplet zone, the (111) surface and the vicinal surface with the slope of p 1 are thermodynamically stable. The vicinal surface is formed by macrosteps with the (111) side surface and the stepped surface with the slope of p 1 .
The characteristics of the profile of the faceted macrostep is classified according to the connectivity of the surface tension [16].

Monte Carlo method
To study the non-equilibrium steady-state with macrosteps, a vicinal surface with the following Hamiltonian and 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), and ∆mu is the driving force of the crystal growth (chemical potential difference between the bulk crystal and the ambient phase). When ∆µ > 0, the crystal grows, whereas when ∆µ < 0, the crystal recedes (evaporates, dissociates, or melts). We required periodic boundary condition in the direction of the mean step-running. In the direction normal to the mean step-running, the surface heights of the lowest side line are connected to the the surface heights at the topmost line by adding the number of elementary steps. As the initial configurations of the surface, the parallel train of steps and the single macrostep are prepared.
In the case of the taking average, the first 2 × 10 8 Monte Carlo steps per site (MCS/site) were discarded, and then averaged over following 2 × 10 8 MCS/site. The snapshots of the simulated surfaces are shown in Fig. 4.

3.2.
Step droplet zone In Fig. 5, the driving force dependence of the average height of the macrosteps and the velocity of the surface [14] are shown.
The average height of the macrosteps are obtained as follows: wherex is selected as the ⟨110⟩ direction (normal to the mean step-running direction),ỹ is the ⟨110⟩ 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. To estimate the growth rate of the surface V , the average surface heighth(t), was calculated, The growth rate of the surface V , is defined as: where t 0 and t max are 2 × 10 8 MCS/site and 4 × 10 8 MCS/site, respectively. By analyzing the obtained values, we found a crossover point ∆µ R /ϵ = 0.005±0.002. As |∆µ| increases, the configuration of the vicinal surface crossover to the two surface co-existent one to the kinetically roughened one. In the two-surface co-existent case, elementary steps attach to or detach from the faceted side surface. To understand the ∆µ dependence of the vicinal surface morphology for |∆µ| < ∆µ R , we considered a step-attachment-detachment model for the time evolution of ⟨n⟩ [14]: 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 (the case of ∆µ R < |∆µ|), whereas when n + > n − (the case of ∆µ f < |∆µ| < ∆µ co , §3.3), ⟨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. Then, in the case of |∆µ| < ∆µ R , after some calculations [14], V is expressed as follows: p 1 (∆µ) = 0.332 + 15.6|∆µ|/ϵ + 4.43 × 10 3 (|∆µ|/ϵ) 2 , where p 1 (∆µ) is the slope of the "terrace surface", v 1 (∆µ) is the mean step velocity of the elementary steps obtained by this Monte Carlo calculations. We show the calculated lines in the Fig. 5 (b) by using Eq. (11) together with Eqs. (12) and (13)  On the other hand, ⟨n⟩ is expressed by z, N m and N step as follows: where N m is the number of the macrosteps, andp = aN step /L is 1/ √ 2. The lines of ⟨n⟩ are shown by the light blue lines in Fig. 5 (a). These lines also reproduce the Monte Carlo results well.

3.3.
Step-faceting zone  When the temperature is low enough in the step-faceting zone, the driving force dependences of the height of the macrostep and the growth velocity of the surface change from the case in the previous subsection. The vicinal surface at equilibrium consists of the (001) surfaces and the (111) surfaces because only these surfaces are thermodynamically stable. The most notable characteristic in the dynamics in the step-faceting zone is the freezing of the motion of the surface in the area near equilibrium. Since the both of the (001) and the (111) surfaces are smooth, there are very few kinks on the surface. Hence, two-dimensional nucleation process is needed to grow the surfaces. However, the mean waiting time of the formation of the critical nucleus becomes exponentially longer as the driving force |∆µ| → 0. The waiting time easily exceed the simulation time near equilibrium. In this manner, the dynamics near equilibrium strongly affected by the finiteness of the system size and time.
Then, by analyzing the obtained values, we found several characteristic values for the driving force in addition to ∆µ R . We show them in Table 1. For |∆µ| < ∆µ f , no nucleation occurs during the simulation time M CS max on the average. Hence, the surface configuration almost freezes and it depends on the initial configuration ( Fig. 6 (a) and (b)).
For ∆µ f < |∆µ| < ∆µ co , the surface grows intermittently in the manner of the single nucleation on the line of the lowest edge of the macrostep.
For ∆µ co < |∆µ| < ∆µ R , the surface grows with the successive detachment of elementary steps from the lowest edge of the macrostep.
To understand the ∆µ dependence of the vicinal surface morphology for ∆µ co < |∆µ| < ∆µ R , we again considered a step-attachment-detachment model for the time evolution of ⟨n⟩. The Eqs. (11) and (14) are also applicable to this case. However, the equations of p 1 (∆µ) and v 1 (∆µ) become complex due to the two-dimensional multi-nucleation on the edge of the macrostep. After some calculations [15], p 1 (∆µ) and v 1 (∆µ) are expressed as follows: Here, ∆µ y (L) is introduced [15] so that the growth velocity of the surface is well explained by the classical nucleation theory qualitatively and quantitatively. Physically, ∆µ y (L) indecates a yield point of the self-detachment of elementary steps from a macrostep. The surface velocity (Eq.(11)) calculated by the Eqs. (15) and (16) is shown in Fig. 7 with black solid lines. The mean height of the macrostep (Eq. 14) are also shown in Fig. 6

Kinetically roughened surface
For ∆µ R < |∆µ|, the surface roughens kinetically. Both of the surface velocity and the mean height of the macrostep obey power law behaviors.

Discussions
It is interesting that the figure for V (Fig. 7) is analogous to Fig. 4(b) in Ref. [25], 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 step attachment-detachment motion on the terrace is analogous to the motion of particles with plastic depinning. However, a common mathematical framework has not yet been obtained.

Conclusion
The height of the faceted macrostep decreases as the absolute value of the driving force of the crystal growth |∆µ| increases.