Theoretical/experimental study of silicon epitaxy in horizontal single-wafer chemical vapor deposition reactors

Article abstract of DOI:10.1149/1.1393391

The main goal of this study is to examine the possibility of using detailed three-dimensional simulations of transport of momentum, energy, and mass in horizontal single-wafer epitaxial silicon reactors in conjunction with relatively simple kinetic models to describe the reactor's performance over the entire range of operating conditions. As the SiHCl₃-H₂ system is a widely used precursor for epitaxial silicon deposition in industrial applications, we have chosen to focus our model development on this system. In the development of the model we have considered the dependence of the gas properties on the gas composition as well as on the temperature. In addition, mass transport due to thermal diffusion has been considered. The accuracy of the simulation model has been examined by comparing the predicted silicon deposition rates and profiles in two commercial chemical vapor deposition (CVD) reactors with the experimentally measured values. A comparison of simulation and experimental results has indicated that a detailed transport model in conjunction with a Langmuir-Hinshelwood type kinetic model for silicon deposition accurately describes the epitaxial silicon deposition rate and deposition profile. In turn, this lumped reaction kinetic model has been used for optimization of commercially available horizontal CVD reactors for epitaxial deposition of silicon.

Full text of DOI:10.1149/1.1393391

1538
Journal of The Electrochemical Society, 147 (4) 1538-1550 (2000)
S0013-4651(99)06-070-X CCC: $7.00 © The Electrochemical Society, Inc.
A Theoretical/Experimental Study of Silicon Epitaxy in Horizontal
Single-Wafer Chemical Vapor Deposition Reactors
Srikanth Kommu,a,b,* Gregory M. Wilson,b,* and Bamin Khomamia,z
aDepartment of Chemical Engineering and the Materials Research Laboratory, Washington University, St. Louis,
Missouri 63130, USA
bMEMC Electronic Materials, Incorporated, Epitaxial Silicon Technology, St. Peters, Missouri 63376, USA
The main goal of this study is to examine the possibility of using detailed three-dimensional simulations of transport of momen-
tum, energy, and mass in horizontal single-wafer epitaxial silicon reactors in conjunction with relatively simple kinetic models to
describe the reactor’s performance over the entire range of operating conditions. As the SiHCl3-H2 system is a widely used pre-
cursor for epitaxial silicon deposition in industrial applications, we have chosen to focus our model development on this system.
In the development of the model we have considered the dependence of the gas properties on the gas composition as well as on the
temperature. In addition, mass transport due to thermal diffusion has been considered. The accuracy of the simulation model has
been examined by comparing the predicted silicon deposition rates and profiles in two commercial chemical vapor deposition
(CVD) reactors with the experimentally measured values. A comparison of simulation and experimental results has indicated that
a detailed transport model in conjunction with a Langmuir-Hinshelwood type kinetic model for silicon deposition accurately
describes the epitaxial silicon deposition rate and deposition profile. In turn, this lumped reaction kinetic model has been used for
optimization of commercially available horizontal CVD reactors for epitaxial deposition of silicon.
© 2000 The Electrochemical Society. S0013-4651(99)06-070-X. All rights reserved.
Manuscript submitted June 14, 1999; revised manuscript received November 15, 1999.
Chemical vapor deposition (CVD) is a process that uses chemi-
purpose of this study is to develop a detailed simulation model of sil-
icon epitaxy in atmospheric horizontal single-wafer CVD reactors.
cally reactive gases to deposit thin solid films on a surface. The use
of these thin solid films is widespread in many industries such as
microelectronics, magnetic materials, optical devices, automobiles,
and ceramics. In the silicon-based semiconductor industry CVD is
used to deposit a wide variety of films including polycrystalline sil-
icon (polysilicon), epitaxial silicon, silicon oxides, and silicon
nitrides. The advantage of the CVD process compared to other depo-
sition processes such as evaporation, sputtering, etc., is that it is
capable of producing thin films of a wide variety of materials with
precisely defined and highly reproducible electrical, optical, chemi-
cal, and mechanical properties.1
Background
The CVD process is schematically shown in Fig. 1. In general,
the highly interconnected processes that take place in a CVD reactor
can be summarized as follows2: convective and diffusive transport of
reactants from the reactor inlet to the wafer surface; chemical reac-
tions in the gas-phase leading to new reactive species and by-prod-
ucts; convective and diffusive transport of the reactants and reaction
products from the gas-phase reactions to the wafer surface; adsorp-
tion of gas-phase species on the wafer surface; surface diffusion of
adsorbed species; heterogeneous surface reactions leading to the for-
mation of a solid film; desorption of gaseous reaction products from
the wafer surface; convective and diffusive transport of reaction
products to the outlet of the reactor.
Different types of reactor chambers are used with different
ranges of operating conditions for the purpose of CVD, e.g., hori-
zontal reactors, vertical reactors, pancake reactors, barrel reactors,
continuous belt reactors, etc. However in recent years, the demand
for silicon wafers that have a large diameter and meet the stringent
specifications of the microelectronics industry is on the rise. This has
led to increased interest in single-wafer reactors that are capable of
producing high quality wafers, one wafer at a time. Hence, in this
study we have focused our modeling efforts on horizontal single-
The CVD process is a complex dynamical process involving
simultaneous mass, momentum, and energy transport, and complex
chemical reactions. Hence, the quality of the films produced will be
determined by the interactions of various transport processes and the
chemical reactions in the CVD chamber, which in turn depend on
process conditions such as flow rates, pressure, temperature, con-
centration of chemical species, reactor geometry, etc. The modern
microelectronics industry also put severe demands on the quality of
films produced by the CVD process. The CVD film deposition tech-
niques in the microelectronics industry need to fulfill general re-
quirements such as high deposition rate, good thickness uniformity
across a wafer (less than 2% thickness nonuniformity), minimum
particulate generation, and economic use of reactants. In addition,
the design and optimization of CVD reactors must fulfill the increas-
ing demands in the semiconductor industry. However, due to the lack
of detailed fundamental models the industrial CVD practitioners are
frequently forced to rely on methods of trial and error as well as on
statistical methods to create purely empirical models of the reactor
behavior. The empirical relations that are produced following these
procedures are difficult to use if the reactants or the reactor geome-
try is changed. On the other hand, mathematical modeling and sim-
ulation provide an excellent economic alternative to trial and error
based experimental techniques. The immediate benefits to be real-
ized are fewer experiments, reduction in waste during experimenta-
tion, and the ability to deal with different reactive species and reac-
tor geometries. In addition, simulation models can be used for opti-
mization of the process conditions and reactor geometry as well as
for design of new processes for high volume production. The main
* Electrochemical Society Active Member.
z
E-mail: bam@poly1.che.wustl.edu
Figure 1. Schematic representation of a typical horizontal CVD reactor.
Downloaded on 2015-06-17 to IP 130.194.20.173 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).

Journal of The Electrochemical Society, 147 (4) 1538-1550 (2000)
1539
S0013-4651(99)06-070-X CCC: $7.00 © The Electrochemical Society, Inc.
wafer epitaxial reactors. In what follows we briefly summarize the
most relevant prior work performed in this area to put into context
our own efforts.
have considered the dependence of the gas properties on the gas
composition as well as on the temperature. In addition, mass trans-
port due to thermal diffusion has been considered. To examine the
accuracy of the simulation model, we have performed experiments in
two commercial CVD reactors. Specifically, deposition rates have
been measured at atmospheric pressure and various wafer tempera-
tures as a function of inlet gas composition for a fixed inlet columet-
ric flow rate as well as a function of inlet gas flow rate for a fixed
inlet gas composition. A detailed comparison of simulation and
experimental results has indicated that a detailed transport model in
conjunction with a Langmuir-Hinshelwood type kinetic model for
silicon deposition accurately describes the epitaxial silicon deposi-
tion process. In turn, this lumped reaction kinetic model has also
been used for optimization of commercially available horizontal
CVD reactors used for epitaxial deposition of silicon.
One-dimensional as well as two-dimensional analyses of trans-
port of mass, momentum, and energy have been used to examine the
dynamics of various CVD reactors with different reaction kinetics.3-5
Although these analyses have provided significant insight on how
various parameters influence the dynamics of the deposition process,
they are inadequate in describing the dynamics of horizontal single-
wafer CVD reactors. This is mainly due to the presence of buoyancy
driven recirculations that are created as a result of temperature gra-
dients present in single-wafer reactors as well as due to a number of
other factors such as edge and entrance effects.6,7 Therefore, a num-
ber of previous investigators have focused their attention on devel-
opment of three-dimensional simulation models of horizontal epitax-
ial CVD reactors. The first of such studies were performed by Mof-
fat and Jenson8,9 where the epitaxial deposition of silicon in a CVD
reactor was simulated by using a reduced set of chemical reactions
for silane decomposition. In general, their studies have demonstrated
that three-dimensional transport effects play an important role in de-
termining the silicon deposition profiles even in the absence of buoy-
ancy driven cells due to the influence of different wall boundary con-
ditions. Kleijn and Hoogendoorn10 have also developed three-dimen-
sional simulation models of horizontal CVD reactor for predicting
GaAs deposition profiles. They demonstrated the relative importance
of convection, diffusion, and thermal diffusion and found that three-
dimensional models can accurately predict the experimental results
under widely varying conditions. However, their study was limited to
deposition processes that are completely transport limited.
Some of the above studies as well as others have demonstrated
that although complex multiple reaction chemistry plays an impor-
tant role in epitaxial CVD, deposition rates and profiles can be pre-
dicted accurately by models based on transport of only one reactive
species.5,8,10-13 Habuka et al.14-16 have used this strategy to model a
horizontal epitaxial silicon reactor. Specifically, they used a detailed
transport model to predict averaged silicon deposition rates for a
SiHCl3 and H2 precursor and found that the nonlinear increase in the
average silicon deposition rate is due to the increase in molecular
weight of the inlet SiHCl3 and H2 gas mixture which can be predict-
ed using a Eley-Rideal type decomposition model for SiHCl3. Habu-
ka et al.16 also compared their model predictions with experiments
in terms of an average silicon deposition rate that takes into account
the effect of wafer rotation. However, this is inappropriate since
wafer rotation only leads to azimuthal averaging but not radial aver-
aging of silicon deposition rates. A more appropriate comparison
would involve measuring deposition rates as a function of inlet gas
composition for a fixed inlet gas flow rate and as a function of inlet
gas flow rate for a fixed gas composition. By following this strategy
one can decouple the effect of gas flow rate and composition on sil-
icon deposition rates and hence perform a detailed comparison of
model predictions and experimental measurements.
Problem Formulation
In order to develop an accurate simulation model that fully de-
scribes the transport and deposition of chemical species in a CVD
reactor one needs to solve the appropriate set of governing equations
under realistic boundary conditions. These governing equations must
accurately describe the gas flow, transport of energy and chemical
species, and the chemical reactions in a CVD reactor. As shown by
earlier studies1,10 the dependence of gas properties on gas composi-
tion and temperature as well as transport due to thermal diffusion
should be considered. Taking all of these factors into consideration,
the governing equations have been arrived at by using the following
realistic assumptions for atmospheric pressure CVD: the gases are
considered to be ideal obeying the ideal gas law and Newton’s law of
viscosity; the gas mixture is assumed to behave as a continuum
(Knudsen number < 0.01); the heat generated/consumed by chemical
reactions is neglected as the reactants are highly diluted; heating due
to viscous dissipation is neglected (Brinkman number < 0.01); pres-
sure variations in the energy equation are neglected as the Mach
numbers are very small; the Reynolds number and the Rayleigh num-
ber are not very large so the gas flow is assumed to be laminar; the
Dufour effect which causes an energy flux in a gas mixture as a result
of concentration gradients, is neglected.
Governing equations.—We are interested in determining the
steady-state silicon deposition rate distribution on a wafer so we have
derived the governing set of equations at steady state. However, as the
governing set of equations and boundary conditions contain several
parameters, it is advantageous to write them in dimensionless form.
In order to make the governing set of equations dimensionless, the
following dimensionless variables are introduced
U
V
T
wA
P
Ref
␳V2
ˆ
ˆ
ˆ
ˆ
ˆ
wA
U ϭ
T ϭ
ϭ
P ϭ
ٌ ϭ ٌL
Twafer
wA,in
␳
␮
␭
CP
ˆ
ˆ
ˆ
␳ ϭ
ˆ
␮ ϭ
␭ ϭ
CP ϭ
␳
␮
␭
Ref
CP,Ref
Ref
Ref
The above summary has clearly shown that three-dimensional
modeling for CVD reactors for transport-limited systems10 as well
as for systems such as silane whose chemistry is well understood8
has been successfully used to describe the reactor dynamics. How-
ever, for systems such as the SiHCl3-H2 mixture whose chemistry of
decomposition is not well understood, the success of 3D modeling
efforts in describing the dynamics of the reactor is uncertain.
Our aim in this study is to examine the predictive capability (i.e.,
in terms of silicon deposition rate and deposition profiles) of detailed
three-dimensional simulation models of horizontal single-wafer epi-
taxial silicon reactors with simple kinetic models for decomposition
of SiHCl3. We have selected to study the SiHCl3-H2 mixture not only
because this system is widely used for growing epitaxial films in
industrial applications but also because the decomposition chemistry
of SiHCl3 is not well understood. Hence, by performing simulations
with various kinetic models and comparing the simulation results
with experimental measurements a more complete picture of SiHCl3
decomposition can be obtained. In the development of the model we
DAB
DAB,Ref
DAT
T
ˆ
ˆ
DAB
ϭ
DA
ϭ
[1]
DAT,Ref
All the reference properties are calculated at a reference temperature
Ref, V is the characteristic velocity of the gas, and L is the charac-
teristic reactor dimension. A summary of all the symbols used is
given in the List of Symbols at the end of this paper. Using the above
dimensionless parameters, the governing equations take the follow-
ing dimensionless forms
T
Mass balance
ˆ
ˆ
ˆ
ٌи␳U ϭ 0
[2]
[3]
Chemical species balance
ˆ
1
1
DT ˆ
ˆ
ٌи
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
ٌи(␳UwA ) ϭ
ٌи(␳DAB ٌwA ) ϩ
ٌT
ˆ
ReиSc
ReиTd
T
Downloaded on 2015-06-17 to IP 130.194.20.173 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).

1540
Journal of The Electrochemical Society, 147 (4) 1538-1550 (2000)
S0013-4651(99)06-070-X CCC: $7.00 © The Electrochemical Society, Inc.
Momentum balance
in atmospheric pressure CVD is given in Table I. In this study we
have considered the effect of all relevant dimensionless groups on
the flow, temperature, and concentration profiles in the reactor, and
by doing so we have fully characterized the reactor dynamics over
the entire range of operating conditions. The results of these studies
are presented in the section on Results and Discussion.
1
2
T
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ ˆ
ˆ
ˆ
ˆ
ˆ
ٌи(␳UU) ϭ
ٌи Ϫ ␮ ٌиU иI ϩ ␮ ٌиU ϩ (ٌиU)
[
]
Re
3
1
ˆ
T Ϫ
Gr
Re2
ˆ
2
ˆ
Reactor geometry.—Figure 2 is a schematic representation of the
ASM Epsilon One single-wafer epitaxial silicon reactor and Fig. 3 is
a the schematic representation of the Applied Materials’ Centura
reactor. These reactors have been used in our experimental studies of
epitaxial deposition of silicon on silicon wafers. The dimensions of
the reactors are shown in the figures, and a description of the com-
putational domain used to simulate the deposition process in these
reactors is given in the Results and Discussion section. In what fol-
lows, a brief description of how silicon is epitaxially deposited on a
wafer in these commercial reactors is given.
In both these reactors, the gases are injected at room temperature
(25ЊC) through the three inlets (inlet 1, inlet 2, and inlet 3) in the gas
injector. The gases then pass through the quartz reactor chamber,
which contains a silicon wafer placed on a graphite susceptor and sur-
rounded by a SiC ring. The wafer is heated radiatively by lamp arrays
placed outside the quartz reactor chamber such that the wafer tem-
perature is maintained a constant (1398 K is a typical wafer temper-
ature). The gases react in the chamber depositing the desired solid sil-
icon film on the surface of the wafer. Finally, the reactant and prod-
uct gases leave the reactor through the outlet, which is fixed at atmos-
pheric pressure. The exterior surfaces of the reactor are cooled with
recirculated air.
ˆ
Ϫ ٌP ϩ
ez [4]
1
ˆ
T Ϫ Ga ϩ 1
2
Energy balance
where
1
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ ˆ
ˆ
CP ٌи(␳UT) ϭ
ٌи(␳ٌT)
[5]
ReиPr
2
g␳Ref L3Twafer
␮
RefCP,Ref
␳RefVL
Re ϭ
Gr ϭ
Pr ϭ
2
␮
␭
Ref
␮
RefTRef
Ref
Twafer
Ga ϭ
[6]
[7]
TRef
␮
RefwA,in
DRTef
␮
Ref
Sc ϭ
Td ϭ
␳Ref DRef
In order to calculate the transport properties of the gas mixture in
the above equations, the transport properties of the constituent species
must be known. From the DIPPR database17 we have obtained the
data for specific heat, viscosity, and thermal conductivity of relevant
gas species as a function of temperature. In turn, we have used poly-
nomial functions to fit these data as functions of temperature (see
Appendix A). The binary diffusion coefficients and the thermal diffu-
sion coefficients were estimated using kinetic theory and the variation
of these properties with temperature have also been determined (see
Appendix B).
It can be seen from the above equations that the flow and the tem-
perature fields in the reactor are determined by four dimensionless
groups, viz., Re, Gr, Pr, and Ga while the species transport is deter-
mined by Re, Sc, and Td. The range of these dimensionless groups
Table I. Important dimensionless groups and their range in
atmospheric CVD reactors.
Dimensionless
group
Definition
Range
Figure 2. Schematic representation of ASM Epsilon One single-wafer epi-
taxial silicon reactor. All dimensions are in centimeters.
Reynolds number
(Re)
Schmidt number
(Sc)
Thermal diffusion
number (Td)
Inertial forces/viscous forces
5-100
0.1-10
0.1-10
Momentum diffusivity/chemical
species diffusivity
Momentum diffusivity/chemical
species diffusivity due to
temperature gradients
Heat transfer peclet Convective heat transfer/
number (PeH) diffusive heat transfer
Mass transfer peclet Convective mass transfer/
1.5-70
0.5-1,000
10-30
number (PeH)
Surface arrhenius
number (Ar)
Grashof number
(Gr)
diffusive mass transfer
Activation energy/potential
energy
Buoyancy forces/viscous forces 100-100,000
Prandtl number (Pr) Momentum diffusivity/thermal
0.3-0.7
1-1.3
diffusivity
Gay-Lussac
number (Ga)
Rayleight number
(Ra)
Temperature difference/reference
temperature
Buoyancy forces/viscous forces 30-70,000
Surface Damkohler Chemical reaction rate/diffusion
number (Da) rate
0.1-500
Figure 3. Schematic representation of Applied Materials’ Centura single-
wafer epitaxial silicon reactor. All dimensions are in centimeters.
Downloaded on 2015-06-17 to IP 130.194.20.173 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).

Journal of The Electrochemical Society, 147 (4) 1538-1550 (2000)
1541
S0013-4651(99)06-070-X CCC: $7.00 © The Electrochemical Society, Inc.
Boundary conditions.—Walls.—No slip boundary condition for
the gas velocity is specified
Two different solvers have been used throughout this study. One
solver is based on the Stone’s method19 and the other is based on the
algebraic multigrid method.20 The solver using Stone’s method gen-
erally proved to be faster than the solver based on the algebraic
multigrid method, however, the multigrid solver which solves the
discretized equations on a series of coarsening meshes proved to be
more robust for cases where the buoyancy term in the momentum
equation is dominant.
ˆ
U ϭ 0
[8]
Different thermal boundary conditions have been used on the walls
of the reactor, i.e.
ˆ
ˆ
ˆ ˆ
ˆ ˆ ˆ ˆ
ˆ
Inf
T ϭ TWall or nиٌT ϭ 0 or nиٌT ϭ h(T Ϫ T ) [9]
Zero total mass flux for each of the chemical species is specified
Experimental
ˆ
ˆ
nи( jAC ϩ jAT ) ϭ 0
[10]
To determine the variation in the silicon deposition rate as a func-
tion of different operating conditions (i.e., temperature of the wafer,
initial concentration of SiHCl3, gas flow rates, and the ratio of gas
flow rates across the three inlets of the reactors and rotation rate of
the wafer) a series of experiments in the commercial CVD reactors
shown schematically in Fig. 2 and 3 were conducted. The following
procedure has been used in performing all the epitaxial silicon depo-
sition experiments. In both reactors, first hydrogen gas was injected
into the reactor through the three inlets (inlet 1, inlet 2, and inlet 3) in
the gas injector at a constant flow rate (60 slm is a typical value of the
gas flow rate). Then the radiative lamps located outside the reactor
were switched on so that the silicon-coated graphite susceptor, sur-
rounded by the SiC ring was heated to about 700ЊC. Then the silicon
wafer is loaded into the reactor and placed on top of the susceptor,
and the power of the lamps is increased to obtain a wafer temperature
of 1125ЊC. Once this protocol was completed, the flow of hydrogen
gas into the reactor was replaced by a gas mixture of SiHCl3 and H2
(a typical value of flow rate is 80 slm of H2 and 20 g/min of SiHCl3).
The gases then pass through the quartz reactor chamber, which con-
tains the silicon wafer with a constant surface temperature of 1125ЊC.
The gases react in the chamber depositing the solid silicon epitaxial
film on the surface of the wafer. Clearly, the deposited layer thickness
depends on the deposition time. In all the experiments, the deposition
time was fixed at 100 s for all different combinations of H2 and
SiHCl3 flows. This leads to a steady-state epitaxial silicon deposition
rates of about 4 ␮m/min. After the deposition step, the flow of SiHCl3
is switched off so that only hydrogen gas flows through the reactor
and the temperature of the wafer is ramped down to about 700ЊC.
Then the wafer is unloaded and the epitaxial silicon thickness is
measured using Fourier transform infrared (FTIR) spectroscopy.21 It
should be noted that after each deposition step, the whole quartz
chamber is etched with HCl to remove any chamber deposits on the
reactor walls. In all the above experiments, reactant and product gases
leave the reactor chamber through the outlet, which is fixed at atmos-
pheric pressure and the exterior surfaces of the reactor are continu-
ously cooled with recirculated air. The results of these experiments
are reported in next section.
Wafer.—No slip boundary condition for the gas velocity is specified
or a wafer rotating speed is specified
ˆ
ˆ
U ϭ 0 or UWafer
[11]
An isothermal boundary condition is employed since the tempera-
ture of the wafer is maintained constant through radiative heating by
the lamps
ˆ
ˆ
T ϭ TWafer
[12]
Based on the stoichiometry and an overall surface reaction leading
to silicon deposition [i.e., SiHCl3(g) ϩ H2(g) r Si(s) ϩ 3HCl(g)],
the mass consumption and generation of each chemical species is
taken into account
ˆ
C
ˆ
T
ˆ
nи( jA ϩ jA ) ϭ RA
[13]
Inlets.—At the inlets the inlet velocity, temperature, and the compo-
sition of the gas mixture are specified, i.e.
ˆ
ˆ
in
U ϭ U
[14]
[15]
[16]
ˆ
ˆ
in
T ϭ T
ˆ
ˆ
CA ϭ CA,in
Outlets.—At the outlets the flow is assumed to be fully developed
ˆ ˆ
nиٌиU ϭ 0
[17]
ˆ ˆ
nиٌT ϭ 0
[18]
[19]
ˆ
ˆ
nиٌCA ϭ 0
Computational Details
Due to the very large temperature gradients present near the sili-
con wafer in the reactor, the body force term in the momentum bal-
ance equation has not been approximated using the Boussinesq
method. Instead the density in the buoyancy term is written as
Results and Discussion
␳ ϭ ␳0 ϩ (␳ Ϫ ␳0)
[20]
Influence of operating conditions on the reactor performance.—
In this section, we discuss the effect of operating conditions on the
dynamics of atmospheric CVD reactors. The purpose of this investi-
gation is to determine the operating conditions that will result in
deposition of a uniform silicon layer.
where ␳ is a reference density calculated at the reference gas tem-
perature0, and the term ␳0g is absorbed into the pressure gradient term
of Eq. 4. In turn, the set of governing equations have been dis-
cretized using the finite volume technique and the SIMPLEC algo-
rithm has been used for pressure corrections.17,18 Specifically, an in-
house developed Fortran code (i.e., for specific reaction kinetics)
coupled with a commercially available software package CFX 4.1c
from AEA Technology, Inc.,17 has been used to perform the compu-
tations. It should be noted that the discretization of the governing
equations results in a set of coupled nonlinear algebraic equations.
In order to reduce the required computational resources, a two-level
iteration technique has been employed to solve these nonlinear equa-
tions. Specifically, an inner iteration loop is used to solve for the spa-
tial coupling of each variable, and an outer loop is used to solve for
the coupling between variables. The convergence criteria for this
iterative procedure requires the error in the continuity equation to
fall below 0.1% of the total mass flow rate entering the system. In
addition, the residuals of all other variables such as velocity, pres-
In general the gases used in atmospheric CVD applications have
very similar Prandtl numbers, thus variations in the Prandtl number
does not play an important role in the proper description of the pro-
cess conditions (see Table I). Similarly the Gay-Lussac number only
varies by a small amount. Therefore, the flow in the reactor is basi-
cally governed by the Reynolds and the Grashof numbers while the
chemical species transport is dictated by Reynolds, Grashoff,
Schmidt, and thermal diffusion numbers. As is the case in most indus-
trial practices, the composition of SiHCl3-H2 gas mixture is chosen
such that SiHCl3 is dilute in H2. Hence it is possible to just solve for
the flow and temperature fields initially and then use the computed
flow and temperature fields to solve the chemical species balance
equations. However, the advantage of using this decoupled technique
in terms of overall central processing unit (CPU) time is very minor,
hence we have solved all the governing equations simultaneously.
sure, and enthalpy are required to be less than 10Ϫ8
.
Downloaded on 2015-06-17 to IP 130.194.20.173 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).

1542
Journal of The Electrochemical Society, 147 (4) 1538-1550 (2000)
S0013-4651(99)06-070-X CCC: $7.00 © The Electrochemical Society, Inc.
Figure 4. Domain discretization for the rectangular reactor geometry. The
total number of finite volume cells equals 8,400. Top: projection of the mesh
along the X-Y midplane; bottom: projection of the mesh along the X-Z
midplane.
Figure 5. Velocity vector and temperature contour plots in the rectangular
reactor geometry (Re ϭ 1, Gr ϭ 1000). (a) Velocity vector, Ux along the X-Z
midplane; (b) velocity vector, Uy along the Y-Z midplane; (c) temperature
contour along the X-Z midplane. The dimension along the Z axis has been
magnified five times. Lines with no arrows correspond to velocities perpen-
dicular to the plane in which plots have been made.
Figure 4 shows a typical mesh used in simulating the dynamics
of a simple rectangular reactor. This simple reactor geometry has
been selected to demonstrate the effect of different dimensionless
groups on the performance of the reactor. The total number of finite
volumes in the domain is equal to 8,400 and typical CPU times for
numerical convergence of the coupled overall mass, momentum, and
energy equations is on the order of 100 s on a Silicon Graphics Indi-
go2 workstation. Using this solution as an initial guess, and coupling
the chemical species balance equations with the conservation equa-
tions for mass, momentum, and energy leads to a converged numer-
ical solution in approximately 100 s. To ensure that the solutions
generated with the mesh shown in Fig. 4 are accurate, selected com-
putations using a much more refined mesh consisting of 200,000
cells were performed. Overall, we did not find a significant differ-
ence between the results from the two meshes, i.e., the maximum
difference between the computed results is less than 2%. Hence, the
results shown below are based on the mesh shown in Fig. 4.
ence the formation of buoyancy induced recirculations. In fact sim-
ulations under conditions where inertial effects are important indi-
cate the buoyancy driven recirculations can be avoided in this rec-
tangular reactor geometry if the Gr/Re2 is maintained below 10.
However, for some conditions which are uncommon for epitaxial
CVD, i.e., relatively high Ra (Ra > 1700-2800)5,10 the flow becomes
The effect of various dimensionless groups on the reactor perfor-
mance has been studied by varying the reactor dimensions, wall
boundary conditions, total pressure, as well as the carrier gas, and its
flow rate. As mentioned earlier, the most important parameters that
govern the reactor performance are Re, Gr, Sc, and Td. In Fig. 5 and
6, the effect of varying Re is illustrated. Specifically, Fig. 5a and b
show the velocity vectors parallel and perpendicular to the mean
flow direction under conditions of Re ϭ 1 and Gr ϭ 1000. As ex-
pected, for a large Grashof number and a small Reynolds number
recirculation cells in planes parallel and perpendicular to the mean
flow direction are observed. Figure 5c depicts the temperature con-
tours along the X-Z midplane plane. As clearly shown by this figure,
the temperature field is uniform above the wafer surface which is a
consequence of the fact that the heat transfer due to conduction is
dominant compared to convection (PeH ϭ 0.7). Figure 5c also shows
that near the edges of the wafer (where the recirculation cell is pre-
sent) the temperature contours have slightly expanded leading to
sections of uniform temperature in the recirculatory regions. Fig-
ures 6a and b demonstrate that as Re is increased (i.e., from Re ϭ 1
to Re ϭ 30), but Gr is held constant (i.e., Gr ϭ 1000) the recircula-
tion cells disappear, i.e., the flow profile is not greatly disturbed by
the sudden heating at the leading edge of the wafer. However, it
should be noted that the slight rearrangement observed in the veloc-
ity vectors is a direct consequence of the change in viscosity as a re-
sult of change in the temperature profile. These figures clearly show
that changes in the ratio of Gr/Re␤ (where 1 Յ ␤ Յ 2) strongly influ-
Figure 6. Velocity vector and temperature contour plots in the rectangular
reactor geometry (Re ϭ 30, Gr ϭ 1000). (a) velocity vector, Ux along the X-Z
midplane; (b) velocity vector, Uy along the Y-Z midplane; (c) temperature
contour along the X-Z midplane. The dimension along the Z axis has been
magnified five times. Lines with no arrows correspond to velocities perpen-
dicular to the plane in which plots have been made.
Downloaded on 2015-06-17 to IP 130.194.20.173 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).

Journal of The Electrochemical Society, 147 (4) 1538-1550 (2000)
1543
S0013-4651(99)06-070-X CCC: $7.00 © The Electrochemical Society, Inc.
unstable, and for these conditions the ratio of Gr/Re2 is not an appro-
priate indicator of the flow kinematics inside the reactor.
sition rate increases with an increase in Re (i.e., PeM is increased).
This is due to the increase in the rate of mass flux of SiHCl3 to the
wafer surface. Moreover, an increase in Re has a pronounced effect
on the deposition rate profile, i.e., with an increase in Re the varia-
tions in the deposition rate from the leading edge to the trailing edge
of the wafer decrease. These results can be rationalized in the fol-
lowing manner. For PeM ϭ 1.2, the residence time of the gas in the
reactor is greater in comparison to higher values of PeM as a result
more of SiHCl3 is converted near the leading edge of the wafer
hence, SiHCl3 is severely depleted down the length of the wafer. The
increase in the deposition rate at the trailing edge of the wafer for
PeM ϭ 1.2 is due to the presence of the buoyancy driven recircula-
tion cell in that location (see Fig. 5a).
In Fig. 8 the effect of variations of thermal diffusion (Td) on sili-
con deposition rate for a fixed Re ϭ 30, Sc ϭ 1.2, and Gr ϭ 1000 is
shown. Figure 8 clearly shows that when Td is small the silicon depo-
sition rate is unaffected by thermal diffusion, however, a significant
increase in Td leads to lower deposition rates. To clearly demonstrate
the effect of thermal diffusion one needs to examine the mass fraction
profiles of SiHCl3 inside the reactor under different conditions. Fig-
ure 9 shows the mass fraction profile of SiHCl3 along the X-Z mid-
plane passing through the center of the wafer in the reactor (for an ini-
tial SiHCl3 mass fraction ϭ 0.72). As shown by this figure, although
the initial mass fraction of SiHCl3 is 0.72, due to thermal diffusion
the SiHCl3 concentration near the surface of the wafer decreases.
This effect can be more clearly demonstrated by considering thermal
diffusion in the absence of surface reaction. Figure 10 shows the sim-
ulation results in the absence of surface reaction. Figure 10 clearly
demonstrates that the SiHCl3 concentration becomes smaller near the
surface of the wafer although no consumption of SiHCl3 at the hot
wafer surface is taking place. This is due to the fact that SiHCl3 mol-
ecules are much heavier than H2 molecules and the temperature gra-
dient between wafer and the surroundings will lead to the transport of
As mentioned earlier, heat transfer due to radiation is not consid-
ered in this study although it can be incorporated. However, we have
chosen not to include radiation effects in the simulations because it
would lead to a tremendous increase in the CPU time required for
the computations. Hence, an attempt has been made to determine an
appropriate wall boundary condition for heat transfer that can quan-
titatively describe the temperature distribution on the reactor cham-
ber surfaces. Specifically, by performing simulations with different
temperature boundary conditions, we have determined that a Robin-
type of boundary condition (see the third boundary condition in
Eq. 9) with a heat transfer coefficient of 50 W mϪ2 KϪ1, is adequate
in predicting the measured temperatures at the reactor walls. Hence,
this boundary condition has been used on the reactor walls for all the
simulations described below.
As the performance of the reactor is determined by Re, Gr, Pr,
Ga, Sc, and Td, we have systematically varied them to determine
their effect on silicon deposition rate on the wafer using a simple
reaction mechanism
SiHCl3 ϩ H2 r Si ϩ 3HCl
[21]
The rate of the above reaction is assumed to be first-order in SiHCl3
concentration. It should be noted that although the above kinetic
model does not fully describe the complex surface reactions that take
place in the deposition process it can be used to qualitatively inves-
tigate the effect of different parameters on silicon deposition rate and
thickness uniformity across a wafer. In all simulations discussed
below, the Damkohler number defined based on the inlet SiHCl3
concentration, height of the reactor, and the rate of surface reaction
was fixed at 100, i.e.
LRSi
Das ϭ
ϭ 100
[22]
C
D
SiHCl3,in SiHCl3ϪH
2
The reason for selecting this Damkohler number is that at higher
Damkohler numbers the rate of mass flux of SiHCl3 to the wafer sur-
face is the controlling factor for determining the silicon deposition
rate compared to the intrinsic kinetics of SiHCl3 decomposition.
Moreover, when the Damkohler number is small (i.e., Das << 1) the
gas flow will have little or no effect on silicon thickness profile on
the wafer. Under these conditions, the silicon deposition rate will be
dictated by the wafer surface temperature uniformity since the depo-
sition process is controlled by reaction kinetics.
Figure 7 shows the effect of Re on silicon deposition rate in the
flow direction along the centerline of the wafer for a fixed Gr ϭ
1000 and a fixed Sc ϭ 1.2. As shown by this figure the silicon depo-
Figure 8. Effect of thermal diffusion on silicon deposition rates in the flow
direction along the centerline of the wafer (Re ϭ 30, Gr ϭ 1000, Sc ϭ 1.2).
The dimension along the Z axis has been magnified five times.
Figure 9. Mass fraction profile of SiHCl3 along the X-Z midplane in the rec-
tangular reactor geometry in the presence of surface reaction (wSiHCl3, in ϭ
0.72). The dimension along the z axis has been magnified five times.
Figure 7. Silicon deposition rate in the flow direction along the centerline of
the wafer as a function of Re (Gr ϭ 1000, Sc ϭ 1.2).
Downloaded on 2015-06-17 to IP 130.194.20.173 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).

1544
Journal of The Electrochemical Society, 147 (4) 1538-1550 (2000)
S0013-4651(99)06-070-X CCC: $7.00 © The Electrochemical Society, Inc.
*
adsorbed molecules). The adsorbed SiCl2 species is most likely
bonded to the wafer surface with the two chlorine atoms pointing to
the gas stream. Therefore, after adsorption one ends up with a sili-
con surface covered with chlorine atoms. This adsorbed layer is then
*
attacked by the hydrogen atoms to reduce SiCl2 to Si and release
HCl gas according to the following reaction
*
SiCl2 ϩ H2(g) r Si(s) ϩ 2HCl(g)
[25]
The rate of adsorption of SiHCl3, Rad, on the silicon wafer surface
(see Eq. 24) is assumed to be first-order in SiHCl3 concentration and
is given by
Figure 10. Mass fraction profile of SiHCl3 along the X-Z midplane in the rec-
tangular reactor in the absence of surface reaction (wSiHCl3,in ϭ 0.72). The
dimension along the Z axis has been magnified five times.
R
ad ϭ Kad(1-X)CSsiHCl3
[26]
SiHCl3 molecules away from the hot wafer surface. Figure 10 also
shows that the concentration becomes even smaller beyond the trail-
ing edge of the wafer, which is due to the huge temperature gradients
in this region. Overall, these simulations clearly demonstrate that
thermal diffusion must be included in the simulation model if accu-
rate predictions of transport of the chemical species and the deposi-
tion rates are desired over the entire range of operating conditions.
In Eq. 26, X denotes the fraction of the sites occupied by SiHCl3,
CSsiHCl3 denotes the SiHCl3 concentration at the surface of the silicon
wafer, and Kad denotes the adsorption rate constant, which is a func-
tion of temperature. The rate of decomposition of the adsorbed
SiCl2 , which is the same as the rate of formation of silicon from
SiCl2 , is assumed to be first order in H2 concentration at the surface
*
*
of the silicon wafer and is given by
Silicon deposition rate measurements and predictions.—An
accurate kinetic model for SiHCl3 decomposition needs to be cou-
pled with the above described transport model to predict silicon
deposition rates as a function of different operating parameters.
There is a reasonable amount of information available regarding
intrinsic kinetic and thermochemical information on chlorosilanes in
the physical chemistry literature. Ho et al.22 have recently performed
a comprehensive review of prior studies23-28 dealing with gas-phase
and surface chemistry of SiHCl3, and proposed a reaction mecha-
nism for the SiHCl3 kinetics. The key observations made by Ho
et al.22 are as follows: the mass flux of SiHCl3 is greatly responsible
for the silicon deposition rate; the deposition rate is insensitive to the
rate of H2 and HCl adsorption/desorption; the concentration of gas-
phase intermediates such as SiCl2 is very small. Although it is pos-
sible to couple the reaction mechanism proposed by Ho et al.22 with
a detailed three-dimensional transport model, it is not practical to
perform reactor design and optimization with such a model since
many simulations have to be performed to deal with the uncertain-
ties present in the surface reaction rates. As an alternative we have
decided to use a lumped reaction kinetic model that can capture the
important observations made by Ho et al.22 First, an attempt was
made to use the simplest possible kinetic model which consists of an
overall surface reaction given by Eq. 21 where the rate of reaction is
assumed to be first order in SiHCl3 and H2 concentrations, i.e.
Rsi ϭ Kr XCHs 2
[27]
where CHs 2 denotes the concentration of hydrogen at the surface of
the silicon wafer and Kr denotes the desorption rate constant, which
is a function of temperature. The fraction of occupied sites at steady
state is obtained from the following relation
∂X
ϭ Rad Ϫ RSi
[28]
∂t
Using Eq. 26, 27, and 28
s
KadC
SiHCl3
X ϭ
[29]
KadCSsiHCl ϩ KrCHs
3
2
and from Eq. 27 and Eq. 29 the silicon deposition rate is obtained
s
KadC
SiHCl3
s
RSi
ϭ
[30]
KadC
KrCHs
SiHCl3
1 ϩ
2
Rsi ϭ Kf[SiHCl3][H2]
[23]
However, using the above reaction kinetics in the simulations did not
lead to a good comparison between model predictions and experi-
mental results (the results are shown in Fig. 13, 14, and 15). It should
be noted that Kf in Eq. 23 is determined using a regression analysis
which is described in more detail later in this section. Clearly, the
simple kinetic model shown in Eq. 23 is incapable of taking into
account one of the key observations made by Ho et al.,22 namely, the
deposition rate is sensitive to the desorption rate of H2. To take this
into account, we have used a Langmuir-Hinshelwood kinetic model
to describe the surface chemistry and in what follows, the approach
used to simplify the reaction chemistry of SiHCl3 decomposition to
the Langmuir-Hinshelwood kinetic model is described. The silicon
surface is assumed to have a fixed number of reactive sites. The reac-
tants, i.e., SiHCl3 molecules, are transported from the gas phase to
the surface of the silicon wafer. SiHCl3 molecules then attack the
unoccupied reactive sites and decompose at a certain rate according
to the following reaction
*
SiHCl3(g) r SiCl2 ϩ HCl(g)
[24]
Figure 11. Domain discretization of the ASM Epsilon One reactor geometry.
Total number of finite volume cells ϭ 43,414. Top: projection of the mesh
along the X-Y midplane; bottom: projection of the mesh along the X-Z mid-
plane.
Equation 24 describes the adsorption reaction mechanism where
SiCl2 is the adsorbed species (we assume that the adsorption is lim-
ited to only one monolayer, and there is no interaction between the
*
Downloaded on 2015-06-17 to IP 130.194.20.173 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).

Journal of The Electrochemical Society, 147 (4) 1538-1550 (2000)
1545
S0013-4651(99)06-070-X CCC: $7.00 © The Electrochemical Society, Inc.
Table II. Rate constants at different temperatures.
Temperature (K)
Kf [m4/(kmol s)]
Kad (m/s)
Kr (m/s)
1373
1398
1423
5.6
7.0
8.6
1.00
1.00
1.00
3.83 ϫ 10Ϫ3
4.23 ϫ 10Ϫ3
4.63 ϫ 10Ϫ3
and Kr. This can be accomplished by experimentally measuring the
silicon deposition rates on the leading edge of the hot wafer. The rea-
son for selecting the leading edge of the hot wafer is that the mass-
transfer resistance is almost negligible at this position due to the neg-
ligible thickness of the mass-transfer boundary layer. So one can
assume that the concentration of SiHCl3 at the leading edge is close
to the initial concentration of SiHCl3. Therefore, by measuring the
silicon deposition rate at the leading edge of the hot wafer at differ-
ent temperatures one can determine the values of Kad and Kr as func-
tions of temperature. However, it is very difficult to experimentally
measure the deposition rate at the leading edge of the hot wafer
because in the commercial reactors used to perform experiments, the
wafer sits on a hot graphite susceptor that is surrounded by a SiC
ring. Due to this limitation we have decided to couple the proposed
kinetics shown by Eq. 30 with a full three-dimensional simulation of
mass, momentum, and energy transport in the reactor to determine
Kad and Kr. Specifically, simulations have been performed for vari-
ous sets of Kr and Kad at a fixed temperature. Figures 11 and 12 show
the domain discretizations used to perform simulations of the ASM
Figure 12. Domain discretization of the Applied Materials’ Centura reactor
geometry. Total number of finite volume cells ϭ 46,320. Top: projection of
the mesh along the X-Y midplane; bottom: projection of the mesh along the
X-Z midplane.
In order to couple this phenomenological kinetic model of silicon
deposition rate to the transport model one needs to determine Kad
Figure 13. Theoretically predicted and experimentally measured silicon deposition rates in the flow direction along the centerline of the wafer as a function of
wafer temperature for a fixed inlet volumetric flow rate of 1.7 ϫ 10Ϫ3 m3/s and a fixed inlet gas composition of wSiHCl3, in ϭ 0.71. Deposition rate is defined
as the ratio of the deposited film thickness to the deposition time. Measurements from the two reactors are indistinguishable. (a) T ϭ 1398 K; (b) T ϭ 1423 K;
(c) T ϭ 1448 K.
Downloaded on 2015-06-17 to IP 130.194.20.173 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).

1546
Journal of The Electrochemical Society, 147 (4) 1538-1550 (2000)
S0013-4651(99)06-070-X CCC: $7.00 © The Electrochemical Society, Inc.
Epsilon One and the Applied Materials Centura reactors, respective-
ly. To investigate the effect of various operating conditions on the
silicon deposition rate a number of simulations have been performed
with these meshes. To ensure that the solutions generated with the
meshes shown in Fig. 11 and 12 are accurate, selected computa-
tions using a much more refined mesh consisting of approximately
200,000 cells were performed. Overall, we did not find a significant
difference between the results from the two meshes, i.e., the maxi-
mum difference between the computed results is less than 2%.
Hence, the results shown below are based on the meshes shown in
Fig. 11 and 12.
The computed and measured silicon deposition rates along the
centerline of the wafer as a function of volumetric flow rate and gas
composition have been used to determine Kr and Kad at a fixed tem-
perature. Specifically, a regression analysis was used to minimize
the error between the simulated results and the experimental meas-
urements as a function of Kad and Kr
should be noted that Kf in the simple kinetic model (see Eq. 23) is
also determined by a similar regression analysis. Table II shows the
adsorption and desorption rate constants and Kf at the three different
wafer temperatures considered. From the values of Kad and Kr shown
in Table II it is evident that the silicon deposition rate is more sensi-
tive to the desorption rate constant compared to the value of the
adsorption rate constant. This observation is fully consistent with the
observations made by Ho et al.22
Figure 13, 14, and 15 show the comparisons between model pre-
dictions and experimental results. The experimental results shown in
Fig. 13, 14, and 15 are typical silicon deposition rates obtained under
normal operating conditions from the ASM and the Centura reactors.
Figure 13 shows the model prediction of silicon deposition rate as a
function of wafer temperature along the centerline of the wafer sur-
face for a nonrotating wafer for a fixed inlet volumetric flow (1.7 ϫ
10Ϫ3 m3/s) and fixed inlet composition (SiHCl3 mass fraction ϭ
0.71). It should be noted that the units of silicon deposition rate in
the plot are microns/minute and the units of RSi are converted to
microns/minute as G(␮m/min) ϭ RSiи6.0и107 MwSi/␳Si. Figure 13
shows that the silicon deposition rate increases with an increase in
the temperature of the wafer although the rate of increase is not very
high in the temperature range considered in this study. Figure 13 also
demonstrates that the Langmuir-Hinshelwood model is capable of
predicting the experimentally observed decrease in silicon deposi-
tion rate from the leading edge to the trailing edge of the wafer at
different temperatures. However, the linear kinetic model is inca-
xϭN
2
Error ϭ
(Gx,model Ϫ Gx,Experimental)
[31]
∑
xϭ1
In Eq. 31, Gx is the silicon deposition rate at any position x along
the centerline of the wafer and N is the total number of points on the
wafer surface at which the deposited layer thickness is measured.
The same procedure was repeated at different wafer temperatures. It
Figure 14. Theoretically predicted and experimentally measured silicon deposition rates in the flow direction along the centerline of the wafer as a function of
inlet gas composition for a fixed inlet volumetric flow rate of 1.7 ϫ 10Ϫ3 m3/s and a fixed wafer temperature of 1398 K. Deposition rate is defined as the ratio
of the deposited film thickness to the deposition time. Measurements from the two reactors are indistinguishable. (a) wSiHCl , in ϭ 0.65; (b) wSiHCl3, in ϭ 0.71;
3
(c) wSiHCl3, in ϭ 0.76.
Downloaded on 2015-06-17 to IP 130.194.20.173 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).

Journal of The Electrochemical Society, 147 (4) 1538-1550 (2000)
1547
S0013-4651(99)06-070-X CCC: $7.00 © The Electrochemical Society, Inc.
pable of predicting this trend. This is due to the fact that the desorp-
tion mechanism is not accounted for in this model. Figure 14 depicts
the model predictions of the silicon deposition rate as a function of
the inlet SiHCl3 composition along the centerline of the wafer sur-
face for a nonrotating wafer under conditions of a fixed inlet volu-
metric flow rate (1.7 ϫ 10Ϫ3 m3/s) as well as a fixed wafer temper-
ature of 1398 K. This figure clearly shows that, as expected, the
silicon deposition rate increases with an increase in the inlet SiHCl3
concentration. Moreover, the proposed Langmuir-Hinshelwood
kinetic model leads to excellent predictions of experimental results
at the three different SiHCl3 concentrations while the simple linear
kinetic model predictions are not very good. Figure 15 depicts the
model prediction of silicon deposition rate as a function of inlet vol-
umetric flow along the centerline of the wafer surface for a nonro-
tating wafer for a fixed inlet composition (SiHCl3 mass fraction ϭ
0.71) and a fixed wafer temperature of 1398 K. As shown by this fig-
ure, an increase in volumetric flow rate leads to an increase in sili-
con deposition rate which essentially means that the rate of transport
of chemical species plays an important role in determining the sili-
con deposition rate. Once again the Langmuir-Hinshelwood type of
kinetic model predicts the experimental results well for the three vol-
umetric flow rates considered. Overall, the above comparisons clear-
ly demonstrate that a detailed transport model coupled with the
Langmuir-Hinshelwood type kinetics is capable of predicting silicon
deposition rates over a wide range of processing conditions and
hence can be used to effectively design and optimize CVD reactors
with the SiHCl3-H2 reactant system.
Although it is not possible to separate exactly the influence of mass
transfer and kinetics on the silicon deposition rate from the above
analysis, the simulations results clearly suggest that the silicon depo-
sition rate in the region around the leading edge of the wafer is more
sensitive to the reaction kinetics and the deposition rate at the trailing
edge is more sensitive to the transport properties. This can be ex-
plained by the fact that near the leading edge the mass-transfer bound-
ary layer is thin, hence the overall deposition rate is governed by the
reaction rate. However, at the trailing edge the thickness of the mass-
transfer boundary layer is significant, hence, the ability to rapidly get
the reactants to the wafer surface becomes the controlling factor.
Effect of Wafer Rotation Rate on the Thickness Uniformity
Our experimental studies indicate that by rotating the wafer and
adjusting the flow rate ratios across the three reactor inlets the thick-
ness uniformity of the deposited silicon can be optimized. We have
also observed that the simulated deposition rates and profiles for
rotating and static wafers are very similar when the wafer is rotated
with a speed in the range of 5-50 rpm (see Fig. 16). This is due to the
fact that in this range of rotation speeds (i.e., typical range used in
industrial reactors) wafer rotation affects the flow field only in the
close proximity of the wafer surface, hence the mass-transfer bound-
ary layer thickness is almost unaffected. As a result, the bulk SiHCl3
Figure 15. Theoretically predicted and experimentally measured silicon deposition rates in the flow direction along the centerline of the wafer as a function of
inlet volumetric flow rate for a fixed inlet gas composition of wSiHCl , in ϭ 0.71 and a fixed wafer temperature 1398 K. Deposition rate is defined as the ratio of
the deposited film thickness to the deposition time. Measurements 3from the two reactors are indistinguishable. (a) 1.1 ϫ 10Ϫ3 m3/s; (b) 1.7 ϫ 10Ϫ3 m3/s; (c)
2.4 ϫ 10Ϫ3 m3/s.
Downloaded on 2015-06-17 to IP 130.194.20.173 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).

1548
Journal of The Electrochemical Society, 147 (4) 1538-1550 (2000)
S0013-4651(99)06-070-X CCC: $7.00 © The Electrochemical Society, Inc.
Figure 16. Comparison of silicon deposition rates in the flow direction along
the centerline of the wafer with and without wafer rotation in the ASM reac-
tor. Volumetric flow rate ϭ 1.7 ϫ 10Ϫ3 m3/s; T ϭ 1398 K; wSiHCl3, in ϭ 0.76.
concentration profile (i.e., SiHCl3 concentration profile at a height
midway from the wafer surface to the top of the reactor) remains the
same irrespective of the wafer rotation speed (see Fig. 17).
Clearly, in the absence of wafer rotation, the deposition rate is
nonuniform both in the direction of the flow (see Fig. 16) as well as
in the transverse direction to the flow (see Fig. 18a). Hence, wafer
rotation is required if good thickness uniformity is to be achieved
since rotation translates any location on the wafer through both high
(leading edge) and low (trailing edge) deposition regimes leading to
azimuthal averaging. Hence, to compare the computed deposition
rates with the experimental ones, the continuous change in the posi-
tion of the wafer due to wafer rotation must be taken into considera-
tion, i.e., the film thickness at a given position on the wafer should
be calculated by integrating the deposition rate over all the positions
that are at the same distance from the wafer center. After integrating
the local deposition rates along different radii, we have found that the
numerical results using Langmuir-Hinshellwood kinetics compared
very well with the experimental measurements (see Fig. 18b). Fig-
ure 18b also shows that for rotation speeds less than 50 rpm the depo-
sition rate as function of radial position from the center is identical
Figure 17. SiHCl3 mass-fraction profile along the X-Y midplane in the ASM
reactor. Volumetric flow rate ϭ 1.7 ϫ 10Ϫ3 m3/s; T ϭ 1398 K; wSiHCl3, in ϭ
0.67. (a) Stationary wafer; (b) wafer rotated at 35 rpm.
for a rotating wafer and a stationary wafer when the deposition rate
has been averaged in the azimuthal direction (thickness nonunifor-
mity is defined as the percent difference between the thickness of the
deposited layer at any location on the wafer with respect to the thick-
ness of the deposited layer at the center of the wafer), confirming the
fact that wafer rotation mainly leads to azimuthal averaging of the
thickness of the deposited layer.
As wafer rotation mainly leads to azimuthal averaging of the
thickness of the deposited layer, other strategies should be explored
to control the silicon thickness uniformity in the radial direction. As
mentioned above, this can be accomplished by adjusting the inlet
flow rate ratios across the wafer. In fact it can be shown that adjust-
Figure 18. Silicon deposition rates as a function of wafer rotation speed. Volumetric flow rate ϭ 1.7 ϫ 10Ϫ3 m3/s, T ϭ 1398 K; wSiHCl3, in ϭ 0.76. Flow rate
ratio across the inlets ϭ 1.5:1. (a) Comparison of computed and experimentally measured silicon deposition rates in the transverse direction to the mean flow
along the centerline on a stationary wafer. (b) Comparison of computed thickness nonuniformity in the radial direction as a function of wafer rotation rate. The
center of the disk corresponds to position 100.
Downloaded on 2015-06-17 to IP 130.194.20.173 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).

Journal of The Electrochemical Society, 147 (4) 1538-1550 (2000)
1549
S0013-4651(99)06-070-X CCC: $7.00 © The Electrochemical Society, Inc.
Washington University assisted in meeting the publication costs of this
article.
ing the inlet flow ratios has a major effect on the shape of the depo-
sition rate profile in the radial direction. Figure 19 compares the
numerical and experimental silicon thickness nonuniformity profiles
obtained when the wafer is rotated with different inlet flow rate
ratios. Figure 19 shows that for a wafer rotation speed of 35 rpm,
when the gas velocity through all the three inlets is maintained the
same a thickness nonuniformity of 7% is obtained with the thickness
at the edge of the wafer being greater than the thickness at the cetner
giving rise to a U-shaped thickness profile. This is very undesirable
because the thickness nonuniformity should generally meet a specif-
fication of less than 2% for an industrial accepted quality epitaxial
silicon wafer. Figure 19 shows that when the gas velocity through
inlet 2 is increased to twice the gas velocity through inlet 1 and inlet
3 an inverted U-shaped thickness profile with a thickness nonunifor-
mity close to 7% is obtained. Figure 19 also shows that for a flow
rate ratio of 1.5:1 between the central to the side inlets an acceptable
thickness nonuniformity (i.e., a less than 2% variation across the
wafer) is obtained. The above simple analysis shows that by adjust-
ing the gas flow rates through the three gas inlets (i.e., varying the
rate of mass transport to different sections of the wafer) one can opti-
mize the epitaxial thickness uniformity on the wafer.
Appendix A
Specific heat as a function of temperature
[J/(kg K)]
Gases
Hydrogen
1.408и104 ϩ 6.33и10Ϫ1 T
Trichlorosilane 5.96и102 ϩ 6.76и10Ϫ2 T
HCl
8.09и102 Ϫ 5.5и10Ϫ2 T ϩ 7.20и10Ϫ5T2 Ϫ 1.34и10Ϫ8T3
The specific heat of the gas mixture is calculated using the following rela-
tionship
CpSiHCl ϪH ϭ wSiHCl CpSiHCl ϩ wH CpH
2
[A-1]
3
2
3
3
2
Viscosity as a function of temperature
(Pa s)
Gases
Hydrogen
3.0и10Ϫ06 ϩ 2.0и10Ϫ08 (T/298) Ϫ 3.0и10Ϫ12 (T/298)2
Trichlorosilane 2.0и10Ϫ7 ϩ 4.0и10Ϫ8 (T/298) Ϫ 5и10Ϫ12 (T/298)2
HCl
4.0и10Ϫ07 ϩ 5.0и10Ϫ8 (T/298) Ϫ 8.0и10Ϫ12 (T/298)2
The viscosity of the gas mixture is calculated using the following relation-
ship:
Summary
n
xi␮
A simulation model for horizontal single-wafer epitaxial silicon
atmospheric CVD reactors has been developed. In the development
of the model the dependence of the gas properties on the gas com-
position and temperature as well as mass transport due to thermal
diffusion have been considered. As SiHCl3 is a widely used precur-
sor for growing epitaxial films in industrial applications the simula-
tion model has been used to study the epitaxial silicon deposition
from a SiHCl3-H2 precursor. To examine the predictive capability of
the simulation model, a selected number of silicon epitaxial deposi-
tion experiments were performed in two commercial CVD reactors.
A detailed comparison of simulation and experimental results has in-
dicated that a detailed transport model that includes mass transport
due to thermal diffusion in conjunction with a Langmuir-Hinshel-
wood type kinetic model for SiHCl3 decomposition accurately de-
scribes the epitaxial silicon deposition process in a broad range of
operating conditions. In turn, this lumped reaction kinetic model has
been successfully used for optimization of commercially available
horizontal CVD reactors.
i
␮
ϭ
mix
n
[A-2]
[A-3]
∑
iϭ1
xjΦij
∑
jϭ1
where
2
1
2
1
2
1
4
Ϫ
Mj
Mi
1
Mi
Mj
␮
i
Φij
ϭ
1 ϩ
1 ϩ
␮
8
j
Thermal conductivity as a function of temperature
(W/m K)
Gases
Hydrogen
Trichlorosilane
HCl
4.7и10Ϫ02 ϩ 5.0и10Ϫ4 (T) Ϫ 6.0и10Ϫ8 (T)2
6.2и10Ϫ03 ϩ 5.0и10Ϫ5 (T) Ϫ 7.0и10Ϫ9 (T)2
1.0и10Ϫ04 ϩ 5.0и10Ϫ5 (T) Ϫ 3.0и10Ϫ9 (T)2
The thermal conductivity of the gas mixture is calculated using the following
relationship
Acknowledgments
n
n
Ϫ1
xi
This work has been supported by a grant from MEMC Electron-
ic Materials, Inc. We are grateful to Tom Torack for performing the
deposition experiments as well as Dr. Robert W. Standley and Dr.
Sriram Ramanathan for helpful comments.
␭ ϭ ␰
␭ixi ϩ (1 Ϫ ␰)
[A-4]
[A-5]
∑
∑
␭
i
iϭ1
iϭ1
where
␰ ϭ 0.312 ϩ 0.325xH Ϫ 0.311xH ϩ 0.469xH
2
2
2
The expressions A-2 to A-4 are obtained from Bretsznadjer.29
Appendix B
Binary diffusion coefficients as a function of
temperature (m2/s)
Gas pairs
Hydrogen-hydrogen
4.0и10Ϫ5 (T/298) ϩ 10Ϫ4 (T/298)2
Ϫ 8.0и10Ϫ6 (T/298)3 ϩ 5.0и10Ϫ4 (T/298)4
Trichlorosilane-hydrogen 7.0и10Ϫ6 (T/298) ϩ 3.0и10Ϫ5 (T/298)2
Ϫ 3.0и10Ϫ6 (T/298)3 ϩ 2.0и10Ϫ4 (T/298)4
HCl-hydrogen
2.0и10Ϫ5 (T/298) ϩ 5.0и10Ϫ5 (T/298)2
Ϫ 3.0и10Ϫ6 (T/298)3 ϩ 7.0и10Ϫ8 (T/298)4
The thermal diffusion coefficients of the gas species are calculated using the
following relationships
P
DSTiHCl
ϭ
MSiHCl MH DSiHCl 2 kSiHCl
ϪH ϪH
3 2 3 3
[B-1]
[B-2]
3
2
MRT
Figure 19. Theoretically predicted and experimentally measured silicon
thickness profiles as a function of inlet gas flow rate ratios in the ASM reac-
tor. Volumetric flow rate ϭ 1.7 ϫ 10Ϫ3 m3/s; T ϭ 1398 K; wSiHCl3, in ϭ 0.67.
P
DHTCl
ϭ
MHClMH DHClϪH kHClϪH
2 2 2
MRT
Downloaded on 2015-06-17 to IP 130.194.20.173 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).

1550
Journal of The Electrochemical Society, 147 (4) 1538-1550 (2000)
S0013-4651(99)06-070-X CCC: $7.00 © The Electrochemical Society, Inc.
where the thermal diffusion ratios (kHCl-H and kSiHCl -H2) are obtained as a
function of the gas species concentrations2and temper3ature from Hirshfelder
et al.30
␳
␮
␭
^
density of gas (kg mϪ3
)
viscosity of gas (kg mϪ1 sϪ1
)
thermal conductivity of gas (W KϪ1
dimensionless parameter
)
List of Symbols
CA
CpA specific heat of chemical species A (J kgϪ1 KϪ1
Da Damkohler number
molar concentration of chemical species A (kmol mϪ3
)
References
)
1. C. R. Kleijn, Ph.D. Thesis, Delft University of Technology, Delft, The Netherlands
(1991).
DAB binary diffusion coefficient of chemical species A in diluent B (m2 sϪ1
)
2. K. F. Jensen, H. Simka, T. G. Mihopoulos, P. Futerko, and M. Hierlemann, in
Advances in Rapid Thermal and Integrated Processing, F. Roozeboom, Editor,
Series E: Applied Sciences, Vol. 318, p. 305, Kluwer Academic Publishers, Ams-
terdam (1995).
3. H. Schlichting, Boundary-Layer Theory, 7th ed., McGraw-Hill, New York (1968).
4. G. Evans and R. Grief, Num. Heat Transfer, 14, 373 (1988).
5. V. J. Ven, G. J. M. Rutten, M. J. Raaymakers, and L. J. Giling, J. Cryst. Growth,
76, 352 (1986).
6. C. Weber, C. Opdorp, and M. Keijser, J. Appl. Phys., 67, 2109 (1990).
7. D. I. Fotiadis, S. Kieda, and K. F. Jensen, J. Cryst. Growth, 102, 441 (1990).
8. H. K. Moffat and K. F. Jensen, J. Cryst. Growth, 77, 108 (1986).
9. H. K. Moffat and K. F. Jensen, J. Electrochem. Soc., 135, 459 (1988).
10. C. Kleijn and C. J. Hoogendoorn, Chem. Eng. Sci., 46, 321 (1991).
11. M. E. Coltrin, R. J. Kee, and J. A. Miller, J. Electrochem. Soc., 131, 425 (1984).
12. M. E. Coltrin, R. J. Kee, and G. H. Evans, J. Electrochem. Soc., 136, 819 (1989).
13. J. Ouzzani, K. Chiu, and F. Rosenberger, J. Cryst. Growth, 91, 497 (1988).
14. H. Habuka, M. Katayama, M. Shimada, and K. Okuyama, Jpn. J. Appl. Phys., 33,
1977 (1994).
DAT
thermal diffusion coefficient of chemical species A in diluent B (kg
m
Ϫ1 sϪ1
)
ez
g
Ga
Gr
Gx
Unit normal vector along z direction
acceleration due to gravity (m sϪ2
Gay-Lussac number
)
Grashof number
silicon deposition rate at any position x on the wafer surface
(␮m minϪ1
heat transfer coefficient (W mϪ2 KϪ1
diffusion flux of chemical species A driven by concentration gradient
(kg mϪ2 sϪ1
diffusion flux of chemical species A driven by temperature gradient
(kg mϪ2 sϪ1
adsorption rate constant (m sϪ1
reaction rate constant (m4 kmolϪ1 sϪ1
desorption rate constant (m sϪ1
characteristic reactor dimension (m)
)
h
)
jAC
)
jAT
)
Kad
Kf
Kr
L
)
)
15. H. Habuka, T. Nagoya, M. Katayama, M. Shimada, and K. Okuyama, J. Elec-
trochem. Soc., 142, 4272 (1995).
16. H. Habuka, T. Nagoya, M. Mayasumi, M. Katayama, M. Shimada, and K. Okuya-
ma, J. Cryst. Growth, 169, 61 (1996).
)
MwA molecular weight of a gas species A (kg kmolϪ1
)
n
P
Pr
R
RA
Rad
Re
Rsi
Sc
T
unit normal vector
pressure (Pa)
17. CFX 4.1c Solver, AEA Technology, CFX International, Oxfordshire, U.K. (1997).
18. S. V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere Publ. Corp.,
Washington, DC (1983).
Prandtl number
19. H. L. Stone, SIAM J. Num. Anal., 7, 104 (1968).
universal gas constant (J kmolϪ1 KϪ1
)
20. R. D. Lonsdale, Int. J. Num. Meth. Heat Fluid Flow, 3, 553 (1993).
21. ECO Software User’s Guide, Nicolet Instrument Corporation, Madison, WI (1994).
22. P. Ho, A. Balakrishna, M. J. Chacin, and A. Thilderkvist, in Fundamental Gas-
Phase and Surface Chemistry of Vapor-Phase Materials Synthesis, M. D. Allen-
dorf, M. R. Zachariah, L. Mountziaris, and A. H. McDaniel, Editor, Vol. 98-23,
p. 117, The Electrochemical Society Proceeding Series, Pennington, NJ (1998).
23. K. L. Walker, R. E. Jardine, M. A. Ring, and H. E. O’Neal, Int. J. Chem. Kinet., 30,
69 (1998).
24. M. D. Su and H. B. Schlegel, J. Phys. Chem., 97, 9981 (1993).
25. M. T. Swihart and R. W. Carr, J. Phys. Chem., A102, 1542 (1998).
26. P. Ho, M. E. Coltrin, J. S. Binkley, and C. F. Melius, J. Phys. Chem., 89, 4647
(1985).
27. M. A. Mendicino and E. G. Seebauer, J. Electrochem. Soc., 140, 1786 (1993).
28. P. Bratu, K. L. Kompa, and U. Hofer, Chem. Phys. Lett., 251, 1 (1996).
29. S. Bretsznajder, Prediction of Transport and other Physical Properties of Fluids,
Pergamon Press, Oxford (1971).
30. J. O. Hirscfelder, C. F. Curtiss, and R. B. Bird, Molecular Theory of Gases and Liq-
uids, John Wiley & Sons, Inc., New York (1967).
surface chemical reaction rate of chemical species A (kg mϪ2 sϪ1
)
rate of adsorption (kmol mϪ2 sϪ1
Reynolds number
)
silicon deposition rate (kmol mϪ2 sϪ1
Schmidt number
)
temperature (K)
thermal diffusion number
ambient gas temperature (K)
Td
TInf
Twafer wafer surface temperature (K)
Twall wall temperature (K)
U
V
wA
velocity of gas (m sϪ1
)
characteristic velocity of gas (m sϪ1
mass fraction of chemical species A
)
wA,in mass fraction of chemical species A at the reactor inlets
fraction of reactive sites on wafer surface occupied by SiHCl3
X
Downloaded on 2015-06-17 to IP 130.194.20.173 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).