Thin-film deposition modeling of hydrogenated amorphous silicon in the afterglow of argon plasma

Accepted Manuscript Thin-Film Deposition Modeling of Hydrogenated Amorphous Silicon in the Afterglow of Argon Plasma S.W. Chau, Q.T. Pham, M.N. Nguyen PII: DOI: Reference:

S0045-7930(13)00366-6 http://dx.doi.org/10.1016/j.compfluid.2013.09.021 CAF 2303

To appear in:

Computers & Fluids

Received Date: Revised Date: Accepted Date:

15 January 2013 16 September 2013 17 September 2013

Please cite this article as: Chau, S.W., Pham, Q.T., Nguyen, M.N., Thin-Film Deposition Modeling of Hydrogenated Amorphous Silicon in the Afterglow of Argon Plasma, Computers & Fluids (2013), doi: http://dx.doi.org/10.1016/ j.compfluid.2013.09.021

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Thin-Film Deposition Modeling of Hydrogenated Amorphous Silicon in the Afterglow of Argon Plasma S. W. CHAU*, Q. T. Pham, M. N. Nguyen Dept. of Mech. Eng., National Taiwan University of Sci. and Tech., Taipei, Taiwan, ROC *Corresponding

author, [email protected]

ABSTRACT: This paper proposes an integrated numerical approach to model the thin-film deposition of amorphous silicon in the afterglow of argon plasma. The investigated plasma enhanced chemical vapor deposition (PECVD) process is governed by gas flow equations, chemical kinetics equations and surface reaction equations. Nineteen species, i.e. Ar, Ar*, Ar+, Si2H6, Si2H4, SiH4, SiH3, SiH2, SiH, H2, H, Si, Si+, SiH3+, SiH2+, SiH+, H2+, H+ and electron (e) are considered, where fifty-two chemical reaction equations are adopted to describe the reactions in space among silane, hydrogen and argon plasma. Continuity, momentum and energy equations are solved to describe the velocity, pressure and temperature field of gas flow in the process chamber, while the density of neutral species is calculated from the corresponding transport equation considering its convection and diffusion effect, as well as the production and destruction due to chemical reactions. The transport equations of charged species are modeled by a drift-diffusion approximation. The ambipolar diffusion is assumed for charged species in the afterglow of argon plasma due to the vanishing of external electrical field. The equation of energy balance for electron is applied to the estimation of the electron temperature in the process chamber. In this study, ten reaction mechanisms of thin film growth are taken into account. The species considered in the deposition process on substrates are silane, the free radicals of silane, hydrogen, hydrogen atom, silicon atom, the species with dangling bond on the substrate surface covered by active and the passive site. All governing equations are discretized by using a finite volume method and solved by an iterative approach. The variables, such as position, velocity and moment, are expressed by Cartesian coordinates, where the non-staggered

variable arrangement is adopted to define field variables. The deferred correction approach is employed to approximate for both diffusion and source terms, where SIMPLE algorithm is used to decouple pressure and velocity. A parallelized numerical scheme is implemented via MPI library. Numerical simulations with different operation parameters are performed to show their influence on the deposition rate. The numerical results indicate the deposition rate exhibits a double-peak characteristic on the substrate and give an average value varying between 5 Å/s and 10 Å/s for the studied PECVD process in the afterglow of argon plasma. Keywords: PECVD, Afterglow, Deposition Rate, Argon Plasma.

Nomenclature A

Area of the substrate

Al

Fitted values of parameters to determine the rate constant

Bl

Fitted values of parameters to determine the rate constant

Cl

Fitted values of parameters to determine the rate constant

ce

Mobility constant of electron

ci

Mobility constant of ions

cp

Specific heat

Da

Ambipolar diffusion coefficient

De

Diffusion coefficient of electron

Di

Diffusion coefficient of ions

Dij

Diffusion coefficient of binary mixture

Dφ

Diffusion coefficient

dx

Dimensionless cell size

da

Distance between two argon orifices

do

Pitch between two plasma orifices

dg

Distance between two gas rings

E

Electron field

Ec

Relative discretization error

Einel,l

Rate constant of inelastic collision in the reaction l

ΔEinel,l

Energy loss of electron due to inelastic collisions in the reaction l

f

Gravitational force

Fi

Flux of species i on wall

ha

Height of chamber

hg

Distance between gas rings and substrate

kb

Boltzmann constant

kl

Rate constant for the reaction l

Kel,j

Elastic collision rate constant

Mj

Molecular weight of species j

M*

Mean molecular weight of three main gases

me

Mass fraction of electron

mj

Mass fraction of main gas j

Nc

Number of cells

ne

Density of electron

ni

Density of ions

nt

Number of particles per unit area

P

Pressure

Pw

Power of plasma source

R

Universal gas constant

Rsource

Production rate

Rsink

Consumption rate

Sφ

Source term

*

T

Dimensionless temperature

Te

Electron temperature

Tg

Gas temperature

Tw

Wall temperature

U

Velocity vector

Vi

Thermal velocity of species i

wx

Width of chamber

wz

Depth of chamber

yi

Molar fraction

αs

Deposition rate

∞

α

Grid-independent

ε

Mean electron energy Field variable

Γe

Electron flux

Γi

Ion flux

γl

Sticking rate constant for the deposition reaction l

η

Viscosity

ηj

Viscosity of species j

κ

Thermal conductivity

κj

Thermal conductivity of species j

μe

Mobility of electron

μi

Mobility of ions

θa

Fraction of surface covered by active sites

θp

Fraction of surface covered by passive sites

ρ

Density

σi

Collision diameter of species i

σij

Binary collision diameter between two species i and j

Ω

Collision integral

1. Introduction

For decades plasma supported processes, such as plasma-enhanced chemical vapor deposition (PECVD), have been widely applied in the manufacturing industries for the purpose of surface cleaning, coating and etching [1]. PECVD, which has been used for various applications, for instance active-matrix liquid crystal displays, linear image sensors for facsimile machines, microelectronics, computer chip, notably transistors and solar cells, plays an important role especially in the industrial production of silicon-based thin films [2]. Both gas-phase and surface reactions are important in this plasma assisted process starting from the dissociation of the feed gas activated by plasma, where electrical energy is first employed to generate a glow discharge (plasma) in the environment, and the energy is then transferred into a gas mixture composed of different types of species, i.e. reactive radicals, electron, ions, neutral atoms and molecules, and other highly excited species. These atomic and molecular fragments later interact with designated substrate, where either etching or coating process occurring on the substrate depends on the nature of chemical interactions. The reactive and energetic species in the gas phase are mainly formed by the collision reaction in the gas phase and hence the gas temperature becomes an essential issue in the PECVD process. The substrate is usually maintained at a temperature range between 25°C and 400°C because low working temperature leads to several desirable properties of thin film, such as superior adhesion, low pinhole density, good step coverage, and high uniformity [3]. In addition, PECVD process often requires reactions between gas-phase precursor components. Therefore the operational pressure typically adopted in the system

is ranging from 0.1 to 10 Torr, which ensures adequate low gas pressure to sustain the discharge process [3]. There are two major types of PECVD devices. The first type is based on a simple direct-current (DC) discharge, which can be readily produced at a few Torr between two conductive electrodes and is suitable for the deposition of conductive materials. Nevertheless, insulating films could rapidly extinguish DC discharge as they are deposited. Furthermore, excitation of capacitive discharge by applying an alternating-current (AC) or radio-frequency (RF) signal between an electrode and the conductive walls of a reactor chamber, or between two cylindrical conductive electrodes facing one another is the second type, which is known as a parallel plate reactor [4]. In RF-PECVD process, two distinct frequencies are employed. The former one is in the low-frequency range, usually around 100 kHz, requires several hundred volts to sustain the discharge. The latter one is in the high-frequency range, at the frequency of 13.56 MHz. The advantages of PECVD process in comparison with the other types of chemical vapor deposition processes are low working temperature, convenient control over film stress, less demand for physical/chemical chamber cleaning, and offering a wide range of deposition material, such as organic (plasma polymerization) and inorganic film (amorphous hydrogenated silicon), silicon dioxide (SiO2), silicon nitride (SiN), silicon carbide (SiC), poly silicon (poly-Si) and diamond like carbon (DLC) [5]. However, limited capacity, wafer loading problem and easily suffering from contamination are typical drawbacks of PECVD. Traditional approach, where the substrate is arranged in the glow region, can result in immediate interaction of substrate with neutral and ionic species composing the plasma. Neutral

species predominantly contribute to space and surface chemical reactions and charged particles are the main reason for positive ion bombardment. The combination of these two reaction types give unique deposition kinetics on the substrate, which is characterized by the region of high monomer fragmentation and a coated film with crosslinked network, as well as high surface energy indicating reduced monomer retention [6]. To relocate the substrate several centimeters downstream from the glow regime, i.e. in the afterglow region of plasma, can easily overcome this problem due to the quick diminishing of, or even absence of, the density and the energy of electrons and ions, which clearly deliver a film of low surface energy. Experimental study on the silicon nitride deposition in the afterglow of Ar-N2 plasma excited by RF power supply was conducted with a gas mixture of SiH4 and Ar as the feed gas [7]. The best deposition rate of 14Å/s for silicon nitride film with a refractive index of 1.95 was found at the substrate temperature of 400°C. The deposition rate between Si and SiN thin film was experimentally compared in the afterglow of N2 and H2 microwave plasma, where the feed gas SiH4 was separately injected into the reacting tube 50 mm in front of the circular substrate [8]. The measurement suggested that the silicon film formation was mainly governed by the SiH4 dissociation and reactions among silane species, which were directly related to the density of H2 plasma. In addition to the conventional model of SiH3 radical contributing to the Si film formation, the experimental result also implied that a model of SiH2 radical should be also essential in the deposition process. Fluorocarbon coatings from hexafluoropropylene oxide with the assistance of RF-PECVD obtained from the glow and afterglow arrangement was investigated [6],

where chemical and morphological features of thin film were extensively examined. The introduction of afterglow approach instead of a glow one could allow the hydrophobic surface to give superhydrophobic behavior. Recently the afterglow approach adopted in a DC glow discharge was proposed to enhance the film deposition rate of carbonitride using a gas mixture of CH4 and N2 with successful result [9]. The available works devoted to the film formation in the afterglow of plasma were mostly based on experimental measurements, and the theoretical study on the afterglow deposition via fully three-dimensional modeling is not available. Gas flow, chemical reactions and plasma dynamics are three major parts affecting the coating performance of PECVD process, which should be all included in the mathematical model and solved in a coupled way in order to perform a thorough analysis of plasma discharge process. In this work, kinetic simulation incorporated with fluid dynamics is adopted to investigate low-pressure and non-equilibrium situations of thin film deposition in the afterglow of RF plasma. This study proposes an integrated numerical model to account for the chemical reaction between plasma and working gases, along with the deposition mechanism on substrate surface, where the influences of processing conditions on the flow characteristics of afterglow-type PECVD process, as well as the growth rate of thin film, are predicted. The modeled plasma reactor is a prototype of afterglow-type PECVD, Fig. 1, where plasma source of hollow-cathode discharge is employed to provide high-density plasma flow targeting uniform thin film growth for large area.

2. Governing Equations

As depicted in Fig. 1 argon plasma is generated in the RF plasma source located above the reaction chamber, and then injected into the reaction chamber through two rows of small orifices connecting the process chamber and discharge space between electrodes. Unlike the small discharge space inside the plasma source subjected to strong and time-varying electric field, the process chamber does not experience any external electric field, where the system pressure, the substrate temperature, the input power of plasma source, and the flow rate of feed gas and argon plasma are all kept constant during the deposition process. A relatively weak electric field induced by the ambipolar diffusion is considered inside the process chamber, and thus a steady state approach can be presumed. Continuity, momentum and energy equations in the steady form are employed to describe the velocity, pressure and temperature field of gas flow in the process chamber. The steady continuity equation of gas flow is given as ∇ ⋅ ( ρU ) = 0 ,

(1)

where ρ is the density and U the velocity vector. The steady momentum equation of gas flow is expressed as U ⋅ (∇ρU ) = −∇P + ∇ ⋅ (η∇U ) + ρf ,

(2)

where P is the pressure, ηthe viscosity and f the gravitational force. The steady energy conservation of gas flow can be written as

U ⋅ ∇( ρc pT ) = ∇ ⋅ (κ∇T ) ,

(3)

where cp is the specific heat, κ the thermal conductivity and T the gas temperature. Because the gas mixture is mainly composed of three main species in the process chamber, i.e. SiH4, H2 and Ar, the contribution of other species to the fluid properties is neglected. The fluid properties cp, ρ, ηand κare calculated by the following equations: c p = ∑ m j c p, j ,

(4)

PM * , RT

(5)

j

ρ=

1

η

=∑ j

mj

ηj

(6)

,

κ = ∑ m jκ j ,

(7)

j

where the subscript j indicates three main gases considered in this study, mj the mass fraction of species j, cp,j the specific heat of species j, R the universal gas constant, ηj the viscosity of species j, κj the thermal conductivity of species j and M* the mean molecular weight of three main gases and defined by the following equation: M* =

1 mj ,

∑M j

(8)

j

where Mj is the molecular weight of species j. Table 1 summarizes the fluid properties of three main gases [10]. A laminar flow is assumed in the process chamber because of small Reynolds number (in the order of 103), while the gas density is determined by the ideal gas law due to constant and small chamber pressure P (in the order of 1 Torr) maintained in the process chamber. In order to describe the behavior of the chemical interactions among reacting gases, transport equations for the considered species, which participate in the deposition process in the process chamber, are solved. The following nineteen species, i.e. Ar, Ar*, Ar+, Si2H6, Si2H4, SiH4, SiH3, SiH2, SiH, H2, H, Si, Si+, SiH3+, SiH2+, SiH+, H2+, H+ and electron (e) are considered. Fifty-two equations of chemical reactions are employed to describe the interactions among SiH4, H2 and Ar plasma in the chemical reactions [11, 12], which are summarized in Table 2. The considered reactions include the electron-impact reactions with argon (3 equations), the metastable reaction (1 equation), the electron-impact reactions with silane (10 equations), the electron-impact reactions with hydrogen (7 equations), other electron-impact processes (6 equations), the neutral radical reactions (9 equations), the excitation transfer and radiactive decay (6 equations), ion-neutral reactions (8 equations) and collisions with disilane reactions (2 equations). The rate constant of reactions adopted in Table 2 is given as the following form:

kl = AlTeBl exp(

− Cl ), Te

(9)

where kl is the rate constant for the reaction l, Te the electron temperature, and Al, Bl and Cl the fitted values of parameters to determine the rate constant. The transport equations of electron and ions are described by the drift-diffusion approximation, where the steady transport equations of electron and ions are given as follows:

∇ ⋅ (ne U) + ∇ ⋅ Γe = Re,source − Re,sink ,

(10)

∇ ⋅ (ni U) + ∇ ⋅ Γi = Ri,source − Ri,sink ,

(11)

where ne is the density of electron and ni the density of ions, Rsource and Rsink the production and consumption rate of species, respectively. The electron flux Γe and ion flux Γi are written as follows: Γ e = − De ∇ n e − μ e n e E ,

(12)

Γ i = − Di ∇ ni − μ i ni E ,

(13)

where De is the diffusion coefficient of electron, Di the diffusion coefficient of ions, μe the mobility of electron, μi the mobility of ions and E the electron field. For neutral species other than charged particles, the steady transport equation is given in the following form:

∇ ⋅ (ni U) + ∇ ⋅ Γi = Ri ,source − Ri ,sink ,

(14)

where ni is the density of neutral species i. The flux Γi of neutral species i is defined as follows: Γ i = Di ∇ ni ,

(15)

where Di is the diffusion coefficient of neutral species i. The diffusion coefficient of neutral species Di is calculated by the following equation:

Di =

∑y

j

j

yj

∑D j

,

(16)

ij

where i represents the neutral species other than three main background gases (Ar, SiH4 and H2) and j stands for three main background gases, yj the molar fraction of three main background gases and Dij the diffusion coefficient of binary mixture between the pair of species i and background gas j. The

diffusion coefficient between two neutral species is determined by the binary mixing theory [13]:

Dij =

⎛ 1 1.8583 × 10−7 1 ⎞⎟ 1 T 3⎜ + , ⎜ M M ⎟ Pσ 2 Ω 760 j ⎠ ij ⎝ i

(17)

where σij is the binary collision diameter between two species and Ω the collision integral. The binary collision diameter is defined as the arithmetic average of collision diameters:

1 2

σ ij = (σ i + σ j ) ,

(18)

where σi is the collision diameter of species i. The dimensionless collision integral Ω is obtained from the following fitted equation: Ω=

1.06036 0.19300 1.03587 1.76474 + * + * + * , T *0.15610 e 0.47635T e1.52996T e 3.89411T

(19)

where T* is the dimensionless temperature defined as

T* =

T

,

(20)

⎛ε ⎞ ⎛ε ⎞⎛ε ⎞ ⎜ ⎟ = ⎜ ⎟⎜ ⎟ , ⎝ K ⎠ ij ⎝ K ⎠i ⎝ K ⎠ j

(21)

⎛ε ⎞ ⎜ ⎟ ⎝K⎠

⎛ε ⎞ where ⎜ ⎟ is the Lennard-Jones parameter. The collision diameter and Lennard-Jones parameters ⎝K⎠ for neutral species considered in this study are given in Table 3 [14]. Because the density of electron and ions are not low enough to be neglected, a considerable space charge is formed as the result of charge separation and the generated polarization field impedes further violation of charge neutrality, where the ambipolar diffusion is adopted to describe the diffusion phenomenon between electron and ions [15]. The ambipolar diffusion coefficient Da required in the transport equations of charged species is expressed as follows:

Da =

μi De + μe Di . μi + μ e

(22)

The diffusion coefficient of ions Di can be calculated by the following equations [15]: Di = μ iT ,

(23)

μi =

ci =

cz =

ci , P

(24)

∑y

z

z

y ∑z c z z

∑y

yj

j

(25)

,

(26)

j

j

∑c

,

j z

(M z + M j ) c zj , = j c j+ 2M z

(27)

where the subscript z indicates the ion species Ar+, SiH3+, SiH2+, SiH+, H2+, H+ and j the three main gases, cz the mobility constant of ion z, ci mobility constant of ions and c zj mobility constant of ion z in the gas j. The mobility constant of ion in its corresponding state gas are summarized in Table 4 [16]. The diffusion coefficient of electron De is calculated with a similar approach: De = μ eTe ,

μe =

ce =

(28)

ce , P

(29)

∑y

j

j

yj

∑c j

,

(30)

j e

where ce is the mobility constant of electron. The mobility constant of electron in different gases are given in the Table 5 [15, 17]. By adopting the assumption of ambipolar diffusion, the electron and ion flux is expressed as the following equation: Γ e = − Da ∇ n e ,

(31)

Γ i = − D a ∇ ni .

(32)

The equation of energy balance for electron is applied to estimate the electron temperature in the chamber [18, 19]:

5 m ⎛5 ⎞ ∇ ⋅ ⎜ εΓ e − ne De∇ε ⎟ = −cΓ e ⋅ E − ∑ ΔEinel,l K inel,l − 3∑ e K el, j kb (Te − T ) , 3 ⎝3 ⎠ l j mj

(33)

3 where ε is the mean electron energy (= k bTe ), ∆Einel,l the energy loss of electron due to inelastic 2 collisions in the reaction l, Kinel,l the rate constant of inelastic collision in the reaction l, Kel,j the rate constant of elastic collision with main gas j, me the mass fraction of electron, mj the mass fraction of main gas j, kb the Boltzmann constant and c the electric charge of an electron. The subscript l indicates the reactions given in Table 3, which result in energy loss of electron, due to inelastic collisions, i.e. (Reactions 1-3, 7-14, 18-19, 21, 22-27 in Table 2). The elastic collision rate constant Kel,j is set equal to 10-13 (m3/s) for three main gases [4, 20]. The electron field E induced in the chamber as a result of ambipolar diffusion is presented as the following equation [15]:

E=

Di − De ∇ne . μi + μe ne

(34)

The energy loss due to inelastic energy collision is expressed as the following equation:

∑ ΔE

inel,l

l

K inel,l = ∑ Cl kl ni ne ,

(35)

where Cl and kl are the parameter and rate constant for the reaction l in Table 3. The species considered in the deposition process on substrates are silane (SiH4), neutral radicals of silane (SiH3, SiH2, SiH), hydrogen (H2), hydrogen atom (H), silicon atom(Si), the species with dangling bond on the substrate surface (SiH-S, Si-S) as active sites and the passive site (Si-B), where

θa and θp denote the fraction of surface covered by active and passive sites, respectively, and θa +θp =1 [3]. Table 6 lists the reaction mechanism of thin film deposition considered in this study. The rate charge of active site is then given as follows [21-23]:

dθ a 1 = − γ1 FSiH3θ a + γ2 FSiH3 (1 − θ a ) − γ 5 FSiHθ a + γ 6 FSiH (1 − θ a ) − γ 9 FHθ a + γ 10 FH (1 − θ a ) , dt n

[

]

(36)

where γl is the sticking rate constant for the deposition reaction l on substrates given in Table 6, nt the

1 number of particles per unit area, Fi the flux of species i to the substrate (= niVi ) and Vi the thermal 4 velocity of species i (=

8 RT ). For a steady-state approach the left-hand-side term vanishes and the πM i

active site θa is consequently obtained as:

θa =

γ 2 nSiH VSiH + γ 6 nSiHVSiH + γ 10 nHVH . (γ 1 + γ 2 ) nSiH VSiH + (γ 5 + γ 6 ) nSiHVSiH + (γ 9 + γ 10 ) nHVH 3

3

(37)

3

3

The deposition rate αs on the glass substrate is then determined by the following expression:

α s = (γ 1FSiH θa + γ 2 FSiH θ p + γ 3 FSiH θa + γ 4 FSiH θ p + γ 5 FSiHθa + γ 6 FSiHθ p + γ 7 FSiθa + γ 8 FSiθ p ) 3

3

2

2

.

(38)

3. Numerical Scheme The integral form of governing equations obtained through a finite volume discretization of governing equations is given as follows:

∫ (φU ⋅ n)dS = ∫ (Dφ∇φ ⋅ n)dS + ∫ Sφ dV , s

s

V

(39)

where the control volume V is defined by the control surface S with the outer normal n, φ the field variable, Dφ the diffusion coefficient and Sφ the source term. The kinetic model is first expressed as a system of partial differential equations, which are further discretized by using the finite volume approach and solved by an iterative approach. The variables, such as position, velocity and moment, are expressed by Cartesian coordinates, where the non-staggered variable arrangement is adopted to define the variables. All variables are stored in the cell center. In order to achieve highly accuracy in the discretized equation, both upwind difference (UD) and central difference (CD) methods are applied to calculate the convention term with the deferred correction approach and the latter one is employed to approximate for both diffusion and source terms [24]. The solution process of the coupling equation system starts from solving the momentum equations, Eq.(2), for temporary velocity components in the computational domain with a given pressure field. The gas temperature is then obtained by solving the energy equation of gas mixture, Eq.(3), and the pressure field is then updated with the help of the continuity equation, Eq.(1), where the temporary velocity field previously calculated from the momentum equations is employed. The density of electron is estimated from the transport equation of electron, Eq.(10), while the density of ions is predicted via the corresponding transport equation, Eq.(11). For neutral species, Eq. (14) is then adopted to

compute their density. The electron temperature is attained by the solution of the energy equation of electron, Eq.(33). The deposition rate is computed from Eq.(38), after the active site of the substrate is determined based on the available density of species on the substrate. The fluid and species properties in space are then updated using the new gas and electron temperature field. The abovementioned procedure defines a complete outer iteration in our calculation. The proposed segregated and iterative scheme is proven to provide sufficient coupling among the differential equation system when decoupling them from each other. A modified SIP (strongly implicit procedure) solver described in [25] is employed to solve each governing equation in turn for attaining specific field variables, such as the density of species, which defines an inner iteration in our calculation. SIMPLE algorithm is used to decouple pressure and velocity as updating the pressure field in an outer iteration, while the proposed numerical scheme is parallelized and implemented via MPI library [26]. 4. Numerical Results Figure 2 gives the schematic illustration of the computational domain and the employed coordinate system. The principle dimensions of investigated process chamber are summarized in Table 7, where ha expresses the height of chamber, hg the distance between gas rings and substrate, da the distance between two argon orifices, do the pitch between two plasma orifices, dg the distance between two gas rings, wx the width of chamber and wz the depth of chamber. Due to the dense arrangement of plasma orifices along the depth direction only part of the process chamber is modeled, where the computational domain with approximate 180000 hexahedral cells is typically adopted to discretize the region only containing two plasma orifices along the depth (z) direction. The boundary conditions of field variables are given in Table 8, where Tg denotes the gas temperature and Tw the wall temperature. At the plasma orifice the density of Ar, Ar+, Ar* and electron are specified, while the flow rate of silane and hydrogen are given at the gas ring, respectively. In this study, the Chemkin software is used to estimate the plasma composition at the plasma orifice, which has been

validated by the corresponding experimental measurement [27]. The species flux on solid wall is determined from its thermal velocity. The processing condition for the standard case is Tg=298 K, Tw=523 K and P=1 Torr, where the electron temperature at the plasma orifice is attained with the value of 1.8 eV for the power of plasma source Pw=750 W. The total flow rate of silane and hydrogen is identical to that of argon (300 sccm) and the ratio of flow rate between silane and hydrogen is 9 to 1. For a typical calculation the convergence is assumed when the residual of all variables is reduced by six orders of its magnitude. Figure 3 depicts the dependence of average deposition rate of substrate (α) on the number of cells (NC) in a half model and the corresponding grid-independent solution α ∞ (=7.89 Å/s) obtained from the Richardson’s extrapolation [28], where solid line denotes the numerical predictions of different grids and the dotted line the grid-independent solution. The average deposition rate is defined as follows:

α=

1 α s dA , A∫

(40)

where A denotes the area of the substrate. Table 9 compares the average deposition rate of a typical case predicted in this study with the estimated deposition rate via numerical approach of a recent reference [14], where the reference case was modeled in the glow arrangement. Two case were both at the same system pressure (1 Torr) and substrate temperate (532 K). Despite of different input power and flow rate, two studies gave the same tendency of change in the deposition rate as the gas pressure and input power varied in the opposite direction. The predicted deposition rates, 2.5 Å/s and 7.89 Å/s, were of the same order of magnitude falling between 1 Å/s and 10 Å/s. The deposition rate of present study was about three times greater than that of the previously published result because of a much higher flow rate of feed gas and an afterglow arrangement [9]. This comparison with previous published numerical result provides a justification of the proposed mathematical model adopted in this study, which is proven appropriate for the investigated problem. Figure 4 shows the relative discretization error Ec of the proposed numerical method denoted by the solid line, as well as that of an ideal second-order method represented by the dotted line. The relative discretization error

is defined as follows: Ec =

α ∞ −α α∞

,

(41)

where dx is the dimensionless cell size. This comparison indicates the approximation accuracy of the proposed numerical method is close to a second-order one. Figures 5 to 8 illustrate the density distributions of considered species on the center plane (z=0) of the process chamber for the standard case, where a zoomed view of computational domain is shown to highlight the species variation near the substrate. Figure 5 shows the density distribution of feed gases SiH4 and H2 inside the chamber, as well as argon, the plasma gas. The densities of background gases obtained in the process chamber are of similar order of magnitude (1015/cm3), while they are at least five-order-of-magnitude larger than other species, which justifies our assumption in determining the effective gas properties by taking only these three gas components into account. An obvious low-density region of silane and hydrogen is found between gas rings and below the argon plasma orifices, which is clearly due to the consumption of feed gases in the chemical reactions with energetic argon plasma. This region is regarded as the major location where the chemical reactions mainly occur. Figure 6 suggests that the density of ion and metastable of argon quickly decays from the plasma orifice, while the electron temperature slowly decreases in the main reaction zone due to the energy gain from the self-induced electrical field in the ambipolar diffusion process. The free radicals as a consequence of dissociation reactions are very important to the enhancement of deposition process on the substrate. Figure 7 depicts the density distribution of neutral radicals, such as Si2H6, Si2H4, SiH3, SiH2, SiH, Si and H. The most dominant components of neutral radicals are Si and SiH3, which are the major contributors in the deposition of thin film. The density of Si and SiH3 is in the order of 1012/cm3. Other neutral radicals play trivial role in the deposition process due to their small densities. The densities of ions dissociated from silane and hydrogen, i.e. SiH3+, SiH2+, SiH+, Si+, H2+ and H+, are represented in Fig. 8, which give the maximum value in space below the plasma orifice and decline toward the substrate and gas outlet.

The major components of ions in the process chamber are H2+, SiH3+, SiH2+ and SiH+, while Si+ and H+ deliver much smaller magnitude of density than other ions. Figures 9 to 13 show the comparisons of predictions between two- and three-dimensional calculations for parametric study with the standard case, such as the system pressure, the input power of plasma source, the chamber height, the temperature of substrate and the distance between plasma orifices. The numerical predictions suggest that two-dimensional results give a minor over-prediction (about 0.5 Å/s) when compared to its three-dimensional counterpart. This is mainly because of a very dense arrangement of plasma orifices is adopted in this study. Figure 9 (a) reveals the deposition rate along the centerline (y=0 and z=0) for three different chamber pressures, where a plateau region with double peaks characterizes the deposition behavior on the substrate. The region of substrate located between two plasma orifices delivers high deposition rate of 12 Å/s for P=1 Torr, while it sharply fall to 5 Å/s for the region outside plasma orifices. The average deposition rate gradually decreases from 8.4 Å/s to 6.2 Å/s as the system pressure grows from 0.5 Torr to 2.0 Torr, Fig. 9(b). The increase of pressure leads to the high collision frequency of electron with other particles, which prevent it from reacting with specific species to create neutral radicals required in the deposition process. The influence of the input power of plasma source on the deposition rate is disclosed in Fig. 10(a) and 10(b), where the deposition rate evidently grow with the input power. The increase in input power employed in the plasma source results in the elevation of electron temperature, where αs rises from 5.7 Å/s to 12 Å/s provided the power is doubled from 500 W to 1000 W. The growth in input power boosts the electron regeneration process and is favorable to the production of neutral radicals essential to the deposition process. However, the film deposition can suffer from non-uniformity due to input power increase as a consequence of high electron penetration close to the substrate without enough time diffusing in space. Fig.11 (a) suggests that severe problem in the uniformity of film may occur for very small chamber height. The raise of chamber height can provide smooth distribution of film thickness and also has positive impact on the average deposition rate, Fig. 11(b), where the chamber height of 12cm can generate an increase of 15% in the average deposition rate than the 4c m

one. Large chamber height can offer more space and time for chemical reaction and diffusion process of species, which significantly improve the evenness of film and also explains the smoothing effect through the height growth. The temperature of substrate is another nontrivial parameter to increase the deposition rate, Figure 12(a). An almost linear dependence of average deposition rate on the temperature of substrate is predicted in the numerical simulations, Fig. 12(b). This is principally owing to the high reaction rate on the substrate surface accompanied with the rise of temperature. The growth of temperature by 200 K can yield an gain of 3 Å/s. Nevertheless, the enhancement also brings some negative influence on the film smoothness. The arrangement of reducing the distance between two plasma orifices can merge the two peaks in the deposition rate curve into a single one, Fig.13 (a), which is disadvantageous from the viewpoint of film uniformity. But the increase in this distance can also suffer from the problem of smoothness because low plasma density in the middle region between plasma orifices. Thus an appropriate distance can be optimized through the numerical simulations. In addition, small distance between plasma orifices can promote the average deposition rate through very high plasma density favoring the deposition process, Fig.13 (b). 5. Conclusions The mathematical model for the afterglow-type PECVD process of argon plasma, including fifty-two reactions in space, and ten reactions on substrate, is proposed, where nineteen species are considered in the three-dimensional modeling with hexahedral grids. All governing equations are discretized by using a finite volume method and solved by an iterative approach. The variables, such as position, velocity and moment, are expressed by Cartesian coordinates, where the non-staggered variable arrangement is adopted to define field variables, where the parallelized numerical scheme is implemented via MPI library. The grid-dependence study of the proposed method is first conducted to verify the prediction accuracy of method and the employed approach is proven close to second-order accurate. The numerical prediction is also validated with the result of a previously published literature. The numerical predictions indicate that argon, hydrogen, and silane are the main background gases in the process chamber. In addition, neutral radicals, followed by ions, have small

species density in the investigated PECVD process. Numerical simulations are carried out to investigate the influence of geometrical and processing parameters on the deposition rate of the substrate. The numerical predictions suggest that two-dimensional calculations give a minor over-prediction (about 0.5Å/s) when compared with three-dimensional results, which give an average value ranging from 5 Å/s to 10 Å/s under the studied conditions. A plateau region with double peaks is found to describe the typical deposition behavior on the substrate for a process chamber equipped with two parallel rows of plasma orifices. The increase of system pressure results in low deposition rate due to the high collision frequency of electron with other particles. The deposition rate improves but the smoothness of film thickness deteriorates with the growth of input power, which is explained by the substantial increase of electron temperature. Large chamber height can improve the evenness of film, as well as the deposition rate. The average deposition rate is found quasi-linearly dependent on the temperature of substrate, where a temperature increase of 200 K can produce an improvement of 3 Å/s in the average deposition rate. The decrease in the distance between plasma orifices brings some negative influence on the film smoothness but apparently enhances the rise of average deposition rate.

References [1] Franz G. Low pressure plasmas and microstructuring technology. Springer; 2009. [2] Bruno G, Capezzuto P, Madan A. Plasma deposition processes of amorphous silicon-based materials. Academic Press; 1995. [3] Fridman A. Plasma chemistry. Cambridge University Press; 2008. [4] Lieberman MA, Lichtenberg AJ. Principles of plasma discharges and materials processing. Wiley; 2005. [5] Park, JH. Deposition of coatings by PECVD. Department of Electrical Engineering, Auburn University; 2003. [6] Intranuovo F, Sardella E, Rossini P, Agostino R, Favia P. PECVD of fluorocarbon coatings from hexafluoropropylene oxide: glow vs. afterglow. Chemical Vapor Deposition 2009,15:95-100. [7] Mitchell RR, Young RM, Partlow WD, Bevan MJ, Chantry PJ, Kline LE. Silicon nitride deposition using N2-rare gas radio-frequency afterglows. IEEE International Conference on Plasma Science 1990. [8] Meikle S, Nakanishi Y, Hatanaka Y. Thinfilm deposition in the afterglow of N2 and H2 microwave plasmas. J. Vac. Sci. Technol. A 1991,9:1051-1054. [9] Grigorian GM, Kochetov IV. Preparation of carbonitride films in the active and afterglow phases of a glow discharge. Plasma Physics Report 2013;39:412-419. [10] Lide DR. Handbook of Chemistry and Physics. Taylor and Francis Group LLC; 2009. [11] Meeks E, Larson RS, Ho P, Apblett C, Han SM, Edelberg E, Aydil ES. Modeling of SiO2 deposition in high density plasma reactors and comparisons of model predictions with experimental measurements. Journal of Vacuum Science and Technology A: Vacuum, Surfaces and Films 1998;16:544-563. [12] Kushner MJ. A model for the discharge kinetics and plasma chemistry during plasma enhanced chemical vapor deposition of amorphous silicon. Journal of Applied Physics1988;

63:2532-2551. [13] Bird RB, Stewart W, Lightfoot EN. Transport Phenomena, Second Edition. Wiley; 2001. [14] Krzhizhanovskaya VV. A virtual reaction for simulation of plasma enhanced chemical vapor deposition. Ph.D Thesis, The Institutional Repository of the University of Amsterdam; 2008. [15] Raizer YP. Gas discharge physics. Springer; 1991. [16] Salabas A, Gousset G, Alves LL. Two-dimensional fluid modelling of charged particle transport in radio-frequency capacitively coupled discharges. Plasma Sources Science and Technology 2002;11:448-465. [17] Perrin J, Leroy O, Bordage MC. Cross-sections, rate constants and transport coefficients insilane plasma chemistry. Contributions to Plasma Physics 1996;36:3-49. [18] Sakiyama Y, Graves DB. Neutral gas flow and ring-shaped emission profile in non-thermal RF-excited plasma needle discharge at atmospheric pressure. Plasma Sources Science and Technology 2009;18:1-11. [19] Hsu CC. Diagnostic Studies and Modeling of Inductively Coupled Plasma. PhD Thesis, University of California; 2006. [20] Denysenko IB, Ostrikov K, Xu S, Yu MY, Diong CH. Nanopowder management and control of plasma

parameters

in

electronegative

SiH4 plasmas.

Journal

of

Applied

Physics

2003;94:6097-6107. [21] Kessels WMM, Van De Sanden MCM, Severens RJ, Schram DC. Surface reaction probability during fast deposition of hydrogenated amorphous silicon with a remote silane plasma. Journal of Applied Physics 2000;87:3313-3320. [22] Amanatides E, Stamou S, Mataras D. Gas phase and surface kinetics in plasma enhanced chemical vapor deposition of microcrystalline silicon: The combined effect of power and hydrogen dilution. Journal of Applied Physics 2001;90:5786-5798. [23] Hsin WC, Tsai DS, Shimogaki Y. Surface reaction probabilities of silicon hydride radicals in SiH4/H2 thermal chemical vapor deposition. Industrial and Engineering Chemistry Research

2002;41:2129-2135. [24] Khosla PK, Rubin SG. A diagonally dominant second-order accurate implicit scheme. Computers & Fluids 1974;2:207-209. [25] Chau SW. Numerical investigation of free stream rudder characteristics using a multi-block finite volume method. PhD Thesis, Universität Hamburg; 1997. [26] Chau SW, Hsu KL. Modeling steady axis-symmetric thermal plasma flow of air by a parallelized magneto-hydrodynamic flow solver. Computers & Fluids 2011;45:109-115. [27] Pham QT. Modeling of gas reaction in argon plasma enhanced chemical vapor deposition process. Master thesis, National Taiwan University of Science and Technology; 2012. [28] Richardson LF. The Approximate arithmetical solution by finite differences of physical problems involving differential equations with an application to the stresses in a masonry dam. Transaction of Royal Society London 1910;210:307-357.

Table 1. Fluid properties of three main species [10].

Gas SiH4 H2 Ar

η cp κ M [kg/(m·s)] [J/(mol·K)] [W/(K·m)] [kg/kmol] 1.08×10-5 42 1.918×10-2 32 -6 8.6×10 21 1.683×10-1 2 -5 -2 2.09×10 20 1.636×10 40

Table 2. Reaction mechanism among Ar plasma, SiH4 and H2.

No.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36

Reaction (Electron - impact reactions with argon) E + Ar → E + + E + Ar → Ar+ → Ar+ + E + Ar* (Metastable reaction) → Ar + Ar* + Ar* (Electron - impact reactions with silane) → SiH4 + E + SiH4 → SiH4 + E + SiH4 → SiH3 + E + SiH4 E + SiH4 → SiH2 + → SiH3+ + E + SiH4 → SiH2+ + E + SiH4 + → SiH + E + SiH4 → Si+ + E + SiH4 E + SiH4 → H2+ + → H+ + E + SiH4 (Electron - impact reactions with hydrogen) → H2 + E + H2 → H2 + E + H2 E + H2 → H2 + → 2H + E + H2 E + H2 → H2+ + E + H → H + + E + H → H+ (Other electron - impact processes) → SiH3+ + E + SiH3 → SiH2+ + E + SiH3 → SiH2+ + E + SiH2 + E + SiH2 → SiH + + E + SiH → SiH+ + E + SiH → Si+ (Neutral radical reactions) → H2 + H + SiH4 SiH3 + SiH3 → SiH2 + → Si + SiH2 + M → SiH + SiH2 + H (Excitation transfer and radiactive decay) → H + Ar* + H2 → SiH3 + Ar* + SiH4 → SiH2 + Ar* + SiH4 * → SiH2 + Ar + SiH3 → SiH + Ar* + SiH2

Ref.

Ar* 2E 2E Ar+ E E H 2H H H2 H2 2H2 SiH SiH2

[11] [11] [11] +

+ + + + + + + +

E

[11]

E E 2E 2E H + 2E H + H +

[11] [11] [11] [11] [11] [11] [11] [11] [11] [11]

2E 2E 2E

E E E E 2E E 2E E H 2E H 2E H SiH3 SiH4 H2 H2 H H H H H

[11] [11] [11] [11] [11] [11] [11]

+

2E

+

2E

+

2E

+

M

+ + + + +

Ar Ar H + Ar Ar

[11] [11] [11] [11] [11] [11] [11] [11] [11] [11]

Ar

[12] [12] [12] [12] [12]

37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52

Ar* + SiH → Si + (Ion-neutral radical reactions) + Ar → Ar+ H2 + Ar+ + H2 → H2 + + + SiH4 → SiH3+ H + + SiH4 → SiH3+ H2 + + SiH4 → SiH3+ SiH2 + Ar + SiH4 → Ar + SiH4 → Ar Ar+ + SiH4 → Ar Ar+ (Other neutral-radical reactions) → Si2H4 SiH + SiH4 Si2H4 + H2 → SiH2 Si2H4 + SiH4 → Si2H6 → SiH4 H + Si2H6 → SiH2 Si + Si2H6 (Collisions with disilane reactions) → SiH4 + E + Si2H6 → Si2H4 + E + Si2H6

H

+

+ + + + + + + +

H2 Ar H2 H2 SiH3 SiH3 SiH2+ SiH+

+ + + + +

H SiH4 SiH2 SiH3 Si2H4

SiH2 H2

+ +

Ar

[12]

+

H

+ + +

H H2 H +

H2

[11] [11] [11] [11] [11] [11] [11] [11] [11] [11] [11] [11] [11]

E E

[12] [12]

Table 3. The collision diameter and Lennard-Jones parameter [14]. Species

σ(Å)

Ar H H2 Si SiH SiH2 SiH3 SiH4 Si2H4 Si2H6

3.44 2.50 2.94 3.84 3.66 3.80 3.94 4.08 4.60 4.82

ε K

(K)

120 30 36 100 96 133 170 207 312 301

Table 4. The mobility constant of ion in gas. j Ar SiH4 H2

c jj+ (Torr·m2 /V·s) Ref. 0.14 0.21 0.88

[16] [16] [16]

Table 5

The mobility constant of electron in gas. j Ar SiH4 H2

cej (Torr·m2/V·s) Ref. 33 10 37

[15] [17] [15]

Table 6. The reaction mechanism of silicon thin film deposition. No.

Reaction

Ref.

1

SiH3 +

Si(S)

→

SiH(S)

+

Si(B)

2

SiH3 +

SiH(S)

→

SiH4

+

Si(S)

3

SiH2 +

Si(S)

→

Si(S)

+

Si(B)

+ H2

[23]

4

SiH2 +

SiH(S)

→

Si(B)

+ SiH(S) + H2

[23]

5

SiH +

Si(S)

→

SiH(S)

+

Si(B)

[21]

6

SiH +

SiH(S)

→

Si(S)

+

Si(B) Si(B)

7

Si

+

Si(S)

→

Si(S)

+

8

Si

+

SiH(S)

→

Si(B)

+ SiH(S)

9

H

+

Si(S)

→

SiH(S)

10

H

+

SiH(S)

→

Si(S)

+ H2

[21] [22]

+ H2

[22] [23] [23] [21]

+

H2

[22]

Table 7. Principal dimensions of process chamber. Dimension da dg wx ha hg do wz Length (cm) 15 30 114 9.5 3.5 0.5 121

Table 8. Boundary conditions of field variables. Field Plasma Variables Orifice

Gas Ring

Substrate Chamber Chamber Outlet Wall ∂Te =0 ∂N

Te

Te,o

−

T

Tg

Tg

Ts

nAr

nAr,o

0

∂nAr =0 ∂N

nSiH4

0

nSiH 4 , g

FSiH 4

nH 2

0

nH 2 , g

FH 2

∂Te =0 ∂N

∂Te =0 ∂N

∂nAr =0 ∂N ∂nSiH 4 =0 ∂N ∂nH 2 =0 ∂N

∂nAr =0 ∂N

−

Tw

FSiH 4

FH 2

where N denotes the normal direction of the boundary.

Table 9. Comparison of predicted deposition rate in afterglow PECVD. Source

Deposition rate α (Å /s)

This study

7.89

[14]

2.5

Processing Dependence Dependence of α on P of α on Pw Conditions P = 1 Torr Pw = 750 W − + G = 300 sccm Ts = 523 K P = 1 Torr Pw =1250 W − + G = 60 sccm Ts = 523 K

where GGas denotes the total flow rate of feed gas (SiH4 and H2).

Fig. 1. Schematic illustration of the studied afterglow PECVD system.

Fig. 2 Definition of principal dimensions describing the process chamber.

Fig. 3. Average deposition rate for different grids.

Fig. 4. Comparison of relative discretization error for average deposition rate.

Fig. 5. Density distribution of main gases.

Fig. 6. Electron temperature profile and density distribution of Ar+ and Ar*.

Fig. 7. Density distribution of neutral radicals.

Fig. 8. Density distribution of ions.

(a)

(b) Fig. 9.

(a) Comparison of deposition rate along the centerline at different P. (b) The dependence of average deposition rate on P.

45

(a)

(b) Fig. 10.

(a) Comparison of deposition rate along the centerline at different Pw. (b) The dependence of average deposition rate on Pw.

46

(a)

(b) Fig. 11.

(a) Comparison of deposition rate along the centerline at different ha. (b) The dependence of average deposition rate on ha.

47

(a)

(b) Fig. 12.

(a) Comparison of deposition rate along the centerline at different Ts. (b) The dependence of average deposition rate on Ts.

48

(a)

(b) Fig. 13.

(a) Comparison of deposition rate along the centerline at different da. (b) The dependence of average deposition rate on da.

49