0018-8646/99/$5.00 (C) 1999 IBM Modeling and simulation methods for plasma processing by S. Hamaguchi Methods used for the modeling and numerical simulation of the plasma processes used in semiconductor integrated-circuit fabrication are reviewed. In the first part of the paper, we review continuum and kinetic methods. A model based on the drift-diffusion equations is presented as an example of a continuum model; the model and associated numerical solutions are discussed. The most widely used simulation method for kinetic modeling is the Particle-In-Cell/Monte-Carlo-Collision (PIC/MCC) method, in which the plasma is modeled by a system of charged superparticles (each of which represents a collection of a large number of ions or electrons) that move in self-consistent electromagnetic fields and collide via given collision cross sections. In the second part of the paper, we review the modeling and simulation of the evolution of surface topography in plasma etching and deposition. 1. Introduction Modeling and numerical simulations of plasma processing can be useful in many ways. An improved understanding of a plasma-processing system can be achieved by comparing predictions from numerical simulations with experimental observations. The optimization of existing processes or the development of new plasma processes that offer better processing results may follow such an understanding. Modeling and simulations based on reliable physical/chemical modeling of a plasma-processing system can significantly reduce the number of associated experiments that otherwise would have to be performed. Additionally, such modeling and simulations can also be used in the computer-aided design (CAD) of new process systems and to optimize manufacturing processes within the framework of existing processing systems. The ultimate goal may be the modeling and simulation of an entire manufacturing system. Although current methods are far from achieving this goal, several physics-based simulation tools have already played important roles in the development of several plasma-processing systems. In this paper we focus on the modeling and simulation of plasma processes used in the fabrication of integrated circuits. Conventional plasma-processing tools are essentially based on parallel-plate discharges [1], in which a voltage is applied to electrodes to break down a gas and sustain the generated discharges. Plasmas generated in this method are usually only partially ionized and have high neutral gas pressures. Because higher ion fluxes toward substrates are generally more desirable in plasma processing, high-plasma-density processing tools [2] with lower gas pressures are beginning to replace conventional plasma tools. Several different methods can be used to generate high-density plasmas for semiconductor applications. For example, oscillating rf electric fields can be used to generate and sustain a high-density plasma without using steady-state magnetic fields [known as an inductively coupled plasma (ICP)] or by using steady-state magnetic fields (known as a helicon plasma because helicon waves are induced along the magnetic field lines). Additionally, electron cyclotron resonance (ECR) can be used to generate a high-density plasma. A neutral-loop discharge (NLD) is a recently developed method in which a plasma is generated by rf fields along a closed magnetic neutral line [3]. The goal of this paper is to present fundamental aspects of plasma-process modeling based on relevant physical/chemical mechanisms. To achieve this goal and because of the limited available space, we are unable to discuss details of the most recent development of plasma-process simulation tools. However, references to such work are cited, and the interested reader is encouraged to examine them. In Section 2, a continuum model is discussed which is more appropriate for highly collisional, low-density/high-pressure discharges, although such a model is also known to be applicable to some high-density/low-pressure plasmas. For less collisional plasmas with possibly non-Maxwellian ion and/or electron distribution functions, we present in Section 3 a particle method as a numerical technique for solving kinetic plasma models. The sheath region of the plasma has different characteristics from the bulk plasma (such as the presence of strong electric fields) and requires different modeling methods, which are discussed in Section 4. Once fluxes from the plasma source are determined from model calculations, one can predict the outcome of materials processing by modeling processes [4] that occur on the substrate, such as etching and deposition; this is discussed in detail in Section 5. 2. Continuum models If the neutral-gas pressure is sufficiently high, collisional processes dominate the dynamics of the plasma. To model such a plasma, drift-diffusion equations are widely used [5-10]. The drift-diffusion equations are "reduced" equations derived from the well-known fluid equations (i.e., mass- and momentum-conservation laws) for ions and electrons under the assumption of high collision frequencies. For example, the momentum-conservation law for electrons may be written as m[SUB]e n[SUB]e ([round dee]v[SUB]e / [round dee]t + v[SUB]e o (1) [nabla]v[SUB]e) = -[nabla]p[SUB]e -en[SUB]e E-m[SUB]e n[SUB]e [nu][SUB]e v[SUB]e , where m[SUB]e is the electron mass, n[SUB]e is the electron density, v[SUB]e is the electron "fluid" (drift) velocity, p[SUB]e is the electron fluid pressure, E is the electric field, and [nu][SUB]e is the electron collision frequency at which the plasma electrons collide with the neutral gas molecules. Electron-ion collisions are negligible compared with electron-neutral collisions in a high-pressure plasma. If the inertia term (i.e., the left-hand side) of Equation (1) is negligibly small, we obtain v[SUB]e = -[nabla]p[SUB]e / m[SUB]e n[SUB]e [nu][SUB]e - eE (2) / m[SUB]e [nu][SUB]e. Equation (2) states that v[SUB]e is determined by the balance among the pressure force, electrostatic force, and drag force. The electron flux [Gamma][SUB]e = n[SUB]e v[SUB]e is then given by [Gamma][SUB]e = -n[SUB]e [muon][SUB]e E-[nabla](D[SUB]e n[SUB]e), (3) where [muon][SUB]e = e/m[SUB]e [nu][SUB]e is the electron mobility and D[SUB]e = k[SUB]B T[SUB]e / m[SUB]e [nu][SUB]e is the electron diffusivity. The equation of state p[SUB]e = n[SUB]ek[SUB]BT[SUB]e was also used. The relation D[SUB]e/[muon][SUB]e = k[SUB]BT[SUB]e/e is known as the Einstein relation. Similarly, the ion drift velocity v[SUB]i and ion flux [Gamma][SUB]i = n[SUB]iv[SUB]i can be expressed as [Gamma][SUB]i=n[SUB]i [muon][SUB]i E - [nabla](D[SUB]i n[SUB]i), (4) where [muon][SUB]i = Z[SUB]i e/m[SUB]i [nu]i is the ion mobility, Z[SUB]i e and [nu]i are the ion charge and ion-neutral collision frequencies, and D[SUB]i = k[SUB]B T[SUB]i / m[SUB]i [nu]i is the ion diffusivity. If the values of the mobilities and diffusivities are given, we need equations only for the electron and ion densities, n[SUB]e and n[SUB]i, and the electric field E to close the system. This can be achieved by using the mass conservation laws and Poisson's equation: [round dee]n[SUB]e / [round dee]t+[nabla] o [Gamma][SUB]e = S[SUB]e , (5) [round dee]n[SUB]i / [round dee]t+[nabla] o [Gamma][SUB]i = S[SUB]i , (6) [Delta][phi]= - e / [epsilon][SUB]0 (Z[SUB]i n[SUB]i - n[SUB]e). (7) In Equations (5) and (6), S[SUB]e and S[SUB]i denote the source terms for the electrons and ions, to be discussed below. In Poisson's equation (7), [phi] denotes the electrostatic potential, i.e., E = -[nabla][phi]. Equations (5)-(7), together with Equations (3) and (4), constitute the drift-diffusion equations for the system. Although we have considered only single-ion species, it is straightforward to extend Equations (5)-(7) to those for the systems of multi-ion species, including negative ions. Sources and sinks of electrons and ions in the bulk plasma are usually due to ionization, recombination, electron attachment, and charge exchange. In the case of electropositive discharges with single-ion species, we may write S[SUB]e = S[SUB]i = k[SUB]i n[SUB]e n[SUB]n - k[SUB]r n[SUB]i n[SUB]e , where k[SUB]i and k[SUB]r are ionization and recombination rate coefficients. (In general, S[SUB]e [not equal to] S[SUB]i, especially for electronegative plasmas. See for example References [10-12].) As is well known, there is a threshold energy E[SUB]th for ionization, and the ionization rate coefficient depends on the electron temperature T[SUB]e. It is often assumed that the electron temperature dependence of the ionization rate coefficient is given by the Arrhenius relation, i.e., k[SUB]i=A[SUB]i exp (-(E[SUB]th / k[SUB]B T[SUB]e)) , (8) where A[SUB]i is a constant. For example, for argon, which is commonly used in plasma processing, A[SUB]i [approximately equal to] 1 x 10**-14 m**3/s. Ions and electrons are also lost to the walls and electrodes of the chamber that contains the plasma. Some electrons are also created at the walls and electrodes due to secondary-electron emission. These effects must be incorporated in the boundary conditions. The other boundary conditions may be given as follows. For typical parallel-plate rf discharges, the potential [phi] at the anode and wall boundaries may be set to zero, whereas the potential at the cathode may be set at the applied voltage, e.g., [phi] = V[SUB]dc + V[SUB]rf cos [omega]t, where [omega] is the applied rf angular frequency and V[SUB]dc and V[SUB]rf are the dc and ac components of the applied voltage. Depending on the physical conditions one wishes to consider, there are a variety of choices for boundary conditions for electron and ion fluxes. For example, we may set [Gamma][SUB]e=k[SUB]r[SUP](s) n[SUB]e - [gamma][Gamma][SUB]i , [Gamma][SUB]i=n[SUB]i [muon][SUB]i E, at the boundaries, where [Gamma][SUB]e, [Gamma][SUB]i, and E are the surface normal components of [Gamma][SUB]e, [Gamma][SUB]i, and E, respectively; k[SUB]r[SUP](s) is the electron surface recombination coefficient; and [gamma] is the secondary-electron emission coefficient. Note that [Gamma][SUB]e ([Gamma][SUB]i) is defined to be positive if the electrons (ions) flow into the boundary, and [gamma] is set to zero if [Gamma][SUB]i < 0. Well before the widespread use of the drift-diffusion equations for plasma-process simulations, equations similar to Equations (5) and (6) were extensively used to model transport of electrons and holes in semiconductor materials [13-15]. Such modeling of semiconductor devices by similar drift-diffusion equations (which are also known as the fundamental semiconductor equations) dates back to the work by Gummel [16] in 1964. Thanks to such extensive early studies, we now have access to accumulated knowledge of mathematical structures and numerical techniques for the drift-diffusion equations. As an example of relevant simulation results, Figure 1 depicts one-dimensional results for parallel-plate helium discharges obtained by Boeuf [17]. The simulation is based on a system of equations similar to Equations (5) and (6). In the figure, spatial profiles of the electric field (E) and ion and electron densities are shown at four different times in an rf cycle. The assumptions made in this simulation are the following: The conducting electrodes are separated by a distance of 4 cm, the gas pressure and temperature are 1 Torr and 300 K, respectively, the driving frequency [omega]/2[pi] is 10 MHz, the driving voltage V[SUB]rf is 500 V with no dc component (V[SUB]dc = 0), and the secondary-electron emission coefficient [gamma] is 0.1. So far we have assumed that the electron temperature profile can be characterized by a given function (or constant). However, as is evident from Equation (8), the electron temperature profile can affect the ionization rate, which in turn affects the plasma profile. To treat the electron temperature more realistically, one can couple the electron-energy transport equation with Equations (5) and (6) through Equation (8). The equation for the electron-energy transport can be expressed [18] as 3/2 [round dee]/[round dee]t (n[SUB]e k[SUB]B T[SUB]e) = (9) -[nabla] o q[SUB]e + Q[SUB]e, where the electron-energy flux q[SUB]e and the rate of electron heat generation Q[SUB]e can be approximated by q[SUB]e= 5/2 [Gamma][SUB]e k[SUB]B T[SUB]e - [Kappa][SUB]e [nabla](k[SUB]B T[SUB]e), with [Kappa][SUB]e being the electron heat conductivity and Q[SUB]e = -e[Gamma][SUB]e o E-n[SUB]e n[SUB]n (10) [SUM] k[SUB]j [epsilon][SUB]j + Q[SUB]ext. j In Equation (10), the first term represents the joule heating, and the second term represents the electron-energy loss due to collisions (such as ionization, excitation, and elastic collisions). Here k[SUB]j and [epsilon][SUB]j denote the collision rate coefficient and energy loss per electron per collision for the collision process denoted by j [e.g., for ionization, k[SUB]j = k[SUB]i, as given in Equation (8)]. The last term of Equation (10) represents an electron heat source associated with external effects, for example, the power-deposition profile in an ICP or ECR plasma. Continuum equations more general than Equations (5) and (6) coupled with Equation (9) have also been used to model various discharges. For example, one may employ the full ion momentum equations instead of the drift-diffusion approximation given by Equation (4). Some earlier one-dimensional simulations for rf glow discharges based on continuum models may be found in References [5-10, 17, 19]. More recently, the continuum equations have been extended to higher dimensions. Two-dimensional simulation results for rf glow discharges are found, for example, in References [11, 12, 20, 21]. Although it is known that the overall structures of plasma discharges are well described by continuum equations such as Equations (5) and (6), these equations do not model the change in the velocity-distribution functions, inasmuch as all of the distribution functions are assumed to be Maxwellian. As Godyak et al. [22-24] have demonstrated experimentally, however, some rf discharges exhibit bi-Maxwellian electron distributions, i.e., distributions with two electron temperatures. Since electron kinetic energies, especially those near ionization threshold energies, affect ionization rates significantly, it is necessary to consider the time evolution of distribution functions if one wishes to model the plasma more accurately. To this end we must use kinetic and hybrid models, as discussed in the next section. 3. Kinetic and hybrid models Kinetic modeling of plasmas determines time evolution of the velocity-distribution functions. For a rarefied gas or plasma, where interparticle correlations are sufficiently weak, the time evolution of distribution functions is known to be governed by the Boltzmann equation [18, 25]. For example, the electron-distribution function f[SUB]e (t, x, v) can be characterized by the Boltzmann equation [round dee]f[SUB]e / [round dee]t + v o [round dee]f[SUB]e / (11) [round dee]x - eE / m[SUB]e o [round dee]f[SUB]e / [round dee]v = ([delta]f[SUB]e / [delta]t)[SUB]coll , where [round dee]/[round dee]x (= [nabla]) and [round dee]/[round dee]v are gradients in the real and velocity spaces, and the right-hand side represents the collision term. The specific form for the collision term ([delta]f[SUB]e / [delta]t)[SUB]coll depends on the collision processes to be considered, but it typically includes ionization, excitation, dissociation, and recombination. For ions, a distribution function may be defined for each excited state of each species, which is also governed by an equation similar to (11). Obtaining numerical solutions to the Boltzmann equation under general conditions is a difficult computational problem. One scheme that can be used to solve the equation exactly is the convective scheme [26-28]. A more widely used alternative method is the particle in cell (PIC) method [29, 30]. In the PIC method, a group of electrons (or ions) are represented by a single simulation particle (superparticle or pseudoparticle), and the motion of each particle is assumed to follow Newton's law. Each superparticle can represent 10**8-10**12 real particles (such as electrons), so that 10000 superparticles--a manageable number on today's engineering workstations--may represent all of the electrons and ions in (a portion of) an actual plasma system. It is also assumed that such superparticles have no interparticle interaction except for very short-range interactions corresponding to collisions: Forces exerted on superparticles originate from self-consistent electromagnetic fields. In this way, only collective motions (except for collisions) are modeled by a PIC simulation. For parallel-plate rf discharges, only the electrostatic field need be taken into account. The equations of motion for the superparticles are therefore given by [xdots][SUB]j = -q[SUB]j / m[SUB]j [nabla][phi], (12) [Delta][phi]=-[rho] / [epsilon][SUB]0 . (13) Here x[SUB]j = x[SUB]j(t) is the position of the jth (electron or ion) superparticle at time t, q[SUB]j and m[SUB]j are its charge and mass, and [rho] is the charge density. In the PIC method, the field quantities such as potential [phi](x, t) and charge density [rho](x, t) are numerically evaluated only at spatial grid points, whereas particles are allowed to occupy any position within the discharge. The right-hand-side terms of Equation (12) are therefore evaluated from the interpolation of [phi] values at neighboring grid points (force weighting), whereas the value of [rho] in the right-hand side of Equation (13) at each grid point is evaluated from charges of neighboring particles (particle weighting). For these weighting processes to be accurate, a sufficient number of simulation particles must be present in each grid cell--thus the name "particle in cell." The boundary conditions for Equations (12) and (13) may be developed as follows. When an electron reaches the boundary, it is assumed to be adsorbed. For an ion, it is also assumed to be adsorbed, but secondary electrons may be emitted with a probability [gamma], depending on the impinging ion energy. One may also assume that some charged particles are reflected by the boundary or subjected to some chemical reactions at the boundary. Choices of boundary conditions depend on the physical conditions of the boundary walls and electrodes. There are many other issues associated with the PIC algorithms which are beyond the scope of this paper, and interested readers are referred to, e.g., References [29, 30]. To simulate highly collisional plasmas using the PIC method, the Monte Carlo collision (MCC) method has been widely adopted [31-36]. By using the knowledge of collision cross sections for all possible collisions and generating a sequence of random numbers, the MCC method can be used to determine whether a specific collision takes place during a given time interval; if it occurs, the particle velocity is revised accordingly. There is an efficient scheme to perform this task, designated as the null collision algorithm; detailed descriptions of the scheme can be found, for example, in References [37-40]. The MCC process is typically performed after the possible new positions of all of the particles are determined in each simulation cycle. Figure 2 schematically shows a sequence of computational processes in a single time-step cycle for a PIC/MCC simulation. Vahedi et al. [41] applied PIC/MCC simulation to simulate argon rf glow discharges and have examined the applicability of bi-Maxwellian electron distribution functions, which are typically observed under conditions of relatively low gas pressure. Prior to their simulations, Godyak et al. measured electron-energy distribution functions [22-24] in rf discharges and postulated that the bi-Maxwellian distributions result from two competing electron-heating mechanisms taking place in the discharge. One heating mechanism is ohmic heating (i.e., collisional heating), which typically takes place in the bulk plasma, and the other is stochastic heating [42, 43], which occurs near the plasma-sheath boundaries. In stochastic heating, electrons gain energy from oscillating rf sheaths. The concept proposed by Godyak et al. is the following: Low-energy electrons are generated by inelastic collisions, such as electron-impact ionizations, in the bulk plasma. Especially if the electron energies are near the Ramsauer minimum, such electrons cannot be thermalized because of the low collision frequency. Therefore, most of them simply oscillate collisionlessly and remain in an abnormally low temperature range. Some of the low-energy electrons, however, reach the bulk-sheath boundary region and can gain energy from the stochastic heating mechanism. High-energy electrons generated in this way are eventually lost to the wall or lose their energy through inelastic collisions, and a steady state with two electron temperatures is maintained. Figure 3 shows electron-energy distribution functions measured by Godyak et al. [22-24]. Figure 4 shows the electron-energy distribution function obtained from one-dimensional PIC/MCC simulations (one dimension in coordinate space and two dimensions in velocity space) by Vahedi et al. [41] under approximately the same conditions as those by Godyak et al.: parallel-plate argon rf (13.56 MHz) discharge with a gap separation of 2 cm, discharge current 2.56 mA-cm**-2, gas pressure 100 mTorr. It is seen that the simulated electron-energy distribution function, where the broken lines represent two temperatures (0.5 eV and 3.0 eV) of the bi-Maxwellian distribution [41], is in good agreement with the experimental results of Figure 3. Another observation made in both experiments and simulations is that the electron-energy distribution varies from a bi-Maxwellian to an ordinary Maxwellian function as the neutral gas pressure is increased [24, 41]. This transition is caused by a transition in the electron-heating mode from a stochastic-heating-dominant regime at low pressures to an ohmic-heating-dominant regime at high pressures. The main advantage of a kinetic simulation over a continuum-model simulation is its usefulness for obtaining detailed information on discharge characteristics, such as velocity-distribution functions. Although we have thus far discussed only electron-distribution functions, distribution functions of ions and radicals can also be simulated using the PIC/MCC method. As we see in the next section, energy and angular distributions of ions and radicals near the wafer surface strongly influence microscopic material-surface processes such as etching and deposition. Therefore, kinetic modeling is always desirable for such process simulations. The major disadvantage of a kinetic simulation is, however, its computational cost. Extending the PIC/MCC method to a system of higher dimensions (e.g., 2D in coordinate space and 3D in velocity space) may be relatively straightforward in terms of numerical methods, but running such simulations with satisfactory accuracy can be extremely expensive. However, some progress has been made in this direction. For recent development in 2D PIC/MCC simulations for materials process applications, see References [44, 45]. To reduce the computational cost associated with kinetic simulations but still retain some of their advantages, several researchers have used hybrid schemes, i.e., combinations of continuum simulations and kinetic simulations [40, 46-53]. For example, Sommerer and Kushner [40] used the "test-particle MC method" to simulate rf discharges, in which the self-consistent fluid equations (for both ions and electrons) were supplemented by particle models of electrons subjected to MCC processes. In this method, the electron-simulation particles are used only as "test particles" to evaluate the electron-distribution function under given ion, electron, and neutral-species densities and electric fields obtained from the fluid equations. The obtained electron-distribution functions are then used to determine various reaction rates (including ionization and excitation rates) and transport coefficients as functions of space and time. Since the electron density n[SUB]e (t, x) is obtained by solving the electron fluid equations, not by counting the number of electron-simulation particles in a unit volume, this is not a PIC simulation. The lack of the particle-weighting process of Figure 2 allows one to use a very small number of electron-simulation particles, typically 1% or less of the equivalent PIC simulation. The advantage of this kinetic-fluid hybrid method is that detailed information on the electron-distribution function is obtained, making it possible to more accurately estimate various reaction rates without adding too much computational cost. Depending on the physical/chemical phenomena to be studied, different types of kinetic-fluid hybrid methods may be adopted. For example, Porteous and Graves [51] employed an electron-fluid ion-particle method to study magnetically confined plasmas such as ECR plasmas with long ion mean free paths. In contrast to the test-particle MC method discussed above, Porteous and Graves used fluid equations only for electrons, and modeled ions completely by particles with MCC. Their particle MCC method was, however, not the PIC/MCC method either, because the ion density was not calculated from the number density of the simulation particles at each instance of time but was averaged over trajectories of many simulation particles. This simulation scheme is essentially the direct-simulation Monte Carlo (DSMC) method, which has been used extensively to model neutral rarefied gas flows [54]. Owing to this trajectory-averaging strategy, fewer particles are needed to obtain sufficient statistics in the DSMC method than in the PIC/MCC method. However, in contrast to PIC/MCC simulations, only time-averaged steady-state solutions can be obtained from DSMC simulations. More recently, the DSMC method coupled with an electron-fluid model has been used to model the dynamics of both ions and neutral species in two-dimensional high-density ICP plasmas [55, 56]. It is relatively easy to implement various collision processes in particle simulations such as the PIC/MCC and DSMC methods, compared to those based on continuum models such as those characterized by the drift-diffusion equations. The major problem of particle simulations is, as mentioned above, the computational cost. The problem is especially severe for two- or three-dimensional systems, and it is often necessary to use high-performance vector or parallel computers to obtain meaningful results. To accurately simulate high-density low-pressure plasmas, where ion and neutral mean free paths are comparable with the chamber dimensions, it may be necessary to use particle simulations such as those discussed above. However, for most practical purposes (e.g., designing plasma-processing tools), it is known that continuum (or fluid) simulations are more feasible because of their reasonably predictive capabilities and lower computational costs, even for high-density-plasma simulations. The real challenge in simulating high-density plasmas with low gas pressures often lies in evaluating the power depositions self-consistently. As mentioned earlier, most high-density plasmas for materials processing applications are generated by means of high-frequency driving sources such as rf coils (for ICP and NLD), rf antennas (helicon-wave and surface-wave plasmas), and microwave waveguides (for ECR plasmas). Therefore, to obtain the power deposition profile self-consistently, one must solve the Maxwell equations together with the equations for the plasma. While the power-deposition profiles [i.e., Q[SUB]ext (x) in Equation (10)] were simply assumed to be represented by given functions in some earlier studies [51, 53], recent simulation studies attempt to solve the coupled Maxwell equations with some simplifying assumptions (such as a cold-plasma approximation for the dielectric tensor for Maxwell's equations) to evaluate the power-deposition profiles [57]. For recent two- and three-dimensional simulations of self-consistent high-density plasma sources for materials-processing applications, the reader is referred, for example, to References [55-66]. 4. Sheath modeling The sheath region in a plasma-processing system usually denotes a gaseous boundary region between the bulk plasma and the electrode, where the electron density is low and strong electric fields are present. The materials processes are strongly influenced by the characteristics of ion- and neutral-species fluxes impinging on the substrate through the sheath. For dry etching and deposition processes, energy and angular distributions of those incoming fluxes are two important factors. As discussed in the previous section, obtaining information on ion- and neutral-species distribution functions from kinetic or hybrid simulations can be very costly. On the other hand, less computationally intensive fluid-type simulations provide no information on the distribution functions. To compromise, one may use less expensive fluid-type simulation to model the global structures (such as ion and electron temperatures and densities) of the plasma and treat the sheath region separately, using a kinetic model. Much effort has been made to measure and compute the ion-distribution functions in sheaths of various discharges. Davis and Vanderslice [67] made the first systematic measurements of the energy distribution of ions striking the cathode and found that the energy distributions vary considerably according to the collisionality (i.e., the ratio of the sheath thickness to the mean free path). Taking into account only charge-exchange collisions, Davis and Vanderslice also proposed a simple model of the ion distribution and succeeded in explaining some of the measured energy distributions. Kushner [68] computed the energy and angular distributions of the impinging ions in rf discharges, using the Monte Carlo (MC) technique for charge-exchange collisions. Thompson et al. presented more comprehensive MC simulations [69] for several different forms of elastic ion-neutral interactions as well as charge-exchange collisions in dc and rf discharges with predefined electric fields. Farouki et al. [70] reported MC simulations with self-consistent electric fields due to space-charge effects, taking into account both elastic scattering and charge-exchange collisions. Under some conditions, analytical expressions for the ion-distribution function can be obtained. Wannier [71] has solved the Boltzmann equation analytically for ions with ion-neutral elastic collisions in the limit of a strong uniform electric field. Lawler [72] has also solved the Boltzmann equation with charge-exchange collisions and obtained an ion-distribution function whose asymptotic limit is Wannier's equilibrium solution. Hamaguchi et al. [73] have obtained analytical expressions for ion-distribution functions in collisional dc sheaths. According to that work, in the limit of weak ion-neutral elastic collisionality with a constant electric field, the normalized ion-energy flux function [Gamma][SUB]en ([eta]) is given by [Gamma][SUB]en ([eta])=-ln(1-[eta]), (14) where [eta] is the normalized kinetic energy of the ion, i.e., [eta] = v[SUB]z**2 / v[SUB]max**2, with v[SUB]max = [square root]v[SUB]B**2 - 2qV/m[SUB]i being the ion velocity resulting from the sheath potential V (<0) and ion charge q in the absence of collisions. The Bohm velocity (or ion sound velocity), v[SUB]B = (k[SUB]B T[SUB]e / m[SUB]i)**1/2, is the initial velocity of the ion when it enters the sheath region from the bulk plasma. Since v[SUB]B is typically small (as k[SUB]B T[SUB]e << -qV), we may write [eta] = m[SUB]i v**2 / (-qV). It is shown in Figure 5 that Equation (14) is an excellent approximation to the self-consistent energy flux function. The exact analytical formula of the ion-energy distribution function in the self-consistent electric field (represented by the solid curve in Figure 5) is given by Reference [73], Equation (104). Figure 6 shows time-averaged ion-distribution functions at the cathode of collisionless rf sheaths obtained by Hamaguchi et al. [74, 75]; the solid lines (upper figures) are from the analytical expressions given below [Equation (15)] and histograms (lower figures) are simulation results [76]. The electric field in the sheath was assumed to oscillate sinusoidally; i.e., E[SUB]0 (x) + E[SUB]1 (x)cos[omega]t, with E[SUB]0 (x) and E[SUB]1 (x) having arbitrary spatial dependences. (The origin was assumed to be located at the plasma-sheath boundary and the electrode was assumed to be located at x = d.) If the ion-transit frequency [omega][SUB]tr = [qE[SUB]0 (d) / m[SUB]i d]**1/2 (i.e., the inverse of the typical time that an ion takes to cross the sheath of width d) is sufficiently small compared to the rf frequency (i.e., [epsilon] = [omega][SUB]tr / [omega] <<[ 1), the time-averaged ion energy flux is given [74] by [Gamma][SUB]en ([eta]) = n[SUB]Iv[SUB]B / (15) 2[pi][square root]([nu][SUB]+ - [square root][eta]) ([square root][eta]-[nu]_) , where n[SUB]I is the ion density at the plasma-sheath boundary (i.e., x = 0) and [nu][SUB]+- = 1 +- qE[SUB]1(d)/m[SUB]i[omega]v[SUB]max**2. The width of the energy spread, in terms of the dimensional quantity[fancyE] = 1/2 m[SUB]iv[SUB]max**2[eta], is given by [Delta][fancyE] = 2qE[SUB]1 (d) / [omega] (-2qV/m[SUB]i)**1/2. It is also known that the ponderomotive force due to the oscillating sheath field exerts a retarding influence on the ion trajectories [74]. 5. Surface process modeling Modeling and numerical simulation of plasma etching and deposition systems can facilitate a better understanding of microscopic phenomena taking place on the substrate. In this section, we discuss the modeling and simulation of the evolution of surface topography on the substrate using information on ion and radical fluxes impinging on the substrate (which may be obtained from plasma simulations discussed above or experimental observations). Using such information, one can simulate the evolution of surface profiles during processing. First we discuss continuum surface models, in which boundaries of materials are represented by continuous surfaces. Although collections of atoms and molecules are actually involved, a continuous-surface representation is still a good approximation at dimensions currently of interest (i.e., 0.1-10 [muon]m, pattern widths) for microelectronic structures. To further simplify the problem, we consider only 2D problems, such as a traverse cross section of infinitely long trenches, where boundaries are expressed by curves. For the purpose of numerical simulation, such curves are discretized and represented by a sequence of nodes connected by line segments. The etching or deposition rates (i.e., "velocities" of the surface) are then evaluated at each nodal point on the basis of the magnitudes and angles of incoming ion and neutral fluxes. If the surface velocity (i.e., "rate") is known at each point on the surface, the equation that governs surface motion may be obtained in the following manner. The x-axis is assumed to be horizontal and the z-axis vertical. The y-coordinate, needed for 3D analyses, is chosen accordingly to form the usual right-hand coordinate system. We assume that the boundary curve of the material is given by an equation of the form [psi](x, z, t) = 0, and that [psi](x, z, t) > 0 (< 0) represents the material (vacuum) side of the boundary. Etching and deposition processes may then be formulated as the time evolution of this surface. Let C=C[SUB]n [nHAT]+C[SUB]t [tHAT] (16) be the velocity vector at point (x, z) of the boundary. Here [nHAT]=1 / [square root][psi][SUB]x**2 + [psi][SUB]z**2 ([psi][SUB]x [psi][SUB]z), [tHAT]=1 / [square root][psi][SUB]x**2 + [psi][SUB]z**2 (-[psi][SUB]z [psi][SUB]x) denote the unit normal and unit tangent vectors ([psi][SUB]x = [round dee][psi]/[round dee]x and [psi][SUB]z = [round dee][psi]/[round dee]z). Note that the normal velocity components C[SUB]n > 0 and C[SUB]n < 0 represent etching and deposition of the material in this formulation. Differentiating the equation [psi][x(t), z(t), t] = 0 with respect to t and substituting Equation (16) yields [psi][SUB]t + C[psi][SUB]x**2 + [psi][SUB]z**2 = 0. (17) Here we write C = C[SUB]n for simplicity. Note that the tangential velocity component C[SUB]t does not appear in Equation (17). The velocity along the curve does not alter its shape, and the motion of the curve is determined only by the normal component C [identity] C[SUB]n. Equation (17) is known as the Hamilton-Jacobi equation if the function C depends only on time t, position (x, z), the unknown function [psi], and its first derivatives ([psi][SUB]x, [psi][SUB]z) [77-80]. In general etching/deposition problems, however, C can be a more complicated function. If the height of the boundary curve z is a "function" of x and one can write z = u(x, t), we may write [psi](x, z, t) = u(x, t) - z and simplify Equation (17) as u[SUB]t + C[square root]1+u[SUB]x**2 = 0. (18) Suppose an ion beam bombards the surface in the negative z direction and the slope of the boundary curve is denoted by p = u[SUB]x = -[psi][SUB]x / [psi][SUB]z, as shown in Figure 7. Then the function f(p) = [square root]C1 + u[SUB]x**2 represents the volume of the material removed from the surface per unit beam flux and unit time and is proportional to the sputtering yield Y--i.e., the flux (number of atoms) sputtered from the surface per single incident beam particle (ion or atom). Therefore, f(p) is sometimes called a flux function [81, 82]. It is important to note that, although motion of the moving boundary always satisfies Equation (17) [or Equation (18)], the converse is not true. In other words, a solution to Equation (17) [or Equation (18)] satisfying given initial-boundary conditions does not necessarily describe physically plausible surface evolution. Indeed, Equation (17) [or Equation (18)] can admit two or more solutions for given initial-boundary conditions, and only one of them is physically meaningful. For an example of nonphysical solutions to Equation (17) [or Equation (18)], the reader is referred to Reference [4]. The problem of spurious solutions to Equation (17) [or Equation (18)] arises from the fact that we have tried to model surface evolution only by geometric motion of a surface (or a curve in the case of two dimensions) with given surface velocities. For example, a curve on a two-dimensional plane can intersect itself to form a loop during its motion, whereas the surface of a solid is always represented by a non-self-intersecting surface, no matter how the shape of the solid evolves. In this sense, there is an essential difference between geometric surface motion and the surface evolution of a material. To construct a system of equations that provides a unique solution for all t > 0 that represents physically meaningful material surface evolution, therefore, additional conditions known as "entropy conditions" must be imposed. Detailed discussion of these conditions is found in References [4, 81]. In many earlier simulations, each node and linear segments representing the material surface were moved in the direction of the specified surface velocity [often without referring to Equation (17) or Equation (18)]. For example, if the surface velocity is a function of the surface gradient p [identity] u[SUB]x (e.g., for the etching-rate function of Figure 9, discussed later), the surface velocity is not well defined at sharp corners (where p is not uniquely defined). In some earlier studies, ad hoc geometric considerations were employed to avoid such uncertainties. Consequently, the simulations suffered from a basic lack of credibility. If the surface velocity C does not depend on second- or higher-order derivatives of [psi] or u, Equation (17) and Equation (18) are of first order and can be solved by the method of characteristics. Several numerical simulation codes have been proposed employing this method. Problems still arise, however, when discontinuities of the surface gradient form (i.e., two characteristics intersect). Some existing models that employ the method of characteristics appear to require subtle geometric "adjustments," such as eliminating or avoiding the formation of nonphysical loops ("delooping") based on geometric hypotheses to deal with such situations (see for example References [83-92]). More recently, a numerical algorithm has been proposed for obtaining a solution for Equation (17) [or Equation (18)] with the entropy conditions, using the same curve discretization (i.e., nodal points connected by linear segments). The algorithm is designated the shock-tracking method, and through its use, a physically plausible solution can be obtained for the propagation of discontinuities of the surface slope (i.e., "shock wave"). For details, the reader is referred to Reference [81]. The shock-tracking method is employed in the simulation code "SHADE," which has been applied to various dry-etching and deposition problems [93-98]. The simulation results discussed next have been obtained using SHADE unless indicated otherwise. Most existing continuum-simulation codes employ some combination or modification of the boundary-moving algorithms mentioned above. Among them, SAMPLE [99], SPEEDIE [100], DEPICT [101], and EVOLVE [102] are examples of the most widely (and commercially) available 2D codes for simulating etching, deposition, and/or photolithography. There is also an entirely different class of methods that can be used to solve Equation (17) or Equation (18), without appealing to the characteristic or shock-tracking methods [103-105]. For example, one can solve Equation (17) [or Equation (18)] directly for [psi] (or u) on the x-z surface (or x axis), using a combination of a Lax-Wendroff-type finite-difference scheme with an upwind (or artificial diffusion) scheme. The "numerical" diffusion inherent in this method stabilizes the solution. The method is especially suited to the case in which the surface velocity depends explicitly on the local curvature and the PDE is of second order. As we discuss shortly, the curvature-dependent surface velocity introduces a "diffusion" effect, so that slope discontinuities are not formed in this case. Here we briefly discuss diffusive flows along the surface of the solid. If the temperature of the solid is sufficiently high (e.g., near its melting point), surface diffusion can cause changes in the surface topography. The surface diffusion is linearly dependent on the surface curvature [Kappa] and smooths out the surface [106]. To incorporate surface diffusion into Equations (17) and (18), we can replace C with C - [nu][Kappa], where -[nu][Kappa] is the surface velocity due to the curvature-driven diffusion ([nu] is typically a small positive constant). Since [Kappa]=-[psi][SUB]xx [psi][SUB]z**2 - 2[psi][SUB]xz [psi][SUB]x [psi][SUB]z + [psi][SUB]zz [psi][SUB]x**2 / ([psi][SUB]x**2 + [psi][SUB]z**2)**3/2 =u[SUB]xx / (1+u[SUB]x**2)**3/2 , Equations (17) and (18) become [psi][SUB]t + C[square root][psi][SUB]x**2 + [psi][SUB]z**2 (19) = -[nu][psi][SUB]xx [psi][SUB]z**2 - 2[psi][SUB]xz [psi][SUB]x [psi][SUB]z + [psi][SUB]zz [psi][SUB]x**2 / [psi][SUB]x**2 + [psi][SUB]z**2 and u[SUB]t + C[square root]1 + u[SUB]x**2 = [nu]u[SUB]xx (20) / 1+u[SUB]x**2. Solving Equation (17) or Equation (19) with the condition [psi](x, z, t) = 0 on a spatial grid is called a level-set method. Such a method makes it possible to simulate complex topological changes that occur during the course of surface evolution more easily than other methods based on line discretization. In the level-set method, however, the problem is solved in a space one dimension higher; e.g., [zeta] = [psi](x, z, t) is obtained from Equation (17) or Equation (19) first, and then the condition [psi](x, z, t) = 0 is used to obtain the solution curve on the x-z plane. Therefore, it is generally more expensive computationally than line-discretization-based methods. We next discuss the estimation of the surface velocity C in Equation (17) or (18). A typical schematic of gas-phase transport of ionic and neutral species is shown in Figure 8. Since the length scales of features in which we are interested are sufficiently smaller than the mean free paths, collisionless transport may be assumed for all gas-phase species. To simplify the following discussion, we also assume that the ion flux is unidirectional and vertical to the substrate surface. The sputtering yield Y depends on the angle [theta] formed by the incident flux and surface normal. Figure 9 shows an example of the angular dependence of an etching-rate function. The dependence is proportional to Y([theta]) cos [theta]. (It is normalized to 1 at [theta] = 0 in Figure 9.) Some etching-rate functions peak at nonzero angles [107], as in the case of the example shown in Figure 9. The velocity distribution function of the re-emitted atoms at position X on the boundary curve typically has a functional form given by f[SUP]R =f[SUB][nu] (|v|, X) cos[SUP][nu]+1 [theta], (21) where [theta] is the angle between the velocity vector v and the surface normal [Ncarat], as indicated in Figure 8. Depending on its angular dependence, the re-emission is designated as an over-cosine, cosine, or under-cosine distribution for [nu] > 0, [nu] = 0, or -1 < [nu] < 0, respectively. Such sputtered atoms have sufficiently low energies that they do not resputter the surface material and can be adsorbed on the surface when they land on other sites on the surface. The probability that an adsorbed atom remains on the surface indefinitely is the sticking coefficient [specialS]. The angular distribution of desorbed atoms is assumed to have the same functional form given in Equation (21). The incoming flux at position x on the surface boundary is related to the velocity distribution function f(v, x) as [specialF][SUP]in (x) = [integral][SUB][ncarat] o v<0[/SUB] ([ncarat] o v)f(v, x) dv, where [ncarat] is the unit vector normal to the surface at x and the integration is over the half velocity space ([ncarat] o v < 0). The outgoing flux is defined in a similar manner, with the integration over the other half space [ncarat] o v > 0. The outgoing flux from the surface, i.e., the flux of re-emitted atoms, is the sum of the desorption flux and resputtering flux: [specialF][SUP]out (x)=(1-[specialS])[specialF][SUP]in (x) (22) + [specialF][SUP]etch (x), where in(x) is the incoming flux of (low-energy) atoms. Therefore, the first term (1 - [specialS])[specialF][SUP]in represents the total desorption flux. The second term represents the outgoing flux due to the ion sputtering (etching), i.e.,[specialF][SUP]etch = Y([theta])[specialF][SUP]ion , where [specialF][SUP]ion is the magnitude of the unidirectional ion flux and [theta] is the angle formed by the surface normal and the vertical direction ([theta] = 0 in the example of Figure 8). Suppose that the only source of neutral atoms is the sputtering of the surface by ion bombardment. Using the re-emission distribution function, Equation (21), we may relate [4] the desorption flux [specialF][SUP]in (x) to the outgoing metal atom flux [specialF][SUP]out (X) at every position X (at which the unit normal vector is denoted as [Ncarat]) on the surface as [specialF][SUP]in (x) = A[SUB][nu]+4 ([nu]+2) / [pi] (23) [integral][SUB][specialB] dsK[SUB][nu] (x, X)[specialF][SUP]out (X), where K[SUB][nu](x, X)=g(x, X) ([Rhat] o [Ncarat])|[Rhat] (24) o [ncarat]|[SUP][nu]+1 / R is the integration kernel, the coefficient A[SUB][nu]+4 is given by A[SUB][nu]+4 = [Gamma]([nu]+3 / 2) [Gamma](1/2) / 2[Gamma]([nu]/2 + 2), and [specialB] denotes the surface boundary. For example, A[SUB]4 = [pi]/4 for the cosine distribution ([nu] = 0). In Equation (24), the function g is the visibility factor (i.e., g = 1 if the points x and X are on the line of sight, and g = 0 otherwise), R = x - X, R = |R|, and [Rhat] = R/R. Substituting Equation (22) into Equation (23) results in an integral equation for [specialF][SUP]in (x) for given source terms [specialF][SUP]etch (x). If the flux [specialF][SUP]in (x) of the incoming material at position x on the boundary curve is known, the deposition rate is given by D = (m/[rho])[specialF][SUP]in, where m and [rho] are the atomic mass and mass density of the deposited material. The net surface velocity is then the difference between the etching rate and the deposition rate. For more details of the derivation of Equation (23), the reader is referred to Reference [4]. To illustrate how such simulations can be applied to actual plasma deposition processes, we consider ionized-magnetron metal sputter-deposition processes [95, 96, 108, 109] as examples. An ionized-magnetron sputter-deposition system consists of a conventional magnetron sputter-deposition system and an inductively coupled plasma generator [108, 109]. Metal atoms (such as Cu) sputtered from the target are ionized by the high-density plasma (typically argon plasmas), which is generated by a radio-frequency induction (RFI) coil connected to an rf generator. The ion bombardment energy on the wafer surface, therefore the amount of resputtering, can be controlled through a bias voltage. The directionality of the ion flux in such a system provides better deposition at the bottom of a trench than conventional magnetron sputtering systems. With an appropriate bias voltage that induces simultaneous resputter deposition on the trench sidewalls, conformal deposition of thin metal liners over high-aspect-ratio trenches can be achieved. Ideally, plasma and sheath simulation discussed in the previous sections should be used to estimate energy and angular distributions of incoming ion fluxes and neutral-species fluxes at the substrate [110]. However, for the sake of simplicity, we assume here that the ion beam incident on the substrate is unidirectional and vertical and the neutral flux has a uniform angular distribution with a 52[degree] cutoff (i.e., collimation) angle due to the finite size of the target [98]. The etching-rate function is assumed to be proportional to the function of Figure 9, and the sticking coefficient [specialS] is assumed to be unity. The total neutral incoming flux [specialF][SUP]in consists of the resputtering flux [specialF][SUP]R defined by the right-hand side of Equation (23) and the direct neutral flux [specialF][SUP]D (i.e., [specialF][SUP]in = [specialF][SUP]R + [specialF][SUP]D). Since [specialF][SUP]D is known, we use Equation (22) (where [specialF][SUP]in is replaced by [specialF][SUP]R + [specialF][SUP]D) and Equation (23) (where the left-hand side is replaced by [specialF][SUP]R) to determine [specialF]R. Figure 10 shows simulation results for the deposition of a metal liner on the walls of trenches having aspect ratio 2.5. (Since Ar sputtering is taking place, we consider Y as the "effective" sputtering yield due to both metal and Ar ions. Namely, for a single ion impinging on the surface, Y metal atoms are sputtered, whether the sputtering is due to the metal ion or other Ar ions.) The ratio of metal ion flux to neutral flux is assumed to be 1:1. The shaded area represents the initial SiO2 trenches, and curves above it represent the profiles of the deposited film. At a sputtering yield of zero, very little deposition is expected. However, at a sputtering yield of 1.0, atoms sputtered from the trench bottom are expected to be redeposited on the vertical walls, predicting good conformal coverage. Figure 11 shows scanning electron microscopy (SEM) photographs of Cu deposition obtained under comparable conditions. The dc magnetron power was 10 kW for a 200-mm target cathode and up to 30 kW for a 300-mm cathode. The RFI coil was connected to a 13.56-MHz rf generator having a maximum power capacity of 3 kW. The bias voltage could be varied up to -200 V. Figure 12 depicts the results of a simulation of the effects of sputtering yield on metal deposition onto two trenches. Each boundary curve represents the metal film profile at equal time intervals, and the shaded area represents the initial trenches. As in Figure 10, the ratio of the ion to neutral fluxes was assumed to be 1:1. The sputtering yield at [theta] = 0 was assumed to be 0, 0.6, and 0.8. The simulation predicts that as the film thickness becomes comparable to the trench width, the atoms sputtered from one side of the trench will be collected on the opposite side, leading to a lateral buildup of resputtered deposit which can eventually result in the closing of the trench. Comparable effects are observed experimentally, as shown in Figure 13. The film depicted in part (a) was deposited with ion energies of 20 eV or less, which can cause a directional deposition without causing major resputtering at the surface. The film depicted in part (b) was deposited at an ion energy of about 80 eV, and that of part (c) at an ion energy of about 120 eV. As can be seen, both (b) and (c) show the effects of increasing sputtering yields. We now briefly discuss another entirely different numerical method for solving surface evolution problems in etching and/or deposition processes. Instead of representing the material surface by a continuous surface, one may represent both the material and incoming fluxes as collections of particles. As in PIC simulations, one can use superparticles, which represent a large number of ions and neutral species. For example, in SIMBAD [111], a commercially available simulation code, thin films of scale lengths of 1 [muon]m are represented by an aggregation of 10**4-10**5 disks in two dimensions. In contrast to PIC or molecular dynamics (MD) simulation, no force is assumed to be exerted on these disks during the gas-phase transport. When the disks reach the substrate, they may sputter the surface materials (also represented by disks), be adsorbed and desorbed with given probabilities (sticking coefficients), or undergo diffusion to minimize the surface chemical potential. One of the advantages of using particle (or disk) models instead of the continuum models discussed above is that microstructures (such as the packing density of particles) of the deposited film can also be taken into consideration. Figure 14 shows a sample calculation of MgF2 film deposition obtained from SIMBAD [112]. A disadvantage of the particle method may be its computational cost, because a very large number of particles must be used to evolve the surface, especially in 3D simulations. The particle method for surface evolution can be made more sophisticated by incorporating appropriate interparticle potentials for each simulation particle. Such an MD simulation represents a more realistic view of the microscopic process. This method makes it possible to estimate microstructure and its stresses with good accuracy. However, MD simulations are computationally intensive and costly [113]. Summary In this paper, we have reviewed some of the most widely used methods for modeling and simulating the plasma processes that are used in integrated-circuit fabrication. As is discussed in Sections 2 and 3, continuum models or particle models, or a combination of both (hybrid models), may be used to describe the dynamics of processing plasmas. The continuum models, such as those based on the drift-diffusion equations, are less time-consuming in terms of numerical computation but have several built-in approximations, such as a Maxwellian distribution for each species. The particle models are, on the other hand, less restrictive in terms of such built-in approximations and are therefore applicable to a plasma process under a wider variety of conditions; and they provide more detailed information on the system. The drawback is that such particle models generally require more computational resources (CPU time and storage) to simulate the plasma dynamics compared to a continuum model. If the plasma is highly collisional, a continuum model and its fluidlike simulation can provide accurate information on its dynamics. However, under certain conditions (for example, when the distribution functions are known to be highly non-Maxwellian), kinetic simulation may be the only feasible choice. After discussing sheath modeling briefly in Section 4, we have reviewed numerical methods to simulate the time evolution of microscopic surfaces under etching and/or deposition conditions. Although it is not presented in this paper, a combination of bulk plasma and sheath simulations with microscopic surface-evolution simulations can provide insights into the problem of controlling microscopic process results by adjusting macroscopic process-tool parameters such as gas pressures and bias voltages. Because of space limitations, we were able to give an account of only the fundamental aspects of plasma-processing modeling and simulation. Sample results of actual plasma processes were presented only to illustrate how such models and simulations can be applied to actual embodiments. The reader is therefore encouraged to refer to the referenced publications and also to the publications quoted therein for more details. Acknowledgment The author thanks S. M. Rossnagel for providing the SEM micrographs of Figures 11 and 13. References 1. B. N. Chapman, Glow Discharge Processes, John Wiley & Sons, Inc., New York, 1980. 2. O. A. Popov, Ed., High Density Plasma Sources, Noyes Publications, Park Ridge, NJ, 1995. 3. T. Uchida, Jpn. J. Appl. Phys. 33, L43 (1994). 4. S. Hamaguchi, in Thin Films: Modeling of Film Deposition for Microelectronic Applications, S. M. Rossnagel and A. Ulman, Eds., Academic Press, Inc., San Diego, 1996, Vol. 22, p. 81. 5. B. E. Thompson and H. H. Sawin, J. Appl. Phys. 60, 89 (1986). 6. D. B. Graves and K. F. Jensen, IEEE Trans. Plasma Sci. 14, 92 (1986). 7. A. D. Richards, B. E. Thompson, and H. H. Sawin, Appl. Phys. Lett. 50, 492 (1987). 8. D. B. Graves, J. Appl. Phys. 62, 88 (1987). 9. M. S. Barnes, T. J. Cotler, and M. E. Elta, J. Appl. Phys. 61, 81 (1987). 10. S.-K. Park and D. J. Economou, J. Appl. Phys. 68, 3904 (1990). 11. N. Sato and Y. Shida, Jpn. J. Appl. Phys. 36, 4794 (1997). 12. T. H. Chung, L. Ming, H. J. Yoon, and J. K. Lee, Jpn. J. Appl. Phys. 63, 2874 (1997). 13. W. Lee, S. E. Laux, M. V. Fischetti, G. Baccarani, A. Gnudi, J. M. C. Stork, J. A. Mandelman, E. F. Crabbe, M. R. Wordeman, and F. Odeh, IBM J. Res. Develop. 36, 208 (1992). 14. W. M. Coughran, Jr., J. Cole, P. Lloyd, and J. K. White, Eds., Semiconductors, Part I (the IMA volumes in mathematics and its applications, Vol. 58), Springer-Verlag, New York, 1994. 15. W. M. Coughran, Jr., J. Cole, P. Lloyd, and J. K. White, Eds., Semiconductors, Part II (the IMA volumes in mathematics and its applications, Vol. 59), Springer-Verlag, New York, 1994. 16. H. Gummel, IEEE Trans. Electron Devices ED-11, 455 (1964). 17. J.-P. Boeuf, Phys. Rev. A 36, 2782 (1987). 18. K. Miyamoto, Plasma Physics for Nuclear Fusion, MIT Press, Cambridge, MA, 1980. 19. K. Okazaki, T. Makabe, and Y. Yamaguchi, Appl. Phys. Lett. 54, 1742 (1989). 20. J.-P. Boeuf, J. Appl. Phys. 63, 1342 (1988). 21. K. Satake, T. Monaka, O. Ukai, Y. Takeuchi, and M. Murata, Jpn. J. Appl. Phys. 36, 4789 (1997). 22. V. A. Godyak and R. B. Piejak, Phys. Rev. Lett. 65, 996 (1990). 23. V. A. Godyak, R. B. Piejak, and B. M. Alexandrovich, IEEE Trans. Plasma Sci. 19, 660 (1991). 24. V. A. Godyak, R. B. Piejak, and B. M. Alexandrovich, Plasma Sources Sci. Technol. 1, 36 (1992). 25. N. A. Krall and A. W. Trivelpiece, Principles of Plasma Physics, McGraw-Hill Book Co., Inc., New York, 1973. 26. W. N. Hitchon, D. J. Koch, and J. B. Adams, J. Comp. Phys. 83, 97 (1989). 27. T. J. Sommerer, W. N. G. Hitchon, R. E. P. Harvey, and J. E. Lawler, Phys. Rev. A 43, 4452 (1991). 28. G. J. Parker and W. N. Hitchon, Jpn. J. Appl. Phys. 36, 4799 (1997). 29. C. K. Birdsall and A. B. Langdon, Plasma Physics Via Computer Simulation, McGraw-Hill Book Co., Inc., New York, 1985. 30. T. Tajima, Computational Plasma Physics, Addison-Wesley Publishing Co., New York, 1989. 31. C. K. Birdsall, IEEE Trans. Plasma Sci. 19, 65 (1990). 32. D. Vender and R. W. Boswell, IEEE Trans. Plasma Sci. 18, 725 (1990). 33. M. Surendra and D. B. Graves, IEEE Trans. Plasma Sci. 19, 144 (1991). 34. V. Vahedi, M. A. Lieberman, M. V. Alves, J. P. Verboncoeur, and C. K. Birdsall, J. Appl. Phys. 69, 2008 (1991). 35. M. V. Alves, M. A. Lieberman, V. Vahedi, and C. K. Birdsall, J. Appl. Phys. 69, 3823 (1991). 36. M. M. Turner and H. B. Hopkins, Phys. Rev. Lett. 69, 3511 (1992). 37. V. Vahedi, G. DiPeso, C. K. Birdsall, M. A. Lieberman, and T. D. Rognlien, Plasma Sources Sci. Technol. 2, 261 (1993). 38. H. R. Skullerud, J. Phys. D 1, 1567 (1968). 39. S. L. Lin and J. N. Bardsley, J. Chem. Phys. 66, 435 (1977). 40. T. J. Sommerer and M. J. Kushner, J. Appl. Phys. 71, 1654 (1992). 41. V. Vahedi, C. K. Birdsall, M. A. Lieberman, G. DiPeso, and T. D. Rognlien, Plasma Sources Sci. Technol. 2, 273 (1993). 42. M. A. Lieberman, IEEE Trans. Plasma Sci. 16, 638 (1988). 43. M. A. Lieberman, IEEE Trans. Plasma Sci. 17, 338 (1989). 44. V. Vahedi, C. K. Birdsall, M. A. Lieberman, G. DiPeso, and T. D. Rognlien, Phys. Fluids B 5, 2719 (1993). 45. V. Vahedi and G. DiPeso, J. Comp. Phys. 131, 149 (1997). 46. Ph. Belenguer and J.-P. Boeuf, Phys. Rev. A 41, 4447 (1990). 47. J.-P. Boeuf and L. C. Pitchford, IEEE Trans. Plasma Sci. 19, 286 (1991). 48. N. Sato and H. Tagashira, IEEE Trans. Plasma Sci. 19, 102 (1991). 49. M. Surendra, D. B. Graves, and G. M. Jellum, Phys. Rev. A 41, 1112 (1990). 50. K. H. Shoenbach, H. Chen, and G. Schaefer, J. Appl. Phys. 67, 154 (1990). 51. R. K. Porteous and D. B. Graves, IEEE Trans. Plasma Sci. 19, 204 (1991). 52. M. Surendra, D. B. Graves, and L. S. Plano, J. Appl. Phys. 71, 5189 (1992). 53. R. A. Stewart, P. Vitello, and D. B. Graves, J. Vac. Sci. Technol. B 12, 478 (1994). 54. G. A. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows, Oxford Science Publishers, Oxford, UK, 1994. 55. D. J. Economou, T. J. Bartel, R. S. Wise, and D. P. Lymberopoulos, IEEE Trans. Plasma Sci. 23, 581 (1995). 56. J. Johannes, T. Bartel, G. A. Hebner, J. Woodworth, and D. J. Economou, J. Electrochem. Soc. 144, 2448 (1997). 57. S. Rauf and M. J. Kushner, J. Appl. Phys. 81, 5966 (1997). 58. D. P. Lymberopoulos and D. J. Economou, IEEE Trans. Plasma Sci. 23, 573 (1995). 59. V. P. Gopinath and T. A. Grotjohn, IEEE Trans. Plasma Sci. 23, 602 (1995). 60. M. M. Turner, Jpn. J. Appl. Phys. 36, 4784 (1997). 61. T. Morimoto, Y. Yasaka, M. Tozawa, T. Akahori, H. Amano, and N. Ishii, Jpn. J. Appl. Phys. 36, 4769 (1997). 62. H. Muta, Y. Ueda, and Y. Kawai, Jpn. J. Appl. Phys. 36, 4773 (1997). 63. R. J. Hoekstra, M. J. Grapperhaus, and M. J. Kushner, J. Vac. Sci. Technol. A 15, 1913 (1997). 64. P. L. G. Ventzek, R. J. Hoekstra, and M. J. Kushner, J. Vac. Sci. Technol. B 12, 461 (1994). 65. A. P. Paranjpe, J. Vac. Sci. Technol. A 12, 1221 (1994). 66. H.-M. Wu, B. W. Yu, M. Li, and Y. Yang, IEEE Trans. Plasma Sci. 25, 1 (1997). 67. W. M. Davis and T. A. Vanderslice, Phys. Rev. 131, 219 (1963). 68. M. J. Kushner, J. Appl. Phys. 58, 4024 (1985). 69. B. E. Thompson, H. H. Sawin, and D. A. Fisher, J. Appl. Phys. 63, 2241 (1988). 70. R. T. Farouki, S. Hamaguchi, and M. Dalvie, Phys. Rev. A 44, 2664 (1991). 71. G. H. Wannier, Bell Syst. Tech. J. 32, 170 (1953). 72. J. E. Lawler, Phys. Rev. A 32, 2977 (1985). 73. S. Hamaguchi, R. T. Farouki, and M. Dalvie, Phys. Rev. A 44, 3804 (1991). 74. S. Hamaguchi, R. T. Farouki, and M. Dalvie, Phys. Rev. Lett. 68, 44 (1992). 75. S. Hamaguchi, R. T. Farouki, and M. Dalvie, Phys. Fluids B 4, 2362 (1992). 76. R. T. Farouki, S. Hamaguchi, and M. Dalvie, Phys. Rev. A 45, 5913 (1992). 77. E. D. Conway and E. Hopf, J. Math. Mech. 13, 939 (1964). 78. M. G. Crandall and P. L. Lions, Trans. Amer. Math. Soc. 277, 1 (1983). 79. M. G. Crandall, L. C. Evans, and P. L. Lions, Trans. Amer. Math. Soc. 282, 487 (1984). 80. M. G. Crandall and P. L. Lions, Math. Comp. 43, 1 (1984). 81. S. Hamaguchi, M. Dalvie, R. T. Farouki, and S. Sethuraman, J. Appl. Phys. 74, 5172 (1993). 82. P. D. Lax, Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves, SIAM Regional Conference Series in Applied Mathematics, Society for Industrial and Applied Mathematics, Philadelphia, 1973. 83. M. J. Nobes, J. S. Colligon, and G. Carter, J. Mater. Sci. 4, 730 (1969). 84. G. Carter, J. S. Colligon, and M. J. Nobes, J. Mater. Sci. 6, 115 (1971). 85. D. J. Barber, F. C. Frank, M. Moss, J. W. Steeds, and I. S. T. Tsong, J. Mater. Sci. 8, 1030 (1973). 86. J. P. Ducommun, M. Cantagrel, and M. Marchal, J. Mater. Sci. 9, 725 (1974). 87. R. Smith, G. Carter, and M. J. Nobes, Proc. Roy. Soc. Lond. A 407, 405 (1986). 88. R. Smith, S. J. Wilde, G. Carter, I. V. Katardjiev, and M. J. Nobes, J. Vac. Sci. Technol. B 5, 579 (1987). 89. I. V. Katardjiev, J. Vac. Sci. Technol. B 7, 3227 (1989). 90. I. V. Katardjiev, G. Carter, and M. J. Nobes, J. Phys. D 22, 1813 (1989). 91. C. W. Jurgensen and E. S. G. Shaqfeh, J. Vac. Sci. Technol. B 7, 1488 (1989). 92. J. C. Arnold, H. H. Sawin, M. Dalvie, and S. Hamaguchi, J. Vac. Sci. Technol. A 12, 620 (1994). 93. S. Hamaguchi and M. Dalvie, J. Vac. Sci. Technol. A 12, 2745 (1994). 94. S. Hamaguchi and M. Dalvie, J. Electrochem. Soc. 141, 1964 (1994). 95. S. Hamaguchi and S. M. Rossnagel, J. Vac. Sci. Technol. B 13, 183 (1995). 96. S. Hamaguchi and S. M. Rossnagel, J. Vac. Sci. Technol. B 14, 2603 (1996). 97. S. Hamaguchi, A. A. Mayo, S. M. Rossnagel, D. E. Kotecki, K. R. Milkove, C. Wang, and C. E. Farrell, Jpn. J. Appl. Phys. 36, 4762 (1997). 98. A. A. Mayo, S. Hamaguchi, J. H. Joo, and S. M. Rossnagel, J. Vac. Sci. Technol. B 15, 1788 (1997). 99. W. G. Oldham, A. R. Neureuther, C. Sung, J. L. Reynolds, and S. N. Nadgeonkar, IEEE Trans. Electron Devices ED-27, 1455 (1980). 100. J. Ignacio, F. Ulacia, and J. P. McVittie, J. Appl. Phys. 65, 1484 (1989). 101. T. Thurgate, IEEE Trans. Computer-Aided Design IC's & Syst. 10, 1101 (1991). 102. T. S. Cale and V. Mahadev, in Thin Films: Modeling of Film Deposition for Microelectronic Applications, S. M. Rossnagel and A. Ulman, Eds., Academic Press, Inc., San Diego, 1996, Vol. 22, p. 175. 103. D. S. Ross, J. Electrochem. Soc. 135, 1235 (1988). 104. D. S. Ross, J. Electrochem. Soc. 135, 1260 (1988). 105. S. Osher and J. A. Sethian, J. Comp. Phys. 79, 12 (1988). 106. Jean Philibert, Atom Movements Diffusion and Mass Transport in Solids, Les Editions de Physique, Les Ulis Cedex A, France, 1991. 107. P. C. Zalm, in Handbook of Ion Beam Processing Technology, J. J. Cuomo, S. M. Rossnagel, and H. R. Kaufman, Eds., Noyes Publications, Park Ridge, NJ, 1989, p. 78. 108. S. M. Rossnagel and J. Hopwood, Appl. Phys. Lett. 63, 3285 (1993). 109. S. M. Rossnagel and J. Hopwood, J. Vac. Sci. Technol. B 12, 449 (1994). 110. M. J. Grapperhaus, Z. Krivokapic, and M. J. Kushner, J. Appl. Phys. 83, 35 (1998). 111. M. J. Brett, S. K. Dew, and T. J. Smy, in Thin Films: Modeling of Film Deposition for Microelectronic Applications, S. M. Rossnagel and A. Ulman, Eds., Academic Press, Inc., San Diego, 1996, Vol. 22, p. 2. 112. R. N. Tait, T. Smy, and M. J. Brett, Thin Solid Films 187, 375 (1990). 113. C.-C. Fang, V. Prasad, R. V. Joshi, F. Jones, and J. J. Hsieh, in Thin Films: Modeling of Film Deposition for Microelectronic Applications, S. M. Rossnagel and A. Ulman, Eds., Academic Press, Inc., San Diego, 1996, Vol. 22, p. 117. Received February 23, 1998; accepted for publication July 29, 1998 Satoshi Hamaguchi Department of Fundamental Energy Science, Graduate School of Energy Science, Kyoto University, Gokasho, Uji-shi, Kyoto 611-0011, Japan (hamaguch@energy.kyoto-u.ac.jp). Dr. Hamaguchi is an Associate Professor of Energy Science at Kyoto University. This paper was written while he was a Research Staff Member at the IBM Thomas J. Watson Research Center in Yorktown Heights, New York. After receiving B.S. and M.S. degrees in physics from the University of Tokyo, Japan, in 1982 and 1984, he came to the United States under a Fulbright scholarship to continue his graduate studies at the Courant Institute of Mathematical Sciences, New York University. Upon graduation in 1988, he joined the Institute for Fusion Studies at the University of Texas at Austin, where, as a Research Fellow, he worked on thermonuclear fusion plasmas until joining IBM in 1990. He joined the faculty of Kyoto University in 1998. His current research interests include the kinetic theory of low-temperature plasmas and of thermonuclear fusion plasmas, as well as the applications of plasma etching and deposition processes in the fabrication of LCDs and integrated circuits. Professor Hamaguchi holds a Ph.D. in physics from the University of Tokyo and a Ph.D. in mathematics from New York University.