Investigation of a Simpliﬁed Mechanism Model for Prediction of Gallium Nitride Thin Film Growth through Numerical Analysis

: A numerical procedure was performed to simplify the complicated mechanism of an epitaxial thin-ﬁlm growth process. In this study, three numerical mechanism models are presented for verifying the growth rate of the gallium nitride (GaN) mechanism. The mechanism models were developed through rate of production analysis. All of the results can be compared in one schematic diagram, and the differences among these three mechanisms are pronounced at high temperatures. The simpliﬁed reaction mechanisms were then used as input for a two-dimensional computational ﬂuid dynamics code FLUENT, enabling the accurate prediction of growth rates. Validation studies are presented for two types of laboratory-scale reactors (vertical and horizontal). A computational study including thermal and ﬂow ﬁeld was also performed to investigate the ﬂuid dynamic in those reactors. For each study, the predictions agree acceptably well with the experimental data, indicating the reasonable accuracy of the reaction mechanisms.


Introduction
Metal organic chemical vapor deposition (MOCVD) is a crucial process [1][2][3][4] for manufacturing compound semiconductor devices.Through the force convection mechanism, thin films are deposited on solid substrates from the gas inlet to the gas outlet.During this process, the reactant vapors enter the chamber from the top of the showerhead and flow vertically down into the reactor.Upon contacting the surface of the susceptor, the reactant vapors flow radially outward above the susceptor and finally exhaust through the outlet.By using MOCVD, a high-quality GaN epitaxial thin film was obtained for the first time on a sapphire wafer [5][6][7].
Design on simulation, which is based on commercial numerical software, is now widely used in the semiconductor industry in the optimization and troubleshooting of epitaxial growth processes that occur on the substrate [8,9].Furthermore, it is commonly employed in reactor design.In the past 20 years, few studies have investigated the GaN mechanism; therefore, the experimental data available for III-nitride MOCVD growth chemistry are extremely limited, particularly the chemical kinetic rate data in the Arrhenius equation.Although, a simplified reaction model was introduced by Safvi [10], and several features of the GaN growth chemistry have been established by experimental and numerical computational studies [10][11][12][13][14][15][16][17][18].Low temperature kinetic effects in MOVPE of III-V compounds and a novel approach to simulation of the group-III nitrides growth method were presented by Yakovlev et al.The growth rate and the solid phase composition are predicted theoretically and compared with available experimental data, even the application in a commercial Axitron reactor [19][20][21].
In the previous year, a very systematic and exhaustive work about a hybrid kinetic Monte Carlo and an exponentially-weighted moving average (EWMA) control algorithm for a multi-scale PECVD growth was developed by Crose et al. [22].In order to improve uniformity in the presence of spatially non-uniform species concentrations, it shows that the thin films can be driven to the desired thickness in four wafer zones.In the same year, Crose et al. developed a multi-scale CFD model to understand the plasma chemistry and transport phenomena in a PECVD reactor [23].The mass, momentum and energy equation have been solved by discretization.All indicate that numerical modeling is a very helpful way for industry application.
By contrast, CFD modeling has recently become a critical tool used for designing and optimizing the performance of an MOCVD reactor.Flow stability maps have been developed on the basis of two-dimensional (2D) and three-dimensional numerical modeling by using possible process parameters.It demonstrates that all typical flow types encountered in a rotating-disc reactor capture the effects of all of the aforementioned process parameters on the flow stability in the reactor [24][25][26][27][28][29][30][31][32][33].
A complicated chemistry method derived from the ab initio quantum chemistry methods developed in 2005 is used for growth rate prediction [15].Such detailed reaction mechanisms, comprising mostly elementary and ring reactions, have the advantage of requiring little or no calibration by using experimental data.However, for a complicated mechanism, a single computer has difficulty in rapidly obtaining a solution, and the calculation time is considerably long.
Although a complicated mechanism provides sufficient accuracy for the calculations, and the computational capacity is extremely high, calculation requires several days, and using a standard computer results in computational difficulties.Parallel computing may be the most effective means for solving this problem, but the cost of this means is too high.This study used the complicated chemical mechanism from previous studies [15] and simplified it through rate of production (ROP) analysis.The optimization reaction mechanism model was then employed as input to a 2D CFD code, enabling the accurate prediction of growth rates for GaN thin film (FLUENT, 2D calculation).This study developed a simplified mechanism that can be easily calculated using a standard computer to predict the GaN growth in the blue-LED process in the high temperature of mass transport zone (around 1050 • C), which was not developed in the past in the literature.It can be used as the key equipment design mechanism in academia and industry.

Materials and Methods
Predicting the GaN growth rate and thin film uniformity at the surface or wafer requires a coupled solution of the conservation of mass, momentum, energy and species equations.In other words, finite-rate homogeneous chemical reactions must be considered.For these purposes, the commercial CFD code ANSYS FLUENT was used.The proposed mechanisms for GaN growth were employed in a series of simulations of two reactors.

Procedure for the Numerical Modeling
Figure 1 illustrates the complete procedure for the numerical modeling.The first step was 0-D computation, including an ROP-optimized analysis for the mechanism model (complete model, reduced model and new model) and general theory.In the next step, we imported the most simplified model (new model) for the 2-D software (Ansys 14.5, Canonsburg, PA, USA), FLUENT, to calculate the growth rate in two examples.Figure 2 briefly illustrates the procedure for the development of the growth rate; this figure describes how we obtained the growth rate for GaN thin film.

Introduction for 0-Dimensional Analysis
In CHEMKIN, the perfectly stirred reactor (PSR) model was used in this study.The contents of a PSR are assumed to be nearly spatially uniform because of the high diffusion rates of forced mixing.In other words, the rate of conversion of reactants to products is controlled by chemical reaction rates and not by mixing processes.Thus, we considered the reactor to be limited by reaction kinetics.Mass transport to the reactor walls was assumed to be infinitely fast.Therefore, the relative importance of surface reactions to gas-phase reactions was determined only by the surface-to-volume ratios of each

Introduction for 0-Dimensional Analysis
In CHEMKIN, the perfectly stirred reactor (PSR) model was used in this study.The contents of a PSR are assumed to be nearly spatially uniform because of the high diffusion rates of forced mixing.In other words, the rate of conversion of reactants to products is controlled by chemical reaction rates and not by mixing processes.Thus, we considered the reactor to be limited by reaction kinetics.Mass transport to the reactor walls was assumed to be infinitely fast.Therefore, the relative importance of surface reactions to gas-phase reactions was determined only by the surface-to-volume ratios of each

Introduction for 0-Dimensional Analysis
In CHEMKIN, the perfectly stirred reactor (PSR) model was used in this study.The contents of a PSR are assumed to be nearly spatially uniform because of the high diffusion rates of forced mixing.In other words, the rate of conversion of reactants to products is controlled by chemical reaction rates and not by mixing processes.Thus, we considered the reactor to be limited by reaction kinetics.Mass transport to the reactor walls was assumed to be infinitely fast.Therefore, the relative importance of surface reactions to gas-phase reactions was determined only by the surface-to-volume ratios of each material and the relative reaction rates.In this section, we introduce the equations that are used in this study.
• Global mass balance equation: m is the mass in the reactor; . m in, out is the mass flows in and out of the system; and .m reaction is the reaction term that includes generation and consumption work in both the gas and surface phases.

• ROP analysis:
In CHEMKIN, ROP analysis can help the user determine the most dominant reaction paths in complicated reaction mechanisms; the general expression for this analysis is: ROP analysis enables the rapid identification of dominant reaction paths.It indicates each path by a percentage on the x-axis (production or consumption), for example 0.1 = 10% and 0.2 = 20%.The y-axis shows the reaction details for different temperatures.R i,production is the production rate for species i; ∑ j R j,production is the total production rate for species i; R i,consumpation is the consumption rate for species i; and ∑ j R j,consumption is the total consumption rate for species i.

Introduction for 2-Dimensional (2-D) Analysis
In FLUENT computation, first, the laminar flow model is used for viscous purposes; second, the species transfer model is used for concentration diffusion; and third, the finite-rate model is used for chemical reaction.The steady-state, pressure-based solver is introduced for this study, so the computational time is the mean to the time for iteration.Thus, the time discretization ANSYS FLUENT is based on finite volume method; by default, ANSYS FLUENT stores discrete values of the scalar φ at the cell centers; however, face values φ f are required for the convection terms and must be interpolated from the cell center values.This is accomplished using an upwind scheme."Upwind" means that the face value φ f is derived from quantities in the cell upstream.The spatial discretization for pressure, momentum, all species and energy is the second-order upwind scheme; when the second-order accuracy is desired, quantities at cell faces are computed using a multidimensional linear reconstruction approach.In this approach, higher order accuracy is achieved at cell faces through a Taylor series expansion of the cell-centered solution about the cell centroid.Thus, when second-order upwind is selected, the face value φ f is computed using the following expression: where φ and ∇φ are the cell-centered value and its gradient in the upstream cell and → τ is the displacement vector from the upstream cell centroid to the face centroid.This formulation requires the determination of the gradient ∇φ in each cell.By the way, second-order upwind is available in the pressure-based and density-based solvers.
The following five equations listed systematically are for 2D computation in FLUENT, including the final calculation of the obtained growth rate of GaN.

• Continuity equation:
The general integral mass conservation equation in differential form is given by: Coatings 2017, 7, 43 → v is the velocity vector; and t denotes time.
• Momentum equation of motion: The general Navier-Stokes equation is: For an incompressible fluid, the left side is the momentum difference of the system, µ is the viscosity coefficient and → f is the body force.
• Energy equation: This equation involves the total work done by pressure and the work done by viscous force; k is the thermal conductivity for fluid; T is temperature; and e is the energy per unit mass for fluid.

• Species transportation equation:
This conservation equation is expressed in the following form: where R is the net rate of production of species by chemical reaction, Y is the local mass fraction of each species and S is the rate of creation by addition from the dispersed phase plus any user-defined sources.
. J is the diffusion flux of species.
• Laminar finite-rate equation and Arrhenius equation: The laminar finite-rate equation is written in the general form as follows: N is the number of chemical species in the system; ν i,r is the stoichiometric coefficient for reactant i in reaction r; ν i,r is the stoichiometric coefficient for product i in reaction r; M i is a symbol denoting species i; k f,r is the forward rate constant for reaction r; and k b,r is the backward rate constant for reaction r.
The forward rate constant (growth rate) is computed using the Arrhenius expression: where A is the pre-exponential factor, B is the temperature exponent, E a is the activation energy for the reaction and R is the universal gas constant.

2D Numerical Assumption
In this study, the following assumptions were made for the 2D model: (1) the gas mixture is treated as a continuum; (2) the steady state is considered; (3) gas behavior is axisymmetric; (4) radiation between the walls is neglected; and (5) incompressible gas behavior is ideal.

Boundary Condition
The simulations described in this section were performed through the commercial fluid dynamics code FLUENT, which employs the finite volume method to solve governing fluid equations of mass, momentum, energy and species conservation laws in a stationary state.The gas-flow conditions in various MOCVD reactors and their chemical conditions were calculated using FLUENT.The boundary conditions for the model are as follows.A no slip condition was imposed on the reactor wall, and a velocity boundary condition was specified as the velocity at the showerhead or nozzle gas inlet, which was obtained from the given volume flow rate.The temperature and species mole fraction were also specified for the inlet.For the outlet, the outflow boundary condition was selected; this condition assumed fully developed flow at the outlet.For the wall, zero velocity, constant temperature, a zero species gradient and zero thermal flux were assumed.

Computation of the Complete Mechanism Model
A mathematical representation was developed for the complete mechanism model through ROP analysis using the 0D software CHEMKIN (version 4.0, Reaction Design, San Diego, CA, USA).Tables 1 and 2 summarize the gas and surface reaction mechanisms [15].The Arrhenius coefficients and activation energy are all temperature sensitive; the small pressure difference does not affect the final conclusion and computation.Consequently, the pressure difference between the derivative and process is ignored.
The gas reaction mechanism considers several aspects, including the pyrolysis reactions for trimethyl gallium (TMG), dimethyl gallium (DMG) and monomethyl gallium (MMG) and the adduct reactions with NH 3 and CH 3 ; interaction among CH 3 , H 2 and H; and finally, the pyrolysis reaction for adducts.TMG and MMG decompose slowly because of their large activation energy barriers; however, the pyrolysis reaction of DMG to MMG is rapid.Therefore, the conversion of TMG to MMG is rapid at high temperatures, and the concentration of MMG is expected to be high at higher temperatures.Therefore, MMG and MMG:NH 3 are expected to be the dominant gallium-containing precursors at high temperatures, whereas TMG and TMG:NH 3 are expected to contribute to deposition at low temperatures.For the surface reaction, Table 2 lists three paths, and Table 3 shows the compositions for chemical compounds on the surface.2.98 0 Table 3.Chemical compositions for the chemical compounds on the surface [15].

Compounds Names Chemical Formula
COMPM1(S) Path 1 displays the ring reaction between MMG and NH 3 .First, the MMG in the gas phase is absorbed by the substrate to form the surface site MMG(S) that then interacts with NH 3 in the gas region and forms COMPM1(S).This procedure is repeated until COMPM5(S) is formed.COMPM5(S) is desorbed by a CH3 radical to form a ring structure RINGM1, which is then attached to a nearby surface Ga (S).Finally, the hydrogen atoms H 2 are eliminated, resulting in GaN formation.
Path 2 for GaN formation starts from adsorption at the surface site GaNH 2 (S).This species can directly be formed from the gas phase reaction of MMG and NH 3 or through the reaction of the adsorbed Ga(S) and adsorbed NH 2 (S).GaNH 2 (S) then absorbs the MMG and NH 3 in the gas phase to form the compound COMPMM1(S).Similar to Path 1, this process is repeated until the elimination of CH4 to produce a ring structure RINGM1(S).Finally, the ring absorbs CH 3 radicals from each Ga atom of the ring followed by the elimination of CH 4 , resulting in GaN formation.
In Path 3, TMG is adsorbed by N(S) to form TMG(s). NH 3 is then adsorbed by TMG(S) to form TCOM1(S).Moreover, TCOM1 can be formed by the adduct TMG:NH 3 in the gas phase, which is directly absorbed by the surface N(S).Similar to the processes in Paths 1 and 2, this procedure repeats until a complex TCOM3(S) is formed.Finally, this complex TCOM3(S) eliminates two molecules of CH 4 to produce GaN.To verify the accuracy of our CHEMKIN software, we incorporated the complete mechanism model presented in Tables 1 and 2 into the solver.Table 4 lists the input parameters for the calculation of the growth rate.Figure 3 illustrates the predicted result compared with the results of previous studies.On the basis of the aforementioned mechanisms, we observe that in a high-temperature zone, this numerical result precisely matches the experimental data.This section indicates that the numerical code CHEMKIN is a reliable software tool, and we used it to describe chemical mechanisms in this study.We describe the simplified procedure in the next section.
Coatings 2017, 7, 43 10 of 23 result precisely matches the experimental data.This section indicates that the numerical code CHEMKIN is a reliable software tool, and we used it to describe chemical mechanisms in this study.We describe the simplified procedure in the next section.

Development of the Complete Mechanism Model
For the LED industry, rapidly and accurately predicting epitaxial uniformity and the growth rate is critical.The accuracy is verified using the aforementioned 0-D calculation.For the complete model, the reaction number of the complete mechanism is up to 69, as listed in Tables 1 and 2. Therefore, we had to simplify the complete model.ROP analysis is used for determining the contribution of each reaction to the net production or destruction rate of a species.It is particularly useful for 0-D systems, where the computational expense for the added calculations is low, and the data can be considered from a large reactor set.Therefore, ROP analysis rapidly identifies dominant reaction paths.
Figure 4 displays the ROP analysis of a complete model and the reaction path at different temperatures for producing a GaN thin film.For a blue-LED process, the critical process temperature is 1323 K.For a process temperature greater than 1300 K, the most dominant reaction paths are Reaction Numbers 13 and 41 in Table 1.The percentage of the pyrolysis reaction for RINGM2(S) is extremely significant and can be up to 90%. Figure 5

Development of the Complete Mechanism Model
For the LED industry, rapidly and accurately predicting epitaxial uniformity and the growth rate is critical.The accuracy is verified using the aforementioned 0-D calculation.For the complete model, the reaction number of the complete mechanism is up to 69, as listed in Tables 1 and 2. Therefore, we had to simplify the complete model.ROP analysis is used for determining the contribution of each reaction to the net production or destruction rate of a species.It is particularly useful for 0-D systems, where the computational expense for the added calculations is low, and the data can be considered from a large reactor set.Therefore, ROP analysis rapidly identifies dominant reaction paths.
Figure 4 displays the ROP analysis of a complete model and the reaction path at different temperatures for producing a GaN thin film.For a blue-LED process, the critical process temperature is 1323 K.For a process temperature greater than 1300 K, the most dominant reaction paths are Reaction Numbers 13 and 41 in Table 1.The percentage of the pyrolysis reaction for RINGM2(S) is extremely significant and can be up to 90%. Figure 5 depicts the schematics for Paths 1 and 2. Path 1 displays the successive reaction between MMG and NH 3 .
The MMG in the gas phase is absorbed by the substrate to form the surface site MMG(S) and interacts with NH 3 until the ring structure RINGM2 is formed.The hydrogen atoms H 2 are then eliminated, resulting in GaN formation.Similar to that in Path 1, the process in Path 2 repeats until the elimination of CH 4 ; a ring structure RINGM2(S) is produced; however, the source for Path 2 is GaNH 2 .The ROP analysis in Figure 4 demonstrates that at high temperatures (1300-1400 K), the reaction for GaNH   The MMG in the gas phase is absorbed by the substrate to form the surface site MMG(S) and interacts with NH3 until the ring structure RINGM2 is formed.The hydrogen atoms H2 are then eliminated, resulting in GaN formation.Similar to that in Path 1, the process in Path 2 repeats until the elimination of CH4; a ring structure RINGM2(S) is produced; however, the source for Path 2 is GaNH2.The ROP analysis in Figure 4 demonstrates that at high temperatures (1300-1400 K), the reaction for GaNH2 in Path 2 almost does not occur.Figures 6-8 depict the ROP analyses for species consumption of TMG, DMG and MMG, respectively.The MMG in the gas phase is absorbed by the substrate to form the surface site MMG(S) and interacts with NH3 until the ring structure RINGM2 is formed.The hydrogen atoms H2 are then eliminated, resulting in GaN formation.Similar to that in Path 1, the process in Path 2 repeats until the elimination of CH4; a ring structure RINGM2(S) is produced; however, the source for Path 2 is GaNH2.The ROP analysis in Figure 4 demonstrates that at high temperatures (1300-1400 K), the reaction for GaNH2 in Path 2 almost does not occur.Figures 6-8 depict the ROP analyses for species consumption of TMG, DMG and MMG, respectively.Figures 6 and 7 display three consumption paths for TMG and DMG in the gas phase; however, at high temperatures, TMG tends to convert to DMG.A low percentage of TMG converts to TMG(S), and DMG tends to decompose into MMG in high-temperature zones.Moreover, as displayed in Figure 8, MMG is consumed in the ring reaction or converts to MMG(S) at the surface.To summarize, we can conclude that nearly no CH 3 exists in the gas phase for Path 2. Thus, Path 2 is not crucial in high-temperature zones.
From Sengupta [15], the upper reaction (ring structure) is the principal procedure for producing the GaN; every pathway is first connected by three Ga sources and three N sources, is then connected by another active point and finally releases the CH 4 to produce the GaN, but this is too complicated to calculate.Therefore, this study assumed a single-direction method of producing GaN: when the Ga source is adsorbed by the surface, the NH 3 in the gas phase is physically adsorbed on the surface; when the energy is sufficient, it releases CH 4 and H 2 to form GaN.
Hence, we deleted Path 2 to simplify the complete model.Aside from the elimination of the reaction path, the ring structure proposed by D. Sengupta [15] was assumed, and the surface reactions  2 reflect ring deposition.Every GaN deposition and formation route is connected to three Ga sources and three N sources, is then connected to active ring species and releases CH 4 gases.This assumption regarding ring deposition generated many surface medium species and led surface reactions to entail excessively complicated processes.Therefore, we assumed that GaN deposition occurs in a single direction.After the Ga source is absorbed by the surface, NH 3 in the gas phase is absorbed by the surface physically.When energy is sufficient, CH 4 and H 2 are released to form GaN.
Table 5 lists the new form of the surface mechanisms in the reduced model (the gas mechanism remains the same as that of the complicated mechanism) for Path 1. Figure 9 provides the final predicted result (including Path 1 and Path 3) of this new form of the mechanisms for the reduced model; the result strongly agrees with that of the complete model.Figures 6 and 7 display three consumption paths for TMG and DMG in the gas phase; however, at high temperatures, TMG tends to convert to DMG.A low percentage of TMG converts to TMG(S), and DMG tends to decompose into MMG in high-temperature zones.Moreover, as displayed in Figure 8, MMG is consumed in the ring reaction or converts to MMG(S) at the surface.To summarize, we can conclude that nearly no CH3 exists in the gas phase for Path 2. Thus, Path 2 is not crucial in high-temperature zones.
From Sengupta [15], the upper reaction (ring structure) is the principal procedure for producing the GaN; every pathway is first connected by three Ga sources and three N sources, is then connected by another active point and finally releases the CH4 to produce the GaN, but this is too complicated to calculate.Therefore, this study assumed a single-direction method of producing GaN: when the Ga source is adsorbed by the surface, the NH3 in the gas phase is physically adsorbed on the surface; when the energy is sufficient, it releases CH4 and H2 to form GaN.
Hence, we deleted Path 2 to simplify the complete model.Aside from the elimination of the reaction path, the ring structure proposed by D. Sengupta [15] was assumed, and the surface reactions from Route 5-Route 12, from Route 20-Route 34 and from Route 39-Route 40 in Table 2 reflect ring deposition.Every GaN deposition and formation route is connected to three Ga sources and three N sources, is then connected to active ring species and releases CH4 gases.This assumption regarding ring deposition generated many surface medium species and led surface reactions to entail excessively complicated processes.Therefore, we assumed that GaN deposition occurs in a single direction.After the Ga source is absorbed by the surface, NH3 in the gas phase is absorbed by the surface physically.When energy is sufficient, CH4 and H2 are released to form GaN.
Table 5 lists the new form of the surface mechanisms in the reduced model (the gas mechanism remains the same as that of the complicated mechanism) for Path 1. Figure 9 provides the final predicted result (including Path 1 and Path 3) of this new form of the mechanisms for the reduced model; the result strongly agrees with that of the complete model.

Mechanism Optimization
We investigated whether the mechanism (reduced model) could be simplified further.First, the adduct reaction in the gas phase is an ambiguous reaction for III/V compound deposition.It consumes Coatings 2017, 7, 43 15 of 23 the source for species III and results in a lower deposition rate for epitaxy thin film [16,17].Karpov [18] claimed that the adduct reaction usually occurs in a low-temperature area (gas inlet) in the MOCVD reactor.Moreover, the adduct compound tends to fully undergo pyrolysis, as described previously, at a high susceptor temperature (1273 K).The utilization of species III can be increased by adjusting certain thermal and flow mechanisms, such as the pressure in the reactor, uniformity of the thermal field, and inlet distribution of III/V gas.The abandoned species III gas can be improved by adjusting these factors.
According to the literature, the adduct compounds in the blue-LED deposition mechanism are TMG:NH 3 , DMG:NH 3 , MMG:NH 3 , DMG:NH 2 and MMG:NH 2 .To investigate the consuming path for each species III gas, an ROP analysis was conducted for TMG, DMG and MMG and is discussed in this section.
As illustrated in Figure 6, in a low-temperature zone (600-800 K), TMG tends to form an adduct with NH 3 .When the temperature is raised to a middle-temperature region (1000-1100 K), the surface adsorption reaction is the main consuming path for TMG.However, in a high-temperature zone, all of the TMG is rapidly consumed by DMG (1100-1400 K).Therefore, when the growth temperature is lower than 1100 K, the main path to the formation of the GaN thin film is entirely related to TMG.With the temperature increasing (1100-1400 K), TMG changes to DMG easily, and the main deposition path changes gradually.
As displayed in Figure 7, similar to TMG, in the high-temperature zone, all of the DMG is consumed by MMG, and the pyrolysis reaction is the main path; however, DMG tends to form an adduct with NH3 in the low-temperature region (<800 K).Moreover, as illustrated in Figure 8, as with TMG and DMG, in the low-temperature region, MMG tends to forming an adduct with NH 3 at a high temperature (900-1400 K) and participates in the ring reaction at the wafer surface to form GaN thin film.
The results of the study revealed that when the process temperature is lower than 800 K, Ga species tend to exist in the form of an adduct in the gas phase.The critical temperature for a blue-LED process is usually higher than 1000 K.However, few adducts exist in such a high-temperature area.Therefore, for a blue-LED process, we can perform the following simplification procedure.In a high-temperature zone, the activation energy is sufficiently small in the reaction from DMG to MMG that we assume that all TMG is consumed to immediately form MMG and that the adduct assumption is negligible.Therefore, for a gas reaction, the mechanism can be reduced to only one chemical reaction, and the reaction for the adduct is deleted; Table 6 presents the results (new model).

Numerical Analysis of the 2-D Model
In this section, we consider the new model mechanism as an input parameter for 2D calculation.Table 7 lists the input parameters and numerous boundary conditions for 2D computation; Figure 11a,b illustrates the reactor shapes for Experiments 1 and 2, respectively.

Numerical Analysis of the 2-D Model
In this section, we consider the new model mechanism as an input parameter for 2D calculation.Table 7 lists the input parameters and numerous boundary conditions for 2D computation; Figure 11a,b illustrates the reactor shapes for Experiments 1 and 2, respectively.Figure 12 is a schematic diagram of the mesh for each reactor.The first simulation experiment was performed for a simple vertical double-nozzle reactor, displayed in Figure 11a, for which the experimental data are provided in Table 7.The average growth rate is attained at 1273 K and a constant wall temperature of 300 K.A general method of providing the flow stability in MOCVD reactors is the use of dimensionless numbers that arise from the Navier-Stokes equations for incompressible flow.The Grashof number (Gr) characterizes the significance of natural convection in the reactor.In the preceding equations, flow stability criteria can be defined with some critical values of the dimensionless groups that represent the ratio of the relative strengths of different forces in the reactor.
The flow and temperature field in Experiment 1 are shown for the two cases.Figure 13a shows the flow field for Case 1, while Figure 14a shows the path lines for Case 2. The corresponding temperature profiles are shown in Figures 13b and 14b, respectively.There is one sidewall recirculation near the Figure 12 is a schematic diagram of the mesh for each reactor.The first simulation experiment was performed for a simple vertical double-nozzle reactor, displayed in Figure 11a, for which the experimental data are provided in Table 7.The average growth rate is attained at 1273 K and a constant wall temperature of 300 K. Figure 12 is a schematic diagram of the mesh for each reactor.The first simulation experiment was performed for a simple vertical double-nozzle reactor, displayed in Figure 11a, for which the experimental data are provided in Table 7.The average growth rate is attained at 1273 K and a constant wall temperature of 300 K.A general method of providing the flow stability in MOCVD reactors is the use of dimensionless numbers that arise from the Navier-Stokes equations for incompressible flow.The Grashof number (Gr) characterizes the significance of natural convection in the reactor.In the preceding equations, flow stability criteria can be defined with some critical values of the dimensionless groups that represent the ratio of the relative strengths of different forces in the reactor.
The flow and temperature field in Experiment 1 are shown for the two cases.Figure 13a shows the flow field for Case 1, while Figure 14a shows the path lines for Case 2. The corresponding temperature profiles are shown in Figures 13b and 14b, respectively.There is one sidewall recirculation near the A general method of providing the flow stability in MOCVD reactors is the use of dimensionless numbers that arise from the Navier-Stokes equations for incompressible flow.The Grashof number (Gr) characterizes the significance of natural convection in the reactor.In the preceding equations, flow stability criteria can be defined with some critical values of the dimensionless groups that represent the ratio of the relative strengths of different forces in the reactor.
The flow and temperature field in Experiment 1 are shown for the two cases.Figure 13a shows the flow field for Case 1, while Figure 14a shows the path lines for Case 2. The corresponding temperature profiles are shown in Figures 13b and 14b, respectively.There is one sidewall recirculation near the reactor wall in Figure 13a.We can clearly see that this may be contributed by the natural convection and the reactor shape.Moreover, Figure 13b illustrates the gradient of the temperature field over the susceptor at an operating pressure of 100 torr.We used the volume flow rate and assumed that the working gases are all ideal gases so that there is a feature P ∝ ρ 2 in a natural convection mechanism.
Coatings 2017, 7, 43 17 of 23 reactor wall in Figure 13a.We can clearly see that this may be contributed by the natural convection and the reactor shape.Moreover, Figure 13b illustrates the gradient of the temperature field over the susceptor at an operating pressure of 100 torr.We used the volume flow rate and assumed that the working gases are all ideal gases so that there is a feature P ∝ ρ 2 in a natural convection mechanism.The natural convection is a special natural phenomenon caused by thermal expansion and contraction that dominates the physical field if the inertial force effect is small.However, in the flow field in these two cases, the recirculation patterns are closer to the susceptor for higher inner feed flow rates.A high total flow rate (12 slm for both cases, inner tube feed of 5 slm for Case 1200 sccm for Case 2) in the reactor leads to recirculation above the susceptor.The recirculation for Figures 13a and  14a is all geometry driven and exists even at room temperature at the same process flow rates.
Removal of the recirculation through a lowering of the flow rate would allow more time for parasitic gas phase reactions, which would lead to lower growth rates.Lower flow rates would also lead to an increase in the temperature upstream of the susceptor, resulting in an increase in the premature decomposition of the reactants.Lower reactor pressures could also decrease the significance of recirculation at the expense of lower growth rates.The temperature profile for Figure 14b in Case 2 reactor wall in Figure 13a.We can clearly see that this may be contributed by the natural convection and the reactor shape.Moreover, Figure 13b illustrates the gradient of the temperature field over the susceptor at an operating pressure of 100 torr.We used the volume flow rate and assumed that the working gases are all ideal gases so that there is a feature P ∝ ρ 2 in a natural convection mechanism.The natural convection is a special natural phenomenon caused by thermal expansion and contraction that dominates the physical field if the inertial force effect is small.However, in the flow field in these two cases, the recirculation patterns are closer to the susceptor for higher inner feed flow rates.A high total flow rate (12 slm for both cases, inner tube feed of 5 slm for Case 1200 sccm for Case 2) in the reactor leads to recirculation above the susceptor.The recirculation for Figures 13a and  14a is all geometry driven and exists even at room temperature at the same process flow rates.
Removal of the recirculation through a lowering of the flow rate would allow more time for parasitic gas phase reactions, which would lead to lower growth rates.Lower flow rates would also lead to an increase in the temperature upstream of the susceptor, resulting in an increase in the premature decomposition of the reactants.Lower reactor pressures could also decrease the significance of recirculation at the expense of lower growth rates.The temperature profile for Figure 14b in Case 2 The natural convection is a special natural phenomenon caused by thermal expansion and contraction that dominates the physical field if the inertial force effect is small.However, in the flow field in these two cases, the recirculation patterns are closer to the susceptor for higher inner feed flow rates.A high total flow rate (12 slm for both cases, inner tube feed of 5 slm for Case 1200 sccm for Case 2) in the reactor leads to recirculation above the susceptor.The recirculation for Figures 13a and  14a is all geometry driven and exists even at room temperature at the same process flow rates.
Removal of the recirculation through a lowering of the flow rate would allow more time for parasitic gas phase reactions, which would lead to lower growth rates.Lower flow rates would also lead to an increase in the temperature upstream of the susceptor, resulting in an increase in the premature decomposition of the reactants.Lower reactor pressures could also decrease the significance of recirculation at the expense of lower growth rates.The temperature profile for Figure 14b in Case 2 is flat (uniform) above the substrate.For Case 1, the temperature profiles are compressed near the center of the substrate, because of the high velocity of the inner feed, while there is an expansion close to the edge of the substrate.The distributions of the thermal and flow fields are all in agreement with the results reported in the literature [10].
Figures 15 and 16 provide the growth results from Cases 1 and 2 in Experiment 1; the predicted growth rates correspond well with the experimental data.For Case 1, the species of Group III appear from the inner nozzle, and those for Group V appear from the exterior nozzle.This reactor configuration implies that the growth rate of the inner wafer is high and decreases from the wafer center to the outer edges.This is due to species depletion occurring as the reactant flows over the susceptor.
Coatings 2017, 7, 43 18 of 23 is flat (uniform) above the substrate.For Case 1, the temperature profiles are compressed near the center of the substrate, because of the high velocity of the inner feed, while there is an expansion close to the edge of the substrate.The distributions of the thermal and flow fields are all in agreement with the results reported in the literature [10].Figures 15 and 16 provide the growth results from Cases 1 and 2 in Experiment 1; the predicted growth rates correspond well with the experimental data.For Case 1, the species of Group III appear from the inner nozzle, and those for Group V appear from the exterior nozzle.This reactor configuration implies that the growth rate of the inner wafer is high and decreases from the wafer center to the outer edges.This is due to species depletion occurring as the reactant flows over the susceptor.For Case 2, where diffusion dominates the transport, the exact placement of the nozzle is inconsequential; the match is inadequate near the wall of the reactor, where the TMG jet directly impinges.This may be attributed to a lack of geometric details pertaining to both positions in the literature.
Experiment 2 was performed for a vertical showerhead reactor, displayed in Figure 11b, for which the initial and boundary conditions are provided in Table 7.The distributions of the thermal and flow fields are all in agreement with the results from the literature in Figure 17 [11].is flat (uniform) above the substrate.For Case 1, the temperature profiles are compressed near the center of the substrate, because of the high velocity of the inner feed, while there is an expansion close to the edge of the substrate.The distributions of the thermal and flow fields are all in agreement with the results reported in the literature [10].
Figures 15 and 16 provide the growth results from Cases 1 and 2 in Experiment 1; the predicted growth rates correspond well with the experimental data.For Case 1, the species of Group III appear from the inner nozzle, and those for Group V appear from the exterior nozzle.This reactor configuration implies that the growth rate of the inner wafer is high and decreases from the wafer center to the outer edges.This is due to species depletion occurring as the reactant flows over the susceptor.For Case 2, where diffusion dominates the transport, the exact placement of the nozzle is inconsequential; the match is inadequate near the wall of the reactor, where the TMG jet directly impinges.This may be attributed to a lack of geometric details pertaining to both positions in the literature.
Experiment 2 was performed for a vertical showerhead reactor, displayed in Figure 11b, for which the initial and boundary conditions are provided in Table 7.The distributions of the thermal and flow fields are all in agreement with the results from the literature in Figure 17 [11].For Case 2, where diffusion dominates the transport, the exact placement of the nozzle is inconsequential; the match is inadequate near the wall of the reactor, where the TMG jet directly impinges.This may be attributed to a lack of geometric details pertaining to both positions in the literature.
Experiment 2 was performed for a vertical showerhead reactor, displayed in Figure 11b, for which the initial and boundary conditions are provided in Table 7.The distributions of the thermal and flow fields are all in agreement with the results from the literature in Figure 17 [11].However, in this experiment, the rotation rate of the susceptor is very high (1200 RPM), and the flow field is unstable.Furthermore, the dominating mechanism for this experiment is rotation-induced inertial force convection.In Figure 18, the average growth rate is attained at temperatures of 1100-1300 K.The predicted growth rates correspond well with the experimental data at each temperature point (Figure 18).Experiment 3 was performed for a horizontal quartz reactor [13], displayed in Figure 19, for which the initial and boundary conditions are provided in Table 7.The average growth rate is attained at temperatures of 600-1400 K. Just likes Experiment 1 and 2, the predicted growth rates (Figure 20) correspond well with the experimental data in a blue-LED process temperature.How to get quick solution in mass-transport temperature is our focal issue; in the low temperature zone, the adduct consumption is deleted for this purpose, so the prediction is not quite correct in that zone, but that is not why we are concerned about it.Table 8 lists the calculation time; it can clearly be seen that the calculation time has dramatically dropped in the comparisons among these three models in the same conditions.However, in this experiment, the rotation rate of the susceptor is very high (1200 RPM), and the flow field is unstable.Furthermore, the dominating mechanism for this experiment is rotation-induced inertial force convection.In Figure 18, the average growth rate is attained at temperatures of 1100-1300 K.The predicted growth rates correspond well with the experimental data at each temperature point (Figure 18).Experiment 3 was performed for a horizontal quartz reactor [13], displayed in Figure 19, for which the initial and boundary conditions are provided in Table 7.The average growth rate is attained at temperatures of 600-1400 K. Just likes Experiment 1 and 2, the predicted growth rates (Figure 20) correspond well with the experimental data in a blue-LED process temperature.How to get quick solution in mass-transport temperature is our focal issue; in the low temperature zone, the adduct consumption is deleted for this purpose, so the prediction is not quite correct in that zone, but that is not why we are concerned about it.Table 8 lists the calculation time; it can clearly be seen that the calculation time has dramatically dropped in the comparisons among these three models in the same conditions.However, in this experiment, the rotation rate of the susceptor is very high (1200 RPM), and the flow field is unstable.Furthermore, the dominating mechanism for this experiment is rotation-induced inertial force convection.In Figure 18, the average growth rate is attained at temperatures of 1100-1300 K.The predicted growth rates correspond well with the experimental data at each temperature point (Figure 18).Experiment 3 was performed for a horizontal quartz reactor [13], displayed in Figure 19, for which the initial and boundary conditions are provided in Table 7.The average growth rate is attained at temperatures of 600-1400 K. Just likes Experiment 1 and 2, the predicted growth rates (Figure 20) correspond well with the experimental data in a blue-LED process temperature.How to get quick solution in mass-transport temperature is our focal issue; in the low temperature zone, the adduct consumption is deleted for this purpose, so the prediction is not quite correct in that zone, but that is not why we are concerned about it.Table 8 lists the calculation time; it can clearly be seen that the calculation time has dramatically dropped in the comparisons among these three models in the same conditions.

Conclusions
A simplified process that used CHEMKIN for modeling GaN growth was demonstrated.Several reaction pathways were calculated through ROP analysis.The activation barriers and vibrational frequencies were employed for calculating rate constants.The new model was used for calculating the growth rate for different reactors through CFD.The predicted growth rates correspond acceptably well with the experimental data.The most simplified mechanism can be used to predict the growth rate of a blue-LED process before an actual process begins; this is useful for LED companies.All of these observations demonstrate that the proposed mechanisms are crucial pathways for the epitaxial growth of GaN.This study provides a useful and predictable mechanism for GaN thin film research on the blue-LED process.The new information from this study is to provide a very usefully simplified mechanism for predicting the blue-LED process at a high temperature of the mass transport zone (around 1050 °C), which was not developed in the past in the literature.In the future, this mechanism will be imported into a human-machine interface in the industry.Before the R/D researchers start their own experiment, they can run the quick numerical computation in the interface, to see if the result is bad or good, saving the gas consumption.It can be provided not only for academia on the GaN mechanism research, but also for the R/D study of the blue-LED knowledge in industry.In spite of some  Table 8.The calculation time in three models.

Conclusions
A simplified process that used CHEMKIN for modeling GaN growth was demonstrated.Several reaction pathways were calculated through ROP analysis.The activation barriers and vibrational frequencies were employed for calculating rate constants.The new model was used for calculating the growth rate for different reactors through CFD.The predicted growth rates correspond acceptably well with the experimental data.The most simplified mechanism can be used to predict the growth rate of a blue-LED process before an actual process begins; this is useful for LED companies.All of these observations demonstrate that the proposed mechanisms are crucial pathways for the epitaxial growth of GaN.This study provides a useful and predictable mechanism for GaN thin film research on the blue-LED process.The new information from this study is to provide a very usefully simplified mechanism for predicting the blue-LED process at a high temperature of the mass transport zone (around 1050 °C), which was not developed in the past in the literature.In the future, this mechanism will be imported into a human-machine interface in the industry.Before the R/D researchers start their own experiment, they can run the quick numerical computation in the interface, to see if the result is bad or good, saving the gas consumption.It can be provided not only for academia on the GaN mechanism research, but also for the R/D study of the blue-LED knowledge in industry.In spite of some

Conclusions
A simplified process that used CHEMKIN for modeling GaN growth was demonstrated.Several reaction pathways were calculated through ROP analysis.The activation barriers and vibrational frequencies were employed for calculating rate constants.The new model was used for calculating the growth rate for different reactors through CFD.The predicted growth rates correspond acceptably well with the experimental data.The most simplified mechanism can be used to predict the growth rate of a blue-LED process before an actual process begins; this is useful for LED companies.All of these observations demonstrate that the proposed mechanisms are crucial pathways for the epitaxial growth of GaN.This study provides a useful and predictable mechanism for GaN thin film research on the blue-LED process.The new information from this study is to provide a very usefully simplified mechanism for predicting the blue-LED process at a high temperature of the mass transport zone (around 1050 • C), which was not developed in the past in the literature.In the future, this mechanism will be imported into a human-machine interface in the industry.Before the R/D researchers start their own experiment, they can run the quick numerical computation in the interface, to see if the result is bad or good, saving the gas consumption.It can be provided not only for academia on the GaN mechanism research, but also for the R/D study of the blue-LED knowledge in industry.In spite of some deviations in the low temperature zone, more verification work is underway through our own experiments.

Figure 1 .
Figure 1.Schematic procedure of the complete numerical modeling.

Figure 2 .
Figure 2. Schematic procedure for the 0-D and 2-D numerical analysis.

Figure 2 .
Figure 2. Schematic procedure for the 0-D and 2-D numerical analysis.

Figure 2 .
Figure 2. Schematic procedure for the 0-D and 2-D numerical analysis.

Figure 3 .
Figure 3. Result of the complete mechanism model vs. literature data [11].
depicts the schematics for Paths 1 and 2. Path 1 displays the successive reaction between MMG and NH3.

Figure 3 .
Figure 3. Result of the complete mechanism model vs. literature data [11].
2 in Path 2 almost does not occur.Figures 6-8 depict the ROP analyses for species consumption of TMG, DMG and MMG, respectively.

Figure. 4 .
Figure. 4. Rate of production analysis for the complete model.

Figure 4 . 23 Figure. 4 .
Figure 4. Rate of production analysis for the complete model.

Figure 6 .
Figure 6.ROP analysis for TMG in the gas phase.Figure 6. ROP analysis for TMG in the gas phase.

Figure 6 .
Figure 6.ROP analysis for TMG in the gas phase.

Figure 7 .
Figure 7. ROP analysis for DMGa in the gas phase.

Figure 8 .
Figure 8. ROP analysis for MMGa in the gas phase.

Figure 7 .
Figure 7. ROP analysis for DMGa in the gas phase.

Figure 6 .
Figure 6.ROP analysis for TMG in the gas phase.

Figure 7 .
Figure 7. ROP analysis for DMGa in the gas phase.

Figure 8 .
Figure 8. ROP analysis for MMGa in the gas phase.Figure 8. ROP analysis for MMGa in the gas phase.

Figure 8 .
Figure 8. ROP analysis for MMGa in the gas phase.Figure 8. ROP analysis for MMGa in the gas phase.

Figure 9 .
Figure 9. Complete model vs. the reduced model.

Figure 9 .
Figure 9. Complete model vs. the reduced model.

Figure 10
Figure 10 provides the final predicted results of the complete model, reduced model and new model and the experimental data of their comparisons.The difference between the new model and the complete model is small.

Figure 10
Figure 10 provides the final predicted results of the complete model, reduced model and new model and the experimental data of their comparisons.The difference between the new model and the complete model is small.

Figure 10 .
Figure 10.Comparison among all mechanism models.

Figure 13 .
Figure 13.Schematic diagram of for Experiment 1, Case 1 (a) flow field and (b) thermal field.

Figure 14 .
Figure 14.Schematic diagram of for Experiment 1, Case 2 (a) flow field and (b) thermal field.

Figure 13 .
Figure 13.Schematic diagram of for Experiment 1, Case 1 (a) flow field and (b) thermal field.

Figure 13 .
Figure 13.Schematic diagram of for Experiment 1, Case 1 (a) flow field and (b) thermal field.

Figure 14 .
Figure 14.Schematic diagram of for Experiment 1, Case 2 (a) flow field and (b) thermal field.

Figure 14 .
Figure 14.Schematic diagram of for Experiment 1, Case 2 (a) flow field and (b) thermal field.

Figure 17 .
Figure 17.Schematic diagram of for Experiment 2 (a) flow field and (b) thermal field.

Figure 17 .
Figure 17.Schematic diagram of for Experiment 2 (a) flow field and (b) thermal field.

Table 2 .
[15]ace-phase mechanism of the complete model[15](the unit of E a is Calorie).

Table 5 .
Surface mechanisms (Path 1 + Path 3) are considered for the reduced model.

Table 5 .
Surface mechanisms (Path 1 + Path 3) are considered for the reduced model.
k = A  e −Ea/RT A n Ea (Cal)

Table 6 .
New model for GaN deposition.

Table 7 .
Process parameters for 2D modeling.
Figure 10.Comparison among all mechanism models.

Table 7 .
Process parameters for 2D modeling.

Table 8 .
The calculation time in three models.

Table 8 .
The calculation time in three models.