A Monte Carlo Method for Low Pressure Radio Frequency Discharges

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.


Introduction
G low 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 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 semiquantitative 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.

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 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 T m is the mean flight time of an electron moving with a velocity v(ε), in which v(ε) is the drift velocity of an electron.where Q t is the total collision cross-section in m 2 , N the gas number density and ε the electron energy in eV.Q t 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 R 1 at the end of each step, where R 1 is a random number uniformly distributed between 0 and 1.When the electron undergoes a collision The nature of the collision is determined in the following way: P 2,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 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 R 2 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: , where E exc is the excitation threshold and E 0 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 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 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.

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.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.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 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.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, 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.

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.

Figure 1 .
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.

Figure 3 .Figure 4 . 2 .
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).

Figure 4 . 4 .
Figure 4.4.Temporal and spatial variations Figure 5.1.electron density as a function of the electric field for the fourth rf cycle.ofposition for the fourth rf cycle.

Figure 5 . 2 .
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.