Numerical Analysis of Gas Flow Instabilities in Simplified Vertical HVPE GaN Reactors

Abstract

This paper investigates the gas flow and the mass transport in simplified axial-symmetric vertical HVPE reactors for the growth of GaN bulk crystals through numerical simulations. We evaluate the relative significance of different flow and transport phenomena in dependence on the direction of gravity. The performed simulations show that buoyancy effects due to density differences between neighboring gas lines are the main factor causing the deformation of laminar flow patterns and the formation of recirculation cells within the growth zone. Baroclinic instabilities have been identified as the source for these phenomena. In contrast, typical vertical temperature gradients show only a minor impact on the stability of the gas flow within the growth zone in the vicinity of the growing crystal. Based on these results, major differences of the species transport in vertical HVPE reactors, where the flow is parallel or anti-parallel to the direction of gravity, referred to as down-flow and up-flow, respectively, are summarized. The performed analysis of the interplay and relative significance of different flow effects in the HVPE environment allows a general recommendation for reactor design and scaling with respect to stable gas flow conditions within the growth zone.

1. Introduction

Hydride vapor phase epitaxy (HVPE) is currently the main growth technique for the fabrication of bulk GaN crystals or the production of single freestanding GaN wafers [1,2,3]. Similar to other VPE methods, gas flow and mass transport phenomena strongly influence the uniformity of the growth rate, material composition and doping in the HVPE process and thus can affect the remaining built-in strain of the crystal [4,5]. Thus, the growth of GaN bulk crystals by HVPE requires a more precise control of the growth process in comparison to the fabrication of thin epitaxial layers. However, the nature of this method leads to a complex interplay between physical properties of the gas flow and the reactor geometry and therefore poses a number of challenges in the reactor design. HVPE of GaN metallic Ga, HCl and NH${}_3$ are precursors, and a mixture of N${}_2$ and H${}_2$ is used as the carrier gas. Initially and within a separate reactor zone, the so-called source zone, HCl reacts with the Ga melt surface under the formation of GaCl according to

$$2\mathrm{Ga}+2\mathrm{HCl}\leftrightarrow 2\mathrm{GaCl}+{\mathrm{H}}_2.$$

(1)

The formation of GaN on the substrate surface occurs within a second temperature zone at higher temperatures, the so-called growth zone, according to the following reaction

$$\mathrm{GaCl}+{\mathrm{NH}}_3\leftrightarrow \mathrm{GaN}+\mathrm{HCl}+{\mathrm{H}}_2.$$

(2)

This reaction takes place close to thermodynamic equilibrium and requires high partial pressures of GaCl and NH${}_3$ to achieve a reasonable growth rate [6]. In order to prevent NH${}_4$ Cl formation at low temperatures and parasitic depositions of GaN on the reactor parts at high temperatures, the precursors should be transported into the reactor by separated gas lines and should only start to mix at a short distance of a few cm away from the substrate. This means that the precursors in the HVPE reactor should be separated before they are nearly simultaneously mixed within the growth zone. These aspects distinguish HVPE strongly from other VPE growth techniques and require a quite unique and complex reactor setup for this growth technique. Therefore, the previous extensive studies of reactor flow dynamics, e.g., for CVD [7,8], cannot answer the questions of how to design the flow dynamics of the HVPE growth process and how to scale it. Indeed, strongly different molar masses of precursors and the components of the carrier gas result in significant density differences between the individual gas lines and their strong interaction in the area between the nozzles and the substrate (see Figure 1) at typical HVPE flow conditions. This can lead to severe flow instabilities due to buoyancy effects, which can make the control of the HVPE growth process extremely difficult.

There are quite a number of numerical studies on MOCVD (metalorganic vapour phase deposition) reactors and related processes of different geometrical concepts [9]. The reactor concepts studied were horizontal and vertical as well as vertical in up-flow or down-flow configurations. Included in former studies of MOCVD reactors, one can moreover find investigations of multi-wafer reactors of industrial scale and the effects of high speed rotation and process parameters on the layer quality and thickness homogeneity. The modelling of the quality and homogeneity of III-nitride epitaxial layers is unambiguously bound to the kinetics of the chemical processes on the surface of the substrate in conjunction with the temperature and flow patterns established in the reactor. The majority of publications on MOCVD reactors, therefore, focus on the further understanding of chemical reactions on the surface of the growing layer. In principle, the MOCVD reactors can be considered in the same way as the HVPE system in our study. The majority the numerical simulations existing deal with effects due to temperature fields, gas velocities and buoyancy; see e.g., [10,11,12,13,14]. However, the systems exhibit narrow distances between the gas inlet and the substrate position compared to an HVPE reactor, where the scaling is much larger. Our study can be seen as complementary to the published knowledge and will add fundamental value to the published studies as the inclusion of baroclinicity was not seen in the existing simulations and will give the additional ability to analyse and more thoroughly explain the effects of buoyancy and the appearance of vortices or instabilities in flow patterns inside systems with different geometries and process conditions. This will help to identify the most critical parameters (either in process parameters or reactor hardware) with regard to the optimization target (homogeneity) for an epitaxial layer or a thicker bulk crystal and sheds additional light on the so far not always well understood appearance of flow instabilities and possible counter measures.

Besides the numerical approach, ab initio theoretical methods are used to study nanoscale phenomena of epitaxial growth [15]. Comprehensive ab initio molecular dynamics [16] in conjunction with experimental validation can give additional important aspects to the understanding of epitaxial growth of III nitrides.

Therefore, the design of HVPE growth process is practically impossible without numerical simulations. Numerous publications provide CFD (computational fluid dynamics) investigations of horizontal [17,18,19,20] as well as vertical [21,22,23] HVPE reactors. Some studies [18,22] emphasize a large impact of gravity and buoyancy effects on the gas flow and mass transport in the HVPE reactor. However, the issues related to the design of the gas mixture and its flow stability for the given reactor geometry are still not investigated in full detail. Furthermore, the studied reactor setups are usually limited by a substrate diameter, which is up to 2 inches, and no numerical studies investigating the scaling aspects of HVPE for larger substrate diameters have been available until now.

This paper analyzes the stability of the gas flow in vertical HVPE reactors with down-flow and up-flow configurations for the growth of 3 inch GaN crystals through numerical simulations. A simplified vertical reactor setup with a conventional, so-called coaxial nozzle geometry is systematically investigated for different gas compositions and flow conditions using numerical simulations with OpenFOAM. The behavior of the gas flow will be described by means of non-dimensional Reynolds and Grashof numbers and their respective relationships. Based on these relationships, it is possible to derive criteria for the flow stability inside the reactor. Through the presented analysis, we are able to propose guidelines for the design of an HVPE reactor allowing the growth of GaN crystals under stable gas flow conditions in the growth zone.

The paper is structured as follows. In Section 2, we present a simplified model of the vertical HVPE reactor together with the boundary conditions and the numerical scheme. In Section 3, we perform a theoretical analysis of flow instabilities and their possible root cause in the expected flow regime by using the baroclinic term from the vorticity equation of the Navier–Stokes equation. Then, in Section 4, we investigate the occurrence and the dependence of flow instabilities with respect to different flow regimes with the help of numerical experiments. Finally, in Section 5, we propose guidelines for the setup of HVPE reactors and their gas inlets based on the previous findings. The presented results can be used as a guideline for the practical design of the vertical growth geometry and the gas composition for HVPE of GaN.

2. Modelling of Gas Flow in an HVPE Reactor

2.1. The Reactor Geometry

For our study, we use a simplified geometry of a vertical HVPE reactor; see Figure 1. This geometry does not consider the Ga source and is limited to the growth zone, including a coaxial nozzle arrangement and a substrate holder. Depending on the direction of gravity, this reactor setup can operate as a down-flow (same direction of gas flow and gravity) or up-flow (opposite direction of gas flow and gravity) HVPE reactor. Four separated coaxial gas lines are used for the transport of the precursors and the carrier gas. The middle gas line for the transport of GaCl is designated as inner run line (IRL). In order to prevent parasitic formation of GaN on the nozzles, a gas line with carrier gas, denoted as the separation line (SL), separates the IRL from the gas line containing ammonia, denoted as the outer run line (ORL). The last carrier gas line labeled as the side wall purge (SWP) helps to reduce parasitic growth of GaN on the outer reactor wall. Based on our experimental experience, the dimensions of the nozzle should allow a reasonable homogeneity of the growth rate at typical temperatures of around 1250 K on the 3-inch substrate placed on the rotated substrate holder at a distance of 80 mm from the nozzles.

2.2. The Fluid Dynamical System

For modelling the gas flow in the inner part of the reactor, we use the compressible Navier–Stokes equations in the following form:

$$\begin{array}{cc}\hfill {\rho }_t+\nabla ·\left(\rho \mathbf{u}\right)& =0,\hfill \\ \hfill {\left(\rho \mathbf{u}\right)}_t+\nabla ·(\rho \mathbf{u}\otimes \mathbf{u})+\nabla p& =\nabla ·\tau +\rho \mathbf{g},\hfill \\ \hfill E_t+\nabla ·\left(\mathbf{u}(E+p)\right)& =\nabla ·\lambda \nabla T+\langle \mathbf{u},\nabla ·\tau \rangle +\rho \langle \mathbf{u},\mathbf{g}\rangle ,\hfill \\ \hfill {\left(\rho Y_i\right)}_t+\nabla ·\left(\rho \mathbf{u}{Y}_i\right)& =\nabla ·D\nabla Y_i,\hfill \end{array}$$

(3)

where $\rho$ is the density, $\mathbf{u}$ is the fluid velocity, E is the total energy, D the diffusion coefficient, $\lambda$ is the thermal conductivity, $Y_i$ is the mass fraction of species i and $\mathbf{g}$ is the vector of gravitational acceleration. We assume the pressure p to satisfy the ideal gas law as

$$p=\rho \frac{R}{M}T,$$

(4)

where T is the temperature, M is the molar mass of the gas mixture and R is the universal gas constant. The viscous stress tensor $\tau$ is given by

$$\tau =\mu \left(\nabla \mathbf{u}+\nabla {\mathbf{u}}^T\right)-\mu \frac{2}{3}\nabla ·\mathbf{u}\mathbf{Id},$$

(5)

where $\mathbf{Id}$ is the $3\times 3$ identity matrix, and $\mu$ is the dynamic viscosity and is assumed to follow a Sutherland relation of the form

$$\mu =\overline{A}\frac{\sqrt{T}}{1+\frac{\overline{T}}{T}},$$

(6)

where the coefficients can be found in Appendix A. The diffusion coefficient D is determined via a constant Schmidt number of $Sc=1$ by the viscosity $\mu$, i.e.,

$$D=\frac{\mu }{\rho Sc}.$$

(7)

We justify this approximation from reference values of Schmidt numbers for representative gas pairs given in [24], where we find that $Sc\approx 1$. Lastly, the thermal conductivity is determined via a modified Eucken formula [25] given by

$$\lambda =\mu C_v\left(1.32+1.77\frac{R}{C_v}\right),$$

(8)

where $C_v$ is the heat capacity at constant volume. In the case of an ideal gas, there is

$$C_p-C_v=R,$$

(9)

defining $C_v$ in terms of the heat capacity at constant pressure $C_p$ and R. $C_p$ is given in polynomial form in Appendix A.

2.3. Boundary Conditions

At the periphery of the computational domain, fixed temperatures were used, which were obtained from a global thermal model of the whole HVPE reactor. The temperature profile at the inner reactor wall is shown in Figure 1. The inlet temperature in all gas lines was set to 1100 K. In each gas line, the gas mixture entered the reactor with a prescribed velocity ${\mathbf{u}}_j$, where $j\in \left\{\mathrm{SWP},\mathrm{ORL},\mathrm{SL},\mathrm{IRL}\right\}$; see also Figure 1. The pressure gradient $\nabla p$ was zero at all boundaries, except at the outlet, where a constant pressure $p=$ 91,192 Pa was used. The gradient of the mass fraction $Y_i$ of each species i was set to zero at all boundaries, except at the inlet j, where a certain mass fraction $Y_{i,j}$ entered the reactor. No-slip conditions for the gas flow velocities were applied at all solid walls. The substrate and substrate holder were rotating around the z-axis with a constant rotation rate of 6 rpm, and the no-slip boundary condition was adjusted here accordingly.

2.4. Numerical Scheme

The model was set up using the open source software OpenFoam v18 [26]. Specifically, we used the solver rhoreactingbuoyantfoam in order to numerically approximate solutions of the governing Equation (3) in their three-dimensional, time-dependent form with respect to the boundary conditions described in Section 2.3. The mesh was generated using the software package Gambit using the Cooper algorithm. The mesh generated then consisted of about 143,000 hexahedral cells. In terms of mesh quality, Gambit offers the metric EQUISIZE SKEW (ES), which is the skewness over the cell size. ES is in the range from 0 to 1, where 0 reflects the best quality and 1 the worst quality. We find that in the mesh used for the calculations, over $90\%$ of the cells were below 0.2, and none were over 0.65.

3. A Connection between Buoancy and Baroclinicity

In order to investigate the stability of the fluid flow in the reactor, we analyzed the so-called vorticity equation derived from the Navier–Stokes system (3) in its non-dimensional form. Under the assumptions that the reference variables scale as

$$u_0=\frac{x_0}{t_0},{\omega }_0=\frac{1}{t_0},{\tau }_0=\frac{{\mu }_0{u}_0}{x_0},$$

(10)

we obtain the vorticity equation in the following form

$${\omega }_t+(\mathbf{u}·\nabla )\omega -(\omega ·\nabla )\mathbf{u}+\omega \nabla ·\mathbf{u}=\frac{1}{\gamma M{a}^2}\frac{\nabla \rho \times \nabla p}{{\rho }^2}+\frac{1}{Re}\frac{\mu }{\rho }\Delta \omega -\frac{1}{Re}\frac{\nabla \rho \times \nabla ·\tau }{{\rho }^2},$$

(11)

where $\omega =\nabla \times \mathbf{u}$ is the vorticity vector and

$$Ma=\frac{u_0}{c_0},Re=\frac{x_0{u}_0{\rho }_0}{{\mu }_0},c_0=\sqrt{\gamma \frac{p_0}{{\rho }_0}}$$

(12)

are the Mach number, Reynolds number and the speed of sound, respectively. $\gamma$ is the isentropic coefficient and can be expressed as the ratios of heat capacities, i.e., $\gamma =\frac{C_p}{C_v}$.

We point out that, in our model, there is $\nabla \times \mathbf{g}=0$, and therefore the gravitational source term does not directly appear in (11). However, the connection to buoancy can be made by considering the relevant fluid regime. We expect the flow regime to be in a range of $Ma\approx {10}^{-3}$. Therefore, we find that in (11), the first term on the right-hand side, the baroclinic term [27], scaling with the inverse of $M{a}^2$ is of great importance. Hence, we now focus on the composition of the baroclinic term.

Due to the analysis presented in the Appendix B, we find that for flows at low Mach numbers, the pressure can be approximated by its hydrostatic part satisfying the relation (A8). Therefore, using (4), (A7), (A8) and (14), we can rewrite the baroclinic term in (11) as follows

$$\frac{1}{\gamma M{a}^2}\frac{\nabla \rho \times \nabla p}{{\rho }^2}\approx \frac{G{r}_M}{R{e}^2}\nabla lnM\times \mathbf{g}-\frac{G{r}_T}{R{e}^2}\nabla lnT\times \mathbf{g},$$

(13)

by using the following dimensionless Grashof numbers:

$$G{r}_T=\frac{{\rho }_0^2{g}_0{x}_0^3\nabla T_0}{{\mu }_0^2{T}_0},G{r}_M=\frac{{\rho }_0^2{g}_0{x}_0^3\nabla M_0}{{\mu }_0^2{M}_0}.$$

(14)

Equation (13) shows the direct connection of the baroclinic term to buoyancy at low Mach numbers due to the appearance of the Grashof numbers. See also [28], where the baroclinic term is used to explain the well-known buoancy-driven Rayleigh–Taylor instability. In fact, Equation (13) also provides the basis to distinguish buoancy effects based on thermal or molar gradients in the fluid. Moreover, since the gravitational acceleration $\mathbf{g}$ is considered constant, the effect of the molar or thermal gradients on the fluid flow can be easily computed by evaluating the cross-product of the respective gradient with the constant vector $\mathbf{g}$. Finally, since changing from the up-flow to the down-flow configuration of the HVPE reactor and vice versa could be understood as just flipping the sign of the gravitational acceleration $\mathbf{g}$, Equation (13) already suggests that the flow in the different reactors might have a different response to the same thermal and molar gradients.

We finally use (13) in (11) to find that

$$\begin{array}{ll}{\omega }_t+(\mathbf{u}·\nabla )\omega & -(\omega ·\nabla )\mathbf{u}+\omega \nabla ·\mathbf{u}\\ & =\frac{1}{Re}\left(\frac{G{r}_M}{Re}\nabla lnM\times \mathbf{g}-\frac{G{r}_T}{Re}\nabla lnT\times \mathbf{g}+\frac{\mu }{\rho }\Delta \omega -\frac{\nabla \rho \times \nabla ·\tau }{{\rho }^2}\right).\end{array}$$

(15)

We conclude from (15) that in the low Mach number flow regime, the behaviour of the vorticity is described by three non-dimensional numbers, i.e., the Reynolds number $Re$ and the two Grashof numbers $G{r}_M$ and $G{r}_T$. We will make use of these non-dimensional numbers and their respective ratios as given in (15) to describe the respective fluid regimes for which we conduct the numerical experiments.

To close this section and as an outlook for different reactor designs, we would like to add some remarks. First, due to the properties of the cross-product, we conclude from Equation (13) that only horizontal molar or thermal gradients contribute to the baroclinic term. Therefore, we focus our analysis in the numerical experiments exactly on those lateral gradients. Moreover, this point of view of buoyancy may also be used in horizontal reactors as, for example, studied in [18,29,30], since it depends only on the interplay of horizontal gradients of molar mass or temperature in the reactor with gravity. Then, the baroclinic term describes the resulting forces on the fluid due to buoyancy. However, the analysis of buoyancy effects through the baroclinic term only makes sense if transport phenomena through advection are relevant in the model, since the baroclinic term is derived from the pressure term in the momentum equation.

4. Numerical Experiments

We perform several numerical experiments to evaluate the flow stability in the down-flow and up-flow reactor. By using simplified gas compositions in the run lines, we propose a novel approach for the analysis of the gas flow stability in vertical HVPE reactors. Starting from the laminar N${}_2$ flow, we add for the growth reasonable amounts of NH${}_3$ or GaCl/H${}_2$ in the ORL or IRL, respectively, and investigate its effect on the stability of the gas flow. We base our investigations on the stability in different fluid regimes on the molar Grashof $G{r}_M$, the thermal Grashof $G{r}_T$ and the Reynolds numbers $Re$. The definition of the respective reference values are given in

Table 1. Reference values used to determine non-dimensional numbers.

Reference Quantity	${\mathit{x}}_0$	${\mathit{g}}_0$	${\mathit{\rho }}_0$	${\mathit{\mu }}_0$	$\mathit{\nabla }\mathit{ln}\left({\mathit{T}}_0\right)$
value	$0.3$ m	$9.81$$\frac{\mathrm{m}}{{\mathrm{s}}^2}$	$1.1$$\frac{\mathrm{kg}}{{\mathrm{m}}^3}$	$5.0\times {10}^{-5}\frac{\mathrm{kg}}{\mathrm{m}s}$	0.1 $\frac{1}{m}$

. Therefore, it should be mentioned at this point that the thermal Grashof number $G{r}_T$ is constant in all simulations because the thermal boundary conditions are not changed, i.e., $G{r}_T\approx 1.3\times {10}^7$.

The reference velocity $u_0$ is defined as the area-averaged mean velocity at the inlets, i.e.,

$$u_0=\frac{{\sum }_j{u}_j{A}_j}{{\sum }_{j=1}{A}_j},$$

(16)

where $A_j$ denotes the surface of the inlet $j\in \left\{\mathrm{SWP},\mathrm{ORL},\mathrm{SL},\mathrm{IRL}\right\}$. In order to define the inlet velocities, $u_j$, we consider reference inlet velocities relative to the inlet velocity $u_{IRL}=\pm 0.172$ m/s, where $"+"$ denotes the up-flow reactor and $"-"$ denotes the down-flow reactor. Then, in order to vary the Reynolds number, we define the inlet velocities by scaling the reference values by a constant factor. According to (16), we allow negative and positive Reynolds numbers in (12) to refer to the down-flow and up-flow reactor, respectively.

Furthermore, we define $\frac{\nabla M_0}{M_0}=\nabla ln\left(M_0\right)$ in terms of the possible four different molar masses at the inlets. Due to the axial symmetric reactor geometry, we obtain three different possible values, i.e.,

$$\nabla ln{\left(M_0\right)}_{j_2}^{j_1}=\frac{ln\left(M_{j_1}\right)-ln\left(M_{j_2}\right)}{x_0},(j_1,j_2)\in \{\left(\mathrm{SWP},\mathrm{ORL}\right),\left(\mathrm{ORL},\mathrm{SL}\right),\left(\mathrm{SL},\mathrm{IRL}\right)\}.$$

(17)

We use each of the above definitions in the respective simulations. Lastly, we propose a definition of unstable flow behavior. We aim for a criterion that reflects the occurrence of vortices and recirculizations in the reactor.

Definition. 

A flow is unstable in a domain Ω over the time $t\in \left[T_0,T_1\right]$ if

$$\mathbf{I}\left(\mathrm{\Omega }\right):={\int }_{T_0}^{T_1}{\int }_{\mathrm{\Omega }}\chi (\mathbf{u},\mathbf{v})dxdt>0,$$

(18)

where

$$\chi (\mathbf{u},\mathbf{v})=\left\{\begin{array}{c}1ifarccos\left(\frac{\langle \mathbf{u},\mathbf{v}\rangle }{\mathbf{u}\mathbf{v}}\right)>{100}^{\circ },\hfill \\ 0\mathit{else},\hfill \end{array}\right.\mathit{with}\mathbf{v}=\left\{\begin{array}{cc}-\mathbf{g}\hfill & \mathit{for}\mathit{Up}-\mathit{Flow}\hfill \\ \mathbf{g}\hfill & \mathit{for}\mathit{Down}-\mathit{Flow}\hfill \end{array}\right..$$

(19)

In other words, the function defined in (19) identifies regions where the flow velocity vector deviates more then ${100}^{\circ }$ from its intended direction, which in the up-flow case is against the vector of gravity and in the down-flow case along the vector of gravity. Accordingly, the function $\mathbf{I}\left(\mathrm{\Omega }\right)$ defined in (18) is able to identify if instabilities occur over a defined time interval in the specified domain $\mathrm{\Omega }$. In order to distinguish different flow phenomena, we split the reactor ${\mathrm{\Omega }}_{HVPE}$ in an outer part and in an inner part, denoted by ${\mathrm{\Omega }}_{Out}$ and ${\mathrm{\Omega }}_{In}$, respectively, with ${\mathrm{\Omega }}_{HVPE}={\mathrm{\Omega }}_{Out}\cup {\mathrm{\Omega }}_{In}$. In Figure 2, we depict the definition of ${\mathrm{\Omega }}_{Out}$ and ${\mathrm{\Omega }}_{In}$ for the down-flow reactor. The definition for the up-flow reactor is analogous following Figure 1.

4.1. Variation of NH${}_3$ in ORL

First we investigated the variation of ${\mathrm{NH}}_3$ in the ORL for the up-flow and the down-flow reactor. Apart from the ammonia in the ORL, only ${\mathrm{N}}_2$ was used in all the other gas lines. Following (17), we used $\nabla ln{\left(M_0\right)}_{\mathrm{ORL}}^{\mathrm{SWP}}$ in (14) to define the molar Grashof number. Consequently, from the molar masses of ${\mathrm{NH}}_3$ and ${\mathrm{N}}_2$, we found that $\nabla ln{\left(M_0\right)}_{\mathrm{ORL}}^{\mathrm{SWP}}>0$, and therefore, $G{r}_M>0$ in this simulation series.

In Figure 3, snapshots of the distribution of $Y_{NH3}$ in a vertical cut of both the up-flow and the down-flow reactor for $Re=\pm 600$ and $G{r}_M=8.5\times {10}^7$ are shown. It is obvious that the gas flow behaves differently for the two setups. In the up-flow configuration, the distribution of $Y_{{\mathrm{NH}}_3}$ is symmetric and follows the general direction of the flow, while in the down-flow configuration, the ammonia also flows into the outer and upper part of the reactor.

We performed multiple simulations at different Reynolds and Grashof numbers and checked for the occurrence of instabilities. Figure 4 shows exemplarily typical regions where the instabilities in the down-flow and up-flow configuration occur for different Re and Gr numbers.

In Figure 4a, the instabilities are present mainly in the outer part of the growth zone below and above the position of the gas nozzles for the down-flow reactor, while in Figure 4c, the instabilities are located exclusively in the upper part between the SWP gas line and the reactor wall. The most critical disturbance for crystal growth of GaN in the down-flow configuration is shown in Figure 4b as the instabilities occur directly in the vicinity of the growing crystal in the inner part of the growth zone as well as in the outer part.

Typical locations of the instabilities in the up-flow configuration are shown in Figure 4d–f. Again, we can distinguish less harmful cases, where the instabilities occur only in the outer part (see Figure 4f) and more severe cases, where the instabilities are located close to the growing crystal in the inner part (see Figure 4d) or even in the inner and outer part (see Figure 4e).

Figure 5 summarizes the results of the variation of NH${}_3$ in the ORL for the down-flow (see Figure 5a) and up-flow (see Figure 5b) configurations in the form of a stability diagram in the $G{r}_M$ / $Re$ over $G{r}_T$ /$Re$ parameter space. For the down-flow reactor, it is obvious that for $G{r}_T/Re<-0.3\times {10}^5$ and independent of the $G{r}_M/Re$ ratio, the flow is always unstable and the instabilities occur mainly in the outer parts; see, for example Figure 4c. For $G{r}_T/Re>-0.3\times {10}^5$ and $G{r}_M/Re>-1.2\times {10}^5$, the flow is stable without any disturbance, while the flow becomes unstable again for $G{r}_T/Re<-0.3\times {10}^5$ and $G{r}_M/Re<-1.2\times {10}^5$ with disturbances in the outer part of the growth zone; see Figure 4a. Figure 4b refers to a fluid regime outside any of the previous described stability boundaries for $G{r}_M/Re$ and $G{r}_T/Re$. Moreover, we find that for a given $G{r}_M$, the stability region can be extended by increasing the Re number.

The stability regions for the up-flow reactor behave differently from the down-flow configuration, as shown in Figure 5b. Instabilities in the outer part of the reactor are found only typically for $G{r}_T/Re>0.5\times {10}^5$; compare Figure 4f. This is slightly higher than in the down-flow configuration, where the transition took place at $G{r}_T/Re<-0.3\times {10}^5$. Opposite to the down-flow reactor, no stability region exists in the up-flow configuration, and the flow is always unstable for $G{r}_T/Re<0.5\times {10}^5$, where instabilities occur mostly in the inner part of the reactor; see Figure 4d. Thereby, it is obvious that there exists a transition region between $0.2\times {10}^5<G{r}_T/Re<0.5\times {10}^5$. Within this region, instabilities occur in the inner and outer part; see Figure 4e.

In the next section, we present the results regarding the stability of the flow with respect to a variation of the gas mixture in the IRL.

4.2. Variation of H${}_2$ in IRL

In the second case study, we considered the variation of ${\mathrm{H}}_2$ in the IRL, where a mixture of GaCl, ${\mathrm{N}}_2$ and ${\mathrm{H}}_2$ was used. In all the other lines, only ${\mathrm{N}}_2$ was applied. We varied the amount of $Y_{{\mathrm{H}}_2,IRL}$ while keeping $Y_{GaCl,IRL}=0.212$ constant until a maximum such that the molar mass of the IRL equaled that of pure ${\mathrm{N}}_2$. Following (17), we use $\nabla ln{\left(M_0\right)}_{IRL}^{SL}$ in (14) to define the molar Grashof number. Consequently, there was $\nabla ln{\left(M_0\right)}_{IRL}^{SL}<0$, and therefore, $G{r}_M<0$ in this simulation series.

In Figure 6, snapshots of the distribution of $Y_{GaCl}$ in a vertical cut of both the up-flow and the down-flow reactors for $Re=\pm 600$ and $G{r}_M=-7.2\times {10}^7$ are shown. Again, the gas flow behaves differently for the two setups. While in the down-flow reactor, the GaCl flows laminarly from the nozzle to the growing crystal, in the up-flow reactor, the GaCl spreads significantly towards the outer part of the reactor after leaving the nozzle.

Again we performed a case study for different Reynolds and Grashof numbers by varying the inlet velocities, and we progressively substituted ${\mathrm{H}}_2$ by ${\mathrm{N}}_2$ in the IRL until there was no molar mass gradient. Figure 7 shows exemplary typical regions where the instabilities in the down-flow and up-flow configuration occur for different Re and Gr numbers. In the down-flow reactor, the instabilities occur only in the outer part of the reactor, while for the up-flow configuration, the flow is disturbed in the inner and outer part depending on the Re and Gr numbers.

Figure 8 summarizes the results of the variation of H${}_2$ in the IRL for the down-flow (see Figure 8a) and up-flow (see Figure 8b) configurations in the form of a stability diagram in the $G{r}_M/Re$ over $G{r}_T/Re$ parameter space. For the down-flow geometry, when the gas mixture in the IRL is varied, there also exists a critical $G{r}_T/Re$ ratio. For $G{r}_T/Re<-0.3\times {10}^5$, the flow is always unstable in the outer part of the growth zone (compare Figure 7a), while for $G{r}_T/Re>-0.3\times {10}^5$, no disturbance of the flow occurs at all in the considered parameter range. Interestingly, the transition takes place almost at the same critical $G{r}_T/Re$ ratio for a variation of the gas mixture in the IRL as in the case where we varied the gas composition in the ORL. However, in the case of the variation of the gas composition of the ORL, instabilities were also found in the inner part for a high $G{r}_M/Re$ ratio, which do not exist in the present case.

For the up-flow reactor (see Figure 8b), there are again three stability regions. For $G{r}_T/Re>0.5\times {10}^5$, instabilities are present mainly in the outer areas (compare Figure 7d), whereas for $G{r}_T/Re<0.3\times {10}^5$, the flow is disturbed in the inner part close to the growing crystals (compare Figure 7b). In between, a stable flow regime exists for $G{r}_M/Re<-1.0\times {10}^5$.

4.3. Analysis of the Instabilities for the Down-Flow Reactor

In order to find the different causes for the observed flow instabilities, we recall the low Mach number approximation of the vorticity Equation (15), i.e.,

$$\begin{array}{ll}{\omega }_t+(\mathbf{u}·\nabla )\omega & -(\omega ·\nabla )\mathbf{u}+\omega \nabla ·\mathbf{u}\\ & =\frac{1}{Re}\left(\frac{G{r}_M}{Re}\nabla lnM\times \mathbf{g}-\frac{G{r}_T}{Re}\nabla lnT\times \mathbf{g}+\frac{\mu }{\rho }\Delta \omega -\frac{\nabla \rho \times \nabla ·\tau }{{\rho }^2}\right).\end{array}$$

(20)

From Equation (20), it is obvious that due to the cross-product, not only the magnitude, but also the direction of the molar and thermal gradients with respect to the direction of gravity are important.

From the thermal boundary conditions shown in Figure 1, we expect stable thermal conditions for the down-flow configuration above the growing crystal because the gas is hottest at the nozzles and coldest at the crystal, while in the upper reactor part from the gas inlet to the gas outlet at the nozzles, the thermal conditions are considered to be unstable as it is coldest at the inlet and hottest at the outlet. For the up-flow configuration, the thermal conditions are reversed, which means it is unstable under the growing crystal and stable along the gas lines.

In addition, we have to consider that horizontal gradients of the molar mass within the different gas lines exist as well as a horizontal temperature difference within each gas line because cold gas enters the hot reactor, leading to a cold jet in each gas line; see Figures 10 and 12. Since the vector of gravitational acceleration $\mathbf{g}$ is vertical and for every vector $\mathbf{w}$ there is $\mathbf{w}\times \mathbf{w}=0$, the baroclinic term is expected to depend more on the horizontal variations of the molar mass and the temperature field.

Figure 9 shows the distribution of the x-component of the baroclinic term for the variation of NH${}_3$ Figure 9a and GaCL/H${}_2$ Figure 9b in the down-flow reactor. The x-component of the baroclinic term acts as a source for the vorticity component, resembling rotations in the y-z-plane. This is also referred to as baroclinic torque [28]. The absolute value of the baroclinic torque represents the acceleration of the fluid, while its sign indicates whether the torque is in the positive (counterclockwise) or negative (clockwise) direction.

We start our analysis by identifying regions in the reactor that show a large amount of baroclinicity. For both cases in Figure 9, we find baroclinic torque to be active at the outer walls of the reactor and beneath the nozzle outlet between the ORL and the SWP, where NH${}_3$ is varied in the ORL, and between the IRL and the SL, where H${}_2$ is changed in the IRL.

In the case of the NH${}_3$, in the ORL, the baroclinic torque accelerates the gas against its intended flow direction (see Figure 9a), reducing the vertical fluid velocity beneath the nozzle and transporting NH${}_3$ towards the outer part of the reactor. Therefore, instabilities in the outer part of the reactor (compare Figure 4a), and in the inner and outer part (compare Figure 4b), can be triggered with an increase in $G{r}_M$. The approximation of the baroclinic term in (13) lets us conclude that higher fluid velocities, i.e., higher Re numbers, might suppress the rising of NH${}_3$ in the outer part of the reactor, as can also be seen in Figure 5a, since with higher fluid velocities, one moves on rays towards the origin in the plot.

On the other hand, in the case of GaCl in the IRL, the baroclinic torque accelerates the gas along its intended flow direction, leading to a stable flow pattern in terms of the dependence on the molar mass gradients. This is further underlined by the results in the Figure 7 and Figure 8a, as no instabilities in the inner part of the reactor above the growing crystal are observed.

Next, we analyze the baroclinicity at the outer reactor walls. The baroclinic torque shown by the white circles in Figure 9 has the same direction for both cases, variation of NH${}_3$ in the ORL (see Figure 9a) and variation of GaCl/H${}_2$ ratio in IRL (see Figure 9b). It is expected that the baroclinic torque accelerates the gas towards the reactor side wall, respectively, the quartz cylinder separating the SWP from the ORL, which will promote the occurrence of instabilities. From Figure 10, it is clear that this baroclinic torque at the outer walls is induced by horizontal thermal gradients. These horizontal gradients originate from the heating of the gas through the reactor wall while flowing downwards towards the mixing zone. This hypothesis is supported by the results from Figure 4 and Figure 7, where we find that the instabilities occur in both cases, i.e., variation of NH${}_3$ in the ORL and variation of GaCl/H${}_2$ ratio in IRL, at the same critical $G{r}_T/Re$ ratio, i.e., $G{r}_T/Re<-0.3\times {10}^5$. However, from a practical point of view. it can be stated that usually in real crystal growth experiments, the Reynolds number is relatively high and therefore should be not in a regime where these instabilities driven by horizontal temperature gradients occur.

4.4. Analysis of the Instabilities for the Up-Flow Reactor

In case of the up-flow reactor, we again turn to the analysis of the effect of the baroclinic torque. Figure 11 shows the distribution of the x-component of the baroclinicity in the up-flow reactor for the variation of NH${}_3$ in the ORL (see Figure 11a) and variation of $GaCl/H_2$ ratio in IRL (see Figure 11b). In an analogy to the down-flow reactor, two regions of high baroclinicity can be identified, namely the outer walls and the area between the nozzle and the growth zone.

From Figure 4 and Figure 7 together with Figure 5 and Figure 8, we find that the instabilities in the outer part correlate with the $G{r}_T/Re$ ratio rather than the $G{r}_M/Re$ ratio and occur especially at high values of $G{r}_T/Re>0.5\times {10}^5$, i.e., at small Re numbers. From Figure 12, it can be concluded that these instabilities are again triggered by horizontal thermal gradients, as in the case of the down-flow reactor. However, one difference from the down-flow reactor is the small dependence of these instabilities on the $G{r}_M$ in the case of NH${}_3$ in the ORL. Looking at the distribution of $Y_{N{H}_3}$ in Figure 3 and the baroclinic torque at the nozzle in Figure 11, we find that the baroclinic torque becomes stronger with an increasing amount of NH${}_3$, pushing the gas flow towards the center. This results in a region with reduced fluid velocity at the outer wall in the mixing zone. This is also in line with the result from Figure 5 that the region with instabilities extends with increasing $G{r}_M$. Increasing the Re number will counteract this effect and suppress the instabilities in this part of the reactor.

Now, we focus on instabilities that occur in the inner part of the up-flow reactor. We investigate the results from the variation of GaCl and H${}_2$ in the IRL first. From Figure 11b, it can be seen that the baroclinic torque at the nozzle outlets caused by molar mass gradients is responsible for an acceleration of the fluid against its intended flow direction. This effect is accompanied by a transport of the gas away from the center of the reactor. This is underpinned by the results from Figure 7 showing that these instabilities can be suppressed by decreasing $G{r}_M$; see also Figure 8b.

Now, we investigate the variation of NH${}_3$ in the ORL; see Figure 11a. In fact, the molar mass gradient between SWP and ORL is always positive, i.e., $\nabla ln{\left(M_0\right)}_{ORL}^{SWP}>0$. However, the molar mass gradient between ORL and SL is always negative, i.e., $\nabla ln{\left(M_0\right)}_{SL}^{ORL}<0$. This results in a locally negative molar Grashof number, i.e., $G{r}_M<0$. Therefore, these instabilities follow the same root cause as in the case of the GaCl and H${}_2$ in the IRL, namely, the acceleration of the fluid due to baroclinic torque induced by molar mass gradients is against its intended flow direction and therefore leads to instabilities. This is again underpinned by the results from Figure 5b, where one finds that the condition of stability on the Reynolds number relaxes as $G{r}_M$ decreases.

5. Consequences for Gas Inlet Designs and Reactor Scaling

In this chapter, we give an outlook on the consequences that we find through our analysis both for the design of stable gas mixtures for the up-flow and the down-flow reactor as well as for further reactor designs.

5.1. Towards Stable Gas Inlet Configurations

By the analysis of our numerical experiments, we conclude that the molar masses of the gases at the different inlets have a strong influence on the stability of the gas flow in the reactor; see Figure 5 and Figure 8. By considering the definition of the molar Grashof number (12) and the molar gradients (17), a stability criterion for the relation of the molar masses at each inlet can be defined.

First, we consider the down-flow reactor. In the case of $M_{ORL}<M_{SWP}$, instabilities may occur if the molar Grashof number gets too large or the fluid velocity is too small. However, in the case of $M_{IRL}\ge M_{SL}$, we find no instabilities. Due to the axisymmetric reactor geometry, we deduce the following stability criterion for the down-flow reactor:

$$M_{IRL}\ge M_{SL}\ge M_{ORL}\ge M_{SWP}.$$

(21)

However, it should be noticed that small violations of the stability criterion (21) might not immediately lead to instabilities as long as they can be suppressed by a sufficiently high fluid velocity; see also Figure 5.

For the up-flow reactor, we find that the condition $M_{ORL}\le M_{SWP}$ produces stable behavior, especially in the outer part of the reactor. However, in this case, unstable behavior beneath the growing crystal is found since $M_{SL}>M_{ORL}$. This corresponds to the instabilities found in the case where $M_{IRL}>M_{SL}$. Therefore, we formulate the stability criterion for the up-flow reactor as

$$M_{IRL}\le M_{SL}\le M_{ORL}\le M_{SWP}.$$

(22)

5.2. Scaling of Vertical HVPE Reactors

Lastly, we concern ourselves with consequences from the analysis of the flow stability for the scaling of vertical HVPE reactors. For this purpose, we use Figure 5 and Figure 8 and create Figure 13. We recall the definitions of the Grashof and Reynolds numbers from (12) and (14), respectively. Since we are concerned with the scaling of the vertical HVPE reactors, it is important to see that

$$\frac{G{r}_M}{Re}=O\left(x_0^2\right)=\frac{G{r}_T}{Re}.$$

(23)

Therefore, if an up-scaling of the vertical HVPE reactors considered here is proposed, then given a certain base configuration that leads to a fluid regime that can be found in Figure 13, one moves this base configuration on a straight line away from the origin, i.e., along the pink arrows in Figure 13. However, moving away from the origin leads to unstable configurations. If in that case a stabilization can not be reached by applying the stability criterions (21) or (22) on the molar masses, a stabilization through higher fluid velocities can be considered since

$$\frac{G{r}_M}{Re}=O\left(u_0^{-1}\right)=\frac{G{r}_T}{Re}.$$

(24)

Therefore, if it is decided to increase the fluid velocity from a given base configuration, then one moves in Figure 13 in the opposite direction than in the case of varying $x_0$, i.e., along the green arrows towards the origin in Figure 13. At least in the down-flow reactor, there are stable configurations. However, given that one starts in a stable configuration and wants to stay stable through a reactor scaling, it is necessary by considering (23) and (24) that

$$u_0=O\left(x_0^2\right),$$

(25)

i.e., the velocity $u_0$ needs to scale with the square of the length scale $x_0$. However, this can be fulfilled only in theory. In practice, one would expect that, on the one hand, one reaches the boundary to turbulent flow at a certain point for large $x_0$. On the other hand, the higher gas velocities for larger reactor sizes would also be connected with much higher costs for the gases.

6. Conclusions

We analyzed two types of HVPE reactors for the growth of GaN bulk crystals regarding their stability with respect to gas inlet compositions and the orientation of the gas flow direction with respect to the gravitational vector. Therefore, three-dimensional, time-dependent numerical simulations have been performed, and the stability of the gas flow has been analyzed in dependence of the Reynolds number Re and the molar as well as the thermal Grashof numbers, $G{r}_M$ and $G{r}_T$, respectively. We identified the baroclinic term as the main contributor of the observed instabilities. The baroclinic term is sensitive to the contributions of molar mass gradients and horizontal temperature gradients. Based on these findings, we have been able to formulate a stability criterion in terms of molar masses for the down-flow reactor and for the up-flow reactor and to give further insight into arguments with respect to the scaling of vertical HVPE reactors.