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.

1. 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.

2. Macroscopic Model: Heat-Flux

In the following, we discuss the macroscopic model, which is based on continuum equations for the heat-flux.

2.1. 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 can be approximated by considering the gas as pure argon (Ar). The reduced heat equation iswhere is the temperature, is the time, and is the internal energy of the argon gas. The parameters are given as being the density of the argon gas, being the thermal conductivity, being the configuration number, and being the gas constant for argon.

(b) The temperature evolution inside the region of solid materials (e.g., inside the silicon carbide crystal, silicon carbide powder, graphite, and graphite insulation) is described by the heat equationwhere is the density of the solid material, is the internal energy, is the thermal conductivity, and 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 conditionwhere is the normal vector of the gas domain, is the radiosity, and is the irradiosity. The irradiosity is determined by integrating along the whole boundary of the gas domain (cf. [5]). Moreover, we havewhere is the radiation, 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.

2.2. 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 of such that, for (with solid, gas) and 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 we derive the following nonlinear equations for the temperature variables,where the time interval is The temperature is given as where represents cylindrical coordinates. For the right-hand sides, we demand and

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.

2.3. 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].

Within the following specific material functions and parameters for the processes, the thermal conductivity is given in W/(m K), the electrical conductivity is given in 1/(Ohm m), the mass density is given in the specific heat is given in J/(K kg), the temperature is given in and the relative gas constant is given in J/(K kg). Further, the emissivity and relative magnetic permeability are given dimensionless.

(i) For the gas phase (argon), we have where

(ii) For graphite felt insulation, we have where

(iii) For the graphite, we have where

(iv) For the SiC crystal, we have

(v) For the SiC powder, we have

The functions are programmed in our flexible software package WIAS-HiTNIHS.

In the next section, we discuss the microscopic model.

2.4. 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 that is sufficiently small to yield accurate kinetics, we solve the following subproblems consecutively for (cf. [14, 15]): where we assume that the operator has a large time scale (macroscopic model) and has a small time scale (microscopic model). Further and 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 where are given linear operators being generators of the -semigroup and 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 where the iterations are

The proof is given in [15].

In the next subsection, we present the methods for the microscopic model.

3. 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 [2023].

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 ; 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 for 1.20 picoseconds, and 3 structures were selected as initial starting geometries at picosecond (trajectory A50), picosecond (trajectory A60), and 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 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 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.

4. 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.

4.1. 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 (cf. [11]). The numerical experiments are performed using the software WIAS-HiTNIHS (cf. [4]) based 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 (with the coordinates and corresponding to the points and in Figure 1) is crucial. On the other hand, in the physical growth experiments, usually only the temperatures and (with the coordinates and corresponding to the points and in Figure 1) are measurable and their difference is often used as an indicator for In Figure 3, we present the temperature differences and As a result of our computations, the temperature difference can only restrictively be used as an indicator for the temperature difference (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 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 and for successive grids. For the error analysis, we apply the following error differences:where and are solutions evaluated at the point which has been computed using the grids and respectively. The elements of the grid are approximately of the elements of the grid 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.

4.2. 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 a preliminary study exploring the applicability of QM/MD simulations for SiC crystal growth. We have first concentrated on the polar C-surface 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.

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 units, and behave similar to “wobbling ” 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, units are oriented mainly perpendicular to the surface, while the more visible dimers do not show such an alignment preference. The surface itself retains the structure of alternating –C units. A new layer of –C units is being deposited with a somewhat inhomogeneous structure containing and units at first, and gradually becoming more homogeneous due to annealing.

5. 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.