Plasma enhanced chemical vapor deposition modeling

Sur/ace and Coatings Technology, 49 (1991) 387—393

387

Plasma enhanced chemical vapor deposition modeling E. Hyman, K. Tsang, I. Lottati and A. Drobot Science Applications International corporation, 1710 Goodridge Drive, McLean, VA 22102 (U.S.A.)

B.

Lane*, R. Post and H. Sawint

Applied Science and Technology, 35 cabot Road, Wohurn, MA 01801 (U.S.A.)

Abstract We are developing a model to simulate the plasma enhanced chemical vapor deposition (PECVD) of thin diamond films. The emphasis to date has been on the development of stand-alone modules to simulate the microwave-induced time-dependent electric and magnetic fields, the generation and energization of plasma electrons in the discharge, the non-equilibrium hydrocarbon chemistry, and the development of a two-dimensional unstructured mesh hydrodynamics solver capable of simulating flow through geometrically realistic reactors. The coupling of the various modules, and the incorporation of a turface chemistry module for the substrate deposition, into a self consistent reactor model is underway. We present some preliminary results from components of a model 2.45 GHz microwave reactor employing H, with 1% CH 4 and operating at a gas pressure of 5.3 x l0~Pa (40 Torr). We have completed an electromagnetic model of the microwave energy deposition in the plasma and calculated the field patterns in the reactor. We have also performed point calculations of the time-dependent electron distribution and of the build-up of atomic hydrogen, the gas temperature, and the resulting generation of CH1, C, H2, and other hydrocarbon radicals. We have also completed a fluid simulation of the flow through the reactor using unstructured mesh techniques. The results we discuss in this paper indicate that careful treatment of non-equilibrium processes in PECVD reactors as well as accurate representation of reactor geometry are essential to a useful simulation capability.

I. Introduction The ability to deposit thin diamond films rapidly onto substrates with a high degree of uniformity using the plasma enhanced chemical vapor deposition (PECVD) technique is a high priority technology goal. It is generally recognized that an improved understanding of the microscopic mechanisms in PECVD reactors and of the sensitivity of the various reactor parameters is needed. The important design issues for PECVD reactors are as follows: efficient coupling of microwave energy to the plasma and to the process gas; efficient transport of activated process gas to the wafer or substrate deposition area; efficient use of the injected gas; uniformity of chemically active species flux across the deposition area. It is desirable to avoid reactor designs that have the following: high microwave electric fields in regions away from the desired plasma formation location, leading to plasma discharge near chamber walls or breakdown of dielectric materials; flow patterns which carry activated species to the reactor walls

*permanent address: MIT, Plasma Fusion Center, Cambridge, MA 02139, U.S.A. ~Permanent address: MIT, Chemical and Electrical Engineering Department, Cambridge, MA 02139. U.S.A.

0257-8972/91/S3.50

or out through pumping ports rather than to the wafer; stagnant or circulating flow patterns above the wafer, buffering the wafer from the desired chemically active species. To understand these issues and to provide input to improved reactor designs we are developing a self-consistent numerical model which simulates each of the essential mechanisms in the PECVD reactor. A physically realistic model requires careful simulation of the electromagnetics, the plasma physics, the neutral gas flow, and the homogeneous and heterogeneous chemistry. Furthermore, the different elementary processes in the reactor are highly interactive; for this reason it is difficult to foresee intuitively the impact of varying one or another reactor parameter. For example, the microwave source induces a complex, geometrically dependent and time varying electric field which ionizes the gas; the resultant build-up of electrons alters the developing electric field distribution. The microwave energy input heats the electrons, and the energetic part of the non-Maxwellian electron energy distribution dissociates the gas, inducing a rise in the gas temperature. The amount of dissociation and heating depends sensitively on the high energy tail of the electron distribution which consequently must be accurately determined. The interaction of the neutral gas with the plasma alters the molecular input stream of H2 to include a substantial

~ 1991

—

Elsevier Sequoia, Lausanne

388

E. Hyman ci a!. / Plasma enhanced CVD modeling

I

affectsThe the hydrocarbon ionization rate and the iselectron distribution. chemistry non-equilibrium and both the flux and the spatial distribution of ap-

~[~7) ~ ~ k.~ v~1

propriate radicals reaching the wafer are very sensitive to the geometrical configuration of the reactor and to the details of the flow configuration through the reactor. We have focused to date on the development of modules to simulate the microwave-induced time-dependent electric and magnetic fields, the generation and energization of plasma electrons in the discharge, the evolution of the molecular and atomic hydrogen gas, the non-equilibrium hydrocarbon chemistry, and the development of a two-dimensional unstructured mesh hydrodynamics solver capable of simulating flow through geometrically realistic reactors. The coupling of these modules, and the incorporation of a surface chemistry module for the substrate deposition, into a self consistent reactor model is underway. In the next section we describe in some detail the generation of the microwave field and the transfer of the field energy to the electrons. This is followed by preliminary model results and our conclusions,

2. Microwave field and plasma generation The absorption of microwaves and the creation of the plasma which transfers energy to the neutral species in the reactor involves the solution of two closely coupled problems. They are (1) the determination the electromagnetic in the cornplex of geometry of the reactorfield and patterns (2) the formation of the electron distribution function. At a pressure of 5.3 x lO~Pa (40 Torr) and a gas temperature in the plasma region greater than 2000 K, the mean free path of an electron with neutrals is approximately 5 x l0~m. During the time an electron gains the average electron energy (approximately 2 eV) typical of the reactors we are modeling, it undergoes around 150 collisions and has a mean displacement of approximately 7 x l0~m. Thus, to an excellent approximation, the heating of the electrons results from the microwave electric fields which are local to the electron’s spatial location. The electron distribution function satisfies the Boltzmann equation. Because an electron undergoes many collisions as it is heated, the distribution function is nearly isotropic and can be well approximated by the zero and first order terms of a spherical expansion, the latter representing a distortion of the distribution function in the direction of the applied field, oscillating at the microwave frequency w. The equation for the electron distribution function is

1

(

component of atomic hydrogen, and this, in turn,

~

v~, t~F 0 2+ w2 -~i;-

2

,.

+ v2V

(-‘—

VF

I

/

(

0 + vV~ V~—~ J \~ ~v 2m I 3VmFo) ~

\Ym

=Lj+Lx~~~~(V

where E 0 is the amplitude of the electric field, e, m~, and v are the electron charge, mass, and velocity respectively, vrn is the electron momentum transfer frequency, F0 is the zero order approximation to the distribution function, V is the bulk fluid velocity, M is the neutral mass, and L1 and L~are the inelastic loss terms that affect the distribution function respectively via ionization and excitation of rotational, vibrational, and electronic levels. The first term on the left represents the electron velocity diffusion due to the cumulative affect of many small angle scatterings of the electron induced by the oscillating electric field. The next two terms give the affect of the divergence of the diffusive and convective fluxes respectively. The last term on the right gives the energy loss due to elastic collisions. Below a critical electric field the ionization rate is exceedingly small and hence the electron density and power deposited per unit volume are also small. Above the critical field the ionization rate and power deposition increase rapidly. If the power deposition is kept approximately constant and equal to the power injected into the reactor the electric field will rapidly adjust to a level close to the breakdown value. power 2fle,The where ~e is dethe posited per unit volume scales as E0 electron density; as the electron density rises the electric field drops. The above considerations provide the necessary prescription for determining the time-dependent evolution of the electric field, the electron density, and the electron energy distribution. Using a set of elastic and inelastic cross-sections, the last two parameters define the time-dependent evolution of the fluid, including the build-up of atomic hydrogen and the rise in the gas temperature as the gas dissociates. In addition to H 2 and H, the Boltzmann calculation monitors the evolution of H,~, H3~,H~,and H— and separately tracks each of the three lowest vibrational levels of H2. The hydrocarbon chemistry is initiated by energetic electrons, but being trace constituents the hydrocarbons do not significantly affect the electron development. It is useful to take advantage of the separation of time scales inherent in this problem. The electron distribution function is established on a time scale of around l0_8 s, the electron density growth occurs over approximately 1 us, and the hydrogen dissociation and hydrocarbon chemistry as well as fluid convection and diffusion occur on a millisecond time scale.

389

E. Hyman ci al. / Plasma enhanced CVD modeling

ELECTRON DENSITY (M3)

3. Results We present calculations from the modules of the PECVD model that have been constructed and tested. The results obtained are designed to identify important physical mechanisms and to determine the regimes where they are critical. The calculation of the electric and magnetic fields was accomplished using SAIC’s MASK code, a general two-dimensional electromagnetic code designed for the study of microwave devices of arbitrary geometrical configuration. The code introduces the electric fields at the input port and allows them to propagate into the reactor, which can include arbitrarily shaped regions of dielectric or conducting bodies. It employs a finite difference representation of the full set of time-dependent Maxwell equations and solves the initial value problem. Figure 1 shows the results of a simulation for a generic reactor in which the plasma is modeled as a spherical shell with a finite conductivity and a finite dielectric constant. Ultimately the plasma model will be replaced by results of Boltzmann calculations over the plasma region. The MASK calculation retains both the field strength and the phase dependence and determines the energy deposition in the target plasma. The results shown are contours of constant field amplitude for the axial (a) and radial (b) electric components in an azimuthally symmetric configuration. The bottom horizontal line is the axis of symmetry. The input

__________

_________

______

0.~

0~

Input Port

0.10 /

Reactor Axis

_______

Plasma Shell

Substrate

Exit Port

____________________ ________ ____________

(b)

~

0.L~

0.10

0.14

0.19

02*

Fig. 1. MASK calculation of(a) the axial and (b) the radial components of the electric field in a model reactor. Shown in the figure are contours of constant amplitude. The reactor axis is the bottom horizontal line. The substrate shelf is at the right and the outlet for the reacting gases is above it. The plasma shell is centered on the axis to the left of the substrate,

iii ~iii~

1018 H—

-

—~

1016 T

.

..

.

~.

-

1010

-

_•_i

—

-

—

.

1012

—

-

~~08

:/

—

j 1

--—.

10

—

-

—

-

----I

.~_

— —

r~l

—-—-~----—

——

.—~..

1o-

—

——

— -

—

.

—

_._.

10

~

-

-4

—

10

-

10

10

TIME (S) Fig. 2. Time development of the electron density at a point in the reactor of high electric field corresponding to the one-point simulation described in the text.

wave is introduced on the left and propagates into the reactor region through a radially expanding transition region. The substrate is on the right in the figure and immediately above it is the outlet for the reacting gases. The figure indicates regions of high and low field concentration which will provide an important tool for reactor design. This calculation determines the energy deposition in the plasma and, hence, the reactor efficiency. When performed self consistently it also predicts the shape of the plasma region and the subsequent coupling to the hydrodynamic calculation. The Boltzmann module calculation simulates the electron and heavy particle evolution at a location within the reactor where the electric field is sufficiently high to create and sustain a plasma. We maintain a constant deposited microwave power and assume the gas pressure is kept constant at 5.3 x l0~Pa (40 Torr). The simulation runs for several milliseconds; beyond which time advection and diffusion effects, not included in this calculation, would become important. The initial conditions are [H 3, [CH 2] = 1.28 T x =1024 4} =2 l.28x 3, gas temperature 300rn-K. Figure shows 1022 electron rnthe density rising rapidly to approximately l0’~m3 in around 1 ns, a result of the very energetic electron distribution at early times. Thereafter, it increases nearly two more decades over about l0-~s, the slower increase being reflective of the less energetic electron spectrum as the electrons give up energy to the various inelastic processes. The increase on a millisecond time scale is associated with the conversion of the gas from the molecular to the atomic state which leads to an increasing fraction of atomic ions (which recombine much more slowly, than molecular ions) and also causes

E. Hyman et a!. / Plasma enhanced CVD modeling

390

Electron distribution 1

I

I

-

1.5 ns

at i

I

-

1

-

:

1 0

:

~ •

100

.

-

•

•

-

t

0

_____

I

8

4

J

I

1 2

I

1 6

20

Energy (eV) (a)

Electron distribution at 55.2x10 3 ~

6

r

I

$

1 01

1

~

~

detailed balance. The hydrocarbon chemistry is mitiated by the electrons which dissociate H2 (and also the

1•

CH~), ing thecausing gas. Figure the release 4 shows of chemical the gas temperature energy and heatas a of time for the simulation described above. In

1

Fig. show the to evolution 12 hydrocarbon plus 5H-,weand H out 3 ms, atof which time the H species

function

1o

an adjustment in the electron distribution function. The evolution of the electron distribution function as determined by the Boltzmann equation is illustrated for these three time regimes in Fig. 3 in which the ordinate scale is arbitrary. In Fig. 3(a) at 1.5 ns we see the very energetic electron spectrum; the average electron energy is around 8 eV. In Fig. 3(b) at 55.2 u~ the average energy has dropped to 2 eV and in Fig. 3(c) at 4.8 ms when atomic hydrogen predominates, the spectrum has changed again, the average electron energy increasing moderately to about 2.7 eV. The evolution of the hydrocarbon species is simulated with a chemistry code that uses the output of the Boltzmann code. The reactions used in the hydrocarbon modelarelistedinTablelalongwiththeconstantsA, b, and E which determine the rate coefficient k according to k = ATh exp( —E/T). The code calculates the rate for each reverse reaction that is not known, using

2 and H 1

0

-

~ r

1

~\‘,‘

F

-

-‘.‘•-.••\~‘

~

1 0

-

after about 2.5 ms with rise in temperature, thermal dissociation of H2 becomes predominant. The build-up of CH3 due to the dissociation of CH4 occurs very early

\1

1

7

j

\~

-

1 0

densities approximately Although the formation of Hare is initiated by the equal. electron dissociation of H2,

.~

\

1 O~~

11

N

0

4

8

1

2

1 6

(eV)

Energy

20

(b)

Electron distribution at 4.8 ms 3

I

I

I

I

r-

1 0

r 101

I

i~FI

I

I

I

~

-

~.

r

1o~

r

1o~

‘r-

1 ~ 7

_________

i

-

k

-

1 o

I

I

-

.~

~

-

(approximately 30 us) but it reacts with itself to form C2H6 and drops to a local minimum before 1 ms. As the temperature increases, however, the CH3 recombination reaction rate decreases and CH3 increases to a new maximum near 2 ms; thereafter it decreases once more as the CH4 becomes exhausted. Acetylene (C2H,) results from the chain of reactions initiated by the formation of C,H6, thence to C2H5 and in turn to C2H4, C2H3 and finally to C2H2 which persists to the end of the simulation. In general, the hydrocarbon species do not have time to reach the equilibrium values that the gas temperature would dictate. Thus, the time between their formation in the plasma and their reaching the substrate determines the densities of the critical radical species reaching the substrate.

-

r -

-

0 (c)

-

I

4

I

__________

8

Ener

1 2

“ .‘

(eV’ ‘

/

I

—

J 1 6

20

Fig. 3. Electron energy distribution for the simulation as in Fig. 2: (a) at a time before inelastic processes reduce the average electron energy: (b) at an intermediate time when the average electron energy is about 2 eV; (c) at a time when dissociation of the H, is nearly complete. The ordinate scale is arbitrary.

E. Hyman et al. / Plasma enhanced CVD modeling

391

‘

TABLE I. Hydrocarbon reactions and rate coefficients 60’ ‘1(m)’°

Reaction

A(l0

H+H+H,—.H,+H, H,+H CH 2-.H+H+H2 4+H—~CH+H2 CH,+H2—.CH4+H CH4 + ~ —~CH,+ H + CH3+H+c~4~CH4+c~4 CH+CH,—.C2H,+H C2H5+H-.CH+CH3 CH+CH,-.C,H4+H, CH~+c~~lCH2+H+c~ CH, + H CH + H2 CH2+CH,—*C,H4+H C2H6 + H —*C,H, + H2 C2H6+CH—IC,H,+CH4

I

s1)~

2.7 x 10” 9 l.5xl0’ 3.6 x lO_20 1.1 x 10.21 (016)33 x l0~ (8.6(2.2 x 1021 1.3 x l0~ 5.Ox 10_il 1.7 x I0~ (18.6(17 x l0~ 6.6 x 10-I I 6.6 x 10-’’ 9.0 x I0_22 9.1 x l0_25 I~:~l.7 x 10’ 2.3 x l0_12 l.7x l0 2.5x I0~° (is6)~.3 X l0~ (~‘~4.3 x l0~ 7.0 x 10” 3.3 x 10 1f86~5.0< io-~ (~~l.lI x l0—° lOx 10’° 2.5x lO~ 3.0 x 1012 5.8x10’’

C2H,+C,H,—~C2H4+C,H6 ~ C2H4+H—*C2H,+H, C2H4+c~~*C2H3+H+~~4 C,H4+CH,-.C,H,+CH4 C,H+H—~C,H2+H, C,H, + ~ —~C2H2+ H + ~ C2H2 + H + ~!~4—lC2H,+ C,H,+H-.C,H+H, C,H+H,-.C,H2+H C,H,+CH,-.C,H+H C2H,+C,H—IC4H,+H C2H2 +v~4_C2H+ H + ~

(18.6)6.6 x l08

b

E(K)

Range(K)

—0.6 0

0 4.84x104

100—5000 2500—8000

3.0 3.0 0 —3.0 0 0 0 0 0 0 3.5 4.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

4.40 x 10’ 3.90 x 10’ 4.45 x l0~ 0 1.34 x l0~ 0 1.61 x l0~ 4.56 x l0~ 0 0 2.62 x 10’ 4.17 x 10 3.43 x l0~ 0 1.56x lQ~ 5.14x 10’ 3.99 ~< l0’~ 4.86 x l0~ 5.59 x 10’ 0 1.61 x l0~ 3.5 x 102 l.19x l0~ l.56x 10’ 0 0 5.38 x l0~

300—2500 300—2500 1500—3000 300—2500 1500—3000 300—1500 1500—2500 1500—3000 300—2500 300—2500 300—2000 300—2000 800—2500 300—1200 700—1500 700—2000 1500—2500 1500—2500 300—1000 300—2500 500—2500 300—500 300—3000 300—3000 >298 300—2500 1500—3500

n denotes the number of reactants.

Temperature (K) 4.0 10

itiiiiiui

iuIIIIIIIrIiIIIy)..J._

3 3.0 10

/

3 2.0 10

-

7

/ / /

-

-

_._._._ -

3 1.0 10

-

_._

._._._~

-

-

0.0

-

IIIIIIIIIIIIIIIIIII

T

hI~11I~

ime ~ms,

Fig. 4. Evolution of the gas temperature for the simulation as in Fig. 2.

Finally, we present preliminary results of a fluid simulation of a generic PECVD reactor, using SAIC’s FUGG code which is capable of performing fluid calculations over arbitrarily complex geometries. The code

employs an unstructured grid allowing extremely fine resolution in critical areas while employing coarser griddinginregionswherequantitiesvaryslowly.The code was designed for the study of flow problems dominated by convection and is presently being modified to incorporate thermal conduction and viscosity effects. In Fig. 6 we show two examples of the code’s triangular gridding capability. In both cases a crosssection of the azimuthally symmetric model reactor is shown, where the left vertical boundary represents the reactor axis. Three gas inlet ports are modeled allowing the gas to enter at the top (in Fig. 6(a) the plenum region above the inlet ports is also modeled). The gas exits through the horizontal boundary at the bottom right.ThesubstratewaferisrepresentedinFig.6(a)by the left half of the lower horizontal boundary and in Fig. 6(b) by the shelf at the lower left. In Fig. 6(a) the variable gridding capability is clearly illustrated and, in particular, the fine gridding needed in the inlet ports is shown. Figure 7 shows results of a fluid calculation for the reactor of Fig. 6(b) in which hydrogen gas enters at 50 m s’ at a pressure of 5.3 x l0~Pa (40 Torr). A heating source of 1.5 kW over a spherical volume of radius 0.035 m, centered on the reactor axis and mid-

E. Hyman et al. / Plasma enhanced CVD modeling

392

HYDROCARBON DENSITIES (M 3) 1024 CH4 1022.

d

~

~-

1 o~°

lnpt*pc*ts

,•

1018~

CH

-

/

3

CH2

~/ ..~

~

/CH

/~\

1016_

0

1

2

3

reactor axis

~-

Time (ms) (a)

~

HYDROCARBON DENSITIES (M

3)

k

‘

1024

-

.

C2 H2 1022.

1 020

/

H

- ~

)~ -4~

1018.•

C~H6 -

1016..

0

-

/ C H .;

~

.~.

..

2 5:

1

••~

1

C H j—_

2 Time (ms)

.

I,,

exit pan

j~Z. C2H3 :c~w-~ C3H3

~

.

~

~

0.10

4 2

3

(b) Fig. 5. Evolution of the hydrocarbon densities for the simulation as in Fig. 2.

Fig. 7. Velocity field for reactor of Fig. 6(b) in which 1-17 enters at a velocity of 50 ms’~ at a pressure of 5.3 x 10’ Pa (40 Torr). A heating source of 1.5 kW to simulate the effect of the plasma is centered on the axis between the inlet port and the wafer in a spherical volume of radius 0.035 m. The ordinate and abscissa dimensions are in meters.

r

~ L Fig. 6. Examples of the FUGG code’s unstructured gridding capability for two model reactors ((a) and (b)).

(b)

E. Hyman ci a!. / Plasma enhanced CVD modeling

way between the inlet port and the wafer, is included to simulate approximately the effect of the plasma source. Shown are velocity vectors for the flow 2.6 ms after the plasma is turned on. Also calculated but not shown are the pressure, density, and temperature fields. While conclusions should be tempered because of the current lack of inclusion of thermal effects in the code and because the results represent a transient pre-steady-state stage, the effects of buoyancy are apparent. The complex vortex flows seen suggest this reactor configuration would be very poor for efficient diamond deposition.

4. Discussion and conclusions We have presented results of a model under development that will permit the simulation of PECVD reactors of arbitrary geometry. The model will be an important tool providing better understanding of the microscopic processes occurring within the reactor, permit parameter studies to identify those parameters which critically affect both the rate of deposition and the uniformity of the deposition over the wafer surface, and ultimately enable the design of improved reactors. We have identified several critical elements in the modeling effort that need to be treated carefully if simulation results are to be meaningful. First, the electromagnetic fields which initiate the plasma formation need to be determined in the realistic reactor geometry, including

393

effects of all metallic, dielectric, and insulator elements actually present, to ensure that the fields are high in the desired plasma formation region but not elsewhere. Second, the coupling of the fields to the plasma electrons, to determine accurately both the time development of the electron density and their energy distribution, is most important for determining the evolution of the rate of hydrogen dissociation and the rise in the gas temperature. This, in turn, critically determines the non-equilibrium hydrocarbon chemistry development, a third area that needs to be carefully modeled. Finally, the flow of hydrocarbon radicals to the wafer is very sensitive to the reactor’s geometrical configuration, its thermal properties, and the location of the plasma relative to the wafer. In conclusion, our results emphasize the highly non-equilibrium and coupled nature of PECVD reactor processes and the strong influence of reactor geometry. A numerical simulation that is useful must address all such issues.

Acknowledgment This research was supported by the U.S. Army Missile Command and sponsored by the Defense Advanced Research Projects Agency (DARPA) under contract DAAHO1-90-C-0279 and by the Independent Research and Development Program of Science Applications International Corporation (SAIC).