Mathematical modeling of cold-wall channel CVD reactors

Journal of Crystal Growth 94 (1989) 131—144 North-Holland, Amsterdam

131

MATHEMA11CAL MODELING OF COLD-WALL CHANNEL CVD REACTORS W.L. HOLSTEIN and J.L. FITZJOHN Engineering Technology Laboratory, Experimental Station, E. I. du Pont de Nemours & Co. (Inc.), Wilmington, Delaware 19898, USA

and E.J. FAHY, P.W. GILMOUR and E.R. SCHMELZER Central Research and Development Department, Experimental Station, El. du Pont de Nemours & Co. (Inc.), Wilmington, Delaware 19898. USA

Received 20 August 1987; manuscript received in final form 1 June 1988

A two-dimensional numerical model employing the method of finite elements is described for mass transfer-controlled chemical vapor deposition in channel reactors. This model includes all relevant transport processes occurring in the two primary dimensions, including thermal diffusion, buoyancy forces, and the diffusion of energy, mass, and momentum in the axial (flow) direction. The model represents conditions occurring in reactors with large aspect ratios (channel width/height) in the absence of longitudinal roll waves, conditions conducive to the best growth uniformity. The model was used to calculate growth rate uniformity for metalorganic chemical vapor deposition (MOCVD) growth of InP from trimethyl indium. Thermal diffusion was found to be a significant term in the mass transfer equation, decreasing growth rate by about 30% near the front of the susceptor and by lesser amounts further downstream. Growth rate was found to be insensitive to susceptor temperature but dependent on the binary diffusion coefficient and thermal diffusion factor of the reactant in the carrier gas. Use of a tilted susceptor allows for conditions to be established for which there is a region of relatively uniform growth rate. Tilt angle must be chosen to balance between high growth rate and uniformity over large area.

1. Introduction Chemical vapor deposition is a key technology in electronic material processing. Excellent recent review articles by Hess et a!. [1] and Jensen [2] provide a good overview of the technology, with a strong emphasis on the processing aspects and recent developments towards understanding the complex transport processes which occur,

channel reactor (fig. 1), which has a planar top wall. For metalorganic chemical vapor deposition (MOCVD) growth, the susceptor is usually made of graphite, into which the wafer is recessed. The top wall is often cooled with water to provide a more regularly defined temperature distribution within the reactor. This reactor design has usually been referred to as a horizontal reactor. We use the name “channel” reactor to emphasis the pres-

1.1. Channel reactors

Chemical vapor deposition (CVD) in the mass transfer-limited regime is carried out in several

~z0~:

~ gas flow is either parallel or perpendicular to the wafer surface. Included in the former type is the 0022-0248/89/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

WATER FLOW

II

1 FLOW

~

~

Fig. 1. Cold-wall channel reactor.

132

W.L. Holstein et at.

/ Mathematical modeling

of cold-wall channel C VD reactors

1~s are

ence of a planar top wall and the use of the

model. Points on the susceptor surface x~,. x~= x/[cos(O)],

(Ia)

1.2. Susceptor tilt

Ys = x tan(9).

(Ib)

The earliest reactors employing susceptor tilt for metalorganic chemical vapor deposition (MOCVD) [3] used a block with a tilt angle 9 to it (fig. 2a). It was soon realized that this design

For small 9, cos(9)

reactor in orientations other than the horizontal,

defined by

ing h between the susceptor and the cold wall varies with axial position x through the relation h

=

h 0

could lead to the formation of vortices upstream and downstream of the susceptor. More recent reactors have eliminated the possibility of vortices in the front of the susceptor by flushing the leading edge of the susceptor (fig. 2b), while allowing for vortices to remain at the rear of the susceptor. This type of design continues to find wide use. A further refinement involves flushing both the leading and trailing edges of the susceptor (fig. 2c), eliminating all vortices, We will be primarily concerned with the reactor design of fig. 2c. Rectangular coordinates x and y are referenced to an origin at the front of the susceptor: x is axial distance parallel to the top wall and y is the direction perpendicular to the top wall. We neglect the dimension perpendicular to the x and y directions in the two-dimensional

I and x~ x. The cell spac-

— .~

tan(O),

(2)

where h0 is the cell spacing at the front edge of the susceptor (x = 0). 1.3. Transport equations The modeling of fluid flow, heat transfer, and mass transfer in CVD reactors is complex because of the numerous physical phenomena occurring simultaneously, the large temperature differences, and the temperature-dependence of the physical properties. Chemical vapor deposition from dilute concentrations of reactants in a carrier gas (e.g. hydrogen) and in the absence of homogeneous gasphase reactions is described by the following equations: conservation of momentum p(vvv)=

—vp+vT+pg,

(3)

a

conservation of energy

GAS FLOW —+

pç(v.vT)=v.(KVT),

______________

(4)

species specific conservation of mass b

pvVwV _-i5~~\

(~‘~‘~

pD(Vw+wcsV ln T), continuity equation

(5)

(6) c

Fig. 2. Different methods in which susceptor tilt can be introduced: (a) tilted block; (b) flushed leading edge; (c) flushed leading and trailing edges.

where p is density, v is gas velocity, p is pressure, ‘r is the viscous stress tensor, g is the~gravitational acceleration vector, Ci,, is specific heat at constant pressure, T is temperature, K is thermal conductivity, D is reactant-carrier gas binary diffusion coefficient, w is mass fraction of the reactant in the carrier gas, and a is reactant-carrier gas ther-

WL. Holstein et at.

/ Mathematical modeling

ma! diffusion factor. The viscous stress tensor for a Newtonian fluid is T=~t[Vv+ (vv)T]

—

~L(V

v)I,

T

(7)

133

of cold-wall channel CVD reactors

calculated momentum and even thermal boundary

layer thicknesses are often of dimensions comparable to or greater than the spacing above the susceptor.

Correspondingly,

experiments

have

where ~sis viscosity, (vv)T is the transpose of Vv, and I is the identity tensor. In addition, the equation of state for an ideal gas holds for CVD, relating density to pressure, temperature, and mean molecular weight M

shown that the temperature profile across the Fe-

p

(8)

dle and back of the susceptor, as indicated by

where R is the gas constant. Terms for viscous dissipation and reversible compression have not been included in the energy

results presented later. In addition, isothermal

conservation equation, since they are negligible

not play a role, analytic models can provide rea-

for conditions typically encountered in CVD. In

sonable fits to experimental data. However, it

addition, the Dufour term, which reflects energy transfer due to a concentration gradient, has been

should be noted that almost all models contain at least one adjustable constant, usually the diffusion

omitted, since dilute reactant concentrations are involved. This term indeed has been shown to be negligible for typical MOCVD growth conditions [4]. Significantly, the thermal diffusion term in the

coefficient, which is often not known beforehand and is therefore generally chosen to yield the best fit to the data. Bloem and Giling [17] have shown that most analytic models yield growth rate pre-

species-specific mass conservation equation (last

dictions of similar form.

=

Mp/RT,

actor is established fairly rapidly [13]. Sheppard [14] presented a model that assumed fully developed velocity and temperature profiles.

These assumptions more closely reflect actual growth conditions, particularly towards the mid-

models have been presented [15,16]. Under conditions where buoyancy forces do

term on the right-hand side of eq. (5)), which

Berkman et al. [7] were the first to include

reflects mass transfer due to a temperature gradi-

susceptor tilt in an analytic model for chemical

ent, cannot be neglected.

vapor deposition (CVD). This model was based on

1.4. Previous models

the stagnant layer assumption. Chen [18] presented an analytic model based on temperature-in-

The earliest analytic models for mass transfercontrolled CVD include the stagnant layer models of Eversteyn et al. [5], Bloem [6], and Berkman et a!. [7] and the boundary layer models of Grove [8],

dependent physical properties, with the exception of a temperature-dependence to the diffusion coefficient.

models do not correctly describe the physical phe-

Recently, Van de Ven et a!. [16] presented an analytic model for predicting MOCVD growth rate uniformity in the presence of a tilted susceptor for the established flow region corresponding

nomena occuring in channel reactors. The stag-

to locations along the middle and back of the

nant-layer models assume a mixed, or turbulent, layer near the cold wall and a stagnant, or laminar, layer near the hot susceptor. These models were

susceptor. They found an optimal tilt angle 9*

Luznetsov and Belyi [9], and Hitchman [10]. These

suggested by the observation of a clear layer in flow visualization studies using smoke particles, an effect resulting from thermopheresis [11,12], rather than the presence of both laminar and turbulent regimes. Reynolds numbers for CVD indicate that flow is always well within the laminar regime. Boundary layer models are unwarranted, They incorrectly assume a developing velocity profile at the leading edge of the susceptor. Also,

“

—

—

tan~’[2 84 D ‘

0 / TM \ v0h~~ T0)

]‘

(9)

where 1~is a reference temperature (300 K), v0 is feed gas velocity calculated at T0, D0 is reactantcarrier gas binary diffusion coefficient at T0, and the mean gas temperature TM is related to the cold wall temperature T~ and the susceptor temperature 7~by TM = (T~+ 7~)/2. The parameter c is the temperature-dependence of the diffusion

134

W. L. Holstein et at.

/

Mathematical modeling of cold-watt channel CVD reactors

coefficient. For tilt angle 9 = 9*, eq. (9) predicts uniform growth rate. For 9<8*, growth rate is predicted to decrease with axial distance; For

complete to date from the standpoint of predicting the complex fluid flow that can occur from longitudinal roll waves. An exact solution is ob-

9

9* growth rate is predicted to increase with

tamed in the two dimensions normal to the direc-

axial distance. The first two-dimensional numerical model of a

tion of primary gas flow for each grid point in the axial direction and a marching procedure is used to obtain solution in the third (axial) direction. This solution method required the neglect of the

>

channel reactor was presented by Secrest et al. [19] and applied to vapour-phase growth of GaAs

through the chloride process. The first numerical model for CVD in a cold-wall channel reactor was presented by Juza and Cermak [20], who applied their model to the deposition of silicon from SiC1

4. The two-dimensional laminar flow model was simplified by neglecting diffusion of mass, energy,

and momentum in the axial direction. Buoyancy forces were not included and both inertial and viscous effects transverse to the main flow direction were omitted. Transport properties were calculated by assuming a dilute gas mixture of reactant gas in the carrier gas. Their work was the first to numerically show the need to include thermal, or Soret, diffusion in the mass transport

axial diffusion of momentum, mass, and energy as well as the assumption of much greater pressure variations in the axial direction than in the other two directions. As a result, the effect of downstream phenomena on upstream processes are not included. Downstream phenomena can influence upstream flow when buoyancy forces are important in both vertical and horizontal reactor orientations [27]. At low flow rates, axial diffusion of mass is also important [16]. The model has recently been modified to include thermal diffusion and susceptor tilt [28]. Rhee et al. [29] have presented a model in three dimensions.

equation, and the general treatment provided a firm basis for future numerical models. Coltrin et

2. Description of the model

al. [21—23] developed a model with similar assumptions for heat and momentum transport for the study of silicon deposition from SiH4. The

Fluid flow and growth were modeled in two dimensions. The best growth uniformity is

mass transfer part of their model is much more

achieved for reactors of large aspect

complicated. Models that include a third dimension are of interest because of the onset of longitudinal roll waves in horizontal channel reactors at where high 2C~h3(z~T)/~tKT, number Ra = gp gRayleigh is gravitational acceleration, h is the spacing

(width/height) in the absence of longitudinal roll waves, conditions where the flow is largely two-dimensional. Buoyancy and diffusive terms in the axial directions can readily incorporated in two-dimensional models, but be their incorporation

between the hot susceptor and the cold upper wall, and ~T is the temperature difference between the hot susceptor and the cold upper wall. Takahashi et al. [24] presented the first such model in three dimensions. Constant fluid properties were assumed, except for variable gas density in the buoyancy term, which is equivalent to the Boussinesq approximation in the classic Bénard problem [25]. Some progress was made in showing the dependence of longitudinal roll waves on buoyancy forces, and as might be expected, roll waves were predicted to lead to very nonuniform growth. The laminar flow model in three dimensions presented by Moffat and Jensen [26] is the most

ratio

in three-dimensional models requires a large increase in computation time. 2.1. Physical and chemical assumptions in the model A model is only as accurate as the assumptions on which it is based. Too often, the limitations inherent in a model are overlooked. These limitations result both from simplifications in the solution procedure and an incomplete knowledge of the chemical phenomena occurring. Even with present computational capabilities, numerous assumptions are imposed by the need to obtain a solution within a reasonable amount of computation time. In addition to the neglect of the third

W.L. Holstein et at.

/ Mathematical modeling of cold-wall channel CVD

dimension, we assumed dilute gas mixtures for the

GAS FLOW

a

calculation of transport properties and fixed boundary conditions. The dilute gas mixture approximation has been used routinely in CVD models. Gas density, viscosity and thermal conductivity are taken as the values of the carrier gas. Typical mole fractions of Group III organometallics in the feed are 0.0001, and the approximation would be excellent if the organometallics were the only addition to the carncr gas. However, Group V hydrides are usually present in mole fractions of 0.001—0.01, which strongly influences gas density. For a phosphine

135

reactors

h ____________

~ ‘

b

Vt

~Lx

0

L~ L5L~

x Fig. 3. (a) Schematic of reactor used in modeling; (b) temperature profiles used for the lower wall.

mole fraction of 0.005 in hydrogen, the gas density is about 8% greater than for pure hydrogen. For arsine at the same concentration, the density is about 19% greater. Physical properties other than density are only weakly affected by arsine or phosphine at typical mole fractions. The dilute concentration assumption also negates the need to invoke the Stefan—Maxwell equations for diffusion of mass. We consider the MOCVD growth of InP from trimethyl indium and phosphine (CH

3)31n + PH3

=

InP + 3 CH~,

The reaction mechanism is only poorly understood. The rate of MOCVD growth of Ill—V compounds has generally been assumed to be controlled by the rate of diffusion of a group III species to the surface [30].

2.2. Selection of boundary conditions Setting boundary conditions is one of the most critical parts of any model. For the present model, only the boundary conditions for momentum transport can be selected with the certainty that they have no influence on the model predictions. The model is set up with lengths downstream and upstream of the susceptor of sufficient size to ensure that entrance and exit boundary conditions do not influence the model predictions (fig. 3a).

The uncertain boundary conditions are those for T and w on the lower wall and that for T on the upper wall,

The no-slip conditions v~= 0 and v~= 0 apply at the reactor walls. A parabolic profile is assumed for v~ far enough upstream to avoid influence from downstream phenomena. Correspondingly v~, = 0 in this entrance region. These assumptions are excellent. For well-designed reactors, parabolic velocity profiles are established fairly rapidly. The entrance length for the establishment of velocity profiles for gases is roughly equal to 0.04 Re h where the Reynolds number Re = pvh/p [31], corresponding to a distance of 0.5—3 cm for typical CVD growth conditions. The exit boundary conditions are set constant at dv~/dx = 0 and v~.= 0. We assume a temperature profile on the upper and lower walls as shown in fig. 3b. The upper wall temperature was set at T~.The susceptor of length L~ is assumed to be at a constant growth temperature 7~.The lower reactor walls are assumed to be at a constant temperature T~, with the exception of transition zones of length L~ at each end of the susceptor, over which a constant temperature gradient is present. We routinely set L~ to be 2 cm, as suggested by rough experimental measurements. The use of temperature transition zones allows for the avoidance of temperature discontinuities in the boundary conditions, which would otherwise occur at the leading and trailing edges of the susceptor. The entrance gas temperature was set at T~. The temperature boundary condition at the exit of the reactor was selected as dT/dx = 0.

The mass fraction of organometallic in the feed stream entering the reactor is set at w0, which assumes that the feed stream has been uniformly

136

W. L. Holstein et al.

/ Mathematical modeling

of cold-wall channel C VD reactors

mixed prior to entry into the reactor. The boundary condition at the reactor exit was selected as dw/dx = 0. The mass fraction is set at w = 0 on the hot susceptor, reflecting the assumption of complete and rapid decomposition of the organometallic at the susceptor surface. The assumption is made that no decomposition of the organometallic occurs on any of the reactor walls, both upper and lower. Thus the boundary condition for these surfaces is represented by zero flux perpendicular to the surface

effective implementation of a finite element solution strategy. However, the continuity equation does not include second derivative terms and thus requires special treatment. Substituting for p from eq. (8) into eq. (6) for constant p and M yields

~D(Vw

in eq. (3), where A is a large constant. Thus a large “artificial stress” was entered in the momentum

+ ~~vT)

•,,

T

=

0,

V •(v/T) = 0. (10) To preserve total mass conservation, we replace vp by the term vp

=

{

V ~tT~

[v

.(v/T)]

}

(11)

where n is a unit vector normal to the surface, These boundary conditions are somewhat disconcerting, since they impose a discontinuity in growth at the edges of the susceptor in the model solution. Better boundary conditions would include kinetically controlled decomposition and the transition from kinetically controlled to mass transfer-controlled growth with increasing temperature on the bottom wall before and after the susceptor. Incorporation of kinetics is hindered by the absence of kinetic data for the decomposition of the organometallics on materials of construction for MOCVD reactors.

equation for any velocity field not satisfying the continuity equation (10). The finite-element method then ensures conservation of mass, ~ ~ /7’~ 0 ~l2

2.3. Method of solution

mechanics [32]. The key dimensionless numbers influencing the solution are the Reynolds number Re = pvh/p., the local Peclet number for heat transfer PeT (local) = vpC~z.~x/K,and the local Peclet number for mass transfer of the organometallic in the carrier gas PeM (local) = v~x/D where Z~xis the grid dimension. Re reflects the relative magnitudes of inertial and viscous forces and is less than 50 for typical growth conditions, indicating that flow is well within the laminar regime. The grid spacings are made small enough to ensure that PeT(local) < 20 and PeM(local) < 100 calculated at the mean gas temperature to ensure that upwinding was not needed for handling the convection terms in the finite element solution. The energy, mass, and momentum equations are all coupled; that is, the solution of one is affected by the results of the others. However, we note that, for dilute concentrations of organome-

Fluid flow, temperature and concentration profiles, and growth rates in the MOCVD reactor were modeled by the finite-element technique [32—34]. All of the terms in eqs. (1)—(6) were included. The solution domain was divided into triangular elements with nodes at corners and midpoints of the line segments. Gas density p can be eliminated by substitution of eq. (8), leaving five unknown variables: v~, vp.. T, w, and p. Temperature, velocity, and mass fraction were represented by biquadratic polynomials. Absolute pressure can be considered constant, but small variations in pressure occur that affect gas flow, We allow for variable pressure and treat it specifically, as discussed below. The equations of conservation of energy, momentum, and species-specific mass have the general second-order derivative terms needed for

“

/

/

in its solution of the momentum equation. For our 7/K suitable. This application, we found A = 10 technique, sometimes referred to as the penalty function method [33], couples the continuity and momentum equations. When incorporated in this form into a finite-element solution strategy, it restricts the computer search to solutions satisfying eq. (12). This method is widely used in the implementation of finite-element analysis to fluid

W. L. Holstein et aL

/ Mathematical modeling of cold-wall channel

tallics in the carrier gas, the energy and momen-

137

C VD reactors

Insufficient experimental data were available to

tum conservation equations are not affected by

calculate the binary diffusion coefficient and ther-

solution of the mass conservation equation. Thus, we can decouple the mass conservation equation and use a two-step, or serial, procedure where the coupled energy and momentum conservation

ma! diffusion factor. Lennard-Jones parameters were estimated as described previously [35]. The binary diffusion coefficient was estimated from the Lennard-Jones parameters through the standard process [36]. The binary diffusion coefficient can be expressed in the form

equations are first solved, and the resulting temperature and flow fields are then used to solve the mass conservation equation. This procedure has been used previously [26]. Simultaneously with the solution of the species specific mass equation, thefrom growth rate across theconservation susceptor was calculated the mass flux

sur face

I

j

of the organometallic to the susceptor aw

=

PD(VW +

~vT)

.n.

(13)

Growth rate G was then calculated from the relations p G= (Ms/MRpS)j,

(14)

D= D

0(P0/P)(T/T0~.

(18)

and7~ c== 300 1.775. For trimethyhndium—hydrogen 2/s with K and P0 = 1 . atm,. D0 = 0.330 cm Thermal diffusion factors were estimated from .

0

‘

-

Lennard-Jones parameters via the simplified relationship for dilute gas mixture derived previously [35]. For dilute concentrations of trimethyl indium in hydrogen, a varies from 0.67 at 50°C to 1.29 at

650°C.

For MOCVD growth of InP, M~= 149.8 g/mol and p 3. The molecular weight of the 5 = 4.79 g/cm limiting reactant trimethyl indium MR is 159.9 g/mol.

where M 5 is molecular weight of the solid, MR 15 molecular weight of the controlling reactant, and Ps

is density of the solid,

For the purpose of calculating growth rates, we assume a feed mole fraction of trimethyl indium of 1.00 X i0~, corresponding to a feed mass fraction w0 of 7.87 x i0~. Growth rate increases linearly with w0.

2.4. Calculation of physical properties

The approximation of dilute concentrations of reactants in the carrier gas was employed. Den-

3. Results

sity, viscosity, and thermal conductivity were assumed to be those of the pure carrier gas. They can be expressed as a function of temperature and pressure by

3.1. Thermal diffusion

p

=

=

p0(1~/T)(p/p0),

(15)

.t~~(T/T~~)

(16)

K= K0(T/1~)”.

(17)

For hydrogen with T0 = 27°C3,(300 and Po = 1 p~K) = 8.89 X i0~ atm, .p0s, = a 8.19 x iO~ g/cm = 0.66, K g/cm 0 = 17.4 J/cm ‘S. K, and b = 0.68. For hydrogen, C, varies from 14.3 J/g’ K at 50°C to 14.8 J/g’ K at 650°C and is independent of pressure over the range of pressures of interest for CVD. It was assumed to be 14.5 J/g’ K for all conditions,

We consider the case of a channel reactor in the normal horizontal orientation for the growth of indium phosphide from trimethyl indium in hydrogen carrier gas at atmospheric pressure. The spacing above the reactor h is selected to be 2.0 cm. The susceptor length is selected to be 10.0 cm, except where noted, with side. 2 cm The temperature tion zones L~ on either susceptortransitemperature is selected as 650°C and the entrance and cold-wall temperatures as 500 C. The flow rate per unit width of the reactor is selected as 0.833 liter (STP)/min ‘ cm, corresponding to a flow of 10 liter (STP)/min in a reactor of 12 cm width. We have shown experimentally that the 2.0

138

W.L. Holstein ci at.

/

Mathematical modeling of cold-walt channel CVD reactors

cm spacing above the susceptor is such that longitudinal roll waves do not occur for these conditions [12]. Other effects from buoyancy forces are also minimal for these conditions [27].

for The 4b, the respectively. above conditions calculated velocity The temperature areand shown temperature inprofiles figs. 4afields and are established rapidly. They are fully established fairly for distances greater thanalmost a thermal entrance length LT [271 of about 2.4 cm from the leading edge of the susceptor. Mass fraction pro-

files for the organometallic normalized to the feed

0.8

E

0

0.6 ~1RM 0.4 DIFFUSION ~

OMITTED

0.2

mass fraction, w/w

0, were calculated for both the inclusion of thermal diffusion (fig. 4c) and its omission (fig. 4d). As expected, with thermal dif fusion included, the mass fraction, of the organometallic near the cold wall increases, even to the extent that it is nearly 40% greater than in theCalculated feed. growth rates along the susceptor are shown in fig. 5 for both the inclusion and omission of thermal diffusion. The neglect of thermal diffusion results in growth rates ranging from about 30% too high at 1 cm from the leading edge of the susceptor to 25% too high at 1 cm from the trailing edge. These results agree with previous results [4,6,20,35,37] indicating that thermal diffu-

-

.

.

-.

a

-~

~

~‘

I

-_______________________________ 2

4

6 X

8

(cm)

Fig. growth ratenoacross the susceptor for = and 0.1 in cm, a reactor susceptor tilt for inclusion atm, 5.h0Calculated =32.0 Q/s =with 0.833 liter(STP)/min.cm, andp w0 = 7.87X10 omission of thermal diffusion.

sion cannot be neglected in modeling CVD. In the remaining discussion, we consider only the results where thermal diffusion was included. 3.2. Sensitivity of results to physical properties and boundary conditions The effects of boundary conditions and finite

-

~

7

0

SUSCEPTOR 50°C

CEPTOR

________________________________________ ~

element grid mesh were examined. Increasing the grid mesh was found to have no effect on the growth rate, except for distances less than 1 cm from the leading and trailing edges of the susceptor. We consider results at the leading and trailing edges unmeaningful an~ay because of the discontinuity in growth imposed by the boundary conditions. Changing L~ over the range of 0—3 cm again influenced growthedges rates only 1 cm of the leading and trailing of thewithin susceptor.

CEPTOR

Decreasing the susceptor temperature from 650 to 6 ____—

d SIJSCEPTOR I Fig. 4. Calculated fields in reactor without suscepior tilt for p = 0.1 atm, h 0 = 3: 2.0 velocity cm, Q/s(a),= 0.883 temperature .liter(STP)/mincm, (b), and normalized fraction of organometallic w/w and w0 mass = 7.87X10 0 for inclusion (c) and omission (d) of effects resulting from thermal diffusion.

600°C resulted in less than a 0.2% decrease in growth rate. The effect of the diffusion coefficient was also studied. As noted, diffusion coefficients must be estimated at present. The values used here are probably accurate to within 20%. For the conditions modeled, a 20% increase in diffusion coeffi-

W.L. Holstein er aL

/ Mathematical modeling

cient produces about an 11% increase in growth rate over the entire axial distance. Thermal diffusion factors are also not known precisely. A 20% increase in thermal diffusion factor decreases growth rate from 5.7% at 1 cm from the leading edge of the susceptor to 5.2% at 1 cm from the trailing edge. We observe that uncertainty in mass transfer properties influences the predictions of this and other models, while uncertainty in the temperature transition zones at the leading and trailing edges of the susceptor have only a small effect. We have not examined the uncertainty originating from lack of complete knowledge of the reactant decomposition kinetics. We conclude, however, that enough uncertainty exists to prevent exact predictions of absolute growth rates a priori. Neverthe-

139

of cold-wall channel CVD reactors

o.E

0.6

100 80

C 0.4 00

0.2

0 x (cm)

Fig. 6. Calculated growth rate across the susceptor for p

atm, h3,and 8=0°—8°(L 0 = 2.0 cm, Q/s = 0.883 liter(STP)/min’cm, X10

5=lOcm)and

w0

=

=

0.1 7.87

100 (Ls=7.Scm).

less, we feel models such as this give excellent

predictions of growth-rate uniformity. As noted, a number of analytic models have been presented, including isothermal models. We used the model to examine the isothermal mass transfer problem by setting T~ and T~ to the mean temperature of 350 °C. By using the diffusion coefficient as an adjustable parameter, we were able to fit the predictions of the full model to within a few percent. Thus it is not surprising that analytic models with an adjustable diffusion coefficient have found some success in fitting curves of growth rate versus axial position, despite their obvious oversimplification of the complex transport processes. In the presence of moderate or strong buoyancy forces [27], simple analytic models fail.

than can also yield regions of uniform growth rate. For typical tilt angles, the region of uniform growth rate increases as tilt angle decreases. But this is accompanied by a decrease in growth rate. These two factors must be balanced in choosing a tilt angle for use. Velocity, temperature, and concentration profiles for the case of 6° tilt are shown in figs. 7a—7c. Due to the decrease in channel height, the velocity continues to increase with axial distance (fig. 7a), in contrast to the case of no tilt (fig. 4a), where a constant velocity profile is established for distances greater than a thermal entrance length 50

______________________________________________

3.3. Susceptor tilt

We consider first growth rate uniformity as a function of tilt angle. Growth rate is plotted as a function of axial position in fig. 6 for susceptor tilt angles of 00_100. Tilt angles greater than about 50 resulted in a minimum in growth rate with axial position. The location of the minimum, and correspondingly, the optimal location for placement of the wafer, are dependent on the tilt angle. the tiltrate angle increases, the For region of uniformAsgrowth moves upstream. longer susceptors than shown in fig. 6, tilt angles less

3, = and Fig. 7. Calculated fields in the reactor for p = 0.1 atm, h0 2.0 cm, Q/s=0.833 liter(STP)/min’cm, 9 = 60: velocity (a), temperature (b), w0=7.87X10 and normalized mass fraction of organometallic w/w 0 (c).

140

W.L. Holstein ci at.

/ Mathematical modeling of cold-wall channel

of about 2.4 cm. For the tilted susceptor reactor, uniform temperature gradients above the susceptor are not established (fig. 7b). This result contrasts with the case of no tilt, where the ternperature gradients above the susceptor are constant for distances greater than about LT (fig. 4b). We discuss further the implications of this result with respect to growth uniformity below. LT

a

CVD reactors

N~

h~

I’ 0 Tw b

T

h

4. Discussion

Vs C

h~

4.1. Thermal diffusion y

The effect of thermal diffusion in cold wall reactors is to decrease the overall mass transfer rate to the hot susceptor surface. For typical reactor geometries and operating conditions with hydrogen or helium carrier gases, thermal diffusion results in decreased growth rate over the entire susceptor with respect to that expected in the absence of thermal diffusion. However, for very high reactant conversions, a point is reached downstream in the reactor where the growth rate is actually higher due to thermal diffusion [37], a result of the decreased consumption of reactant upstream. In all cases, thermal diffusion results in decreased total deposition rate integrated over the length of the susceptor. 4.2. Established flow with no tilt

Numerical models provide insight into the transport processes and the sensitivity of growth rates to process conditions. We have already noted the significance of thermal diffusion. The model confirms that temperature profiles are established fairly rapidly (fig. 4b). Velocity profiles are also established rapidly (fig. 4a). Towards the back of the susceptor, the temperature and velocity profiles are fully established. A comparison between the isothermal entrance region and the fully established region over the middle and rear of the susceptor is shown in fig. 8. Temperature T is constant (50°C) in the entrance region. The profile of velocity v~ is parabolic, which also results in a parabolic profile for total mass flow pv~..In the established flow region over

~

—~

eva Fig. 8. Comparison of profiles for temperature (a). velocity v~. (b), and mass flow pv~ (c) at the reactor entrance region (— — —) and in the established flow region ( ).

the middle and rear of the susceptor, the temperature varies from 50°C at the cold wall to 650°C at the hot wall. The temperature gradient is nearly constant, deviating slightly because of the temperature dependence of the thermal conductivity. The velocity profile is nearly parabolic, deviating slightly because of the temperature dependence of viscosity, but is greater than for the entrance region by nearly a factor of two owing to gas expansion. The total mass flow profile is greatly altered by the inverse-temperature dependence of density, eq. (8). Thus, as the gas flows through the reactor, most of it flows through near the cold wall at relatively low temperature. Noting further from fig. 4c that the mass fraction of the organometallic is smallest where the temperature is highest, we conclude that as the organometallic flows through the reactor, only a very small fraction of it is being heated to high temperature at any one time. 4.3. Optimal tilt angle

The present results on susceptor tilt are fundamentally different from the predictions of Van de Ven et al.’s analytic model [16], eq. (9). For the conditions analyzed here (h0 2.0 cm, 7~ 300 =

=

W. L Holstein et aL

__________

,~-WAFER

/ Mathematical modeling

_______

l[sUSCEPTORi~

Fig. 9. Reactor with variable cell spacing for growth rate uniformity over large area.

2/s, c = 1.775),

K, v0 7.6 cm/s, D0 = 0.330 cm eq. (9) yields an optimum tilt angle 9* of 6.2°. Our results differ in that we conclude that a single optimal tilt angle cannot be defined for growth =

rate uniformity alone. For any tilt angle and for a susceptor of there properis length, a position canrate be found where a minimum in growth uniformity with position, and correspondingly a region of uniform growth. For lower tilt angles,

141

of cold-wall channel CVD reactors

susceptor surface below, assuming constant susceptor temperature, in order to assess the magnitude of the variability in heat flux. The heat flux from the susceptor surface q is equal to the combined heat fluxes by conduction ~ and radiation q~, q=q~+q~. (19) Ignoring radiation to the susceptor, the heat flux by radiation is q~ ea’T~, (20) where T 5 is susceptor temperature, e is emissivity and the Stefan—Boltzmann constant a’ flux 5.67 by X 2~K4. For small 8, the heat 1012 J/s~cm conduction is =

=

q~ —K(dT/dy)~... =

(21) where K is thermal conductivity of the gas, and the temperature gradient (dT/dy) is evaluated at the susceptor surface y5, defined by eq. (ib). The thermal entrance length LT [27], the axial distance from the leading edge of the susceptor where convective effects influence heat transfer, is about 2.4 cm for the conditions analyzed here. For distances greater than about LT, heat transfer occurs only by conduction and radiation. For 0° tilt, q~is greatest at x 0 and decreases through the thermal entrance region, finally becoming about constant for positions on the susceptor greater than LT. The temperature gradient above the wafer is about constant for x> LT, and correspondingly, the total heat flux from the susceptor q is about constant. For a uniformly heated susceptor, we expect the wafer temperature to be about constant for x > LT. For the case of susceptor tilt, heat flux from the susceptor is not uniform for x > LT. Rather, heat flux increases with axial distance due to the increase in thermal conduction resulting from the decrease in cell spacing h between the hot suscep55,

the location of the minimum occurs further downstream, and the axial length over which uniform growth occurs is larger. However, the growth rate is reduced for constant reactant feed gas concentration. Thus there is no single optimal tilt angle from the standpoint of growth rate uniformity alone. Tilt angle must be chosen to balance between high growth rate and growth rate uniformity over large axial distance, We have considered the case of a tilted susceptor. For small tilt angles, the results also apply to reactors with tapered top walls [16]. For both the traditional tilted susceptor [3] and tapered top wall [16] designs, the constraint of planarity of the susceptor and top wall ultimately limits the reactor performance. Improved growth rate uniformity over large surface areas can be obtained by forming the susceptor or top wall to a preferred non-planar shape. Since the wafers are planar, this is most easily accomplished by varying the shape of the top wall (fig. 9). Numerical models such as the one used here can help to determine the ideal shape of such a reactor.

=

htorcan lead to problems in controlling thefrom susceptemperature due to variable heat flux the

tor and the cold wall. We can estimate the change in q with axial position. The heat loss due to radiation is constant with axial position. We take e to 2.bes. about 0.6.LT, ForweT5can650 °C, q~ 2.46 For x> approximate eq. J/cmby (21)

susceptor surface with axial position. We present approximate calculations of heat flux from the

qc

4.4. Heat flux from susceptor

The use of a reactor with a variable cell spacing

=

=

K(T~— Tw)/h,

=

(22)

142

W.L Holstein et aL

/ Mathematical modeling

where T~ is susceptor temperature, T~ is wall temperature, and thermal conductivity evaluated at the mean temperature TM. For

cold

K is T 5= 50°C and hydrogen carrier gas, K

650°C, T~ 2.85 x iO~ J/cm’ s’ K at the mean temperature of 350°C. We compare two positions on the susceptor for 6° tilt in the region where good growth rate uniformity is observed in fig. 7, x 4 cm and x 8 cm. From eq. (2), h 1.58 cm at x 4 cm and h 1.16 8 cm. At xs; at4 2 s cm andatq x 3.54 J/cm2 xcm, 8q~ cm, 1.08 q~ J/cm 1.47 J/cm2. s and q 3.93 J/cm2 s. The heat flux due to thermal conduction is about 36% greater at x 8 cm than at x 4 cm. Correspondingly, total heat flux is about 11% greater at x 8 cm than at x 4 cm. The increased heat flux can lead to a decrease in susceptor temperature with axial distance for x> LT, in contrast to 0° tilt, where heat flux is constant for x > LT and uniform temperature is expected. In making the calculations with the model, we have assumed that the susceptor temperature was constant. In actuality, variable heat flux can lead to variations in susceptor temperature. We have shown above that changes in susceptor temperature have little effect on growth rate, but when doping levels or composition are dependent on growth temperature, tilted susceptors might be more prone to problems in dopant or compositional uniformity. =

=

=

=

=

=

=

=

=

‘

=

=

=

=

addition, the results of numerical models provide insight into the physical processes that occur. There is no single susceptor tilt angle which will lead to optimal growth rate uniformity in CVD reactors. Tilt angle must be chosen to provide a balance between high growth rate and growth rate uniformity over large area. Computer-aided design of pew reactor geometries can yield reactors exhibiting higher reactant utilization efficiencies than possible with susceptor tilt alone.

=

=

=

‘

of cold-wall channel CVD reactors

Acknowledgments

=

=

This work was supported by BT & D Technologies, Ltd. We thank the Directors of Research and Development of BT&D Technologies for permission to publish this paper. We are grateful to R.H. Moss and A.D. Nelson for helpful discussions throughout the course of this work. We thank T.R. Anderson for insightful discussions on tilted susceptor reactors. The assistance of A.L. Squyres in development of the computer graphics is also gratefully acknowledged.

List of symbols a b

5. Conclusions

c

Analytic models provide a conceptual understanding of the transport processes but cannot realistically incorporate all of the complex phenomena occurring simultaneously. Numerical models provide a more complete representation of the transport processes. Even with the best of present numerical models, there is too much uncertainty in boundary conditions, growth chemistry, and the values of physical constants to exactly predict growth rates a priori. These limita-

Ci,, D D e0

tions are much less serious for predicting relative growth rate with position or process conditions. Suchgrowth modelsrate are uniformity of considerable value in predicting and reproducibility. In

g g

Exponent in temperature dependence of viscosity Exponent in temperature dependence of thermal conductivity Exponent in temperature dependence of diffusion coefficient Specific heat (J/g’ K) Diffusion coefficient (cm2/s) 2/s) Diffusion Emissivitycoefficient at 1~(cm Gravitational acceleration (980 cm/s2) Gravitational acceleration vector (cm/ ~2)

G h h 0 I

jK

Growth rate (nm/s) Cell spacing (cm) Cell spacing at leading edge of susceptor (cm) Identity tensor 2’ s) Mass fluxconductivity (g/cm Thermal (J/cm. s~K)

W.L Holstein et aL / Mathematical modeling of cold-wall channel C VD reactors

L

5 LT, L M MR

143

Susceptor length (cm) Thermal entrance lengthlength (cm) for botTemperature transition torn wall (cm) Molecular weight of gas (g/mol) Molecular weight of limiting reactant (g/mol)

a’ ‘r

Stefan—Boltzmann 2~K4) constant (5.67 x 10- 12 J/sstress -cm tensor (g/cm2. s) Viscous

Molecular weight of solid (g/rnol) Pressure (atm) Reference (1 atm) Total heat pressure flux (J/cm2. s) Heat flux by conduction (J/cm2 s) Heat flux by radiation (J/cm2. s) Volumeric flow rate at standard temperature and pressure (liter(STP)/min) Gas constant (82.06 atm crn3/mol. K) Reactor width (cm) Absolute temperature (K)

PeT (local) Local Peclet number for heat transfer

Dimensionless numbers

PeM(local)Local Peclet number for mass transfer (= vL~x/D)

M 5 p p0 q

q~ q~ Q

R s T

(

Ra Re

vpC,,z.lx/K)

2C~h3(L~T)/~iK)

Rayleigh number Reynolds number(= (=gp pvh/~t)

‘

References [1] D.W. Hess, K.F. Jensen and T.J. Anderson, Rev. Chem. Eng. 3 (1985) 97. 121 K.F. Jensen, Chem. Eng. Sci. 42 (1987) 923. [3] SJ. Bass, J. Crystal Growth 31(1975)172.

T 0

T~ T5

v v

w

w0 x x~ I.~x y

Reference temperature (300 K) Cold wall and entrance temperature (K) Susceptor temperature (K) Gas velocity (cm/s) Gas velocity vector (cm/s) Gas velocity at T0 (cm/s) Gas velocity in axial direction (cm/s) Gas velocity in height direction (cm/s) Mass fraction of reactant Mass fraction of reactant in feed gas Axial distance (cm) Value of x on susceptor surface (cm) Grid spacing (cm) Distance in direction perpendicular to top wall (cm) Value of y on susceptor surface (cm)

c/k 0

9* A p p5 a0

and R. Pollard, J. Electrochem. Soc. 131

(1984) 2917. [5] F.C. Eversteyn, P.J.W. Severin, C.H.J. van den Brekel, and H.L. Peek, J. Electrochem. Soc. 117 (1970) 925. [61J. Bloem, J. Electrochem. Soc. 117 (1970) 1397. [7] S. Berkman, VS. Ban and N. Goldsmith, in: Heteroepitaxial Semiconductors for Electronic Devices, Eds. G.W. Cullen 264281.and CC. Wang, (Springer, New York, 1978) pp. [8] A.S. Grove, Physics and Technology of Semiconductor

Thermal diffusion factor Lennard Jones energy parameter (K) Reactant utilization efficiency Susceptor tilt angle (deg) Defined by eq. (1) (deg)

Devices (Wiley, New York, 1967). [9] F.A. Luznetsov and VI. Belyi, J. Electrochem. Soc. 117 (1987) 785. [10] M.L. Hitchman, J. Crystal Growth 48 (1980) 394. [11] L.J. Giling, in: Crystal Growth of Electronic Materials, Ed. E. Kaldis (Elsevier, Amsterdam, 1985) pp. 71—85. [12] W.L. Holstein and J.L. Fitzjohn, J. Electrochem. Soc., submitted. [13] L.J. Giling, J. Electrochem. Soc. 129 (1982) 634. [14] W.H. Shepherd, J. Electrochem. Soc. 112 (1965) 988. [15] F.C. Rundle, Intern. J. Electron 24 (1968) 405. [16] J. van de Ven, G.M.J. Rutten, M.J. Raaijmakers and L.J. Giling, J. Crystal Growth 76 (1986) 352. [17] J. Bloem and L.J. Giling, in: Current Topics in Materials Science, Vol. 1, Ed. E. Kaldis (North-Holland, Amsterdam, 1977) pp. 214—342. [18] K. Chen, J. Crystal Growth 70 (1984) 64. [19] B.G. Secrest, W.W. Boyd and D.W. Shaw. J. Crystal Growth 10(1971)251.

Penalty function constant (1/K) Viscosity (g/cm s) 3) Gas density (g/cm 3) Solid density (g/cm 3) Gas density at 7~and (g/cm (A) Lennard-Jones distancePo parameter

[20] J. Juza and J. Cermak. J. Electrochem. Soc. 129 (1982) 1627. Coltnn, R.J. Kee and J.A. Miller. in: Proc. 9th [21] M.E. Intern. Conf. on Chemical Vapor Deposition, Eds. M. Robinson, C.H.J. van den Brekel, G.W. Cullen and J.M. Blocher, Jr. (Electrochemical Society, Penningion, NJ, 1984) pp. 31—43.

Greek symbols

a

[41J.P. Jenlunson

144

W.L Holstein et aL

/

Mathematical modeling of cold-wall channel CVD reactors

122] ME. Coltrin, R.J. Kee and J.A. Miller, J. Electrochem. Soc. 131 (1984) 425. [23] ME. Coltrin, R.J. Kee and J.A. Miller, J. Electrochem. Soc. 133 (1986) 1207. [24] R. Takahashi, Y. Koga and K. Sugawara, J. Electrochem. Soc. 119 (1972) 1406. [25] S. Chandrasekhar, Hydrodynamics and Hydromagnetic Stability (Clarendon, London, 1961). [26] H. Moffat and K.F. Jensen, J. Crystal Growth 77 (1986) 108. [27] W.L. Holstein and J.L. Fitzjohn, J. Crystal Growth 94 (1989) 145. [28] H.K. Moffat and K.F. Jensen, J. Electrochem. Soc. 135 (1988) 459. (29] 5. Rhee, J. Szekely and O.J. Ilegbusi, J. Electrochem. Soc. 134 (1987) 2552. (30] D.H. Reep and S.K. Ghandi, J. Electrochem. Soc. 130 (1983) 675.

[31] H. Schlichting, Boundary-Layer Theory (McGraw-Hill, New York, 1968). [32] AJ. Baker, Finite Element Computational Fluid Mechaflies (McGraw-Hill, New York, 1983). [33] O.C. Zienkiewicz and Y.K. Cheung, The Finite Element Method in Structural and Continuum Mechanics (McGraw-Hill, London, 1967). [34] G.E. Myers, Analytical Methods in Conduction Heat Transfer (McGraw.Hill, New York, 1971). [35] W.L. Holstein, J. Electrochem. Soc. 135 (1988) 1788. [36] R.C. Reid and T.K. Sherwood, The Properties of Gases and Liquids (McGraw-Hill, New York, 1966). [37] J. Ouazanni, K-C. Chiu and F. Rosenberger, in: Proc. 10th Intern. Conf. on Chemical Vapor Deposition, Ed. G.W. Cullen (Electrochemical Society, Pennington, NJ, 1987) pp. 165—174.