- Email: [email protected]
Numerical study of SiC CVD in a vertical cold-wall reactor
Computational Materials Science 24 (2002) 520–534 www.elsevier.com/locate/commatsci
Numerical study of SiC CVD in a vertical cold-wall reactor A.N. Vorob’ev a,*, S.Yu. Karpov a, M.V. Bogdanov a, A.E. Komissarov a, O.V. Bord a, A.I. Zhmakin b, Yu.N. Makarov c b
a Soft-Impact Ltd., P.O. Box 33, 194156 St. Petersburg, Russia A.F. Ioffe Physical Technical Institute, Russian Academy of Sciences, Polytechnicheskaya 26, 194021 St. Petersburg, Russia c STR Inc., P.O. Box 70604, Richmond, VA 23255, USA
Received 7 August 2001; received in revised form 29 January 2002; accepted 1 February 2002
Abstract The conventional heat and mass transport model is extended to describe silicon cluster formation in the gas phase and is employed for a numerical analysis of SiC chemical vapor deposition in a commercial vertical rotating disc reactor. The model is verified by comparing the computed growth rate with available experimental data. The growth rate is studied as a function of precursor flow rates varied in a wide range of values. It is found that the growth rate is limited by the gas mixture depletion in silicon atoms due to homogeneous nucleation. The secondary phase formation on the growing surface is analyzed. The SiC growth window depending on the precursor flow rates is calculated, and a significant influence of the homogeneous nucleation on the window width is shown. The model results predict that the Si/C ratio on the wafer can considerably differ from that at the reactor inlet. 2002 Elsevier Science B.V. All rights reserved. Keywords: Chemical vapor deposition; Cluster; Growth rate; Nucleation
1. Introduction The superior physical properties of silicon carbide (SiC) as compared with silicon make it a perspective material for high power and hightemperature electronics [1]. This semiconductor possesses extremely high thermal conductivity and high breakdown voltage. The wide bandgap of SiC provides a low leakage current of the p–n junction
*
Corresponding author. Tel.: +7-812-554-45-70; fax: +7-812326-61-94. E-mail address: [email protected] (A.N. Vorob’ev).
even at high-temperatures. These and other advantageous properties allow us to hope that the industry could replace high-power Si transistors, thyristors and rectifiers used in high-voltage applications by analogous SiC devices in the future. Silicon carbide technology has made tremendous strides in the last few years, with a variety of encouraging device and circuit demonstrations. Nevertheless, there are a number of problems to be resolved before SiC devices become of any substantial commercial value. In particular, highvoltage bipolar devices require thick epitaxial layers with low doping concentration and, hence, a high carrier lifetime. For getting a reasonable yield
0927-0256/02/$ - see front matter 2002 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 7 - 0 2 5 6 ( 0 2 ) 0 0 2 2 0 - 3
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
of these devices, it is necessary to increase the growth rate of epitaxial layers, to achieve a good surface morphology and to provide a high uniformity and controlability of the doping level. Two latter factors are found to depend significantly on the Si/C ratio in the gas mixture supplied onto the growth surface. Chemical vapor deposition (CVD) is widely used for producing SiC films. This growth method turned out to be suitable for meeting the requirements mentioned above. At the same time, there is ongoing competition of different epitaxial setups in reaching a point which can be called production suited and there has been no final decision as to which setup is more advantageous. At present, four main processes are available: • • • •
cold-wall horizontal CVD cold-wall vertical low pressure CVD hot-wall horizontal CVD high-temperature CVD (HTCVD)
SiC growth in reactors of the first two types is carried out at 1400–1500 C and the growth rate for epitaxial layers is of the order of 2–3 lm/h [2]. This means that an active layer of 50 lm takes about 18–25 h to grow. The growth process in reactors of the third and fourth types differs considerably from that in standard CVD. First, there is significant sublimation of the graphite walls occurring at higher growth temperatures (1800– 2300 C) [3]. Second, the growth rates obtained by using hot-wall horizontal reactor CVD, especially HTCVD, vary from a few tens of a micrometer per hour to 0.5 mm/h [4]. The present paper focuses on the numerical study of SiC growth in cold-wall vertical low pressure CVD. There are a few related problems in the performance of this technology. We will emphasize two of them. First, since the carrier gas in reactors of this type goes downward through the inlet at the top to the outlet at the bottom, the flow can become unstable due to natural convection generated by the vertical temperature gradient. As a result, the radial uniformity of the SiC epilayer is violated, and undesirable impurities from the reactor side walls are incorporated into the crystal. Second, spontaneous condensation of silicon
521
vapor is possible under certain operating conditions. This is supported by the fact that an irradiating layer is observed above the substrate in a cold-wall vertical reactor [2]. The gas-phase nucleation depends on the silane flow rate and results in losses of silicon and in a considerable decrease in the growth rate. Besides, the gas-to-particle conversion of silicon influences the Si/C ratio and, consequently, affects the properties of the growing SiC epilayer. Normally, growth conditions are chosen in such a way as to avoid both the vortex formation and the gas-phase nucleation. The aim of this paper is to get an insight into the complex mechanism of SiC growth in a vertical cold-wall reactor and to reveal several important features of the process. Particular attention is given to the growth rate dependence on precursor flow rates and to the effect of silicon vapor condensation in the gas. 2. Model 2.1. Mass transport The gas flow in a cold-wall reactor is essentially subsonic. Under these conditions, the Navier– Stokes equations for low Mach number flow are found to be most suitable for description of gas dynamics and reactive species transport. By neglecting the compressibility effects, steady-state conservation equations modeling a Ns -species mixture motion may be presented in the following vector form [5]: Continuity r qV ¼ 0
ð1Þ
Momentum r qVV þ pI^ s^ ðq q0 Þg ¼ 0
ð2Þ
Energy r ðqVh þ qÞ ¼ 0
ð3Þ
Species r ðqVcs þ Js Þ ¼ W_ s dsc W_ c ; s ¼ 1; 2; . . . ; Ns 1
ð4Þ
522
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
State Ns X cs qRT ¼ P0 ¼ const; Ms s¼1
ð5Þ
where q is the mixture density, V is the gas mixture velocity, p is the dynamic or gauge pressure, h is the mixture enthalpy, cs is the mass fraction of the s-th species, s^ is the viscous stress tensor, q is the heat flux, Js is the mass diffusion flux of the s-th species, W_ s is the generation/consumption rate of the s-th species due to homogeneous chemical reactions, W_ c is the consumption/generation rate of condensible species due to homogeneous nucleation/evaporation, dsc is the Kronecker symbol, P0 is the static pressure, R is the gas constant, T is the gas mixture temperature, Ms is the molar weight of the s-th species, I^ is the unity tensor, q0 is the gas mixture density at the reactor inlet. Eqs. (1)–(5) are closed by the well-known constitutive relations for the molecular transfer fluxes s^, q and Js [5]. 2.2. Homogeneous chemistry Normally, silane (SiH4 ) and propane (C3 H8 ) are used as precursors highly diluted by hydrogen
carrier gas. The temperature strongly influences the pathways of thermal decomposition of these precursors, producing different gaseous products involved in SiC growth. There is a large body of information on SiH4 and C3 H8 chemical transformation in the literature (e.g., see [6,7]). The most detailed data on the gas-phase reaction rates are collected in [7]. In this paper, pyrolysis of the reactants considered is simulated by 83 elementary reactions and by about 30 decomposition products. It should be noted that most of the species have very low concentrations and do not influence Si and C atom transport to the growing surface. Besides, it is impossible to use all of these reactions in 2D and 3D computations because of the computer resource limitation. So, we carried out a preliminary 1D simulation presetting the axial temperature profiles typical for cold-wall CVD. This allows us to reduce significantly the initial reaction set and to identify nine major species important for accurate prediction of the growth rate. The results obtained for the temperature variation from 300 K at the inlet to 1670–1780 K at the substrate are given in Tables 1 and 2. The reaction rate coefficients krþ and kr appear in the expression for W_ s in Eq. (4):
Table 1 Gas-phase reactions r 1 2 3 4 5 6
Reactions
a 10
H þ H þ M ! H2 þ M C3 H8 ! C2 H5 þ CH3 C2 H5 þ H ! CH3 þ CH3 CH4 þ H ! CH3 þ H2 SiH4 ! SiH2 þ H2 SiH2 ! Si þ H2
9:200 10 1:698 1016 1:000 1011 2:200 101 6:671 1029 1:060 1014
b
c
Units
0:600 0:000 0:000 3:000 4:795 0:880
0:000 4:263 104 0:000 4:396 103 3:188 104 2:261 104
m6 kmol2 s1 s1 m3 kmol1 s1 m3 kmol1 s1 s1 s1
Forward reaction mechanism rate coefficients in form krþ ¼ aT b exp ðc=T Þ, 300 K 6 T 6 1780 K.
Table 2 Gas-phase reactions r 1 2 3 4 5 6
Reactions HþHþM H2 þ M C3 H8 CH3 þ C2 H5 C2 H5 þ H CH3 þ CH3 CH4 þ H CH3 þ H2 SiH4 SiH2 þ H2 SiH2 Si þ H2
a
b 13
2:714 10 3:194 103 4:158 106 1:835 102 2:868 1022 1:045 107
0:284 1:804 0.791 3:485 4:203 0.388
c
Units 4
5:210 10 8:215 102 3:363 103 3:544 103 6:279 103 2:391 103
Reverse reaction mechanism rate coefficients in form kr ¼ aT b expðc=T Þ, 300 K 6 T 6 1780 K.
m3 kmol1 s1 m3 kmol1 s1 m3 kmol1 s1 m3 kmol1 s1 m3 kmol 1 s1 m3 kmol1 s1
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
W_ s ¼ Ms
Nhr X
þ m sr msr
r¼1
kr
Ns Y qk M k k¼1
m Ns Y qk kr Mk k¼1 !
þ
krþ
mkr
;
ð6Þ
where Nhr is the total number of homogeneous reactions, mþ kr and mkr are the stoichiometric coefficients for the r-th reaction, and qk ¼ qck is the density of the k-th species. Thus, for simulating cold-wall SiC CVD, six gas-phase reactions and the following minimum set of species are chosen: H2 , H, C3 H8 , C2 H5 , CH4 , CH3 , SiH4 , SiH2 and Si. 2.3. Heterogeneous chemistry Generally, modeling of the chemical interaction between gaseous reactants and a solid surface involves a set of heterogeneous reactions. Surface reactions are characterized by complicated reaction mechanisms including a number of steps, and each surface reaction depends on the partial pressures of the gaseous species, the rate constants of each individual reaction step as a function of temperature and surface properties. However, there is usually little or no detailed information on the individual reaction step and the rate constants. Therefore, one of the steps is frequently assumed to be rate-limiting whereas all the other steps are taken to be in equilibrium (e.g., see [8]). With these assumptions, an expression for the overall reaction rate can be derived. A comparison of the predicted and measured growth rates is needed to show the validity of the assumptions and to evaluate unknown constants. In practice, this approach provides the availability of free parameters which should be adjusted to reproduce a definite set of experimental data. But the adjusted parameters are often found to be entirely unusable for the prediction of growth under different conditions. This is a crucial problem in modeling the heterogeneous chemistry for CVD. An alternative approach is based on the idea that the transport of species from the gas phase to the growing epitaxial layer is the stage limiting growth rate. The approach is feasible provided
523
that the characteristic time of all surface processes is much smaller than that of the gas-phase diffusion. In this case, near-equilibrium between the gaseous species and the crystal is established at the growth surface. Since the species transport in a cold-wall reactor is very slow and the growth temperature is high enough (1450 C), we assume that the collision and subsequent reaction of the gas molecules with the heated substrate surface are limited by mass transport. In this consideration, all the heterogeneous reactions run under nearly equilibrium conditions. In accordance with the Gibbs phase rule, we can write Nsr mass action law equations for all species involved in the surface reactions: Ns Y w mir xi ¼ Kr ;
r ¼ 1; 2; . . . ; Nsr
ð7Þ
i¼1
Here Nsr ¼ Ns þ Ncon Ne , Ncon and Ne are the number of condensed phases formed on the growth surface and the total P number of elements, respecNs tively, xwi ¼ ðqwi =Mi Þ= i¼1 qwi =Mi ¼ pi =P0 , qwi and pi are the molar fraction, the density and the partial pressure of the i-th species on the surface, respectively; mir is the stoichiometric coefficient of the i-th species in the r-th surface reaction, Kr is the equilibrium constant of the r-th surface reaction, determined from the thermodynamic properties of the relevant species. The total species flux onto the surface includes the convective ðqwi V w Þ and diffusion ðJiw Þ terms and is defined as: qwi V w þ Jiw ¼ Mi
Nsr X
mir w_ r ;
i ¼ 1; 2; . . . ; Ns ;
ð8Þ
r¼1
where V w is the Stephan velocity, w_ r is the rate of the r-th surface reaction. To close the set of Eqs. (7) and (8), the normalizing condition should be added: Ns X i¼1
xwi ¼
Ns X
pi =P0 ¼ 1
ð9Þ
i¼1
The set of non-linear Eqs. (7)–(9) must be solved with respect to ðNs þ Nsr þ 1Þ-variables on the solid surface: xwi , w_ r and V w .
524
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
The deposition rate is then found as: Vg ¼
Ns X
Mi
i¼1
Nsr X
mir w_ r =qsurf ;
ð10Þ
r¼1
where qsurf is the surface material density. The solution to the problem (7)–(9) is not influenced by the choice of reactions associated with the mass action laws in (7), provided they are independent. For instance, to describe SiC deposition on a SiC substrate, the following independent reactions can be chosen: CH4 þ Si $ 2H2 þ SiCðsolidÞ 2H $ H2 C3 H8 $ CH3 þ C2 H5 CH4 þ H $ CH3 þ H2 C2 H5 þ H $ 2CH3 SiH4 $ SiH2 þ H2 SiH2 $ Si þ H2 2.4. Gas-phase nucleation As mentioned above, a key feature of the process considered is homogeneous nucleation of silicon. There is a large amount of information on the subject. At present, there are two principal approaches to simulate a gas flow coupled with condensation. One is well-known and is based on the classical theory suggested by Volmer [9], and developed by Becker and D€ oring [10], Zeldovich [11] and Frenkel [12]. Despite the obvious drawbacks and a number of theoretical attacks, the classical theory is adequate in itself to predict the nucleation rate consistent with experimental data for many important processes under quasi-stationary conditions. However, nucleation is not the only stage of condensation. Nuclei can grow in a supersaturated vapor and, generally, their growth rate varies with the cluster size. Analyzing experimental data, Hill [13] concluded that the growth rate is independent of the radius for small clusters, and it is sufficient to employ equations for the first
three moments of the size-distribution function for modeling the drop evolution after the nucleation. This consideration significantly reduces the set of governing equations. The other approach using the discrete-sectional model has been lately applied to the problem discussed (e.g., see [14]). The method considers two parts of the simulated particle-size range: discrete and sectional. The discrete part handles the detailed interactions of monomers and small clusters up to a predetermined maximum cluster size. The upper-size range of the distribution is approximated by a small number of sections. However, the discrete-sectional model has been applied mostly to 1D flow of binary mixture because of the limited computational power [15,16]. By now, this approach has been hardly feasible to 2D and 3D flow of a multi-component gas mixture coupled with condensation of a species. Estimations show that under the operating conditions in the vertical cold-wall reactor, the drop size is a few orders of magnitude less than the mean free path of molecules in the vapor. Therefore, clusters can be considered as a pseudo-gas of heavy particles with the size-independent growth rate. Since the gas flow is slow, a steady-state sizedistribution of nuclei becomes achievable. In addition, our preliminary estimates and further 2D computations within the employed model show that such important mechanism as coagulation can be neglected under operating conditions typical for CVD of SiC in the considered reactor. The numerical concentration of clusters here does not exceed 1016 m3 , whereas their size is small (aver). Then, the estimation of age radius of 30–50 A mean free path of the particles gives 0.2 m and larger, while the characteristic size of the reactor is 0.1 m. So, the probability of collision between two neighboring clusters to generate a larger cluster is negligible. This results primarily from the fact that the precursors are highly diluted by the carrier gas, and the process is operated under low total pressures (50 Torr). In this situation, the approach suggested in [13] is suited to simulate the gas-phase condensation of silicon in a non-isothermal flow of a multi-component gas mixture. An approximate kinetic equation can be derived from either the Focker–Plank equation or
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
525
following simple considerations. Let us introduce the size distribution function (f) which is normalized as Z 1 4 qc ¼ pq0Si r3 f ðrÞ dr; ð11Þ 3 r
1 ðT Þ is supersaturation. Within the where s ¼ pSi =pSi assumptions mentioned above, the rate of mass exchange between the clusters and the gas phase can be simulated by using the Hertz–Knudsen equation; hence
where r is the particle radius, qc is the density of the particle pseudo-gas, q0Si is the liquid density. By neglecting the diffusion, this function can vary due to three main mechanisms: the particle motion in the gas flow, the gas-phase nucleation and the condensation/evaporation of gas molecules on/from the cluster surface. The balance of these mechanisms is expressed as o r ð Vc f Þ þ ð12Þ r_ f ¼ Idðr r Þ; or
r_ ¼
where Vc is the transport velocity of the particle pseudo-gas, r is the critical nucleus radius, I is the formation rate of critical-sized nuclei, r_ is the cluster-growth rate after nucleation. The process of steady-state nucleation is considered to be a flow of constant net current of nuclei through a steady distribution of nucleus sizes. In the other words, the net rate at which g-sized nuclei grow to ðg þ 1Þ-sized nuclei is equal to the net rate, at which g-sized nuclei grow from ðg 1Þ-sized ones, so that the concentration of nuclei of size g is constant. With this simplifying concept, the nucleation rate (I) is given by the Becker–D€ oring theory with Bernard’s correction factor [10]: sffiffiffiffiffiffiffiffiffiffiffiffiffiffi q2Si 2rSi NA3 4 rSi NA I¼a 0 exp p 3 RT qSi pMSi3 ð13Þ
r2 ð1 3g2=3 þ 2g1 Þ ; where a is the condensation coefficient, rSi is the surface tension coefficient of liquid silicon, NA is the Avogadro number, MSi is the silicon molar weight. The critical nucleus radius, r , and the number of molecules contained in the critical nucleus, g , are calculated from a maximum of the net change in free energy necessary for reducing the vapor pressure to the flat-film saturation value and for creating a liquid drop [13] r ¼
2rSi MSi ; q0Si RT ln s
4 NA q0Si 3 r ; g ¼ p 3 MSi
ð14Þ
1 a pSi pSi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi; 0 qSi 2pRT =MSi
ð15Þ
1 where pSi is the partial pressure of silicon and pSi is the flat-film saturation pressure. With account of Eqs. (13)–(15), the evolution equations for the first three moments of the sizedistribution function, Xn , are derived by multiplying Eq. (12) by rn and by integrating over r from r to 1
r ðVc Xn Þ ¼ rn I þ n_rXn1 ; where Xn is defined as Z 1 rn f ðrÞ dr Xn ¼
n ¼ 0; 1; 2;
ð16Þ
ð17Þ
r
Since a large temperature gradient is found near the heated substrate in a vertical cold-wall CVD reactor, the thermophoretic force can make the particle path deviate significantly from the actual gas flow streamlines [17–19]. The basic mechanism of thermophoresis is that a suspended particle in a temperature gradient experiences a net force in the direction of decreasing temperature. This force arises because molecules impacting the particle on opposite sides have different average velocities due to the temperature gradient. The thermophoretic velocity is derived theoretically and accounted for as [18] l Vc ¼ V KT rT ; ð18Þ qT where l is the gas mixture viscosity, KT is the thermophoretic coefficient. The value of this constant depends on the molecular flow regime. A measure of this regime is the Knudsen number––the ratio of the mean free path length of the gas molecule to the particle radius. For large Knudsen’s numbers, KT of 0.55 is predicted by the theory describing the motion of a single particle in the temperature gradient [17]. As follows from literature (e.g., see [18]), this value is chosen as a basic one and varied from 0.0 to 0.9 to bring the
526
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
computations into agreement with available experiment. In this study, the value of KT ¼ 0:2 appears most appropriate for the conditions encountered in the process considered. Transport of clusters is described by their continuity equation r ðqc Vc Þ ¼ W_ c ;
ð19Þ
where 4 W_ c ¼ pq0Si ðr3 I þ 3_rX2 Þ 3
ð20Þ
Expression (20) is introduced into Eq. (4) to provide the consumption or generation of a condensible species due to homogeneous nucleation or evaporation of clusters. Since the precursors SiH4 and C3 H8 are highly diluted by H2 under typical operating conditions, the corresponding terms associated with spontaneous condensation were neglected in Eqs. (1)–(3).
3. Results and discussion 3.1. Gas dynamics The CVD system elaborated by EMCORE Corp. (New Jersey) has been close to treat SiC growth. Precursors highly diluted by hydrogen carrier gas can be supplied into the reactor using either the whole area of the reactor inlet or a hole/ holes of a smaller diameter. The process is operated under the following conditions [2]: the total pressure of 50–300 Torr, the hydrogen flow of 15–30 slm, the susceptor rotation rate of 750 rpm, and the wafer temperature of 1725 K. The flow in this type of reactors is said to occur either in the forced convection mode, where streamlines are curved but non-circulating above the susceptor, or in a mixed convection mode, with undesirable vortices caused by the heated susceptor. The mode transition can be roughly estimated by using the well-known dimensionless flow parameters, namely, the Reynolds (Re) and Grashof (Gr) numbers [20]. For small values of the Gr=Re3=2 ratio, there is no flow recirculation above the rotating susceptor. However, the inlet gas velocity is an additional parameter that enters the
problem through the specification of the flow stability boundary. Obviously, some computations are needed to refine the boundary. An experimental and numerical analysis of the flow pattern as a function of the total pressure was carried out in [2]. It was shown that the mode transition occurs in the narrow range of 200–250 Torr for the total flow rate of 30 slm and the rotation rate of 750 rpm. With increasing total pressure, the flow becomes unstable and essentially 3D. Since most experimental data on the growth rate reported in [21] concern the total flow rate of 15 slm and the rotation rate of 750 rpm, we have found numerically the transition between the forced and mixed convection modes with the pressure variation. Our computations show that the flow remains to be stable up to 75 Torr. Fig. 1 gives an example of the velocity vector and temperature distribution at the total pressure of 50 Torr when the top inlet area is wholly used for supplying the gas into the reactor. One can see that there is a steep temperature gradient near the wafer heated to 1725 K and the flow is non-circulating everywhere, except the zone close to the reactor outlet, which is inessential to the epilayer growth. The typical mass concentration distribution of precursors and major species are shown in Figs. 2
Fig. 1. Stable flow pattern and temperature distribution.
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
527
central hole at the inlet. In this case, the inlet parameters were recalculated in such a way to keep the same mass fluxes of precursors and hydrogen supplied into the reactor. For this purpose, we varied the inlet density at the constant velocity over whole inlet area. As an illustration, let us consider the following typical operating conditions: • Total pressure ¼ 50 Torr • Flow rates: H2 ¼ 15 slm, C3 H8 ¼ 6:6 sccm, SiH4 ¼ 17 sccm • Inlet temperature ¼ 300 K Fig. 2. Mass concentration distribution of propane and methane.
and 3. It follows from the figures that propane and silane are completely decomposed near the heated substrate, producing CH4 and SiH2 . These species and atomic silicon provide SiC growth on a wafer placed at the center of the holder top. Fig. 4 illustrates the formation of a large vortex above the substrate at the total pressure of 100 Torr. Elevation of the total pressure to 200 Torr does not change the general features of the flow pattern, while the species mixing is intensified. The same analysis was carried out when the precursors go into the reactor only through the
Fig. 3. Mass concentration distribution of silane and silylene.
Calculated inlet parameters: 1. Whole inlet area is used: Inlet velocity ¼ 0:121148 m/s Mass fractions: H2 ¼ 0:972939, C3 H8 ¼ 0:009418, SiH4 ¼ 0:017643 Mixture density ¼ 0:005414 kg/m3 2. Central hole is used: Inlet velocity ¼ 0:121148 m/s At outer ring: Mass fractions: H2 ¼ 1:0, C3 H8 ¼ 0:0, SiH4 ¼ 0:0 Mixture density ¼ 0:005275 kg/m3 At central hole: Mass fractions: H2 ¼ 0:668616, C3 H8 ¼ 0:115333, SiH4 ¼ 0:216051 Mixture density ¼ 0:007675 kg/m3 Thus, the mixture density can rise then almost 1.5 times at the central hole. Note that the transition of the gas flow from the laminar regime to the unstable one is somewhat delayed. When the radius ratio of the central hole and the reactor inlet is 1:4, the flow remains to be non-circulating to 150 Torr at the propane and silane flow rates of 6.6 and 17 sccm, respectively. The flow stabilization is due to the fact that the increase in the gas mixture density at the reactor center tends to suppress vortices generated by the steep temperature gradient near the wafer. Analyzing the preliminary results, we limited the range of investigations and mostly concentrated on the mass transport under typical gas dynamic conditions: the total pressure of 50 Torr,
528
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
Fig. 5. Growth rate vs propane flow rate at SiH4 ¼ 17 sccm. The hole at the inlet is used for supplying SiH4 and C3 H8 . Full circles are for the experiment. Fig. 4. Unstable flow pattern and temperature distribution.
the susceptor rotation rate of 750 rpm, the wafer temperature of 1725 K, and the hydrogen flow rate of 15 slm. The propane and silane flow rates are allowed to vary in a wide range. 3.2. Verification of the model The model was verified by comparing the computations with the available experimental data on the SiC growth rate as a function of the input propane and silane flow rates. There is no detailed information on what inlet flange was used for supplying precursors into the growth chamber at the operating parameters mentioned above. For this reason, we carried out two sets of computations corresponding to the precursor supply through the central hole and the whole inlet area. The growth rate variation with the propane flow rate is plotted in Fig. 5. The precursors are supplied through the central hole of a radius four times smaller than that of the reactor inlet and the silane flow rate is kept at 17 sccm. As seen from the figure, the computations fit well the experiment. The sloping portion of the curve corresponds to regimes where silane is in excess and the growth rate is controlled by the propane flow rate. The growth rate begins to saturate at the input
carbon-to-silicon ratio C/Si 0.6 and can no longer be raised by elevating the propane flow rate. Normally, the saturation is due the transition from the propane-limited mode to the silane-limited and is observed at C/Si 1. We speculate that the earlier saturation in Fig. 5 is related to homogeneous nucleation in the gas. The main species assumed to be responsible for this phenomenon is silicon produced by the decomposition of silylene SiH2 Si þ H2 (see the reaction in Table 1). The silicon gasto-particle conversion depletes the vapor in the reactive species needed for heterogeneous surface reactions of the film growth. Support for the speculation comes from experimental evidence on the irradiating layer observed through the reactor viewport and attributed to black-body radiation scattering by silicon clusters [2]. The computed distribution of silicon cluster density over the substrate is shown in Fig. 6. It reproduces well the shape and the location of the irradiating layer in [21]. The result obtained suggests the following mechanism of the gas-phase depletion. The clusters affected by the thermophoretic force are accumulated in a thin layer located over the susceptor and cannot reach the growing surface. The gas flow takes the clusters away in the radial
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
529
Fig. 7. Growth rate vs propane flow rate at SiH4 ¼ 17 sccm. The whole inlet is used for supplying SiH4 and C3 H8 . Full circles are for the experiment. Fig. 6. Typical Si cluster density distribution.
direction (see Fig. 1), and, hence, some of the Si atoms stored in the clusters are lost to the SiC growth on the wafer. Note that the model neglecting silicon vapor condensation predicts saturation at C/Si P 1.0 and significantly higher growth rates [21]. A similar growth rate dependence corresponding to the use of the whole inlet section for supplying the precursors is presented in Fig. 7. A comparison of the curves in Figs. 5 and 7 shows that a higher propane flow rate of 5 sccm (C/ Si 0.8) is required to get the maximum growth rate in the latter case. This is due to the fact that some amount of precursors goes between the water-cooled wall and the susceptor as far as the outlet and does not reach the wafer. Thus, a careful choice of the inlet flange design can provide a higher precursor utilization and, hence, a higher growth rate. The growth rate is shown in Fig. 8 as a function of the silane flow rate at the constant propane flow rate of 6.6 sccm. The solid and dashed curves correspond to the supply of precursors through the central hole and the whole inlet section, respectively. As seen from the theoretical plots, a deviation of the growth rate dependence from the
Fig. 8. Growth rate vs silane flow rate at C3 H8 ¼ 6:6 sccm. Precursors are supplied through the hole (solid) and whole inlet (dashed). Full circles are for the experiment.
linear that occurs at the silane flow rate of 10 sccm. This is due to the onset of nucleation supported by the sharp increase of the cluster density over the substrate. In addition, the computations predict a decrease of the growth rate at high SiH4 flow rates. This behavior is associated with the greater Si vapor supersaturation and, therefore, with a higher gas-to-particle conversion of silicon.
530
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
Note that the theoretical curves shown in Fig. 8 are in poor agreement with the experimental data compared to the case of propane flow rate variation. In particular, the predicted growth rate begins to decrease at the SiH4 flow rate greater than 15 sccm, while the measured growth rate does not exhibit such a behavior. In our opinion, this discrepancy is related to an extremely high sensitivity of the growth rate to the detailed pathway of silicon nucleation. Indeed, the growth rate is directly affected by the material losses due to the gasto-particle conversion. The cluster nucleation, in turn, depends very non-linearly on the silicon vapor supersaturation. Therefore, small deviation of the real growth conditions from those set in computations may result in the considerable discrepancy between the theoretical predictions and experimental observations. In our case, the effect of unaccounted details of the growth process is found to be insignificant for the accurate prediction of the growth rate behavior at the propane flow rate variation. However, this effect is remarkable when the SiH4 flow rate is varied. On the other hand, the comparison of the computational and experimental data shows that the model suggested catches the major physical mechanisms underlying the CVD of SiC, as a whole. Therefore, we apply further the model to the computational analysis of the process by varying the precursor flow rates in a wider range. 3.3. Effect of precursor flow rate variation on the growth rate For getting a better insight into the growth rate behavior, we carried out computations, varying the precursor flow rate in the following ranges: C3 H8 ¼ 1:5–8:0 sccm, SiH4 ¼ 4:0–60:0 sccm. The growth rate variation with the propane flow rate is shown in Fig. 9(a) and (b) at several constant SiH4 flow rates. As seen from the figures, the growth rate initially rises at low propane flow rates and then reaches a maximum. This level increases in the range of the silane flow rate from 4.0 to 17 sccm. The higher is the silane flow rate, the larger is the propane flow rate at which the onset of the growth rate saturation occurs (see Fig. 9(a)). The maximum level and the saturation onset depend
Fig. 9. Growth rate vs propane flow rate at several constant SiH4 flow rates: (a) 4–17 sccm, (b) 17–60 sccm. The whole inlet is used for supplying SiH4 and C3 H8 .
inversely on the silane flow rate at SiH4 > 17 sccm (see Fig. 9(b)). This is due to the fact that the depletion of the gas mixture in silicon atoms involved in the growth process is enhanced by the silicon gas-to-particle conversion at higher silane flow rates. The growth rate is plotted in Fig. 10 as a function of the silane flow rate at several constant propane flow rates. The plateaus demonstrate the growth rate saturation at a low propane flow rate from 1.5 to 3.3 sccm where C3 H8 controls the process. At C3 H8 > 3:3 sccm, the plateaus disap-
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
Fig. 10. Growth rate vs silane flow rate at several constant C3 H8 flow rates. The whole inlet is used for supplying SiH4 and C3 H8 .
pear, and the growth rate dependence on silane flow rate peaks at SiH4 17 sccm. The subsequent decrease of the growth rate in both cases results from the gas-phase nucleation of silicon. At C3 H8 < 3:3 sccm, this occurs differently after the gas mixture becomes fairly depleted in silicon and the carbon-containing species becomes in excess. Note that an additional increase of the propane flow rate from C3 H8 ¼ 5:5 sccm and higher does not practically influence the maximum growth rate. For visualizing more clearly the trends discussed above, a 3D growth rate surface is plotted in Fig. 11 as a function of C3 H8 and SiH4 flow rates. It follows from the figure that the growth rate cannot exceed a critical value. The limitation is caused by the drastic silicon condensation in the vapor at high silane flow rates and means that it is impossible to provide a desirable increase in the growth rate via the precursor flow rate variation. Thus, the high flow rates of silane result in a greater loss of silicon atoms involved in SiC growth, and in spite of the inlet Si/C ratio being much more unity, the growth rate is then controlled by SiH4 . 3.4. Gas-phase nucleation effect on the Si/C ratio on the growth surface It is well-known that Si/C ratio plays a crucial role for many properties of SiC crystals. In par-
531
Fig. 11. 3D surface of growth rate vs SiH4 and C3 H8 flow rate.
ticular, the doping by various impurities, the concentration of carbon vacancies, surface morphology including step bunching, etc. are affected by the surface Si/C ratio which may differ significantly from the Si/C ratio at the reactor inlet. Since the Si/C ratio controls the adsorption layer composition on the growth surface, it should be directly related to the collision frequencies of molecules containing silicon and carbon atoms, respectively. So, we calculate the Si/C ratio as Pn Si fi bi pi ; ð21Þ Si=C ¼ Pi¼1 n C i¼1 fi bi pi where fiSi and fiC are the number of Si and C atoms in a molecule of the i-th species, respectively, 1=2 bi ¼ ð2pmi kT Þ and mi are the Hertz–Knudsen factor and the molecular mass of the i-th species, respectively, pi is the partial pressure of the i-th species. Normally, crystal growth is operated under excess of either silicon or carbon species. The growth rate is then controlled by supplying the species with a lower inlet concentration. Near the growth surface, the concentration of such species additionally decreases due to a more efficient incorporation into the crystal, as compared to the species in excess. For this reason, the actual Si/C ratio on the growth surface can differ significantly from that at the reactor inlet. Obviously, gas-phase processes like homogeneous nucleation
532
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
Fig. 12. Si/C ratio on the wafer and growth rate vs SiH4 flow rate at C3 H8 ¼ 6:6 sccm. No gas-phase nucleation is accounted for.
not contributing to the crystallization may alter additionally the Si/C ratio on the growth surface. The experimental [2] and predicted growth rates and Si/C ratio are plotted in Figs. 12 and 13, respectively, as a function of silane at the constant propane flow rate of 6.6 sccm with and without accounting for nucleation. Other operating parameters had their typical values: the hydrogen carrier gas flow rate of 15 slm, the total pressure of 50 Torr, and the growth temperature of 1450 C. One can see saturation of the deposition rate with
Fig. 13. Si/C ratio on the wafer and growth rate vs SiH4 flow rate at C3 H8 ¼ 6:6 sccm. Gas-phase nucleation is accounted for.
silane flow rate increase, usually observed in the experiments. Model results without the account of the gas-phase nucleation give a clear saturation of the dependence at the point where the inlet Si/C ratio becomes close to unity (Fig. 12). The predicted growth rates are higher than the experimental values. When the gas-phase nucleation is taken into account (Fig. 13), the consistence of the computed and experimental growth rates is better. What is important is that the growth rate saturation starts earlier. Note that the local values of the Si/C ratio over the wafer are also less than unity under the conditions with the inlet Si/C ratio larger than unity. The local Si/C ratio is predicted to decrease with silane supply at high inlet Si/C ratios. This effect results from the fairly high silane flow rate. To get a general pattern of the Si/C ratio near the growth surface, we made computations varying the precursor flow rates in the following ranges: C3 H8 ¼ 1:5–8:0 sccm, SiH4 ¼ 4:0–60:0 sccm. The Si/C behavior is illustrated in Fig. 14. Note that the proportion between silicon and carbon atoms here is an abrupt function and is either much less or much more than unity. This results from the accepted assumption of the near-equilibrium growth conditions at temperatures as high as 1450 C. The linear portion of the curve for Si=C ¼ 1 separates the SiC growth regimes when the silane flow rate is insufficiently high for the onset of nucleation. The Si/C ratio is then con-
Fig. 14. Si/C ratio near the growth surface vs precursor flow rates.
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
trolled completely by that at the inlet. The silicon gas-to-particle conversion manifests itself as the curve slope. The area of Si=C < 1 becomes larger than the area of Si=C > 1 due to the drastic depletion of the vapor in silicon near the growth surface. For this reason, the higher silane flow rate is unable to raise the Si/C ratio at SiH4 > 17 sccm. 3.5. Secondary phase formation The formation of parasitic (secondary) phases on the SiC surface (solid graphite and either liquid or solid silicon) is an undesirable effect leading to degradation of the material properties. Here we suggest a special criterion to predict whether a condensed phase is formed on the growth surface. We assume a secondary phase arisen to be in equilibrium both with the SiC crystal and with the common gas phase. Then, in addition to the heterogeneous chemical reaction producing crystalline SiC, we consider the reaction giving rise to the secondary phase. Let W_ SiC be the rate of a heterogeneous reaction producing a SiC crystal and W_ SP the rate of the secondary phase generation on the surface. These reaction rates can be calculated using the mass balance equations for the growth surface coupled with the mass action law equations corresponding to heterogeneous equilibrium vapor-crystal-secondary phase equilibrium. Depending on the rates of these reactions, the following cases are possible: 1
4. At W_ SP < 0 and W_ SiC < 0, SiC crystal evaporates congruently. The boundaries between these cases can be found from the equations W_ SP ¼ 0 and W_ SiC ¼ 0. The analysis of the secondary phase formation was made for the following variation of precursor flow rates: C3 H8 ¼ 0:5–8:0 sccm, and SiH4 ¼ 4:0–60:0 sccm. Fig. 15 shows by solid lines the boundaries between the zones where the secondary phase formation and/or growth of the SiC single crystal are predictable. For comparison, the similar boundaries computed with the neglect of the silicon condensation in the gas are shown by dashed lines. It is seen that without nucleation there would be a very narrow growth window where single-crystal SiC could be obtained. The silicon nucleation alters the Si/C ratio over the substrate and, hence, broadens the growth window significantly. It should be noted that the shape of the growth window is rather complex and is controlled totally by the silicon clustering. On the other hand, the growth rate is predicted to peak at SiH4 17 sccm (see Fig. 10 or 11). So, according to Fig. 15, the propane flow rate is allowed to vary in the range 3:3–5:5 sccm to avoid the parasitic phase formation on the growing surface. It can be also concluded from the comparison of the growth window with that computed with neglect of silicon clustering that homogeneous nucleation indirectly influences SiC crystal quality, and, undoubtedly,
1. At W_ SP < 0 and W_ SiC > 0, a single SiC crystal grows on the surface without secondary phase formation. 2. At W_ SP > 0 and W_ SiC < 0, the growth surface becomes completely covered by a secondary phase. 3. At W_ SP > 0 and W_ SiC > 0, simultaneous deposition of both SiC and secondary phase occurs on the growth surface. Then the surface fraction occupied by the secondary phase can be esti mated from the ratio as W_ SP = W_ SP þ W_ SiC .
1
In more detail, the criterion for a secondary phase formation, its verification and examples of application to various growth systems will be published in a separate paper.
533
Fig. 15. Secondary phase map on the growth surface.
534
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
should be accounted for when optimizing the SiC CVD in a cold-wall vertical reactor.
4. Summary An advanced model of heat and mass transport coupled with spontaneous condensation of silicon in the gas is suggested to simulate SiC CVD in a cold-wall rotating-disc vertical reactor. The comparison of computed growth rates with available experimental data shows a reasonable agreement under typical operating conditions. It is found that homogeneous nucleation considerably influences the CVD process. In particular, it is impossible to provide an increase in the growth rate over some maximum value by varying the precursor flow rates due to a drastic silicon gas-phase condensation at high silane flow rates. The analysis shows that the silicon-to-carbon ratio near the growth surface can radically differ from that at the reactor inlet. In turn, this changes the proportion between the silicon and carbon atoms on the wafer and increases the range of silane and propane variation where the SiC single epilayer is expected to grow without secondary phase formation or morphology degradation. The model allows us to predict the optimal precursor flow rates, ensuring the maximum growth rate at relatively broad growth window.
References [1] A. Itoh, H. Matsunami, IEEE Electron Dev. Lett. 16 (1995) 280.
[2] R. Rupp, Yu.N. Makarov, H. Behner, A. Wiedenhofer, Phys. Stat. Sol. (b) 202 (1997) 281. [3] O. Kordina, C. Hallin, A. Ellison, A.S. Bakin, I.G. Ivanov, A. Henry, R. Yakimova, M. Touminen, A. Vehanen, E. Janzen, Appl. Phys. Lett. 69 (10) (1996). [4] O. Kordina, C. Hallin, A. Henry, J.P. Bergman, I. Ivanov, A. Ellison, N.T. Son, E. Jenzen, Phys. Stat. Sol. (b) 202 (1997) 321. [5] Yu.E. Egorov, A.I. Zhmakin, Comput. Mat. Sci. 11 (1998) 204–220. [6] C.D. Stinespring, J.C. Wormhoudt, J. Crystal Growth 87 (1988) 481–493. [7] M.D. Allendorf, R.J. Kee, J. Electrochem. Soc. 138 (3) (1991) 841–852. [8] M.E. Coltrin, R.J. Kee, G.H. Evans, J. Electrochem. Soc. 136 (3) (1989) 819–829. € bersattigten [9] M. Volmer, A. Weber, Keimbildung in u Gebilden 119 (1926) 277–301. [10] R. Becker, W. D€ oring, Kinetische Behandlung der Keim€bers€attigten D€ampfen, Annalen der Physik 24 bildung in u (1935) 719–752. [11] J. Zeldovich, J. Exp. Theor. Phys. (U.S.S.R.) 12 (1942) 525. [12] J. Frenkel, Kinetic Theory of Liquids, Dover, New York, 1969. [13] Ph.G. Hill, J. Fluid. Mech. 25 (3) (1965) 593–620. [14] J.J. Wu, R.C. Flagan, J. Colloid Interf. Sci. 123 (2) (1988) 339–352. [15] S. Yu, Y. Yoon, M. M€ uller-Roosen, I.M. Kennedy, Aerosol Sci. Technol. 28 (1998) 185–196. [16] S.-Y. Lu, H.-C. Lin, C.-H. Lin, J. Crystal Growth 200 (1999) 527–542. [17] L. Talbot, R.K. Cheng, R.W. Schefer, D.R. Willis, J. Fluid Mech. 101 (4) (1980) 737–758. [18] K.L. Walker, F.T. Geyling, S.R. Nagel, Thermophoretic deposition of small particles in the modified chemical deposition (MCVD) process, J. Am. Ceram. Soc. 63 (1980) 552–558. [19] D.I. Fotiadis, K.F. Jensen, J. Crystal Growth 102 (1990) 743–761. [20] G. Evans, R. Greif, J. Heat Transfer 109 (1987) 928–935. [21] A.N. Vorob’ev, Yu.E. Egorov, Yu.N. Makarov, A.I. Zhmakin, A.O. Galyukov, R. Rupp, Mat. Sci. Eng. B 61–62 (1999) 172–175.
Numerical study of SiC CVD in a vertical cold-wall reactor A.N. Vorob’ev a,*, S.Yu. Karpov a, M.V. Bogdanov a, A.E. Komissarov a, O.V. Bord a, A.I. Zhmakin b, Yu.N. Makarov c b
a Soft-Impact Ltd., P.O. Box 33, 194156 St. Petersburg, Russia A.F. Ioffe Physical Technical Institute, Russian Academy of Sciences, Polytechnicheskaya 26, 194021 St. Petersburg, Russia c STR Inc., P.O. Box 70604, Richmond, VA 23255, USA
Received 7 August 2001; received in revised form 29 January 2002; accepted 1 February 2002
Abstract The conventional heat and mass transport model is extended to describe silicon cluster formation in the gas phase and is employed for a numerical analysis of SiC chemical vapor deposition in a commercial vertical rotating disc reactor. The model is verified by comparing the computed growth rate with available experimental data. The growth rate is studied as a function of precursor flow rates varied in a wide range of values. It is found that the growth rate is limited by the gas mixture depletion in silicon atoms due to homogeneous nucleation. The secondary phase formation on the growing surface is analyzed. The SiC growth window depending on the precursor flow rates is calculated, and a significant influence of the homogeneous nucleation on the window width is shown. The model results predict that the Si/C ratio on the wafer can considerably differ from that at the reactor inlet. 2002 Elsevier Science B.V. All rights reserved. Keywords: Chemical vapor deposition; Cluster; Growth rate; Nucleation
1. Introduction The superior physical properties of silicon carbide (SiC) as compared with silicon make it a perspective material for high power and hightemperature electronics [1]. This semiconductor possesses extremely high thermal conductivity and high breakdown voltage. The wide bandgap of SiC provides a low leakage current of the p–n junction
*
Corresponding author. Tel.: +7-812-554-45-70; fax: +7-812326-61-94. E-mail address: [email protected] (A.N. Vorob’ev).
even at high-temperatures. These and other advantageous properties allow us to hope that the industry could replace high-power Si transistors, thyristors and rectifiers used in high-voltage applications by analogous SiC devices in the future. Silicon carbide technology has made tremendous strides in the last few years, with a variety of encouraging device and circuit demonstrations. Nevertheless, there are a number of problems to be resolved before SiC devices become of any substantial commercial value. In particular, highvoltage bipolar devices require thick epitaxial layers with low doping concentration and, hence, a high carrier lifetime. For getting a reasonable yield
0927-0256/02/$ - see front matter 2002 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 7 - 0 2 5 6 ( 0 2 ) 0 0 2 2 0 - 3
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
of these devices, it is necessary to increase the growth rate of epitaxial layers, to achieve a good surface morphology and to provide a high uniformity and controlability of the doping level. Two latter factors are found to depend significantly on the Si/C ratio in the gas mixture supplied onto the growth surface. Chemical vapor deposition (CVD) is widely used for producing SiC films. This growth method turned out to be suitable for meeting the requirements mentioned above. At the same time, there is ongoing competition of different epitaxial setups in reaching a point which can be called production suited and there has been no final decision as to which setup is more advantageous. At present, four main processes are available: • • • •
cold-wall horizontal CVD cold-wall vertical low pressure CVD hot-wall horizontal CVD high-temperature CVD (HTCVD)
SiC growth in reactors of the first two types is carried out at 1400–1500 C and the growth rate for epitaxial layers is of the order of 2–3 lm/h [2]. This means that an active layer of 50 lm takes about 18–25 h to grow. The growth process in reactors of the third and fourth types differs considerably from that in standard CVD. First, there is significant sublimation of the graphite walls occurring at higher growth temperatures (1800– 2300 C) [3]. Second, the growth rates obtained by using hot-wall horizontal reactor CVD, especially HTCVD, vary from a few tens of a micrometer per hour to 0.5 mm/h [4]. The present paper focuses on the numerical study of SiC growth in cold-wall vertical low pressure CVD. There are a few related problems in the performance of this technology. We will emphasize two of them. First, since the carrier gas in reactors of this type goes downward through the inlet at the top to the outlet at the bottom, the flow can become unstable due to natural convection generated by the vertical temperature gradient. As a result, the radial uniformity of the SiC epilayer is violated, and undesirable impurities from the reactor side walls are incorporated into the crystal. Second, spontaneous condensation of silicon
521
vapor is possible under certain operating conditions. This is supported by the fact that an irradiating layer is observed above the substrate in a cold-wall vertical reactor [2]. The gas-phase nucleation depends on the silane flow rate and results in losses of silicon and in a considerable decrease in the growth rate. Besides, the gas-to-particle conversion of silicon influences the Si/C ratio and, consequently, affects the properties of the growing SiC epilayer. Normally, growth conditions are chosen in such a way as to avoid both the vortex formation and the gas-phase nucleation. The aim of this paper is to get an insight into the complex mechanism of SiC growth in a vertical cold-wall reactor and to reveal several important features of the process. Particular attention is given to the growth rate dependence on precursor flow rates and to the effect of silicon vapor condensation in the gas. 2. Model 2.1. Mass transport The gas flow in a cold-wall reactor is essentially subsonic. Under these conditions, the Navier– Stokes equations for low Mach number flow are found to be most suitable for description of gas dynamics and reactive species transport. By neglecting the compressibility effects, steady-state conservation equations modeling a Ns -species mixture motion may be presented in the following vector form [5]: Continuity r qV ¼ 0
ð1Þ
Momentum r qVV þ pI^ s^ ðq q0 Þg ¼ 0
ð2Þ
Energy r ðqVh þ qÞ ¼ 0
ð3Þ
Species r ðqVcs þ Js Þ ¼ W_ s dsc W_ c ; s ¼ 1; 2; . . . ; Ns 1
ð4Þ
522
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
State Ns X cs qRT ¼ P0 ¼ const; Ms s¼1
ð5Þ
where q is the mixture density, V is the gas mixture velocity, p is the dynamic or gauge pressure, h is the mixture enthalpy, cs is the mass fraction of the s-th species, s^ is the viscous stress tensor, q is the heat flux, Js is the mass diffusion flux of the s-th species, W_ s is the generation/consumption rate of the s-th species due to homogeneous chemical reactions, W_ c is the consumption/generation rate of condensible species due to homogeneous nucleation/evaporation, dsc is the Kronecker symbol, P0 is the static pressure, R is the gas constant, T is the gas mixture temperature, Ms is the molar weight of the s-th species, I^ is the unity tensor, q0 is the gas mixture density at the reactor inlet. Eqs. (1)–(5) are closed by the well-known constitutive relations for the molecular transfer fluxes s^, q and Js [5]. 2.2. Homogeneous chemistry Normally, silane (SiH4 ) and propane (C3 H8 ) are used as precursors highly diluted by hydrogen
carrier gas. The temperature strongly influences the pathways of thermal decomposition of these precursors, producing different gaseous products involved in SiC growth. There is a large body of information on SiH4 and C3 H8 chemical transformation in the literature (e.g., see [6,7]). The most detailed data on the gas-phase reaction rates are collected in [7]. In this paper, pyrolysis of the reactants considered is simulated by 83 elementary reactions and by about 30 decomposition products. It should be noted that most of the species have very low concentrations and do not influence Si and C atom transport to the growing surface. Besides, it is impossible to use all of these reactions in 2D and 3D computations because of the computer resource limitation. So, we carried out a preliminary 1D simulation presetting the axial temperature profiles typical for cold-wall CVD. This allows us to reduce significantly the initial reaction set and to identify nine major species important for accurate prediction of the growth rate. The results obtained for the temperature variation from 300 K at the inlet to 1670–1780 K at the substrate are given in Tables 1 and 2. The reaction rate coefficients krþ and kr appear in the expression for W_ s in Eq. (4):
Table 1 Gas-phase reactions r 1 2 3 4 5 6
Reactions
a 10
H þ H þ M ! H2 þ M C3 H8 ! C2 H5 þ CH3 C2 H5 þ H ! CH3 þ CH3 CH4 þ H ! CH3 þ H2 SiH4 ! SiH2 þ H2 SiH2 ! Si þ H2
9:200 10 1:698 1016 1:000 1011 2:200 101 6:671 1029 1:060 1014
b
c
Units
0:600 0:000 0:000 3:000 4:795 0:880
0:000 4:263 104 0:000 4:396 103 3:188 104 2:261 104
m6 kmol2 s1 s1 m3 kmol1 s1 m3 kmol1 s1 s1 s1
Forward reaction mechanism rate coefficients in form krþ ¼ aT b exp ðc=T Þ, 300 K 6 T 6 1780 K.
Table 2 Gas-phase reactions r 1 2 3 4 5 6
Reactions HþHþM H2 þ M C3 H8 CH3 þ C2 H5 C2 H5 þ H CH3 þ CH3 CH4 þ H CH3 þ H2 SiH4 SiH2 þ H2 SiH2 Si þ H2
a
b 13
2:714 10 3:194 103 4:158 106 1:835 102 2:868 1022 1:045 107
0:284 1:804 0.791 3:485 4:203 0.388
c
Units 4
5:210 10 8:215 102 3:363 103 3:544 103 6:279 103 2:391 103
Reverse reaction mechanism rate coefficients in form kr ¼ aT b expðc=T Þ, 300 K 6 T 6 1780 K.
m3 kmol1 s1 m3 kmol1 s1 m3 kmol1 s1 m3 kmol1 s1 m3 kmol 1 s1 m3 kmol1 s1
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
W_ s ¼ Ms
Nhr X
þ m sr msr
r¼1
kr
Ns Y qk M k k¼1
m Ns Y qk kr Mk k¼1 !
þ
krþ
mkr
;
ð6Þ
where Nhr is the total number of homogeneous reactions, mþ kr and mkr are the stoichiometric coefficients for the r-th reaction, and qk ¼ qck is the density of the k-th species. Thus, for simulating cold-wall SiC CVD, six gas-phase reactions and the following minimum set of species are chosen: H2 , H, C3 H8 , C2 H5 , CH4 , CH3 , SiH4 , SiH2 and Si. 2.3. Heterogeneous chemistry Generally, modeling of the chemical interaction between gaseous reactants and a solid surface involves a set of heterogeneous reactions. Surface reactions are characterized by complicated reaction mechanisms including a number of steps, and each surface reaction depends on the partial pressures of the gaseous species, the rate constants of each individual reaction step as a function of temperature and surface properties. However, there is usually little or no detailed information on the individual reaction step and the rate constants. Therefore, one of the steps is frequently assumed to be rate-limiting whereas all the other steps are taken to be in equilibrium (e.g., see [8]). With these assumptions, an expression for the overall reaction rate can be derived. A comparison of the predicted and measured growth rates is needed to show the validity of the assumptions and to evaluate unknown constants. In practice, this approach provides the availability of free parameters which should be adjusted to reproduce a definite set of experimental data. But the adjusted parameters are often found to be entirely unusable for the prediction of growth under different conditions. This is a crucial problem in modeling the heterogeneous chemistry for CVD. An alternative approach is based on the idea that the transport of species from the gas phase to the growing epitaxial layer is the stage limiting growth rate. The approach is feasible provided
523
that the characteristic time of all surface processes is much smaller than that of the gas-phase diffusion. In this case, near-equilibrium between the gaseous species and the crystal is established at the growth surface. Since the species transport in a cold-wall reactor is very slow and the growth temperature is high enough (1450 C), we assume that the collision and subsequent reaction of the gas molecules with the heated substrate surface are limited by mass transport. In this consideration, all the heterogeneous reactions run under nearly equilibrium conditions. In accordance with the Gibbs phase rule, we can write Nsr mass action law equations for all species involved in the surface reactions: Ns Y w mir xi ¼ Kr ;
r ¼ 1; 2; . . . ; Nsr
ð7Þ
i¼1
Here Nsr ¼ Ns þ Ncon Ne , Ncon and Ne are the number of condensed phases formed on the growth surface and the total P number of elements, respecNs tively, xwi ¼ ðqwi =Mi Þ= i¼1 qwi =Mi ¼ pi =P0 , qwi and pi are the molar fraction, the density and the partial pressure of the i-th species on the surface, respectively; mir is the stoichiometric coefficient of the i-th species in the r-th surface reaction, Kr is the equilibrium constant of the r-th surface reaction, determined from the thermodynamic properties of the relevant species. The total species flux onto the surface includes the convective ðqwi V w Þ and diffusion ðJiw Þ terms and is defined as: qwi V w þ Jiw ¼ Mi
Nsr X
mir w_ r ;
i ¼ 1; 2; . . . ; Ns ;
ð8Þ
r¼1
where V w is the Stephan velocity, w_ r is the rate of the r-th surface reaction. To close the set of Eqs. (7) and (8), the normalizing condition should be added: Ns X i¼1
xwi ¼
Ns X
pi =P0 ¼ 1
ð9Þ
i¼1
The set of non-linear Eqs. (7)–(9) must be solved with respect to ðNs þ Nsr þ 1Þ-variables on the solid surface: xwi , w_ r and V w .
524
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
The deposition rate is then found as: Vg ¼
Ns X
Mi
i¼1
Nsr X
mir w_ r =qsurf ;
ð10Þ
r¼1
where qsurf is the surface material density. The solution to the problem (7)–(9) is not influenced by the choice of reactions associated with the mass action laws in (7), provided they are independent. For instance, to describe SiC deposition on a SiC substrate, the following independent reactions can be chosen: CH4 þ Si $ 2H2 þ SiCðsolidÞ 2H $ H2 C3 H8 $ CH3 þ C2 H5 CH4 þ H $ CH3 þ H2 C2 H5 þ H $ 2CH3 SiH4 $ SiH2 þ H2 SiH2 $ Si þ H2 2.4. Gas-phase nucleation As mentioned above, a key feature of the process considered is homogeneous nucleation of silicon. There is a large amount of information on the subject. At present, there are two principal approaches to simulate a gas flow coupled with condensation. One is well-known and is based on the classical theory suggested by Volmer [9], and developed by Becker and D€ oring [10], Zeldovich [11] and Frenkel [12]. Despite the obvious drawbacks and a number of theoretical attacks, the classical theory is adequate in itself to predict the nucleation rate consistent with experimental data for many important processes under quasi-stationary conditions. However, nucleation is not the only stage of condensation. Nuclei can grow in a supersaturated vapor and, generally, their growth rate varies with the cluster size. Analyzing experimental data, Hill [13] concluded that the growth rate is independent of the radius for small clusters, and it is sufficient to employ equations for the first
three moments of the size-distribution function for modeling the drop evolution after the nucleation. This consideration significantly reduces the set of governing equations. The other approach using the discrete-sectional model has been lately applied to the problem discussed (e.g., see [14]). The method considers two parts of the simulated particle-size range: discrete and sectional. The discrete part handles the detailed interactions of monomers and small clusters up to a predetermined maximum cluster size. The upper-size range of the distribution is approximated by a small number of sections. However, the discrete-sectional model has been applied mostly to 1D flow of binary mixture because of the limited computational power [15,16]. By now, this approach has been hardly feasible to 2D and 3D flow of a multi-component gas mixture coupled with condensation of a species. Estimations show that under the operating conditions in the vertical cold-wall reactor, the drop size is a few orders of magnitude less than the mean free path of molecules in the vapor. Therefore, clusters can be considered as a pseudo-gas of heavy particles with the size-independent growth rate. Since the gas flow is slow, a steady-state sizedistribution of nuclei becomes achievable. In addition, our preliminary estimates and further 2D computations within the employed model show that such important mechanism as coagulation can be neglected under operating conditions typical for CVD of SiC in the considered reactor. The numerical concentration of clusters here does not exceed 1016 m3 , whereas their size is small (aver). Then, the estimation of age radius of 30–50 A mean free path of the particles gives 0.2 m and larger, while the characteristic size of the reactor is 0.1 m. So, the probability of collision between two neighboring clusters to generate a larger cluster is negligible. This results primarily from the fact that the precursors are highly diluted by the carrier gas, and the process is operated under low total pressures (50 Torr). In this situation, the approach suggested in [13] is suited to simulate the gas-phase condensation of silicon in a non-isothermal flow of a multi-component gas mixture. An approximate kinetic equation can be derived from either the Focker–Plank equation or
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
525
following simple considerations. Let us introduce the size distribution function (f) which is normalized as Z 1 4 qc ¼ pq0Si r3 f ðrÞ dr; ð11Þ 3 r
1 ðT Þ is supersaturation. Within the where s ¼ pSi =pSi assumptions mentioned above, the rate of mass exchange between the clusters and the gas phase can be simulated by using the Hertz–Knudsen equation; hence
where r is the particle radius, qc is the density of the particle pseudo-gas, q0Si is the liquid density. By neglecting the diffusion, this function can vary due to three main mechanisms: the particle motion in the gas flow, the gas-phase nucleation and the condensation/evaporation of gas molecules on/from the cluster surface. The balance of these mechanisms is expressed as o r ð Vc f Þ þ ð12Þ r_ f ¼ Idðr r Þ; or
r_ ¼
where Vc is the transport velocity of the particle pseudo-gas, r is the critical nucleus radius, I is the formation rate of critical-sized nuclei, r_ is the cluster-growth rate after nucleation. The process of steady-state nucleation is considered to be a flow of constant net current of nuclei through a steady distribution of nucleus sizes. In the other words, the net rate at which g-sized nuclei grow to ðg þ 1Þ-sized nuclei is equal to the net rate, at which g-sized nuclei grow from ðg 1Þ-sized ones, so that the concentration of nuclei of size g is constant. With this simplifying concept, the nucleation rate (I) is given by the Becker–D€ oring theory with Bernard’s correction factor [10]: sffiffiffiffiffiffiffiffiffiffiffiffiffiffi q2Si 2rSi NA3 4 rSi NA I¼a 0 exp p 3 RT qSi pMSi3 ð13Þ
r2 ð1 3g2=3 þ 2g1 Þ ; where a is the condensation coefficient, rSi is the surface tension coefficient of liquid silicon, NA is the Avogadro number, MSi is the silicon molar weight. The critical nucleus radius, r , and the number of molecules contained in the critical nucleus, g , are calculated from a maximum of the net change in free energy necessary for reducing the vapor pressure to the flat-film saturation value and for creating a liquid drop [13] r ¼
2rSi MSi ; q0Si RT ln s
4 NA q0Si 3 r ; g ¼ p 3 MSi
ð14Þ
1 a pSi pSi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi; 0 qSi 2pRT =MSi
ð15Þ
1 where pSi is the partial pressure of silicon and pSi is the flat-film saturation pressure. With account of Eqs. (13)–(15), the evolution equations for the first three moments of the sizedistribution function, Xn , are derived by multiplying Eq. (12) by rn and by integrating over r from r to 1
r ðVc Xn Þ ¼ rn I þ n_rXn1 ; where Xn is defined as Z 1 rn f ðrÞ dr Xn ¼
n ¼ 0; 1; 2;
ð16Þ
ð17Þ
r
Since a large temperature gradient is found near the heated substrate in a vertical cold-wall CVD reactor, the thermophoretic force can make the particle path deviate significantly from the actual gas flow streamlines [17–19]. The basic mechanism of thermophoresis is that a suspended particle in a temperature gradient experiences a net force in the direction of decreasing temperature. This force arises because molecules impacting the particle on opposite sides have different average velocities due to the temperature gradient. The thermophoretic velocity is derived theoretically and accounted for as [18] l Vc ¼ V KT rT ; ð18Þ qT where l is the gas mixture viscosity, KT is the thermophoretic coefficient. The value of this constant depends on the molecular flow regime. A measure of this regime is the Knudsen number––the ratio of the mean free path length of the gas molecule to the particle radius. For large Knudsen’s numbers, KT of 0.55 is predicted by the theory describing the motion of a single particle in the temperature gradient [17]. As follows from literature (e.g., see [18]), this value is chosen as a basic one and varied from 0.0 to 0.9 to bring the
526
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
computations into agreement with available experiment. In this study, the value of KT ¼ 0:2 appears most appropriate for the conditions encountered in the process considered. Transport of clusters is described by their continuity equation r ðqc Vc Þ ¼ W_ c ;
ð19Þ
where 4 W_ c ¼ pq0Si ðr3 I þ 3_rX2 Þ 3
ð20Þ
Expression (20) is introduced into Eq. (4) to provide the consumption or generation of a condensible species due to homogeneous nucleation or evaporation of clusters. Since the precursors SiH4 and C3 H8 are highly diluted by H2 under typical operating conditions, the corresponding terms associated with spontaneous condensation were neglected in Eqs. (1)–(3).
3. Results and discussion 3.1. Gas dynamics The CVD system elaborated by EMCORE Corp. (New Jersey) has been close to treat SiC growth. Precursors highly diluted by hydrogen carrier gas can be supplied into the reactor using either the whole area of the reactor inlet or a hole/ holes of a smaller diameter. The process is operated under the following conditions [2]: the total pressure of 50–300 Torr, the hydrogen flow of 15–30 slm, the susceptor rotation rate of 750 rpm, and the wafer temperature of 1725 K. The flow in this type of reactors is said to occur either in the forced convection mode, where streamlines are curved but non-circulating above the susceptor, or in a mixed convection mode, with undesirable vortices caused by the heated susceptor. The mode transition can be roughly estimated by using the well-known dimensionless flow parameters, namely, the Reynolds (Re) and Grashof (Gr) numbers [20]. For small values of the Gr=Re3=2 ratio, there is no flow recirculation above the rotating susceptor. However, the inlet gas velocity is an additional parameter that enters the
problem through the specification of the flow stability boundary. Obviously, some computations are needed to refine the boundary. An experimental and numerical analysis of the flow pattern as a function of the total pressure was carried out in [2]. It was shown that the mode transition occurs in the narrow range of 200–250 Torr for the total flow rate of 30 slm and the rotation rate of 750 rpm. With increasing total pressure, the flow becomes unstable and essentially 3D. Since most experimental data on the growth rate reported in [21] concern the total flow rate of 15 slm and the rotation rate of 750 rpm, we have found numerically the transition between the forced and mixed convection modes with the pressure variation. Our computations show that the flow remains to be stable up to 75 Torr. Fig. 1 gives an example of the velocity vector and temperature distribution at the total pressure of 50 Torr when the top inlet area is wholly used for supplying the gas into the reactor. One can see that there is a steep temperature gradient near the wafer heated to 1725 K and the flow is non-circulating everywhere, except the zone close to the reactor outlet, which is inessential to the epilayer growth. The typical mass concentration distribution of precursors and major species are shown in Figs. 2
Fig. 1. Stable flow pattern and temperature distribution.
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
527
central hole at the inlet. In this case, the inlet parameters were recalculated in such a way to keep the same mass fluxes of precursors and hydrogen supplied into the reactor. For this purpose, we varied the inlet density at the constant velocity over whole inlet area. As an illustration, let us consider the following typical operating conditions: • Total pressure ¼ 50 Torr • Flow rates: H2 ¼ 15 slm, C3 H8 ¼ 6:6 sccm, SiH4 ¼ 17 sccm • Inlet temperature ¼ 300 K Fig. 2. Mass concentration distribution of propane and methane.
and 3. It follows from the figures that propane and silane are completely decomposed near the heated substrate, producing CH4 and SiH2 . These species and atomic silicon provide SiC growth on a wafer placed at the center of the holder top. Fig. 4 illustrates the formation of a large vortex above the substrate at the total pressure of 100 Torr. Elevation of the total pressure to 200 Torr does not change the general features of the flow pattern, while the species mixing is intensified. The same analysis was carried out when the precursors go into the reactor only through the
Fig. 3. Mass concentration distribution of silane and silylene.
Calculated inlet parameters: 1. Whole inlet area is used: Inlet velocity ¼ 0:121148 m/s Mass fractions: H2 ¼ 0:972939, C3 H8 ¼ 0:009418, SiH4 ¼ 0:017643 Mixture density ¼ 0:005414 kg/m3 2. Central hole is used: Inlet velocity ¼ 0:121148 m/s At outer ring: Mass fractions: H2 ¼ 1:0, C3 H8 ¼ 0:0, SiH4 ¼ 0:0 Mixture density ¼ 0:005275 kg/m3 At central hole: Mass fractions: H2 ¼ 0:668616, C3 H8 ¼ 0:115333, SiH4 ¼ 0:216051 Mixture density ¼ 0:007675 kg/m3 Thus, the mixture density can rise then almost 1.5 times at the central hole. Note that the transition of the gas flow from the laminar regime to the unstable one is somewhat delayed. When the radius ratio of the central hole and the reactor inlet is 1:4, the flow remains to be non-circulating to 150 Torr at the propane and silane flow rates of 6.6 and 17 sccm, respectively. The flow stabilization is due to the fact that the increase in the gas mixture density at the reactor center tends to suppress vortices generated by the steep temperature gradient near the wafer. Analyzing the preliminary results, we limited the range of investigations and mostly concentrated on the mass transport under typical gas dynamic conditions: the total pressure of 50 Torr,
528
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
Fig. 5. Growth rate vs propane flow rate at SiH4 ¼ 17 sccm. The hole at the inlet is used for supplying SiH4 and C3 H8 . Full circles are for the experiment. Fig. 4. Unstable flow pattern and temperature distribution.
the susceptor rotation rate of 750 rpm, the wafer temperature of 1725 K, and the hydrogen flow rate of 15 slm. The propane and silane flow rates are allowed to vary in a wide range. 3.2. Verification of the model The model was verified by comparing the computations with the available experimental data on the SiC growth rate as a function of the input propane and silane flow rates. There is no detailed information on what inlet flange was used for supplying precursors into the growth chamber at the operating parameters mentioned above. For this reason, we carried out two sets of computations corresponding to the precursor supply through the central hole and the whole inlet area. The growth rate variation with the propane flow rate is plotted in Fig. 5. The precursors are supplied through the central hole of a radius four times smaller than that of the reactor inlet and the silane flow rate is kept at 17 sccm. As seen from the figure, the computations fit well the experiment. The sloping portion of the curve corresponds to regimes where silane is in excess and the growth rate is controlled by the propane flow rate. The growth rate begins to saturate at the input
carbon-to-silicon ratio C/Si 0.6 and can no longer be raised by elevating the propane flow rate. Normally, the saturation is due the transition from the propane-limited mode to the silane-limited and is observed at C/Si 1. We speculate that the earlier saturation in Fig. 5 is related to homogeneous nucleation in the gas. The main species assumed to be responsible for this phenomenon is silicon produced by the decomposition of silylene SiH2 Si þ H2 (see the reaction in Table 1). The silicon gasto-particle conversion depletes the vapor in the reactive species needed for heterogeneous surface reactions of the film growth. Support for the speculation comes from experimental evidence on the irradiating layer observed through the reactor viewport and attributed to black-body radiation scattering by silicon clusters [2]. The computed distribution of silicon cluster density over the substrate is shown in Fig. 6. It reproduces well the shape and the location of the irradiating layer in [21]. The result obtained suggests the following mechanism of the gas-phase depletion. The clusters affected by the thermophoretic force are accumulated in a thin layer located over the susceptor and cannot reach the growing surface. The gas flow takes the clusters away in the radial
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
529
Fig. 7. Growth rate vs propane flow rate at SiH4 ¼ 17 sccm. The whole inlet is used for supplying SiH4 and C3 H8 . Full circles are for the experiment. Fig. 6. Typical Si cluster density distribution.
direction (see Fig. 1), and, hence, some of the Si atoms stored in the clusters are lost to the SiC growth on the wafer. Note that the model neglecting silicon vapor condensation predicts saturation at C/Si P 1.0 and significantly higher growth rates [21]. A similar growth rate dependence corresponding to the use of the whole inlet section for supplying the precursors is presented in Fig. 7. A comparison of the curves in Figs. 5 and 7 shows that a higher propane flow rate of 5 sccm (C/ Si 0.8) is required to get the maximum growth rate in the latter case. This is due to the fact that some amount of precursors goes between the water-cooled wall and the susceptor as far as the outlet and does not reach the wafer. Thus, a careful choice of the inlet flange design can provide a higher precursor utilization and, hence, a higher growth rate. The growth rate is shown in Fig. 8 as a function of the silane flow rate at the constant propane flow rate of 6.6 sccm. The solid and dashed curves correspond to the supply of precursors through the central hole and the whole inlet section, respectively. As seen from the theoretical plots, a deviation of the growth rate dependence from the
Fig. 8. Growth rate vs silane flow rate at C3 H8 ¼ 6:6 sccm. Precursors are supplied through the hole (solid) and whole inlet (dashed). Full circles are for the experiment.
linear that occurs at the silane flow rate of 10 sccm. This is due to the onset of nucleation supported by the sharp increase of the cluster density over the substrate. In addition, the computations predict a decrease of the growth rate at high SiH4 flow rates. This behavior is associated with the greater Si vapor supersaturation and, therefore, with a higher gas-to-particle conversion of silicon.
530
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
Note that the theoretical curves shown in Fig. 8 are in poor agreement with the experimental data compared to the case of propane flow rate variation. In particular, the predicted growth rate begins to decrease at the SiH4 flow rate greater than 15 sccm, while the measured growth rate does not exhibit such a behavior. In our opinion, this discrepancy is related to an extremely high sensitivity of the growth rate to the detailed pathway of silicon nucleation. Indeed, the growth rate is directly affected by the material losses due to the gasto-particle conversion. The cluster nucleation, in turn, depends very non-linearly on the silicon vapor supersaturation. Therefore, small deviation of the real growth conditions from those set in computations may result in the considerable discrepancy between the theoretical predictions and experimental observations. In our case, the effect of unaccounted details of the growth process is found to be insignificant for the accurate prediction of the growth rate behavior at the propane flow rate variation. However, this effect is remarkable when the SiH4 flow rate is varied. On the other hand, the comparison of the computational and experimental data shows that the model suggested catches the major physical mechanisms underlying the CVD of SiC, as a whole. Therefore, we apply further the model to the computational analysis of the process by varying the precursor flow rates in a wider range. 3.3. Effect of precursor flow rate variation on the growth rate For getting a better insight into the growth rate behavior, we carried out computations, varying the precursor flow rate in the following ranges: C3 H8 ¼ 1:5–8:0 sccm, SiH4 ¼ 4:0–60:0 sccm. The growth rate variation with the propane flow rate is shown in Fig. 9(a) and (b) at several constant SiH4 flow rates. As seen from the figures, the growth rate initially rises at low propane flow rates and then reaches a maximum. This level increases in the range of the silane flow rate from 4.0 to 17 sccm. The higher is the silane flow rate, the larger is the propane flow rate at which the onset of the growth rate saturation occurs (see Fig. 9(a)). The maximum level and the saturation onset depend
Fig. 9. Growth rate vs propane flow rate at several constant SiH4 flow rates: (a) 4–17 sccm, (b) 17–60 sccm. The whole inlet is used for supplying SiH4 and C3 H8 .
inversely on the silane flow rate at SiH4 > 17 sccm (see Fig. 9(b)). This is due to the fact that the depletion of the gas mixture in silicon atoms involved in the growth process is enhanced by the silicon gas-to-particle conversion at higher silane flow rates. The growth rate is plotted in Fig. 10 as a function of the silane flow rate at several constant propane flow rates. The plateaus demonstrate the growth rate saturation at a low propane flow rate from 1.5 to 3.3 sccm where C3 H8 controls the process. At C3 H8 > 3:3 sccm, the plateaus disap-
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
Fig. 10. Growth rate vs silane flow rate at several constant C3 H8 flow rates. The whole inlet is used for supplying SiH4 and C3 H8 .
pear, and the growth rate dependence on silane flow rate peaks at SiH4 17 sccm. The subsequent decrease of the growth rate in both cases results from the gas-phase nucleation of silicon. At C3 H8 < 3:3 sccm, this occurs differently after the gas mixture becomes fairly depleted in silicon and the carbon-containing species becomes in excess. Note that an additional increase of the propane flow rate from C3 H8 ¼ 5:5 sccm and higher does not practically influence the maximum growth rate. For visualizing more clearly the trends discussed above, a 3D growth rate surface is plotted in Fig. 11 as a function of C3 H8 and SiH4 flow rates. It follows from the figure that the growth rate cannot exceed a critical value. The limitation is caused by the drastic silicon condensation in the vapor at high silane flow rates and means that it is impossible to provide a desirable increase in the growth rate via the precursor flow rate variation. Thus, the high flow rates of silane result in a greater loss of silicon atoms involved in SiC growth, and in spite of the inlet Si/C ratio being much more unity, the growth rate is then controlled by SiH4 . 3.4. Gas-phase nucleation effect on the Si/C ratio on the growth surface It is well-known that Si/C ratio plays a crucial role for many properties of SiC crystals. In par-
531
Fig. 11. 3D surface of growth rate vs SiH4 and C3 H8 flow rate.
ticular, the doping by various impurities, the concentration of carbon vacancies, surface morphology including step bunching, etc. are affected by the surface Si/C ratio which may differ significantly from the Si/C ratio at the reactor inlet. Since the Si/C ratio controls the adsorption layer composition on the growth surface, it should be directly related to the collision frequencies of molecules containing silicon and carbon atoms, respectively. So, we calculate the Si/C ratio as Pn Si fi bi pi ; ð21Þ Si=C ¼ Pi¼1 n C i¼1 fi bi pi where fiSi and fiC are the number of Si and C atoms in a molecule of the i-th species, respectively, 1=2 bi ¼ ð2pmi kT Þ and mi are the Hertz–Knudsen factor and the molecular mass of the i-th species, respectively, pi is the partial pressure of the i-th species. Normally, crystal growth is operated under excess of either silicon or carbon species. The growth rate is then controlled by supplying the species with a lower inlet concentration. Near the growth surface, the concentration of such species additionally decreases due to a more efficient incorporation into the crystal, as compared to the species in excess. For this reason, the actual Si/C ratio on the growth surface can differ significantly from that at the reactor inlet. Obviously, gas-phase processes like homogeneous nucleation
532
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
Fig. 12. Si/C ratio on the wafer and growth rate vs SiH4 flow rate at C3 H8 ¼ 6:6 sccm. No gas-phase nucleation is accounted for.
not contributing to the crystallization may alter additionally the Si/C ratio on the growth surface. The experimental [2] and predicted growth rates and Si/C ratio are plotted in Figs. 12 and 13, respectively, as a function of silane at the constant propane flow rate of 6.6 sccm with and without accounting for nucleation. Other operating parameters had their typical values: the hydrogen carrier gas flow rate of 15 slm, the total pressure of 50 Torr, and the growth temperature of 1450 C. One can see saturation of the deposition rate with
Fig. 13. Si/C ratio on the wafer and growth rate vs SiH4 flow rate at C3 H8 ¼ 6:6 sccm. Gas-phase nucleation is accounted for.
silane flow rate increase, usually observed in the experiments. Model results without the account of the gas-phase nucleation give a clear saturation of the dependence at the point where the inlet Si/C ratio becomes close to unity (Fig. 12). The predicted growth rates are higher than the experimental values. When the gas-phase nucleation is taken into account (Fig. 13), the consistence of the computed and experimental growth rates is better. What is important is that the growth rate saturation starts earlier. Note that the local values of the Si/C ratio over the wafer are also less than unity under the conditions with the inlet Si/C ratio larger than unity. The local Si/C ratio is predicted to decrease with silane supply at high inlet Si/C ratios. This effect results from the fairly high silane flow rate. To get a general pattern of the Si/C ratio near the growth surface, we made computations varying the precursor flow rates in the following ranges: C3 H8 ¼ 1:5–8:0 sccm, SiH4 ¼ 4:0–60:0 sccm. The Si/C behavior is illustrated in Fig. 14. Note that the proportion between silicon and carbon atoms here is an abrupt function and is either much less or much more than unity. This results from the accepted assumption of the near-equilibrium growth conditions at temperatures as high as 1450 C. The linear portion of the curve for Si=C ¼ 1 separates the SiC growth regimes when the silane flow rate is insufficiently high for the onset of nucleation. The Si/C ratio is then con-
Fig. 14. Si/C ratio near the growth surface vs precursor flow rates.
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
trolled completely by that at the inlet. The silicon gas-to-particle conversion manifests itself as the curve slope. The area of Si=C < 1 becomes larger than the area of Si=C > 1 due to the drastic depletion of the vapor in silicon near the growth surface. For this reason, the higher silane flow rate is unable to raise the Si/C ratio at SiH4 > 17 sccm. 3.5. Secondary phase formation The formation of parasitic (secondary) phases on the SiC surface (solid graphite and either liquid or solid silicon) is an undesirable effect leading to degradation of the material properties. Here we suggest a special criterion to predict whether a condensed phase is formed on the growth surface. We assume a secondary phase arisen to be in equilibrium both with the SiC crystal and with the common gas phase. Then, in addition to the heterogeneous chemical reaction producing crystalline SiC, we consider the reaction giving rise to the secondary phase. Let W_ SiC be the rate of a heterogeneous reaction producing a SiC crystal and W_ SP the rate of the secondary phase generation on the surface. These reaction rates can be calculated using the mass balance equations for the growth surface coupled with the mass action law equations corresponding to heterogeneous equilibrium vapor-crystal-secondary phase equilibrium. Depending on the rates of these reactions, the following cases are possible: 1
4. At W_ SP < 0 and W_ SiC < 0, SiC crystal evaporates congruently. The boundaries between these cases can be found from the equations W_ SP ¼ 0 and W_ SiC ¼ 0. The analysis of the secondary phase formation was made for the following variation of precursor flow rates: C3 H8 ¼ 0:5–8:0 sccm, and SiH4 ¼ 4:0–60:0 sccm. Fig. 15 shows by solid lines the boundaries between the zones where the secondary phase formation and/or growth of the SiC single crystal are predictable. For comparison, the similar boundaries computed with the neglect of the silicon condensation in the gas are shown by dashed lines. It is seen that without nucleation there would be a very narrow growth window where single-crystal SiC could be obtained. The silicon nucleation alters the Si/C ratio over the substrate and, hence, broadens the growth window significantly. It should be noted that the shape of the growth window is rather complex and is controlled totally by the silicon clustering. On the other hand, the growth rate is predicted to peak at SiH4 17 sccm (see Fig. 10 or 11). So, according to Fig. 15, the propane flow rate is allowed to vary in the range 3:3–5:5 sccm to avoid the parasitic phase formation on the growing surface. It can be also concluded from the comparison of the growth window with that computed with neglect of silicon clustering that homogeneous nucleation indirectly influences SiC crystal quality, and, undoubtedly,
1. At W_ SP < 0 and W_ SiC > 0, a single SiC crystal grows on the surface without secondary phase formation. 2. At W_ SP > 0 and W_ SiC < 0, the growth surface becomes completely covered by a secondary phase. 3. At W_ SP > 0 and W_ SiC > 0, simultaneous deposition of both SiC and secondary phase occurs on the growth surface. Then the surface fraction occupied by the secondary phase can be esti mated from the ratio as W_ SP = W_ SP þ W_ SiC .
1
In more detail, the criterion for a secondary phase formation, its verification and examples of application to various growth systems will be published in a separate paper.
533
Fig. 15. Secondary phase map on the growth surface.
534
A.N. Vorob’ev et al. / Computational Materials Science 24 (2002) 520–534
should be accounted for when optimizing the SiC CVD in a cold-wall vertical reactor.
4. Summary An advanced model of heat and mass transport coupled with spontaneous condensation of silicon in the gas is suggested to simulate SiC CVD in a cold-wall rotating-disc vertical reactor. The comparison of computed growth rates with available experimental data shows a reasonable agreement under typical operating conditions. It is found that homogeneous nucleation considerably influences the CVD process. In particular, it is impossible to provide an increase in the growth rate over some maximum value by varying the precursor flow rates due to a drastic silicon gas-phase condensation at high silane flow rates. The analysis shows that the silicon-to-carbon ratio near the growth surface can radically differ from that at the reactor inlet. In turn, this changes the proportion between the silicon and carbon atoms on the wafer and increases the range of silane and propane variation where the SiC single epilayer is expected to grow without secondary phase formation or morphology degradation. The model allows us to predict the optimal precursor flow rates, ensuring the maximum growth rate at relatively broad growth window.
References [1] A. Itoh, H. Matsunami, IEEE Electron Dev. Lett. 16 (1995) 280.
[2] R. Rupp, Yu.N. Makarov, H. Behner, A. Wiedenhofer, Phys. Stat. Sol. (b) 202 (1997) 281. [3] O. Kordina, C. Hallin, A. Ellison, A.S. Bakin, I.G. Ivanov, A. Henry, R. Yakimova, M. Touminen, A. Vehanen, E. Janzen, Appl. Phys. Lett. 69 (10) (1996). [4] O. Kordina, C. Hallin, A. Henry, J.P. Bergman, I. Ivanov, A. Ellison, N.T. Son, E. Jenzen, Phys. Stat. Sol. (b) 202 (1997) 321. [5] Yu.E. Egorov, A.I. Zhmakin, Comput. Mat. Sci. 11 (1998) 204–220. [6] C.D. Stinespring, J.C. Wormhoudt, J. Crystal Growth 87 (1988) 481–493. [7] M.D. Allendorf, R.J. Kee, J. Electrochem. Soc. 138 (3) (1991) 841–852. [8] M.E. Coltrin, R.J. Kee, G.H. Evans, J. Electrochem. Soc. 136 (3) (1989) 819–829. € bersattigten [9] M. Volmer, A. Weber, Keimbildung in u Gebilden 119 (1926) 277–301. [10] R. Becker, W. D€ oring, Kinetische Behandlung der Keim€bers€attigten D€ampfen, Annalen der Physik 24 bildung in u (1935) 719–752. [11] J. Zeldovich, J. Exp. Theor. Phys. (U.S.S.R.) 12 (1942) 525. [12] J. Frenkel, Kinetic Theory of Liquids, Dover, New York, 1969. [13] Ph.G. Hill, J. Fluid. Mech. 25 (3) (1965) 593–620. [14] J.J. Wu, R.C. Flagan, J. Colloid Interf. Sci. 123 (2) (1988) 339–352. [15] S. Yu, Y. Yoon, M. M€ uller-Roosen, I.M. Kennedy, Aerosol Sci. Technol. 28 (1998) 185–196. [16] S.-Y. Lu, H.-C. Lin, C.-H. Lin, J. Crystal Growth 200 (1999) 527–542. [17] L. Talbot, R.K. Cheng, R.W. Schefer, D.R. Willis, J. Fluid Mech. 101 (4) (1980) 737–758. [18] K.L. Walker, F.T. Geyling, S.R. Nagel, Thermophoretic deposition of small particles in the modified chemical deposition (MCVD) process, J. Am. Ceram. Soc. 63 (1980) 552–558. [19] D.I. Fotiadis, K.F. Jensen, J. Crystal Growth 102 (1990) 743–761. [20] G. Evans, R. Greif, J. Heat Transfer 109 (1987) 928–935. [21] A.N. Vorob’ev, Yu.E. Egorov, Yu.N. Makarov, A.I. Zhmakin, A.O. Galyukov, R. Rupp, Mat. Sci. Eng. B 61–62 (1999) 172–175.