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

Authors: Jürgen Geiser; Stephan Irle
Journal: Mathematical Problems in Engineering
 (2009)

---

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


---

## Body

## 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(2.1)ρg∂tUg−∇⋅(κg∇T)=0,(2.2)Ug=zArRArT,where T is the temperature, t is the time, and Ug 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, zAr being the configuration number, and RAr 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(2.3)ρs∂tUs−∇⋅(κs∇T)=f,(2.4)Us=∫0Tcs(S)dS,where ρs is the density of the solid material, Us is the internal energy, κs is the thermal conductivity, and cs 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(2.5)κg∇T⋅ng+R−J=κs∇T⋅ng,where ng 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(2.6)R=E+Jref,E=σϵT4(Stefan-Boltzmannequation),Jref=(1−ϵ)J,where E is the radiation, Jref 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 𝒯=(ω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,(2.7)ρm∫ωm,i(Um(Tn+1)−Um(Tn))rdx−Δtn+1∫∂ωm,iκm(Tn+1)∇Tn+1⋅nωm,irds=Δtn+1∫ωm,ifmrdx,where the time interval is Δtn+1=tn+1−tn.
The temperature is given as Tn+1=T(tn+1,x),
where x represents cylindrical coordinates. For the
right-hand sides, we demand fs:=f≥0 and fg=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:rmin=0,rmax=8.4cm,zmin=0,zmax=25.0cm,
the coil rings’ dimensions: rmin=4.2cm,rmax=5.2cm,zmin=0,zmax=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/m3,
the specific heat csp is given in J/(K kg),
the temperature T is given in K, and the relative gas constant RAr is given in J/(K kg).
Further, the emissivity ϵ and relative magnetic permeability μ are given dimensionless.(i) For the gas  phase (argon), we have(2.8)κAr(T)={1.8391410−4T0.800404T≤500,−7.128738+6.61028810−2T−2.44083910−4T2+4.49763310−7T3−4.13251710−10T4+1.51446310−13T5500≤T≤600,4.194410−4T0.671118600≥T,  where σc,Ar=0.0,ρAr=3.7310−3,μAr=1.0,zAr=3/2,RAr=2.08130810−2.(ii) For graphite felt  insulation, we
have(2.9)κIns(T)={8.17510−2+2.48510−4TT≤1473,−1.1902102+0.346838T−3.997110−4T2+2.283010−7T3−6.4604710−11T4+7.254910−15T51473≤T≤1873,−0.7447+7.510−4T1873≥T, where ϵIns=0.2,σc,Ins(T)=2.45102+9.8210−2T,ρIns=170.00,μIns=1.00,csp,Ins=2100.00.(iii) For the  graphite, we have(2.10)κGraphite(T)=37.715exp(−1.9610−4T),ϵGraphite(T)={0.67T≤1200,3.752−7.43610−3T+6.416310−6T2−2.336610−9T3−3.083310−13T41200≤T≤2200,0.792200≥T, where σc,Graphite=104,ρGraphite=1750.0,μGraphite=1.0,csp,Graphite(T)=1/(4.411102T−2.306+7.9710−4T−0.0665).(iv) For the SiC crystal, we haveκSiC-C(T)=exp(9.892+(2.498102)/T−0.844ln(T)),ϵSiC-C=0.85,σc,SiC-C=105,ρSiC-C=3140.0,μSiC-C=1.0,csp,SiC-C(T)=1/(3.91104T−3.173+1.83510−3T−0.117).(v) For the SiC powder, we haveκSiC-P(T)=1.45210−2+5.4710−12T3,ϵSiC-P=0.85,σc,SiC-P=100.0,ρSiC-P=1700.0,μSiC-P=1.0,csp,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[tn,tn+1] that is sufficiently small to yield accurate
kinetics, we solve the following subproblems consecutively for i=0,2,…,2m (cf. [14, 15]):(2.11)∂ci(t)∂t=Aci(t)+Bci−1(t),withci(tn)=cn,i=1,2,…,j,(2.12)∂ci(t)∂t=Aci−1(t)+Bci(t),withci+1(tn)=cn,i=j+1,j+2,…,m, where we assume
that the operator A has a large time scale (macroscopic model) and B has a small time scale (microscopic model).
Further c0(tn)=cn,c−1=0, and cn 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 spaceX:(2.13)∂tc(t)=Ac(t)+Bc(t),0<t≤T,c(0)=c0, where A,B,A+B:X→X are given linear operators being generators of
the C0-semigroup and c0∈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 𝒪(τn2m),
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.


## 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(2.1)ρg∂tUg−∇⋅(κg∇T)=0,(2.2)Ug=zArRArT,where T is the temperature, t is the time, and Ug 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, zAr being the configuration number, and RAr 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(2.3)ρs∂tUs−∇⋅(κs∇T)=f,(2.4)Us=∫0Tcs(S)dS,where ρs is the density of the solid material, Us is the internal energy, κs is the thermal conductivity, and cs 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(2.5)κg∇T⋅ng+R−J=κs∇T⋅ng,where ng 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(2.6)R=E+Jref,E=σϵT4(Stefan-Boltzmannequation),Jref=(1−ϵ)J,where E is the radiation, Jref 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 𝒯=(ω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,(2.7)ρm∫ωm,i(Um(Tn+1)−Um(Tn))rdx−Δtn+1∫∂ωm,iκm(Tn+1)∇Tn+1⋅nωm,irds=Δtn+1∫ωm,ifmrdx,where the time interval is Δtn+1=tn+1−tn.
The temperature is given as Tn+1=T(tn+1,x),
where x represents cylindrical coordinates. For the
right-hand sides, we demand fs:=f≥0 and fg=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:rmin=0,rmax=8.4cm,zmin=0,zmax=25.0cm,
the coil rings’ dimensions: rmin=4.2cm,rmax=5.2cm,zmin=0,zmax=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/m3,
the specific heat csp is given in J/(K kg),
the temperature T is given in K, and the relative gas constant RAr is given in J/(K kg).
Further, the emissivity ϵ and relative magnetic permeability μ are given dimensionless.(i) For the gas  phase (argon), we have(2.8)κAr(T)={1.8391410−4T0.800404T≤500,−7.128738+6.61028810−2T−2.44083910−4T2+4.49763310−7T3−4.13251710−10T4+1.51446310−13T5500≤T≤600,4.194410−4T0.671118600≥T,  where σc,Ar=0.0,ρAr=3.7310−3,μAr=1.0,zAr=3/2,RAr=2.08130810−2.(ii) For graphite felt  insulation, we
have(2.9)κIns(T)={8.17510−2+2.48510−4TT≤1473,−1.1902102+0.346838T−3.997110−4T2+2.283010−7T3−6.4604710−11T4+7.254910−15T51473≤T≤1873,−0.7447+7.510−4T1873≥T, where ϵIns=0.2,σc,Ins(T)=2.45102+9.8210−2T,ρIns=170.00,μIns=1.00,csp,Ins=2100.00.(iii) For the  graphite, we have(2.10)κGraphite(T)=37.715exp(−1.9610−4T),ϵGraphite(T)={0.67T≤1200,3.752−7.43610−3T+6.416310−6T2−2.336610−9T3−3.083310−13T41200≤T≤2200,0.792200≥T, where σc,Graphite=104,ρGraphite=1750.0,μGraphite=1.0,csp,Graphite(T)=1/(4.411102T−2.306+7.9710−4T−0.0665).(iv) For the SiC crystal, we haveκSiC-C(T)=exp(9.892+(2.498102)/T−0.844ln(T)),ϵSiC-C=0.85,σc,SiC-C=105,ρSiC-C=3140.0,μSiC-C=1.0,csp,SiC-C(T)=1/(3.91104T−3.173+1.83510−3T−0.117).(v) For the SiC powder, we haveκSiC-P(T)=1.45210−2+5.4710−12T3,ϵSiC-P=0.85,σc,SiC-P=100.0,ρSiC-P=1700.0,μSiC-P=1.0,csp,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[tn,tn+1] that is sufficiently small to yield accurate
kinetics, we solve the following subproblems consecutively for i=0,2,…,2m (cf. [14, 15]):(2.11)∂ci(t)∂t=Aci(t)+Bci−1(t),withci(tn)=cn,i=1,2,…,j,(2.12)∂ci(t)∂t=Aci−1(t)+Bci(t),withci+1(tn)=cn,i=j+1,j+2,…,m, where we assume
that the operator A has a large time scale (macroscopic model) and B has a small time scale (microscopic model).
Further c0(tn)=cn,c−1=0, and cn 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 spaceX:(2.13)∂tc(t)=Ac(t)+Bc(t),0<t≤T,c(0)=c0, where A,B,A+B:X→X are given linear operators being generators of
the C0-semigroup and c0∈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 𝒪(τn2m),
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.


## 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
temperatureTt;
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 atTt=2000K 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 Tt.
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
ratev 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 Section2.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 Tss=T(rsource,zsource)−T(rseed,zseed) (with the coordinates (rsource,zsource)=(0,0.143) and (rseed,zseed)=(0,0.158),
corresponding to the points Tsource and Tseed in Figure 1) is crucial. On the other hand,
in the physical growth experiments, usually only the temperatures T(rbottom,zbottom) and T(rtop,ztop) (with the coordinates (rbottom,zbottom)=(0,0.028) and (rtop,ztop)=(0,0.173),
corresponding to the points Tbottom and Ttop in Figure 1) are measurable and their
difference Tbt=T(rbottom,zbottom)−T(rtop,ztop) is often used as an indicator for Tss.
In Figure 3, we present the temperature differences Tss and Tbt.
As a result of our computations, the temperature difference  Tbt can only restrictively be used as an indicator
for the temperature difference  Tss (cf. the discussions in [5, 9]).Figure 3
Transient
results for the temperature differencesTbt and Tss.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(rbottom,zbottom) and T(rtop,ztop) for successive grids. For the error analysis,
we apply the following error differences:(4.1)ϵabs=|T˜j+1(r,z)−T˜j(r,z)|,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.Table 1
Computations on
different grids for the errors analysis with absolute differences (cf. (4.1)).
GridGrid point(0,0.028) (Tbottom)Grid point(0,0.173) (Ttop)LevelNumber of nodesSolutionT[K]Absolute differenceT[K]SolutionT[K]Absolute differenceT[K]015322408.112813.291230172409.781.672812.781.012912902410.350.572811.790.4933642252410.460.112811.600.19The 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 Figure4, 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 with23017 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.Tt 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
Figure2. 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 C2 units, and behave similar to “wobbling C2” 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, C2 units are oriented mainly perpendicular to the
surface, while the more visible Si2 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 C2 and Si2 units at first, and gradually becoming more
homogeneous due to annealing.


## 4.1. Macroscopic Model: Simulation of the Temperature Field in the Apparatus

For the numerical results, we apply the parameter
functions in Section2.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 Tss=T(rsource,zsource)−T(rseed,zseed) (with the coordinates (rsource,zsource)=(0,0.143) and (rseed,zseed)=(0,0.158),
corresponding to the points Tsource and Tseed in Figure 1) is crucial. On the other hand,
in the physical growth experiments, usually only the temperatures T(rbottom,zbottom) and T(rtop,ztop) (with the coordinates (rbottom,zbottom)=(0,0.028) and (rtop,ztop)=(0,0.173),
corresponding to the points Tbottom and Ttop in Figure 1) are measurable and their
difference Tbt=T(rbottom,zbottom)−T(rtop,ztop) is often used as an indicator for Tss.
In Figure 3, we present the temperature differences Tss and Tbt.
As a result of our computations, the temperature difference  Tbt can only restrictively be used as an indicator
for the temperature difference  Tss (cf. the discussions in [5, 9]).Figure 3
Transient
results for the temperature differencesTbt and Tss.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(rbottom,zbottom) and T(rtop,ztop) for successive grids. For the error analysis,
we apply the following error differences:(4.1)ϵabs=|T˜j+1(r,z)−T˜j(r,z)|,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.Table 1
Computations on
different grids for the errors analysis with absolute differences (cf. (4.1)).
GridGrid point(0,0.028) (Tbottom)Grid point(0,0.173) (Ttop)LevelNumber of nodesSolutionT[K]Absolute differenceT[K]SolutionT[K]Absolute differenceT[K]015322408.112813.291230172409.781.672812.781.012912902410.350.572811.790.4933642252410.460.112811.600.19The 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 Figure4, 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 with23017 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.Tt 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
Figure2. 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 C2 units, and behave similar to “wobbling C2” 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, C2 units are oriented mainly perpendicular to the
surface, while the more visible Si2 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 C2 and Si2 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.



---

*Source: 716104-2009-02-03.xml*