Kinetic Monte Carlo model for homoepitaxial growth of Ga 2 O 3

We developed a kinetic Monte Carlo (KMC) model for the homoepitaxy of β -Ga 2 O 3 . It comprises adsorption, diffusion, and desorption and reﬂects the structure of β -Ga 2 O 3 with its two kinds of atoms: Ga and O. The knowledge gained from metal organic vapour phase experiments (MOVPE) experiments combined with AFM and TEM characterisation was used for the setup of rules and activation energies for the various surface processes. We performed a set of runs for the growth on ﬂat and vicinal (100) surfaces. The nucleation on the ﬂat surface requires a minimum ratio of the impingement rate of O 2 and Ga. The behavior at different substrate temperatures was similar in experiment and in simulation. At high temperatures, we observe the formation of large islands whereas at low temperatures small islands are formed. The growth rate is increasing with decreasing temperature. On a vicinal surface (6 ◦ ) different growth modes have been observed when using different desorption energies. Low desorption energy (high desorption rate) leads to step bunching, intermediate to step growth, and high energy (low desorption rate) to nucleation on terraces with a ﬁnal conﬁguration similar to step bunching.


I. INTRODUCTION
β-Ga 2 O 3 is a material of great interest for different applications.With a bandgap about 4.8 eV and an estimated electrical break down field of 8 MV cm −1 it is an excellent candidate for power electronics and optoelectronic applications [1,2].Bulk crystals of β-Ga 2 O 3 can be grown from melt with high perfection and thus substrates of high quality are available for homoepitaxial growth of device-grade layers.
Despite the great progress in depositing epitaxial layers of high quality, the understanding of the kinetic processes is still very limited.One attempt to understand the growth dynamics is by means of kinetic Monte Carlo (KMC) methods.The aim of this paper is to set up an appropriate statistical model for KMC computation and investigate the basic dynamics of the model.The guidelines are the results of our various MOVPE experiments and their analysis by means of AFM and TEM.
Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
Using the model as described in the next section we are able to recover experimental growth in a semiquantitative manner.Certain parameters are varied to see the influence on the growth dynamics.

II. MODEL
We use an on-lattice KMC model, which reflects the structure of Ga 2 O 3 .The unit cells of Ga 2 O 3 are forming a grid of cells labeled as shown in Fig. 1.For sake of simplicity, we omit the term "unit" in front of cell from now on.Each cell has a unique neighbor in the lateral and diagonal directions.Please note that we use periodic boundary conditions in x and y directions.In x direction, the boundary might include a step of one cell.This allows to compute a stepped surface.In particular, we have two steps because there are two stable planes per unit cell in b-c direction (see Fig. 3).The neighbors at the boundary for the flat and vicinal surface can be seen in Fig. 1.
Every cell contains eight positions for gallium atoms and 12 positions for oxygen atoms.The labeling of the atoms can be seen in Fig. 2, where two cells are shown from the 010 (b axis) and 001 (c axis) directions.The Ga atoms are labeled from a1 to a8 and the O atoms from a9 to a20.Atoms in the neighboring cells are labeled additionally with an e, w, n, s, t, b for east (+x direction) west (−x direction), north (+y direction), south (−y direction), top (+z direction), and bottom (−z direction), respectively.
An atomic position i in cell m can be occupied or unoccupied.In both cases, there are two possibilities: s i = 2: atom is in bulk and cannot move, s i = 1: atom can move.s i = −1: position is empty and cannot be occupied (gas phase), s i = 0: position can be occupied.For every position we compute the number of bonds to bottom (N bond bottom ), in layer (N bond layer ), and to top (N bond top ).We distinguish three types of bonds with respect to their energy.Bonds in a * direction (e.g., a8-a19) have the FIG. 1. Numbering of cells in the x-y plane.The numbering is given for the first layer (l = 1).For other layers, one has to add l * Nx * Ny to the given numbers.Periodic boundary conditions are applied in x and y directions.The neighboring cells are numbered in blue color.In x direction, we allow to have a step of a unit cell.In this case, the neighboring cells at the left and right boundary are given by the outer expression in the black box.
energy E I , bonds in the b-c plane (e.g., a8-a9, a8-a18) have the energy E II , and the others the energy E III .The total binding energy for an atom at a position will be computed and stored.This is the basis for calculating rates, which will be discussed later in detail.At this point, we like to present only the general outline for setting and choosing events.In the beginning, a list of all possible events will be created by scanning over all cells.The scheme is quite flexible and different kinds of events can be defined.For every event n, the transition rate R n and the initial positions and final positions are stored.The total rate at time t is given by R total (t ) = n max (t ) 1 R n (t ), where n max (t ) is the total number of events at time t.The events include adsorption, diffusion, and desorption.The rate for adsorption is given by an impingement rate (R ads ) and a local sticking coefficient S: R n (t ) = S (t )R ads .For diffusion and desorption the rate is given by R n (t ) = ωe −E n /k B T , where E n , k B , and T are the activation energy, the Boltzmann constant, and the temperature, respectively.The frequency ω is set to the Debye frequency (see below).In every iteration step (time t), one event n 1 out of the complete list is chosen by a random number 0 < r 1 < 1: R total (t ) .r 1 and r 2 are independent random numbers, generated by the package proposed by L'Ecuyer et al. [14].After choosing and performing the event, the cells of the event itself and all neighboring cells are updated concerning the possible events.Technically, this is done by deleting all entries for these cells from the list of events, compute the possible events for all these cells, and store them.The core parts of the code like event choosing and random number generation are based on the code used in Refs.[15,16].

III. KINETICS
As already mentioned, in this paper we will consider the MOVPE process.Gallium is provided via the precursor triethylgallium (TEGa) and the growth temperature is between 750 • C and 850 • C. Oxygen is provided by a flux of O 2 at pressures of some millibar.Only the growth on the (100) plane will be considered, because for no other plane a step flow mode was observed so far.
STEM HAADF images of the surfaces reveal that the terraces are always a B1 or B2 plane (see Fig. 3).These planes have a low surface energy and do not exhibit relevant relaxation according to density functional theory (DFT) calculations [17].Figure 3 shows a sketch of two steps.The plane at the step edge is (201), which has also a low surface energy according to DFT calculations.However, DFT also reveals that the (201) plane exhibits significant relaxation of the top most layer [17].This relaxation is not considered in Fig. 3.

A. General considerations for energies and rates
The main ingredients for KMC are the activation energies for diffusion processes.Currently, no data are available from DFT calculations but a guess for the magnitude can be derived from the investigations of Schewski et al.They analysed the surface structure by AFM and STEM and gave a phenomenological description of the observed twin formation [10].Moreover, they adopt a model by Bales [18] to obtain a quantitative description by a mean field approach.The best fit of the experimental and numerical results for the nucleation density versus terrace width was obtained for a diffusion coefficient of D = 7 × 10 −13 m 2 /s.This diffusion coefficient is based on using one (artificial) particle per cell traveling in c direction.At this point, we assume that the diffusion time is dominated by the diffusion of a gallium and an oxygen atom on the B plane.Therefore we use this diffusion coefficient as a rough estimation for the energy barriers for the diffusion of a single gallium and oxygen atom in c direction on a B plane.
Assuming the Debye frequency ω D to be the characteristic one for the diffusion process we can write where ( We will take this value as a rough estimation for the activation energy of moving in c direction.
FIG. 3. Surface planes as observed in the STEM analysis of grown layers by step flow mode [17].Only the relevant atom labeling is given.Please note that in contrary to Fig. 2 the positive c direction is going to the right.
The bond energy between a Ga and O atom on the B1 or B2 plane is about 0.33 eV resulting from the surface energy of 0.6 J/m 2 computed by DFT [17].However, this value was calculated from a bulk and a slab with surface and thus the bond energy of 0.33 eV is for the situation in bulk.A single Ga atom on the B1 or B2 plane has no further bond except the one to the bottom oxygen atom and therefore, the binding energy for this bond might be much larger than in bulk.We will use this considerations later for defining the activation energies of clusters.
In the experiments a growth rate of about 3 nm/min was observed for a step flow growth.This means that about 2.5 unit cells will be added in a minute in a * direction (a a * = 11.88Å).Or inversely, one layer (half unit cell) is grown in about 12 s.

Adsorption of gallium
In recent experiments triethylgallium (TEGa) was used as the precursor for Ga.A theoretical study based on quantum mechanical calculations predicts the decomposition of triethylgallium to the final product GaH [20].This was also observed in an experimental study of trimethylgallium (TMGa) in reaction with H 2 [21].The additional hydrogen is required to saturate the leaving methane.In the case of triethylgallium, no additional hydrogen is required because ethylene is leaving and not ethane.When using TMGa layer by layer growth of Ga 2 O 3 could only be observed when using water instead of oxygen in the atmosphere [5].With TEGa all experiments have been performed with oxygen.
As a conclusion of the various results we expect GaH to adsorb on the surface.The hydrogen atom can also act as a catalysator for the adsorption of an oxygen molecule.As a result an OH group might be formed.In order not to complicate things too much, we neglect the hydrogen in the current version of KMC and consider only Ga and O atoms.
We set different sticking coefficients for the tetrahedrally and octahedrally coordinated Ga atoms.For the first type (a1-a4), we set N bond bottom is the number of bonds to oxygen atoms at bottom and might be one or two for atoms a1-a4.For the second type (a5-a8), we set In this case, we have only one bond to a oxygen atom at bottom but might have also bondings to oxygen atoms in the same layer [(100) plane] indicated by N bond layer .The more oxygen atoms are present for bonds to the arriving gallium atom the higher the probability to adsorb.For positions a6 and a8, we apply a special rule.In case that for a6 not only the oxygen atom a13, where a6 is going to stick to, is present but also its two neighbors in north and south direction we set S = 0.5.The same holds for a8 and the oxygen atoms a15.The higher sticking coefficient ensures a sufficient adsorption rate at the edges with respect to those on the B1 and B2 plane.The higher probability can be also expected from the fact that the oxygen atoms for the adsorption of a6 and a8 has less bonds than the ones for the adsorption of a5 and a7.Details have to be clarified by DFT calculations.

Adsorption of oxygen
We assume that the O 2 molecule will adsorb on a Ga atom like it is shown in Fig. 4. Firstly, a peroxidelike configuration will be formed [see configuration (b)].Experimental results and quantum-chemical calculations showed the possibility of a spontaneous reaction of Ga and O 2 to form a cyclic peroxide GaO 2 [22].The distance between the oxygen atoms in the peroxide (1.35 Å) is only slightly larger than in the O 2 molecule (1.21 Å).Probably, for the second step [see configuration (c)] there might be an activation necessary, which might be supported by hydrogen.The oxygen atoms are moving to their FIG.4. Adsorption of an oxygen molecule on a Ga atom (here a5).This is just a schematic view to describe the rules for KMC.Probably in reality hydrogen is involved at least in process (b) and one of the oxygen atoms will form an OH bond.
positions in the layer of the Ga atom.This process can happen on the Ga atoms a5, a6, a7, and a8: For the adsorption of O 2 on Ga atoms a5, a7, a6, and a8 there are different possibilities for the oxygen atoms to go in step (b) in Fig. 4. The possibilities depend on the already existing bonds of the Ga atom to oxygen atoms in the same layer (e.g., for a5: a12, a12s, a19, and a19s).If no oxygen atom is bound to the Ga atom, there are in total six possibilities and for each we have S = 1/6S * .If one oxygen atom is bound then there are three possibilities for the adsorbing oxygen atoms, i.e., S = 1/3S * for each possibility.If there are already two oxygen atoms bound to the Ga atom then there is only one possibility with S = S * .On a2 and a4, the peroxide structure can remain because the structure in bulk similar but with a larger distance of the oxygen atoms (2.9 Å).Therefore adsorption should be easy and we set S = 1.0 for the adsorption of a15 and a15s on a2 as well as for a13 and a13n on a4.On the Ga atoms a1 and a3 we do not allow adsorption of O 2 .

C. Diffusion
We allow Ga atoms and O atoms to move.Ga atoms can move in the same layer (a1-a8, a4-a5, a3-a6, a2-a7), a layer down, or a layer up.O atoms can move in the same layer (a16-a15b, a9-a18, a12-a19, a13-a14, a11-a20, a10-a17), or a layer down.In addition, on the B planes also clusters of Ga and O can move.This is based on the assumption that lateral bonds weaken the bond to the bottom layer and thus the activation energy to move a small cluster is not much higher than for the movement of a single Ga atom.
The activation energies for the process are related to the bond energies between neighbor Ga and O atoms.Bond energies are set as follows: For Ga atoms a1-a4, all bonds have energy E III .For Ga atoms a5-a8, bonds to bottom and top O atoms have energy E I and for bonds in the same layer E II .

Diffusion of gallium atoms
From Fig. 2(b), one can see that the diffusion will be anisotropic in b and c directions.In ±b direction, a Ga atom can easily diffuse from one position to the same position in the neighboring cell.In ±c direction, a Ga atom has a longer distance to move to the same position in the neighboring cell.
Firstly, we define the diffusion in the same layer.All Ga atoms can move in the ±b direction (north ↔ south), as, e.g., a1→ a1n and a1→ a1s.The activation energy for diffusion A is given by E bond is the sum of bond energies to nearest neighbors of oxygen and E Ga edge is an extra energy if the Ga atoms are in vicinity of an edge(for a2 or a7 an edge is defined if a1 on left or right side is occupied.For a4 or a5, it is defined if a3 on left or right side is occupied).
Diffusion in ±c direction (east ↔ west) is more complex.In the following, we list the possible movements and the related extra penalty E extra for these diffusions with respect to Eqs. ( 5) and ( 6

Diffusion of oxygen atoms
Oxygen atoms in the layer of Ga atoms a5-a8 can diffuse around these atoms with an activation energy given by E bond .
In particular, we have the following movements:

Diffusion of clusters
We consider the following clusters: Ga 2 O 2 , GaO 4 , GaO 3 , GaO 2 , and GaO.Explicitly, the possible clusters are shown in Fig. 5 with their activation energy E A for diffusion.The energies are given as fraction of the energy between a single Ga or O atom and the B1 (or B2) plane, which is E I ≈ 1.6 eV.In the bulk Ga and O are bond by about 0.33 eV.For the cluster, we assume intermediate values.

D. Desorption
It is expected that the suboxide Ga 2 O is desorbing from the surface, because quadrupole mass spectrometer (QMS) measurements during MBE growth revealed the formation of the volatile gallium suboxide Ga 2 O. [23].Voigt et al. observed a Ga 2 O flux and also analysed the growth rate as a function of the temperature.Based on these results they deduce that Ga 2 O is desorbing.Thus we take this desorption into account and set the activation energies depending on the configuration as follows: The expressions for E 1 des are given in Table I.

IV. CALCULATIONS
We performed runs with different impingement rates of oxygen and different desorption rates of Ga 2 O.In the first set, we change the relation between impingement of Ga and that of O 2 .This will influence the diffusion of Ga atoms on the terraces and thus, it is important for the mean diffusion length.Also the second variation, i.e., that of desorption rate, will influence the mean diffusion length on the surface.Because for desorption two gallium and one oxygen atom are needed, the relationship between desorption rates and surface morphology is not obvious on a first glance.The various parameters for the runs are listed in Table II A. Flat surface We start with a flat B2 plane and the crucial point is to form a critical nucleus on this plane.The growth of the nucleus is in competition with desorption and diffusion.Critical means that an agglomeration of atoms is more likely to grow than to shrink.Please keep in mind that desorbing molecules are gallium suboxides containing three atoms and also diffusion might include more than one atom but up to five in the case of clusters.Firstly, we discuss the results for the runs A2a-A2d.We define the ratio between oxygen and Ga flux as R ads O 2 /Ga = R ads O 2 /R ads Ga .In the case of the low ratio of R ads O 2 /Ga = 10 (A2a) no growth is observed (runs were performed up to 12 s in a domain of 20 × 32 cells).Ga 2 O is desorbing before a critical nucleus could be created.The desorption can be significantly lowered by increasing E Ga des from 1.0 eV (A2a) to 1.5 eV (A3a).Thus nuclei can be formed after short time and we observe a growth (see Fig. 6).However, there is evidence from experiments that desorption takes place.Nucleation can be also enforced by in increasing the O 2 to Ga ratio.We discuss now the runs A2b-A2d in Fig. 6.The larger R ads O 2 /Ga the faster nucleation and growth.The faster nucleation the higher is the density of nuclei and consequently the smaller the 2D structures are.In cases A2c and A2d, we have already several nuclei in the runs with the domain size 20 × 32 so that we do not see a difference in the growth rate when using a domain of size 64 × 64.All runs have in common that the islands grow in height only up to the next B plane.This means that the entire layer is filled before a new critical nucleus is created.This behaviour can be clearly observed in Fig. 6 by the short plateaus at 4, 8, and 12 Ga atoms per unit cell area.This behavior of filling up a layer was also observed in experiments.
In the case of R ads O 2 /Ga = 20 (A2b), it takes some seconds before a critical nucleus is created.The structures are much larger than in the cases discussed above.In order to minimise the self-interaction of the occurring 2D structure with itself, we performed runs with a domain of size 40 × 100.Even with this domain size we observe only one nucleus per B plane.In particular, we performed four runs with different initialisations of the random generator.The results can be seen in Fig. 6.The variation in nucleation time and transient behavior is much larger than in the other cases, because we have only one nucleus and thus no statistical average.
After nucleation an elongated island is forming (see right top figure in Fig. 8.The top of this island is a B1 plane.Because the area of the 2D island becomes much larger than a critical nucleus, such a nucleus can be created before the B1 plane is completed.In this case no plateau occurs in the growth rate at 4 Ga atoms per unit cell area (see Fig. 6).Please note that due to the large size of the island compared to the domain size we have finite size effects.
For comparison we also show the growth rate for the step growth (case A2b step).As we already have a step no initialization is needed and we have almost a straight line in Fig. 6, i.e. a constant growth rate.
Summing up the results for a flat surface: the higher the O 2 impingement rate the higher the growth rate is.If the O 2 impingement rate is too low for the given impingement rate of Ga no critical nucleus is formed.Forming of a critical nucleus can be also induced by reducing the desorption rate by choosing a larger value for E Ga des .For run A3a in Fig. 6 all parameters are the same as for A2a, where no critical nucleus could be achieved within 12 s, except that E Ga des is 1.5 eV instead of 1.0 eV.The plateau is not as well expressed as for the other cases but it is still there.
The behavior in the results of the KMC calculations corresponds with that observed in experiments on a flat (100) plane (details are given in Ref. [24]).Thanks to the excellent cleavage such planes can be easily prepared with high perfection.Experiments have been performed for three different substrate temperatures.Every run has been stopped after a specific time after starting exposure of Ga (10,15,20,25,30, 35, and 40 s).The flux rates were 6.07 × 10 −3 and 0.02 mol min −1 for TEGa and for O 2 , respectively.This is the same as used in experiments with steps [17].Please note that flux in the experiment is measured at the outlet for the precursors TEGa and O 2 , whereas the flux in KMC is given by the flux to the surface.The flux at the surface is unknown in the experiment.Because the transport of TEGa and O 2 through the argon atmosphere is different also the ratios at the inlets and at the surface will be different.In Fig. 7, the surface morphology can be seen from AFM images after 10, 20, and 25 s.Please note that every image is the result of another experiment.This has to be taken into account when comparing the results with those of the numerical calculations.
For R ads O 2 /Ga = 20 (A2b) we performed runs at the three temperatures of the experiment.The results are shown in Fig. 8. Here, we also present the coverage observed in the shots of the experiments (black curves).The time of the experiment was scaled by a factor of 0.1 to be consistent with the KMC results.In addition, the starting point was shifted by 0.5 s (on the scale of Fig. 8).The growth in the experiment was about a factor 10 slower than in the numerical runs, which might be mainly due to the choice of impingement rates in the simulation.As mentioned above in the experiment the FIG. 6. Number of Ga atoms per unit cell area (bc) as a function of time for different runs on a flat surface.At the initial time, the number is defined as 0. Bluish lines represent runs with the domain size of 20 × 32 cells.Red lines represent runs with the domain sizes 64 × 64 and 40 × 100 for A2c and A2b, respectively.For case A2b and domain size 40 × 100, four runs with different initialization for the random generator have been performed resulting in quite different times for the onset of nucleation.For comparison, we also show the dynamics for the step growth of case A2b (A2b step).For an increased desorption energy growth on a flat surface is also observed for the case R ads O 2 = 50 s −1 (A3a).impingement rate is unknown and cannot directly be related to the incoming flux of precursors.Later we will discuss the results for smaller impingement rates which exhibit similar differences in growth for the three temperatures (750 • C, 800 • C, and 850 • C).
In Fig. 8, the curves for the experimental growth rates agree quite well with those of the KMC computations.Please keep in mind that layer growth is a statistical process and every experimental run the situation might be slightly different.One drawback of KMC is that currently the domain is much smaller than the analysed area of experiment (500 × 500 nm 2 ).The large domain (40 × 100) corresponds to an area of 23.2 × 30.4 nm 2 (size is indicated in the lower left image of Fig. 7 by a white box).Therefore, the surface morphology can be compared only up to a certain limit.The morphologies of the KMC runs are shown on top in Fig. 8.For T = 850 • C, one long island is evolving.The island is shown at 2.5 s.For the same coverage there are many smaller islands in the case of T = 800 • C and 750 • C. At the time of 2.5 s there already many overlapping 2D islands and we are not far away from full coverage.The island in the run at T = 850 • C has about 70 atoms in b direction and about 15 atoms in c direction, which means a relation of 21.0 to 8.7 nm, i.e., about 2.4:1.The elongation in the experimental results is more pronounced.Because we have periodic boundary conditions in our calculations the 2D island in the run for T = 850 • C starts to interact with itself.The domain size should be increased but there are restrictions due the cpu time.For the runs with 20 × 32 cells the average time step for an was about, 1 × 10 −8 , 2.5 × 10 −8 , 6 × 10 −9 , and 2.7 × 10 −10 s for A2a, A2b, A2c, and A2d, respectively.For runs A2d (20 × 32) and A2c (64 × 64) we needed about 58 days on a Xeon Gold 6134 cpu (3.2 GHz) for 4 s real time.For run A2b (20 × 32) a total cpu time of 35 h was required for 16 s real time.For runs with size 20 × 32, the cpu time per KMC update is 150-350 μs.For larger domains (64 × 64 and 40 × 100), the cpu time per update is much larger (up to 4000 μs) because of the much longer entry list of possible events.
Both in numerical and experimental runs there is only a small difference in the growth behavior at T = 750 • C and 800 • C.But the results at T = 850 • C differs significantly from those at the two other temperatures.A much longer time for forming a critical nucleus is needed, which can be seen both in KMC calculations (Fig. 8) and experiments (Fig. 7).The delay in forming a critical nucleus is of advantage when going to step flow mode because nucleation on the terrace is unwanted.
We discuss the 2D island at T = 850 • C in more detail.As mentioned above it is elongated in the b direction by a factor of about 2.4.Such a behaviour might be caused by anisotropic diffusion.Considering only the diffusion of single Ga atoms we have ratios between 2.5 and 2.8 for the number FIG. 8. Number of Ga atoms per unit cell area (bc) as a function of time for runs of A2b (R ads O 2 /Ga = 20) at different temperatures.Blue and red lines represent results from KMC runs with different number of cells as denoted in the upper left corner.Black lines are from coverage measurements of experimental results [24].These results are scaled by a factor of 1/10 in time to be comparable with the numerical results.For On top the morphologies are shown at 2 s for all temperatures and 4.5 s for T = 850 • C (visualisation made by OVITO [25]).b and c axes are vertical and horizontal, respectively.The size of the domain is 30.4 × 20.3 nm 2 (40 × 100).The Ga atoms a8 (top of B1 plane) are in light red, the Ga atoms a6 (top of B2-plane) are in blue, all other Ga atoms are in grey. of hoppings in b direction versus c direction.Taking also into account the diffusion of the clusters the ratio is less than one, i.e., more Ga is moving in average in c direction than in b direction.Thus, a difference in diffusion cannot explain the elongation.But anisotropy could be also apparent in the attachment dynamics.It is easier to stick to an island in b direction than in c direction.Thus, growth velocity is larger in 010 direction than in 001 direction.Please note that this happens for the attachment at an edge.The statement cannot be generalised for a comparison of attachment on {010} and {001} planes.
At this point we like to discuss the influence of cluster diffusion.As just mentioned cluster diffusion changes the relation between diffusion in b and c direction.More important is the impact on nucleation.A Ga atom can easily move on a B plane.Once oxygen atoms are attached to the Ga atom, the activation energy for a motion of the Ga atom alone will be much higher.Without a cluster diffusion we would easily have a nucleation on the B plane.Thus, cluster diffusion is essential to meet the growth dynamics in the experiments.
In runs B (see Table II), we use slightly higher energy barriers and a reduced impingement rate to provoke a slower growth rate.As already mentioned in the experiments a slower growth rate was observed than in the numerical runs A. Because of the complex kinetics it is not possible to make an easy scaling of all the rates: not all relations could be kept constant.Therefore it cannot be expected that exactly the same growth dynamics will occur on a different time scale.Nevertheless, the results of runs B are quite similar to those of runs A2b as can be seen in Fig. 9.At T = 850 • C, no critical nucleus is created for a domain of 20 × 32.For the larger domain a delay of more than 20 s is observed before growth is starting.A similar 2D island is formed like in the run A2b.For T = 800 • C and 750 • C, there is nearly no difference in the time evolution of the Ga adatoms for the runs with the small domain 20 × 32 and those with the large domain 40 × 100.This is related to the fact that the structures at these temperatures are much smaller than in the run for T = 850 • C.

B. Step flow
The formation of 2D island growth could be suppressed by using vicinal surfaces.Experimentally, this has been realised by providing (100) substrates with a miscut of 6 o towards the [001] direction [10,17].Because the structure of our KMC is based on unit cells we need to have at least two steps in the domain.A miscut of 6 o means a terrace containing 10 unit cells and thus we choose a domain of 20 × 32 cells.In order to obtain step flow, the diffusion towards the step edges should be fast enough compared to adsorption and desorption.In a first step, we check the influence of the desorption rate on the growth rate and surface morphology.We used three different extra energies E Ga des for desorption (see Table II, runs A1-A3).The growth rates are 0.18, 0.6, and 1.2 Ga atoms/(unit cell × s) for A1a, A2a, and A3a, respectively (see Fig. 10).As described in the last section, the desorption energy depends on the configuration of Ga 2 O on the surface.In the following, we use the energy for a free standing a5-a12-a5n (or a7-10-a7s) configuration as a reference.This means the desorption rates are R des Ga 2 O are 2.4 × 10 7 , 1.4 × 10 5 , and 7.7 × 10 2 s −1 for runs A1a, A2a, and A3a, respectively.In addition we performed runs with extreme values, i.e., E Ga des = 5 eV (run A4a, R des Ga 2 O = 1.5 × 10 −13 s −1 ) and E Ga des = 0 eV (run A0a, R des Ga 2 O = 4.1 × 10 9 s −1 ).The minimal sticking coefficient for Ga atoms is 0.15 (4 Ga atoms with S = 0.1 and 4 Ga atoms with S = 0.2 in a unit cell) which yields the right grey line in Fig. 10.Once there are already oxygen atoms on the surface the sticking coefficient for Ga becomes larger.Therefore the line for run A2b corresponds almost with the one for the minimal sticking coefficient though we have considerable desorption [in average about 0.17 Ga atoms per unit cell area (bc) and second].
Not only the growth rates are different but also the evolving surface structure.In Fig. 11, we show the surface morphology at 10 and 20 s.For an intermediate extra desorption energy (run A2a, E Ga des = 1.0 eV) we observe step flow mode.For a high extra desorption energy, we obtain growth on the terraces.The desorption rate is very low (0.03 per adsorption for A3a) or zero (for A4a) and this leads to an agglomeration of gallium and oxygen on the terraces and step flow mode is suppressed.However, eventually a step of one unit cell is forming (see top right of Fig. 11).Thus, in a post-mortem analysis of an experiment one will find a similar structure as in the case of step bunching.Please note, currently we can only place a step of one unit cell in the computational domain.
Thus we cannot answer the question if steps higher than one unit cell will develop by time.
A similar behavior has been observed in a recent series of experiments [12].Applying a higher chamber pressure and thus lowering the desorption leads to a surface structure with step heights of several unit cells and wide terraces.
In the case of low E Ga des (high desorption rate), we obtain step bunching (A1a and A0a).This can be explained as follows: if the desorption rate is too large any fluctuation at the edges will trigger an instability.Once one terrace becomes slightly wider by fluctuation of the edges more Ga atoms can desorb and less Ga atoms will reach the upper edge.Consequently, the terrace will become even wider.In other word, the process is self-accelerating and we observe step bunching.
The last point we like to discuss is the influence of the oxygen impingement rate on the surface morphology.Increasing the rate by a factor two with respect to the already discussed case A2a does not change the picture (A2b): we still have a step flow growth.Increasing the ratio further to R ads O 2 /Ga = 40 (A2c) leads to a growth on the terraces.There is no chance for the gallium atoms to move to the edges before oxygen is adsorbed and stick them together.On the other hand, if the oxygen impingement rate is small the growth rate is small.For most of the cases shown in Fig. 10, we performed runs with different initialisations of the random generator and different domain sizes 20 × 32 and 20 × 64.The variation is quite small and no significant difference can be observed.

V. CONCLUSIONS
We have developed a kinetic Monte Carlo model, which takes into account details of the kinetics during the homoepitaxy of Ga 2 O 3 .We have investigated the growth on a flat (100) surface and on a (100) plane with a miscut of 6 o .The energy barriers for diffusion were chosen in the range, which was determined by a 1D mean field approach adapted to experimental results of 2D island growth [10].We varied the relation between Ga and O 2 impingement rates as well as the extra energy for Ga 2 O desorption.Both have a significant influence on growth rate and surface morphology.
On a flat surface an increase of the ratio R ads O 2 /Ga of O 2 to Ga impingement rate leads to an increase in the growth rate.Also important is that a larger ratio shortened the time of creating a critical nucleus significantly.For a comparison with our experimental runs, we performed numerical calculation for different temperatures.At the substrate temperature of T = 850 • C a large island is formed after a certain time of creating a critical nucleus.For lower temperatures (T = 800 • C and 750 • C) the creation of a critical nucleus occur much faster and many small islands are formed.The same behavior is observed in experiments.In both experiments and numerical calculations, a 3D growth is not observed: all islands grow together until forming the next layer of a B plane.
On a vicinal surface with a 6 o miscut an increase of the ratio O 2 to Ga impingement rate beyond 40 leads to a growth on terraces.The oxygen atoms bind to the gallium atoms on the terraces and the probability creating a critical nucleus is high compared to the probability that gallium or small cluster can reach the edge.Thus, in this case a step flow mode is not possible.Also the desorption rate has a significant influence on the growth kinetics.There is only a small window where step flow growth takes place.If the desorption is too low nucleation and growth on the terraces will take place.However, in contrary to the above case with large R ads O 2 /Ga , we observe a step of a full unit cell at later stage, i.e., the result is the same as in the case of step bunching.This is in accordance with the result of recent experiments [12].If the desorption rate is high any fluctuation of the edge results in an instability and we obtain step bunching.In an intermediate range, we get a step flow growth as wanted for producing layers of high quality.
The results of a large set of calculations with different parameters give evidence that there is only a small parameter window where a step flow growth with a sufficient growth rate will take place.This correlates with the experience from a large number of MOVPE experiments.

FIG. 2 .
FIG. 2. Atom labeling in the Unit cell.(a) c-a * plane.Please note that the crystallographic c axis is in the minus x direction .(b) b-a * plane.For a better viewing, we also show one neighboring cell.
and T are the activation energy, the Boltzmann constant, and temperature, respectively.a a * = 11.88Å and a b = 3.04 Å are the lattice constants in a * and b directions.The Debye temperature was determined by thermal transport measurements: θ D = 145(5) K [19].Together with the temperature T = 850 • C used in the above mentioned experiments this yields E A ≈ 1.55 eV.

FIG. 5 .
FIG. 5. Configurations for diffusion of clusters.For every configuration the activation energy E A is given.E A might be increased by additional bondings in the layer of the cluster.The possibilities are shown for the Ga atom a5.The same possibilities are set for a7.

FIG. 7 .
FIG. 7. MOVPE growth on a flat (100) plane for three different substrate temperatures.AFM images of the results of four runs stopped after 10 s [(a), (e), and (i)], 15 s [(b), (f), and (j)], 20 s [(c), (g), and (k)], and 25 s [(d), (h), and (l)].b and c axes are vertical and horizontal, respectively.The size of the domain is 500 × 500 nm 2 .For comparison in the bottom left image the size of the computational domain (40 × 100) is shown as a white box.The residual mean square (RMS) value of the roughness is defined as 1 n

FIG. 9 .
FIG. 9. Number of Ga atoms per unit cell area (bc) as a function of time for runs of case B at different temperatures.Blue and red lines represent runs with the domain size of 20 × 32 cells and 40 × 100, respectively.For T = 750 • C and 850 • C, runs with different initialization of the random generator have been performed.For T = 850 • C, in total four runs have been performed.In one case, no nucleation have been observed until 30 s.For all temperatures the morphologies are shown on top for the same coverage of 0.45 Ga adatoms per unit cell area: (a) T = 750 • C, (b) 800 • C, and (c) 850 • C.

FIG. 10 .
FIG. 10.Number of Ga atoms per unit cell area (bc) as a function of time.At the initial time the number is defined as 0. The results for three different E Ga des are shown.Bluish lines represent runs with the domain size of 20 × 32 cells.Red lines represent runs with the domain size 20 × 64.For all cases we performed runs with different initialisations of the random number generator and show the average curves including the RMS errors.For comparison also the case A2b is shown, where R ads O 2 is twice as for A2a.The lower dotted black line represents the growth rate for the minimal sticking coefficient for Ga atoms, the upper dotted black line represents the growth rate for S = 1.

FIG. 11 .
FIG. 11.Surface morphologies for different runs at 10 s (left) and 20 s (right).The desorption rate of Ga 2 O is increasing from top to bottom.The Ga atoms a8 (top of B1 plane) are in light red, the Ga atoms a6 (top of B2 plane) are in blue, all other Ga atoms are in grey.Oxygen atoms have smaller radii and are in red.Visualisation made by OVITO [25].
) (i.e., E dif A = E bond + E extra and E dif A = E bond + E Ga edge + E extra ):

TABLE I .
Activation energies E 1 des for desorption of Ga.

TABLE II .
Parameter for runs.Energies E are given in eV and rates R in s −1 .The values in the top table are valid for all runs A an B.