Investigation of a Simplified Mechanism Model for Prediction of Gallium Nitride Thin Film Growth through Numerical Analysis

Abstract

A numerical procedure was performed to simplify the complicated mechanism of an epitaxial thin-film 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 simplified reaction mechanisms were then used as input for a two-dimensional computational fluid 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 flow field was also performed to investigate the fluid dynamic in those reactors. For each study, the predictions agree acceptably well with the experimental data, indicating the reasonable accuracy of the reaction mechanisms.

1. 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.

2. 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.

2.1. 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.

2.2. 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:

$$\frac{\mathrm{d}m}{\mathrm{d}t}={\dot{m}}_{\mathrm{in}}-{\dot{m}}_{\mathrm{out}}+{\dot{m}}_{\mathrm{reaction}}$$

(1)

m is the mass in the reactor; ${\dot{m}}_{\mathrm{in},\mathrm{out}}$ is the mass flows in and out of the system; and ${\dot{m}}_{\mathrm{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:

$${\mathrm{ROP}}_{\mathrm{i},\mathrm{production}}=\frac{R_{\mathrm{i},\mathrm{production}}}{{\sum }_{\mathrm{j}}{R}_{\mathrm{j},\mathrm{production}}}$$

(2)

$${\mathrm{ROP}}_{\mathrm{i},\mathrm{consumption}}=\frac{R_{\mathrm{i},\mathrm{consumption}}}{{\sum }_{\mathrm{j}}{R}_{\mathrm{j},\mathrm{consumption}}}$$

(3)

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; ∑_jR_j,production is the total production rate for species i; R_{i,consumpation} is the consumption rate for species i; and ∑_jR_{j,consumption} is the total consumption rate for species i.

2.3. 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:

$${\mathsf{\varphi }}_{\mathrm{f}}=\mathsf{\varphi }+\nabla \mathsf{\varphi }·\overrightarrow{\mathsf{\tau }}$$

(4)

where φ and ∇φ are the cell-centered value and its gradient in the upstream cell and $\overrightarrow{\mathsf{\tau }}$ 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:

$$\nabla ·\overrightarrow{v}=0$$

(5)

$\overrightarrow{v}$ is the velocity vector; and t denotes time.

• Momentum equation of motion:

The general Navier–Stokes equation is:

$$\mathsf{\rho }\frac{\partial e}{\partial t}+\mathsf{\rho }(\overrightarrow{v}·\nabla \overrightarrow{v})=-\nabla p+\mathsf{\mu }\nabla \overrightarrow{v}+\mathsf{\rho }·\overrightarrow{f}$$

(6)

For an incompressible fluid, the left side is the momentum difference of the system, μ is the viscosity coefficient and $\overrightarrow{f}$ is the body force.

• Energy equation:

$$\mathsf{\rho }\frac{\partial e}{\partial t}+\mathsf{\rho }(\overrightarrow{v}·\nabla e)=-p(\nabla ·\overrightarrow{v})+\nabla (k\nabla T)$$

(7)

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:

$$\frac{\partial }{\partial t}(\mathsf{\rho }Y)+\nabla (\mathsf{\rho }\overrightarrow{v}Y)=-\nabla ·\dot{\mathrm{J}}+R+S$$

(8)

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. $\dot{\mathrm{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:

$${\displaystyle \sum _{i=1}^N{{\mathsf{\nu }}^{\prime }}_{i,r}{M}_i}\underset{k_{b,r}}{\overset{k_{f,r}}{\rightleftarrows }}{\displaystyle \sum _{i=1}^N{{\mathsf{\nu }}^{″}}_{i,r}{M}_i}$$

(9)

N is the number of chemical species in the system; ${{\mathsf{\nu }}^{\prime }}_{\mathrm{i},\mathrm{r}}$ is the stoichiometric coefficient for reactant i in reaction r; ${{\mathsf{\nu }}^{″}}_{\mathrm{i},\mathrm{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:

$$k_{\mathrm{f},\mathrm{r}}={AT}^n{\mathrm{e}}^{-E\mathrm{a}/RT}$$

(10)

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.

2.3.1. 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.

2.3.2. 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.

3. Results and Discussion

3.1. 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).

Table 1. Gas-phase mechanism for the complete model [15] (the unit of E_a is Calorie). TMG, trimethyl gallium; DMG, dimethyl gallium; MMG, monomethyl gallium.

$\mathit{k}=\mathit{A}{\mathit{T}}^{\mathit{n}}{\mathbf{e}}^{-\mathit{E}\mathbf{a}/\mathit{RT}}$										A	n	Ea
G1	TMG	=	DMG	+	CH3					1.00 × 1047	−9.18	76,996
G2	DMG	=	MMG	+	CH3					7.67 × 1043	−9.8	34,017
G3	MMG	=	Ga	+	CH3					1.68 × 1030	−5.07	84,030
G4	TMG	+	NH3	→	TMG:NH3					2.28 × 1034	−8.31	3115
G5	TMG	+	NH3	→	DMG:NH2	+	CH4			1.70 × 104	2	19,969
G6	DMG	+	NH3	→	DMG:NH3					4.08 × 1031	−7.03	3234
G7	DMG	+	NH3	→	MMG:NH2	+	CH4			5.30 × 105	1.56	20,744
G8	MMG	+	NH3	→	MMG:NH3					7.95 × 1024	−5.21	2094
G9	MMG	+	NH3	→	GaNH2	+	CH4			8.10 × 105	1.3	17,722
G10	NH3	+	CH3	→	NH2	+	CH4			3.31 × 103	2.51	9859
G11	CH3	+	H2	→	CH4	+	H			1.20 × 1012	0	12,518
G12	TMG	+	H	→	DMG	+	CH4			5.00 × 1013	0	10,036
G13	DMG	+	H	→	MMG	+	CH4			5.00 × 1013	0	10,036
G14	TMG:NH3	→	MMG	+	2CH3	+	NH3			1.33 × 1044	−8.24	77,791
G15	CH3	+	H	+	M	→	CH4	+	NH3	2.40 × 1022	−1	0
G16	2CH3	=	C2H6							2.00 × 1013	0	0
G17	2H	+	M	=	H2	+	M			2.00 × 1016	0	0

and

Table 2. Surface-phase mechanism of the complete model [15] (the unit of E_a is Calorie).

Path 1$,\mathit{k}=\mathit{A}{\mathit{T}}^{\mathit{n}}{\mathbf{e}}^{-\mathit{E}\mathbf{a}/\mathit{RT}}$								A	n	Ea
1	MMG	+	N(S)	→	MMG(S)			1.16 × 105	2.98	0
2	MMG(S)	→	MMG	+	N(S)			1.12 × 1014	0.55	107,673
3	NH3	+	MMG(S)	→	COMPM1(S)			3.35 × 107	3.33	0
4	COMPM1(S)	→	NH3	+	MMG(S)			5.70 × 1013	−0.16	8146
5	MMG	+	COMPM1(S)	→	CH4	+	COMPM2(S)	1.23 × 1010	3.22	23,446
6	NH3	+	COMPM2(S)	→	COMPM3(S)			3.35 × 107	3.33	0
7	COMPM3(S)	→	NH3	+	COMPM2(S)			5.70 × 1013	−0.161	8146
8	MMG	+	COMPM3(S)	→	CH4	+	COMPM4(S)	1.23 × 1010	3.22	23,446
9	NH3	+	COMPM4(S)	→	COMPM5(S)			3.35 × 107	3.33	0
10	COMPM5(S)	→	NH3	+	COMPM4(S)			5.70 × 1013	−0.161	8146
11	COMPM5(S)	→	CH4	+	RINGM1(S)			1.23 × 107	3.22	23,446
12	Ga(S)	+	RINGM1(S)	→	RINGM2(S)	+	N(S)	3.35 × 107	3.33	0
13	RINGM2(S)	→	3H2	+	3GaN(B)	+	Ga(S)	3.68 × 109	2.05	59,610
Path 2$,\mathit{k}=\mathit{A}{\mathit{T}}^{\mathit{n}}{\mathbf{e}}^{-\mathit{E}\mathbf{a}/\mathit{RT}}$								A	n	Ea
14	CH3	+	Ga(S)	→	MMG(S)			1.76 × 109	1.39	0
15	MMG(S)	→	CH3	+	Ga(S)			4.54 × 1013	0.0346	79,480
16	NH2	+	Ga(S)	→	NH2(S)			3.17 × 108	1.83	0
17	GaNH2	+	N(S)	→	GaNH2(s)			2.27 × 106	2.247	0
18	GaNH2(S)	→	GaNH2	+	N(S)			4.83 × 1013	0.614	83,881
19	COMPMM1(S)	→	CH4	+	GaNH2(S)			1.49 × 1011	0.609	25,950
20	MMG	+	GaNH2(S)	→	COMPMM1(S)			1.16 × 105	2.98	0
21	NH3	+	COMPMM1(S)	→	COMPMM2(S)			3.35 × 107	3.33	0
22	COMPMM2(S)	→	CH4	+	COMPMM3(S)			1.49 × 1011	0.609	25,950
23	MMG	+	COMPMM3(S)	→	COMPMM4(S)			1.16 × 105	2.98	0
24	NH3	+	COMPMM4(S)	→	COMPMM5(S)			3.35 × 107	3.33	0
25	COMPMM5(S)	→	CH4	+	RINGM1(S)			1.49 × 1011	0.609	25,950
26	NH2(S)	→	NH2	+	Ga(S)			1.45 × 1014	0.09	59,786
27	COMPMM1(S)	→	MMG	+	GaNH2(S)			1.00 × 1014	0.55	42,819
28	COMPMM2(S)	→	NH3	+	COMPMM1(S)			5.70 × 1013	−0.1	8146
29	COMPMM4(S)	→	MMG	+	COMPMM3(S)			1.00 × 1014	0.55	42,819
30	COMPMM5(S)	→	NH3	+	COMPMM4(S)			5.70 × 1013	−0.1	8146
31	Ga	+	N(S)	→	Ga(S)			1.00 × 1011	1.5	0
32	Ga(S)	+	NH2(S)	→	GaNH2	+	Ga(S)	1.00 × 1025	0	0
33	Ga(S)	→	Ga	+	N(S)			1.00 × 1013	0	45,168
34	6CH3	+	RINGM2(S)	→	COM1(S)			7.55 × 107	2.31	0
35	COM1(S)	→	6CH3	+	RINGM2(S)			1.00 × 1013	0.71	45,506
36	COM1(S)	→	6CH4	+	3GaN(B)	+	Ga(S)	4.00 × 1012	0	49,675
Path 3$,\mathit{k}=\mathit{A}{\mathit{T}}^{\mathit{n}}{\mathbf{e}}^{-\mathit{E}\mathbf{a}/\mathit{RT}}$								A	n	Ea
37	TMG	+	N(S)	→	TMG(S)			1.16 × 105	2.98	0
38	NH3	+	TMG(S)	→	TCOM1(S)			3.35 × 107	3.33	0
39	TCOM1(S)	→	CH4	+	TCOM2(S)			1.49 × 1011	0.609	32,785
40	Ga(S)	+	TCOM2(S)	→	TCOM3(S)	+	N(S)	3.35 × 107	3.33	0
41	TCOM3(S)	→	2CH4	+	GaN(B)	+	Ga(S)	1.49 × 1011	0.609	49,675
42	TMG(S)	→	TMG	+	N(S)			1.12 × 1014	0.55	49,675
43	TCOM1(S)	→	NH3	+	TMG(S)			5.70 × 1013	−0.161	11,922
44	TMG:NH3	+	N(S)	→	TCOM1(S)			1.16 × 105	2.98	0
45	TCOM1(S)	→	TMG:NH3	+	N(S)			1.12 × 1014	0.55	49,675
46	TCOM1(S)	→	2CH3	+	MMG(S)	+	NH3 +N(S)	1.12 × 1014	0.55	10,7673
47	MMGNH3	+	N(S)	→	COMPM1(S)			1.16 × 105	2.98	0
48	COMPM1(S)	→	MMG:NH3	+	N(S)			1.12 × 1014	0.55	107,673
49	MMG:NH3	+	COMPM1(S)	→	CH4	+	COMPM3(S)	1.23 × 1010	3.22	23,446
50	MMG:NH3	+	COMPM3(S)	→	CH4	+	COMPM5(S)	1.23 × 1010	3.22	23,446
51	MMG:NH3	+	GaNH2(S)	→	COMPMM2(S)			1.16 × 105	2.98	0
52	MMG:NH3	+	COMPMM3(S)	→	COMPMM5(S)			1.16 × 105	2.98	0

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₃ and CH₃; interaction among CH₃, H₂ 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₃ are expected to be the dominant gallium-containing precursors at high temperatures, whereas TMG and TMG:NH₃ are expected to contribute to deposition at low temperatures. For the surface reaction, Table 2 lists three paths, and

Table 3. Chemical compositions for the chemical compounds on the surface [15].

Compounds Names	Chemical Formula
COMPM1(S)	NH3·MMG(S)
COMPM2(S)	Ga·NH2·MMG(S)
COMPM3(S)	NH3·Ga·NH2·MMG(S)
COMPM4(S)	Ga·NH2·Ga·NH2·MMG(S)
COMPM5(S)	NH3·Ga·NH2·Ga·NH2·MMG(S)
RINGM1(S)	NH2·Ga·NH2·Ga·NH2·Ga(S)
RINGM2(S)	(S)NH2·Ga·NH2·Ga·NH2·Ga(S)
COMPMM1(S)	MMG·GaNH2(S)
COMPMM2(S)	NH3·MMG·GaNH2·Ga(S)
COMPMM3(S)	NH2·Ga·NH2·Ga(S)
COMPMM4(S)	MMG·NH2·Ga·NH2·Ga(S)
COMPMM5(S)	NH3·MMG·NH2·Ga·NH2·Ga(S)
TCOM1(S)	NH3·TMG(S)
TCOM2(S)	NH2·DMG(S)
TCOM3(3)	(S)NH2·DMG(S)
COM1(S)	RINGM2(S)·CH3 complex

shows the compositions for chemical compounds on the surface.

Path 1 displays the ring reaction between MMG and NH₃. First, the MMG in the gas phase is absorbed by the substrate to form the surface site MMG(S) that then interacts with NH₃ 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₂ are eliminated, resulting in GaN formation.

Path 2 for GaN formation starts from adsorption at the surface site GaNH₂(S). This species can directly be formed from the gas phase reaction of MMG and NH₃ or through the reaction of the adsorbed Ga(S) and adsorbed NH₂(S). GaNH₂(S) then absorbs the MMG and NH₃ 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₃ radicals from each Ga atom of the ring followed by the elimination of CH₄, resulting in GaN formation.

In Path 3, TMG is adsorbed by N(S) to form TMG(s). NH₃ is then adsorbed by TMG(S) to form TCOM1(S). Moreover, TCOM1 can be formed by the adduct TMG:NH₃ 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₄ to produce GaN. To verify the accuracy of our CHEMKIN software, we incorporated the complete mechanism model presented in Table 1 and Table 2 into the solver.

Table 4. Process parameters for the 0-D model [11].

Volume Flow Rate	Values
Standard liter per min	–
TMGa flow rate (µmol/min)	7.48 × 10−4
NH3 flow rate (Slm)	2
H2 flow rate (Slm)	10
Parameter	value
Initial pressure (torr)	140
Susceptor temperature (K)	600–1400
Reactor initial temperature (K)	300
Gas inlet temperature (K)	300

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.

3.2. 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 Table 1 and Table 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₃.

The MMG in the gas phase is absorbed by the substrate to form the surface site MMG(S) and interacts with NH₃ until the ring structure RINGM2 is formed. The hydrogen atoms H₂ are then eliminated, resulting in GaN formation. Similar to that in Path 1, the process in Path 2 repeats until the elimination of CH₄; a ring structure RINGM2(S) is produced; however, the source for Path 2 is GaNH₂. The ROP analysis in Figure 4 demonstrates that at high temperatures (1300–1400 K), the reaction for GaNH₂ in Path 2 almost does not occur. Figure 6, Figure 7 and Figure 8 depict the ROP analyses for species consumption of TMG, DMG and MMG, respectively.

Figure 6 and Figure 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₃ 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₄ 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₃ in the gas phase is physically adsorbed on the surface; when the energy is sufficient, it releases CH₄ and H₂ 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 CH₄ 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₃ in the gas phase is absorbed by the surface physically. When energy is sufficient, CH₄ and H₂ are released to form GaN.

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

$\mathit{k}=\mathit{A}{\mathit{T}}^{\mathit{n}}{\mathbf{e}}^{-\mathit{E}\mathbf{a}/\mathit{RT}}$										A	n	Ea (Cal)
1	MMG	+	N(S)	→	MMG(S)					1.16 × 105	2.98	0
2	MMG(S)	→	MMG	+	N(S)					1.12 × 1014	0.55	107,673
3	NH3	+	MMG(S)	→	COMPM1(S)					3.35 × 107	3.33	0
4	COMPM1(S)	→	NH3	+	MMG(S)					5.70 × 1013	−0.16	8146
5	COMPM1(S)	→	CH4	+	COMPM2(S)					1.23 × 1010	3.22	23,446
6	COMPM2(S)	→	H2	+	GaN(B)	+	N(S)			3.68 × 109	2.05	59,610
7	TMG	+	N(S)	→	TMG(S)					1.16 × 105	2.98	0
8	TMG(S)	→	TMG	+	N(S)					1.12 × 1014	0.55	49,675
9	NH3	+	TMG(S)	→	TCOM1(S)					3.35 × 107	3.33	0
10	TCOM1(S)	→	NH3	+	TMG(S)					5.70 × 1013	−0.161	11,922
11	TCOM1(S)	→	CH4	+	TCOM2(S)					1.49 × 1011	0.609	32,785
12	TCOM2(S)	→	2CH4	+	GaN(B)	+	N(S)			1.49 × 1011	0.609	49,675
13	TMGNH3	+	N(S)	→	TCOM1(S)					1.16 × 105	2.98	0
14	TCOM1(S)	→	TMGNH3	+	N(S)					1.12 × 1014	0.55	49,675
15	TCOM1(S)	→	2CH3	+	MMG	+	NH3	+	N(S)	1.12 × 1014	0.55	107,673
16	MMGNH3	+	N(S)	→	COMPM1(S)					1.16 × 105	2.98	0
17	COMPM1(S)	→	MMGNH3	+	N(S)					1.12 × 1014	0.55	107,673

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.

3.3. 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 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₃, DMG:NH₃, MMG:NH₃, DMG:NH₂ and MMG:NH₂. 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₃. 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₃ 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. New model for GaN deposition.

Gas Reaction$,\mathit{k}=\mathit{A}{\mathit{T}}^{\mathit{n}}{\mathbf{e}}^{-\mathit{E}\mathbf{a}/\mathit{RT}}$								A	n	Ea (Cal)
1	TMG	→	MMG	+	C2H6			1.00 × 1047	−9.18	76,996
2	TMG	+	N(S)	→	TMG(S)			1.16 × 105	2.98	0
3	TMG(S)	→	TMG	+	N(S)			1.12 × 1014	0.55	49,675
4	NH3	+	TMG(S)	→	TCOM1(S)			3.35 × 107	3.33	0
5	TCOM1(S)	→	NH3	+	TMG(S)			5.70 × 1013	−0.161	11,922
6	TCOM1(S)	→	3CH4	+	GaN(B)	+	N(S)	1.49 × 1011	0.609	49,675
7	MMG	+	N(S)	→	MMG(S)			1.16 × 105	2.98	0
8	MMG(S)	→	MMG	+	N(S)			1.12 × 1014	0.55	107,673
9	NH3	+	MMG(S)	→	COMPM1(S)			3.35 × 107	3.33	0
10	COMPM1(S)	→	NH3	+	MMG(S)			5.70 × 1013	−0.16	8146
11	COMPM1(S)	→	H2	+	GaN(B)	+	N(S) + CH4	3.68 × 109	2.05	59,610

presents the results (new model).

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.

3.4. 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. Process parameters for 2D modeling.

Volume Flow Rate	2D Experiment 1, Case 1	2D Experiment 1, Case 2	2D Experiment 2	2D Experiment 3
Standard Condition
TMGa flow rate (Slm)	7.48 × 10−4	7.48 × 10−4	3 × 10−3	6.7 × 10−4
NH3 flow rate (Slm)	2	2	6	1
H2 flow rate (Slm)	10 (inner: 5 slm)	10 (inner: 0.2 slm)	7.1	5
Real Condition
TMGa flow rate (lm)	6.24 × 10−3	6.24 × 10−3	1.8 × 10−2	6.27 × 10−4
NH3 flow rate (lm)	16.7	16.7	36	0.9345
H2 floe rate (lm)	83.5	83.5	42	0.4672
Parameter	2D Experiment 1, Case 1	2D Experiment 1, Case 2	2D Experiment 2	2D Experiment 3
Initial pressure (torr)	100	100	140	780 torr
Susceptor temperature (K)	1273	1273	1100–1300	600–1400
Wall temperature (K)	300	300	300	300
Operating density (kg/m3)	0.0241	0.0241	0.0665	0.1882
Initial temperature (K)	300	300	300	300
Rotation rate (RPM)	0	0	1200	0

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 Figure 13b and Figure 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 ∝ ρ² 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 Figure 13a and Figure 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].

Figure 15 and Figure 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].

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. The calculation time in three models.

Complete Model	Reduced Model	New Model
10–12 h	4–6 h	0.2–0.5 h

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.

4. 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.