Macro- and Microsimulations for a Sublimation Growth of SiC Single Crystals

The numerous technical applications in electronic and optoelectronic devices, such as lasers, diodes, and sensors, demand high-quality silicon carbide (SiC) bulk single crystal for industrial applications. We consider an SiC crystal growth process by physical vapor transport (PVT), called modified Lely method. We deal with a model for the micro- and macroscales of the sublimation processes within the growth apparatus. The macroscopic model is based on the heat equation with heat sources due to induction heating and nonlocal interface conditions, representing the heat transfer by radiation. The microscopic model is based on the quantum interatomic potential and is computed with molecular dynamics. We study the temperature evolution in the apparatus and reflect the growth behavior of the microscopic model. We present results of some numerical simulations of 
the micro- and macromodels of our growth apparatus.


Introduction
The motivation for this study comes from the technical demand to simulate a crystal growth apparatus for SiC single crystals.The single crystals are used as a high-valued and expensive material for optoelectronics and electronics cf. 1 .We concentrate on a deterministic model for simulating crystal growth; alternative models are discussed with comprehensive probabilistic modeling see 2 .
The silicon carbide SiC bulk single crystals are produced by a growth process through physical vapor transport PVT , called modified Lely method.The modeling for the thermal processes within the growth apparatus is done in 3, 4 .In this paper, we propose one step more in the modeling of the macroscopic and microscopic parts.The idea is to exchange results from the macroscopic to the microscopic scale to obtain a feedback to control the growth process of the SiC bulk.Here the benefits are an acceleration of solving interactive growth processes of the crystal with their underlying temperature in the apparatus.Using only standard codes, which are decoupled, a simple parameter exchange of temperature and pressure in the deposition region cannot resolve the growth problem accurately.We propose a first framework of a combined model, which is based on the authors' knowledge of a novel work and a first approach to a coupled solver method.

Macroscopic Model: Heat-Flux
In the following, we discuss the macroscopic model, which is based on continuum equations for the heat-flux.

Mathematical Model
The underlying equations of the model are given as follows.
a In this work, we assume that the temperature evolution inside the gas region Ω g can be approximated by considering the gas as pure argon Ar .The reduced heat equation is where T is the temperature, t is the time, and U g is the internal energy of the argon gas.The parameters are given as ρ g being the density of the argon gas, κ g being the thermal conductivity, z Ar being the configuration number, and R Ar being the gas constant for argon.b The temperature evolution inside the region of solid materials Ω s e.g., inside the silicon carbide crystal, silicon carbide powder, graphite, and graphite insulation is described by the heat equation where ρ s is the density of the solid material, U s is the internal energy, κ s is the thermal conductivity, and c s is the specific heat.The equations hold in the domains of the respective materials and are coupled by interface conditions, for example, requiring the continuity for the temperature and for the normal components of the heat flux on the interfaces between opaque solid materials.On the boundary of the gas domain, that is, on the interface between the solid material and the gas domain, we consider the interface condition where n g is the normal vector of the gas domain, R is the radiosity, and J is the irradiosity.The irradiosity is determined by integrating R along the whole boundary of the gas domain cf. 5 .Moreover, we have where E is the radiation, J ref is the reflexed radiation, is the emissivity, and σ is the Boltzmann radiation constant.
The density of the heat source induced by the induction heating is determined by solving Maxwell's equations.We deal with these equations under the simplifying assumption of an axisymmetric geometry, axisymmetric electromagnetic fields, and a sinusoidal time dependence of the involved electromagnetic quantities, following 6 .The considered system and its derivation can be found in 3, 4, 7 .
In this paper, we focus on the discretization and material properties, which are important for realistic simulations.Our underlying software tool WIAS-HiTNIHS cf. 4 allows us a flexibility in the grid generation and for the material parameters.
In the next section, we describe the used discretization.

Discretization
For the discretization of the heat equation diffusion equation , we apply the implicit Euler method in time and the finite volume method for the space discretization cf. 3, 4, 8 .We consider a partition T ω i i∈I of Ω such that, for m ∈ {s, g} with s solid, g gas and i ∈ I, ω m,i : ω i ∩ Ω m defines either a void subset or a nonvoid, connected, and open polyhedral subset of Ω.By integrating the corresponding heat equation 2.1 or 2.3 over ω m,i , we derive the following nonlinear equations for the temperature variables, where the time interval is Δt n 1 t n 1 − t n .The temperature is given as T n 1 T t n 1 , x , where x represents cylindrical coordinates.For the right-hand sides, we demand f s : f ≥ 0 and f g 0.
More details of the discretization and of dealing with the interface conditions are presented in 3, 4, 9, 10 .
In the next section, the properties of the materials in the crystal growth apparatus are described.

Material Properties
For the technical realization of the apparatus, we implement the axisymmetric geometry given in 11 , which is presented in Figure 1.Furthermore, the properties of the materials are specified in 3, 9, 12 .

Mathematical Problems in Engineering
Within the following specific material functions and parameters for the processes, the thermal conductivity κ is given in W/ m K , the electrical conductivity σ c is given in 1/ Ohm m , the mass density ρ is given in kg/m 3 , the specific heat c sp is given in J/ K kg , the temperature T is given in K, and the relative gas constant R Ar is given in J/ K kg .Further, the emissivity and relative magnetic permeability μ are given dimensionless.
The functions are programmed in our flexible software package WIAS-HiTNIHS.
In the next section, we discuss the microscopic model.The growth apparatus' dimensions: r min 0, r max 8.4 cm, z min 0, z max 25.0 cm, the coil rings' dimensions: r min 4.2 cm, r max 5.2 cm, z min 0, z max 14.0 cm.

Coupling Method for Macroscopic and Microscopic Models: Operator Splitting
Often simple coupling via the parameters e.g., target-temperature and growth velocity of the bulk is enough for the problem.
Here we propose a new idea of coupling the model equations together, on the one hand the heat equations and on the other hand the kinetic equations for molecules.
For a first idea, we deal with abstract operators, which include the heat-and the kinetics equations.
Using our two standard codes of the macro-and micromodels, we could implement a coupled model, by a so-called iterative operator-splitting method.Such a proposed method couples the two physical processes of the thermal situation in the growth apparatus and their important geometrical differences at the deposition layer with the kinetic molecular model.The benefits are a numerical algorithm, that exchanged the underlying operators of the thermal situation and the kinetic molecular situation, which are computed by each software code independently and coupled via an iterative solver step; see a detailed coupling analysis in 13 .
In the following algorithm, an iteration method is proposed, with fixed splitting discretization step-size τ.
Due to the underlying multiscale problem of kinetics and heat processes, we have to solve fine time scales of kinetic equations and coarse time scales for heat equations.On a time interval t n , t n 1 that is sufficiently small to yield accurate kinetics, we solve the following subproblems consecutively for i 0, 2, . . ., 2m cf.14, 15 : where we assume that the operator A has a large time scale macroscopic model and B has a small time scale microscopic model .Further c 0 t n c n , c −1 0, and c n are initialization and starting conditions.
In the following, we give an overview to the accuracy of the method, which is given in the convergence and the rate of the convergence.
Theorem 2.1.Let us consider the abstract Cauchy problem in a Banach space X: where A, B, A B : X → X are given linear operators being generators of the C 0 -semigroup and c 0 ∈ X is a given element.Then the iteration process 2.11 -2.12 is convergent.The rate of convergence is of higher order and given as O τ 2m n , where the iterations are i 1, 3, . . ., 2m 1.
The proof is given in 15 .
In the next subsection, we present the methods for the microscopic model.

Microscopic Model: Quantum Chemical Molecular Dynamics (QM/MD) of SiC Condensation (Methodology)
The density-functional tight-binding DFTB method is employed as the quantum interatomic potential in our molecular dynamics MD simulations, using atomic and diatomic parameters obtained from density functional theory; see 16 .DFTB is an approximate density functional theory method based on the tight binding approach, and utilizes an optimized minimal LCAO Slater-type all-valence basis set in combination with a two-center approximation for Hamiltonian matrix elements.Parameter sets for Si-Si and Si-C were taken from 17 .Energies and gradients are evaluated direct on the fly during the dynamic simulation.As in our previous simulations of carbon cap 18 and subsequent nanotube formation 19 on the C-and Si-faces of SiC 000-1 surfaces during sublimation evaporation, we have not included charge-or spin-polarization in the present work.Further, we will consider in a next model electrokinetic effect on heat transfer in parallel-plate microchannels, hydrodynamic focusing effects, and nanoeffect as done in 20-23 .
For time propagation we employed a velocity Verlet integrator with a time step of 1.209 fs 50 atomic units and used a Nose-Hoover chain thermostat to generate a canonical ensemble for target temperature T t ; see 24 .The thermostat was employed uniformly in the reaction system.Regarding the atomistic structure of the employed surface model systems, we have chosen the C-face of the same square SiC 000-1 slab unit cell as in our previous study, 19 consisting of two SiC layers terminated by hydrogen atoms to mimic bulk effect in the direction away from the surface.Periodic boundary conditions were employed with a unit cell size of 1000 Å in the direction perpendicular to the surface and 16.0 Å and 15.4 Å in the other two surface-directions to achieve two-dimensional slab periodicity.The geometry optimized structure of this surface model is shown in Figure 2.
During MD simulations, the movements of hydrogen terminating atoms were frozen.Using such an approach, we have effectively introduced a steep temperature gradient from the deepest bulk-side SiC layer to the atoms lying above on the surface.The slab model was then annealed at T t 2000 K for 1.20 picoseconds, and 3 structures were selected as initial starting geometries at t 0.60 picosecond trajectory A50 , t 0.72 picosecond trajectory A60 , and t 0.86 picosecond trajectory A80 .In the vicinity of the surface, 10 SiC molecules were randomly distributed in the gas phase.Since these molecules are nonbonded to the surface, they are subsequently thermostated at T t .Gas phase molecules approaching the surface will experience immediate cooling, which will drive the condensation process during these simulations.
In the microscopic model, we can derive the growth rate v of the seed surface in dependence on temperature and pressure.Based on these growth rates, we can adapt the geometry for the macroscopic model.Such modification helps to give more accurate temperature differences in the macroscopic model and understand the growth process.
In the next section, we present results of our numerical experiments.

Numerical Experiments
We present in the following our macro-and microscopic simulations, where the microscopic simulations take into account the target temperature of the macroscopic model.

Macroscopic Model: Simulation of the Temperature Field in the Apparatus
For the numerical results, we apply the parameter functions in Section 2.3.We consider the geometry shown in Figure 1, using a constant total input power of 10 kW cf. on the software package pdelib cf. 25 which uses the sparse matrix solver PARDISO cf. 26 .We compute the coupled system consisting of the heat equations and Maxwell's equations.For the growth process, the temperature difference T ss T r source , z source − T r seed , z seed with the coordinates r source , z source 0, 0.143 and r seed , z seed 0, 0.158 , corresponding to the points T source and T seed in Figure 1 is crucial.On the other hand, in the physical growth experiments, usually only the temperatures T r bottom , z bottom and T r top , z top with the coordinates r bottom , z bottom 0, 0.028 and r top , z top 0, 0.173 , corresponding to the points T bottom and T top in Figure 1 are measurable and their difference T bt T r bottom , z bottom −T r top , z top is often used as an indicator for T ss .In Figure 3, we present the temperature differences T ss and T bt .As a result of our computations, the temperature difference T bt can only restrictively be used as an indicator for the temperature difference T ss cf. the discussions in 5, 9 .
The further computations are based on the stationary case, dealing with 2.1 by discarding the terms with a time derivative.For this case, the results are virtually equal to the one in the transient case with t > 15000 seconds.For the stationary results, we focus on the error analysis for the space dimension by applying the grid refinement.The solutions for the heat equation are computed at the points T r bottom , z bottom and T r top , z top for successive grids.For the error analysis, we apply the following error differences: abs where T j r, z and T j 1 r, z are solutions evaluated at the point r, z which has been computed using the grids j and j 1, respectively.The elements of the grid j 1 are approximately 1/4 of the elements of the grid j.The results are presented in Table 1.
The result of the refinement indicates the reduction of the absolute difference as it is demanded for the convergence of the discretization method.The method is stabilized in the presented refinement by reducing the differences.In Figure 4, the temperature field is presented for the stationary case.The temperature increases from the bottom to the middle of the graphite pot, and decreases from the middle to the top of the graphite pot.

Microscopic Model: Atomistic QM/MD Simulations of SiC Condensation on the C-face of Si(000-1)
The total time of the three condensation simulations was 24.02 picoseconds.This is admittedly a time too short for the study of crystal growth, which would ideally require annealing simulations on the order of several 100 nanoseconds, but this study is focusing on the initial stages of SiC aggregation and tries to identify key features in the condensation process.As such, this is at present rather preliminary study exploring the applicability of QM/MD simulations for SiC crystal growth.We have first concentrated on the polar Csurface of SiC 0001 since it has a maximum of dangling bonds with highest reactivity.The Si-face and other nonpolar surfaces are much less reactive; see 27 .Since our seed crystal surface slab model contains only two SiC layers, we are also unable to address the issue of polymorphism at present, although it should be noted that our model system rather resembles the cubic 3C than the hexagonal polytypes.T t was chosen as 2000 K for all simulations, and representative snapshots from trajectory A50 are given in Figure 5.We find that under the present conditions with a relatively high density of SiC gas molecules, several of them attach very quickly to the surface 2 after 0.10 picosecond .Also, SiC molecules can react with each other to form dimers, preferably with C-C bonds.Eventually, an average of 5.3 SiC molecules become attached to the surface in the three simulations, with the other molecules being lost to the vacuum layer.
Once attached, the Si atoms on the surface prove to be highly mobile, as their bond radius is larger than the case of carbon, and the binding energies are lower 18 .The carbon atoms on the surface tend to form C 2 units, and behave similar to "wobbling C 2 " entities that we had observed for high-temperature simulations of pure carbon; see 28 .It seems from our simulations at this stage that the system tries to reach an equilibrium with a constant number of C-C in the new layers, and that the Si atoms are more isolated, becoming occasionally attacked by a C2 dimer.In particular, C 2 units are oriented mainly perpendicular to the surface, while the more visible Si 2 dimers do not show such an alignment preference.The surface itself retains the structure of alternating Si-C units.A new layer of Si-C units is being deposited with a somewhat inhomogeneous structure containing C 2 and Si 2 units at first, and gradually becoming more homogeneous due to annealing.

Summary
We have presented a model for the heat transport inside a technical apparatus for crystal growth of SiC single crystals.We introduce the heat equation and the radiation of the apparatus and the coupled situation of the different materials.The equations are discretized by the finite volume method and the complex material functions are embedded in this method.Transient and stationary results are presented leading to some information about the processes within the technical apparatus.We present numerical results for the stationary case to support the accuracy of our solutions.We also presented atomistic quantum chemical molecular dynamics QM/MD simulations based on the density-functional tight-binding DFTB method for initial reactions of gaseous SiC on the polar C-face of SiC 000-1 .In our future work, we concentrate on further implementations and numerical methods for a crystal growth model and use kinetic data obtained from more accurate microscopic model simulations in the simulation of the heat transport.Once longer and a larger number of trajectories are obtained in our microsimulations, it will be possible to deduct an accurate QM/MD-based estimate for the bulk growth, in dependence on the temperature to our macrosimulations.This data will then enter the iterative solution of the heat and kinetics equations of the coupled macroscopic and microscopic models.

Figure 2 :
Figure 2: Optimized geometry of the C-face of the 000-1 SiC surface as initial starting point for QM/MD simulations.Blue spheres correspond to silicon atoms, purple spheres correspond to carbon atoms, and white spheres correspond to hydrogen atoms terminating the slab model in bulk direction.The model is the unit cell used in periodic boundary calculations with infinite surface extension.

Figure 3 :
Figure 3: Transient results for the temperature differences T bt and T ss .

TFigure 4 :
Figure 4: Temperature field for the apparatus simulated for the stationary case with 23017 nodes.

MathematicalFigure 5 :
Figure 5: Simulation of the addition of 10 SiC atoms on the C-face of the 000-1 SiC surface from Figure 2. Blue spheres correspond to silicon atoms, purple spheres correspond to carbon atoms, and white spheres correspond to hydrogen atoms terminating the slab model in bulk direction.Times are given in picoseconds ps , indicating that the moment the snapshots were taken during the dynamics simulations.

Table 1 :
Computations on different grids for the errors analysis with absolute differences cf.4.1 .