# Surface chemistry models for GaAs epitaxial growth and hydride cracking
using reacting flow simulations

Malik Hassanaly  malik[dot]hassanaly[at]nrel[dot]gov Computational Science
Center, National Renewable Energy Laboratory  Hariswaran Sitaraman
Computational Science Center, National Renewable Energy Laboratory  Kevin L.
Schulte  Chemistry & Nanoscience Department, National Renewable Energy
Laboratory  Aaron J. Ptak  Chemistry & Nanoscience Department, National
Renewable Energy Laboratory  John Simon  Chemistry & Nanoscience Department,
National Renewable Energy Laboratory  Kevin Udwary  Kyma Technologies, Inc.
Jacob H. Leach  Kyma Technologies, Inc.  Heather Splawn  Kyma Technologies,
Inc.

###### Abstract

Hydride Vapor Phase Epitaxy (HVPE) is a promising technology that can aid in
the cost reduction of III-V materials and devices manufacturing, particularly
high-efficiency solar cells for space and terrestrial applications. However,
recent demonstrations of ultra fast growth rates ($\sim$ 500 $\mu$m/h) via
uncracked hydrides are not well described by present models for the growth.
Therefore, it is necessary to understand the kinetics of the growth process
and its coupling with transport phenomena, so as to enable fast and uniform
epitaxial growth. In this work, we derive a kinetic model using experimental
data and integrate it into a computational fluid dynamics simulation of an
HVPE growth reactor. We also modify an existing hydride cracking model that we
validate against numerical simulations and experimental data. We show that the
developed growth model and the improved cracking model are able to reproduce
experimental growth measurements of GaAs in an existing HVPE system.

###### pacs:

72.80.Ey,68.43.Mn,47.11.-j,81.15.Gh

††preprint: AIP/123-QED

##  I Introduction

The following article has been accepted by Journal of Applied Physics. After
it is published, it will be found at Journal of Applied Physics.

III-V solar cells have been used for decades in space-based applications King
_et al._ (2006); Hoheisel, Philipps, and Bett (2010) due to their high
specific power and radiation resistance. There is an increasing push to lower
the cost of these light, flexible, high-efficiency devices so that they can be
used in terrestrial applications, including in transportation, building
integration, and consumer electronics. Entry into these markets requires
driving down the cost of III-V material manufacturing, and the recent re-
emergence of hydride vapor phase epitaxy (HVPE) is one possible pathway to do
so Simon  _et al._ (2019). HVPE uses relatively low-cost precursors, uses them
efficiently, and recently showed growth rates for GaAs and GaInP over 500 and
200 $\mu m/h$ McClure  _et al._ (2020); Metaferia  _et al._ (2019)
respectively, which could have a significant effect on overall costs Horowitz
_et al._ (2018). It is the use of uncracked hydride precursors, in this case,
AsH3 and PH3, that enables these high growth rates using HVPE Schulte  _et
al._ (2018). The active hydrogen released from the hydride molecule at the
growth surface is theorized to scavenge chemisorbed Cl atoms that otherwise
act to inhibit growth Hollan and Durand (1979). This hydride-enhanced
mechanism significantly reduces the kinetic barrier to growth, allowing for
rates exceeding 500 monolayers/s at moderate growth temperatures ($\sim 900$
K) that are conducive to the formation of other alloys, e.g. GaInP, needed for
standard solar cell device structures. Previous work by DenBaars  _et al._
(1986) showed that AsH3 could spontaneously crack at quartz surfaces following
the global surface reaction

$\ce{AsH3->\frac{1}{4}As4+\frac{3}{2}H2},$ (1)

which they modeled with a one-step chemical mechanism. A method to inhibit
arsine cracking at the substrate is to increase the velocity of the carrier
gas, thereby limiting the residence time of arsine in the hot portion of the
reactor. However, larger velocities can adversely affect external gas sources
mixing and lead to a non-uniform growth at the substrate, requiring careful
design of the growth system. Furthermore, since surface growth is faster, the
coupling with gas transport is stronger in hydride-enhanced HVPE growth. In
this context, a simulation-based understanding that couples the chemical
kinetics of arsine cracking and GaAs deposition with transport phenomena is
paramount for the optimization of the overall manufacturing process.

Existing growth models for the deposition of GaAs in the literature only
consider deposition from arsenic vapor (As2 or As4) Shaw (1975); Kangawa  _et
al._ (2002), but not from uncracked AsH3. In this work, a one-step kinetic
model for the deposition of GaAs from uncracked arsine and gallium chloride is
developed and describes the following global reaction

$\ce{AsH3(g)+GaCl(g)\rightleftharpoons GaAs(s)+HCl(g)+H2(g)}.$ (2)

The main contributions of the paper are: a) the development of a kinetic model
for surface growth of uncracked arsine; b) uncertainty bounds on the growth
model parameters are provided using a Bayesian calibration approach; c) an
update of the arsine cracking model of DenBaars  _et al._ (1986); d) the
cracking and growth models are integrated into a computational fluid dynamics
(CFD) solver.

The surface growth model is developed and validated by modeling the reaction
kinetics and the fluid dynamics in a lab-scale reactor designed for a 2-inch
substrate, and for which experimental growth data is available. The rest of
the paper is organized as follows. Section II presents the lab-scale reactor
simulated and the dataset used to develop the kinetic model. The numerical
method is described in Sec. III and the calibration of the cracking model from
experimental data is shown in Sec. IV. In Sec. V, the finite rate kinetic
model for GaAs deposition is developed. Concluding remarks are provided in
Sec. VI.

##  II Configuration and experimental dataset

###  II.1 D-HVPE reactor

The configuration simulated in this work (see illustration in Fig. 1) is from
experimental studies by Schulte  _et al._ (2018) that used a Dynamic-Hydride
Vapor phase epitaxy (D-HVPE) reactor where growth rates with uncracked arsine
were also reported. These experiments were conducted with a Dynamic-Hydride
Vapor phase epitaxy (D-HVPE) reactor Young  _et al._ (2013) which addresses
manufacturing challenges posed by HVPE. Highly efficient thin-film solar cells
typically consist of several layers, each one made of a different III-V
compound Biefeld, Koleske, and Cederberg (2015). When the growth rate of the
solar cell is low enough and the reactor equilibrates quickly, the substrate
can be kept idle in the same chamber while different external gas are injected
over time. With a higher growth rate and large residence time of gases - such
as in HVPE deposition reactors - this strategy can lead to non-abrupt
interfaces between the layers, which is detrimental to solar cell efficiency.
In the D-HVPE, the aforementioned issue is avoided by physically moving the
substrate to a different chamber each time a new layer needs to be grown. In
order to prevent inter-chamber flow contamination, the chambers are separated
by vertical gas curtains created by a flow of H2 (see illustration in Fig. 1).

Figure 1: Schematic of the D-HVPE reactor. Solid arrows indicate flows, while
dashed arrows point to different elements of the reactor. An external furnace
maintains part of the reactor wall temperature at $1073.15$ K (area surrounded
by ) and the rest at $923.15$ K (area surrounded by ).

The D-HVPE reactor considered here has been experimentally Metaferia  _et al._
(2019); Schulte  _et al._ (2018, 2019) and numerically McClure  _et al._
(2020); Schulte  _et al._ (2019) shown to provide high growth rate and
efficiency, while potentially enabling in-line manufacturing.

The metal-halide is generated in-situ by injecting HCl into a boat/reservoir
of liquid Ga. The boat is externally heated at $1073.15$ K to drive the
kinetics of the HCl + Ga -> GaCl reaction. The boat conversion efficiency was
previously estimated to be 90% Schulte  _et al._ (2016). It outputs a mixture
of gaseous gallium chloride GaCl and residual HCl that were not converted. A
center port (called inlet A in Schulte  _et al._ (2018)) introduces a mixture
of AsH3 and H2 in the chamber. The reactants are then deposited on a substrate
externally heated at $923.15$ K. In Fig. 1, the external heating is
represented by the red and blue dashed boxes.

###  II.2 Dataset

The dataset used for validation of our model is that presented in Schulte  _et
al._ (2018) where growth rates at the center of the substrate were measured
under a variety of flow conditions. The volume flow rates at standard pressure
and temperature for all the cases considered are shown in Tab. 1. The
experiments were conducted at a pressure of $0.829$ atm and with external gas
injected at ambient temperature.

Table 1: Volume flow rates in standard centimeter cubes per minute (sccm) for each case simulates. “Top Chamber" denotes the top boundary of the chamber that contains the substrate. “Adjacent Chamber" denotes the top boundary of the chamber that does not contain the substrate. [sccm] | Center Port | Top Chamber | Adjacent Chamber | Curtains
---|---|---|---|---
Case |  H2 |  AsH3 |  H2 |  HCl |  GaCl |  H2 |  H2
Case 1 | 1000 | 44.5 | 8200 | 21 | 9 | 10000 | 2500
Case 2 | 1600 | 44.5 | 7600 | 21 | 9 | 10000 | 2500
Case 3 | 2200 | 44.5 | 7000 | 21 | 9 | 10000 | 2500
Case 4 | 2500 | 33.5 | 7500 | 20.8 | 7.2 | 10000 | 2500
Case 5 | 2500 | 33.5 | 7500 | 21 | 9 | 10000 | 2500
Case 6 | 2500 | 33.5 | 7500 | 21.2 | 10.8 | 10000 | 2500
Case 7 | 2500 | 33.5 | 7500 | 21.6 | 14.4 | 10000 | 2500

An important observation from the experiment was the near-zero growth of GaAs
when AsH3 was recirculated into the reactor (through inlet B in Fig. 1)
similar to GaCl rather than direct injection through the center port (inlet A
in Fig. 1). The larger residence time of AsH3 in this case led to increased
cracking at the reactor walls before reaching the substrate. This hypothesis
is further bolstered by the success of growth rate models that considered a
deposition pathway only via cracked arsine Schulte  _et al._ (2019). Under the
conditions where all the arsine has cracked and in the range of operating
conditions investigated, an equilibrium kinetic model Schulte  _et al._ (2016)
indeed predicts growth rates of at most $2$ $\mu$m/h (Fig. 2) which is at
least one order of magnitude lower than growth rates measured in Schulte  _et
al._ (2018). Here, the As4 pathway is inhibited by the low residence time of
external gas sources and the low substrate temperature. Therefore, the growth
rate available in Schulte  _et al._ (2018) is particularly suited to deduce a
kinetic model for the AsH3-pathway for epitaxial growth.

Figure 2: GaAs(s) growth rate via As4 reaction pathway at $923.15$ K for
various partial pressure of GaCl.

##  III Numerical method

To derive a kinetic model for the surface reaction via uncracked hydrides,
computational fluid dynamics (CFD) simulations of the experimental
configuration are conducted. The kinetic model depends on parameters that are
then calibrated to match experimental growth rates. The fluid solver used is
described in App. A, the implementation of surface reactions is explained in
App. B and the transport and thermodynamic parameters are given in App. C.

The boat injector is omitted from the domain to simplify the geometry.
Instead, the boat injector is modeled using an injector spread over the top of
the chamber aligned with the substrate, assuming uniform composition of GaCl,
HCl and H2. The volume flow rates for the top chamber are available in Tab. 1.
Additionally, to reduce the range of scales resolved, the inner region of the
center port is not meshed and its exit is treated as an inlet boundary
condition. Because AsH3 is mostly exposed to walls inside the center port,
this is where it mostly cracks. The boundary conditions at the exit of the
center port need specific treatment to incorporate this effect and are
described in App. D.

A hexahedral-dominant mesh shown in Fig. 3 is used with refinement near the
injector exit and the platter, where large gradients of velocity and
composition are observed. To ensure that the results presented are not
influenced by numerical errors, two different meshes are used by increasing
the resolution in the chamber that contains the substrate. The coarse mesh
contains about $0.5$ million computational cells and the finer mesh consists
of $1$ million cells. The smallest mesh size in all three directions is $0.75$
mm and $0.47$ mm for the coarse and fine meshes, respectively.

Figure 3: Illustration of the mesh used. The blue parts denote the center port
and the substrate. Note the refinement near the center port exit and near
substrate, where gradients are the largest.

##  IV Arsine cracking model

As mentioned in Sec. I, the cracking of hydrides (here AsH3) critically
impacts the growth rate of GaAs layers. Therefore, a cracking kinetic model
needs to be included with other surface reactions and coupled with the
transport equations to adequately predict deposition rates. A global kinetic
model for the decomposition of hydrides at quartz surfaces (the walls of the
D-HVPE reactor considered here are also made of quartz) of AsH3 diluted within
a H2 carrier gas was developed in DenBaars  _et al._ (1986). One of the
experiments reported in DenBaars  _et al._ (1986) included measurement of AsH3
concentration after it has passed through a $5$ cm diameter and $12.5$ cm long
heated quartz tube at a speed of $1.25$ cm/s. The kinetic rate of the global
reaction (Eq. 1) was modeled in the Arrhenius form as

$K_{cracking}=A_{cracking}exp\left(-\frac{Ea_{cracking}}{R_{u}T}\right),$ (3)

where $A_{cracking}$ and $Ea_{cracking}$ are constant kinetic coefficients for
pre-exponential factor and activation energy, $R_{u}$ is the universal gas
constant, and $T$ is the surface temperature. The values reported for
$K_{cracking}$ by DenBaars  _et al._ (1986) were computed by assuming with
respect to volumetric arsine concentration, given by

$K_{cracking}=\frac{1}{\Delta t}log\left(\frac{C_{0}}{C_{f}}\right),$ (4)

where $\Delta t$ is the residence time of the gas in the quartz tube, $C_{0}$
and $C_{f}$ are the volumetric molar concentrations of AsH3 at the inlet and
outlet of a tube through which a mixture of arsine and hydrogen flowed. While
the expression of $\Delta t$ is not explicitly given, the value of
$K_{cracking}$ provided in DenBaars  _et al._ (1986) suggests that it was
calculated as $\frac{L}{U}$, where $L$ is the length of the tube and $U$ is
the streamwise velocity. The expression chosen for $K_{cracking}$ neglected
the effect of surface-to-volume ratio which depends on the geometry of the
tube. For instance, if the tube had a smaller diameter, and under the same
conditions as the experiments – in particular same volume flow rates –, the
apparent $K_{cracking}$ would be larger than the one seen in the experiments.
The importance of including geometry considerations in the kinetic constants
was also recognized elsewhere Harrous  _et al._ (1988). Furthermore, since
$K_{cracking}$ is the rate constant for a surface reaction it should have
units of m.s-1. DenBaars  _et al._ (1986) factored in the area-to-volume ratio
(units of $m^{-1}$) within $K_{cracking}$ that gave rise to a first-order rate
constant with units of $1/s$. The authors reported a value for $Ea_{cracking}$
= 34$~{}kcal/mol$ by varying the temperature of the tube but did not propose a
value for $A_{cracking}$. One of the main contributions of this work is the
calibration of $A_{cracking}$ along with a correction to the expression for
$K_{cracking}$ that accounts for the surface-to-volume ratio.

One can estimate the cracking of AsH3 in a cylindrical tube using a steady one
dimensional advection-reaction equation (justification is provided in App. D).
The mass fraction of arsine $Y_{\ce{AsH3}}$ along the streamwise direction
follows the relation

$\pi R^{2}\frac{d}{dx}\left(\rho UY_{\ce{AsH3}}\right)=-K_{cracking}\rho
Y_{\ce{AsH3}}2\pi R,$ (5)

where $U$ is the streamwise velocity and $R$ is the tube radius. Assuming that
the density $\rho$ varies only due to temperature variations, one can derive
an approximate relation for the ratio $C_{f}/C_{0}$ as

$K_{cracking}=\frac{UR}{2L}log\left(\frac{C_{0}}{C_{f}}\right).$ (6)

The new pre-log factor can be thought of as a residence timescale rescaled
with a non-dimensional ratio of catalytic surface per unit volume. With that
rescaling, the cracking kinetic coefficient now includes the effect of tube
radius and the expression can be directly used to obtain a value for
$A_{cracking}$ by fitting the experimental data to the 1D model.

To incorporate modeling and experimental uncertainty in the calibration of
$A_{cracking}$, a Bayesian calibration approach is adopted Braman, Oliver, and
Raman (2013); Bell  _et al._ (2019). The output of that procedure is a
posterior probability density function (PDF) of $A_{cracking}$ given the
experimental data available and is computed as

$P(A_{cracking}|D)\sim P(D|A_{cracking})P(A_{cracking}),$ (7)

where $P(A_{cracking}|D)$ is the posterior PDF, $P(A_{cracking})$ is a prior
PDF (describes the knowledge of $A_{cracking}$ without any experimental
measurement) and $P(D|A_{cracking})$ is the likelihood of the data. An
uninformative prior $P(A_{cracking})\sim\mathcal{U}(10^{4},2\times 10^{7})$ is
chosen. Following Braman, Oliver, and Raman (2013), a Gaussian likelihood is
adopted,

$P(D|A_{cracking})\sim\frac{1}{(2\pi\sigma)^{N_{d}/2}}exp\left[-\frac{1}{2\sigma^{2}}\sum_{i=1}^{N_{d}}(x_{i}(A_{cracking})-D_{i})^{2}\right],$
(8)

where $N_{d}$ is the number of data points, $x_{i}$ is the model prediction
for the $i^{th}$ experiment, $D_{i}$ is the $i^{th}$ experimental value, and
$\sigma$ is a hyperparameter which describes the experimental errors. Since
$\sigma$ is unknown here, it is fit as part of the Bayesian calibration
procedure with prior $\mathcal{U}(0,10)$, along with $A_{cracking}$ Braman,
Oliver, and Raman (2013). It was found that the choice of prior had little
influence on the calibrated parameters (shown in App. E. The results of the
Bayesian calibration are shown in Fig. 4.

Figure 4: Top: experimental measurement ( ) of AsH3 cracking at different
temperatures, 1-step cracking model prediction with $A_{cracking}=7.5\times
10^{5}$ m.s-1 ( ) corresponding to the largest posterior probability and
uncertainty bounds $\pm\sigma$ ( ). Bottom: corner plot of the 2-dimensional
posterior PDF.

##  V Gallium arsenide growth model

###  V.1 Model functional form

The surface reaction is modeled as an irreversible forward reaction

$\ce{AsH3(g)+GaCl(g)->GaAs(s)+HCl(g)+H2(g)}.$ (9)

The assumption of irreversibility is justified by its large equilibrium
constant which suggests that the backward reaction is unlikely McClure  _et
al._ (2020); Grüter  _et al._ (1989). A finite rate chemistry model is
developed hereafter for the forward reaction. An infinitely fast model is also
constructed and the results are reported in App. F.

The finite rate chemistry model for the reaction at the substrate is
constructed by assuming the Arrhenius form
$K_{reac}=A_{reac}T^{\beta_{reac}}exp(-Ea_{reac}/R_{u}T)$, where $A_{reac}$,
$\beta_{reac}$ are kinetic constants and $Ea_{reac}$ is the activation energy.
In the experiments conducted in Schulte  _et al._ (2018), it was found that
for the temperature range $[873-953~{}K]$, the dependence of the epitaxial
growth rate with temperature was negligible. Therefore, $K_{reac}$ is assumed
to be independent of temperature over similar temperature ranges studied in
this work.

The objective of this section is to find a value of $K_{reac}$. The CFD model
is used to simulate the D-HVPE reactor with different values of $K_{reac}$.
The resulting growth rates are obtained by computing the molar flux of AsH3 at
the substrate, which is equal to the number of moles of GaAs produced. The
density of GaAs(s) is then used to compute the growth rate obtained via CFD,
which is then compared to the ones reported in the Schulte  _et al._ (2018).
The growth rates are measured at the center of the substrate.

###  V.2 Model calibration

Eleven values of $K_{reac}$ are chosen within the range [$45$ m.s-1,$310$
m.s-1]. For each value of $K_{reac}$, the seven cases shown in Tab. 1 are
simulated. In absence of uncertainty, a total of $77$ numerical simulations
would be needed for calibration. In the present case, the optimal value of
$K_{reac}$ depends on the value of $A_{cracking}$ obtained by fitting the
cracking model. To account for uncertainty in the value of $A_{cracking}$, the
calibration procedure is conducted for ten values of $A_{cracking}$ chosen in
the range [$6.5\times 10^{5}$ m.s-1,$8.5\times 10^{5}$ m.s-1], which span the
support of its posterior distribution. The function that related $K_{reac}$
and $A_{cracking}$ is then convolved with the posterior distribution of
$A_{cracking}$ to deduce the PDF of $K_{reac}$. The overall procedure is
illustrated in Fig. 5 and requires a total of $770$ runs.

Figure 5: Illustration of the workflow of the calibration procedure.

Due to the large number of simulations required, the calibration is done with
the coarse grid ($0.5$ million computational cells). The fine mesh ($1$
million computational cells) is used afterwards to verify grid convergence for
the estimated values of $K_{reac}$.

###  V.3 Results

The steady-state distribution of temperature and arsine mass fraction for case
3 (Tab. 1) with $K_{reac}=222.8$ m.s-1 and $A_{cracking}=8.5\times 10^{5}$
m.s-1 is illustrated in Fig. 6. It can be seen that since the flow inside the
center port is heated to a temperature close to $1073.15$ K, the center port
injects a hot flow at the substrate (indicated by the arrow). The contour of
AsH3 (bottom) shows the effect of the flow curtain which avoids contamination
of the adjacent chamber with arsine (left arrow). At the substrate, AsH3 is
consumed, which results in a non-zero normal gradient (right arrow). Similar
features can be observed for GaCl.

Figure 6: Top: steady-state contour of temperature in the D-HVPE reactor for
Case 3. Middle: steady-state contour of mass fraction of AsH3 in the D-HVPE
reactor for Case 3. Bottom: steady-state contour of mass fraction of GaCl in
the D-HVPE reactor for Case 3. Arrows point to regions of interest.

For each $(K_{reac},A_{cracking})$ sample pair, the average error in growth
rate $\langle\varepsilon\rangle$ over the seven cases considered is computed.
The results are shown in Fig. 7. It can be seen that there exists an optimal
value of $K_{reac}$ which leads to the best match of the growth rates. In the
particular case considered ($A_{cracking}=8.5\times 10^{5}$ m.s-1),
$K_{reac}=222.8$ m.s-1 minimizes $\langle\varepsilon\rangle$. To improve the
estimate of the value of $K_{reac}$ that minimizes the error, interpolation
between data points is done with radial basis functions, where hyperparameters
are optimized using the Scikit-Learn library Pedregosa  _et al._ (2011). An
example of the interpolation result is shown in Fig. 7.

Figure 7: Average growth rate prediction error across the seven cases
investigated as a function of the kinetic constant $K_{reac}$ with
$A_{cracking}=8.5\times 10^{5}$ m.s-1. Simulation results ( ) and interpolated
results with radial basis functions ( ) are displayed.

Over the range chosen $K_{reac}$, it can be observed from Fig. 7 that there
exists a value that can serve, at least, as a global minimizer of the
discrepancy with experimental results. It is shown here, that the local
minimum is likely a global minimum. We first argue that the relationship
between the rate of growth $R_{G}$ and the kinetic coefficient $K_{reac}$ is
monotonic and increasing. Physically, this assertion is justified by the fact
that increasing the kinetic rate should increase the rate of production of
GaAs. We also propose a numerical argument by plotting the difference between
the rate of growth obtained in the simulations ($R_{G,sim}$) and in
experiments ($R_{G,exp}$) for all $A_{cracking}$ and all $K_{reac}$ values.
The results are shown in Fig. 8 (top). As expected, the rate of growth is
monotonically increasing with $K_{reac}$ for all $A_{cracking}$ values. Next,
we record the minimum and maximum error over all the cases for each
$A_{cracking}$ and $K_{reac}$. Because of the monotonically increasing
relation between $K_{reac}$ and $R_{G}$, if the maximum error
$R_{G,sim}-R_{G,exp}$ at the lowest $K_{reac}$ value is negative, then
$\forall K_{reac}<45$ m.s-1, $\langle\varepsilon\rangle$ would be larger than
for $K_{reac}=45$ m.s-1 . In Fig. 8 (right), it can be seen that the maximum
absolute value of the error is almost negative for all $A_{cracking}$ at
$K_{reac}=45$ m.s-1. The same argument can be invoked for the minimum error at
$K_{reac}=310$ m.s-1 which is positive. We therefore argue that the minimum
shown is likely a global minimum.

Figure 8: Left: average difference between the rate of growth obtained from
simulations $R_{G,sim}$ and the rate of growth obtained from experiments
$R_{G,exp}$ plotted against $K_{reac}$ for all $A_{cracking}$. Right: minimum
( ) and maximum ( ) difference of the rate of growth obtained from simulations
$R_{G,sim}$ compared to rate of growth obtained from experiments $R_{G,exp}$
plotted against $K_{reac}$ for all $A_{cracking}$.

The growth rates computed with $K_{reac}=222.8$ m.s-1 and
$A_{cracking}=8.5\times 10^{5}$ m.s-1 are plotted against experimental values
in Fig. 9. A reasonable agreement with experiments can be observed.
Furthermore, the trends are consistent with the experiments and physical
intuition. As the volume flow rate of H2, $Q_{a}$ through the center port
increases (left), less AsH3 cracks which leads to larger reaction rates. The
exponential dependence of the amount of cracked AsH3 with the velocity (Eq.
17) should result in a super-linear increase of growth rate with the flow rate
of arsine through the center port. This is the trend observed in the
experiments and the simulations. This observation is in sharp contrast with
the linear increase of the growth rate with increasing concentration of GaCl
(bottom). To ensure that the calibrated $K_{reac}$ is only minimally affected
by the numerical errors, the fine grid (1 million computational cells) is used
to verify that the results are grid converged. The fine and coarse grids
predict similar growth rates and therefore the baseline grid is sufficient for
accurate predictions.

Figure 9: Growth rate of GaAs comparison between experimental data reported in
Schulte  _et al._ (2018) ( ), the calibrated model used with the 0.5 million
cells grid ( ) and the calibrated model used with the 1 million cells grid ( )
for $A_{cracking}=8.5\times 10^{5}$ m.s-1. Results are overlayed with
calibrated model obtained with the 0.5 million cells grid for all other
$A_{cracking}$ values ( ). Top: Case 1-3 where growth rate is plotted against
$Q_{a}$, the volume flow rate of H2 through the center port. Bottom: Case 4-7
where growth rate is plotted against $P_{GaCl}$, the partial pressure of GaCl.

The spatial distribution of the growth rate at the substrate is compared to
experimental data in Fig. 10. A reasonable agreement can be observed for case
3 (strong flow through the center port). A breakdown of the axial symmetry can
be observed for case 2 (weaker flow through the center port) both in
simulation and experimental data. The apparent breakdown of symmetry typically
appears at low center port velocity, such as the one of Case 2. The flow near
the substrate is influenced by a vertical component that originates from the
center port, and by a horizontal component, the crossflow that connects the
two deposition chambers. The horizontal flow can be observed in Fig. 6 (middle
and bottom). When lower center port velocities are used, the relative effect
of the horizontal crossflow may increase and disturb the boundary layer near
the substrate, thereby leading to a breakdown of growth symmetry.
Interestingly, the same symmetry breakdown seems to occur in the experimental
results, albeit to a smaller extent. The discrepancy may be due to a slightly
earlier center port jet breakdown in the simulations due to either numerical
diffusion, or errors in the center jet exit boundary conditions. Finally, we
note that in other experiments reported for the same reactor Schulte  _et al._
(2019), a similar breakdown of symmetry may appear at different conditions.
This result suggests that the model can be used not only to assess the overall
growth rate of the substrate but also its spatial uniformity, which is
critical for the efficiency of solar cells.

Figure 10: Spatial distribution of growth rate at the substrate for case 2
(left) and case 3 (right) between experimental data (top) and simulation data
(bottom) with $A_{cracking}=8.5\times 10^{5}$ m.s-1 and $K_{reac}=222.8$
m.s-1.

As explained in Sec. V.2, the calibration procedure is repeated for ten
different values of $A_{cracking}$ which span the support of its posterior
distribution obtained in Sec. IV. As a result, one obtains the response curve
of the optimal $K_{reac}$ as a function of $A_{cracking}$, which is shown in
Fig. 11 (top). As $A_{cracking}$ increases, the amount of AsH3 that reaches
the substrate decreases. Therefore, to keep constant the growth rates, the
optimal $K_{reac}$ is a monotonically increasing function of $A_{cracking}$.
This trend can be observed in Fig. 11. It can however be seen that the
increase occurs by steps which is a consequence of the discretization of
$A_{cracking}$. Given the dependence of $K_{reac}$ with $A_{cracking}$, one
can sample $A_{cracking}$ using a Markov-chain Monte Carlo (MCMC) approach and
construct a PDF of $K_{reac}$. The results are shown in Fig. 11 (bottom). The
standard deviation of the PDF of $K_{reac}$ can then be used to extract
uncertainty bounds on the calibrated value of the kinetic parameter. We found
$K_{reac}=205\pm 5$ m.s-1

Figure 11: Top: calibrated value of $K_{reac}$ for each value of
$A_{cracking}$ considered. Bottom: constructed PDF of $K_{reac}$.

##  VI Conclusions

In this work, a comprehensive modeling procedure was developed for the
simulation of growth rate via uncracked hydrides in a D-HVPE reactor. The AsH3
cracking model was adjusted for the geometry and a surrogate model for the
cracking through the injectors was developed. Multiple simulations for the
surface reactions at the substrate were used to calibrate a simple kinetic
model for the reaction of uncracked arsine. The proposed model matches
experimental growth rate values and spatial distributions for a range of
typical HVPE operating conditions.

As part of the calibration, an uncertainty quantification procedure was
conducted to include uncertainty in the AsH3 cracking model parameters. While
the influence of other uncertain parameters like numerical CFD errors was
shown to be negligible, uncertainty due to the interpolation during the
calibration procedure will be included as future work.

Since III-V solar cells typically require the growth of different layers, the
same procedure will be applied to develop kinetic models that correspond to
each one of the compounds. Finally, the validated kinetic model will be used
in conjunction with a CFD software to design scaled-up reactors able to
maximize growth rate while minimizing growth non-uniformity.

###### Acknowledgements.

Fruitful discussions with Milo Parra-Alvarez and Olga Doronina are gratefully
acknowledged. This work was authored by the National Renewable Energy
Laboratory, operated by Alliance for Sustainable Energy, LLC, for the U.S.
Department of Energy (DOE) under Contract No. DE-AC36-08GO28308. This work was
supported by the Air Force Research Laboratory (IAG-19-02103). This research
was performed using computational resources sponsored by the U.S. Department
of Energy’s Office of Energy Efficiency and Renewable Energy and located at
the National Renewable Energy Laboratory. The views expressed in the article
do not necessarily represent the views of the DOE or the U.S. Government. The
U.S. Government retains and the publisher, by accepting the article for
publication, acknowledges that the U.S. Government retains a nonexclusive,
paid-up, irrevocable, worldwide license to publish or reproduce the published
form of this work or allow others to do so, for U.S. Government purposes.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A Fluid solver

A variable density low-Mach number solver based on the OpenFOAM framework is
used Hassanaly  _et al._ (2018) to model the flow inside the reactor. The
solver has been successfully applied in several other studies Hassanaly  _et
al._ (2020); Tang and Raman (2021) including configuration with reactive phase
transformations such as the one considered here Koo  _et al._ (2017). The low-
Mach number assumption allows for the use larger time-steps on the order of
fluid convection time-scales($\sim$ $1$ ms) in all the simulations, without
the need to resolve acoustic time-scales. Thermodynamics variables (specific
heat, enthalpy) and transport properties (viscosity, thermal conductivity,
species diffusivity) used in the model are obtained from open source library,
PelePhysics, designed for complex combustion simulations Sitaraman  _et al._
(2021).

The momentum transport equation given by:

$\frac{\partial\rho\bm{u}}{\partial t}+\nabla\cdot(\rho\bm{u}\bm{u})=-\nabla
p+\nabla\cdot\bm{\overline{\tau}},$ (10)

is solved, where the pressure $p$ is computed to ensure mass conservation and
$\bm{\overline{\tau}}$ is the stress tensor.

The gas-phase species are transported using the equation

$\frac{\partial\rho Y_{i}}{\partial
t}+\nabla\cdot(\rho\bm{u}Y_{i})=\nabla\cdot(\rho D_{i}\nabla
Y_{i})+\dot{\omega_{i}},$ (11)

where
$Y_{i}~{}\in~{}\\{Y_{\ce{AsH3}},Y_{\ce{As4}},Y_{\ce{GaCl}},Y_{\ce{HCl}}\\}$
are species mass fractions, $\rho$ is the gas mixture density, $\bm{u}$ is the
flow velocity and $D_{i}$ are the individual species diffusivities. The mass
fraction of the carrier specie $Y_{\ce{H2}}$ is obtained by enforcing the
unity mass fraction sum constraint.

The energy equation is solved in terms of a global temperature among all
species given by:

$\frac{\partial\rho T}{\partial
t}+\nabla\cdot(\rho\bm{u}T)=\nabla\cdot(\lambda\nabla T)-\rho(\nabla
T)\cdot\sum_{i=1}^{N}Cp_{i}Y_{i}\bm{V_{i}}+\dot{\omega_{HR}},$ (12)

where $T$ denotes temperature, $\lambda$ is the thermal conductivity, $Cp_{i}$
are the species heat capacities, $\bm{V_{i}}$ are the species diffusion
velocities, $N$ is the number of species, and $\dot{\omega_{HR}}$ denotes the
heat release rate.

A Smagorinsky Smagorinsky (1963) turbulence model is used and is mostly active
in the shear layer near the center-port exit, while the flow near the
substrate is nearly laminar. For the scalar turbulence, a turbulent Schmidt
number $Sc_{t}=0.72$ and a turbulent Prandtl number $Pr=0.7$ were used Tang
_et al._ (2019). The turbulence model is useful for stabilization purposes in
the shear layers and all the results in this work are obtained when the
reactor reaches steady-state.

##  Appendix B Numerical implementation of finite-rate surface reactions

The surface reactions are implemented in a fashion similar to Maestri and
Cuoci (2013) by imposing a non-zero flux of species at the walls. Formally,
the reaction source term for species $i$ is written as

$\dot{\omega_{i}}=\nabla\cdot\Omega_{i},$ (13)

where $\Omega_{i}$ is a field defined at the cell faces which is null
everywhere in the domain, except at the walls where surface reactions occur.
$\Omega_{i}$ is computed using the kinetic model and by assuming the surface
concentration of the species equals that of the volumetric concentration at
the computational cell immediately adjacent to the face. The same procedure is
used to define the heat release from the walls. The surface reactions are
coupled with the scalar transport equation using a Strang splitting method
Strang (1968), except when infinitely fast chemistry is used (See App. F).

The surface reactions in this work involve a conversion of species in the gas
phase to solid phase products at the catalytic surfaces. Therefore, mass may
leave or enter the gas phase at these surfaces. Similar to other cases where a
phase change occurs like with spray evaporation Heye, Raman, and Masri (2013)
or soot formation Chong  _et al._ (2018), the pressure equation needs to be
augmented with a source term for mass. The mass source term $\dot{m}$ is
simply computed as

$\dot{m}=\sum_{i=1}^{N}\dot{\omega_{i}},$ (14)

where $N$ is the number of species.

##  Appendix C Thermodynamic and transport parameters

The thermodynamic parameters (enthalpy and calorific capacity) are obtained
from polynomial fits that relate them to the local temperature. The
coefficients of the polynomial fits are obtained from Oehlschlaeger  _et al._
(2009). Missing properties of GaCl were substituted with the properties of
GaH. Note that although the polynomial coefficient for the calorific capacity
and enthalpy were substituted, the molar weight (used to compute density) was
not modified. Similarly, the transport parameters (species diffusivity, heat
conductivity, dynamic viscosity) are computed at every grid point using the
species concentrations and temperature. The transport parameters are computed
using the following molecular properties: their Lennard Jones potential well
depth, collision diameter, dipole moment, polarizability, and their rotational
relaxation collision number. The molecular properties are obtained from Kee
_et al._ (1999). Again, since the properties of some species were not
available, they had to be substituted with other ones. In particular, the
molecular properties of GaCl and As4 were substituted with that of GaH and
As2.

##  Appendix D Cracking model validation and center port surrogate model

The two assumptions used to derive the cracking model namely the 1D advection-
reaction approximation and the sole dependence of density on temperature are
addressed. The model is compared to three-dimensional reacting flow
simulations. The validation of these assumptions will further bolster the
correction applied to the cracking kinetic coefficients.

The 1D transport of enthalpy in a cylindrical tube neglecting axial diffusion
can be approximated as

$\frac{d}{dx}(\rho UCpT)\pi R^{2}dx=h(T_{w}(x)-T(x))2\pi Rdx,$ (15)

where $\rho$ is the local density, $Cp$ is the local heat capacity, $T(x)$ is
the local temperature, $T_{w}(x)$ is the wall temperature and $R$ is the tube
radius and $h$ is the convective heat transfer coefficient. The convective
heat transfer coefficient for a fully developed pipe flow can be obtained from
Nusselt number $Nu=3.66$ Bergman  _et al._ (2011). The heat transfer
coefficient can be written as $h\approx Nu\lambda/(R)$ where $\lambda$ is the
thermal conductivity. Assuming a near-constant heat capacity and neglecting
the addition of mass due to surface reactions

$\frac{dT}{dx}=\frac{2h}{R\rho_{0}U_{0}Cp_{0}}(T_{w}(x)-T(x)),$ (16)

where $\rho_{0}$, $U_{0}$, and $Cp_{0}$ are the fluid properties at the pipe
inlet. The same approach for the transport of mass fraction of AsH3 leads to

$\frac{dY_{\ce{AsH3}}}{dx}=\frac{-2Y_{\ce{AsH3}}}{R\rho_{0}U_{0}}A_{cracking}\,\,e^{\frac{-Ea_{cracking}}{R_{u}T_{w}(x)}}\frac{\rho_{0}T_{0}}{T(x)},$
(17)

where $\rho(x)$ was assumed to be solely a function of temperature. Given a 1D
discretization, Eqs. 16 and 17 can be solved with an exponential type of
integrator such as a Rosenbrock integrator.

In Fig. 12 the results for a 3D reacting flow simulation are shown for
conditions similar to the ones of DenBaars  _et al._ (1986). The wall
temperature is held fixed at $723.15$ K and $773.15$ K, the cylinder radius is
$25$ mm, its length is $125$ mm, the pressure is $1$ atm, and the inlet is
made of 50% AsH3 and H2 by volume. The mesh size is $1.5$ mm radially and $4$
mm axially.

Figure 12: Left: contour of mass fraction of AsH3 in a cylinder that
replicates the experiments DenBaars  _et al._ (1986). Flow is from top to
bottom. Right: surface averaged longitudinal profile of mass fraction of AsH3
for a wall temperature of $723.15$ K for the 3D model ( ) and the 1D model ( )
and for a wall temperature of $773.15$ K for the 3D model ( ) and the 1D model
( ).

The profile of $Y_{\ce{AsH3}}$ exhibits variation within a cross-section due
to the relatively larger diameter. Nevertheless, the 1D model successfully
approximates the longitudinal profile of $Y_{\ce{AsH3}}$. In turn, it suggests
that the correction for the kinetic rate for arsine cracking is valid.

The same procedure can be applied to the center port of the HVPE reactor
studied in this work. In this case, the cylinder has a $2$ mm radius and its
length is $308$ mm. The center port is even more suited to the 1D
approximation given its larger length to radius ratio. The wall temperature is
held fixed at $1023.15$ K, and the pressure is $0.829$ atm. The mesh size is
$0.12$ mm radially and $4$ mm axially. The results are shown in Fig. 13. Since
the aspect ratio of the cylinder is smaller than that of DenBaars  _et al._
(1986), the profile of AsH3 is almost one dimensional (left). In turn, the 1D
model accurately reproduces the longitudinal profile of AsH3.

Figure 13: Left: contour of mass fraction of AsH3 in a cylinder that
replicates the center port. The longitudinal direction was rescaled by a
factor 0.032 for plotting purposes. Flow is from top to bottom. Right: surface
averaged longitudinal profile of mass fraction of AsH3 for an outer-wall
temperature of $1023.15$ K without taking into account the heat transfer
through the quartz for the 3D model ( ) and the 1D model ( ), and by taking
into account the heat transfer through the quartz for the 3D model ( ) and the
1D model ( ).

In the experiments reported in Schulte  _et al._ (2018), the center port has
an inner radius $R_{in}=2$ mm and an outer radius of $R_{out}=4.5$ mm. One can
also incorporate the conjugate heat transfer between the quartz annulus and
the inner flow. The inner-wall temperature can be obtained by writing a
balance of heat flux at the inner-wall boundary

$h(T(x)-T_{iw}(x))+\frac{\lambda_{quartz}}{\Delta R}(T_{w}(x)-T_{iw}(x))=0,$
(18)

where $\lambda_{quartz}$ is the heat conductivity of the quartz, $T_{iw}$ is
the inner-wall temperature, $T_{w}$ is the outer-wall temperature and $\Delta
R=R_{out}-R_{in}$. An algebraic expression can be obtained for $T_{iw}$, and
Eq. 16 becomes

$\frac{dT}{dx}=\frac{2h}{R\rho_{0}U_{0}Cp_{0}\left(1+\frac{h\Delta
R}{\lambda_{quartz}}\right)}(T_{w}(x)-T(x)),$ (19)

and Eq. 17 becomes

$\frac{dY_{\ce{AsH3}}}{dx}=\frac{-2Y_{\ce{AsH3}}}{R\rho_{0}U_{0}}A_{cracking}e^{-Ea_{cracking}/(R_{u}T_{iw}(x))}\frac{\rho_{0}T_{0}}{T(x)},$
(20)

where $T_{iw}(x)=T(x)+\frac{T_{w}(x)-T(x)}{1+h\Delta R/\lambda_{quartz}}$.

The 1D model is compared to a 3D model run with CFD, where the heat transfer
through the solid is written assuming constant thermodynamic and transport
properties

$\frac{\partial}{\partial
t}(\rho_{quartz}T)=\frac{1}{Cp_{quartz}}\nabla\cdot(\lambda_{quartz}\nabla
T),$ (21)

where $\rho_{quartz}$ is the density of the solid assumed to be equal to
$2200$ kg.m-3, $Cp_{quartz}$ is the heat capacity of the solid set at $705$
J.kg-1.K-1, and $\lambda_{quartz}$ is equal to $1.46$ W.m-1.K-1. The coupling
between the solid and fluid equations is performed via outer-Picard iterations
Hassanaly  _et al._ (2018). The thermal coupling at every face of the inner-
wall boundary is done by dynamically setting the inner-wall temperature to

$T_{iw}=\frac{T_{s}\frac{\lambda_{s}}{\Delta
s}+T_{f}\frac{\lambda_{fluid}}{\Delta f}}{\frac{\lambda_{s}}{\Delta
s}+\frac{\lambda_{fluid}}{\Delta f}},$ (22)

where $T_{s}$ (respectively $T_{f}$) is the temperature in the solid
(respectively fluid) domain that is radially adjacent to the inner-wall face,
$\lambda_{s}$ (respectively $\lambda_{f}$) is the heat transfer coefficient at
that same location, $\Delta s$ (respectively $\Delta_{f}$) is the distance
between that location and the inner-wall face center. The choice of the inner-
wall temperature boundary condition equates both the temperature and the heat
flux at the inner-wall interface. The mesh size is the same as the one used
without the conjugate heat transfer, and the results are shown in Fig. 13
(right) for an outer-wall temperature of $1023.15$ K. The 1D model and the 3D
model reasonably agree on the longitudinal profile of AsH3. Because the inner-
wall temperature is lower than without the conjugate heat transfer, less AsH3
cracks through the center port. Similar to the results shown in Fig. 12, since
the 1D model neglects radial diffusion, it overestimates the amount of AsH3
immediately adjacent to the walls, thereby slightly overestimating the amount
of AsH3 cracked.

Overall, the 1D model is a reasonable representation at lower computational
costs for modeling transport and chemistry within the center port and it is
proposed to be used as a surrogate for the flow inside the center port to
obtain the boundary conditions at the exit of the injector. The boundary
conditions for velocity can be obtained from mass conservation assuming that
it solely depends on the density at the exit of the tube which, in turn, is
assumed to be a sole function of temperature. The mass fraction of As4 at the
center port exit can be computed by integrating Eq. 20 over the length of the
tube.

In the CFD model, the surrogate model is implemented by first collecting the
value of temperature (transported according to Eq. 12) at the outer-walls of
the center port, at every timestep. The temperature values are averaged
azimuthally to construct a one-dimensional representation of the outer-wall
temperature of the center port. The one-dimensional wall temperature profile
is then directly used in the advection-reaction model of the center port. The
boundary conditions at the injector exit are subsequently evaluated by
integrating the 1D ordinary differential equation systems (Eqs. 19 and 20) for
the scalars and temperature. This procedure allows incorporating longitudinal
gradients of temperature along the center-port walls into the center-port exit
boundary condition. Although the method is repeated at every timestep, the
overhead computational cost is negligible compared to momentum and scalar
transport.

##  Appendix E Effect of prior on the Bayesian calibration

In Bayesian analysis, the posterior distribution depends on the prior, which
in turn may have an effect on the PDF of the calibrated $K_{reac}$. In the
results reported in the paper, the prior was chosen to be
$\mathcal{U}(10^{4},2\times 10^{7})$ for $A_{cracking}$ and
$\mathcal{U}(0,10)$ for $\sigma$ . In Fig. 14, the results are reported with a
narrow prior $\mathcal{U}(10^{5},10^{6})$ for $A_{cracking}$ and
$\mathcal{U}(0,1)$ for $\sigma$; and a wide prior $\mathcal{U}(1,2\times
10^{11})$ for $A_{cracking}$ and $\mathcal{U}(0,100)$ for $\sigma$. As can be
observed, the prior has little effect on the posterior distribution of
$A_{cracking}$.

Figure 14: Black line: Two-dimensional posterior PDF for prior used
($\mathcal{U}(10^{4},2\times 10^{7})$ for $A_{cracking}$ and
$\mathcal{U}(0,10)$ for $\sigma$). Blue line: Two-dimensional posterior PDF
for narrow prior ($\mathcal{U}(10^{5},10^{6})$ for $A_{cracking}$ and
$\mathcal{U}(0,1)$ for $\sigma$). Red line: Two-dimensional posterior PDF for
wide prior ($\mathcal{U}(1,2\times 10^{11})$ for $A_{cracking}$ and
$\mathcal{U}(0,100)$ for $\sigma$).

##  Appendix F Infinitely fast kinetic growth model

An infinitely fast kinetic model is implemented following the hypothesis that
the reactor operates in a transport-limited regime. A mass flux boundary
condition is imposed for the species transport equations at the substrate
using the reactant composition at adjacent computational cells. Molar
stoichiometry is used to determine excess and deficient reactants at the
surface and surface gradients are subsequently set to match the desired
consumption rate.

The growth rate results are reported for Case 2 and Case 3 in Tab. 2. While
the order of magnitude of the growth rates obtained with the infinitely fast
chemistry is reasonable and has the correct trend compared to experiments, it
leads to growth rates about 2 times larger (third row) than the ones observed
in the experiments (second row). We conclude that a simple transport-limited
description of the reaction cannot accurately explain the growth rate of GaAs
via uncracked hydrides.

Table 2: Growth rates of Case 2 and 3 obtained with infinitely fast chemistry. |  $R_{G}$ [$\mu m/h$] for Case 2 |  $R_{G}$ [$\mu m/h$] for Case 3
---|---|---
Schulte  _et al._ (2018) | 28.8 | 63.2
Infinitely fast | 59.1 | 125.1

## References

  * King  _et al._ (2006) R. R. King, C. M. Fetzer, D. C. Law, K. M. Edmondson, H. Yoon, G. S. Kinsey, D. D. Krut, J. H. Ermer, P. Hebert, B. T. Cavicchi, _et al._ , “Advanced III-V multijunction cells for space,” in  _2006 IEEE 4th World Conference on Photovoltaic Energy Conference_ , Vol. 2 (IEEE, 2006) pp. 1757–1762.
  * Hoheisel, Philipps, and Bett (2010) R. Hoheisel, S. Philipps, and A. Bett, “Long-term energy production of III–V triple-junction solar cells on the Martian surface,” Progress in Photovoltaics: Research and Applications 18, 90–99 (2010).
  * Simon  _et al._ (2019) J. Simon, K. L. Schulte, K. A. Horowitz, T. Remo, D. L. Young, and A. J. Ptak, “III-V-based optoelectronics with low-cost dynamic hydride vapor phase epitaxy,” Crystals 9, 3 (2019).
  * McClure  _et al._ (2020) E. L. McClure, K. L. Schulte, J. Simon, W. Metaferia, and A. J. Ptak, “Gaas growth rates of 528 $\mu$ m/h using dynamic-hydride vapor phase epitaxy with a nitrogen carrier gas,” Applied Physics Letters 116, 182102 (2020).
  * Metaferia  _et al._ (2019) W. Metaferia, K. L. Schulte, J. Simon, S. Johnston, and A. J. Ptak, “Gallium arsenide solar cells grown at rates exceeding 300 $\mu$m/h by hydride vapor phase epitaxy,” Nature communications 10, 1–8 (2019).
  * Horowitz  _et al._ (2018) K. A. Horowitz, T. W. Remo, B. Smith, and A. J. Ptak, “A techno-economic analysis and cost reduction roadmap for III-V solar cells,” Tech. Rep. (National Renewable Energy Lab.(NREL), Golden, CO (United States), 2018).
  * Schulte  _et al._ (2018) K. L. Schulte, A. Braun, J. Simon, and A. J. Ptak, “High growth rate hydride vapor phase epitaxy at low temperature through use of uncracked hydrides,” Applied Physics Letters 112, 042101 (2018).
  * Hollan and Durand (1979) L. Hollan and J. Durand, “Fast growth in GaAs VPE at low temperature and high partial pressures,” Journal of Crystal Growth 46, 665–670 (1979).
  * DenBaars  _et al._ (1986) S. DenBaars, B. Maa, P. Dapkus, A. Danner, and H. C. Lee, “Homogeneous and heterogeneous thermal decomposition rates of trimethylgallium and arsine and their relevance to the growth of gaas by mocvd,” Journal of Crystal Growth 77, 188–193 (1986).
  * Shaw (1975) D. W. Shaw, “Kinetic aspects in the vapour phase epitaxy of III–V compounds,” Journal of crystal growth 31, 130–141 (1975).
  * Kangawa  _et al._ (2002) Y. Kangawa, T. Ito, Y. Hiraoka, A. Taguchi, K. Shiraishi, and T. Ohachi, “Theoretical approach to influence of As2 pressure on GaAs growth kinetics,” Surface science 507, 285–289 (2002).
  * Young  _et al._ (2013) D. L. Young, A. J. Ptak, T. F. Kuech, K. Schulte, and J. D. Simon, “High throughput semiconductor deposition system,” (2013), uS Patent App. 13/895,190.
  * Biefeld, Koleske, and Cederberg (2015) R. M. Biefeld, D. D. Koleske, and J. G. Cederberg, “The Science and Practice of Metal-Organic Vapor Phase Epitaxy (MOVPE),” in  _Handbook of Crystal Growth_ (Elsevier, 2015) pp. 95–160.
  * Schulte  _et al._ (2019) K. L. Schulte, W. Metaferia, J. Simon, and A. J. Ptak, “Uniformity of GaAs solar cells grown in a kinetically-limited regime by dynamic hydride vapor phase epitaxy,” Solar Energy Materials and Solar Cells 197, 84–92 (2019).
  * Schulte  _et al._ (2016) K. L. Schulte, J. Simon, N. Jain, D. L. Young, and A. J. Ptak, “A kinetic model for GaAs growth by hydride vapor phase epitaxy,” in  _2016 IEEE 43rd Photovoltaic Specialists Conference (PVSC)_ (IEEE, 2016) pp. 1930–1933.
  * Harrous  _et al._ (1988) M. Harrous, L. Chaput, A. Bendraoui, M. Cadoret, C. Pariset, and R. Cadoret, “Phosphine and arsine decomposition in CVD reactors for InP and InGaAs growth,” Journal of Crystal Growth 92, 423–431 (1988).
  * Braman, Oliver, and Raman (2013) K. Braman, T. A. Oliver, and V. Raman, “Bayesian analysis of syngas chemistry models,” Combustion Theory and Modelling 17, 858–887 (2013).
  * Bell  _et al._ (2019) J. Bell, M. Day, J. Goodman, R. Grout, and M. Morzfeld, “A bayesian approach to calibrating hydrogen flame kinetics using many experiments and parameters,” Combustion and Flame 205, 305–315 (2019).
  * Grüter  _et al._ (1989) K. Grüter, M. Deschler, H. Jürgensen, R. Beccard, and P. Balk, “Deposition of high quality GaAs films at fast rates in the LP-CVD system,” Journal of crystal growth 94, 607–612 (1989).
  * Pedregosa  _et al._ (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, “Scikit-learn: Machine Learning in Python,” Journal of Machine Learning Research 12, 2825–2830 (2011).
  * Hassanaly  _et al._ (2018) M. Hassanaly, H. Koo, C. F. Lietz, S. T. Chong, and V. Raman, “A minimally-dissipative low-Mach number solver for complex reacting flows in OpenFOAM,” Computers & Fluids 162, 11–25 (2018).
  * Hassanaly  _et al._ (2020) M. Hassanaly, Y. Tang, S. Barwey, and V. Raman, “Data-driven Analysis of Relight variability of Jet Fuels induced by Turbulence,” Combustion and Flame 225, 453–467 (2020).
  * Tang and Raman (2021) Y. Tang and V. Raman, “Large eddy simulation of premixed turbulent combustion using a non-adiabatic, strain-sensitive flamelet approach,” Combustion and Flame 234, 111655 (2021).
  * Koo  _et al._ (2017) H. Koo, M. Hassanaly, V. Raman, M. E. Mueller, and K. Peter Geigle, “Large-eddy simulation of soot formation in a model gas turbine combustor,” Journal of Engineering for Gas Turbines and Power 139 (2017).
  * Sitaraman  _et al._ (2021) H. Sitaraman, S. Yellapantula, M. T. H. de Frahan, B. Perry, J. Rood, R. Grout, and M. Day, “Adaptive mesh based combustion simulations of direct fuel injection effects in a supersonic cavity flame-holder,” Combustion and Flame 232, 111531 (2021).
  * Smagorinsky (1963) J. Smagorinsky, “General circulation experiments with the primitive equations: I. The basic experiment,” Monthly weather review 91, 99–164 (1963).
  * Tang  _et al._ (2019) Y. Tang, M. Hassanaly, V. Raman, B. Sforzo, and J. Seitzman, “A comprehensive modeling procedure for estimating statistical properties of forced ignition,” Combustion and Flame 206, 158–176 (2019).
  * Maestri and Cuoci (2013) M. Maestri and A. Cuoci, “Coupling CFD with detailed microkinetic modeling in heterogeneous catalysis,” Chemical Engineering Science 96, 106–117 (2013).
  * Strang (1968) G. Strang, “On the construction and comparison of difference schemes,” SIAM journal on numerical analysis 5, 506–517 (1968).
  * Heye, Raman, and Masri (2013) C. Heye, V. Raman, and A. R. Masri, “LES/probability density function approach for the simulation of an ethanol spray flame,” Proceedings of the Combustion Institute 34, 1633–1641 (2013).
  * Chong  _et al._ (2018) S. T. Chong, M. Hassanaly, H. Koo, M. E. Mueller, V. Raman, and K.-P. Geigle, “Large eddy simulation of pressure and dilution-jet effects on soot formation in a model aircraft swirl combustor,” Combustion and Flame 192, 452–472 (2018).
  * Oehlschlaeger  _et al._ (2009) M. A. Oehlschlaeger, J. Steinberg, C. K. Westbrook, and W. J. Pitz, “The autoignition of iso-cetane at high to moderate temperatures and elevated pressures: Shock tube experiments and kinetic modeling,” Combustion and flame 156, 2165–2172 (2009).
  * Kee  _et al._ (1999) R. Kee, G. Dixon-Lewis, J. Warnatz, M. Coltrin, J. Miller, and H. Moffat, “Transport: a software package for the evaluation of gas-phase, multicomponent transport properties,” Chemkin Collection (1999).
  * Bergman  _et al._ (2011) T. L. Bergman, F. P. Incropera, D. P. DeWitt, and A. S. Lavine,  _Fundamentals of heat and mass transfer_ (John Wiley & Sons, 2011).