Research Article  Open Access
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 highquality 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 highvalued 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: HeatFlux
In the following, we discuss the macroscopic model, which is based on continuum equations for the heatflux.
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 WIASHiTNIHS (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 righthand 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 WIASHiTNIHS.
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., targettemperature 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 socalled iterative operatorsplitting 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 stepsize
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 densityfunctional tightbinding (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 Slatertype allvalence basis set in combination with a twocenter approximation for Hamiltonian matrix elements. Parameter sets for SiSi and SiC 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 Sifaces of SiC(0001) surfaces during sublimation evaporation, we have not included charge or spinpolarization in the present work. Further, we will consider in a next model electrokinetic effect on heat transfer in parallelplate 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 NoseHoover 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 Cface of the same square SiC(0001) 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 surfacedirections to achieve twodimensional 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 bulkside 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 WIASHiTNIHS (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 CFace of Si(0001)
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 Csurface of SiC (0001) since it has a maximum of dangling bonds with highest reactivity. The Siface 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 CC 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 hightemperature 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 CC 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 densityfunctional tightbinding (DFTB) method for initial reactions of gaseous SiC on the polar Cface of SiC(0001). 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/MDbased 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.
References
 St. Müller, R. C. Glass, H. M. Hobgood et al., “Progress in the industrial production of SiC substrates for semiconductor devices,” Materials Science and Engineering B, vol. 80, no. 1–3, pp. 327–331, 2001. View at: Publisher Site  Google Scholar
 X. Yao, B. He, H. Wang, and X. Zhou, “Numerical simulation of dendrite growth during solidification,” International Journal of Nonlinear Sciences and Numerical Simulation, vol. 7, no. 2, pp. 171–176, 2006. View at: Google Scholar
 O. Klein and P. Philip, “Transient numerical investigation of induction heating during sublimation growth of silicon carbide single crystals,” Journal of Crystal Growth, vol. 247, no. 12, pp. 219–235, 2003. View at: Publisher Site  Google Scholar
 P. Philip, Transient numerical simulation of sublimation growth of SiC bulk single crystals. Modeling, finite volume method, results, Ph.D. thesis, WeierstrassInstitute for Applied Analysis and Stochastics, Berlin, Germany, 2003, Report No. 22.
 O. Klein, P. Philip, and J. Sprekels, “Modeling and simulation of sublimation growth of SiC bulk single crystals,” Interfaces and Free Boundaries, vol. 6, no. 3, pp. 295–314, 2004. View at: Google Scholar  Zentralblatt MATH  MathSciNet
 J. Rappaz and M. Swierkosz, “Modelling in numerical simulation of electromagnetic heating,” in Modelling and Optimization of Distributed Parameter Systems, pp. 313–320, Chapman & Hall, New York, NY, USA, 1996. View at: Google Scholar  Zentralblatt MATH  MathSciNet
 O. Klein and P. Philip, “Correct voltage distribution for axisymmetric sinusoidal modeling of induction heating with prescribed current, voltage, or power,” IEEE Transactions on Magnetics, vol. 38, no. 3, pp. 1519–1523, 2002. View at: Publisher Site  Google Scholar
 J. Geiser, “${R}^{3}$$T$: radioactiveretardationreactiontransportprogram for the simulation of radioactive waste disposals,” Institute for Scientific Computation, Texas A&M University, College Station, Tex, USA, April 2004. View at: Google Scholar
 O. Klein and P. Philip, “Transient temperature phenomena during sublimation growth of silicon carbide single crystals,” Journal of Crystal Growth, vol. 249, no. 34, pp. 514–522, 2003. View at: Publisher Site  Google Scholar
 O. Klein and P. Philip, “Transient conductiveradiative heat transfer: discrete existence and uniqueness for a finite volume scheme,” Mathematical Models & Methods in Applied Sciences, vol. 15, no. 2, pp. 227–258, 2005. View at: Google Scholar  Zentralblatt MATH  MathSciNet
 M. Pons, M. Anikin, K. Chourou et al., “State of the art in the modelling of SiC sublimation growth,” Materials Science and Engineering B, vol. 6162, pp. 18–28, 1999. View at: Publisher Site  Google Scholar
 O. Nilsson, H. Mehling, R. Horn et al., “Determination of the thermal diffusivity and conductivity of monocrystalline silicon carbide (300–2300 K),” High TemperaturesHigh Pressures, vol. 29, no. 1, pp. 73–79, 1997. View at: Publisher Site  Google Scholar
 J. Geiser, “Iterative operatorsplitting methods with higherorder time integration methods and applications for parabolic partial differential equations,” Journal of Computational and Applied Mathematics, vol. 217, no. 1, pp. 227–242, 2008. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 J. F. Kanney, C. T. Miller, and C. T. Kelley, “Convergence of iterative splitoperator approaches for approximating nonlinear reactive problems,” Advances in Water Resources, vol. 26, no. 3, pp. 247–261, 2003. View at: Publisher Site  Google Scholar
 I. Faragó and J. Geiser, “Iterative operatorsplitting methods for linear problems,” International Journal of Computational Science and Engineering, vol. 3, no. 4, pp. 255–263, 2007. View at: Publisher Site  Google Scholar
 D. Porezag, Th. Frauenheim, Th. Köhler, G. Seifert, and R. Kaschner, “Construction of tightbindinglike potentials on the basis of densityfunctional theory: application to carbon,” Physical Review B, vol. 51, no. 19, pp. 12947–12957, 1995. View at: Publisher Site  Google Scholar
 E. Rauls, Z. Hajnal, P. Deák, and Th. Frauenheim, “Theoretical study of the nonpolar surfaces and their oxygen passivation in 4H and 6HSiC,” Physical Review B, vol. 64, no. 24, Article ID 245323, 7 pages, 2001. View at: Publisher Site  Google Scholar
 S. Irle, Z. Wang, G. Zheng, K. Morokuma, and M. Kusunoki, “Theory and experiment agree: singlewalled carbon nanotube caps grow catalystfree with chirality preference on a SiC surface,” The Journal of Chemical Physics, vol. 125, no. 4, Article ID 044702, 5 pages, 2006. View at: Publisher Site  Google Scholar
 Z. Wang, S. Irle, G. Zheng, M. Kusunoki, and K. Morokuma, “Carbon nanotubes grow on the C face of SiC (0001) during sublimation decomposition: quantum chemical molecular dynamics simulations,” The Journal of Physical Chemistry C, vol. 111, no. 35, pp. 12960–12972, 2007. View at: Publisher Site  Google Scholar
 C.C. Chang and R.J. Yang, “Hydrodynamic focusing effect on twounmixedfluid in microchannels,” International Journal of Nonlinear Sciences and Numerical Simulation, vol. 9, no. 3, pp. 213–220, 2008. View at: Google Scholar
 J. Fan, L. Wang, and L. Cheng, “Electrokinetic effects on flow and heat transfer in parallelplate microchannels,” International Journal of Nonlinear Sciences and Numerical Simulation, vol. 8, no. 3, pp. 335–345, 2007. View at: Google Scholar
 J. Fan, L. Wang, and L. Cheng, “Forced convection in rectangular microchannels: electrokinetic effects,” International Journal of Nonlinear Sciences and Numerical Simulation, vol. 8, no. 3, pp. 359–374, 2007. View at: Google Scholar
 J.H. He, Y.Q. Wan, and L. Xu, “Nanoeffects, quantumlike properties in electrospun nanofibers,” Chaos, Solitons & Fractals, vol. 33, no. 1, pp. 26–37, 2007. View at: Publisher Site  Google Scholar
 M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford University Press, New York, NY, USA, 1989.
 J. Fuhrmann, Th. Koprucki, and H. Langmach, “pdelib: an open modular tool box for the numerical solution of partial differential equations. Design patterns,” in Proceeding of the 14th GAMM Seminar on Concepts of Numerical Software, W. Hackbusch and G. Wittum, Eds., Kiel, Germany, January 2001. View at: Google Scholar
 O. Schenk and K. Gärtner, “Solving unsymmetric sparse systems of linear equations with PARDISO,” Future Generation Computer Systems, vol. 20, no. 3, pp. 475–487, 2004. View at: Publisher Site  Google Scholar
 J. Pollmann, P. Krüger, and M. Sabisch, “Atomic and electronic structure of SiC surfaces from abinitio calculations,” Physica Status Solidi (B), vol. 202, no. 1, pp. 421–445, 1997. View at: Publisher Site  Google Scholar
 S. Irle, G. Zheng, M. Elstner, and K. Morokuma, “Formation of fullerene molecules from carbon nanotubes: a quantum chemical molecular dynamics study,” Nano Letters, vol. 3, no. 4, pp. 465–470, 2003. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2009 Jürgen Geiser and Stephan Irle. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.