Electromagnetic Modeling of the Propagation Characteristics of Satellite Communications Through Composite Precipitation Layers Science and Technology, 8 (2003) 47-54 © 2003 Sultan Qaboos University A Monte Carlo Method for Low Pressure Radio Frequency Discharges Lahouaria Settaouti and Abderrahmane Settaouti Department of Electronics, University of Sciences and Technology, P.O.Box 1505 Oran, EL-M’NAOUAR, Algeria, Email: settaouti@yahoo.fr طريقة مونت كارلو لتفريغ تردد راديو تحت ضغط منخفض هوارية ستاوتى وعبدالرحمن ستاوتى إن االهتمام بالتفريغ تزيد الراديو إزداد كثيراً في األعوام األخيرة وذلك ألهميته في تقنيات اإللكترونيات الصغيرة : خالصة م بخصوصية التفريغ في غاز موجب التكهرب الذي يكثر استعماله في التطبيقات التكنولوجية في هذا الدقيقة وخاصة اإلهتما البحث نقدم بعض الخصوصيات المفصلة للبالزما تردد الراديو الموجودة عن طريق محاكاة مونت كارلو في غاز موجب . التكهرب ABSTRACT: There is increasing interest in glow discharges because of their importance to a large number of application fields, like the microelectronics industry, flat plasma display panel technology, the laser and light industry and analytical spectrochemistry. To improve the capabilities of rf glow discharges, a good understanding of the discharge physics is highly desirable. The typical calculated results include the radio frequency (rf) voltage, the electrical field distribution, the density of argon ions and electrons, the electron energy distribution function and information about the collision processes of the electrons with the Monte Carlo model. These results are presented throughout the discharge axis and as a function of time in the rf cycle. Moreover, we have investigated how many rf cycles have to be followed before a periodic steady state is reached. KEYWORDS: Monte Carlo Simulation, Glow Discharge, Radio Frequency, Plasma. 1. Introduction Glow discharges are being used in many fields of application. They serve extensively as plasma processing devices in microelectronics, such as for ion etching, thin film deposition, and plasma treatment of surfaces, and they also find application as atomization - excitation - ionization sources in analytical chemistry (Bogaerts et al 1995), (Choi et al 1992), (Wang et al 1997). A wide range of parameter space can be chosen to operate the plasma processes. By parameters here we mean those quantities that can be varied by changing not only the position of some external "knobs" (for example, pressure, current or voltage applied, and frequency), but also other quantities that cannot be varied as easily, such as electrode spacing and electrode area ratio. To attain better results in these application fields, a quantitative understanding of the glow discharge is required. The application of gas discharges in processing of materials for microelectronics has reached the point where it is not possible to achieve improvements purely on an empirical basis. Theoretical investigations have recently reached a new level that promises to be a sufficiently good foundation for predictive models (Nanbu, 2000), (Yonemura et al 2001). However, in order to accomplish that goal, further experimental studies of the discharges and a wide range of data are required. On the other hand, plasma technological devices and processes are complex and it is necessary to have constant monitoring. A thorough physical understanding of such discharges, in particular of the spatial dependencies of plasma parameters as charged particle 47 LAHOUARIA SETTAOUTI and ABDERRAHMANE SETTAOUTI densities and fluxes, electric field and electron temperature, in particular in the plasma and electrode sheath, is therefore desirable. Theoretical study of glow discharge in the past has been mostly qualitative or semi- quantitative in nature. Rigorous approach to the problem requires functional solutions of the discharge equations which are difficult to obtain in analytic forms. An alternative is to use a large computer and solve the discharge equations numerically by means of some algorithms. Consequently, there has been significant interest in developing models for these plasma based processes (Gogolides, 1997), (Denpoh et al. 2000). A validated process model has several benefits ranging from providing an understanding of mechanisms and in-process control to aid in equipment design. A current challenge in designing plasma processing tools is the development of computer models of rf discharges that can accurately describe non-equilibrium charged particle transport and plasma chemistry, yet execute quickly enough to make more realistic multidimensional simulations feasible. As in most areas of computational plasma and discharge physics, it is possible to describe the physical processes in the discharge with a variety of approaches. Generally one uses either a fluid treatment (Franklin, 1999), a kinetic scheme (Loffhagen et al 2002), or a hybrid combining aspects of the fluid and kinetic schemes (Garrigues et al. 2000). Briefly, kinetic schemes solve for at least some of the species velocity distribution functions in both time and space whereas fluid models assume distribution functions. Kinetic schemes are therefore intrinsically higher dimensional than fluid models since kinetic schemes involve solutions in velocity space as well as physical space and time. All other things being equal, kinetic schemes are generally much more costly computationally than fluid models because they provide more information. The Monte Carlo simulation has become increasingly important as a simulation tool particularly in the area of low temperature plasma physics. The motivation for the development of the computer simulation described in this paper has its basis in the need to quantitatively understand microelectronic plasma processing. 2. Description of the model Monte Carlo methods as applied to gas discharge problems involve evaluating the percentage of a given species of particles emanating from a given source, after experiencing energy loss and gain and terminating in defined categories. Its technique is realistic because a large number of particles is followed from the source throughout their life history. The fundamental physical concept is the mean free path or the mean free flight time, with the collision equations being formulated subject to the conditions such as the number density, and electric field intensity. The computing time for the Monte Carlo technique depends upon the number of test electrons released from the space of the interval and the number of collisions occurring while each electron travels the distance from the cathode to the anode. A simple group of electrons, typically 100-500, is randomly seeded between the electrode plates. The initial distributions of simulation particles are assumed to be Maxwellian of 5 eV for electrons. The three dimensional motion of electrons between two successive collisions under the electric field is solved numerically with use of Euler method. Electron flight time between two successive collisions is determined by the electron collision cross-section with Ar neutral particles. The scattering after all collisions is assumed to be isotropic. The secondary electron emission coefficient due to electrons flowing to the two electrodes is assumed to be zero. The distribution of the electric field is calculated from the solution of the one dimensional Poisson's equation with the use of the electron and ion densities at the end of each time step. Electron trajectories using the Monte Carlo method (Helin et al 1996, Settaouti et al 1999, Nanbu, 2000, Settaouti et al., 2001) for an rf parallel plate discharge in an electropositive gas (Ar) while oscillating the applied electric field providing a time and spatially dependent electric field over many rf cycles are computed. The electrode system is assumed to consist of parallel plates with a diameter that is larger than the electrode separation. Every electron, during its transit in the gas, performs a succession of free flights punctuated by elastic or inelastic collisions with atoms of gas defined by collision cross-sections. During the successive collisions 48 A MONTE CARLO METHOD FOR LOW PRESSURE RADIO FRQUENCY DISCHARGES for every electron, certain information (velocity, position, etc.) is stored in order to calculate, from appropriate sampling methods, transport coefficients and macroscopic coefficients. The probability of collision and the nature of collision are simulated by comparison with computer generated random numbers. The probability of collision in the time step ∆T is       ∆− −= mT T P exp1 (1) Tm is the mean flight time of an electron moving with a velocity v(ε), ( ) ( )εε vtNQ Tm . 1 = (2) in which v(ε) is the drift velocity of an electron. where Qt is the total collision cross-section in m2, N the gas number density and ε the electron energy in eV. Qt is a function of the electron energy. At t = 0 an electron observes a free flight time with a randomly selected angle of entry depending on the cosine distribution. The collision is simulated by comparing P with R1 at the end of each step, where R1 is a random number uniformly distributed between 0 and 1. When the electron undergoes a collision 1exp1 RT T P m ≥               ∆− −= (3) The nature of the collision is determined in the following way: P2,j is the probability that collision process j takes place, j = 1, 2, 3,... n, including momentum transfer, excitation, and ionization collisions. This leads to t i j Q Q P =2 ∑ = 12 jP Nj PPPP 222,21,2 ≤≤≤ jj PPPRPPP ,22,21,221,22,21,2 ....++≤≤++ − (4) determines the j-th collision process. The sum of the fractional probabilities is equal to unity, and the interval [0,1] is divided into segments with lengths corresponding to these fractional probabilities. A new random number R2 between 0 and 1 is generated, and the interval into which this random number falls, determines the type of collision that occurs. The new energy and direction after the collision depends upon the type of collision: for excitation, the excitation threshold energy E is given by: excEEE −= 0 , where Eexc is the excitation threshold and E0 is the electron energy before collision. For ionization, the total energy before collision is divided between the primary (original) electron and the secondary electron created in the ionization collision. For elastic collisions, the new kinetic energy of the electron is calculated by     −−= )cos1(210 χεε M m (5) which is deduced from the hard sphere model. χ is the scattering angle of the electron after collision, where m and M are the masses of electron and argon atom, respectively. After a collision the angles are determined by isotropic distribution. The assumption of isotropic scattering is a hypothesis which is consistent with the first order theory if, at the same time the total collision cross section is replaced by the momentum transfer cross section. The total number of electrons in the gap increases over many orders of magnitude, hence scaling is necessary to limit the number of simulation particles. When, the total number of simulation particles, exceeds the maximum allowable number of simulation particles, which is specified in the program input data, a statistical subroutine is introduced to choose a new group of larger particles to represent the old larger group 49 LAHOUARIA SETTAOUTI and ABDERRAHMANE SETTAOUTI of smaller particles. The subroutine contains a weighing of velocity distribution of the old group, so that the new group is equivalent in phase space to the old group. In this simulation, the steady state is defined as that in which the electron and ion densities at the center of the bulk region stop increasing, when the initial densities are at least one order of magnitude smaller than the final densities. 3. Results and discussion We present typical results of numerical simulations using a Monte Carlo method of rf plasmas sustained by an rf external source of 13.56 MHz and 200 V applied to the electrode at z = 0, in an argon model gas. The electrode system is assumed to consist of parallel plates with a diameter that is larger than the electrode separation. Electrode separation and gas pressure used are 4 cm and 0.3 Torr, respectively. The electron collision cross sections set of Ar employed is that reported by Delcroix (1960), Christophorou (1971), Huxley et al (1974) and Vossen et al (1978). It is seen from Figure 1 that the electric field in the first cycle oscillates with a period corresponding to the rf cycle. Figure 2 shows that the mean electron energy in the first rf cycle is maximum when the electric field for this cycle is maximum, minimum where the electric field is minimum, and is almost spatially constant with small oscillation with the rf external field. 0.2 0.4 0.6 0.8 1.0 0.8 1.6 2.4 3.2 4.0 -100 0 100 E (V/cm) Di sta nc e ( cm ) Time (cycle) 0.2 0.4 0.6 0.8 1.0 0.8 1.6 2.4 3.2 4.0 0 15 30 ε (eV) D ist an ce (c m ) Time (cycle) Figure 1. Temporal and spatial variations Figure 2. Temporal and spatial variations of the of the electric field for the first rf cycle. mean energy of electrons for the first rf cycle. The electric field in the first cycle oscillates with a period corresponding to the rf cycle. The electron density is low at the beginning of the avalanche growth but as time increases, the enhancement becomes greater, and oscillates with the change in the direction of the electric field (Figure 3). When the applied electric field is sufficiently large, the electron accelerates to a high enough kinetic energy that makes inelastic collisions between the electron and the background gas occur. As the electron proceeds from the cathode to the anode, driven by the applied electric field, it undergoes many ionizing collisions in each of which an ion and another electron are formed. Thus, we obtain an expanding cloud of electrons traveling toward the anode, known as electron avalanche, and a cloud of ions, almost stationary in the time scale of motion of the electron avalanche, remaining behind, and the space charge field distortion begins. 50 A MONTE CARLO METHOD FOR LOW PRESSURE RADIO FRQUENCY DISCHARGES 0.2 0.4 0.6 0.8 1.0 0.8 1.6 2.4 3.2 4.0 1 ne (cm -3) Distance (cm) Tim e ( cyc le) 0 1 2 3 4 -60 -30 0 30 60 E (V/cm) 16T/25 17T/25 18T/25 19T/25 Distance (cm) Figure 3. Temporal and spatial variations Figure 4.1. Electric field in argon discharge of the electron density for the first rf cycle. as a function of position for different times, and T is the rf cycle (third cycle). 0 1 2 3 4 -200 -100 0 100 200 300 E (V/cm) T/25 2T/25 4t/25 Distance (cm) 0 1 2 3 4 -1000 0 1000 E (V/cm) 18T/25 19T/25 20T/25 21T/25 Distance (cm) Figure 4.2. Electric field in argon discharge Figure 4.3. Electric field in argon discharge as a function of position for different times, as a function of position for different times, and T is the rf cycle (fourth cycle). and T is the rf cycle (fourth cycle). The space charge field distortion begins in the third cycle (Figure 4.1), and in the beginning of the fourth cycle (Figure 4.2), the field behind and ahead of the avalanche is enhanced, while in the bulk of the avalanche it is weakened by the space charge field. Figures 4.3 and 4.4 show, the spatial and temporal evolution of the electric field during the fourth rf cycle, when the plasma is established. There are two distinct regions with respect to the flow of electron energy in a parallel 51 LAHOUARIA SETTAOUTI and ABDERRAHMANE SETTAOUTI plate rf discharge; the bulk plasma and the sheaths. The electric field for the discharge is symmetric about the discharge mid-plane (In the parallel plate gap there are the sheath region, the bulk plasma and the sheath region). We also see from the Figures 4.3 and 4.4 that the electric field in the sheath regions, adjacent to the electrodes, oscillates with significant amplitude with period corresponding to that of rf cycle. Most of rf external electric field is absorbed in the sheath regions. As the cycle advances, electrons are pushed out of the right sheath and into the left one. The movement of electrons causes a modulation of sheath width. The sheaths are quite wide causing the discharge to behave capacitively. The apparent sheath thickness is 0.3 cm, (The sheath width was defined as the distance from the electrode to the edge region of the bulk plasma) and most of the voltage difference is concentrated in the sheaths. In the center and close to the electrodes the temporal modulation of electric fields is fairly well represented by a sinusoid with little dc bias.. The central bulk plasma is characterized by a weak electric field (few V/cm). At the times an electrode is cathodic, the plasma near the electrodes resembles the cathode fall of dc discharge in that there is a large voltage drop between the electrode and the bulk plasma in the center of the discharge. The calculated electric field (temporal and spatial variations) and sheath thickness compared with the measured values (Sato et al 1990), (Guo et al 1996), show good agreement. 0.2 0.4 0.6 0.8 1.0 0.8 1.6 2.4 3.2 4.0 -2000 -1000 0 1000 2000 E (V/cm) di sta nc e ( cm ) time (cycle) 0 1 2 3 4 100 101 102 103 104 105 106 107 108 109 1010 Distance (cm) ne (cm -3) Figure 4.4. Temporal and spatial variations Figure 5.1. electron density as a function of the electric field for the fourth rf cycle. of position for the fourth rf cycle. At the beginning of the growth of the avalanche, the electron density is low, and as time increases, the enhancement becomes greater and the electron density oscillate with the change in the direction of the electric field. When the plasma is established, electrons are contained in the plasma in part by the space charge fields. Depending on the phase of the applied potential, electrons move toward either electrode. Their movement affects the local field causing electrons to pile in the bulk sheath interface. The charged particle density profiles for argon are shown in Figures 5.1 and 5.2. From Figure 5.2 we see the positive ion density (due to the ionization) as a function of position. The central bulk plasma is characterized by a high density of electrons and positive ions, and the approximately equal densities of electrons and positive ions. The small peaks around the two sheath edges appear in the electron density. The positive ion density form an envelope around electrons, which move in and out of the sheath region. The agreement between, 52 A MONTE CARLO METHOD FOR LOW PRESSURE RADIO FRQUENCY DISCHARGES the calculated charged particle density compared with the measured values (Sato et al 1990), (Guo et al 1996), is good. In Figure 6, the profiles of the electron mean energy as a function of the distance show the different natures of the three regions in the discharge. The electron mean energy is low through the plasma bulk region and the relatively large modulations take place in sheath regions. The maximum in the electric field in the sheath regions create a maximum in the average electron energy shown in Figure 6. Here the secondary electron emission was assumed to be zero. Indeed the maximum in the electric field is close to the maximum in energy, but not exactly at the same position in space or time because we do not assume local equilibrium with the field. The average electron energy is low in time and space in the bulk region, when the electric field is minimum. 0 1 2 3 4 100 101 102 103 104 105 106 107 108 109 1010 Distance (cm) ne (cm -3) 0 1 2 3 4 0 2 4 6 8 10 ε (eV) Distance (cm) Figure 5.2. Positive ion density as Figure 6. The mean energy of electrons a function of position for the fourth as a function of position for the fourth rf cycle. rf cycle. 4. Conclusion We have presented typical numerical results obtained using the Monte Carlo simulation code, and have shown the detailed properties of rf plasmas in argon gas which is sustained by an rf external source of 13.56 MHz and 200 V. The simulation results show that there are two distinct regions with respect to the flow of electron energy in a parallel plate rf discharge, the bulk plasma and the sheaths. The electric field is maximum in the sheath and minimum in the bulk region. The electric field in the sheath changes almost linearly from both electrodes to the edge region of the bulk plasma. The movement of electrons causes a modulation of sheath width. The sheaths are quite wide causing the discharge to behave capacitively. The discharge parameters are symmetric about the discharge mid-plane because we assumed that the geometry consists of parallel plates with a diameter that is larger than the electrode separation.. The model is in good qualitative agreement with experimental measurements. 53 LAHOUARIA SETTAOUTI and ABDERRAHMANE SETTAOUTI References BOGAERTS, A., GIJBELS, R. 1995. The Role of Fast Argon Ions and Atoms in the Ionization of Argon in a Direct Current Glow Discharge: A Mathematical Simulation. J. Appl. Phys. 78(11): 6427-6431. CHOI, S. W., LUCOVSKY, G., BACHMANN, K. J. 1992. Remote Plasma Enhanced Chemical Vapor Deposition of GaP with in situ Generation of Phosphine Precursors. J. Vac. Sci. Technol. B10(3): 1070-1073. CHRISTOPHOROU, L. G., 1971. Atomic and Molecular Radiation Physics. Wiley, New York. DELCROIX, J. L., 1960. Introduction to the Theory of Ionized Gases. Wiley , New York. DENPOH, K., NANBU, K. 2000. Self Consistent Particle Simulation of Radio Frequency CF4 Discharge: Effect of Gas Pessure. Jpn. J. Appl. Phys. 39(Part. 1., 5B): 2804-2808. FRANKLIN, R. N. 1999. The Fluid Model of the Positive Column of a Discharge with Negative Ions at Low Pressure Joining Plasma and Sheath. J. Phys. D: appl. Phys. 32: L71-L74. GARRIGUES, L., HERON, A., ADAM, J. C., BOEUF, J. P. 2000. Hybrid and particle in Cell Models of a Stationary Plasma Thruster. Plasma Sources Sci. Technol. 9: 219-226. GOGOLIDES, E. 1997. A Synthetic Approach to RF Plasma Modeling Verified by Experiments: Demonstration of a Predictive and Complete Plasma Simulator. Jpn. J. Appl. Phys. 36(Part. 1., 4B): 2435-2442. GUO, X. M., ZHOU, T. D., PAI, S. T., 1996. Experimental Study of Spatial Distribution of Ar Glow Discharge Plasma. Phys. Plasmas, 3(10): 3853-3857. HELIN, W., ZULI, L., DAMING, L. 1996. Monte Carlo Simulation for Electron Neutral Collision Processes in Normal and Abnormal Discharge Cathode Sheath Region. Vacuum. 47(9): 1065- 1072. HUXLEY, L. G. H. and CROMPTON, R. W., 1974. The Diffusion and drift of Electrons in Gases. Wiley, New York. LOFFHAGEN, D., SIGENEGER, F., WINKLER, R. 2002. Study of the Electron Kinetics in the Anode Region of a Glow Discharge by a Multiterm Approach and Monte Carlo Simulations. J. Phys. D: Appl. Phys. 35: 1768-1776. NANBU, K. 2000. Probability Theory of Electron-Molecule, Ion-Molecule, Molecule-Molecule, and Coulomb Collisions for Particle Modeling of Materials Processing Plasmas and Gases. IEEE Trans. Plasma Science. 28(3): 971-990. SATO, A. H., LIEBERMAN, M. A., 1990. Electron-Beam Probe Measurement of Electric Fields in rf Discharges. J. Appl. Phys., 68(12): 6117-6124. SETTAOUTI, A., MIMOUNI, A., 1999. Monte Carlo Simulation of Electron Swarms in SF6 in Uniform Electric Fields. The Third International Conference on Computational Aspects and Their Applications in Electrical Engineering, ''CATAEE'99''. October 19-20, 1999, Philadelphia University, Jordan. SETTAOUTI, A., AOUABDI, L., 2001. Monte Carlo Simulation of Electron Swarms in Nitrogen in Uniform Electric Field. 4th Jordanian International Electric & Electronics Engineering Conference (JIEEEC'2001), April 16-18, 2001, Amman, Jordan. VOSSEN, J. L., KERN, W.,1978. Thin Film Processes. Academic Press. INC. San Diego (USA). WANG, D. Z., DONG, J. Q. 1997. Kinetics of low pressure rf discharges with dust particles. J. Appl. Phys. 81(9): 38-42. YONEMURA, S., NANBU, K. 2001. Electron Energy Distribution in Inductively Coupled Plasma of Argon. Jpn. J. Appl. Phys. 40(Part. 1, 12): 7052-7060. Received 15 July 2002 Accepted 17 April 2003 54 A Monte Carlo Method for Low Pressure Radio Frequency Discharges Lahouaria Settaouti and Abderrahmane Settaouti Department of Electronics, University of Sciences References