The aim of this work is to better understand the basic physical phenomena occurring in the growth reactors by making use of the numerical modeling. This present thesis will be structured as follow.

Chapter 2 is devoted to provide an overview of SiC material, its growth processes and modeling. In the first part, its crystal structure, polytypes, properties, and its applications in power electronic devices will be given. Secondly, the bulk growth processes will be described and the significant contribution of the process modeling will be discussed regarding the process development. The purposes of this thesis will then be detailed.

In Chapter 3, an introduction to the Finite Element Method (FEM) which is the basis of the numerical modeling used throughout this thesis will be given. This method is capable of solving partial differential equations (PDEs). It has been proved to be a powerful tool in the modeling of electromagnetic problems of irregular geometry. This is perfectly adapted to our case in which the induction heating of complex reactor geometry has been used. An overview of the method will be first given. Then we will provide an example of using FEM for modeling the heat transfer in steady state. Since the PDEs associated with the growth processes are generally nonlinear, an approach for solving nonlinear system will be discussed at the end of the chapter.

Chapter 4 is dedicated to the modeling of induction heating and heat transfer. This chapter is divided into four main sections. First, the physical model including induction heating and heat transfer will be given. The boundary conditions for electromagnetic and heat transfer will also be given in the second section. In the third section, we will discuss the model treating SiC as a semi-transparent media. Its spatial and temporal effect on the growth process and the limitations of such model will be also discussed. The last section focuses on the modeling results and validations. This includes the computed I-V characteristic and measurement, temperature distribution, and the effect of coil position.

In Chapter 5, the modeling of mass transport that couples the vapor species transport and the fluid flow in the PVT system will be first given. This includes again the physical model and boundary conditions. Such coupled model make used of the temperature profile computed from the induction heating and heat transfer model. Then special attention will be paid to the chemical reactions at the solid-vapor interfaces. Both 'comparative' and 'coupled' approaches will be

i

General information

This PhD thesis was conducted at the Laboratoire des Matériaux et du Génie Physique (LMGP) which is a joint research unit of CNRS and Grenoble INP. The laboratory is a research laboratory in materials science and materials for biomedical engineering and it is organized into 3 research teams: Crystal Growth (XTO), Interface between Materials and Biological Matter (IMBM), and Thin films, Nanomaterials and Nanostructures (FM2N).

The current thesis was financially supported by the European Union Seventh Framework

Programme under an acronym NetFISiC (Training Network on Functional Interfaces for SiC).

The main scientific of the project objective is to provide Silicon carbide material (of various polytypes) with improved and adequate functional interfaces for getting a step forward in electronic devices performance. Research efforts were dedicated to solve the problems faced by important devices like MOSFET and Schottky diodes. Besides, some fundamental research was performed both on the growth aspect and on new and innovating devices. Applications in high temperature, high power and hash environment were targeted.

The NetFISiC consortium is made of 12 European partners, including 3 companies, coming This thesis was elaborated within the WP1. A close collaboration with another ESR recruited at the same laboratory and within the same work-package was made. The collaboration was focused on developing a better comprehension of the growth of SiC by PVT process using the complete experiments, characterization, and modeling approach. The topics that were tackled ii are development of the process, interaction between vapor species with the crucible, foreign polytype nucleation, and nitrogen doping. The experimental results and characterization are given in Chapter 4 and 5.

Apart from the NetFISiC network, collaborations were made with Laboratoire de Science et Ingénierie des Matériaux et Procédés (SIMaP) for the thermodynamic calculations (results are given in Chapter 5) and within LMGP for the combined experiments, characterization, and numerical modeling approach for better visualization and understanding the growth of SiC by Top Seeded Solution Growth (TSSG) process (results are given in Chapter 6).

In total, eight conference proceedings including one oral presentation, and one journal paper were obtained. In addition, four journal papers are in preparations.

Acknowledgement

It is a great pleasure to acknowledge the many people for their help, support, and all kind of effort throughout the period of my PhD thesis.

First of all, I wish to express my deepest gratitude to my supervisor Prof. Jean-Marc Dedulle and my co-supervisor Dr. Didier Chaussende for giving me the opportunity to work with an interesting topic and the opportunity to collaborate with them. From the very first day until the last moment of my thesis, they gave me not only many helpful suggestions for improving my works, but they also shown me the way of becoming a real scientist. With his kindness and clear explanations, Jean-Marc really paved me the way for understanding and advancing the field of process modeling. From the many discussions we had, Didier undoubtedly helped and impressed me by his real experience and perspectives as a crystal grower. Special thanks would go to Dr. Gabriel Ferro. Through his sympathetic questions and advices, Gabriel certainly comforted me and assisted me to improve my presentation skill. I will keep all their guidance and effort as an invaluable gift they gave to me.

Second, I would like to thank Dr. Elisabeth Blanquet and Prof. Christian Chatillon for their very kind help, collaboration, and discussion throughout the work related to the thermodynamic calculations. Special thanks go to Prof. Michel Pons for his kindness and many fruitful discussions as one of the experts in the process modeling of SiC. I would really appreciate to thank to Dr. Frederic Mercier for his gentle explanations and discussions concerning the modeling of the TSSG process. I would like also to thank to the current and former members of the Cristallogenese group, Roland Madar, Thierry Ouisse, Odette Chaix, Eirini Sarigianaidou, Pierre Bouvier, Joseph La Manna, Lucile Parent-Bert, Julien Lefebure, Sven Tengeler, Lu Shi, Martin Seisse, Le Thi Mai Hoa, Hoang-Long LE-TRAN, and Demitrios Zevgitis for their help, support, and friendships. We shared many great times both inside the laboratory and many events arranged outside. For me, they have been all like my family. Special thanks would go to my colleague Dr. Nikolaos Tsavdaris and PhD student Yun Ji Shin. Apart from the daily scientific discussions and their contributions to the many results in my thesis, they have been my good friends. We have been through both happiness and anxiousness moments. I am honestly pleased to collaborate and be friend with them. I would also appreciate to thank to all the members of Laboratoire des Materiaux et du Genie Physique for their kindness and hospitality. Their warm greeting created a generous working environment. Exclusively thanks would go to the former director Dr. Bernard Chenevier and the current director Dr. Franz Bruckert of the lab. The technical assistance for computing resource from Arnaud Robin is also greatly acknowledged.

I am sincerely grateful to the consortium of a European project NetFISiC (young researcher and their suporvisors) for sharing their views and experience in the field of SiC material and related techniques during the common shared training events and conferences which built me a strong scientific background in the research field. Besides, the friendships created within the network are incalculable.

I should thank to the many people writing papers as well as textbooks which I used as the references. Without them, I would not have enough information and idea to complete this thesis.

Finally, it is difficult to exaggerate my appreciation to my parents and my brother.

Throughout my life, they always support me, teach me, encourage me, as well as love me.

Chapter 1 Introduction

Silicon carbide (SiC) is a wide bandgap semiconductor material. Its properties make it adapted for electronic devices requiring high temperature, high frequency, high power, and/or hash environment operation. SiC exists in different kind of polytypes which are coming from a one-dimensional polymorphism. The difference in stacking sequences leads to different properties of each polytype, such as, bandgap, electrical breakdown voltage. This broadens the area of SiC application. The most technologically important polytypes are 3C, 4H, 6H, and 15R-SiC.

For SiC single crystal growth, the Physical Vapor Transport (PVT) process has shown its ability for the growth of high quality and large size crystals. The demonstration of 6 inch diameter SiC wafers has been achieved. "Defect-free" wafer was also claimed to be realized. However, only the 4H and 6H-SiC polytypes are commercially available. There is no 3C or 15R-SiC wafers with adequate quality for electronic devices owing to the difficulties in getting the stable growth conditions for these two polytypes. Besides, the liquid phase route, especially the Top Seeded Solution Growth (TSSG) process has brought back a huge attention concerning research and development in the past ten years. The "lower" growth temperatures used in such process and the "closer" to thermodynamic equilibrium conditions are expected to provide better control of defects generation. Though lower than PVT, the growth temperatures are still high, which makes the process very difficult to control. Due to the high temperatures (typically about 2000 °C), there is no way to measure the process parameters. The best approach to access such information is to use the modeling which serves as a tool for "process characterization". The proper model is necessary to pave the way first for a better understanding and second for a better monitoring of the growth process.

discussed. Those approaches are characterized by the strength of the coupling between mass transport modeling and thermodynamics. In addition to the effect of the thermodynamic databases in the full process modeling using the former approach, the capability and the extension possibility of the latter approach will be discussed. Finally, general modeling results and validations will be presented.

Chapter 6 contributes to the modeling of the TSSG process. Similar to the two former chapters, physical model and boundary conditions will be first discussed. This is followed by the study of the effects of the process parameters, namely, crystal rotation speed and operating frequency, on the growth of SiC crystals. In the last part of this chapter, a combined numerical and analytical model will be presented. This model has been developed for getting a comprehensive visualization of the interaction between the fluid flow and the step flow that will affect the surface morphology of the crystal. The validation of this model is included.

The last chapter will give the general conclusions of this thesis. Open issues and the perspective for the future work will also be discussed.

5

Chapter 2

Silicon Carbide

Growth processes and modeling Silicon carbide (SiC) is a wide bandgap semiconductor material. Its potential is far beyond silicon and very highly desirable for electronic devices, especially for high power, high temperature, and high frequency applications [1]. In this chapter, SiC crystal structure and polytypes will be described followed by its basic properties and electronic applications. Then, an overview of the growth processes for producing SiC single crystals, focusing on the bulk growth technique, will be given. Finally, the importance of numerical modeling for SiC growth processes and the aims of this thesis will be discussed.

SiC crystal structure and polytypes

Silicon carbide is a solid compound consisting of 50% of carbon and 50% of silicon.

Thus the stoichiometric coefficient for both silicon and carbon are unity. Concerning the structure of carbon and silicon atoms in SiC, all silicon (carbon) atoms are surrounded by four carbon (silicon) atoms forming a tetrahedron (Fig. 2.1). The nearest neighbors to any Si or C atom in the crystal are similar. The distance between same type of atom is approximately 3.08 Å [2] with small variation depending on the polytype. The bonding between Si and C atoms of SiC composes of primary covalent and partly ionic (88% covalent and 12% ionic) where the Si atoms are positively charged.

In the crystallization process, SiC establishes a special case of polymorphism called Due to the fact that SiC does not exist in liquid phase as can be seen from the phase diagram of Si-C system (Fig. 2.6), the conventional method used to grow Si and other semiconductors based on the solidification from liquid phase cannot be applied. Alternative growth techniques are then necessary to realize the growth of SiC. Concerning the bulk growth process, a variety of growth techniques has been studied, such as, the Lely method, Modified-Lely method, High Temperature Chemical Vapor Deposition (HTCVD), Continuous Feed-PVT (CF-PVT), and Top Seeded Solution Growth (TSSG) processes. Those techniques can be grouped into two main routes: vapor phase route and liquid phase route. Overview on the growth of SiC single crystal for both routes focusing on the new emerging processes was discussed in [11].

Growth from vapor phase

One of the big steps in the advancement of the SiC growth process can be attributed to the invention of the Lely growth method in 1955 [12]. The schematic representation the process design and the temperature profile can be shown in Fig. 2.7a. The SiC powder source was put between the wall of graphite crucible and the porous graphite leaving the empty space at the center of the crucible. The temperature distribution is adjusted so that the center of the cylindrical growth chamber has a minimum temperature. Thus at elevated temperature, the vapor species sublimed from the source will travel toward the center. The porous graphite will then behave as a center for the nucleation. Due to the fact that the growth is performed close to the equilibrium and the growth can be assimilated to spontaneous nucleation, the growing crystals have very low defect density and are almost free of micropipe. The main drawback of this growth method is the limited and random size of the crystals.

In 1978, Tairov and Tsvetkov [13] developed the so called modified Lely growth method also known as Physical Vapor Transport (PVT) or seeded sublimation method. The concept is to sublime the SiC powder source at high temperature, the Si and C containing gas species are then transported by the physical mean (mainly diffusion) and re-crystallize at the SiC single crystal seed having a slightly lower temperature than the source. The schematic representation of this growth method is shown in Fig. 2.7b. The growth is typically carried out in a quasi-close low porosity graphite crucible in argon ambient with low pressure and in the temperature ranges 1800-2600 °C. Compared to the Lely method, the extensive spontaneous nucleation along the inner crucible wall can be suppressed using this method. The main challenge is to control the temperature and the gas species concentration fields, which has been achieved through the support of the process modeling [14][15][16][17][18][19][20]. The control of the growth at the seed can then be achieved and the control of polytypes is possible to some extent [21][22][23][24][25][26][27][START_REF] Schulze | Controlled Growth of 15R-SiC Single Crystals by the Modified Lely Method[END_REF]. Very recently, a complete macroscopic approach, including experiments, characterization, and modeling, to the nucleation and propagation of foreign polytype inclusions has been reported [START_REF] Tsavdaris | Macroscopic Approach to the Nucleation and Propagation of Foreign Polytype Inclusions during Seeded Sublimation Growth of Silicon Carbide[END_REF]. The PVT technique has attracted exceptional research efforts and is the most mature process to date. Its abilities are to produce, for example, large size boules (150 mm in diameter) with low dislocation density (EPD=2912 cm -2 , TSD=302 cm -2 , BPD=324 cm -2 ) [30], and micropipe-free single crystalline ingots [31]. The batch nature of the process, however, results in the intrinsic limitations of the PVT process. The evolution of the geometrical parameters, namely polycrystalline powder properties and ingot dimension, contribute to the excessive changes in the growth behavior [32].

With the aims to grow high quality and high purity SiC at high growth rate, Kordina et al [33] presented the HTCVD process as an alternative method to grow SiC boules. The HTCVD geometry is considered as a vertical type growth reactor where the precursors like SiH 4 and C 3 H 8 diluted in a carrier gas used as Si-source and additional C-source, respectively, are fed upward through inlets placed at the bottom of the crucible and the exhaust at the top (Fig. 2.7c). Even if the technique resembles the Chemical Vapor Deposition (CVD) technique, the growth process considerably differs from the CVD process. The chemistry of the HTCVD process requires a separated inlet for Si-and C-containing precursors and an adapted temperature distribution. The purpose is to control a 'SiC cluster cloud' formed by the homogeneous nucleation acting as a 'virtual SiC solid source' for the sublimation step [34]. The main advantages of this technique are the continuous supply of the material, the direct control of Si/C ratio, and the excellent control of purity and doping. On the other hand, there are several difficulties using HTCVD process. The surface morphology of the grown crystal is very sensitive to the growth condition, especially during the temperature ramp up. Moreover, the premature cracking of SiH 4 can lead to either parasitic depositions of Si at the cold parts of the reactor or a blocking of the inlet [33]. Since the last ten years, the liquid phase route for the growth of SiC single crystal has attracted again attention. The potentials of using the liquid phase techniques for the growth of bulk SiC crystal were discussed [40]. In order to avoid the thermodynamic problems related to the stoichiometric melt, the use of non-stoichiometric solutions using, for instant, Si as a solvent and C as a solute has become an alternative way toward the growth of SiC from liquid phase.

Examples of the liquid phase growth set-ups for SiC are shown in Fig. 2.8. On the perspectives of large diameter and reasonable boule length, the TSSG method seems to be the most feasible approach due to the intrinsic limitations of the other two methods concerning especially upscaling [40]. In the TSSG process, the graphite crucible does not only act as a container, but also as the carbon source. The carbon dissolved from the crucible walls will be transported by both diffusion and convection to the SiC seed crystal attached to the seed holder where the temperature is relatively lower. There are some difficulties using this method. First, the manipulation of liquid silicon at high temperature is not very easy due to the high Si vapor pressure and the extreme reactivity with any container material (other than graphite and SiC).

Second, several convective phenomena, such as buoyancy convection, Marangoni convection, forced convection, and electromagnetic convection, can possibly occur depending on the process conditions. Thus the control of the fluid convection in the melt is one of the main challenging issues. This will, in combination with the control of temperature profile, directly affect the growth front and its stability. In fact these problems are associated with the low C solubility in liquid Si as the latter implies using high temperature to increase the solubility limit and/or introduce some convection to enhance the mixing. An alternative approach is to use another solvent such as adding a metal to Si to increase the C solubility. However, another problem will arise concerning the purity of the crystals and especially the contamination from the solvent.

Apart from such difficulties, the liquid phase route shows many excellence advantages: the elimination of harmful defect such as micropipes [41,42], the enlargement of the crystal diameter and length [43], and the high structural quality [44].

Process modeling

In the previous section, it can be seen that many challenging topics in the growth process are related to the precise control of the physico-chemical parameters and their distribution in the Thanks to the pioneering works on the process modeling of SiC single crystal growth of Hofmann et al [45] and Pons et al [46], many research groups started to work in parallel to develop the models and/or codes for the mass transport, chemical reactions, and mechanical stress in the crystal [21,. Nowadays, the modeling of heat and mass transfer has reached the state of the art [START_REF] Chen | Progress in modeling of fluid flows in crystal growth processes[END_REF][START_REF] Chen | Optimization of the design of a crucible for a SiC sublimation growth system using a global model[END_REF]. On the other hand, different approaches and models have been proposed for the chemical reactions, especially at the solid-vapor interfaces [21,[START_REF] Segal | Transport phenomena in sublimation growth of SiC bulk crystals[END_REF][START_REF] Chen | Kinetics and modeling of sublimation growth of silicon carbide bulk crystal[END_REF][START_REF] Danielsson | Predicted nitrogen doping concentrations in silicon carbide epitaxial layers grown by hot-wall chemical vapor deposition[END_REF][START_REF] Karpov | Analytical model of silicon carbide growth under freemolecular transport conditions[END_REF][START_REF] Nishizawa | Growth and Doping Modeling of SiC-CVD in a Horizontal Hot-Wall Reactor[END_REF].

Those models require a certain degree of connection, depending on the model, with the thermodynamics. The model for predicting the dislocation density has also been developed based on the thermal stress calculations [17] while the most recent achievement is the implementation of Alexander-Haasen (AH) model to describe the dynamics of plastic deformation in SiC where the effects of the cooling time on the final dislocation density, residual stress, and stacking faults can be studied [53].

The reliability of the modeling results can be strongly dependent on the database used in the computation. Those include, for example, the thermal conductivity, electrical conductivity, and emissivity for the solid components, and thermo-chemical data for the vapor species.

Actually, this issue may become very problematic, especially for the process modeling of SiC. This is caused by the lack (or not very precise value) of the materials properties data under such high temperature normally used in the growth. As a consequence, this factor defines whether the model can be considered as to provide quantitative results or qualitative trends [11]. In addition to such issue, some other issues such as the simplification of the physical system as well as the choice of the boundary conditions can affect the modeling to some extent. These issues are caused by the fact that it is impossible to model the entire universe. We can model only the physical phenomena that we are interested in, in a defined space (and time). The other things outside will be represented by the boundary conditions of the system. Thus, it is quite important to keep in mind that the absolute or exact solutions are not very realistic and cannot be expected from the process modeling even if the quantitative results are obviously the best results we would like to know. However, providing that the databases are not "completely wrong", the numerical modeling has a strong potential to give the qualitative trend which is very useful for the development and understanding of the process. In some cases, the modeling results can even be considered on a quantitative footing.

Purposes of this thesis work

Besides many successes in applying the numerical modeling to describe the SiC growth processes, there are some strong major challenges which are not considered in a satisfactory way even in the most mature process like PVT. There are (at least) three main open questions for better understanding of the PVT process [START_REF] Chaussende | Open Issues in SiC Bulk Growth[END_REF].

i.) How to describe correctly the heterogeneous reactions at the solid-vapor interfaces?

ii.) Is there any link between chemistry and crystal polytypes?

iii.) How the doping incorporation can be described?

Concerning the TSSG process, the control of the fluid convection and the growth front stability are the main issues toward the development of high quality SiC crystal, especially for electronic devices applications. Thus in this thesis work, we are aiming to use the numerical modeling as a tool in order to tackle those challenging fundamental and technological issues, both for the industrially used PVT process and for the emerging TSSG process. The numerical modeling results will be compared and discussed either with our experimental results if available or with the data available in the literature.

It is worth noting that in this thesis, we are not aiming to develop the code for the numerical modeling since the commercial software, such as, COMSOL multiphysics [80], CFD- Reactor [84], or the open source such as, OpenFOAM [85], are available. They all can treat modeling in process engineering, either specifically for the SiC growth (Virtual Reactor) or for general (coupled) physics phenomena. In this thesis, the numerical modeling was performed using COMSOL Multiphysics software. Such software package is working based on the Finite Element Analysis (FEA). Thus, an overview of the Finite Element Method (FEM) will be given in the following chapter including an example for solving the heat transfer in steady state.

ACE [81], ANSYS Fluent [82], Abaqus [83], Virtual

Chapter 3 Finite Element Method (FEM)

The SiC growth processes involve many physical and chemical phenomena [1]. The system of equations contains unknown multivariable functions as well as their partial derivatives; the whole process system is thus a system of partial differential equations (PDEs). In fact, PDEs are naturally governing science and engineering applications due to the complex balance equations, such as transport of mass, momentum, species, and energy. The whole class of these transport phenomena usually expressed in the integral equations over the domain results in the PDEs in the continuum approximation. Thus, the system to be solved is in reality continuous and the exact solutions can be solved mainly in the oversimplified situation. With aim to approach the true continuum solutions without oversimplification, several numerical methods of discretization have been proposed for solving such PDEs, for example, the Finite Difference Method (FDM), the Finite Volume Method (FVM), and the Finite Element Method (FEM). Even if the approaches in solving PDEs are different, these methods share important characteristics.

First the continuous domain is discretized into a set of discrete subdomain called element in which their behavior is specified by the finite number of parameters. Moreover the approximated functions on each element meet the conditions of continuity between various elements of the field. Advantages and disadvantages of using these numerical methods are summarized in Table 3.1 [2]. The FDM method is relatively easy to be programmed but its ability is limited to a simple geometry and regular grid (element). The FEM is much better to treat correctly the interface problems, even if this method is less trivial for fluid dynamics than FVM. Since the crystal growth is associated mainly with the problem of boundary conditions, the FEM method is more adapted in this current PhD topic. Throughout this thesis, the FEM will be used as a tool for approximating the solutions to the boundary value problems of the PDEs system. An overview of the FEM and examples of solving heat transfer problem using this method will be presented in this chapter [3][4][5][6].

  u A    u B 
The unkn

    0 p Lu u A          d w d w . ( 3.4)
If (3.4) is satisfied for all w, then the differential equations (3.1) must be satisfied at all points in the domain. In fact by inserting the approximation (3.3) into (3.1), we obtain the residual or error of the differential equation. Thus Eq. (3.4) can be considered as a weighted integral of such residuals and this approximation may be called method of weighted residual. This method aims at minimizing the residual by multiplying with the test functions and integrates over the domain.

If (3.4) is satisfied for any w, the residual will approach zero and the approximated solution will approach the exact solution. The accuracy of this method is dependent on the choice of the test function. One of the frequently used methods is the Galerkin method. This method assume the test functions to be similar to the shape functions such that w j = ϕ j . Another example of the method to obtain the integral form is the Ritz method. In such case the problem is characterized by the energy function and the energy minimization technique is used. The integral forms obtained by two techniques are indeed equivalent. Since the energy function is not always known, the weighted residual method is more generally used.

To illustrate the transformation procedure, we may consider the system that is governed by the Poisson's equation

  0 p u p Lu u A 2       . (3.5)
Equation (3.4) can be rewritten as

  0 p u 2       d w . (3.6)
It is possible to perform integration by parts (or other methods to shift the derivatives) on Eq.

(3.4). From the relation

  u v v u uv         , (3.7)
we can re where K ji is called the stiffness matrix and f j is called the load vector. The superscript e defines the element (locally based) matrix. For the N DOF degree of freedom, (3.14) is equivalent to a system of N DOF equations. Such system will then be solved by appropriate numerical methods.

     Then app      to (3.8), w     w
          i u i i an be rearran       d i j ve integral f e procedure s        
            i u i i         d i u n i ct
     d j p n      d j p . approximat iscrete syste        
The size of the matrix is (N DOF ×N DOF ) which depends on the number and order of elements geometrically used to discretize the domain. The larger system size requires the more powerful computer for reasonable computation time.

Example of heat transfer in steady state

This example relates to the calculation of the temperature distribution in the steady state.

The diffusion equation of heat can be written by the formulation representing the conservation of energy

  Q T k      (3.16)
with the Neumann type boundary condition written by

  amb T T h T k n T k          n , (3.17)
where T is the temperature (field variable), k is the thermal conductivity, Q is the heat density, and h is the exchange coefficient. The weighted integral for residual of Eq. (3.13) can be written:

               d Q d T k j j   . (3.18)
From the relation (3.7), (3.18) is simplified to

                     d Q d T k d T k j j j    . (3.19)
Applying the divergence theorem (3.9) to (3.19), we obtain

                  d Q d n T k d T k j j j    . (3.20)
The boundary condition (3.17) can be directly subtract to Eq. (3.20) to achieve

                    d hT d Q hTd d T k amb j j j j     . (3.21)
This is a general formulation for the volume and surface integrals which can be also applied to the case where a domain is composed of materials with different physical properties (Fig. 3.5)

where the continuity relations are: .23) In this case, the integral forms (3.20) become

2 1 T T  and (3.22)     0 2 2 2 1 1 1           n T k n T k . ( 3
                                d Q nd T k nd T k d T k j ext ext j j j 2 1 2 1 2 1 2 1     . (3.24)
Due to the continuity relations at the interface (3.23), the integral of the area along the interface

(ᴦ 1 ᴖᴦ 2 ) is canceled, 0 2 1 2 2 2 1 1 1 2 1                           d n T k d n T k nd T k j j j    . (3.25)
It can be seen that FEM is a method that meets the continuity of flow at the interfaces of the materials with different physical properties. Eq. (3.20) remains valid and the surface integral can be applied on the border of the study area (ᴦ 1ext ᴗᴦ 2ext ) to obtain

    2 1 Next task volume/s condition    v e e k
Then a li element The next step is to define the stiffness matrix K ji of the system. Considering an example of six discrete mixed type (triangle and rectangle) elements shown in Fig. 3.6, each element has its own identifying number and nodal connections corresponding to the coefficients of coupling between the nodes. Concerning the case that the properties can be found in the global coordinates, the stiffness matrix/load vector can be then written as a summation of the global matrix/vector written for each element (Fig. 3.7a). Each separated shaded area corresponds to either a single coefficient or submatrix of type K ji if more than one quantity is considered at nodes. The resulting matrix/vector assembly in the global nodes numbering is shown in Fig. 3.7b with the N DOF × N DOF = 9×9 overall matrix dimension. This leads to a linear system with N DOF unknowns (T value at the nodes) and N DOF equations of the forms

  T T which re        v e i for all j  Fig. 3.5 D      d T k j  k is to discr surface by n (3.17) to (3      e j d T k  inear combin    1 i i i T  , sults in the f      e e j k     DOF N , 1  Domains of       
      j i ji S T K  . ( 3.30)
The stiffness matrix and load vector can be also called the matrix and vector integrant, respectively. These integrants are calculated numerically using the Gaussian quadrature rule which approximates the definite integral by replacing the integral of the function by a weighted sum of function values at specific points in the domain:

         k k k k k e e f d z y x f     , , , , . (3.31)
Weighting function (ω k ) and Gauss points (α,β,γ) are known on all types of element: 2D (triangle and quadrilateral), and 3D (tetrahedral, prism, pyramid, hexahedral).

Matrices generated by the discretization PDEs can be very large. However, most of the coefficients of this matrix are zero (see for example in Fig. 3.7b). This type of matrix is called a "sparse matrix". Methods for solving linear systems of sparse matrices are thus used to limit the number of coefficients to be stored in the memory. The matrices are then represented on the basis of knowledge of the position of the non-zero coefficients (shaded area in Fig. 3.7b). The non-zero coefficients are confined within a so-called matrix band. It can be seen that the matrix is symmetric thus only the matrix elements in the upper half of the band above the diagonal are stored in the memory. In fact, the numbering of the nodes is really important to reduce the bandwidth of the matrix therefore reducing the memory needed for the storage of coefficients. The field of materials development and engineering processes is in reality involved in the coupled physico-chemical phenomena. Such coupled phenomena practically occur in the industrial environment, for example, electromagnetism and heat transfer in the induction heating process, and the fluid dynamics and the chemical species transport in reacting flow process.

These kinds of problems are necessary to be considered in a coupled regime. The coupled phenomena are normally implemented using the two mains techniques: strong coupling and weak coupling. The former one is performed by simultaneously solve all of the coupled equations within a single resolution matrix. This method provides more accurate results with a higher cost of computational time and memory. The latter one is done by solving the equations sequentially. The overall system may be decomposed into several subsystems. This method can benefit from the computational speed using the specialized solver but with some risk concerning an instability issue. Giving the magneto-thermo-hydrodynamics problem as an example of the coupled phenomena, the number of the unknown per node is 9 (3 components of magnetic vector potential, scalar electric potential, temperature, 3 velocity components, and pressure). The matrix system at each nodal point for both strong and weak coupling techniques can be seen in Fig. 3.9.

For the system of NBN nodes, the size of global matrix system calculated from the number of degrees of freedom square (N DOF ) 2 =(9×NBN) 2 in the case of strong coupled is much larger than the weak coupled one (4 2 +1 2 +4 2 )×NBN 2 . If the physical properties are linear, the system of (3.30) can be solved either by direct methods (MUMPS, PARDISO, SPOOLES) based on LU decomposition (LU factorization) or iterative methods (conjugated gradient, GMRES, etc.).

A nonlinear system

If the system (3.30) is nonlinear, for example the thermal conductivity is dependent on the temperature, the equation of heat transfer or operator L in (3.1) is nonlinear. This case is also a general case for fluid mechanics. In such case, the nonlinear solver based on the Newton's method is used. Thus the solution will be obtained iteratively where the initial conditions and the maximum number of iterations must be specified. The formulation representing the residual in the nonlinear system with L+1 number of iteration can be written as This technique is known to slow down the calculation but gives a more secure convergence. In fact, the method using the Jacobian matrix is not the only possibility but a fixed point method can be also used to directly iterate the system (3.29). However the Jacobian matrix is used to ensure the good convergence resolution method. In some industrial software, the Jacobian matrix is either symbolically or numerically assembled.

  0 1   L u F . ( 3
In the scope of this thesis, we limit the modeling only to steady state calculations.

However, it is worth noting one of the important features of the transient FEM calculations. The term that contains a time derivative is developed either by using an implicit or explicit finite difference. The former one is better to avoid the stability problems while the latter is better in case of a large deformation. A system of the type (3.30) can be formed while the solution at the initial time and the time step has to be specified. The transient problem can be also a nonlinear problem, thus the iterations technique of the nonlinear system has to be imposed.

In summary, the FEM method has a strong capability to solve the PDEs governing the physical phenomena associated with the growth processes in a complex geometry. Such process characterization tool, however, requires proper definition of the physical systems and the boundary conditions. We will consider step by step the following systems related to the SiC growth processes. The induction heating and heat transfer system will be discussed in the next chapter. Then the system of mass transport concerning the PVT method will follow. Finally, the full TSSG process will be discussed.

Chapter 4

Induction heating and heat transfer

The thermal conditions in growth system were studied by means of coupled induction heating and heat transfer in the entire region of furnace. All the calculations were conducted by Finite Element Method (FEM) in a two-dimensional axisymmetric space-dimension since the growth reactor is shaped axisymmetric. In this chapter, the physical model and the boundary conditions are presented, followed by the effect of the semi-transparency of SiC. Finally, the modeling results and validations will be discussed.

Physical model

Induction heating

The practical method for high temperature heating of crucible is by electromagnetic induction. In this RF induction heating regime, the electromagnetic phenomenon is described by the Maxwell's equations which provide the distribution of electromagnetic field. The induced eddy currents are generated according to Lenz's law. The Ohmic resistance caused by the current dissipates energy and the crucible is consequently heated by the Joule's effect. The differential form of Maxwell equations are written by [1] ,

t       B E (4.1) , t f       D J H (4.2) , f     D (4.3) . 0    B (4.4)
If low frequency (f < 1 MHz) is used in the experiments [2], the currents and the electromagnetic field vary slowly and the dimension of the growth reactor is small compared to the wavelength.

In this case, the Maxwell's equations for quasi-static approximation under the assumption that 0    t D is valid. The elimination of the displacement current density term reduces the computational cost. Under the assumption that the magnetic field is time-harmonic with an angular frequency ω, there is however no computational cost for including the displacement current in Ampere's law. Thus the Maxwell-Ampere's law (4.2) can effectively be used [3].

From (4.4), the magnetic field B can be formulated by the magnetic vector potential A such that

A B    , (4.5)
From Eq. (4.1), the relation becomes

. 0             t A E (4.6)
The expression for electric field can be written in terms of a vector potential and scalar potential in such a way that it provides a zero curl, thus

. t V      A E (4.7)
The current density is proportional to the force per unit charge f and the electrical conductivity.

In the system heated by induction, the force that drives the charge to produce current is a current caused by electromagnetic force and the externally generated current density J e . Thus the free current density can be written

  . e e f J B v E J f J         (4.8)
The total current density J can thus be written as

  e d f t J D B v E J J J           . ( 4.9)
From (4.5), (4.7), and (4.8), the Maxwell-Ampere's law (4.2) for linear materials can be written

  e r r t V t t V J A v A A A                                          0 0 1 . (4.10)
This forms an equation of two potentials: scalar potential V and vector potential A. In fact, these two potential are not uniquely defined from the electric (4.7) and magnetic (4.5) fields since it is possible to impose extra conditions on V and A providing that E and B are not affected. For example, substituting a set of potentials

t V V      ~and (4.11)     A A ~.
(4.12) will give the same electric and magnetic field

  . t V t t V                        A A E (4.13)   A A B          (4.14)
Thus, it is necessary to add the constraint (gauge transformation) to ensure the unique solution.

For a particular gauge transformation of the form

    V i , providing that V i     A A ~ and (4.15) 0 ~ V , ( 4.16)
equation (4.10) can be reduced to the magnetodynamic equation for the vector potential:

    e r r i J A v A A                   0 2 0 1 , (4.17)
The work done per unit time per unit volume can be expressed as

. t t          A A J E (4.18)
By assuming the relaxation time for heat dissipation is much shorter than the time scale for heat conduction, the volume density of the generated heat power can be computed from the mean value of (4.18). Thus, power delivered per unit volume can be formulated by

  * 2 2 1 A A J E        Q . ( 4
      . th p p Q T k T C T t C             v (4.20)
The local heat generation in unit volume Q th is obtained from the electromagnetic induction heating calculated by Eq. (4.19). Eq. (4.20) is written based on the major assumption that the heat production and consumption (exothermic and endothermic) according to chemical reactions in the gas mixture does not contribute much to the heat transfer for the sublimation process studied in this thesis. Moreover, the quasi-steady state is assumed during crystal growth process.

Therefore, the governing stationary state equation for thermal energy simplifies to

    th Q T k T p C        v  . ( 4.21)
The temperature distribution profile in the system can be solved from Eq. (4.21). The first term on the left hand side is related to the translational heat transfer. Normally the solid parts in the growth system do not move so this term becomes zero. In the gas phase especially in the growth chamber, the term for heat transfer by convection does not vanish. The computation of the convective heat flux requires the knowledge of the mean velocity of the vapor species, which can be obtained only from the modeling of fluid dynamics and mass transport. Thus, the iteration steps are needed for the computation. The basic iteration procedure is associated with the following steps:

1. Guessing the initial value of the fluid velocity (typically we impose the case of no convective heat transfer so the fluid velocity is zero), the temperature distribution can be solved using the induction heating and heat transfer modules 2. Using the temperature profile computed from the 1 st step, the fluid velocity can now be solved using the reacting flow module 3. Performing the calculation as done in the 1 st step of the procedure using the updated fluid velocity computed from the 2 nd step

The iteration process will be done until the convergence of the calculations is obtained; usually after several iteration steps. However if the heat transfer by radiation is taken into accounted, it becomes the dominant contribution at the solid-vapor interface. The similar conclusion was reported from the estimation of the heat transfer quantities [4]. The temperature distribution in the growth chamber without and with convective heat transfer is shown in Fig. 4.1a and Fig. 4.1b, respectively. The convective heat transfer modifies slightly the temperature distribution in the gas phase. The upward convection flow pushes the isotherms toward the seed crystal. The changes in the isothermal shape when taking convective heat transfer into account is in agreement with [5] resulting in the steeper axial temperature gradient toward the seed crystal.

The temperature at the solid-vapor interfaces are very slightly modified as can be seen for example for the radial temperature profile along the SiC seed crystal (Fig. 4.2) and overall temperature gradient is not affected by the gas convection in this work. This is in contrast to [5] due to the fact that the gas velocity was amplified in order to numerically demonstrate the effect of the convective heat transfer in such work. In the scope of this thesis, only the heterogeneous reactions at the solid-vapor interface are considered while the chemical reactions are neglected in the gas mixture. Since the exclusion of convective heat transfer does not affect critically the computational accuracy, the convective heat transfer can be neglected in order to reduce the computational cost. It is worth noting that neglecting the convective heat transfer is one of the major assumptions in many numerical simulation works [2,4,6,7].

Electromagnetic boundary conditions

In the simulation of open boundary problem, some mathematical transformation is needed. A set of boundary conditions must be imposed for the magnetic model. The zero vector potential is applied on the limiting box corresponding to infinity distance from the coil. Moreover, due to the axisymmetric design of the reactor, a zero normal induction field is defined at the symmetry axis.

Heat transfer boundary conditions

A set of boundary conditions are formulated for the heat transfer problem. First, the zerotemperature gradients is imposed at the symmetry axis. Second, the thermal flux boundary conditions related to the radiation are fixed at boundaries of the cavities. The net inward flux is calculated from the difference between the irradiation and the radiosity written as

  J G T k       n , (4.22)
where G is the total arriving flux or the irradiation. The generalized equation for the irradiative flux including the mutual irradiation m G coming from other surfaces and the ambient irradiation amb G can be written as

4 amb amb m amb m T F G G G G      . ( 4
            d J G m 4 r r n r n  , (4.24)
where the local radiosity is projected onto point x when computing the heat flux arrived from x .

The second term represent the ambient irradiation where amb F is the ambient view factor. This function describes the proportion of view from each point covered by ambient conditions. Its value can vary from zero to unity at all point depending on the geometry. Thus the ambient view The thermal flux boundary conditions associated with the radiative phenomenon of the reactor outer wall are implemented by assuming that the surrounded ambients behave like a blackbody where its emissivity and reflectivity equal to unity and zero, respectively. Moreover, the mutual irradiation equals to zero and the ambient view factor equals to unity. Therefore, the boundary condition for the surface-to-ambient radiation becomes

    4 4 T T T k amb e         n . (4.29)

Semi-transparency of SiC

In general an opaque body was assumed for all materials used in the growth system within a fixed range of wavelength

  2 1 ,    
. This means that the radiative heat transfer was taken into account only at the solid-vapor boundaries where the non-reflected waves are completely absorbed and emitted. However, some materials exhibit a semi-transparency behavior. In such case the absorption process does not take place only at the solid-vapor interface but parts of the energy were also absorbed and transmitted in the material medium itself. This phenomenon has been shown to have a strong influence in the crystal growth process [8].

Concerning the process modeling for the growth of SiC, the effects of semi-transparency were reported [5,9,10]. The overall temperature field is slightly affected by the semi-transparency while the local effect on the crystal and its vicinity is significant. First, the influence on the temperature distribution in the crystal may affects directly the thermal stress and hence the dislocation density. Moreover the temperature distribution close to the crystal-vapor interface affects the shape of the growing crystal [10]. In this section, we are considering the method to treat SiC as a semi-transparency material and then we will discuss the effects on the growth process and the validation of this model.

In order to compute the heat exchange in the semi-transparent medium, the band-energy method is used [8]. The schematic representation of the band-energy method is shown in Fig. 4.4. The system domain containing the semi-transparent medium can be decomposed into two enclosures. The first enclosure corresponds to the range of wavelength 1  where SiC is opaque so the heat fluxes emitted from surface (1) interact with the surface (2). The second enclosure corresponds to the wave in a range 2  where the SiC is transparent so the heat fluxes emitted from (1) radiation can be w

  Q r n   Q r   Q r n
where q general contribut the opaqu to be e SiC rystal are better described by the model with opaque SiC. Thus, it is very important to discuss about the validity and limitations of the semi-transparent model. First of all, we are considering the monochromatic plane wave of wavelength λ having an initial intensity I 0 and propagating through the material of thickness h. In the simplest case where the reflection and dispersion are neglected, the wave transmitted through the material will be attenuated due to the absorption.

     T k n      T k n      T k n 1
The relation between the initial intensity and the outgoing intensity can be written [8]  

z I I h       exp 0 , (4.34)
where the absorption coefficient   can be correlated to the skin depth (the distance that it take to reduce the amplitude by factor of 1/e) by

    1  . (4.35)
We may define the spectral optical thickness for which its magnitude is dimensionless such as

h h       1   . (4.36)
The material is opaque if

1    whereas it becomes transparent if 1   
for the wavelength λ. Regarding the band-energy method for treating SiC as a semi-transparent body, we considered that SiC is opaque for some range of wavelengths (such as 1  ) but transparent for others (such as 2  ). Moreover, the transmitted wave is assumed to travel through to crystal without interaction with the crystal itself. Such assumption is accurate only if the ranges of wavelength in which   is of the same order as unity is negligible. This is the first limitation of the model. The second consideration is the dependency of the absorption coefficient of SiC on the temperature, polytype, wavelength, and the doping level [11][12][13][14][15]. There are some discrepancies in those experimental observations so it is not very easy to draw any strong conclusion here. For the optical band, the absorption coefficient having very low values was mostly observed in a relatively pure crystal (low doping concentration) while increasing when the charge carrier concentration increases [12]. Most of the non-zero value of   varies from few cm -1 to several thousand cm -1 . So let's consider here the two extreme cases, 1    cm -1 (lower limit) and 55 1000    cm -1 (upper limit). The SiC crystal will become transparent if 1    meaning that the crystal thickness should be at least one order of magnitude lower than   , i.e., thinner than 0.1 cm and 0.0001 cm for the lower and upper limits, respectively. If the SiC crystal used as a seed is a bit less than 0.1 cm (as used in a standard growth process), the upper limit predicts that SiC is always opaque. On the other hand, the lower limit predicts that the transparency should affect only the initial stage of the growth (few hundred microns) while the crystal becomes opaque after. This requires also that the seed and the newly grown layer must be relatively pure.

Otherwise SiC crystal is always opaque.

The more simplified consideration can be analogous to the poor conductor ( ε

  

). In such case, the skin depth is independent of the frequency and can be calculated from

[1]     2  . (4.37)
It is found that the skin depth is very low in SiC (few tens of microns). Moreover this relation shows that the skin depth is reduced if the electrical conductivity increases. Since the electrical conductivity can be increased by increasing the doping concentration, the crystal can become opaque in a highly doped condition. This is the case in our PVT growth system since the (nitrogen) doping level reaches high 10 18 to low 10 19 cm -3 even in an unintentional doping experiment. Another argument for the invalidity of the semi-transparent model in our PVT process can possibly relate to the temperature dependence of the energy bandgap in semiconductors. In most of the semiconductor materials, the bandgap decreases as the temperature increases. The use of very high temperature during the PVT growth process can drastically reduce the bandgap in SiC crystal. As a consequence, SiC crystal becomes an opaque body with respect to the radiative heat fluxes at the growth temperature. For these reasons, we will assume that the SiC is an opaque body in this thesis.

Modeling results and validations

In order to validate the model, the modeling results have to be compared with the experimental results. Since the computation is mainly performed for the steady state of the growth, then the best case should be the comparison between the modeling results and the experimental observation during the growth experiment. However, this is not always very easy due to the limitation in the real time measurement especially for the validation of the induction heating and heat transfer discussed in this chapter. The creation of an open space for measurement or an introduction of the external probe would substantially disturb the system. Thus in this section we are aiming to verify the modeling with the real time measurements if they are possible otherwise the after growth observations will be used to compared. The modeling results include the I-V characteristic of the induction coil and the effect of coil position. A list of physical properties of materials used in the simulations can be found in appendix A.

I-V characteristic and measurement

In the set-up for SiC crystal growth, the control of power supply and geometry of the reactor is very crucial, for example, for the optimization of the power consumption and controlling the thermal gradient in the growth chamber where the mass transport takes place. The topological connection of the power generator, capacitors, and the reactor is shown in Fig. 4.10a.

The equivalent RLC parallel circuit with resistance in series with inductor is drawn in Fig. 4.10b.

In such RLC circuit, there are two types of power losses: active losses due to the resistive components in the system and the reactive losses due to the inductance. The apparent power S (VA) defined as a measure of alternating current power is the vector summation of active P (Watts) and reactive Q (VAR) powers, and it is a product of root-mean-square of circuit voltage and current written by

2 2 Q P I V S rms rms    . (4.38)
Since the growth of bulk SiC needs rather high temperature of around 1700-2300 °C, the power generator must be able to provide enough active power to heat the crucible up to the such temperature. Moreover, the capacitance is necessary to be integrated in order to generate enough reactive power consumed by the inductor. By mean of an adaptation of the impedance, the optimum design of the coil can be obtained when the impedance of the equivalent R-L-C circuit is near the impedance of the power generator so that the system functions as close as possible to the resonance condition. This optimal design can be obtained by modifying the geometric parameters of the coil such as, the number of turns, the diameter and the height, and the capacitors (frequency). Since the circuit topology of our reactor is a parallel circuit, the parallel Their values and distributions are dependent on the geometry. In fact, this is one of the reasons why the crystals growers face the problem called "reactor dependency". With the same set of input or controlled parameters, the observed results may be differing to one another due to the effect of reactor design. This will lead to the discrepancies concerning the interpretation of the experimental results as well as the conclusions. On the contrary, the numerical modeling of the growth process grants us the possibility to enter the physico-chemical parameters (Fig. 4.15). Then those parameters can be compared directly with the observations, making used of analytical or theoretical approaches. Such comparisons are geometry independent, thus the problem of reactor dependent can be minimized. This point clearly shows the potential of the numerical modeling in the growth processes.

Effect of coil position

The adjustment of the coil position is one of the simplest ways to control the temperature gradient inside the growth chamber. This is much easier than dealing with the crucible geometry.

In order to study the effect of the coil position on the temperature distribution, three coil positions are simulated and compared. The reference coil position is chosen to be similar to the geometry shown in Fig. The curvatures (or shapes) of these streamlines follow the distribution of the isothermal contours.

In the growth cavity the average temperature gradient computed by the temperature difference between top part of the source (T source ) and the seed crystal (T seed ) at the symmetry axis is almost 9 K/cm. If the coil position is shifted upward by ∆z, the hot point shifts upward by almost double the distance of ∆z. The temperature distribution and the temperature gradient streamlines are modified mainly in the powder source. First, the temperature gradient streamlines direct both upward and downward directions in the powder source (Fig. 4.16b-c). Moreover, T source decreases faster than T seed (Fig. 4.16d), thus the average temperature gradient decreases from almost 9 to 6.5 and 4.5 K/cm by shifting up the coil position by 1 and 2 cm, respectively. [20], the combination of these two mechanisms contributed to the recrystallization phenomena on the top part of the source. The last possible mechanism could be the back flow of the SiC vapor species caused by the re-sublimation of the seed. This will be discussed in details in the next chapter. As a consequence, the strongly graphitized area is expected to be located at the bottom of the crucible for the standard coil position. When the coil position is shifted up by 2 cm, the strongly graphitized area locates close to the powder-vapor interface. In this case, part of the sublimed species will diffuse through the porous powder to top part of the powder and transport to the seed (and other relatively cold points) while another part will diffuse down to the bottom of the crucible due to the temperature gradient. Since the sublimation rate is high on the top part of the powder and there is no supplied SiC vapor from the bottom part, the powder becomes graphitized in this area. This can limit the growth of SiC due to two possible effects. First, the growth rate will be decreased due to the reduction of the temperature gradient. Furthermore, the efficiency of the growth concerning the powder consumption may be decreased since some amount of the powder (downward diffusion) does not contribute to the crystal growth. This effect of increasing the coil position is confirmed by the experiments where the growth rate decreases by more than a factor of three, i .e. from around 200 μm/h to 60 μm/h. The theoretical analysis of the mass transport in the powder charge was proposed based on the Darcy's law where the evolution of the powder concerning the porosity and graphitization degree in a longterm growth of bulk SiC crystal can be predicted [21,22].

Conclusion

The coupled induction heating and heat transfer were implemented in the growth furnace.

Their physics represented by the electromagnetic equations, heat transfer equations, and momentum and continuity equations were numerically solved with the defined boundary conditions. The heat transfer by conduction was shown to be negligible and SiC crystal was treat as an opaque body. This coupled model provides the temperature distribution in the reactor which is necessary for the mass transport modeling. Two types of modeling results and validations were shown. First, similar I-V tendency between numerical and measurements was obtained. Second, a direct comparison of the diffusion path of SiC species in the source powder and the temperature gradient streamlines was made regarding the effect of the coil position.

Chapter 5

Mass transport

The mass transport of the vapor species in the gas mixture is coupled between the species transport and the fluid flow making use of the temperature profile computed from the induction heating and heat transfer as describe in the previous chapter. In this chapter, the physical model is presented followed by the boundary conditions concerning both comparative and coupled approaches. The effect of the choice of thermodynamic databases using the former approach will be discussed. Then the chemical model concerning the interaction between vapor species and graphite crucible will be given. Finally the general modeling results and validations will be presented.

Physical model

In order to model the fluid flow, it is necessary to justify whether the flow regime can be described by the standard continuum approach or the statistical mechanics. Concerning the dimensionless analysis, the Knudsen number is normally used to characterize the flow regime and is defined by

C L Kn   , ( 5
    F τ I v v v              p t , (5.3)
where F is the volume force and τ is the viscous stress tensor written by

     I v v v τ          3 2 T .
(5.4)

In order to take into account the conservation of mass, the mass continuity equation is included:

  0       v   t . ( 5.5)
Moreover, the conservation for the i th species in a gas mixture is given by the species balance equation:

  i i i i R t             v j . ( 5.6)
Here, i R denotes the reaction rate of the i th species. The total diffusive mass flux i j is composed of two terms. The first term F i j is the flux due to a gradient of concentration described by Fick's law. The second term, the thermodiffusion flux T i j is the flux due to the temperature gradient described by the Soret's effect. Thus the total diffusive flux becomes

T i F i i j j j   . (5.7)
Concerning the interaction of species in the multicomponent diffusion, the governing equations for the flow behavior are better described by the Maxwell-Stefan equations written as 70

T T D D T i j N j ij i i       d j 1  , (5.8)
where ik D ~ is the multicomponent diffusion coefficient which can be computed using the binary diffusion coefficient by the matrix transformation written as follows:

          i k kj i ik i k kj i j i ij j i adjB D adjB D x x   (5.
               Q l l l j j j j j j j P x P x 1 1 g g d , (5.10)
where j g denotes the external force acting on species j. The coefficient T i D refers to the thermal diffusion coefficient. In this thesis, we use the Maxwell-Stefan model as the governing equations for the flow behavior since it is more comprehensive than the Fick diffusion model. The former takes into account the interaction between the gas species while the latter one does only for the interaction between gas species and the bulk material (solvent). In fact, Fick's model is the limit case of the Maxwell-Stefan model when the gas mixture is highly diluted. In the standard PVT set up, the crucible is not completely closed. First, graphite itself is not a tight material because of its porosity. Second, the joint between different graphite parts are possible leak that allows exchanges between the inner part of the crucible and the reactor's chamber. This argument is verified since the growth rate is a function of the ambient pressure. Moreover by introducing the nitrogen gas flow outside the growth cavity one can obtain even heavily nitrogen doped SiC crystal. Thus we assume that the internal pressure (inside the growth cavity) equals to the external pressure P in (5.10) due to the exchange of the gas. Thus the total pressure of the gas mixture in the growth cavity is constant. Moreover, there is no external force acting on the vapor species. The diffusional driving force (5.10) is now only a function of mole fraction gradient j x  . If the thermal diffusion coefficient is very low, the Soret's effect can be also neglected. This is the case for the growth of SiC by sublimation method. We now rewrite the mass flux (5.8) as

j N j ij i i x D     1  j . (5.11)
It is more convenient to consider the molar flux instead of the mass flux. By using the relation

M M x i i i   , (5.12)
where M is the mean molar mass defined by

   N i i i M x M 1 , (5.13)
the diffusive molar flux can be written

j N j ij i i i x D M x M     1  j .
(5.14)

Boundary conditions

The molar flux at the reacting boundaries N i is described by the summation of two contributions:

v j N i i M ρx M i i   . (5.15)
The first term is the diffusive molar flux defined in (5.14). The second term is the convective molar flux. At the solid-vapor interface, the large volume change caused by the sublimation or condensation reactions at the interfaces produces the advective flow or Stefan's flow velocity:

n v              1 1 N i i i N M PM RT , (5.16)
where M is the mean molar mass. Here, the summation includes only the reactive species since the inert gas does not participate in the heterogeneous reactions. It can be seen from both (5.15) and (5.16) that there is a coupling between the fluid flow and the mass transfer since the velocity of the fluid at the reacting boundary has to be computed from the net species fluxes. At the nonreacting boundaries, such as at the graphite walls where the deposition does not occur, all the velocity components are set to zero. It is worth noting that the contribution of the buoyancy force in the fluid dynamics of the vapor phase is much smaller than the effect of Stefan flow in the PVT process under the practical growth conditions. This follows the fact that such growth process is typically carried out in low pressure. The buoyancy force is then very small. The macroscopic flow of the reactive vapor species owing to the sublimation/condensation processes determines the flow pattern since the velocity magnitude in the growth cell taking into account the Stefan flow (10 -3 m/s) can be three order of magnitude higher than the velocity field caused by purely buoyancy flow (10 -6 m/s) [2]. Thus by taking the Stefan flow into account, the buoyancy force can be neglected.

We assume that there is no homogeneous vapor phase reaction. Only the heterogeneous reactions at the solid-vapor interface are considered. The boundary molar flux was formulated by the Hertz-Knudsen model

  RT M P P N N N i i i i i des i ads i i    2 1      , (5.17)
where N i is the net molar flux of the i th vapor species participating in the heterogeneous reactions, P i is the partial pressure, P i * is the equilibrium partial pressure, M i is the molar mass, R is the gas constant, T is the temperature, γ i and α i are the sticking and evaporation coefficients, respectively. The transport of atom or molecule from the vapor phase to the solid phase is associated with several processes, including, the adsorption, diffusion, desorption, and incorporation [3]. Here, the sticking coefficient refers to the probability that an atom in the gas phase adsorbs, diffuses, and finally incorporates to the bulk solid by attachment to the kink site at a step edge. The reverse process is assumed for the definition of the evaporation coefficient where the atom finally desorbs from the surface. The schematic representation of the Hertz-Knudsen model is shown in Fig. 5.2. A thin Knudsen layer exists between the solid SiC crystal and the vapor as mentioned in [4] where its thickness is approximately one of the mean free path of the molecules [5]. Within such layer, two different molecular streams of adsorbed ( ads

  n v                     RT Px x D RT Px P P N i j N j ij i i i i i i 1 ~, (5.18)
where n is the normal vector to the surface, and

RT M i i   2 1  . (5.19)
The condensation/sublimation rate is proportional to the difference between the actual partial pressure and the equilibrium partial pressure. If the species partial pressure is higher than its equilibrium pressure at the temperature of the crystal surface, the vapor state is supersaturated and the net molar flux becomes positive (growth). On the opposite, the vapor is in the undersaturated state if the species partial pressure is lower than its equilibrium pressure, and the flux becomes negative (sublimation). Eq. (5.18) can also be reformulated in terms of the unknown variable for solving the species transport problem (mole fraction) using the relation

x i =P i /P:   n v                     RT Px x D RT Px P P x N i j N j ij i i i i i i 1 ~. (5.20)
The sublimation growth of SiC is usually carried out in a reduced pressure of a neutral gas (Ar) in a semi-closed graphite crucible. In such growth environment, the SiC source powder decomposes into a Si-rich gas phase, graphite, and the remaining SiC. Due to this non-congruent nature in sublimation of SiC, the three-phase heterogeneous equilibrium "SiC-C-vapor" can be assumed at the source-vapor interface. The equilibrium partial pressure of the vapor species can be determined uniquely as a function of temperature. The three main Si-C vapor species are Si, Si 2 C, and SiC 2 [15][16][17] as identified and quantified by high temperature mass spectrometry using Knudsen effusion cell. With such a Si-rich gas phase, the crystallization process leads to stoichiometric SiC crystals. It is worth nothing that SiC can accept a slight deviation from stoichiometry, mainly in terms of C deficiency [18], but this deviation from stoichiometry is rather small and will not be considered.

From the process modeling point of view, the implementation of the stoichiometry condition is not unique as different approaches can be found. Gao et al. [9] do not consider any stoichiometry constraint at growth interface. In the case of the close space sublimation method [19], the stoichiometric sublimation of the source is imposed. Furthermore, as the sublimation rate and the condensation (growth) rate are assumed to be identical, meaning that no accumulation of either Si or C is allowed, the condensation is thus stoichiometric. By using the two-phase heterogeneous equilibrium "SiC-vapor" [20] for the silicon carbide surface, there exists one thermodynamic degree of freedom, thus the vapor phase composition cannot be determined uniquely. The equation of stoichiometric incorporation was then suggested to close the set of equations. However, the forced convection induced by the boundary Stefan flow and the multicomponent diffusion effects were neglected in this model by assuming a dilute system in which the concentration of carrier gas prevails over the reactive species. The modification of this model to include these two effects was presented [4]. In such case, however, they consider one-dimensional mass transport for the axial mass fluxes while the radial ones are determined from the concrete model of mass exchange between the growth chamber and the external environment. Ma et al. [10] assumed that reactive vapor species in the growth chamber mix at perfect stoichiometric condition and the growth rate can be determined from the one-dimensional mass transport-kinetics.

In the PVT growth of SiC, the mole fraction of silicon is higher than the one of carbon in the vapor phase. The Si/C ratio close to the crystal surface should be larger than unity for a wide range of growth temperatures. From (5.18), assuming that the sticking coefficient of all species is unity, one would expect a higher global silicon flux than the carbon flux, resulting in a silicon excess in the condensed phase, i.e. in the crystal. By introducing the evaporation coefficients extracted from high temperature mass spectrometry studies [16,21] in the modeling of SiC bulk growth [22], the total pressure of the reactive vapor species over solid SiC were found to decrease. In fact, the imbalanced silicon and carbon fluxes at the growing crystal should still be occurring because the stoichiometric condition was not imposed. It is very important to take into account the stoichiometric incorporation for the process modeling of SiC single crystal growth in order to realize the actual growth mechanism and understand the process. Two main approaches used to tackle such problem will be presented and discussed: the comparative approach and the coupled approach.

Comparative approach

In this approach, we will consider the growth cavity as a whole system consisting of vapor phase (Si, Si 2 C, SiC 2 , Ar) and the condensed phase including SiC due to the presence of SiC seed and SiC powder source, and C due to the presence of graphite crucible and the graphitized powder source. This allows us to prior assume the three-phase heterogeneous equilibrium "SiC-C-vapor". The temperature dependent vapor phase composition is obtained from [23]. The evaporation coefficients extracted from high temperature mass spectrometry during the sublimation of "SiC-C" diphasic mixture [16,21] shall be used at all the solid-vapor interfaces. We will present the three models to take into account the stoichiometric incorporation of Si and C atoms into SiC crystal grown from PVT method. These models are built based on the two methods: the atomic fluxes balance, and the adaptable sticking coefficient(s). The terms comparative approach means that both numerical and thermodynamic modeling are performed and compared. The appropriate boundary conditions for stoichiometric condition are discussed in terms of both kinetics and thermodynamic mechanisms.

Atomic fluxes balance (AFB)

To strain the stoichiometric incorporation of Si and C atoms at the crystal-vapor interface, one possibility is to balance the atomic silicon and carbon fluxes:

2 2 2 2 2 2 SiC C Si C SiC C Si Si Si N N F N N N F         . (5.21)
This condition can be formulated by using the weak constraint. With this method, an extra unknown called the Lagrange multiplier has been added to the system of unknowns to ensure the constraint residual equation. From equation (5.21), the weak constraint residual equation is of the forms:

0 2 2       SiC C Si Si C Si N N N F F  .
(5.22)

Adaptable sticking coefficients

We are developing the second method by considering the decomposition of the molar flux at the boundary. From Eq. (5.21), the atomic fluxes balancing equation can be rewritten:

2 2 2 2 2 2 2 2 2 2 SiC SiC C Si C Si C SiC SiC C Si C Si Si Si Si F F                  
, where

(5.23)

  i i i i P P x      . (5.24)
The first assumption for this method is that the sticking coefficient represents the maximum kinetically limited deposition rate. We propose here two models to take the stoichiometry problem into account. The first method is to modify the sticking coefficient of Si, and Si 2 C

species but keeping the one of SiC 2 unchanged. From now on we will call this model "Mod2".

This model is based on the consideration that the species flux providing more C than Si atoms should be kept at its highest deposition rate. In other words, SiC 2 is considered as a rate determining species. Thus γ SiC2 is the unmodified value. Then, Si and Si 2 C species fluxes which provide mainly Si should be reduced to match the carbon flux. This can be done by modifying their sticking coefficients to have the lower values, meaning that these two species deposit into the growing crystal at the rate lower than their maximum kinetically limited rate in order to obey the stoichiometry condition. Rearrangement of equation (5.23) using the modified sticking

coefficients * Si  and * 2 C Si  gives 2 2 2 2 * * * SiC Si C Si Si SiC SiC Si          .
(5.25)

From Eq. (5.25), one unknown modified sticking coefficient is remained undefined. In order to close the equation system, another assumption is necessary. Here we assume that the ratio of the initial experimentally observed values and the modified ones are equal:

Si C Si Si C Si     2 2 * *  . (5.26)
The modified sticking coefficient of Si species can be determined from (5.25) and the one of Si 2 C from (5.26). These modified sticking coefficients are dependent on boundary molar fluxes of the reactive species.

The second model is to modify the sticking coefficient of Si species but keep the other two constant. We will call this model "Mod1". The assumption for this method is to keep the carbon supply species at their maximum kinetically limited rates. In this case, assumption (5.26) is not necessary and the modified sticking coefficient of Si species is expressed by

Si C Si C Si SiC SiC Si      2 2 2 2 *    .
(5.27)

Numerical modeling

The growth system is similar to the one shown in Fig. 4.12. The temperature distribution in the growth chamber was calculated as described in the previous chapter. The temperature on the symmetry axis is 2300 °C at the seed and 2340 °C at the source. The vapor species sublimed from the hot source powder will become supersaturated and deposit on the colder seed side.

When computing the mass transport, all methods discussed in the previous section provide stoichiometric incorporation of Si and C atoms into SiC crystal. Without the stoichiometry condition, the Si/C ratio of the solid would be around 2.15 which is much larger than the limit of experimental values [18]. Thus the stoichiometry condition must be taken into account. The distributions of the species mole fraction along the symmetry axis close to the seed are plotted in The summation of the species mole fraction must be unity. Thus small deviations of the SiC 2 mole fraction in Mod2 and the deviations of Si 2 C and SiC 2 in Mod1 from the initial system without constraint are observed (Fig. 5.3) even if their sticking coefficients are kept unchanged in the method considered. The Si/C ratio distribution along the symmetry axis is plotted in Fig. 5.4. The distributions of Si/C ratio look very similar when using the atomic fluxes balance compared to the system without constraint. However, their gradients change in the area very close to the seed surface as a result of the changes of mole fraction gradient of the vapor species as shown in Fig. 5.3. With the adaptable sticking coefficients method, the vapor phase becomes richer in Si, not only along the symmetry axis but also everywhere in the growth cavity because less Si flux contributes to the condensed phase.

The radial distribution of the total species molar fluxes obtained from all models are plotted in Fig. 5.5. If we assume a simple heterogeneous reaction of the forms xSi(gas)+ySi 2 C(gas)+zSiC 2 (gas) ↔ mSiC (solid), (5.28) where x, y, z, and m are the stoichiometric coefficients at the crystal-vapor interface, the stoichiometric coefficients can be estimated from Fig. In this case, it does not matter if T 2 or T 1 is used to plot since ΔT is constant. When T 2 is fixed at 2200 °C, the order of each species contributions is similar to the first case at low ΔT (Fig. 5.8b). However, at high ΔT in which τ Si2C is larger than τ SiC2 , τ Si becomes negative meaning that the Si vapor is produced from SiC instead of consuming. Further increase of ΔT, i.e. for ΔT=300 °C, leads to the formation of Si in the liquid phase. However, such extreme conditions where the net Si flux leads to the emission of Si vapor species or the precipitation of extra condense phase are out of the scope of this work and will not be considered in this model. By considering the heterogeneous reaction (5.30), the species partial contributions

  2 2 2 SiC C Si i i      
are plotted in Fig. 5.9 as a function of seed temperature divided by temperature difference (T 2 /ΔT). At high T 2 /ΔT ratio, the species partial contribution of Si and Si 2 C are quite similar. As T 2 /ΔT decreases, the partial contribution of Si and SiC 2 decreases while the one of Si 2 C increases. The partial contribution of Si 2 C is approaching SiC 2 and the one of Si is approaching zero. The general trend is that the partial contribution of SiC 2 is always much larger than Si at this temperature range.

Discussion

The atomic fluxes balance method improves the initial system's boundary conditions.

This method does not change the Stefan flow, but adapts the diffusive flux just near the seed to obtain stoichiometry. In contrast, both mole fraction and mole fraction gradient are affected when using the adaptable sticking coefficients method. Thus, the question arises for which method corresponds more to the physics of the crystal growth?

The chemical description of the non-congruent condensation from a complex vapor phase

is still an open issue. The recent paper from Fujimoto et al. [25] is an interesting approach as it allows a comprehensive understanding of vapor phase composition variations in P-T-x phase diagram. It was pointed out that, the SiC source powder thermally decomposes into Si-C vapor with a composition of v 1 and solid carbon. As the temperature is lower at the seed side, vapor v 1 solidifies to form SiC single crystal and the vapor phase composition becomes v 2 which is richer in Si compared to v 1 . From Fig. 5.6, both methods show that the Si/C ratio is higher on the seed side than the source side. Thus both methods seem to be acceptable regarding the analysis of phase diagram discussed in such report.

In this section, it is worth noting also that the similar thermodynamic database [23] was used in both numerical modeling and thermodynamic calculations to ensure the compatibility when comparing those two sets of results. The new question arises concerning the effect of the choice of the thermodynamic database on the full process modeling. This will be discussed in the following section.

Effect of choice of thermodynamic databases

While developing a numerical model for the global process simulation, a careful analysis and choice of the databases is crucial. To date different thermodynamic databases are available for describing the vapor phase in equilibrium in the Si-C chemical system. They are differing by the methods used for their measurements. The first attempt for the gaseous phase study of the Si-C system was done by Drowart et al. [15] using mass spectrometry. The composition of the gas phase and enthalpies of formation for the main species were determined. Some other works related to the determination of equilibrium partial pressure of vapor species in Si-C system were proposed. Thermodynamic analysis of equilibrium processes in the gas phase during silicon carbide sublimation done by Lilov [26] has shown that the composition and the stoichiometry of the gas phase were strongly dependent on temperature. High-temperature mass spectrometry with a multiple Knudsen-cell furnace was used in [16,17]. The enthalpies of the formation for Si, Si 2 , Si 3 , Si 2 C, and SiC 2 molecules were determined. The SGPS database for Si-C system was compiled by SGTE group [27] using data from [START_REF] Schulze | Controlled Growth of 15R-SiC Single Crystals by the Modified Lely Method[END_REF]. This database can be found in commercial software, such as FactSage [START_REF] Tsavdaris | Macroscopic Approach to the Nucleation and Propagation of Foreign Polytype Inclusions during Seeded Sublimation Growth of Silicon Carbide[END_REF]. Another database origin is the JANAF table [23]. The aim of this section is to assess the different thermodynamic databases available and compare them regarding the global modeling of SiC growth by PVT. The calculated results namely, growth rates and Si/C ratios in front of the seed crystal are compared and discussed.

Equilibrium partial pressure in front of SiC-C diphasic mixture

The temperature dependent equilibrium partial pressure of the three main vapor species in "SiC-C-vapor" exists entirely in the growth chamber reduces the thermodynamic degree of freedom of the heterogeneous process. This prevents the variation of the state variable at the growth interface. Thus the kinetics parameters, such as the sticking coefficients, have to be modified in order to realize the stoichiometric growth of SiC as we have done so far. In fact, the problem of sticking coefficient still remains as an open issue. They are very complex and should include energetic and kinetic terms, i. e., temperature, velocity, surface roughness [6,32,33].

The experimental values will be necessary to improve the accuracy of the model. However, an alternative approach to avoid dealing with the evaporation and sticking coefficient as well as cross-checking the thermodynamic databases will be presented in the following section.

Coupled approach

Apart from the problem of boundary conditions for considering the stoichiometry of SiC in the growth process modeling, there are also other strong issues which are not considered in a satisfactory way. This is the case for instance for the thermochemistry of the Si-C system and its coupling to mass transfer. Recently, it has been reconsidered and a better understanding of the solid-vapor interface chemistry was reported [25,34]. But two essential questions are still open: i) how to treat the thermochemistry of the solid SiC crystal? ii) How the latter could be coupled to the numerical modeling of the process? In this section, we present a method to take into account the stoichiometric incorporation of vapor species. More importantly, the proposed method gives an assessment to the chemistry of solid SiC crystal, i.e. the activities of both Si and C atoms in the crystal during the PVT growth process. This approach employs the fully coupling between numerical modeling and the thermodynamic calculations.

Method and calculations

In contrast to the comparative method, this coupled method does not prior assume the three-phase heterogeneous equilibrium (SiC-C-vapor with variance = 0). Thus at the crystalvapor interface, the equilibrium partial pressure cannot be computed directly from the temperature of the solid phase (SiC-gas with variance = 1). The extra thermodynamic degree of freedom is compensated by the condition for SiC stoichiometry via the consideration of SiC solid solution (s, s). By treating SiC as a solid solution, following the work done in [34] the Si and C components can be considered in the non-stoichiometric Si x C y compound according to the equilibrium Si x C y = xSi (in SiC (s, s) ) + yC (in SiC (s, s) ).

(5.31)

Since the SiC crystal is nearly stoichiometric [18], x and y are close to unity and the equilibrium constant for reaction (5.31) is written by

  SiC C Si p a a a T K   1 .
(5.32)

This equilibrium constant is calculated from the standard free enthalpies obtained from the JANAF thermochemical database [23] by the following expression

    T K RT T G p r ln     (5.33)
The activity of the pure components equal to unity by the definition of the standard state considering that there is only small deviation from the stoichiometric composition. The

2 2 2 2 2 2 2 2 2 2 2 2 * * * C Si SiC C Si C Si Si Si Si C Si C Si SiC C Si C Si C Si Si Si Si C Si C Si P P P P P P                  . ( 5.40)
In this model, we assume that the condensation coefficients of all vapor species are SiC  representing the mean condensation mechanism. The same assumption is applied for the evaporation coefficients so the mean evaporation coefficient is SiC  . Using the assumption that SiC SiC

  

, the solution of (5.40) is independent of the condensation and evaporation coefficients. Thus from the relations (5.32), and (5.37)-(5.39), Eq. (5.40) is rewritten

A a K K a K K a K Si C Si Si C Si Si Si    1 2 2 4 2 1 3 1 2    , ( 5.41)
where the differences in Si and C fluxes condensing on the seed surface is denoted by

2 2 2 2 C Si C Si C Si C Si Si Si P P P A       .
(5.42)

Eq. (5.41) can be rearranged in forms of the quadratic equation:

  0 2 2 4 2 1 2 3 1 2      C Si Si Si C Si Si K K Aa a K K K    .
(5.43)

The analytic solutions of Eq. (5.43) are

     C Si Si C Si C Si Si Si K K K K K K K K A A a 2 2 2 3 1 2 4 2 1 3 1 2 2 2 4              . ( 5.44)
Since all parameters are real and the discriminant

    2 2 4 2 1 3 1 2 2 4 C Si C Si Si K K K K K A        is
always positive, the two roots are real and unequal. It can be further proved that one root is positive and another one is negative. The negative value of activity has no physical meaning.

Thus there exists only one positive value of Si activity that verifies the steady-state of growth.

The corresponding C activity can be computed using Eq. (5.32). One may note that by adding the Si 2 and/or Si 3 species in the vapor phase, the quadratic equation (5.43) becomes either cubic or quartic equations which should still be possible to solve analytically or numerically. However, these additional species should not affect the calculated Si activity since the vapor pressure of Si is around two orders of magnitude higher than Si 2 and Si 3 even in the extreme case at SiC-Si phase limit [34]. Thus in general growth conditions, the error from excluding these two species should results in less than 1% error of the calculated activity.

Recall the Hertz-Knudsen formulation for the condensation-sublimation at the seed surface (5.20), we know from (5.37)-(5.39), (5.42) and (5.44) that the equilibrium partial pressure is a function of the partial pressure of all vapor species. Thus Eq. (5.20) may be written as proposed for the growth under equilibrium condition [35]. The results showed that the concentration of different type of both Schottky and Frenkel defects are dependent on the partial pressure of silicon vapor. It may be possible to apply the results from our process modeling concerning the equilibrium partial pressure of Si to such analysis in order to study the behavior of the defects under the real growth conditions.

    n n                             1 1 3 2 1 1 3 2 1 ) , , ( ) , , (
N i i i i i i j N j ij i i i i x x x f P x M M x x D RT Px x x x f P x     . ( 5
The capability of this model in predicting the state (Si-or C-rich) of the growing crystal should be able to correlate to the occurrence of cubic or hexagonal polytype with the growth conditions. A quantitative LEED study suggested that a cubic stacking sequence is induced by the surface reconstruction with the excess Si [36]. This was found only when the surface of 4H-SiC (0001) is Si enriched during the annealing process. It is still unclear if this mechanism is applicable for the growth process where much higher temperature is used. Normally at such high temperature, the 3C polytype is thermally instable and transform to the hexagonal polytype due the solid state transformation caused by the imperfections of the crystal [37]. Birnie et al. [18] found the predominant 3C-SiC grains when equilibrate the 3C-SiC grains in Si-rich environment while 6H-SiC grains are predominant in the C-rich environment. The interesting point is that such experiments were done at very high temperature (2400 °C). This might suggest that the mechanism proposed by [36] overcomes the solid state transformation. The growth of cubic SiC prefers the low growth temperature typically below 2000 °C to prevent the polytype transformation. Using such low temperature, growth of 3C-SiC on 3C-SiC can be stabilized using the separated Si container [38]. This extra Si source increases the Si partial pressure by 25-30%. Keeping low temperature at the seed, very large ΔT (500 °C) and very low pressure (10 -3 -10 -4 mbar) were used in [39] and the 3C-SiC can be grown on both 6H(0001) and 3C(100)

substrates. All the experimental results discussed here suggest that the cubic polytype of SiC is preferable for the growth under Si rich condition. Thus, we may propose a guideline for the stabilization of SiC polytype according the chemistry at the growing surface predicted from our model. First, the growth environment providing the crystal chemistry closer to the SiC-Si phase limit, i.e. high ΔT and/or low pressure, should result in higher probability for the nucleation of cubic polytype. The crystal imperfections such as stress and defects should be also minimized in order to avoid the polytypic transformation. Moreover, the growth under low ΔT and/or high pressure bringing the crystal chemistry closer to SiC-C phase limit should be preferable for stabilizing the hexagonal polytype.

Another model was proposed for the transformation of SiC polytypes during the growth of an epitaxial layer [40]. This may serve as a basis to control the amount of C V N to stabilize the preferential SiC polytype in the growth experiment.

To fabricate n-and p-type SiC epilayer, nitrogen and aluminum are generally used as a dopant, respectively. The spectroscopic evidences from Choyke [43] and Davis et al. [44] shown that the nitrogen atom occupies the C site in SiC crystal while the aluminum atom substitute the Si site [45]. Larkin et al. [46] proposed that the Si/C ratio in the CVD growth reactor strongly affect the doping incorporation for epilayer grown on 6H-SiC Si-face. Nitrogen and carbon compete for the C sites while aluminum and Si compete for the Si sites of the growing SiC epilayer. This is the well-known 'site-competition' model. Even though, it is still not clear why this model is failed to describe the behavior on the C-face. The incorporation mechanism of nitrogen, aluminum, and boron was also proposed by Kimoto et al. [45]. Such mechanism involves the structure of the growing surface under Si-or C-rich condition. The surface coverage is differed depending on whether it is Si-or C-face. The C coverage increases with increasing C/Si ratio on Si face thus preventing the incorporation of N atoms. On the other hand, the surface on C-face is believed to be terminated with C atoms independent on the C/Si ratio in the range of 1-2. Thus the site-competition model does not effectively work for the growth on C-face.

Concerning the CVD process, the amount and ratio of injected precursor gas are controllable.

However this is not the case for PVT since the compositions are always enriched with Si by considering only the vapor phase itself. Indeed the surface chemistry of the growing crystal was completely unknown. Thus in order to quantify the doping incorporation, the information predicted from our model could be correlated with the surface coverage. For example, the growth conditions bringing the surface chemistry far from SiC-C phase boundary should multiply the competitiveness of the nitrogen providing higher nitrogen incorporation rate. Nevertheless, this requires the extension of the model to include N into the system.

Three examples concerning the point defects formation, occurrence of cubic or hexagonal polytype, and the quantitative doping incorporation showed the capability of the newly proposed modeling approach. The knowledge of the chemistry of the growing SiC crystal provides one step further for understanding and developing the growth process. Moreover, it is possible to extend the proposed coupled approach to consider the interaction between the vapor species and the graphite container.

Interaction between vapor species and graphite crucible

In seeded sublimation growth process (PVT) of SiC single crystals, graphite container is generally used. Despite the small number of gaseous species commonly considered (Si, Si 2 C, SiC 2 ), the heterogeneous reactions between vapors and solids are very complex. The shape of the crystal growth front does not only depend on the isothermal shape but also on the global mass transport and chemical reactions considering for example, the heterogeneous adsorption/desorption processes and the interaction with the graphite crucible or tantalum crucible [47]. To optimize the growth process for obtaining high quality crystal, the shape of the growing crystal becomes an important issue. For example, the changes in the curvature of the crystal shape can induce polytype destabilization [48].

Several reports mentioned the important role of the graphite crucible as a carbon provider for the crystal growth [2,[49][50][51][52]. The strongest supporting report for this issue was the work done by Herro et al. [49]. With the starting amount of 88% of 13 C in the source powder, the crystal consisted of only 40% of 13 C after growth. Numerical modeling results concerning the interaction between the vapor species and the graphite wall were presented [2]. However, no detail of the calculation was given. In this section, we are considering the chemical reactions of the vapor species with the graphite wall. We will focus on the effects of additional mass transport to the growing crystal and the resulting crystal growth front shape.

Method and calculations

In the previous section, we consider the equilibrium condition between the solid solution of SiC and the vapor phase for the growing SiC crystal. In this case we may apply the same procedure but replace SiC solid solution crystal by C solid representing the graphite wall.

Concerning the etching process, elementary processes should include the element of the material to etch. In this study we assume that the excess Si vapor does not incorporate into the graphite wall but interact with it and produce SiC 2 species. This is supported by the fact that what is experimentally observed on most of the crucible walls is graphite corrosion and not SiC In this case the standard state defines a C =1, the equilibrium partial pressure of SiC 2 vapor species in front of the graphite surface can be written as a function of the equilibrium partial pressure of Si and its equilibrium constant. In order to close a set of equations, a mass balance equation for stoichiometric SiC (5.21) will be replaced by another mass balance equation considering that the silicon atoms will not be incorporated into the graphite wall

0 2    SiC Si Si N N F . ( 5
       , (5.50)
The equilibrium partial pressure of Si can be obtained:

2 5 2 2 * C Si Si C Si C Si Si Si Si K P P P        .
(5.51)

Since the equilibrium partial pressure of SiC 2 species can be computed from the equilibrium partial pressure of Si, these equilibrium partial pressures can be subtracted back into the boundary condition at the graphite wall which is in fact similar to (5.20). Then a closed set of equations as shown in (5.45) can be formulated at the C-vapor interface. The important point to consider in this case is the mean sticking coefficient C  at the graphite surface. We should now define this value as an etching coefficient. Even if the solution of (5.51) is independent on the etching coefficient, the value of this coefficient will determine the reactivity of the graphite container with the vapor species. The crucible with good quality should be inert to the reactive gases and the etching coefficient becomes zero. In this section, we consider a small deviation from the perfect inert graphite wall. Thus we set a very low value ( C  =0.0005) of the etching coefficient in order to investigate the effect of the etching of the crucible in the process modeling.

Results and discussion

The schematic representation of the growth chamber used in the experiment and modeling is shown in Fig. 5.17a. The two main graphite parts consist of, the graphite guide (part A), and the seed holder (part B). For the modeling without graphite etching, only the deposition of polycrystalline SiC is allowed on the graphite walls. On the other hand, when using the model integrated with the graphite etching, both the deposition of polycrystalline and the consumption in the area close to the surface. The latter case was confirmed using grazing incidence X-ray diffraction. SEM observations were conducted in order to investigate the surface effect of the heat treatment. However, we could not observe any remarkable difference in surface morphology of the graphite parts in both heat treated and non-heat treated ones. We can thus conclude that all the evolutions observed in our case are only coming from the interaction between the graphite wall surfaces and the gaseous species and not from the evolution of the graphite itself upon annealing. Our graphite is perfectly stable with time upon heat treatment at high temperature.

General modeling results and validations

The typical temperature distribution in the growth chamber and the fluid velocity are shown in Fig. 5 species is directed from the source to the seed. On the other hand, the partial pressure of argon gas shows the reverse accumulation behavior. The diffusive flux of argon is directed from the seed to the source which opposes its convective flux. This is corresponding to the boundary condition of the total zero deposited argon flux at all boundaries. The argon gas following the convective flow of the sublimed species will accumulate in front of the seed. Thus the concentration of argon in front of the seed is higher than the powder and a reversed concentration gradient is created. The highest sublimation rate is located at the lateral part and the lowest is located at the center. The back flows of the Si-rich species are also observed at the center of the source powder. This may result in the recrystallization in the powder source at its central part which may block the path ways of the sublimated species in the bulk powder below the surface.

The high sublimation rate in the lateral part of the powder can lead to the fast graphitization rate since the sublimation of SiC is non-congruent. vapor interface results in the non-stoichiometric growth. The stoichiometric incorporation of the vapor species has been taken into account by both comparative and coupled approaches. The comparative approach requires a prior assumption that the three phase heterogeneous equilibrium exists entirely in the growth chamber. Such approach has to deal with the kinetic parameters (the sticking coefficients) since the thermodynamic degree of freedom is reduced. On the other hand, the coupled approach, treating SiC as a solid solution, requires a strong coupling with thermodynamics of the crystal growth. This model has been shown to be more comprehensive for the growth process for two main reasons. First good agreement between modeling and experiments were achieved concerning the growth/etching areas, growth rate, and interface shape of the crystal. Second, it really opens the way to have an access to the chemistry of SiC crystal.

The physico-chemical parameters of the solid phase, activities of Si and C in the SiC compound, can be linked to controllable parameters: temperature and pressure. Further correlations with the observable parameters such as, point defects density, stable polytypes, and doping concentration are possible. However, this requires an extension of the model.

Chapter 6 Top Seeded Solution Growth

So far, we have been focused on the process modeling of SiC growth from the vapor phase. In this chapter, growth from liquid phase will be considered as we will discuss the modeling of the Top Seeded Solution Growth (TSSG) process. Physical model and boundary conditions will be discussed first. Then the modeling results concerning the effects of the process parameters including crystal rotation speed and operating frequency will be given. Finally, we will present the combined numerical and analytical modeling to describe the interaction between fluid flow and the step flow at the growth front and its validation.

Physical model and boundary conditions for induction heating and fluid dynamics

The typical configuration of the TSSG reactor is shown in Fig. 6.1. Pure liquid silicon used as a solvent is placed in the high density graphite crucible. Similar to the PVT process described in the 4 th chapter, the induction heating is employed as a method to heat the graphite crucible. Almost all the sets of equation and the numerical procedure are the same. However, there are two main differences in the TSSG model. First, we assume the incompressible flow of fluid in order to save the computational cost. The fluid density is now constant so the continuity equation (5.5) becomes

0    v  . (6.1)
The divergence of the fluid velocity vanishes, so the momentum conservation equation (5.3) for the steady state calculation becomes

          F v v I F I v v                     T p p . ( 6
 g g F F 0 0 buo T T         . (6.3)
The Boussinesq approximation is accurate when   1 0   T T  . This is the case for the TSSG process for the growth of SiC. If the crucible is heated from the side wall, the fluid will flow up along the side wall and the convection pattern caused by the buoyancy force will be similar to Fig. 6.2a. On the other hand, the fluid will flow up at the center of liquid if the crucible is heated from the bottom and the flow directions will be reversed compared to Fig. 6.2a. The intensity of the buoyancy convection can be estimated from one of the dimensionless numbers, the Rayleigh number computed from the product of the Grashof number and the Prandtl number. The Grashof number is the ratio of the buoyancy to the viscous forces:

2 3 2 2 3     TL g TL g Gr     . (6.4)
The Prandtl number is the ratio of momentum diffusivity to the thermal diffusivity defined as

k C p      Pr . (6.5)
Thus, the Rayleigh number can be written

k C TL g TL g Gr Ra p     3 2 3 Pr       . (6.6)
Two critical values (Ra c1 , Ra c2 ) of the Rayleigh number can be used to characterize the stability of the flow. The transition from no flow to steady flow occurs at Ra c1 =1708 [5]. The transition from steady flow to the unsteady and time-dependent flow occur at Ra c2 =40000 [6]. From (6.6), it can be seen that the buoyancy convection can be controlled by the liquid height and the axial temperature gradient.

Forced convection

The rotation of either seed or crucible directly induces the fluid movement in the azimuthal direction. The force convection caused by the crystal rotation is implemented via the boundary condition for the azimuthal fluid velocity at the seed-liquid interface:

r seed    v . (6.7)
The steady state liquid motion in a rotating liquid is two-dimensional with respect to the coordinate axis rotating with liquid [7]. The appearance of such two-dimensional flow to an observer staying in a stationary frame with respect to the rotating fluid is a set of spiral streamlines rising or falling with a common axis of a symmetry axis of the rotation. In a twodimensional axisymmetric configuration used in the simulations, the flow in the azimuthal direction (out-of-plane) is uncoupled with the radial and axial directions (in-plane) so such spiral movement cannot be observed. However, the flow in the rz-plane induced by the azimuthal flow is still correctly described. To exemplify, the crystal acts as a centrifugal fan drawing in the fluid axially, transmit angular momentum to it in the boundary layer and then expelling it tangentially [8]. Such mechanism causes the flow pattern shown in Fig. 6.2b where the back flow is caused by the presence of the crucible wall. In this case, the Reynolds number defined as a ratio of inertial force to viscous,

 L stir v Re  , ( 6.8)
can be used to describe the stability of the flow. For a flow in a pipe of diameter d, the laminar flow occurs when Re d < 2000 while the turbulent flow occurs when Re d > 4000 [9]. In between those two critical values, a transition flow where both laminar and turbulent flows are possible depending on other factors such as the smoothness of the flow. The forced convection can be controlled by adjusting the rotation speed and the radius of the seed (or crucible).

Marangoni convection

The Marangoni convection appears when the surface tension of the interface, in this case the liquid-gas interface, is dependent on the temperature such as

    0 0 T T T       . (6.9)
where χ is the Marangoni coefficient describing temperature derivative of the surface tension.

The governing equation describing the Marangoni effect at the interface is written by

s T s            s n , (6.10)
where s is the tangential vector along the interface. Such equation specifies that the shear stress on the surface is proportional to the temperature gradient on such surface. The shear force in the plane tangential to the interface can be written as .11) If the free surface does not deform, v n =0. From (6.10) and (6.11), we obtain the relation

                                                                   n s n n s n s s s n s s n s n n v v 1 0 v 2 v v v v v 2 0 1 s n . ( 6
n v s        s T . ( 6.12)
We are taking into account the Marangoni force via a set of boundary conditions expressing the balancing of the temperature gradient on the surface with the shear stress (6.12). Such boundary condition is implemented by a weak form finite element formulation. The weak contribution The Marangoni convection is quite difficult to control. Even though, a special type of crucible configuration was implemented in the numerical modeling study with aim to minimize such effect on the SiC crystal growth [1].

0 v s                             d n d s T ( 6

Electromagnetic convection

Electromagnetic stirring is normally the side effect of the heating used in experiments.

The movement of charges in the magnetic field induces the Lorentz force. This is the case for our growth system in which the magnetic field and the current flow in liquid are presented due to the induction heating. In such case, the volume force computed from the time-averaged Lorentz force can be written

  * 2 1 B J F   real EM . (6.16)
It is possible to estimate the effect of electromagnetic force on the fluid velocity [11] from the screen parameter

2 0 2 L f R liquid w     , ( 6
.17)

and the Alfven velocity,

  2 1 0 max     B U A . (6.18)
The average velocity in bulk liquid can be computed from

4 1 6 . 0   w A R U U . (6.19)
In general, electromagnetic force density is highest at the area close to the crucible wall due to the skin effect. Thus, the fluid velocity should be much higher near the crucible wall. An example of the fluid convection caused by the electromagnetic force is shown in Fig. 6.2d. The electromagnetic convection is in many cases undesirable since it is difficult to control.

Physical model and boundary conditions for mass transfer

In contrast to the PVT modeling, the electromagnetic computation is coupled with the fluid dynamics. Indeed the electromagnetic, heat transfer, and fluid dynamics are coupled in the TSSG model while the modeling of mass transfer will be computed using the results from such coupled calculations. This is possible since both solutal Marangoni and solutal buoyancy convection can be neglected in this study as mentioned in 6.1. Moreover, the transport of diluted species is considered in this model. The governing equation describing the carbon transport is the steady state convection-diffusion equation written by

0 2       C C D u . (6.20)
In the growth of SiC using the TSSG process, the crucible is not only a container but also the carbon source, meaning thus that the crucible-liquid interface does not only induce mechanical stirring but is also the place of a complex reaction/dissolution of solute. The equilibrium concentration of carbon is imposed on both crucible-liquid and crystal-liquid interfaces The crystallization takes place when S > 0, while the dissolution takes place when S < 0. The incorporation of solute in the solid crystal should depend also on the interface velocity. However, if the supersaturation is very small and the growth rate is low, it is reasonable to assume a thermodynamic quasi-equilibrium at the solid-liquid interface. The growth rate can be computed form the diffusive flux of carbon atom normal to the seed surface: This equation specifies that the growth rate is proportional to the supersaturation and inversely proportional to the diffusion boundary layer thickness (δ D ). The distribution of carbon concentration in front of the seed (Fig. 6.10b) shows that the supersaturations are almost constant for all the rotation speed. On the other hand, the solute boundary layer thickness decreases when increasing the rotation speed. This is also the case for the area close to the seed periphery. Thus, in these systems, the growth rate was enhanced by the reduction of the boundary layer thickness due to the seed rotation. Note that this model does not consider the interface attachment kinetics which allows the concentration at the interface to be higher than the equilibrium value [13].

  n     C D M R SiC SiC g  . ( 6

Effect of frequency

The use of an external magnetic field has been recently proposed and theoretically investigated [14][15][16]. In [15], the authors separated both heating and electromagnetic mixing by computing respectively a resistive heating to control the temperature distribution and an AC field generated by an induction coil, the latter affects only the fluid velocity. With a specific design of the growth configuration [15], i. e., the combination of using one-turn coil placed above liquid surface and the low frequency (typically 50 Hz) was found to be more effective for carbon transport enhancement. This theoretical case is very promising but its technical implementation is not industrially straight forward. Thus, in this section, we are considering a technically relevant single coil system where the electromagnetic field is used for both heating and mixing.

We will mainly focus on the effect of coil frequency on the fluid flow and mass transport. The crystal temperature was kept constant at 1700˚C. The coil currents were adjusted to get the suited temperature when varying the frequency from 5 to 100 kHz. The distribution of Lorentz force density in liquid for different coil frequency is shown in Fig. 6.12. The Lorentz force is concentrated only at side wall of the container for 100 kHz while it extends to the inner part of liquid when the frequency decreases. For further analysis, the maximum value of the Lorentz force density in liquid and the Reynolds number computed from the average velocity and Alfven velocity are plotted as a function of the coil frequency in Fig. 6.13a and Fig. 6.13b, respectively.

The maximum Lorentz force density and the Reynolds number increase as the frequency decrease. azimuthal and radial components. If the circular seed is mounted to the graphite rod at the symmetry axis, the schematic representations of tangential flow is shown in Fig. 6.17. By considering the seed at rest, we define the fluid velocity relative to the seed velocity (Δv) which will be called from now on as the fluid flow or the fluid velocity. In order to investigate the interaction between step flow and fluid flow, the parallel and antiparallel flows are defined. The former refers to the case where the fluid flow and step flow have the same direction while the latter refers to the case where the flows oppose to each other. Fig. 6.17a represents the fluid flow with only azimuthal component (Δv ϕ >0, Δv r =0) while Fig. 6.17b is the one with only radial component (Δv ϕ =0, Δv r >0). For ϕ component, the positive and negative fluid velocities are defined for clockwise and counterclockwise directions, respectively. For r component, the positive value refers to the flow out of the seed center, while it is negative for the opposite flow direction. The boundary line is drawn to separate the surface area having the parallel or antiparallel flow. In Fig. 6.17a the boundary line is the horizontal line crossing the seed center. Defining this horizontal line parallel to the off-cut direction as a reference line, the boundary line will be deviating ξ degree from the reference line if the fluid flow has both azimuthal and radial components (Fig. 6.17c). This angle of deviation is derived from the trigonometry consideration:

                                        0 , v , v 0 , v , v tan v v tan 1 1 seed fluid seed r fluid r r z r z r z r z r r , (6.27)
where r is a distance from the seed center, and z is the position in liquid below the seed-liquid interface. This formula is verified for Δv r =0 where ξ is 0º (Δv ϕ >0) or 180º (Δv ϕ <0), and for Δv ϕ =0 where ξ is 90º (Δv r >0) or -90º (Δv r <0). The positive and negative values of ξ represent the rotation of the boundary line in clockwise and counterclockwise directions, respectively. In fact, the question arises in which z position in liquid (z fluid ) should be considered when using Eq. (6.27). If the ratio between the solute (δ D ) and momentum (δ v ) boundary layer thicknesses is much less than unity, the near solid-liquid interface can be considered as a motionless fluid and the mass transport in this region is governed only by the diffusion process [18]. This criterion may allow us to define the velocities at the diffusion layer thickness as representative for fluid velocities near the seed. The ratio between solute and momentum layer thicknesses is related to the Schmidt number (Sc) as [19] Since the model has been verified concerning the growth rate, it should be possible to use other calculated parameters for the qualitative analysis. This is the case for the angle of deviation ξ defined in Eq. (6.27). For the case of high rotation speed of 100 rpm in counterclockwise direction, the values of fluid velocity for both radial and azimuthal components relative to the seed are plotted along the radius (Fig. 6.21). Both radial and azimuthal components of fluid velocity increase along the seed radius. In such case, the solute boundary layer thickness is almost constant δ D =0.2 mm along the seed radius. The ξ value at r=6 mm is calculated to be around 28º. The phase parameter diagram for these growth conditions is plotted in Fig. 6.22 for both and (01-10) facets where the A-D areas are separated by the boundary lines. The surface of the as grown crystal using 4H-SiC (0001) seed crystal 4º off-axis towards [11][12][13][14][15][16][17][18][19][20] is compared with the diagram. The seed rotation speed was set to 100 rpm in counterclockwise direction (clockwise as observe from above the liquid surface). Each image was taken in the area r>4 mm where its azimuthal location on the crystal surface is defined by the angle θ relative to the [11][12][13][14][15][16][17][18][19][20] for -85º and 175º. However, if we neglect the radial flow such as the diagram for the purely azimuthal flow shown in Fig. 6.20, the behavior of macrosteps formation become poorly described. For example, for θ = -40º, Fig. 6.20 shows the only antiparallel flow while macrosteps close to (10-10) plane were experimentally observed. Another example is the case for θ=140º, Fig. 6.20 predicts only the parallel flows for both and (01-10) whereas only the macrosteps predominantly close to (01-10) plan were observed. Moreover, the deviations for the cases of -85º and 175º will become even larger. Thus the effects of the fluid flow in both components tangential to the seed have to be considered since they provide better description of the macrosteps orientations. The occurrence of the macrosteps increases the surface roughness and results in poor crystal quality. The improvement of the surface morphology roughness can be achieved by decreasing the area of parallel flow, such as changing the center of rotation at the seed [22].

However the formation of macrosteps may give some benefit owing to the dislocation conversion mechanism [23]. Thus the results from the combined numerical and analytical modeling show the possibility to use the phase parameters to control and optimize the quality of the SiC crystal grown from the TSSG process regarding the macrosteps formation.

Conclusion

We have discussed the modeling of the TSSG process in this chapter where the special attention has been paid to the contribution of the four main convections on the fluid dynamics.

The process parameters, namely operating frequency and crystal rotation speed, have been investigated. The electromagnetic convection is the main contribution governing the growth process at the experimental conditions using low frequency (14 kHz) while the combined

General conclusions and perspectives

In this thesis, we made use of the numerical modeling as a tool for "process characterization". This tool was able to provide an assessment of the physico-chemical parameters and their distributions inside the growth chamber. Such information is not possible or at least extremely difficult to obtain during the experiment. This is the case especially for the growth of SiC using PVT and TSSG techniques since the experiments were carried out in a quasi-closed graphite crucible, often called as a "black box" process. Thus the modeling approach can be used as an alternative way to make a link between the controllable parameters and the observable results providing better monitoring and better understanding of the growth processes.

With the modeling of induction heating and heat transfer, the temperature distribution in the growth reactor can be obtained. This is one of the crucial information for the crystal growers since the temperature gradient creates a driving force for the crystal growth. The model was verified by the similar tendency of the I-V curves, and the correlation between the diffusion path of the SiC vapor species and temperature gradient streamlines in the powder source. Concerning more details of the model, the heat transfer by gas convection was shown to be negligible due to the dominant radiative heat transfer at the high growth temperature. The semi-transparency of SiC was shown to strongly affect the temperature profile in the vicinity of the crystal and the temperature gradient. However, such effect was not taken into account in the modeling of our reactor for three main reasons. First, the crystals may become opaque since the experimental observations reveal that the doping level is quite high in the crystals even in the case of unintentionally doping experiments. Second, the present band energy model for taking into account a semi-transparency behavior is somewhat simplified where the interactions between the transmitted wave and the crystal were neglected. Finally, the discrepancies in an absorption coefficient dependent on the wavelength, temperature, doping level, and polytype can lead to completely different conclusions. These points should be examined in more details, possibly with a proper experiment and measurement design for the materials properties related issue.

The mass transport model combined the fluid dynamics and the species transport models.

In the PVT growth, special attention was paid to the boundary conditions and the chemical model at the solid vapor interfaces. Two approaches were proposed to take into account the stoichiometry condition at the growing surface. The comparative approach was able to provide the results respecting such condition. However, the first approach was not very convenient for further development of the model since it dealt with the kinetic parameters, sticking coefficient, which is very complex and difficult to determine and verify. The coupled approach provides better comprehension of the growth process since the prior assumption for the three-phase heterogeneous equilibrium was not required. The stoichiometry condition was successfully respected with additional information on the crystal chemistry. Modeling results were validated with good agreement with the experiments. More importantly, this approach really opened the way to get access to the physico-chemical parameters in the crystal, the activities of Si and C atoms in SiC solid solution, during the growth process and their dependency on the process parameters. The correlations between the state of the growing crystals and the occurrence of cubic or hexagonal polytype and defects density were made. The correlation with the doping concentration was also suggested as a possible extension of the model. Those results were discussed base on the evidences found in the literature while the real mechanisms are still unclear. Thus, the extension of the present model and further studies on the mechanisms responsible for any valid correlations are required.

In the full modeling of the TSSG process, the induction heating, heat transfer, and fluid dynamics were coupled. The heat transfer by the fluid convection was shown to be more efficient than the diffusion since the temperature profile in liquid clearly followed the convection pattern.

Such convection pattern was governed by the four main contributions: buoyancy convection, forced convection, Marangoni convection, and electromagnetic convection. The effects of the process parameters, namely, crystal rotation speed and operating frequency were studied. Four main parameters could be monitored during the growth: temperature, velocity field, concentration profile, and supersaturation. The electromagnetic convection was the main contribution governing the growth process using the low frequency while the combined buoyancy and Marangoni became dominant at high frequency. At the growth conditions, the enhancement of the growth rate could be obtained by increasing in the rotation speed of the seed.

This effect is mainly due to the reduction of the boundary layer thickness in front of the seed.

The calculated growth rate as a function of rotation speed exhibit good agreement with the experiments while some deviation at very high rpm was found. With the dimensionless analysis, such flow should be already in the turbulent regime. This is the case also for the growth using very low frequency (5 kHz). Thus, it will be necessary in the future to reconsider the model in the transient and/or three-dimensional system for higher accuracy. In addition, the combined numerical and analytical modeling was developed granting a comprehensive visualization of the interaction between the fluid flow and the step flow that affects the surface morphology of the crystal. The phase parameter was introduced to correlate the macrostep formation with the fluid flow. Such model was successfully validated with the experimental results. Apart from those achievements, the TSSG process is still needed to be optimized to fully utilize its capability for growing high quality SiC crystal. This will include the control of the fluid convection and suppression of the growth instabilities. The basic mechanism for the interaction at the crucibleliquid interface also remains to be studied. This will possibly provide more appropriate boundary condition, rather than the equilibrium concentration, to be implemented in the mass transport modeling.

The present thesis really emphasizes on developing a better understanding of the SiC growth processes. The complete experiment, characterization, and process modeling approach will be a key to obtain a full description of the processes. This obviously requires further extension of the model as well as additional experiments and characterizations.

   N i p i p i C x C 1 J/(kg•K). (A.3)
For specific heats of all gas species, the data are obtained from NIST-JANAF Thermochemical Tables. The tabulated functions have been extracted from these discrete temperature dependent data sets in order to use in the simulation. The thermal conductivity of solid materials is written in forms of analytic function depending on the temperature. In the gas phase, the Mason-Saxena-Wassiljewa expressions have been used. This MSW method gives the accuracy of ±6-8%. The thermal conductivity of the gas mixture using MSW mixing rules is in the forms of

      N i N j ij i i i x k x k 1 1 , (A.4)
where Wassiljewa function ( ij  ) can be calculated by Mason and Saxena expression written as

2 4 1 2 1 2 1 1 1 8 1                                          i j j i j i ij M M k k M M . (A.5)
For each gas species, the thermal conductivity can be calculated from the Chapman-Enskog's formulae based on the kinetic theory of gases assuming Lennard-Jones potentials:

k i i i M T k     2 3 10 6318 . 2  W/(m•K). (A.6)
The collision integral ( k  ) represented the deviation from rigid sphere behavior is defined by

                i i i i k T T T T .
              1 3 32 3 , (A.8)
where d is the equivalent average diameter of the powder granules.

A.2.2 Fluid flow properties

The fluid flow problem needs the density and the dynamic viscosity of the gas mixture. The former case is the same property used in the heat transfer model. The latter one, the dynamic viscosity has been calculated utilizing again MSW mixing rules:

      N i N j ij i i i x x 1 1   . (A.9)
In this case, Wassiljewa function for dynamic viscosity is in the forms of

2 4 1 2 1 2 1 1 1 8 1                                        i j j i j i ij M M M M   . (A.10)
Chapman-Enskog's formulae based on the kinetic theory of gases assuming Lennard-Jones potentials are implemented for calculating the dynamic viscosity of the i th species:

       2 5 10 441 . 8 i i i T M kg/(m•s), (A.11)
where

k    
as defined in Eq. (A.7).

A.2.3 Mass transport properties

In this mass transport model, the diffusion coefficients are required. There are two types of diffusion coefficients according to two different mass transport mechanisms considered here: thermal diffusion and isothermal diffusion. The thermal diffusion coefficients of all gas species are very small in PVT growth. Their values can be found in an order of 10 -10 (kg m -1 s -1 ) in our growth process. Thus, the isothermal diffusion mechanism dominates. The multicomponent Maxwell-Stefan diffusion coefficients are used in the mass transfer equations. In our case, we assume the low-density gas mixtures condition in which Maxwell-Stefan diffusivities can be replaced with the binary diffusivities for the species pairs presented in the mixture. Thus the binary diffusivities will be used as the input data to compute multicomponent Fick diffusivities using matrix transformation. The diagonal terms of Maxwell-Stefan diffusivity matrix defined by ii D representing the diffusion of species through itself are equal to unity. The off-diagonal terms ij D representing the mutual-diffusion of species i at infinite dilution through species j equal to the diffusion of species j at infinite dilution of species i. Thus, the binary diffusivity matrix is

symmetric   ji ij D D 
. The binary diffusion coefficients ij D for the gas mixture are given by

                 ij D ij j i ij T P M M T D 2 3 4 1 1 10 876 . 5 m 2 /s, ( A12
)
where In the TSSG process, the diffusion coefficient of C in liquid Si is set to

  2 j i ij      . The function     ij D T is defined by                     ij ij ij ij

A.3 Properties of various materials used in the modeling

In some cases, the materials venders provide the property data of their products.

However, those data do not cover the entire range of our experimental conditions, i.e., the lack of data at high temperature. Thus it is necessary to extrapolate those properties and/or compare the given properties with the models presented in the literatures in order to utilize the appropriate ones as well as to improve the accuracy of the model and minimize the error that may occur.

Physical properties of materials used in the simulation are listed as follow.

Insulator

Electrical conductivity:



Graphite crucible

Electrical conductivity: Marangoni coefficient:

                     
4 10 5 . 2      Si K m N 

Vapor species

Species

Molar mass (g/mol) Dans le but d'obtenir des solutions réalistes sans simplifications excessives, plusieurs méthodes numériques de discrétisation des équations aux dérivées partielles ont été proposés, par exemple, la méthode des différences finies (MDF), la méthode des volumes finis (MVF) et la méthode des éléments finis (MEF). Même si les approches de résolution des équations aux dérivées partielles sont différentes, elles partagent des caractéristiques communes importantes. Tout d'abord le domaine continu est discrétisé en un ensemble de sous-domaines discrets appelés éléments dans lesquels leurs lois de comportement sont définis par un nombre fini de paramètres.

σ i (Å) ε i /k B (K) D i T (kg/(m•s)) Si 28
De plus, les fonctions approchées sur chaque éléments doivent respecter les conditions de continuité entre les différents éléments du domaine. Les avantages et les inconvénients de l'utilisation de ces méthodes numériques sont synthétisés dans le tableau 1. La méthode des différences finies est relativement facile à programmer, mais sa capacité est limitée à des géométries simples et une grille régulière. La méthodes des éléments finis permet de traiter de façon naturelle les conditions aux interfaces entre deux milieux différents, même si cette méthode est moins triviale pour la mécanique des fluides que la méthode des volumes finis. Étant donné que la croissance des cristaux est principalement associée à des phénomènes aux interfaces et donc mathématiquement à des problèmes de conditions aux limites, la méthode des éléments finis est plus adaptée pour le sujet de la thèse de doctorat. Tout au long de cette thèse, un logiciel commercial basé sur la méthode des éléments finis servira comme un outil de résolution du système d'équations aux dérivées partielles, associé à un ensemble de conditions aux limites régissant le procédé de croissance de SiC.

Chauffage par induction et transfert de chaleur

La distribution de température dans la cavité de croissance a été étudiée en résolvant les équations décrivant le chauffage par induction couplées aux transferts dans l'ensemble des régions du four. Tous les calculs ont été menées par méthode des éléments finis (FEM) dans l'hypothèse d'une géométrie axisymétrique en deux dimensions, puisque le réacteur de croissance peut-être considéré de symétrie axiale. Le chauffage par induction électromagnétique est une méthode pratique pour le chauffage haute température du creuset. Dans ce régime de chauffage par induction RF, le phénomène électromagnétique est décrit par les équations de Maxwell qui permettent de calculer la distribution de champ électromagnétique. L'équation pour résoudre le problème électromagnétique peut être réduite à l'équation de la magnétodynamique en potentiel vecteur :

    e r r i J A v A A                   0 2 0 1 . (1)
Les courants de Foucault induits sont générés selon la loi de Lenz. Les courants induits circulant dans les parties résistives (conductrices) provoque une dissipation d'énergie nommée effet Joule.

Le creuset est chauffé en conséquence par effet de Joule, la densité de puissance fournie peut être formulée par :

  * 2 2 1 A A J E        Q . (2)
La densité de pertes Joule servira comme source dans l'équation de transfert de chaleur. Trois mécanismes sont considérés : la conduction, la convection et le rayonnement. L'équation de l'énergie utilisée pour calculer la distribution de température dans le réacteur de croissance est de la forme :

    . Q T k T C p         v (3)
En outre, la conservation de la i ème espèce dans un mélange gazeux est donnée par l'équation de bilan des espèces :

  0         i i v j . (8)
Le flux molaire N i aux frontières réactives est décrit par la somme de deux contributions :

v j N i i M ρx M i i   . (9)
Le premier et second terme sont respectivement les flux molaires diffusif et convectif. À l'interface solide-vapeur, le changement de volume important causé par les réactions de sublimation ou de la condensation produit un flux d'advection, nommé aussi vitesse de Stefan et s'écrit de la façon suivante :

n v              1 1 N i i i N M PM RT . (10)
Ici, la somme comprend uniquement les espèces réactives, car le gaz inerte ne participe pas aux réactions hétérogènes. Les équations ( 9) et (10)

  SiC C Si p a a a T K   1 . ( 12
)
La constante d'équilibre K 1 p est écrite pour x et y proche de l'unité, puisque le cristal de SiC est proche de la stoechiométrie. L'activité des composés stoechiométriques est égale à l'unité par définition de l'état standard, donc a SiC = 1. La constante d'équilibre K 1 p est calculée à partir des enthalpies libres provenant de la base de données thermochimique JANAF. La formation de trois espèces gazeuses principales à l'interface cristal-vapeur peut être attribuée aux réactions hétérogènes suivantes :

Si (in SiC (s, s) ) = Si(g) with   Si Si p a P T K * 2  , (13) 2Si
(in SiC (s, s) ) + C (in SiC (s, s) ) = Si 2 C(g) with   C Si C Si p a a P T K 2 * 3 2  , (14)
Si (in SiC (s, s) ) + 2C (in SiC (s, s) ) = SiC 2 (g) with

  2 * 4 2 C Si SiC p a a P T K  . (15)
L'équation de bilan massique exprimant l'incorporation stoechiométrique d'atomes de Si et C dans le cristal en croissance, ferme l'ensemble des équations : En ce qui concerne la description quantitative de la croissance des marches, nous introduisons le "paramètre de phase" :

2 2 2 2 2 2 C Si C Si C Si C Si Si N N N N N       . ( 16
         sin , (17)

Conclusion

Dans le cadre de cette thèse, nous avons utilisé la modélisation numérique comme un outil pour la caractérisation de procédés. Cet outil a été en mesure de fournir une évaluation des paramètres physico-chimiques et leurs distributions dans la cavité de croissance permettant d'améliorer le contrôle du procédé de croissance et de mieux les comprendre.

Avec la modélisation du transfert de chaleur et des phénomènes électromagnétiques, la distribution de température dans le réacteur de croissance peut être obtenue. C'est l'une des informations fondamentales pour les fabricants de cristaux car le gradient de température créé une force motrice permettant la croissance des cristaux. Le modèle de transport de masse est combinée avec la mécanique des fluides et le transport des espèces gazeuses. Dans le cadre de la croissance en phase vapeur (PVT), une attention particulière a été portée sur les conditions aux limites et le modèle chimique aux interfaces solide vapeur. En utilisant la modélisation numérique avec une approche couplée et la thermodynamique, la condition de stoechiométrie a été respecté avec succès, donnant des informations supplémentaires sur la chimie des cristaux.

Les résultats de la modélisation ont été validés et un bon accord avec les expériences a été obtenu.

Plus important encore, cette approche a vraiment ouvert la voie pour accéder aux paramètres physico-chimiques dans le cristal et plus précisément les activités des atomes de Si et de C en solution solide de SiC, durant la phase de croissance et de leur dépendance sur les paramètres du procédé. Des corrélations entre la qualité des cristaux et l'apparition de polytype cubique ou hexagonal et la densité de défauts ont été proposés. La corrélation avec la concentration de dopage a également été suggérée comme une possible extension du modèle. Ces résultats ont été discutés sur la base des articles trouvés dans la bibliographie, tandis que les mécanismes réels ne sont pas encore clairs. Ainsi, l'extension du modèle actuel et d'autres études sur les mécanismes responsables des corrélations sont requis.

Pour la modélisation complète du procédé de croissance en phase liquide (TSSG), il est nécessaire de coupler les phénomènes électromagnétiques, les transferts thermiques et la mécanique des fluides. On a étudié les effets des paramètres du procédé tel que la vitesse de rotation de cristal. Quatre paramètres principaux peuvent être étudiés en cours de la croissance : la température, le champ de vitesse, le profil de concentration et la sursaturation. La convection électromagnétique à basse fréquence est la principale contribution régissant le procédé de croissance. Pour des conditions de croissance donnée, on peut améliorer la vitesse de croissance en augmentant la vitesse de rotation du substrat. Le calcul de la vitesse de croissance en fonction de la vitesse de rotation est en bon accord avec les expériences.

De plus, l'association de la modélisation analytique et numérique a été développée permettant une meilleure compréhension de l'interaction entre l'écoulement du fluide et la croissance par avancée de marches qui affecte la morphologie du cristal. Le paramètre de phase vise à établir une corrélation entre la formation des macro-marches avec l'écoulement du fluide.

Ce modèle a été validé avec succès avec les résultats expérimentaux. En delà de ces avancés, le procédé de croissance TSSG doit être amélioré afin d'obtenir une croissance de monocristaux de SiC de haute qualité.

186 Il devra permettre le contrôle de la convection du fluide et la suppression des instabilités de croissance. Le mécanisme de base décrivant l'interaction à l'interface liquide-creuset reste à étudier, afin de donner des conditions aux limites plus appropriées (plutôt que d'imposer concentration d'équilibre) dans le cadre de la modélisation du transport de masse.

La présente thèse contribue au développement d'une meilleure compréhension des procédés de croissance de SiC. L'utilisation simultanée de l'expérience, la caractérisation et la modélisation numérique est un facteur clé pour obtenir une description complète des procédés. Il est nécessaire de poursuivre les extensions au modèle existant en associant des expériences et des caractérisations additionnelles.

Summary of the thesis

Silicon Carbide (SiC) is one of the most desirable materials for power electronic devices. The development of the growth process, to achieve larger size and higher quality is on the way. Even if modeling has proved its ability to assist the optimization of the growth processes, there are still some strong issues which are not considered in a satisfactory way. This thesis aims to use the modeling tools to tackle those challenging fundamental and technological issues on both industrially used PVT and emerging TSSG processes.

In the PVT process, special attention is paid to the physical chemistry at the solid-vapor interfaces. Especially, we investigated the way to model the stoichiometric crystallization of SiC from a non-congruent vapor as the Hertz-Knudsen model was shown to be not adapted. We thus considered SiC as a solid solution using coupled mass transport modeling and thermodynamics. This approach gives an assessment to the chemistry of the SiC crystal which can be linked to the controllable parameters. Such correlations may serve as a basis to control the points defect density, stable polytypes, and doping concentration. Concerning the TSSG process, the effects of crystal rotation speed and operating frequency are studied. The electromagnetic convection is the main contribution governing the growth process using low frequency while the combined buoyancy and Marangoni convections become dominant at high frequency. In the experimental conditions using low frequency, the crystal growth rate could still be enhanced by increasing the rotation speed. The phase parameter is also introduced using the combined fluid dynamics and analytical modeling. This provides a comprehensive visualization of the interactions between fluid flow and step flow and a guideline to improve the surface morphology of the crystal.

Résumé de la thèse en français

Le carbure de silicium (SiC) est l'un des matériaux les plus prometteurs pour les dispositifs électroniques de puissance. Même si la modélisation a prouvé sa capacité à assister le développement du procédé de croissance des cristaux de SiC, de nombreux aspects ne sont toujours pas décrits de façon satisfaisante. Cette thèse a pour but d'utiliser les outils de modélisation pour étudier les phénomènes fondamentaux ayant lieu dans la croissance massive du SiC, que ce soit en phase vapeur (PVT) ou en phase liquide (TSSG).

Pour la méthode PVT, une attention particulière est portée sur la physico-chimie de l'interface solide-vapeur. Pour simuler la cristallisation stoechiométrique du SiC à partir d'une phase gazeuse noncongruente, phénomène mal décrit par le modèle d'Hertz-Knudsen, nous avons considéré le SiC comme une solution solide en utilisant la modélisation couplée du transfert de masse et de la thermodynamique. Cette approche donne une évaluation de la composition du cristal de SiC qui peut être liée à des paramètres contrôlables. De telles corrélations peuvent servir de base pour contrôler la densité de défauts ponctuels, la stabilité des polytypes et la concentration de dopage. Dans le cas du procédé TSSG, les effets de la vitesse de rotation du cristal et de la fréquence du champ magnétique sont étudiés. La convection électromagnétique est la principale contribution régissant les mouvements de fluide à basse fréquence alors que la convection naturelle et l'effet Marangoni deviennent dominants à haute fréquence. Dans les conditions expérimentales utilisant les basses fréquences, la vitesse de croissance du cristal pourrait encore être accrue en augmentant la vitesse de rotation. Une modélisation analytique, couplée aux calculs de dynamique des fluides a permis de décrire les interactions entre le flux de liquide et la direction de propagation des marches de croissance à la surface du cristal. Un paramètre de phase a été introduit comme critère d'apparition de macromarches, à l'origine d'instabilités morphologiques.

from 7

 7 different countries. In total, 12 early stage researchers (ESR) and 4 experienced researchers (ER) were recruited and trained with the network. The consortium is divided into 3 technical work-packages: o WP1 (Material growth and related aspects) is dedicated to the development of less mature polytypes than 4H (3C and 15R). o WP2 (Characterization of material and functional interfaces) is in charge of studying the properties of the materials/surface and interfaces. o WP3 (Devices and demonstrators) is in charge of the electrical testing and fabrication of the targeted devices.

  Fig. 2.1 T

Fig. 2

 2 Fig. 2.7 S Lely met

  Fig. 2.9 D the contr

  Fig. 3.1. In o n u in the app

  Fig. 3.2 T

N

  are the number of volume element and surface element, respectively. The number of degrees of freedom DOF N equals to the number of node (NBN) multiplied by the number of unknown per node. Equations (3.13) can be written as a set of linear equations of the form

  Fig. 3.6 D

  Fig. 3.8.

  Fig. 3.10

Fig. 4

 4 Fig. 4.1 D heat tran isotherm

Fig. 4

 4 Fig. 4.3 C

Fig. 4

 4 Fig. 4.4 S

  Fig. 4.5 T at the init

  Fig. 4.10 and react

4 . 14 .

 414 The others are shifted up either by 1 or 2 cm. In all cases, the temperature on top of the crucible is fixed at 2200 °C. The isothermal lines and the temperature gradient streamlines are shown in Fig.4.16. For the reference coil position, the hot point locates at the bottom corner of the crucible containing SiC powder source (Fig.4.16a). The temperature gradient streamlines direct upward from the hot point to the cold point (top parts of the crucible).

Fig. 4

 4 Fig. 4.16 crucible plotted a

  Fig. 4.16 ar t of the hot n computed vations. Onc art of the po he powder w finally bec art of the po anisms. First esulting from the vapor sp p part where

. 1 ) 1 . 2 .

 112 where is λ the molecular mean free path and L C is the characteristic length of the fluid flow. The flow behavior is classified by four different flow regimes [1]: Continuum flow (Kn < 0.01) Slip flow (0.01 > Kn > 0.1) 3. Transition flow (0.1 > Kn > 10) 4. Free molecular flow (Kn > 10).

Fig. 5

 5 Fig. 5.1 K the grow be the d paramete

  Fig. 5.2 S thin Knu

Fig. 5 . 3 .

 53 Fig. 5.3. The overall species distribution looks very similar for the model without the stoichiometry condition (no constraint) and the one with the atomic fluxes balance. However, in the region very close to the seed surface, there are the distortions of the mole fraction curve when using atomic fluxes balance. The mole fraction gradients of Si rich species (Si, Si 2 C) decrease while the one of the C rich species (SiC 2 ) increases. At the seed surface, the mole fractions of all reactive species are unchanged compared to the case without the stoichiometry condition. Thus from (5.20), only the diffusive fluxes are adjusted in order to respect the stoichiometric constraint but not the convective ones and the Stefan's velocity remains unaffected.

Fig. 5

 5 Fig. 5.3 D the symm interface

5 . 5 .

 55 Fig. 5.5 T stoichiom

  Fig. 5.10 C system (d) SGPS

39 )

 39 equilibrium state in SiC with excess C (SiC-C) gives a C =1 with the minimum value of a Si while a Si =1 with minimum value of a C for SiC with excess Si. In the non-stoichiometry domain of SiC compound (just before the precipitation of either C or Si at the surface), the Si and C activities vary between SiC-C and SiC-Si phase limits at a given temperature. The formation of the three main vapor species can be attributed to the following reactions: Si (in SiC (s, s) ) = Si(g),(5.34) 2Si (in SiC (s, s) ) + C (in SiC (s, s) ) = Si 2 C(g),(5.35) Si (in SiC (s, s) ) + 2C (in SiC (s, s) ) = SiC 2 (g).(5.36)The corresponding equilibrium constant for reactions (5.34)-(5.36) are respectively Since the standard state defines a SiC =1, a C can be written in terms of a Si (5.32) and the equilibrium partial pressure of each vapor species is a function of a Si and its equilibrium constant (temperature dependent).In order to couple the thermodynamic mechanisms with the mass transport modeling we consider separately the chemistry of the solid phases at the source and seed surfaces. On the source surface, the SiC compound sublimes non-congruently and probably the source will be soon a mixture of SiC-C. Thus a C =1 and the equilibrium partial pressure of the vapor species is depending only on the temperature. On the seed surface, the incident mass fluxes sublimed from the source are transported by diffusion and convection to the seed surface. A partial condensation due to the lower temperature occurs and this condensation naturally produces a SiC-Si system due to the large excess incident Si flux. Then the condensed system sublimes more Si than C which means that the excess Si can be re-sublimed toward the source. The growing crystal remains closely stoichiometric and the amount of back flow or the re-sublimation fluxes depend on both the temperature of the seed and the incident fluxes. Thus in a steady-state regime for the growth at the surface of SiC seed, a mass balance equation (5.21) expressing the stoichiometric condition for the incorporation of Si and C atoms is imposed. From (5.17), the SiC stoichiometry condition (5.21) becomes

  Fig. 5.14 distribute condition

  Such model was suggested based on the time variation of the concentration of carbon vacancies C V N in a transition layer. Averaged data concerning the concentration of C V N [41] reveals that C V N in the 4H, 15R, 6H, and 3C-SiC are 7.3×10 20 ,15.1×1020 , 16.3×1020 , and 33.6×10 20 cm -3 , respectively. The introduction of carbon vacancies due to an excess of Si results in the compression of the crystal lattice. Thus a cubic structure of the growing layer is more favorable than the hexagonal structure in term of energy. If the C deficiency is dominant in the crystal, our model predicts the carbon vacancies to be maximized in the growth condition close to SiC-Si phase limit. Such growth condition should be favorable for the formation of cubic polytype as discussed earlier. Thus our model is in very good agreement with the vacancy model of the heteropolytype epitaxy of SiC[42]. The extra information accessible by our model is the variation of the activities with the growth conditions.

  deposition. Thus a heterogeneous (C-vapor) reaction occurring at the graphite wall is assumed to be Si(g) + 2C(s) = SiC 2 (

  Fig. 5.17 show the show the the graph

  .21. The highest temperature at the bottom corner of the growth chamber brings about the highest sublimation rate and large Stefan flow in front of such area. The main part of the fluid velocity is directed from the source to the seed. Part of the fluid goes toward the bottom part of the graphite guide due to the Stefan flow caused by the formation of the polycrystal. The distribution of partial pressure of all the vapor species are shown in Fig. 5.22. The accumulation of the SiC species has a highest value close to the powder source. Thus the diffusive flux of SiC

Fig. 5 . 21 Fig

 521 Fig. 5.21

. 2 )

 2 Fig. 6.1 S section fo

. 13 )

 13 replaces the condition on the tangential component of fluid velocity from the Slip boundary condition where ϕ is the test function for the tangential component of fluid velocity at the surface. On the liquid surface, if the temperature is higher at the crucible-liquid interface compared to the area near the seed crystal, the Marangoni force will cause the movement of fluid from the crucible wall toward the seed. The Marangoni convection for such case is shown in Fig.6.2c. In order to determine the intensity of the flow affected by the Marangoni force, the Reynolds-Marangoni number can be used:

  Ceq is the mole fraction of carbon. Similar to the PVT process, the dissolution and crystallization areas can be evidenced by the supersatuation defined from the deviation of the concentration from equilibrium

  Fig. 6.3 contours

Fig. 6

 6 Fig. 6.12 (b) 25 kH

  Fig. 6.17 arrows re flows are parallel a

  direction. The formation of macrosteps (the dark lines formed by the accumulations of several microsteps) is correlated to the phase parameter diagram. The parallel flow for (01-10) facet (area A) give rises to the macrosteps formation predominantly close to (01-10) plane. The parallel flow criteria in both (01-10) and (10-10) facets (area B) result in possibility to form the zigzag-shape macrosteps while the main feature is still governed by the dominant phase parameter. The predominant macrosteps close to (10-10) plane are observed in the C area where the parallel flow condition is fulfilled only for (10-10) facet. In area D, only the antiparallel flow is predicted and the macrosteps should not be formed. The phase parameter diagram shows good agreement in correlating the flow conditions with the macrostep formations with some deviation

Fig. 6

 6 Fig. 6.22 behavior analytica taken on

  buoyancy and Marangoni dominates at high frequency (100 kHz). The enhancement of the growth rate at the experimental conditions is possible by increasing the rotation speed of the seed. The growth rate dependent on the rotation speed obtained from the calculations is successfully validated with the experiments. The use of very low frequency (5 kHz) and very high rotation speed (200 rpm) may cause the time-dependent and/or turbulent flow and results in the instability of the growth process which is not preferable for the growth even if it provides higher growth rate. Moreover, the combined numerical and analytical modeling has been developed to give a comprehensive visualization of the interactions between the fluid flow and the step flow that have effects on the surface morphology of the growing crystal. The phase parameter has been introduced to correlate the macrostep formation with the fluid flow. Such model has been successfully validated with the experiments. 145 Chapter 7

  Fig. 1 La

  Fig. 5 (a) résistives

∆T

  Fig. 7 (a expérime

  où θ est l'angle défini dans un cercle unité et sert à définir la position sur le cristal.  est l'angle entre l'orientation des marches et la direction de désorientation où les valeurs positives et négatives indiquent que l'angle tourne respectivement dans les sens horaire et antihoraire. Un calcul trigonométrique nous donne l'angle de déviation ξ : la distance depuis le centre du germe et z est la position dans le liquide sous l'interface liquide-germe. Le paramètre de phase varie entre -1 et 1, où les valeurs positives et négatives correspondent respectivement aux écoulements parallèles et antiparallèles. La valeur devient nulle si le la direction d'avancée de marches et la direction de l'écoulement du fluide sont perpendiculaires. Pour un train de marches orienté dans la direction [11-20], α devient nul. En ce qui concerne la symétrie d'ordre 6, représentant le système hexagonal de SiC (Fig. 16 a), α est égal à pour les facettes (01-10) -30º et les facettes (10-10) 30º. Dans les expériences de croissance, un germe avec un angle de quelques degrés de désorientation (0001) vers la direction [11-20] a été utilisé fournissant l'homoepitaxy par l'intermédiaire de la croissance par avancée de marches. Pour la croissance de monocristaux de SiC sous condition de sursaturation faible, la vitesse de croissance moyenne est plus grande sur les plans {11-2n} que sur les plans {1-10n}. C'est le cas pour les conditions typiques de croissance pour le procédé TSSG effectué dans des conditions proches de l'équilibre thermodynamique. En outre les marches à la surface sont concentrées surtout pour le procédé de croissance en solution. Donc au lieu d'avoir une seule marche et un train de marches droit et régulier, une morphologie possible de la surface peut être constitué de recouvrement de marches et devenir des marches en zig-zag avec des fronts de marche dans les directions (10-1n) et (01-1n) considérant que la forme cristalline est limitée par les fronts de marche lents. Dans ce travail, seul le modèle bidimensionnel est considéré. Ainsi, les surfaces de recouvrement (facettes) dans les directions [10-1n] et [01-1n] peuvent être projetées sur le plan c (0001). Cela nous permet d'envisager l'orientation des marches sur le plan c tel que les plans (10-10) et (01-10) représentent respectivement les facettes (10-1n) et (01-1n). Les résultats pour des écoulements azimutaux et radiaux positifs sont indiquées respectivement sur la Fig. 16 b et Fig. 16 c. Quatre domaines A-D sont séparées par des lignes de transition. Les conditions d'écoulements parallèles sont remplies pour les deux facettes {1-100} quand la valeur de Ψ de ces deux facettes sont positifs (zone B). Les facettes (10-10) ou (01-10) permettent de compléter respectivement, les résultats pour un écoulement parallèle dans les zones C et A. Seul un écoulement antiparallèle a été obtenu dans la zone D pour les deux facettes {1-100}.

  Fig. 16 interactio circulaire la direct respectiv

  pour les deux facettes (10-10) et (01-10), où les zones de A à D sont séparés par des lignes de transition. La surface du cristal obtenu à partir d'un germe 4H-SiC (0001) avec 4º de désorientation vers [11-20] est comparée avec le diagramme. La formation de macro-marches (les lignes sombres formées par accumulation de plusieurs micro-marches) est en corrélation avec le diagramme des paramètres de phase. Le critère d'écoulement parallèle pour la facette (zone A) (01-10) permet la formation de macro-marches principalement à proximité du plan (01-10). Le critère d'écoulement parallèle pour les deux facettes (01-10) et (10-10) (zone B) et provoque la possibilité de former des macro-marches en zig-zag tandis que la caractéristique principale est toujours régie par le paramètre de phase dominant. Les macro-marches plus nombreuses à proximité du plan (10-10) sont observées dans la zone C où la condition d'écoulement parallèle est remplie uniquement pour les facettes (10-10). Dans la zone D, seulement un écoulement antiparallèle est prévu et la morphologie de surface est beaucoup plus lisse que dans les autres zones. Le diagramme de paramètre de phase montre une bonne corrélation entre les conditions d'écoulement et la formation des macro-marches.

Table 3 .1 3.1 Over

 3

	T
	seeking f
	a domain
	are linear

   and i are the Lennard-Jones coefficients of the i th species. Moreover, in this work, we have also taken into account the influence of source powder porosity on its thermal conductivity. The heat transfer by radiation in the powder bulk is integrated with the heat transfer by conduction. Thus the effective thermal conductivity is estimated by

	k	eff		k	gas		e	T	d	  solid k
					16145 . 1 14874 . 0	.7732 0 52487 . 0 exp	exp	43787 16178 . 2 2	,	(A.7)
	where	i T		kT		i

*

. i

  est le flux molaire net de la i ème espèce gazeuse participant à des réactions hétérogènes, P i est la pression partielle P i * est la pression partielle d'équilibre, M i est la masse molaire, R est la constante des gaz parfait, T est la température, γ i et α i sont respectivement les coefficients de collage et d'évaporation. La pression partielle de l'espèce réactive peut être calculée par la modélisation globale du transfert de chaleur et du transport de masse. Les trois phases en équilibre hétérogène (SiC-C-gaz avec variance=0) ne sont pas pris en compte dans ce travail, donc la pression partielle d'équilibre ne peut pas être calculée directement à partir de la température de la phase solide (SiC-gaz avec variance = 1) parce que la SiC a été considérée comme une solution solide (s, s). Les atomes de Si et de C peuvent être considérés dans le composé non-stoechiométrique Si x C y comme à

	où N
	poudre sublime produisant les trois espèces gazeuses principales : Si, Si 2 C, SiC 2 . Ces espèces
	gazeuses seront transportés à la fois par diffusion et convection. Ensuite l'espèce gazeuse sera à
	la fois condensé (ou adsorbé) et sublimé sur le substrat. Ainsi, le taux de dépôt net à l'interface
	de cristal-vapeur du substrat est la différence entre les flux adsorbés et désorbés. A l'interface le
	flux molaire de l'espèce réactive est calculé à partir du modèle de Hertz-Knudsen et peut être
	écrit de la façon suivante :
	N	i		N	ads i		N	des i		 	i	P i			i	P i			M  2 1 i	RT	(11)

montrent qu'il existe un couplage entre l'écoulement du fluide et le transfert de masse puisque la vitesse du gaz à la frontière de la réaction doit être calculé à partir du flux net des espèces.

Le gradient de température est calculé à partir de la température à la surface libre de la poudre de SiC et la température du substrat. La température à la surface de la poudre est supérieure à celle du substrat et créé ainsi une force motrice pour la croissance du SiC. À haute température, la i l'équilibre : Si x C y = xSi (in SiC (s, s) ) + yC (in SiC (s, s) ) with

Good agreements were obtained. The numerical modeling was also shown to be able to assist the crystal grower to get access to the physico-chemical parameters and minimize the reactor dependent issue.

According to the definition of the Stefan flow, an additional or removal of species should influence the mean flow of the fluid next to the interface. With the stoichiometry conditions, the net total fluxes at the boundary are changed. As a consequence, the Stefan velocity and the convective fluxes are affected. The method of adaptable sticking coefficients is in agreement with this argument but not the atomic fluxes balance one. In the latter case, only the diffusive fluxes are affected but not the convective fluxes. Further argument concerning the weakness of such method can be discussed as followed. We assumed that the sticking coefficient is identical to the evaporation coefficient obtained from the high temperature mass spectrometer studies of the SiC sublimation process. In such process, the SiC powder would soon become SiC-C diphasic mixture. Thus the obtained evaporation coefficients are associated with the sublimation from SiC-C but not from pure SiC. It was shown that the evaporation coefficients decrease due to the increase in the carbon content at the SiC powder grains [21]. As a result, a set of sticking coefficients from our assumption may not be appropriated even for the initial system without constraint. The adaptable sticking coefficient method allows us to take into account the stoichiometry condition with is not included from the evaporation coefficient taken from the evaporation of SiC-C diphasic mixture itself. Thus, the adaptable sticking coefficients method appears to be a better description concerning the kinetics mechanism of the crystal growth than the atomic fluxes balance method.

The analysis of thermodynamic calculation results aims at answering to the question: which of Mod2 or Mod1 models is more appropriated? By comparing between the species partial contributions obtained from the global modeling in Fig. 5.6 and the species partial contribution in Fig. 5.9, the species partial contribution obtained from Mod1 are quite similar to the results obtained from the thermodynamic calculations at high T 2 /ΔT. When decreasing T 2 /ΔT, the numerical results start to deviate from the thermodynamics results. However, only Mod1 can preserve the general trend where the partial contribution of SiC 2 is much larger than Si. Thus, Mod1 gives the best agreement with the thermodynamic calculations for all the growth conditions. We may note that, by changing the growth conditions (T 1 and ΔT) using Mod1 for the numerical modeling, the partial contributions of Si and SiC 2 can be decreased, while the one of Si 2 C can be increased following the tendency shown in Fig. 5.9.

Appendix A

Physical properties

A.1 Electromagnetic properties

The physical properties related to the electromagnetic problem are magnetic permeability, electrical current density, electrical conductivity, and angular frequency. The permeability of all materials is set equal to that one of the air (1.2566×10 -6 H/m). Thus the relative permeability of all material equals to unity. The electrical current density is dependent on the desired growth temperature. The electrical conductivity varies in each material and is dependent on the temperature.

A.2 Properties of the gas phase

A.2.1 Thermal properties

Several models have been presented for the determination and estimation of physical properties of materials used in the growth process. For the gas phase, the simulations were performed assuming an ideal gas. Thus the gas density is written by

In the mixture of N species (growth cavity), the mean molar mass can be calculated from

Similarly, the specific heat of gas for a mixture of N species is calculated according to the molar fraction of each species reads

List of symbols

The symbols used in this thesis will be listed in the table below. The vector quantities are shown in bold. Some symbols that are not listed here can be found in the text.

Symbols