Modeling of Low-temperature Plasma

Plasma is ionized gas. Hence, it consists of ions and electrons, as well as neutral species. The ionization degree of plasma varies from 100% (fully ionized plasma) to

very low values (e.g. 10⁻⁴–10⁻⁶; weakly ionized plasma). Besides the space plasma, the laboratory plasma is divided into two main groups, which are the high-temperature plasma (fusion plasma), and the low-temperature plasma (gas discharges plasma). Moreover, two sub-groups of gas discharge plasma are classified depended on its working gas pressure, which are thermal equilibrium plasma and non-thermal equilibrium plasma [Lieberman and Lichtenberg, 1994]. The efficient energy exchange between the plasma species due to many collisions occur for high-pressure plasma. Thermal equilibrium discharge is typically used for applications where heat is required, such as for cutting, spraying, welding. On the other hand, for low gas pressure plasma, different temperatures of the plasma species due to its inefficient energy transfer. Non-thermal equilibrium plasma is typically used for applications where heat is not desirable. In recent years, this field of non-thermal equilibrium plasma applications has rapidly expanded due to its non-equilibrium aspect of the plasma. The latter sub-group of gas discharge plasmas is also the second subject of this dissertation.

Some important operating parameters for obtaining different non-equilibrium conditions are briefly summarized as follows [Economou, 2000]:

z The chemical input of working gas and its corresponding gas pressure z The imposed external electromagnetic field structure

z The configurations of plasma chamber and electrodes z The temporal behavior (e.g. pulsing the plasma)

One can realizes that the non-equilibrium plasma conditions are strongly influenced by any one or many of these summarized parameters. For example, the representation of the parameter space in plasma etching is shown in Fig. 1.1. It clearly illustrates that the above-mentioned operating parameters (in the top block) determine the following key plasma properties (listed in the middle block) including the electron velocity distribution function (EVDF), the space and time variation of electron, etc.

Finally, these properties may dominate the figures of merit (listed in the bottom block) including the etching (or deposition) rate, uniformity, etc. Therefore, a computer simulation code should use such many parameters as inputs to help optimize the expected non-equilibrium plasma conditions easily and understand the underlying physics.

Dimensionality of plasma reactor simulations ranges from zero-dimensional to three-dimensional. Low dimensional simulations, such as zero-dimensional, e.g., [Font et. al.,1998], [Meeks and Shon, 1995], and [Deshmukh. and Economou, 1992], and one-dimensional models, e.g., [Midha and Economou, 2000], [Kline et. al. 1989], and [Nedelea and Urbassek, 2004], are best choice in handling very complicated chemistry [Meeks et. al., 1997]. However, two-dimensional simulations can address

the important aspect of reaction uniformity across the wafer radius, e.g., [Shon and Lee, 2002], [Shon et. al.,1999], and [Shon et. al.,1998]. Three-dimensional

simulations are useful for studying azimuthal asymmetries in the reactor due to non-axisymmetric power deposition, or non-axisymmetric gas inlets and pumping ports, e.g., [Kushner et. al.,1996], and [Kushner, 1997].

In general, there are three kinds of plasma simulation approaches; the first is the fluid model, the second is the kinetic model and the third is the hybrid model. In addition, Maxwell’s equations for electromagnetic fields also need to be solved self-consistently coupled with the plasma densities and currents from plasma simulations. For the fluid model, one need to solve the equations, which are derived after taking the moments of the Boltzmann equation with some assumptions regarding, e.g., [Meyyappan, 1994], and [Gogolides and Sawin, 1992]. They are species continuity equation, species momentum equation and species energy equation.

Related publications of fluid model could be found in numerous articles, e.g., [Ventzek et. al.,1993], [Lymberopoulos et. al.,1995], [Bukowski and Graves, 1996], and [Ventzek et. al.,1994], and are not repeated here. Unlike fluid model, kinetic approach yields the particle distribution functions as an output of the simulation.

Especially, it is more accurate than fluid model at low pressures when the species mean free path is comparable to or longer than a characteristic length scale or for

highly non-equilibrium situations. However, Kinetic approach is computationally intensive as compared to fluid model. One of the well-established kinetic approaches is the Particle-In-Cell with Monte-Carlo Collisions (PIC-MCC) method, e.g., [Birdsall, 1991], and [Vahedi et. al.,1993]. In the past two decades, PIC-MCC method has long

been used to study the nonlinear kinetic problems in space and laboratory plasma physics. For self-consistent treatment of the plasma and the background gas, Nanbu combined the Direct Simulation Monte Carlo method (DSMC) with PIC-MCC, e.g., [Nanbu, 2000], and [Serikov et. al.,1997].

Each time step in the PIC-MCC consists of four major steps: charge extrapolation, force interpolation from the solution of the Maxwell’s equations, particle movement, and Monte-Carlo collisions. Briefly speaking, based on the particle positions, charges are assigned to each mesh point and current densities are assigned to the faces between the mesh points. Maxwell's equations are then solved to compute the electric and magnetic fields on the grid. The force on the particles is obtained from the fields at these gird points by interpolation based on the particle position. Particles are then moved according to Newton's law. Particle collisions are handled stochastically in a Monte Carlo module in-between field adjusting time steps. The details of the PIC-MCC method will be given in chapter 4.

As mentioned earlier, most conventional PIC-MCC, e.g., [Birdsall and Langdon,

1991], and [Birdsall, 1991], a structured mesh is usually construed for the

computational domain. However, until very recently, there have been few developments of electrostatic PIC method using a three-dimensional unstructured mesh, mostly designed for thruster plume simulations due to their complicated computational geometry. A hybrid PIC-DSMC code using unstructured mesh, called AQUILA, which has been developed by [Santi et al., 2003] on hall thruster plume simulation. They obtained the improved current density results from better the unstructured mesh resolution. AQUILA uses finite element method to discretize Poisson’s equation with electrons from Boltzmann relation, and then uses Newton-Raphson method to solve the non-liner resulting matrix. [Spirkin et al., 2004]

has also developed a three-dimensional Particle-In-Cell code on a unstructured tetrahedral mesh with finite volume method. This PIC code was applied to the simulation of the flow inside the segmented micro-channel of a directional-retarding potential analyzer. Results show the flow characteristics of the ions and electrons and provide estimates of the collected current by the micro-plate.

Parallel PIC-MCC method has been previously studies by various researchers using different schemes, e.g., [Seidel et al.,2002], [Kawamura et al.,2000], [Decyk, 2002], [Walker,1991], [Walker, 1990], [Lee and Azari, 1992], [Akarsu et al.,1996], [Decyk and Norton, 2004], and [Liewer and Decyk, 1989], since it is the most

computationally demanding compared with other models. Most parallel PIC-MCC schemes, e.g., [Kawamura et al.,2000], [Seidel et al.,2002], and [Lee and Azari, 1992]

employ a Eulerian decomposition scheme in which just paralleling the particle processing without paralleling the field solver since the field solver can be a small percentage of the work load especially when FFT methods are used. In this report, for a fixed number of grid points, the speedup just for this parallel particle processing became more linear with increasing particle number on 2 and 4 CPU symmetric multiprocessor (SMP) machines and on a distributed network of workstations (NOW).

In the past, there have been very few studies concerning on developing dynamic load-balancing technique for particle-based PIC-MCC code, e.g., [Seidel et al.,2002] , [Decyk and Norton, 2004], and [Liewer and Decyk, 1989]. In Seidel’s work [Seidel et al.,2002], he has proposed a method in which the code will dynamically repartition

the computational domain and intends to balance the workload among processors under the framework of structured mesh. Decyk et al. have developed a new algorithm just for PIC method on concurrent processors with distributed memory, which named the general concurrent PIC algorithm (GCPIC). In this algorithm, the physical domain of particle simulation was divided into sub-domains, which are equal to the number of processors. The sub-domain can be re-created to keep the processor loads of particle computations balance (dynamic load balancing) during the transient

period of the simulation, which was called primary decomposition. Again, each sub-domain may have equal numbers of particles, however, unequal numbers of grid numbers. Thus, GCPIC used secondary decomposition to divide the physical domain into number of processors equal sub-domains with equally number of grid points under MIMD paradigm for computing field solver efficiently. However, these reviews showed that parallel PIC methods are not suitable using the SIMD paradigm.