Numerical Modelling of the Czochralski Growth of β-Ga 2 O 3

Our numerical modelling of the Czochralski growth of single crystalline β-Ga2O3 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 β-Ga2O3 crystal. In order to perform the stress analysis, we measured the thermal expansion coefficients and the elastic stiffness coefficients in two samples of a β-Ga2O3 crystal grown at IKZ. Additionally, we analyse published data of β-Ga2O3 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.


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 (Pr = c p µ/λ, 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 (µ = 0.1 Pa • 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.

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: β-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 {e i }.The axes e i are related to the axes of the crystallographic reference system according to e 1 a * , e 2 b 2 and e 3 c, where 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.
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 K.

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: 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 λ(T) = λ 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: 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:

Thermal Expansion
The thermal expansion (β ij ) 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 × 10 −5 Å/K, 1.00 × 10 −5 Å/K, 1.775 × 10 −5 Å/K, and 2 × 10 −5 • /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 K and computing the tensor coefficients we obtain: Thus, there is a large coupling between the directions e 1 and e 3 .In order to understand the large coefficient β 13 , we inspect the relevant equation from [30]: β 13 consists of the algebraic addition of the three slopes da/dT, dc/dT, and dβ/dT, each modified by angle functions and room temperature reference values.However, the dβ/dT containing term is around 50 times as large as the da/dT one and around 25 times as large as the dc/dT one.In the end, β 13 ≈ −0.5 dβ/dT, which leads to the large absolute value of β 13 and which also causes the negative sign of β 11 .Only in the case that da/dT would be more than 3.3 times larger than value measured then β 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 (90 1) × (010) × (001) and (101) × (010) × ( 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 K, where the data was evaluated so far, we have 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).

Elastic Stiffness Coefficients
The same samples were used to determine the elastic stiffness coefficients (c ij ) 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): 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 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: 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.

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 kg • m −3 .From their three data points at temperatures 2073 K, 2123 K, and 2173 K, one can derive a thermal expansion of β = (1.3 ± 0.4) × 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: 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.

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 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:
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 α = 1 × 10 −5 m −1 and S = 1 × 10 −5 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 W/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 λ = 4.4 W • K −1 • 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.0K − 232.1 K log[80 m −1 z − 46.44] and 900 K − 1850 K 50 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.

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 W • m −1 λ i (i = x, z, y) has been computed from the tensor (λ ij ) by appropriate rotations.λ x and λ y are used in the 3D computations.λ r is used in the axisymmetric calculations and was computed as the average of λ x and λ 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 The curve in magenta is from the Ansys-cfx run (3D calculation).The temperature at the top is higher because we used the boundary from a CrysMAS run without melt convection.The red and blue curves are the results of the runs with CGsim and CrysMAS, 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 (µ = 0.1 Pa • 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 α = 1 × 10 −5 m −1 and S = 1 × 10 −5 m −1 is shown (the refraction index was set to 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 α = 400 m −1 and S = 1 × 10 −3 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 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 µ = 0.1 Pa • 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 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 β 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 band 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 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.

Figure 1 .
Figure 1.Unit cell of β-Ga 2 O 3 (left) and standard coordinate system (right).Ga atoms in green (large) and oxygen ones in red (small).The figure is drawn with VESTA [27].

Figure 2 .
Figure 2. Sketch of the geometry and the computational steps.Below: the software used.

Figure 3 .
Figure 3. Two cases of different crystal orientations.

Figure 4 .
Figure 4. Temperature profile for Case 1 along the axis for runs with three different software tools.The curve in magenta is from the Ansys-cfx run (3D calculation).The temperature at the top is higher because we used the boundary from a CrysMAS run without melt convection.The red and blue curves are the results of the runs with CGsim and CrysMAS, respectively.

Figure 5 .
Figure5.Temperature fields in melt, crystal, and seed for the two cases of different growth directions, as observed by 3D calculations with Ansys-cfx.The temperature profile at the center of the crystal is plotted on the right-hand side.Please note that the temperature scales for melt and crystal are given seperately to ensure a more detailed view.Because the temperature varies slightly along the crystal-melt interface, the minimum temperature in the melt and the maximum one in crystal are not identical but overlap.

Figure 6 .
Figure 6.Temperature profiles in the crystal for different runs (isolines every 10 K): (A) CGsim computation with α = 1 × 10 −5 m −1 for Case 1; (B) CGsim computation with opaque crystal for Case 1; (C) CGsim computation with α = 400 m −1 for Case 1; (D) CGsim computation with α = 400 m −1 for Case 2; (E) CFX computation with α = 400 m −1 for Case 2; and (F) CFX computation with transparent crystal for Case 2. Please note that thermal conductivity is unique in the horizontal direction for the CGsim calculations, whereas, in the 3D calculation, it depends on the direction.

Figure 7 .
Figure 7.Comparison of interface shapes.The thick line with colours indicates the preset interface for the 3D calculations, and the colours exhibit the temperatures along this line from the computation for Case 1.The dotted line is the isoline of T = 2081 K.We used axisymmetric calculations with CGsim to compute the crystal-melt interface.The black line is the interface for the same case as the 3D one.The dotted line in magenta is the interface for the case without melt convection.

Figure 8 .Figure 9 .
Figure 8. Von Mises stress for Case 1 (left) and Case 2 (right).The thermal expansion coefficients of [22] have been used.The thin black lines indicate the orginal size of the crystal.The displacements due to the different temperatures are visualized by the surface with the colours of the von Mises stress (deformation magnification factor, (left): 100; (right): 20).