Study on the optimization of the deposition rate of planetary GaN-MOCVD films based on CFD simulation and the corresponding surface model

Metal-organic chemical vapour deposition (MOCVD) is a key technique for fabricating GaN thin film structures for light-emitting and semiconductor laser diodes. Film uniformity is an important index to measure equipment performance and chip processes. This paper introduces a method to improve the quality of thin films by optimizing the rotation speed of different substrates of a model consisting of a planetary with seven 6-inch wafers for the planetary GaN-MOCVD. A numerical solution to the transient state at low pressure is obtained using computational fluid dynamics. To evaluate the role of the different zone speeds on the growth uniformity, single factor analysis is introduced. The results show that the growth rate and uniformity are strongly related to the rotational speed. Next, a response surface model was constructed by using the variables and the corresponding simulation results. The optimized combination of the matching of different speeds is also proposed as a useful reference for applications in industry, obtained by a response surface model and genetic algorithm with a balance between the growth rate and the growth uniformity. This method can save time, and the optimization can obtain the most uniform and highest thin film quality.

Metal-organic chemical vapour deposition (MOCVD) is a key technique for fabricating GaN thin film structures for lightemitting and semiconductor laser diodes. Film uniformity is an important index to measure equipment performance and chip processes. This paper introduces a method to improve the quality of thin films by optimizing the rotation speed of different substrates of a model consisting of a planetary with seven 6-inch wafers for the planetary GaN-MOCVD. A numerical solution to the transient state at low pressure is obtained using computational fluid dynamics. To evaluate the role of the different zone speeds on the growth uniformity, single factor analysis is introduced. The results show that the growth rate and uniformity are strongly related to the rotational speed. Next, a response surface model was constructed by using the variables and the corresponding simulation results. The optimized combination of the matching of different speeds is also proposed as a useful reference for applications in industry, obtained by a response surface model and genetic algorithm with a balance between the growth rate and the growth uniformity. This method can save time, and the optimization can obtain the most uniform and highest thin film quality.
2018 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.

Introduction
GaN-based III-V compound semiconductor thin films are important third-generation semiconductor materials, which are widely used in the manufacture of blue-violet light-emitting diodes, semiconductor lasers and high-frequency, high-power electronic equipment [1][2][3][4][5][6]. Metal-organic chemical vapour deposition (MOCVD) is an important technology for growing GaN thin films due to its advantages such as the large epitaxy area, good reproducibility and high deposition rate [7][8][9][10][11].
Film thickness uniformity is an important index to measure the quality of film growth. In general, the more uniform the film thickness, the more uniform the distribution of components, and the higher the quality of the thin film. The quality of thin film is related to many parameters, such as the gas inlet flow rate, reactor pressure, rotation speed of the base and temperature of the substrate. In order to obtain uniform composition and film thickness, gas flows are required to be laminar. In such a laminar flow, vortices do not occur near the substrate, the gas flow in the direction parallel to the substrate has a uniform temperature field in the direction vertical to the substrate with a high temperature difference, and the particle concentration near the substrate should be as uniform as possible [12][13][14][15][16]. In general, coefficient variation of less than 5% is considered to be important for usable thin films. However, owing to the actual flow state in the reaction chamber, the temperature field and chemical composition of the reaction chamber cannot be visualized by the equipment operator, which introduces unknowns in the process of growth and adjustment.
Fortunately, computational fluid dynamics (CFD) is highly useful for the development of MOCVD equipment. Dimitrios et al., in an early study of the cavity flow state [17][18][19], investigated the particle flow of TiO in a vertical spray reaction chamber, and using the fluid model calculation results and particles in the reaction chamber of the control flow graph experiment, obtained the flow state to determine the feasibility of the numerical simulation. In addition, Theodoropoulos and others investigated the chemical mechanism model [20,21] to study the reaction chamber and coupled field, and proposed a model for the influence of different process parameters on the deposition rate under different conditions. With the improvement of GaN-MOCVD growth mechanism models, scholars such as Mitrovic studied the influence of different process parameters on deposition rate and optimized the design of reaction chambers [22][23][24][25][26][27][28][29][30][31].
At present, the planetary MOCVD chamber can adequately solve the problem relating to the thickness uniformity of large volume and large capacity film deposition. Planetary MOCVD is conducive to large size, large modulus and production of thin films, because it can be operated at low rotational speeds; this avoids large centrifugal forces that lead to film delamination. However, there is little research on the influence of the coupling between different rotation speeds on the flow rate and deposition rate. However, a few research groups have theoretically studied the method to adjust the speed of the different regions to find the optimum deposition rate. It is a difficult problem to adjust the uniformity of the film because of the complex coupling field. The simplest method is to adjust film deposition rate by rotating speed.
This study uses the rotation speeds of the substrate and wafers as variable parameters, and other conditions remain unchanged. The quality of MOCVD films was simulated and optimized, and the variation coefficient of thin film thicknesses was used as a measure of the uniformity and the quality of thin film growth. We present the rotation speeds of substrate and wafers as the input variables, the deposition rates as outputs, and the variation coefficient of the deposition rates are used as a response. First, 100 sets of input and output values are obtained by the design of experiment (DOE) method, and then the response surface model is constructed; finally, optimal results are obtained by using a response surface model (RSM) and a genetic algorithm.

Governing equations
GaN film growth occurs in a state of low pressure, laminar flow and steady state during the growth of MOCVD. Its mathematical model includes the conservation of mass, momentum, energy and individual species.
The conservation equations and conservation equation of momentum are provided [32][33][34]: where ρ and v are the density and velocity vectors of the gas mixture, respectively, and p and g are the pressure and gravitational acceleration, respectively.τ is the shear stress tensor expressed bȳ where µ is the viscosity, and I is the unit tensor. The energy equation is given by where C p is the specific heat capacity at a constant pressure, k is the thermal conductivity and T is the temperature. In addition, H i , H i,0 , M i and J i represent the molar enthalpy, enthalpy of formation, molar mass and diffusion flux of the species, respectively. We employ a finite rate translation model to simulate the growth of GaN, which can be described by where w i is the mass fraction of species i and R i j is the net volumetric rate of creation of species i during reaction j.
The density of the gas mixture is defined by where P op , M w and R are the operating pressure, molecular weight of the gas and universal gas constant, respectively.

Reaction mechanisms
At present, the typical GaN-MOCVD is carried by the carrier gas ammonia (NH 3 ), carrying the metalorganic source precursor Ga(CH 3 ) 3 (TMG) and nitride through the top of the uniquely designed spray mixing device into the reaction chamber. From detailed studies by other researchers, it is known that the complex gas will undergo complex gas and surface chemical reactions in the MOCVD reaction chamber and produce a large number of intermediate products. This paper uses the existing GaN growth reaction mechanism, including five gas-phase reactions and five surface reactions to simulate the growth [20,21,31,35], as shown in tables 1 and 2, respectively. To evaluate the uniformity of film thickness, the coefficient of variation is introduced. It is an absolute value reflecting the degree of data dispersion, and measures the degree of variation of the observed data.   Table 2. Gas-surface interactions in metal-organic vapour phase epitaxy (MOVPE) of GaN. S denotes a free surface site and s = 1 denotes a unity sticking coefficient at zero coverage. wherex and x i are the average deposition rate and the growth rate of the ith location of the area, respectively, S i and S are the area of the ith location and the area of the zones between a and b, respectively and C is the coefficient of variation.

Geometry description
The simulation was carried out in an independently developed conceptual reactor, containing seven 6inch wafers for GaN-MOCVD, as shown in figure 1.   through vertical pipes. After the gas mixture flows through the heated substrate, these two components chemically react to produce monomer molecules, which collide, condense, crystallize and grow into a molecular cloud of particles to eventually form a heterostructured thin film on the substrate surface. The molecular cloud that is not deposited on the substrate leaves the reactor with the unused reactants, by-products and carrier gas through the exhaust flow channel on the sidewall. The turntable is rotated in three areas, with A being the rotation speed of the substrate, B the central wafer rotation speed and C the rotation speed of the surrounding six wafers, as labelled in figure 1a. In order to ensure uniform deposition rate, the substrate and the seven wafers rotate at different speeds. The deposition rate is affected by many factors, such as the gas inlet velocity, substrate rotation speed and pressure in the reactor. In addition, the substrate is large, which creates difficulties in achieving uniform film thickness. Therefore, 12 points are selected on the substrate to study the speed regularity in different regions of the turntable. Average deposition rate along the circles is adopted to measure the thin film thickness.

Boundary conditions
As shown in figure 2, a hexahedral mesh and local refinement are used in areas where the flow rate and temperature gradient are large. We ensure numerical accuracy of the flow field in the reactor by maintaining the same mesh density in the intersection. Before we start the experiment, a grid sensitivity analysis is performed for the MOCVD reactor model to check its grid independency. The uniformity of film thickness is not a dimensioned value, and it can be indicated by the coefficient of variation. Generally, the smaller the coefficient of variation, the smaller the degree of numerical dispersion, and the more uniform the film thickness.
The system was modelled using the Fluent (v. 14.0) fluid dynamics software. The convergence criteria were that the residual value of the energy is less than 10 −6 , and the other variables are less than 10 −4 .
The following boundary conditions were imposed:   (4) The graphite plate has good thermal conductivity and the temperature difference within the susceptor is very small, considering that the temperature is constant. (5) A pressure-outlet boundary condition is adopted such that the given outlet static pressure is 0 Pa.

Simulation results under the reference conditions
The results in figure 3 present the flow field, temperature field and deposition rate of the MOCVD reaction chamber under the conditions given in table 1. Figure 3a shows the velocity profiles of three regions at different rotation speeds. The flow velocity at the inlet, thermal buoyancy and the flow field of each turntable of the substrate are shown in figure 3b. It can also be seen that the inner ring of the entrance acts to encourage gas deposition over the substrate, and the outer ring acts to stabilize the flow field. Figure 3c shows the temperature contours in the reaction chamber. The temperature boundary layer, which has a large temperature gradient, is within 5 mm of the substrate surface. The temperature distribution outside the boundary layer is relatively uniform, in the range of 1100-1273 K, which is suitable for the growth and deposition of thin films.  the central axis of the reactor. Combined with the two diagram observations, it is known that in the upper part of the substrate, because reaction G3 does not need activation energy, Ga(CH 3 ) 3 and NH 3 will quickly form adducts (G3) in the reaction chamber. The adduct is filled in the range of distance entrance 1 mm. When the adduct is near the high-temperature boundary layer, Ga(CH 3 ) 3 :NH 3 decomposes backward to generate Ga(CH 3 ) 3 and NH 3 (G4), and the maximum concentration is generated at a distance of 5 mm from the inlet, where the temperature is about 400 K. However, the undecomposed Ga(CH 3 ) 3 :NH 3 is diffused to the high-temperature area and removed CH 4 to produce an amino group Ga(CH3) 2 :NH 2 (G5). The maximum concentration occurred at the distance entrance 12.5 mm, where the temperature is about 900 K. The amino compounds continue to deposit here to produce GaN films. At the same time, when the temperature rises to the activation energy barrier, which can be crossed over in reactions G1 and G2, Ga(CH 3 ) 3 will be further decomposed, and Ga(CH 3 ) 2 and GaCH 3 will be generated near the substrate. However, the Ga(CH 3 ) 3 concentration of the product of direct decomposition of GaCH 3 is at least two orders of magnitude lower than the product of its addition reaction, Ga(CH3) 2 :NH 2 . The deposition rate of the substrate film is shown in figure 4h; however, it can be seen that the deposition rate is very uneven, with a film thickness variation coefficient of 5.94%. Therefore, we studied the speed of each region in order to obtain a more uniform deposition rate. increase in the speed, and the flow-field distribution on the base still keeps more uniform and straight, as shown in figure 5. When other benchmark conditions are unchanged, and when the rotation speeds of region B are 10, 15 and 20 r.p.m., the rotation speeds of different bases have little influence on the flow field in the reaction chamber. The distribution diagram of the flow field at the three speeds is basically similar, and the distribution of the flow field is even more uniform above the base, as shown in figure 6.

Influence
When other benchmark conditions are unchanged, and when the rotation speeds of region C are 30, 40 and 50 r.p.m., the rotation speeds of different bases have certain influence on the flow field of the reaction chamber. Three-speed flow-field distribution is slightly different in the intersection region of substrates A and C with the increase in the rotational speed fluctuations. With the increase of speed and cyclotron eddy above substrate C, it will have a certain impact on the film deposition. Above the base, a relatively uniform distribution flow is maintained, as shown in figure 7.
The different rotational speeds of the substrates of the three regions have little effect on the temperature field of the reaction chamber, which is basically consistent with that shown in figure 3c. (For further details, refer to electronic supplementary material, figures S1-S3.) Because the rotation speed does not increase significantly, it does not affect the temperature field in speed range.   In the three regions and at different rotational speeds, the concentration distribution of substances in the concentration boundary layer is Ga(CH 3 ) 2 :NH 2 , GaCH 3 and Ga(CH 3 ) 2 , and Ga(CH 3 ) 3 and Ga(CH 3 ) 3 :NH 3 disappear in the cavity. That is, the speed of changing the substrates of different regions does not change the rule that the Ga(CH 3 ) 3 and the NH 3 plus-pyrolysis paths dominate. The main growth precursors of the film growth are still Ga(CH 3 ) 2 :NH 2 and GaCH 3 .

Influence of rotational speed on species transport
In summary, the characteristics of planetary MOCVD involve material deposition at low rotational speeds. Therefore, changing the substrate speed in a relatively low range will not change the overall temperature field, flow field and component distribution. It will also not affect the dominance of addition and pyrolysis paths in MOCVD. It only affects the concentration change of the Ga gas component near the substrate rotation region, which further affects the change in the film deposition rate.

Influence of rotational speed on deposition rate
In the case of other conditions being unchanged, the deposition rate changed as the rotation speeds increase, as shown in figure 11. Figure 11a,b shows the deposition rate dependence on the rotation speed of region A. As the speed increases, the growth rate increases slightly, while the growth uniformity becomes slightly variable. This shows that the effect of substrate rotation speed on deposition rate and composition distribution is very small. Figure 11c,d shows the deposition rate dependence on the rotation speed of region B. The surrounding wafers are symmetrically distributed and rotate at the same speed, and the upper gas flow includes the central wafer to drive the uniform spread of the air distribution. As the speed increases, the growth rate and the growth uniformity again increase by a small amount.  spread around the substrate but also promotes an overall uniform airflow. It shows that the influence of central wafer rotation speed on the deposition rate distribution is relatively large. Figure 11 shows that in the intersection of the inner ring and the outer circle, due to the mismatch of the speed and the flow field, the deposition rate of the connecting parts in the two regions is significantly changed.
We propose that the problem of film thickness uniformity control is converted to the question: how to adjust the rotation speeds X to minimize the coefficient of variation of film uniformity.
Min C(Y): Y = f (X) 5 ≤ x 1 ≤ 20; 20 ≤ x 2 ≤ 50; 20 ≤ x 3 ≤ 50. (5.1) where, x i is the different zone speed, f (X) is the function between the inlet velocity and the average deposition rate, and C(Y) is the function of film thickness variation coefficient. In order to obtain the most uniform film thickness, the variation coefficient of the 12 average deposition rates must be minimized. If the other conditions remain unchanged, the optimal result is obtained by adjusting the values of A, B and C, as shown in table 4.
Numerical simulation takes a lot of time; thus, complete simulation using all the possible range of values of rotation speed will be a very large and time-consuming project. Therefore, in this study an   RSM method was used. One hundred sets of data were simulated by adjusting the values of A, B and C, and then these 100 sets of inputs (A, B and C) and outputs (12 deposition rates) were used to construct the input and output mapping relationship. Then 10 sets of new data were randomly selected to verify the accuracy of the model. Finally, a genetic algorithm was used to optimize the response surface model to find the most uniform (minimum coefficient of variation) result.

Construction and analysis of RSM model
The model has 12 deposition rate results; however, we can select just one of the results of the analysis due to the other results being similar. Figure 12a shows the normal probability plot of residuals, in which most of the points are in the normal distribution, and the fitting results of the model are very good. Figure 12b shows the effect of a single variable on the outcome: B value is the most influential, followed by A, with C having the least impact. The substrate rotation speed, A, has a small numerical change range, and is acting on the gas flow in the whole reactor, so that the change is relatively stable. The rotation speed of the central wafer, B, causes a vortex generated by rotation to spread throughout the reaction chamber, so that the local vortex drives the airflow distribution in the whole chamber to drive the reactant diffusion around, and thus has the greatest impact on deposition rate. The rotation speed, C, of the six circularly distributed wafers, and the induced vortices are symmetrical, resulting in the least impact. Thus, the influence of different rotating speeds on the flow field leads to the change of deposition rate of substrate film. Figure 13 shows the three-dimensional response surface model diagram, which illustrates the effect of the interaction of two parameters on the result. It not only shows the single factor effect, which is included in figure 12b, but also shows the interaction of two factors of A, B and C. From a single factor A surface alone, as the rotation speed increases, the fluctuation trends of A and figure 12b are in agreement; the same is true for factors B and C. The surface showing the A and C coupled effect is not stable, leading to the difference in deposition rate. However, the surfaces of B and C are coupled, and the A and B coupled effect is relatively stable.

Predication performance of RSM model and optimization result
Eight sets of new data were randomly selected, as shown in the electronic supplementary material. The RSM-predicted results were compared with the CFD simulation results and the average deposition rate error is no more than 5%. Therefore, the RSM-predicted and simulated values are consistent in both the numerical values and the changing trends, and for the subsequent optimization of the coefficient of variation, it is highly dependent on the accuracy of trend. Therefore, the RSM model can be used as the corresponding model for the input and output of the MOCVD reactor, which is used for the following optimization process.
Based on the corresponding relation of the input and output, the genetic algorithm is used to find the parameter value of the minimum coefficient of variation of the 12 deposition rates. The optimization results are shown in figure 14. The optimal rotation speeds are A = 40 r.p.m., B = 18 r.p.m. and C = 24 r.p.m., which provide a variation coefficient of 1.87%.

Conclusion
In this study, the MOCVD planetary and wafer rotation speeds were optimized for uniform film thickness, by both numerical simulation and using an RSM model. The important conclusions from this study are as follows.
(1) When the inlet flow rate, turntable temperature and internal pressure are stable, the influence of different speed regions on the temperature field is not obvious, but there is a certain disturbance in the speed and species, further affecting the deposition rate. Especially with the increase in the speed of region C, the disturbance is more obvious. (2) The accuracy of the RSM model is high, such that the predicted value and simulation value can be used to map the relationship between inputs and outputs, and it can be used to optimize the process parameters. The optimum rotation speed of different regions is obtained by using the optimized model, which provides a theoretical basis for obtaining the best film deposition rate. The coefficient of variation of film thickness across the entire substrate is reduced from 5.00 to 1.87%. (3) A model with a large number of grids is used to construct the RSM model, which can investigate the influence of each process parameter on the result, significantly reducing simulation time and cost.