Beyond ab initio reaction simulator: an application to GaN metalorganic vapor phase epitaxy

To develop a quantitative reaction simulator, data assimilation was performed using high-resolution time-of-flight mass spectrometry (TOF-MS) data applied to GaN metalorganic vapor phase epitaxy system. Incorporating ab initio knowledge into the optimization successfully reproduces not only the concentration of CH$_4$ (an impurity precursor) as an objective variable but also known reaction pathways. The simulation results show significant production of GaH$_3$, a precursor of GaN, which has been difficult to detect in TOF-MS experiments. Our proposed approach is expected to be applicable to other applied physics fields that require quantitative prediction that goes beyond ab initio reaction rates.

(*Electronic mail: kusaba@riam.kyushu-u.ac.jp) (Dated: 24 October 2022) To develop a quantitative reaction simulator, data assimilation was performed using high-resolution time-of-flight mass spectrometry (TOF-MS) data applied to GaN metalorganic vapor phase epitaxy system. Incorporating ab initio knowledge into the optimization successfully reproduces not only the concentration of CH 4 (an impurity precursor) as an objective variable but also known reaction pathways. The simulation results show significant production of GaH 3 , a precursor of GaN, which has been difficult to detect in TOF-MS experiments. Our proposed approach is expected to be applicable to other applied physics fields that require quantitative prediction that goes beyond ab initio reaction rates.
Crystal growth technology of III-nitride semiconductors has been actively developed to fabricate optical and electronic devices 1-6 . Reactor simulators that solve for heat and mass transfer, chemical reactions, and phase transitions have played a role in these developments [7][8][9][10][11][12] . The reaction model and its kinetic parameters are key components of reactor simulators, especially for the reaction system for chemical vapor deposition. It is generally quite difficult to determine the complete kinetics of complex gas-phase reaction networks using experiments alone. Thus, such reaction systems have been studied by ab initio calculations, e.g., density functional theory (DFT) and transition state theory (TST) [13][14][15][16][17][18][19][20][21] . Increased semiconductor device performance has recently necessitated the use of simulators for quantitative process optimization that goes beyond obtaining a qualitative understanding [22][23][24][25] . However, the predictive performance of ab initio reaction rates is often inadequate for such quantitative requirements.
The aim of this study is to realize quantitative simulation by data assimilation for a reaction model with rate constant parameters for the GaN metalorganic vapor phase epitaxy (MOVPE) system. A complex reaction network generally involves many parameters, and experimental measurements are costly and provide limited data. Under these circumstances, two measures are taken to ensure that the solution is not indefinite or an overfit. First, a reduced reaction model is adopted that excludes many radical reactions and consists of fewer reactions than conventional models 21 . This reduced model has been demonstrated to reproduce the latest experimental fact on the reaction pathway that the removal of the methyl group in trimethylgallium (TMG) decomposition is caused by reactions with NH 3 26 . Second, the optimization is performed with the strategy that not only focuses on minimizing simulation errors, but also maintains the theoretical basis from ab initio calculations. In addition, the reaction pathway becomes a fundamental consideration in selecting a solution. This approach is used to develop a reliable and quantitative reaction simulator in this study.
Experimental data were obtained by high-resolution timeof-flight mass spectrometry (TOF-MS) [26][27][28] as the detection intensity of the individual molecules in MOVPE environments. This state-of-the-art mass spectrometric technique can distinguish between NH 2 and CH 4 , which is not possible using conventional quadrupole mass spectrometry (QMS) 29 . Indeed, this study relies on the CH 4 concentration data converted from its intensity data. Details of the data conversion are presented in Appendix A. Experimental measurements were performed dozens of times at different temperatures, but only data with kinetic information were used for data assimilation.
The input-output relationship of the simulator is concisely expressed asĉ (1) Here, the outputĉ(x) is the concentration distribution, the input T set is the heater temperature (an experimental condition), the simulation parameter k is the reaction rate constant, and the simulator f solves the ordinary differential equations of the reaction kinetics for our reduced reaction model shown in Table I. The reaction rate constant is determined using DFT: B3LYP/6-311G(d,p) and TST 21 . The modified Arrhenius format given below is used to model the temperature dependence of the ab initio rate constant and subsequently tuned to better reproduce the experimental data: Here, the parameters E act and A are the activation energy and preexponential factor obtained from ab initio calculations, and q E and q A are tuning coefficients for these parameters. Note that E act is not equal to the simple DFT activation energy because of the effects of TST. More details on the simulator can be found in Appendix B. A multiobjective optimization was performed to determine the k (i.e., q E and q A ) that maintains the theoretical basis obtained from ab initio calculations and reproduces the experimental data well. The ab initio ratios of E act were used as the theoretical basis. The concentration of CH 4 (an impurity precursor) was used as a measure of the simulation performance because CH 4 is a stable molecule and thus, the corresponding detection data are reliable and can be converted to concentration data by a reasonable scheme. Our optimization problem is expressed as follows: The first objective function is based on the relative error in the CH 4 concentration between the experimental (c CH 4 ) and simulation (ĉ CH 4 ) results at the detection position x d , where i is the index of the experimental conditions. In this case, only T set was varied. The second objective function is the variance of q E , where j is the reaction index and N is the total number of considered reactions. The smaller the variance is (i.e., the closer the coefficients q E j are to each other), the closer the ratios of the modified activation energies (q E j E act j /q E k E act k ) are to the ab initio ratios (E act j /E act k ). Multiobjective optimization determines the best trade-off solutions (i.e., Pareto solutions) rather than just one optimal solution. Evolutionary algorithms are widely used for finding Pareto solutions to multiobjective optimization problems. A well-known algorithm called the fast elitist nondominated sorting genetic algorithm (NSGA-II) was employed in this study 30 .
Although the parameters are softly constrained by the second objective function, the reaction model consisting of 29 reactions still has 58 parameters. To avoid difficulties in interpreting the tuning results, hard constraints are also imposed by requiring q E to be the same for each reaction group (Groups G1-G5, R7, and R14 in Table I). That is, the chemical equations belonging to the same group have the same reaction partners (H 2 or NH 3 ) and byproducts (CH 4 , NH 3 , or H 2 ) and differ only in the substituents (−CH 3 , −NH 2 , or −H) not involved in the reaction. Thus, these equations are expected to have the same degree of error in the ab initio E act . In addition, as the preexponential factor does not have as large an influence as the activation energy, the same q A was assumed for all reactions for simplicity. Figure 1(a) shows the Pareto solutions obtained from 30,000 function evaluations in the optimization. Here, a marker represents a set of tuning coefficients, i.e., q E for each reaction group and q A . These approximate Pareto solutions have the smallest variance of q E among the many possible sets of coefficients resulting in the same error. The solutions correspond to a trade-off between the two objective functions. That is, the larger the variance of q E is, the lower the root mean squared percentage error (RMSPE) of CH 4 concentration is. However, even with a solution that emphasizes lowering the variance of q E , the RMSPE is only less than 25%. In Fig. 1(b), the horizontal axis is q E and the vertical axis is the RMSPE. A marker in Fig. 1(a) corresponds to 7 markers aligned horizontally in Fig. 1(b) at the same RMSPE level. A solution with an RMSPE of approximately 25% has an almost identical q E of approximately 0.85 for each reaction group, that is, this solution maintains the ab initio ratios of E act . The solutions that emphasize reducing the RMSPE have different q E values for each reaction group. For example, the solution with an RMSPE of 15.0% presented in Fig. 1(a) has a set of q E s at a 15.0% RMSPE level in Fig. 1(b) (i.e., 0.795, 0.885, 0.846, 0.849, 0.848, 0.850 and 0.846 for G1, G2, G3, G4, G5, R7 and R14, respectively). To reduce the RMSPE, the q E for Group G1 needs to be relatively decreased, and the q E for Group G2 needs to be relatively increased. Note that q E E act for all reactions is an absolute value that has been decreased from the ab initio E act to obtain these Pareto solutions. In addition, q A exhibited fluctuations within 20% corresponding to discontinuous changes in q E . Figure 2 is a comparison of the simulation performance between the parameter sets tuned by the optimization and the original, untuned parameter set. The difference in the observed CH 4 concentration corresponds to the difference in the experimental conditions T set . Simulations were performed at each T set using a number of tuned parameter sets and the original parameter set. The results obtained using the tuned parameter sets are plotted in different colors depending on the variance of q E . The parameter sets with the largest and smallest q E variance among the parameter sets shown in Fig. 1 are hereafter referred to as the experimental and ab initio limits, respectively. That is, the experimental and ab initio limits are the solutions with the lowest error in the CH 4 concentration and that best preserve the ab initio E act ratios, respectively. The closer the markers in Fig. 2 are to the proportionality line, the better the simulation performance is. Using the original parameter set resulted in a considerable underestimate of the CH 4 concentration. Lowering E act by using the tuning factors shown in Fig. 1(b) promoted the TMG decomposition reactions (Groups G1, G2, R7, and R14) and increased CH 4 production up to near the observed values. It is remarkable that sufficiently quantitative simulations have been performed that represent a significant improvement over the original, even in the case of the ab initio limit.
Next, we discuss how a solution should be selected among the Pareto solutions. It is nontrivial to realize a proper balance between the two trade-off objective functions. Thus, we consider this problem from a different point of view than that considered for the optimization. . It can be seen that Reaction R8 dominates over Reaction R1 in both the original and ab initio limits. Reaction R1 becomes dominant in the experimental limit. On the other hand, it was reported that in a TOF-MS observation using a heavy hydrogen D 2 supply instead of H 2 , CH 4 or CH 3 D was dominantly detected for the condition with or without NH 3 , respectively 26 . This experimental result indicates that the main pathway involves TMG decomposition reactions with NH 3 (i.e., Group G2) and not reactions with H 2 (i.e., Group G1). In addition, this reaction pathway was qual- itatively supported by an ab initio study 21 . Thus, the ab initio limit solution should be selected to be consistent with the experimentally and theoretically supported reaction pathways. Note that the knowledge of the reaction pathway was not included in the optimization, although the results between the ab initio and experimental limits (indicated by color gradations) are trade-offs in terms of the two objective functions.
As a result of our quantitative reaction simulations, it is worth focusing on the concentration of GaH 3 because it is considered a precursor of the GaN film. Figure 4(a) shows the concentration distribution of GaH 3 at T set = 600 [°C]. It can be seen that GaH 3 is produced in a significant yield under TOF-MS conditions and is the main precursor considering that the TMG supply is 0.029 mol/m 3 . Note that Ga-related products are difficult to detect in TOF-MS experiments under high-temperature conditions, for some reasons 28 . That is, the complementary use of TOF-MS observations and dataassimilated reaction simulations can enable a more detailed study of gaseous compositions than previously possible. On the other hand, other precursor candidates GaH, GaCH 3 , and GaNH 2 were little produced. Nevertheless, even the rates of reactions producing these precursors are considered meaningful in the ab initio limit.
An experimental protocol for better tuning can be discussed in terms of the deviation of the simulated value.  Fig. 2. Immediately after the decomposition reactions begin and when the reactions are almost complete, the deviation of the simulation is small, then the observation data are relatively uninformative as shown in Fig. 4(c). In the transient or kinetic domain, the deviation of the simulation is large, and then the observation data are relatively informative as shown in Fig. 4(b). For more informative observation data, as T set rises, the detection position should be adjusted to the inlet side or the flow velocity should be adjusted faster. Using this protocol to collect more experimental data would enable more reliable rate constants to be determined.
In conclusion, data assimilation of the gas-phase reaction rate in GaN MOVPE was performed by multiobjective optimization. Simple optimization to reproduce the experimental CH 4 concentration did not yield reaction rates that reproduce known major reaction pathways. An approach was proposed to resolve this difficulty in which rate constants are tuned by imposing appropriate constraints based on ab initio calculations. The developed quantitative reaction simulator is a powerful tool for GaN MOVPE research. In addition, the proposed tuning approach is expected to be applicable to other applied physics fields that require quantitative prediction that goes beyond ab initio reaction rates.

ID
Reaction The CH 4 intensity resulting from Ga(CH 3 ) 3 decomposition under a H 2 atmosphere in the absence of NH 3 (see ). Thus, the CH 4 concentration is considered zero in the low-temperature domain and to be three times the Ga(CH 3 ) 3 supply in the equilibrium domain because the full decomposition reaction is Ga(CH 3 ) 3 + 3 H 2 −−→ GaH 3 + 3 CH 4 . The CH 4 intensity ex-hibits a similar temperature dependence for the case of a H 2 atmosphere containing NH 3 , which corresponds to the system on which we performed data assimilation (note that the CH 4 intensity decreases again at higher temperatures to react with NH 3 , but this behavior is out of the scope of our reduced reaction model) 26 . If we assume that the intensity and concentration are proportional to each other in the kinetic domain, the data conversion can be expressed as follows: c CH 4 (T set ) = 3c 0 TMG I(T set ) − I min I max − I min . (A1) Here, c CH 4 is the converted CH 4 concentration, c 0 TMG is the concentration of supplied Ga(CH 3 ) 3 , I is the CH 4 intensity, I max corresponds to the completely decomposed situation, and I min corresponds to the completely undecomposed situation.