Global model simulations of low-pressure oxygen discharges

Abstract

We use a global model (volume averaged) to study plasma discharges in molecular oxygen gas in the 1-100 mTorr pressure range. This model determines densities of positive ions O2+ and O+ , negative ion O- , electrons, ground state O2 and O atoms, and metastables O2(a¹deltag) and O(¹D), and electron temperature as function of gas pressure and input power, for a cylindrical discharge. We apply the model to O2 discharges and the results are compared to a particle-in-cell simulation (PIC), experimental data and a volume-averaged global model developed at the University of California at Berkeley. We find that the total positive ion density increases with pressure at low pressures (up to approximately 30 mTorr), and decreases at higher pressures. The electronegativity decreases with increased power and increased pressure as predicted by the global models presented in the literature. The predictions for electron temperature are also in agreement with these models. However, there is a discrepancy betweeen these global models and PIC simulations and experimental data, for 20 and 40 mTorr cases, concerning electronegativity calculations. PIC simulations yield much higher electronegativities. There are strong indications that this is due to the assumption of Maxwellian electron energy distribution functions in the global model, while in the PIC simulations this is clearly not the case.

Global model; Low pressure oxygen discharges; Particle-in-cell simulation

REGULAR ARTICLES

Global model simulations of low-pressure oxygen discharges

Geraldo Roberson^I; Marisa Roberto^I; John Verboncoeur^II; Patrick Verdonck^III

^IDepartamento de Física, ITA, S. José dos Campos, SP, 12228-900, Brazil

^IIDepartment of Nuclear Engineering, University of California, Berkeley, CA, 94720, USA

^IIIIMEC Kapeldreef 75, 3001 Leuven, Belgium

ABSTRACT

We use a global model (volume averaged) to study plasma discharges in molecular oxygen gas in the 1-100 mTorr pressure range. This model determines densities of positive ions O₂⁺ and O ⁺ , negative ion O ^- , electrons, ground state O₂ and O atoms, and metastables O₂(a¹Dg) and O(¹D), and electron temperature as function of gas pressure and input power, for a cylindrical discharge. We apply the model to O₂ discharges and the results are compared to a particle-in-cell simulation (PIC), experimental data and a volume-averaged global model developed at the University of California at Berkeley. We find that the total positive ion density increases with pressure at low pressures (up to approximately 30 mTorr), and decreases at higher pressures. The electronegativity decreases with increased power and increased pressure as predicted by the global models presented in the literature. The predictions for electron temperature are also in agreement with these models. However, there is a discrepancy betweeen these global models and PIC simulations and experimental data, for 20 and 40 mTorr cases, concerning electronegativity calculations. PIC simulations yield much higher electronegativities. There are strong indications that this is due to the assumption of Maxwellian electron energy distribution functions in the global model, while in the PIC simulations this is clearly not the case.

Keywords: Global model; Low pressure oxygen discharges; Particle-in-cell simulation

I. INTRODUCTION

Oxygen discharges have been applied to numerous applications in plasma processing such as ashing of photoresist, removing polymer films and oxidation or deposition of thin film oxides. Several models have been proposed in order to study these discharges. A macroscopic model of parallel-plate RF oxygen discharge was made by Misium [1] and applied to low-pressure O₂ plasmas. Lichtenberg, Vahedi and Lieberman [2] developed a macroscopic analytical model for electronegative plasma discharges comparing the results with particle-in-cell simulations. Lee, Graves, Lieberman and Hess [3] studied the gas-phase kinetics and plasma chemistry of high density oxygen discharges in the low and high pressure regimes in cylindrical geometry. A global model of high-density plasmas discharges in molecular gases using an inductive cylindrical discharge was proposed by Lee and Lieberman [4], showing that for molecular gases, the electron temperature is no longer only a function of pressure and the reactor geometry, but also strongly depends on the plasma composition. Shibata, Nakano and Makabe [5] investigated the space-time structure of O₂ RF discharges by using the relaxation continuum model. They showed that atomic oxygen, formed by dissociative electron impact in O₂, plays an important role in the RF structure through the detachment process. The effect of the electrode material was taken into account by the same authors in reference [6].

In this article, we examine a volume-averaged (global) model of the plasma chemistry of O₂ gas and compare this model to the particle-in-cell simulation (PIC) [7], the global model developed by Lieberman's group at the University of California at Berkeley [2-4,8,9,10-13] and experimental data [14]. PIC simulation presents some advantages over a volume-averaged model. While the global model assumes a Maxwellian distribution function for the particles in order to calculate the rate coefficients, the PIC model calculates these distributions, which are not exactly Maxwellian. From PIC, we can obtain the pressure-dependent profiles of each species. However, the PIC simulations take a very long time to run, especially for molecular gases and higher pressure cases ( > 100 mTorr), because it requires many particles to obtain sufficient statistics. As a consequence, much more information can be obtained. Despite the assumptions of the global model, including a Maxwellian distribution function for electrons and spatial uniformity of the density profiles with the pressure, it can give an indication of how one parameter depends on another and which reactions among the species are important.

II. THE GLOBAL MODEL

We assume a cylindrical chamber of radius R and length L. A steady flow Q of neutral species is introduced through the inlet. The power is assumed to be deposited uniformly into the plasma bulk. We are using a simplified global model described in reference [4]. The assumptions are listed below.

All densities are assumed to be volume averaged;

A maxwellian electron energy distribution function is assumed;

The discharge gas and ion temperatures are assumed to be constant, at 600 K, irrespective of the discharge conditions;

For an electronegative discharge, the electron density n_e is assumed to be uniform throughout the discharge except near the sheath edge. The negative-ion density n _- is assumed to be parabolic, dropping to zero at the sheath edge, the positive-ion density is n ₊ = n_e+ n _- , with n ₊ =n_e=n _{+ s} at the sheath edge. An appropriate interpolation between electroposive and electronegative profiles is used, as described in reference [4].

We have considered the set of the reactions shown in table I, which includes the neutral species O₂ , O, the positive ions and O ⁺ , the negative ion O ^- and the metastables O₂(a¹Dg) and O(¹D). The rate coefficients are in references [2,3,4,8,9,10]. We followed reference [11] to decide which reactions would be important to include in this model. The density of the metastable singlet delta state O₂(a¹D_g) is roughly 2-11% of the total O₂ density in the pressure range of 1-100 mTorr. This state is mainly created by electron impact excitation and lost by electron impact dissociation.

Table II

The rate of generation and annihilation of each species is determined by the plasma chemistry. For each neutral and charged species, a particle balance equation is generated which accounts for all dominant generation and destruction processes and diffusion losses to the wall. Steady state particle balance equations are written for O₂, O, O(¹D), O₂ (a¹Dg), , O+ and O ^- . The general form of the particle balance equations is flow in + rate of generation = rate of annihilation + pumping loss - as shown in equations (1) to (4).

The particle balance equations for neutral species are given by :

The flow through the reactor has been characterized by a nominal residence time t. This residence time, and the pumping rate k_p = 1/t, has been determined from the flow-rate, pressure, and the reactor volume. The O₂ ''Source'' term in equation (1) is calculated as follows: 4.4836x10¹⁷ x Q in units of molecules/s, where Q is the flow rate in cubic centimeters per minute at STP.

For our model, we consider the production of the ion only through electron-neutral collisions. We consider O ^- as the only negative ion species in the plasma and further we assume that no negative ions are lost to the walls. The justification for this assumption is that the negative ions are essentially trapped in the bulk plasma because of the high positive potential of the plasma with respect to the chamber walls. The particle balances for , O ^- and O ⁺ are written as

Since the bulk of the plasma is essentially neutral, we have

The power balance equation, which equates the change in electron thermal energy to the absorbed power P_abs minus power losses due to elastic and inelastic collisions:

where V is the volume of the discharge, and P_abs can be written as

where r is the number of positive ion species generated in the system, i.e., for O₂, r=2, for generation of both and O ⁺ . P_ev is the electron energy loss due to to all electron-neutral collision processes in the volume, P_ew is the electron energy loss to the walls and P_iw is the ion energy loss to the walls. The total power lost in the volume is

where n_iz,i is the total ionization frequency for production of ith ion (i=and O ⁺ ) and E_L,i is the total energy lost by collision processes in the volume. Equation (11) can be put in the form, as in reference [4]

where u_B,i (i = and O⁺ ) is the Bohm velocity given in a more general form u_B,i = (eT_e (1+ a)/(m_i(1+ag)))^{1 / 2}, a = n_-/n_e and g = T_e/T_i.K_rec,i are the rate coefficients for ion - ion recombination (k₁₇ and k₃₄) and V is the volume of the reactor. A_eff,i is the effective area for ion loss, which is given as

The expression for the ratio of the sheath edge density n_is to the bulk average density n_i is the same as used in reference [4], which is given by

and

J₁ is the cylindrical Bessel function and D_a,i is the ambipolar diffusion coefficient for species i. For the volume averaged discharge

where D_i = eT_i / (m_in_i), where n_i is the total momentum transfer collision frequency.

The energy E_T,i = E_iw + E_ew + E_L,i, where E_iw = 3T_e [4], which is the ion kinetic energy lost per ion lost to the walls and E_ew = 2T_e [4] is the electron kinetic energy lost per electron lost to the walls. E_L,i is the collisional energy lost due to all electron-neutral collision processes in the volume, defined as

where N_s,i is the neutral species number that generates the i^th ion. For O ⁺ , N_s,i = 3 (O₂, O and O(¹D) and for O_{2 +} , N_s,i = 1 (O₂). We do not take into account the generation of from O₂ (a₁Dg) because this generation is lower than this from ground state. The parameter n_iz,ij is the ionization frequency to generate the i^th ion from neutral species j and n_iz,i is the total ionization frequency. The sum over k includes all collision processes that do not produce positive ions and E_all is the threshold energy for these processes The n_elas,j is the elastic frequency collision. For the present study the system of first-order differential equations is allowed to reach a steady state.

III. RESULTS AND DISCUSSIONS

We apply the global model to a cylindrical discharge in a stainless steel chamber with radius R=2.52cm and length L=5 cm in order to compare these results to PIC simulations in reference [7]. We assume an applied power of 20, 30, 40 and 50 W, a flow rate Q=10 sccm and a neutral gas temperature of T_g = 600 K. Since the densities of the atomic oxygen and the metastable oxygen molecule are largely determined by their interactions with the wall, the chamber material and geometry have a decisive influence on the properties of the discharge. The wall sticking coefficients for wall reactions are in reference [9] and they are shown in Table I.

Figures 1a and 1b show O^- , electron, O⁺ and ion density dependence on the neutral pressure for 20 W and 50 W of applied power. One can see the same trends as those in references [9] and [12]. For low pressure the concentration of negative ions is lower than the electron concentration. However, at low power and high pressure, the effect of negative ions becomes more important. The reason can be found in the reactions of generation and loss of these ions. The generation is mainly through electron impact collisions with oxygen molecules in the ground state and metastable state (reactions 9 and 26 from Table I). Figs. 2a and 2b show the rate coefficients (k) and reaction rates (R), which is given by R=rate coefficient x density, related to generation of O^- (i.e. R9 = k9n_e n_O2, R26 = k26n_e n_O2 and R31 = k31n_e n_O2). Fig. 3 shows the same for loss of O^- . The loss terms are electron impact detachment, ion-ion neutralization and detachment by O (reactions 16, 17, 23 and 34 from Table I), since the negatively charged species are trapped in the plasma due to the high potential barrier at the reactor walls. One can see from Fig. 2 that reactions 26 and 9, except for very low pressures ( ~ 1 Torr) are dominant for the generation of O^- . At constant pressure, the negative ion density increases slightly with decreasing power. For low power the generation of O^- is higher, because the molecular oxygen density is higher due to the lower dissociation rate (see Fig. 5a). The loss rate of O^- is also lower for low power (see Fig. 3), since the overall positive ion density decreases with decreasing power (see Fig. 1b) and the generation of O₂ is higher than O for low power.

The behavior of the electron density as a function of gas pressure has been discussed in references [4] and [12]. According to them, it has to do with the balance between the effective loss area for positive ions and the energy loss per created electron-ion pair (Eq. 17). The effective loss area is a decreasing function of pressure since the ions are better confined at higher pressures. The energy loss factor increases with increasing pressure. The positive ion density, and hence the electron density, increases initially but decreases at a higher pressure. Thus the energy loss factor dominates the effective loss area at higher pressures. These trends can be verified in Fig. 1a.

The effect of input power and pressure on the positive ion densities (O₂+ and O ⁺ ) is presented in Fig. 1b. From this figure, we see that the total positive ion density increases with increasing input power and decreasing pressure (the result is more pronounced for O ⁺ than ). At a fixed pressure, the ion density increases as more power is applied to the system. The explanation is related to the behavior of the electron temperature as is shown in Fig. 4a, which is weakly dependent on power and decreases with increasing pressure. The electron temperature is determined by multiple particle balance equations which are strongly coupled. Therefore, T_e is a function of both input power and pressure, even though the dependence on power is much weaker than the pressure dependence. Besides, as pressure increases, the ion-ion neutralization becomes more important what contributes to the decrease of the overall positive ion density at higher pressures.

Figure 4b shows the electronegativity a as a function of neutral gas pressure. The electronegativity as predicted by the global model described in references [9] and [12], decreased with increasing power and increased with increasing pressure. As shown in Fig. 1b we obtained the same trends for the studied range of powers and pressures. These trends are directly explained from Fig. 1a.

Figures 5a and 5b show the neutral species, O₂ and O, and the metastables, O₂ (a¹Dg) and O(¹D) in the discharge. From Fig. 5a, we can see that for low power the oxygen molecule is the dominant neutral species. However, as power increases oxygen molecules are dissociated and therefore the O atoms become the dominant species. As pressure increases, the O density decreases for a fixed power, because the reaction with the walls (reaction 22) decreases as pressure increases. Fig. 6 shows this trend. One can see on this figure the rate coefficients (k21, k22, k28, k29 and k36) and reaction rates (R21, R22, R28, R29 and R36), which are the rate coefficients multiplied by the density of the lost species. The generation of O from reaction 36 decreases as pressure increases what can also explain the behavior of the oxygen atom concentration. As shown in reference [9], lower pressure favors the creation of the metastable oxygen atom O(¹D) ( < 10 mTorr) but this creation decreases with increasing pressure. The metastable O₂ (a₁Dg) increases as pressure increases as can be seen in Fig. 5b. It has a significant role in capacitive rf discharges [6]. According to results from literature [3] in the generation and loss reactions of O(¹D), one can see that the main generation reaction is the excitation of ground-state oxygen atoms, and the main loss mechanism is diffusion loss to the walls at low pressure and ionization at high pressures.

Hence, our global model presents trends which are similar to the global model presented in the literature. A global model which includes transitions from electropositive to electronegative regions, using profiles with non spatial uniformity has been developed by Lieberman's group at Berkeley [13]. They considered different kinds of profiles for different pressures. A result of this model is presented in Table III. It also shows the results obtained by S. Kim et al. [14] (20 mTorr case), experimental results obtained by Stoffels et al. [15] using an RF parallel-plate reactor and the results from the PIC simulation [7] for 20 and 40 mTorr cases. The PIC code has used the planar one-dimensional model to study discharges between parallel plates at 13.56 MHz. For electronegativity, a = n_{O -} /n_e, the average values of the densities were used. g For our global model a_avg is the same as a. The results are compared for an applied power of 20W.

As we can see from Table III, the results for average a from PIC simulation and Stoffels' experimental results are in a reasonable agreement for 20 and 40 mTorr and 20 W cases. However, the global model (from reference [13] and figure 1b) shows lower values of electronegativity. It seems from the global model that as pressure increases the electron density decreases [8,9,12]. It seems that the choice of realistic density profiles is crucial to obtain a better estimation of the electronegativity for low pressure cases. As shown from PIC simulations [7], for low pressure case, the profiles are not flat in the center of the discharge and the eedf's are not Maxwellian, as assumptions described in reference [12]. Reference [13] has modified the profiles, but still considered the eedf Maxwellian. A comparison between the curves for electron density (n_e), obtained by Stoffes [15] and obtained by Kim [13], shows different behaviors for them. While the experimental results show an ascendant curve up to ~ 60 mTorr and with low values of n_e for 10,20 mTorr, the global model shows few variations of n_e for these pressures. The behavior of electron density with pressure depends on the balance between the effective loss area for positive ions and the energy lost by electrons due to the collisions. They depend strongly on the choice of the density profiles and the eedf, which is used to calculate the rate coefficients. It seems at least for low pressures cases, an eedf non-Mawelliana could be used in the improved global model [13].

Table IV shows a comparison between the global model obtained by Gudmundsson [12] and our results, for R=15.2cm, L=7.6 cm and Q=50 sccm for 100, 300 and 500W of absorbed power. Our global model is similar to that of Lee and Lieberman [3]. Gudmundsson [12] has made some improvements including more reactions and self-consistent calculations of the plasma potential.

From Table IV we can see relatively good agreement between the results. The electronegativity obtained by Gudmundsson [12] is in good agreement for 300 W, less for 500 W and 100 W. The electron temperature in general is quite comparable. As we are considering different set of reactions we expect to find only the same trends when we compare both models. These results show how it is important to use a kinetic method such as PIC, which presents realistic profiles and can represent a non-Maxwellian electron energy distribution function, if one wants to obtain good quantitative results. However, the PIC method takes quite a long time to run, especially for pressures higher than 100 mTorr. For a 1GHz clock computer, it takes more than 15 days to run a 200 mTorr case considering the set of reactions presented in reference [7]. The global model is much faster; it takes a few minutes to reach the steady state. However, the general trends of most physical parameters can be adequately predicted by the global model presented in this paper. It also indicates which physical-chemical reactions may be neglected, very important data for reducing the computational load for more complete models such as PIC.

IV. CONCLUSIONS

A volume-averaged global model of electronegative oxygen discharges was developed. The model uses a power balance equation to account for energy deposited into the plasma and lost via collisions and particle flux. The particle densities are modeled via rate equations estimated from collision cross sections assuming Maxwellian electron energy distribution functions. Our results have the same trends as the ones obtained through simulations with a similar global model. We also looked at the results from an improved global model, experimental results and a PIC simulation and all these show that the presented model yields the correct trends and relatively good quantitative results. Hence it may give good indications for which reactions should be included and which may be neglected when performing other - e.g. PIC - simulations. We also notice that the global models do not yield good agreement when compared to PIC simulation and experimental data, for the studied cases. This may be due to the choice of density profiles which do not agree well with experimental data. However, the improved global model [13] modified these profiles and we notice despite this the PIC simulation and experimental data still do not agree with the results of reference [14]. It seems that to calculate the proper rate coefficients it is necessary to use the non-Maxwellian electron energy distribution function, at least for the low pressure cases. In spite of this the global models can still get many trends right at low pressures, although the quantitative agreement with experiments and more complete simulations is worse at low pressures.

Acknowledgements

The financial support of CAPES and FAPESP is acknowledged. One of the authors (M. Roberto) would like to thank J.T. Gudmundsson for helpful discussions.

[14] S. Kim, M.A. Lieberman, A.J. Lichtenberg, Personal communication.

Received on 12 April, 2007

• [1] G.R. Misium, J. Vac. Sci. Technol. A 8, 1642 (1990).
• [2] A.J. Lichtenberg, V. Vahedi, and M.A. Lieberman, J. Appl. Phys. 75, 2339 (1994).
• [3] C. Lee, D.B. Graves, M.A. Lieberman, and D.W. Hess, J. Electrochem. Soc. 141, 1546, (1994).
• [4] C. Lee, M.A. Lieberman, J. Vac. Sci. Technol. A 13 , 368 (1995).
• [5] M. Shibata, M. Nakano, and T. Makae, J. Appl. Phys. 77, 6181 (1995).
• [6] M. Shibata, M. Nakano, and T. Makae, J. Appl. Phys. 80, 6142 (1996).
• [7] M. Roberto, M. J. Verbonceour, P. Verdonck, and E. Cizzoto, Proceedings of SBMICRO 2006. ECS Transactions, 4, 1, 563, (2006).
• [8] J.T. Gudmundsson, A.M. Marakhtanov, K.K. Patel, V.P. Gopinath, and M.A. Lieberman, J. Phys. D: Appl. Phys. 33, 1323 (2000).
• [9] J.T. Gudmundsson, I.G. Kouznetsov, K.K. Patel, and M.A. Lieberman. J.Phys.D: Appl. Phys. 34, 1100, (2001).
• [10] T. Kimura, A.J. Lichtenberg, and M.A. Lieberman, Plasma Sources Sci. Technol. 10, 430, (2001).
• [11] J.T. Gudmundsson, Report RH-17-2004, University of Iceland.
• [12] J.T. Gudmundsson, J. Phys. D: Appl. Phys. 37, 2073, (2004).
• [13] S. Kim, M.A Lieberman, A.J. Lichtenberg, and J.T. Gudmundsson. J. Vac. Sci. Technol. A, 24, 6, 2025, (2006).
• [15] E. Stoffels, W.W. Stoffels, D. Vender, M. Kando, G.M.W. Kroesen, and F.J. Hoog, Physical Review E, 51, 3, 2425 (1995).

Publication Dates

• Publication in this collection
  17 July 2007
• Date of issue
  June 2007

History

• Received
  12 Apr 2007