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

Abstract

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 Ω_g can be approximated by considering the gas as pure argon (Ar). The reduced heat equation is𝜌_g𝜕_𝑡𝑈_g−∇⋅(𝜅_g𝑈∇𝑇)=0,_g=𝑧_Ar𝑅_Ar𝑇, (2.1)(2.2)where 𝑇 is the temperature, 𝑡 is the time, and 𝑈_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, 𝑧_Ar being the configuration number, and 𝑅_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𝜌_s𝜕_𝑡𝑈_s−∇⋅(𝜅_s𝑈∇𝑇)=𝑓,_s=∫^𝑇₀𝑐_s(𝑆)𝑑𝑆, (2.3)(2.4)where 𝜌_s is the density of the solid material, 𝑈_s is the internal energy, 𝜅_s is the thermal conductivity, and 𝑐_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𝜅_g∇𝑇⋅𝐧_g+𝑅−𝐽=𝜅_s∇𝑇⋅𝐧_g, (2.5)where 𝐧_g 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 have𝑅=𝐸+𝐽_ref,𝐸=𝜎𝜖𝑇⁴𝐽(Stefan-Boltzmannequation),_ref=(1−𝜖)𝐽, (2.6)where 𝐸 is the radiation, 𝐽_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.

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 𝑚∈{s,g} (with s solid, g 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,𝜌_{𝑚𝜔𝑚,𝑖}(𝑈_𝑚(𝑇^𝑛+1)−𝑈_𝑚(𝑇^𝑛))𝑟𝑑𝑥−Δ𝑡^𝑛+1_{𝜕𝜔𝑚,𝑖}𝜅_𝑚(𝑇^𝑛+1)∇𝑇^𝑛+1⋅𝐧_{𝜔𝑚,𝑖}𝑟𝑑𝑠=Δ𝑡^𝑛+1_{𝜔𝑚,𝑖}𝑓_𝑚𝑟𝑑𝑥, (2.7)where the time interval is Δ𝑡^𝑛+1=𝑡^𝑛+1−𝑡^𝑛. The temperature is given as 𝑇^𝑛+1=𝑇(𝑡^𝑛+1,𝑥), where 𝑥 represents cylindrical coordinates. For the right-hand sides, we demand 𝑓_s∶=𝑓≥0 and 𝑓_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.

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

Figure 1 

The growth apparatus’ dimensions: 𝑟min=0, 𝑟max=8.4cm, 𝑧min=0, 𝑧max=25.0cm, the coil rings’ dimensions: 𝑟min=4.2cm, 𝑟max=5.2cm, 𝑧min=0,𝑧max=14.0cm.

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³, the specific heat 𝑐_sp is given in J/(K kg), the temperature 𝑇 is given in K, and the relative gas constant 𝑅_Ar is given in J/(K kg). Further, the emissivity 𝜖 and relative magnetic permeability 𝜇 are given dimensionless.

(i) For the gas phase (argon), we have𝜅_Ar(𝑇)=1.8391410⁻⁴𝑇^0.800404𝑇≤500,−7.128738+6.61028810⁻²𝑇−2.44083910⁻⁴𝑇²+4.49763310⁻⁷𝑇³−4.13251710⁻¹⁰𝑇⁴+1.51446310⁻¹³𝑇⁵500≤𝑇≤600,4.194410⁻⁴𝑇^0.671118600≥𝑇, (2.8) where 𝜎_𝑐,Ar=0.0,𝜌_Ar=3.7310⁻³,𝜇_Ar=1.0,𝑧_Ar=3/2,𝑅_Ar=2.08130810⁻².

(ii) For graphite felt insulation, we have𝜅_Ins(𝑇)=8.17510⁻²+2.48510⁻⁴𝑇𝑇≤1473,−1.190210²+0.346838𝑇−3.997110⁻⁴𝑇²+2.283010⁻⁷𝑇³−6.4604710⁻¹¹𝑇⁴+7.254910⁻¹⁵𝑇⁵1473≤𝑇≤1873,−0.7447+7.510⁻⁴𝑇1873≥𝑇, (2.9) where 𝜖_Ins=0.2,𝜎_𝑐,Ins(𝑇)=2.4510²+9.8210⁻²𝑇,𝜌_Ins=170.00,𝜇_Ins=1.00,𝑐_sp,Ins=2100.00.

(iii) For the graphite, we have𝜅_Graphite(𝑇)=37.715exp(−1.9610⁻⁴𝜖𝑇),_Graphite(𝑇)=0.67𝑇≤1200,3.752−7.43610⁻³𝑇+6.416310⁻⁶𝑇²−2.336610⁻⁹𝑇³−3.083310⁻¹³𝑇⁴1200≤𝑇≤2200,0.792200≥𝑇, (2.10) where 𝜎_c,Graphite=10⁴,𝜌_Graphite=1750.0,𝜇_Graphite=1.0,𝑐_sp,Graphite(𝑇)=1/(4.41110²𝑇^−2.306+7.9710⁻⁴𝑇^−0.0665).

(iv) For the SiC crystal, we have 𝜅_SiC-C(𝑇)=exp(9.892+(2.49810²)/𝑇−0.844ln(𝑇)),𝜖_SiC-C=0.85,𝜎_c,SiC-C=10⁵,𝜌_SiC-C=3140.0,𝜇_SiC-C=1.0,𝑐_sp,SiC-C(𝑇)=1/(3.9110⁴𝑇^−3.173+1.83510⁻³𝑇^−0.117).

(v) For the SiC powder, we have 𝜅_SiC-P(𝑇)=1.45210⁻²+5.4710⁻¹²𝑇³,𝜖_SiC-P=0.85,𝜎_c,SiC-P=100.0,𝜌_SiC-P=1700.0,𝜇_SiC-P=1.0,𝑐_sp,SiC-P=1000.0.

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 [𝑡^𝑛,𝑡^𝑛+1] that is sufficiently small to yield accurate kinetics, we solve the following subproblems consecutively for 𝑖=0,2,…,2𝑚 (cf. [14, 15]):𝜕𝑐_𝑖(𝑡)𝜕𝑡=𝐴𝑐_𝑖(𝑡)+𝐵𝑐_𝑖−1(𝑡),with𝑐_𝑖(𝑡^𝑛)=𝑐^𝑛,𝑖=1,2,…,𝑗,𝜕𝑐_𝑖(𝑡)𝜕𝑡=𝐴𝑐_𝑖−1(𝑡)+𝐵𝑐_𝑖(𝑡),with𝑐_𝑖+1(𝑡^𝑛)=𝑐^𝑛,𝑖=𝑗+1,𝑗+2,…,𝑚, (2.11)(2.12) where we assume that the operator 𝐴 has a large time scale (macroscopic model) and 𝐵 has a small time scale (microscopic model). Further 𝑐₀(𝑡^𝑛)=𝑐^𝑛,𝑐₋₁=0, 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 𝐗∶𝜕_𝑡𝑐(𝑡)=𝐴𝑐(𝑡)+𝐵𝑐(𝑡),0<𝑡≤𝑇,𝑐(0)=𝑐₀, (2.13) 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 𝒪(𝜏_𝑛^2𝑚), where the iterations are 𝑖=1,3,…,2𝑚+1.

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 [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 𝑇_𝑡; 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.

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.

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 𝑇_𝑡=2000K for 1.20 picoseconds, and 3 structures were selected as initial starting geometries at 𝑡=0.60 picosecond (trajectory A50), 𝑡=0.72 picosecond (trajectory A60), and 𝑡=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 𝑇_𝑡. 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 10kW (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 𝑇_ss=𝑇(𝑟_source,𝑧_source)−𝑇(𝑟_seed,𝑧_seed) (with the coordinates (𝑟_source,𝑧_source)=(0,0.143) and (𝑟_seed,𝑧_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 𝑇(𝑟_bottom,𝑧_bottom) and 𝑇(𝑟_top,𝑧_top) (with the coordinates (𝑟_bottom,𝑧_bottom)=(0,0.028) and (𝑟_top,𝑧_top)=(0,0.173), corresponding to the points T_bottom and T_top in Figure 1) are measurable and their difference 𝑇_bt=𝑇(𝑟_bottom,𝑧_bottom)−𝑇(𝑟_top,𝑧_top) is often used as an indicator for 𝑇_ss. In Figure 3, we present the temperature differences 𝑇_ss and 𝑇_bt. As a result of our computations, the temperature difference 𝑇_bt can only restrictively be used as an indicator for the temperature difference 𝑇_ss (cf. the discussions in [5, 9]).

Figure 3 

Transient results for the temperature differences 𝑇bt and 𝑇ss.

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 𝑡>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 𝑇(𝑟_bottom,𝑧_bottom) and 𝑇(𝑟_top,𝑧_top) for successive grids. For the error analysis, we apply the following error differences:𝜖_abs=|𝑇^{0𝑥002𝑑𝑐}_𝑗+1(𝑟,𝑧)−𝑇^{0𝑥002𝑑𝑐}_𝑗(𝑟,𝑧)|, (4.1)where ^0𝑥0𝑇^02𝑑𝑐_𝑗(𝑟,𝑧) and ^0𝑥0𝑇^02𝑑𝑐_𝑗+1(𝑟,𝑧) are solutions evaluated at the point (𝑟,𝑧) which has been computed using the grids 𝑗 and 𝑗+1, respectively. The elements of the grid 𝑗+1 are approximately 1/4 of the elements of the grid 𝑗. The results are presented in Table 1.

Table 1 

Computations on different grids for the errors analysis with absolute differences (cf. (4.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.

Figure 4 

Temperature field for the apparatus simulated for the stationary case with 23017 nodes.

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.

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.

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₂ units, and behave similar to “wobbling C₂” 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₂ units are oriented mainly perpendicular to the surface, while the more visible Si₂ 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₂ and Si₂ 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.