Numerical Modelling of the Czochralski Growth of β-Ga₂O₃

Abstract

Our numerical modelling of the Czochralski growth of single crystalline β-Ga${}_2$O${}_3$ crystals (monoclinic symmetry) starts at the 2D heat transport analysis within the crystal growth furnace, proceeds with the 3D heat transport and fluid flow analysis in the crystal-melt-crucible arrangement and targets the 3D thermal stress analysis within the β-Ga${}_2$O${}_3$ crystal. In order to perform the stress analysis, we measured the thermal expansion coefficients and the elastic stiffness coefficients in two samples of a β-Ga${}_2$O${}_3$ crystal grown at IKZ. Additionally, we analyse published data of β-Ga${}_2$O${}_3$ material properties and use data from literature for comparative calculations. The computations were performed by the software packages CrysMAS, CGsim, Ansys-cfx and comsol Multiphysics. By the hand of two different thermal expansion data sets and two different crystal orientations, we analyse the elastic stresses in terms of the von-Mises stress.

1. Introduction

Oxide single crystals are the prerequisite for many applications, as e.g., hosts for lasers, frequency doublers in lasers, sensors exploiting the piezoelectric effect, and substrates for functional layers of oxides. In particular, β-Ga${}_2$O${}_3$ as a transparent semiconducting oxide has recently become of great interest for the design of new types of electronic devices [1]. Bulk crystals of this material can act as a substrate not only for homoepitaxy (see e.g., [2]), but also for epitaxy of mixed systems of (In${}_x$Ga${}_{1-x}$)${}_2$O${}_3$ [3] or (Al${}_x$Ga${}_{1-x}$)${}_2$O${}_3$ [4]. One way of growing single crystals of β-Ga${}_2$O${}_3$ is the Czochralski method, where the crystal is pulled from a melt hosted in an iridium crucible. Heat is induced in the crucible by electro-magnetic heating. Single crystals have been succesfully grown by this method [5,6,7,8]. The major issue associated with growing Ga${}_2$O${}_3$ single crystals from the melt is its thermal instability at high temperatures leading to a strong decomposition. Applying a new approach [5] allowed for growing two-inch diameter crystals with a weight of up to 1 kg. Additionally, free carrier absorption within a growing semiconducting crystal causes problems with heat removal from the growth interface, leading to its inversion from convex to concave and, consequently, to a spiral formation [6]. Growing large crystal size enhances both decomposition (larger melt volume) and free carrier absorption (larger optical thickness of the crystal), and also thermal stresses, which may affect the structural quality, such as dislocation formation, twinning and cracking. Numerical calculations can contribute to a better understanding and eventually lead to improvements of the growth process.

The growth process of β-Ga${}_2$O${}_3$ includes all of the challenges for numerical modelling of oxide crystal growth, namely high melting point temperature causing radiation as an important heat transport mechanism, semitransparency of crystal depending on the doping, time dependent melt flow, formation of facets, and anisotropy of several properties. The latter three points usually require a 3D calculation and, thus, a high amount of computational resources. At first glance, for computing melt convection, a local calculation including the melt and maybe crucible and crystal is enough because the three -dimensionality results only from the flow instability. Oxide melts have a Prandtl number larger than 1 ($\mathrm{Pr}=c_p\mu /\lambda$, $c_P$ specific heat, μ dynamic viscosity, λ thermal conductivity), which means a significant heat transport by convection. In contrast to typical semi-conductor problems, the melt flow is time-dependent but far away from turbulence. Some researchers looked for the onset of the flow instability in oxide melts [9], but the realistic conditions are far beyond the ones considered there. For more realistic Grashof numbers, 3D computations have been performed but only for the melt applying simple boundary conditions at its side, bottom and top [10,11]. The bottom was set adiabatic and a constant temperature was set at the side. At the melt/crystal interface, there was melting point temperature and the free melt surface at the top was also assumed to be adiabatic. Enger et al. [11] performed computations for a melt of Prandtl number Pr = 10.5 and, later, Jing et al. [10] performed runs for Pr = 13.6, as it is the Prandtl number of LiNbO${}_3$ melt. The computations comprised buoyancy and Marangoni convection as well as crystal rotation and exhibited a time-dependent behaviour of the convection. The temperature distribution on the free melt surface showed a structure, which changed in time. The pattern depends on the rotation rate, as was found by Jing et al. However, there are strong restrictions in these computations because of the adiabatic assumption at the free surface of the melt. In real systems, a heat exchange by radiation takes place, adding another non-linearity (the radiation heat flux depends on $T^4$) the non-linearity of the Navier–Stokes equations governing the melt convection. Therefore, Tsukada et al. introduced a coupling of a global axisymmetric calculation with a local 3D one [12]. In the local calculation comprising melt, crucible, and crystal, the time-dependent flow was computed using temperature boundary conditions from the global axisymmetric calculations. Afterwards, the time-averaged velocity field is used in the global calculation for updating the temperature field and the melt/crystal interface shape. In this study, we did not make any back coupling from 3D to axisymmetric calculations but included the complete inner part of the furnace in our 3D calculations. However, we used a rather high viscosity ($\mu =0.1\mathrm{Pa}·\mathrm{s}$), which yields a nearly time-independent flow. The main focus is on the influence of the crystal orientation on the thermal stress. More detailed investigations of the melt convection are part of upcoming research.

The importance of the internal radiation was already shown by Xiao and Derby [13] and, later, by Galazka et al. [14]. The additional heat transport compared to an opaque crystal leads to conical interface deep into the melt. Later, Kobayashi et al. investigated the effect of different absorption coefficients on the thermal stress in the Czochralski growth of a LiNbO${}_3$ crystal [15]. Budenkova et al. applied the discrete ordinate method [16] to compute the interface shape including facet formation [17]. For simplicity, they only consider a central facet of the Bi${}_{12}$GeO${}_{20}$ (BGO) crystal. A much more elaborated algorithm with a kinetic operator to compute the interface motion was applied by Weinstein and Miller in a fully 3D calculation [18]. However, the internal radiation was neglected. Nevertheless, the experimentally observed influence of crystal rotation rate and pulling velocity on the facet formation [19] was also found by their computations.

In our study, we assumed a completely transparent crystal. This is a reasonable approach if the β-Ga${}_2$O${}_3$ crystal is electrically insulating, e.g., by doping with Mg. For comparison, calculations of the temperature field have also been performed for an opaque and a semi-transparent crystal. Considering facets requires a different kind of numeric approach, which was not part of this study. The key point of this study is the anisotropy of the thermal conductivity, the thermal expansion, and the elasticity coefficients. Thus, for the same process, thermal stress might be different for different crystal orientation with respect to the pulling direction. Such effect has been studied some time ago for the Czochralski growth of Gd${}_2$SiO${}_5$ [20]. However, in that study, the thermal field was taken from the estimation for a different material. In this paper, we will present global axisymmetric and 3D calculations of the thermal field that enter the computation of the thermal stress.

2. Material Properties

2.1. Properties of β-Ga${}_2$O${}_3$ Crystal

β-Ga${}_2$O${}_3$ crystallizes in the monoclinic space group C2/m (see unit cell in Figure 1). The lattice constants and the monoclinic angle are as observed by experiments and calculations:

	Exp. [21]	Exp. [22]	Calc. [23]
a	12.23	12.23	12.14
b	3.04	3.04	3.14
c	5.80	5.81	5.85
β	103.7	103.85	103.7

β-Ga${}_2$O${}_3$ exhibits perfect cleavage parallel {100} and good cleavage parallel {001}. Therefore, all physical properties reported in the following refer to a Cartesian reference system {${\mathbf{e}}_i$}. The axes ${\mathbf{e}}_i$ are related to the axes of the crystallographic reference system according to ${\mathbf{e}}_1\parallel {\mathbf{a}}^{*}$, ${\mathbf{e}}_2\parallel \mathbf{b}\parallel 2$ and ${\mathbf{e}}_3\parallel \mathbf{c}$, where ${\mathbf{a}}^{*}$ denotes the first basis vector of the reciprocal crystallographic system. The relation between the axes of the crystallographic reference system and the Cartesian ones is depicted on the right side of Figure 1.

Selected isotropic properties of β-Ga${}_2$O${}_3$ are given as follows:

Property	Symbol	Value	Reference
Density	ρ	$5950\mathrm{kg}/{\mathrm{m}}^3$	[24]
Specific heat	$c_p(@300\mathrm{K})$	$560\mathrm{J}·{\mathrm{kg}}^{-1}·{\mathrm{K}}^{-1}$	[6]
	$c_p(@1500\mathrm{K})$	$720\mathrm{J}·{\mathrm{kg}}^{-1}·{\mathrm{K}}^{-1}$	[6]
Melting point temperature	$T_m$	$(2093\pm 20)\mathrm{K}$	${}^1$ [7]
	$T_m$	$2066\mathrm{K}$	${}^2$ [25]
Latent heat	$\Delta H$	$4.6\times {10}^5\mathrm{J}/\mathrm{kg}$	as YAG [13]
Refraction index	$\mathcal{R}$	1.9	[26]
Absorption coefficient	α	$400{\mathrm{m}}^{-1}$	${}^3$ [6]
${}^1$ Measurement during growth experiment with a two-color Raytek pyrometer; ${}^2$ DTA in air; ${}^3$ Crystal with electron concentration $4\times {10}^{16}{\mathrm{cm}}^{-3}$.

There is some discrepancy in the melting point temperature—measurements have been done in different ways. In all our computations, we use $T_m=2080\mathrm{K}$.

2.1.1. Thermal Conductivity

Recently, several measurements of thermal conductivity have been performed. The thermal diffusivity and capacity along [010] direction of Czochralski-grown β-Ga${}_2$O${}_3$ crystals were measured by Netsch–Gerätebau GmbH (Selb, Germany). At $T=1500$ K, they obtained: $\lambda =4.4\mathrm{W}·{\mathrm{m}}^{-1}·{\mathrm{K}}^{-1}$ [6]. Guo et al. measured the thermal conductivity in four crystallographic directions to get the full matrix for the monoclinic system [28]. The measurements were performed in the range of 200 K–495 K and the measured data were approximated by $\lambda \left(T\right)={\lambda }_0{T}^m$. Measurements have been also performed on samples prepared from crystals grown at IKZ by the Czochralski method using the 2ω-method [29]. For room temperature, the results can be summarized as:

		[6]	[28]	[29]
		λ (W·m−1·K−1)	λ (W·m−1·K−1)	λ (W·m−1·K−1)
a	[100]	-	$11\pm 1$	$11\pm 1$
b	[010]	21	$27\pm 2$	$29\pm 2$
c	[001]	-	$15\pm 2$	$21\pm 2$
	$\left[\overline{2}01\right]$	-	$13\pm 1$	-

Because the measurement of Guo et al. [28] is the only one providing the full data set for the monoclinic system, we use their data extrapolated to $T=2000$K:

$$\left({\lambda }_{ij}\right)=\left(\begin{array}{ccc}1.23& 0& 0.38\\ 0& 2.11& 0\\ 0.38& 0& 1.63\end{array}\right)\mathrm{W}·{\mathrm{m}}^{-1}·{\mathrm{K}}^{-1}.$$

2.1.2. Thermal Expansion

The thermal expansion $\left({\beta }_{ij}\right)$ was recently derived from lattice constants determined at different temperatures between 300 K and 700 K employing powder X-ray diffraction [22]. According to the graphics in [22], the slopes of the linear temperature plots above 300 K are $1.85\times {10}^{-5}$Å/K, $1.00\times {10}^{-5}$Å/K, $1.775\times {10}^{-5}$Å/K, and $2\times {10}^{-5\circ }$/K for the lattice constants a, b, c, and the monoclinic angle β, respectively. Using the Cartesian reference system as declared above, the thermal expansion tensor coefficients can be derived from the unit cell parameters as it is written, e.g., in [30]. Extrapolating the data to $T=2000\mathrm{K}$ and computing the tensor coefficients we obtain:

$$\left({\beta }_{ij}\right)=\left(\begin{array}{ccc}-3.4& 0& -9.9\\ 0& 3.4& 0\\ -9.9& 0& 3.1\end{array}\right)\times {10}^{-6}{\mathrm{K}}^{-1}.$$

Thus, there is a large coupling between the directions ${\mathbf{e}}_1$ and ${\mathbf{e}}_3$. In order to understand the large coefficient ${\beta }_{13}$, we inspect the relevant equation from [30]: ${\beta }_{13}$ consists of the algebraic addition of the three slopes $\mathrm{d}a/\mathrm{d}T$, $\mathrm{d}c/\mathrm{d}T$, and $\mathrm{d}\beta /\mathrm{d}T$, each modified by angle functions and room temperature reference values. However, the $\mathrm{d}\beta /\mathrm{d}T$ containing term is around 50 times as large as the $\mathrm{d}a/\mathrm{d}T$ one and around 25 times as large as the $\mathrm{d}c/\mathrm{d}T$ one. In the end, ${\beta }_{13}\approx -0.5\mathrm{d}\beta /\mathrm{d}T$, which leads to the large absolute value of ${\beta }_{13}$ and which also causes the negative sign of ${\beta }_{11}$. Only in the case that $\mathrm{d}a/\mathrm{d}T$ would be more than 3.3 times larger than value measured then ${\beta }_{11}$ would be positive.

Additionally, we studied the thermal expansion of β-Ga${}_2$O${}_3$ employing a more accurate method, namely high-resolution dilatometry. To this end, two rectangular parallelepipeds with edges parallel $\left(90\overline{1}\right)\times \left(010\right)\times \left(001\right)$ and $\left(101\right)\times \left(010\right)\times (\overline{10}.03)$, respectively, and edge lengths between 8 mm and 10 mm were cut from a Czochralski-grown undoped single crystal. At the highest temperature $T=650\mathrm{K}$, where the data was evaluated so far, we have

$$\left({\beta }_{ij}\right)=\left(\begin{array}{ccc}4.7& 0& -0.17\\ 0& 8.3& 0\\ -0.17& 0& 8.5\end{array}\right)\times {10}^{-6}{\mathrm{K}}^{-1}.$$

More details on experiment and data analysis will be reported elsewhere. Both data of Orlandi et al. [22] and our own were used in the computations. The much smaller coupling in our own data is evident, which has a drastic impact on the results of the von Mises stress (see below).

2.1.3. Elastic Stiffness Coefficients

The same samples were used to determine the elastic stiffness coefficients $\left(c_{ij}\right)$ by a combination of the conventional plate-resonance technique and of resonant ultrasonic spectroscopy (RUS) [31]. RUS experiments were performed in the temperature range of 100 K to 1600 K. Currently, the spectra collected close to room temperature have been analysed, yielding (for better visualization, the symmetric lower half of the tensor has been left empty):

$$\left(c_{ij}\right)=\left(\begin{array}{cccccc}23.8& 13.0& 15.2& 0& -0.4& 0\\ & 35.9& 7.8& 0& 0.2& 0\\ & & 34.6& 0& 1.9& 0\\ & & & 4.9& 0& 0.6\\ & & & & 9.1& 0\\ & & & & & 10.7\end{array}\right)\times {10}^{10}\mathrm{N}/{\mathrm{m}}^2.$$

For the purpose of comparison, we refer to our computations employing the software GULP [32], which contains a potential model with Coulomb interaction, Buckingham potential and core shell spring-like interactions. All parameters were taken from Blanco et al. [23], who used a fitting with respect to the experimental lattice geometry, bulk modulus, and static and high-frequency dielectric constants. Please note that thermal energy is not included in this kind of calculation:

$$\left(c_{ij}\right)=\left(\begin{array}{cccccc}28.5& 13.5& 13.5& 0& -1.3& 0\\ & 40.0& 9.0& 0& 0.8& 0\\ & & 37.6& 0& -3.5& 0\\ & & & 5.0& 0& 2.2\\ & & & & 7.3& 0\\ & & & & & 9.3\end{array}\right)\times {10}^{10}\mathrm{N}/{\mathrm{m}}^2.$$

The anisotropy of the computed elasticity tensor is similar to the experimentally observed anisotropy, e.g., $c_{22}>c_{33}>c_{11}$ and $c_{66}>c_{55}>c_{44}$. However, as expected, the calculation in the athermal limit overestimates the longitudinal elastic stiffnesses by about 12%. More details on our experimental work regarding the elastic behavior of β-Ga${}_2$O${}_3$ will be published elsewhere.

2.2. Properties of Ga${}_2$O${}_3$ Melt

Thus far, not much data for the Ga${}_2$O${}_3$ melt have been published. Mass density has been measured some time ago by Dingwell [33] to be $4870\mathrm{kg}·{\mathrm{m}}^{-3}$. From their three data points at temperatures 2073 K, 2123 K, and 2173 K, one can derive a thermal expansion of $\beta =(1.3\pm 0.4)\times {10}^{-5}$K${}^{-1}$. It is in the same range as the one for YAG [13], which we have originally taken for our calculations. Furthermore, we assume a moderate Prandtl number of 4, which is quite common for oxide melts. Because we like to keep the melt convection stable over time, we have chosen a rather large viscosity. The thermal conductivity is then defined by the viscosity and the Prandtl number. In particular, we used the following values:

density	ρ	$6000\mathrm{kg}·{\mathrm{m}}^{-3}$
specific heat	$c_p$	$800\mathrm{J}·{\mathrm{kg}}^{-1}·{\mathrm{K}}^{-1}$
thermal expansion	β	$1.8\times {10}^{-5}{\mathrm{K}}^{-1}$
Prandtl number	Pr	4
dynamic viscosity	μ	$0.1\mathrm{Pa}·\mathrm{s}$
thermal conductivity	λ	20 $\mathrm{W}/\mathrm{m}\mathrm{K}$

We like to make the following note: the thermal conductivity is probably higher than in realitiy as a consequence of the rather large viscosity and the fixed Prandtl number. This will enhance the heat transport in the melt by conduction. The heat transport by advection is reduced because the flow velocity is smaller due to the large viscosity.

3. Computation

A sketch of the furnace geometry is given on the left-hand side of Figure 2. The heat is induced in the iridium crucible by the electric current with a frequency of $f=30\mathrm{kHz}$ in the copper-made coil. This coil outside the insulation is water cooled, and we apply a fixed temperature of 300 K. All computations have been performed for crystal of 2″ diameter with a length of 91 mm plus a seed of 37 mm length. At this growth stage, the melt height was 62 mm in reference to the crucible bottom. The growth rate was 1.5 mm/h and the rotation rate of the crystal was 5 rpm.

In order to treat the anisotropies in the crystal, we have to perform 3D calculations and have to at least include the inner part where the radiative heat exchange with crystal takes place. In particular, we used the inner part including most all insulations as sketched in Figure 2 (middle part). Thermal boundary conditions in terms of temperatures and the distribution of heat sources in the iridium crucible are taken from an axisymmetric calculation using CrysMAS [34]. The final computation of thermal stress is done by applying the software tool comsol. For the latter step, we took only the crystal including seeds, applying the temperatures obtained in the 3D calculation with Ansys-cfx.

For comparison, various calculations have been done with CrysMAS and the software CGsim (http://www.str-soft.com/). Regarding the numerical schemes of the different software tools, all tools used—CrysMAS, CGsim, and Ansys-cfx—apply finite volume schemes to solve the transport equations. The anisotropy of the thermal conductivity can be handled as follows:

• CrysMAS and CGsim: different values for radial and vertical directions;
• Ansys-cfx: different values for x-, y-, and z-directions;
• comsol: full matrix.

In all runs, we solved the Navier–Stokes equation directly and did not use a turbulence model.

Both CrysMAS and CGsim solve the inverse problem of finding the required heater power, however in a different manner. CrysMAS uses the temperature difference at the three-phase junction ($T_m$ versus current temperature) to adjust the heater power. CGsim applies $T_m$ along the entire interface and computes the latent heat required to fulfill this constraint. At the three-phase junction, the difference between the growth rate computed from the required latent heat and the given one is used for adjusting the heater power. In 3D, the computation of the inverse problem is not applicable because iterative relaxation of the temperature field takes too much time. However, the heat source distribution as obtained from CrysMAS is multiplied by a factor, which was adapted by hand from run to run.

As already mentioned, radiation is an important mechanism of heat transport. It is handled differently in the three software packages used. In cavities filled with transparent material, the wall-to-wall radiation is computed. The view factors are computed once before the run, and the computer time needed for the radiative heat transport is negligible. In CrysMAS, such a cavity can contain both gas and solid materials. Thus, in our case, a completely transparent crystal cavity for radiation is given by the crystal/melt interface, the crystal/holder boundary, the melt interface and the inner part of the furnace. In CGsim, only the gas is treated as transparent. The crystal can be treated as semi-transparent, solving the radiative transport by the discrete ordinate method [35]. Transparency can be achieved by choosing small absorption and scattering coefficients (here, we chose $\alpha =1\times {10}^{-5}{\mathrm{m}}^{-1}$ and $\mathcal{S}=1\times {10}^{-5}{\mathrm{m}}^{-1}$). However, the solid material still has a refractive index, which makes the problem different from a transparent cavity of gas and crystal together. Ansys-cfx uses the Monte Carlo method to compute the movements of photons. This is the most general approach but also the most time-consuming.

Different types of grids have been used in the different software tools:

• CrysMAS: Unstructured triangles for the entire geometry—in regions with convection, a structured rectangular grid is used for solving the Navier–Stokes equation (not applied in this paper).
• CGsim: Every domain can have an unstructured triangle or structured rectangular grid. At interfaces of domains, grids do not need to match, but values are interpolated. Only at the crystal melt must interface grids match. For the view factor calculation, an independent grid definition is used.
• Ansys-cfx: We used a fully structured grid. In all computations presented in this paper, we had a total amount of 1,330,756 nodes. In the crystal, the number of grid elements was 140,676. In order to save computational time, for the Monte Carlo ray tracing, the number of grid points can be reduced. Treating the crystal as transparent, we reduced the number of elements for ray tracing to 19,810 in the crystal.

In the axisymmetric runs, the interface shape can be computed. Both CGsim and CrysMAS can adapt the crystal/melt interface shape to the isoline of the melting-point temperature. The node points on the interface are moved in the vertical direction, and the rest of the grid in adjacent domains is adjusted. In CrysMAS, the drawback is the fact that the view factors are not recomputed and, thus, might not fit to the new shape of the interface. This does not held for CGsim because here the discrete ordinate method is used for computing the heat transport by radiation. In 3D, the computation of the interface shape in our problem is currently beyond our capabilities. In fact, in order to include facet formation, we need to use a different program than Ansys-cfx (see [18]). As we know the crystal shape from the experiment, we used this geometry for our 3D computations. At the crystal/melt boundary, we apply an additional heat flux ($1333\mathrm{W}/{\mathrm{m}}^2$) for latent heat release, corresponding to a growth rate of 1.5 mm/h.

The computation of the 3D thermal stresses were performed by using the software comsol Multiphysics (v. 5.2, COMSOL Multiphysics GmbH, Göttingen, Germany), i.e., by using the modules for linear elasticity and for heat transfer in solids, which were coupled by temperature and thermal expansion. The temperature data for all outer boundaries (crystal–melt interface, surface of cylinder, cone and seed) were taken from global modelling of the inner part of the crystal growth equipment (middle of Figure 2) and set as a Dirichlet-type condition. Hence, the temperature distribution within the crystal was recomputed by using the same thermal conductivity as in the global modelling. For the mechanics, the displacements are set to zero at the ring boundary of the seed end wall where the crystal is fixed to the holder. The reference temperature is the temperature at the top of the seed. The Comsol procedure uses the finite element method, and we selected elements with quadratic shape function. The tensorial material properties such as thermal conductivity, thermal expansion coefficients and elasticity coefficients were treated by rotation operation in order to meet the crystal orientation.

The overall strategy is the following: in the first step (left side of Figure 2), the temperature field is computed by a global simulation within a real growth furnace and chamber. For the calculations, measured temperature values at different points of the furnace were used as monitoring points.

In the first attempt, we performed a calculation using CrysMAS without melt convection and with a uniform thermal conductivity $\lambda =4.4\mathrm{W}·{\mathrm{K}}^{-1}·{\mathrm{m}}^{-1}$. From this calculation, we extract the temperatures at the boundaries of the 3D domain (middle of Figure 2). At the insulation (domain in magenta), the data from CrysMAS is directly used (interpolation to the grid used with Ansys-cfx). At the (virtual) outer boundary of the gas domain (blue), we approximate the temperature profile by $1239.0\mathrm{K}-232.1\mathrm{K}log[80{\mathrm{m}}^{-1}z-46.44]$ and $900\mathrm{K}-1850\mathrm{K}50{\mathrm{m}}^{-1}r$ at the side and the top, respectively. r is the radius and z the vertical coordinate measured from the bottom of the inhousing. Then, we performed the 3D calculations as described in the next section. The 3D calculations of the thermal field were very time-consuming. Therefore, we did not update the boundary conditions by axisymmetric computations including melt convection, as we recently performed with CGsim. In the end, we used the temperature distribution at the crystal and seed boundary as input for the computation with Comsol. The thermal computations was repeated and afterwards the thermal stress was calculated.

4. Results

The main question we address is the influence of the crystal orientation on the von Mises stress. We used two orientations of the crystal: in Case 1, the growth direction is the b-axis ((010) plane in the growth direction) as in all of our growth experiments by the Czochralski method; in Case 2, the growth direction is the a*-axis ((100) plane in the growth direction). This orientation has not been used in the Czochralski growth but was observed in the vertical Bridgman growth [25]. The two orientations are shown in Figure 3. The thermal conductivities in $\mathrm{W}·{\mathrm{m}}^{-1}·{\mathrm{K}}^{-1}$ are:

	Vertical	Horizontal
	${\mathit{\lambda }}_{\mathit{z}}$	${\mathit{\lambda }}_{\mathit{x}}$	${\mathit{\lambda }}_{\mathit{y}}$	${\mathit{\lambda }}_{\mathit{r}}$
Case 1	2.11	1.07	1.79	1.43
Case 2	1.23	2.61	1.63	1.87

${\lambda }_i$ ($i=x,z,y$) has been computed from the tensor $\left({\lambda }_{ij}\right)$ by appropriate rotations. ${\lambda }_x$ and ${\lambda }_y$ are used in the 3D computations. ${\lambda }_r$ is used in the axisymmetric calculations and was computed as the average of ${\lambda }_x$ and ${\lambda }_y$.

In the first step, we computed the temperature field. We compared the runs with different software tools (see Figure 4). Please note that, as mentioned above, the 3D calculation was performed with temperature boundary conditions obtained from a CrysMAS run without melt convection. Therefore, the temperature at the top of the domain for Ansys-cfx (see Figure 2) is about 100 K higher than in the case with melt convection. Thus, the temperature at the top of the seed is also higher in the case of the 3D calculation (about 20 K). The other point to be mentioned is that, in CGsim, the radiation within the crystal is computed by the discrete ordinate method with absorption and scattering coefficients of $\alpha =1\times {10}^{-5}{\mathrm{m}}^{-1}$ and $\mathcal{S}=1\times {10}^{-5}{\mathrm{m}}^{-1}$, respectively.

In Figure 5, the temperature profiles in melt, crystal and seed are shown for the two different crystal orientations. Only small differences occur, which can be seen in more detail in the temperature profile along the axis (plot on the right-hand side). In the melt, the advective heat transport can be clearly seen. Using the high viscosity ($\mu =0.1\mathrm{Pa}·\mathrm{s}$), one large convection roll is established, transporting the heat from the lower part, where the induced heat in the crucible is the highest, towards the melt surface and the crystal–melt interface. We considered buoyancy convection and a crystal rotation of 5 rpm. The resulting convection causes a nearly uniform temperature on the crystal–melt interface. As already mentioned, the interface shape was fixed in the 3D calculations and, thus, the temperature at the interface is a result of the computations. The temperature variation is about 3 K (see also Figure 7).

A more detailed analysis of the temperature in crystal and seed for different runs is given in Figure 6. Four calculations have been performed with CGsim (figures with grey background). In Figure 6A, the result of the run with $\alpha =1\times {10}^{-5}{\mathrm{m}}^{-1}$ and $\mathcal{S}=1\times {10}^{-5}{\mathrm{m}}^{-1}$ is shown (the refraction index was set to $\mathcal{R}=1.4$.), i.e., for a quasi-transparent crystal. In complete contrast, Figure 6B presents the result for an opaque crystal. Because there is no heat transport by radiation, all of the heat from the interface has to be removed by heat conduction, resulting in a much higher temperature gradient than in (A). For a semi-transparent crystal with $\alpha =400{\mathrm{m}}^{-1}$ and $\mathcal{S}=1\times {10}^{-3}{\mathrm{m}}^{-1}$ (C), the gradient is still higher than in case (A) but strongly reduced in comparison to (B). As already noticed from Figure 5, there is not so much difference between Case 1 and Case 2, i.e., for the two different crystal orientations (compare Figure 6C, D). Going to a 3D computation for the semi-transparent crystal, there is the difference that, in the axisymmetric calculation by CGsim, the temperature at the interface is fixed to $T_m=2080\mathrm{K}$, whereas, in the 3D calculation, the latent heat release at the interface is fixed. Therefore, in the latter case (Figure 6E), the maximum temperature is 2110 K, with a considerable temperature gradient in the lower part of the cone. In (F), there is again the case with a transparent crystal. The total temperature difference in crystal and seed is 140 K as compared to 200 K in the case with a semi-transparent crystal and seed (E). Thus, the conductive heat transport is about 1.4 times larger for semi-transparent crystal and seed.

In the 3D calculation, the interface shape is fixed as taken from the typical shape in the growth experiments. This means that there is a temperature variation along the interface. When melt convection with $\mu =0.1\mathrm{Pa}·\mathrm{s}$ is applied, the temperature difference is about 3 K, as indicated by the colours along the interface in Figure 7. The interface shape obtained in the axisymmetric calculation using CGsim is given by the black line in Figure 7. It is much deeper compared with the shape used in the 3D calculation. It is also much deeper than the isoline of $T=2081\mathrm{K}$ (black dotted line) as obtained in the 3D computation. One has to keep in mind that the interface shape has an influence on the convection, and this is vice versa on the heat transport in the melt. However, without convection, the interface is much flatter as computed by CGsim (dotted line in magenta). This underlines the necessity to include melt convection to recover the interface shape as obtained by the real growth process.

The thermal stress was computed only for the temperature fields observed by the 3D calculations employing Ansys-cfx, as indicated in the sketch of Figure 2. We used both the thermal expansion coefficients taken from [22] and from our measurements. The results are completely different due to the strong coupling between [001] and [100] direction, in the case of the coefficients from [22]. Figure 8 shows the results using the thermal expansion coefficients from [22]: in Case 1, the crystal is shrinking in the radial direction and extending in the growth direction yielding in a high von Mises stress in the central part of the interface. In Case 2, the strong coupling via ${\beta }_{13}$ leads to a skew growth of the crystal, which is visible by the magnified geometry deformation on the right-hand side of Figure 8. The influence of the strong coupling is also visible in the high anisotropy of the stress distribution. As a result, the maximum of von Mises stress is higher in Case 1 (about 8.7 MPa) as compared to Case 2 (about 5.1 MPa).

The field distribution of displacement and stress is very different when using the thermal expansion coefficient measured by us (see Figure 9): the coupling via the off-diagonal thermal expansion is rather weak; therefore, in both cases, the crystal geometry nearly maintains its axisymmetry. However, a significant anisotropy in the von Mises stress can be observed. Firstly, we discuss Case 1. the blue stripe of minimal stress at the bottom of the crystal is in the a* direction. The main cleavage plane, the (100) plane, is perpendicular to this direction, and, thus, feels a significant change of stress from its boundary towards its centre.

In Case 2, we have the b- and c-axis in the (100)-plane perpendicular to the growth direction. On the crystal–melt interface, we observe high stress in one direction going up to ≈6.4 $\mathrm{MPa}$ close to the periphery. This is the direction of the c-axis. In the b-axis, the variation from the tip (low) towards the rim (higher) is much smaller.

5. Conclusions

For the first time, we show 3D computation of the thermal stress in the Czochralski growth of β-Ga${}_2$O${}_2$ single crystals. The computations are based on the thermal field obtained by the sequence of an axisymmetric calculation of the entire furnace and 3D calculations of the inner part. The key point in the thermal calculations was the use of the anisotropic thermal conductivity obtained by measurements. As an important ingredient for the computation of thermal stress, we measured the thermal expansion coefficients and the elasticity coefficients. We observed a significant difference of thermal expansion coefficients with respect to those obtained by Orlandi et al. [22]. The two different data sets result in significantly different distributions of stress in the crystal. We also observed a remarkable difference in the stress distribution for the crystal orientations investigated. In the configurations investigated, the maximum von Mises stress is lower when growing in the direction of the b-axis (i.e., the (010) plane is normal to the pulling direction) than when growing in the a* direction (i.e., the (100) plane is normal to the pulling direction).

We are now moving to a more systematic investigation of the influence of crystal orientation on the stress. Furthermore, there are several challenges for further steps. The semi-transparency and its influence on stress has to be studied in more detail. The real interface shape has to be used for the 3D calculations of stress. Moreover, the facets have to be taken into account, and the melt convection has to be analysed going towards a more realistic value of viscosity.