PHYSICAL REVIEW E VOLUME 49, NUMBER 4 APRIL 1994 Monte Carlo simulation of electron swarms at low reduced electric fields M. Yousfi, A. Hennad, and A. Alkaa Universit'e Paul Sabatier, Centre de Physique Atomique, 118, Route de Narbonne, 31 062 Toulouse Cedex, France (Received 2 August 1993) A Monte Carlo method, based on the classical null collision technique, is developed and applied for electron swarm simulation in weakly ionized gases at low E/N values (E being the electric field and iV the background gas density). The dominant low energy collisional processes (elastic, inelastic, and su-perelastic collisions) are properly considered. The influence of the thermal motion of background gas (target not assumed at rest) on electron kinetics, which can play an important role at very low E/N, is also taken into account. At zero field and low E /N conditions, the influence of low energy collision processes on space and time evolution of electron swarm parameters is emphasized in atomic (He) and molecular (H20) gases. Then, the present Monte Carlo method is adapted to the low E/N cases in strongly electronegative gases such as SF6 to overcome the problem of the vanishing of seed electrons by attachment processes which can stop the simulation. In that case, an additional fictitious ionization process with constant collision frequency is considered to obtain hydrodynamic electron swarm parameters in SF6 at low E /N. PACS number(s): 52.20.Fs, 51.50.+ v I. INTRODUCTION Swarm characteristics (distribution functions and transport coefficients) of electrons moving in a weakly or partially ionized gas under the influence of an accelerating electric field can be numerically obtained either from the direct Boltzmann equation solution or from Monte Carlo simulation. Solution of the Boltzmann equation in a nonhydrodynamic regime involving time and in particular space variation of distribution functions constitutes a numerical challenge, while the Monte Carlo method is easier to develop in hydrodynamic as well as nonhydrodynamic regimes. However, the well known Monte Carlo drawback concerns the prohibitive computing time. But, due to progress in computer technology, calculation time is less and less a constraint. Monte Carlo methods can thus be considered as well adapted for numerical simulations of electron swarm motion in gas discharges whatever the geometry (one, two, or three dimensions), the regime (transient or steady, homogeneous or not), and the applied electric field. Concerning the applied electric field E/N (N being the gas density), there are some specific problems in Monte Carlo simulation which depend, in particular, on E /N magnitude. At high E/N values, the large electron amplification due to ionizing collisions obviously enhances the number of electrons treated with Monte Carlo method and can considerably increase computing time. Such a problem—not emphasized in this paper—could be solved by using appropriate weighting techniques. At low E/N values, in the case of electron-molecule collisions, the first problem concerns the correct treatment of collision kinematics. This must include the motion of projectile electrons as well as the thermal motion of target molecules. The motion consideration of both projectile and target particles is important for low impact energy, i.e., at low E/N values when energies of impinging electrons and target molecules are practically of the same order of magnitude. This is what happens at low electron energy not only during elastic collisions but also during inelastic and superelastic collisions involving rotational and vibrational energy levels of target molecule. For these cases, the correct treatment of the collision kinematics (i.e., without assuming molecule at rest) can avoid some errors physically not acceptable in swarm parameter determination at low E /N such as an electron mean energy lower than gas energy. These kinds of non-physical Monte Carlo results must also be avoided or minimized because electron swarm data can be used either as input data in fluid equations for discharge modeling or in methods for fitting electron-molecule collision cross sections by unfolding swarm parameters. At low E/N, there is another problem which can appear in the case of highly electronegative gases such as SF6. The strong electron attachment occurring at low energy (lower than 0.5 eV for SF," ion formation by electron impact) can absorb enough initial seed electrons to stop Monte Carlo code. The purpose of this paper is to present a Monte Carlo method available whatever E/N, but which is more adapted for low E /N values (from 0 Td up to some tens of Td). In this method, described in Sec. II, the energy exchanged between electrons and molecules during elastic, inelastic, and superelastic impacts is properly taken into account. Then, the cross sections used and the analysis of corresponding results, showing the validity of the present methods at low E/N, are given in Sec. Ill in the case of electropositive (He) and slightly electronegative gases (H20). Finally, in the particular case of a strongly electronegative gas (SF6), electron swarm parameters are calculated using an improved Monte Carlo method based on an additional fictitious ionization process in order to avoid the problem of the vanishing of seed electrons due to the high electron attachment 1063-651X/94/49(4)/3264(10)/$06.00 49 3264 © 1994 The American Physical Society 49 MONTE CARLO SIMULATION OF ELECTRON SWARMS 3265 efficiency in SF6 at low E/N. It is to be noted that a symmetric idea (based on fictitious attachment process) was already applied to N2 at high E/N strength by Li, Pitch-ford, and Moratz [1] in order to reduce the number of simulated electrons and therefore the computing time. II. MONTE CARLO METHOD The electron transport in a gas under the influence of an electric field E can be simulated with the help of a Monte Carlo method from an initially great number of seed electrons. These primary electrons are treated one by one from their creation until their disappearance out of the domain of the simulation or by specific collisional processes (e.g., attachment). Every electron, during its transit in the gas, performs a succession of free flights punctuated by elastic, inelastic, or superelastic collisions with molecules 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, distribution functions and transport coefficients. The simulation is stopped when all the primary electrons as well as the secondary electrons (created, for example, by ionization) are treated. The flow chart of Monte Carlo method described hereafter is shown in Fig. 1. As it can be seen, after definitions of the simulation parameters, the gas, and the initial conditions, it is necessary to know first the time of free flight. A. Time of free flight f aigh, The time of free flight is calculated by using the null collision method initially developed by Skullerud [2] for simulation of ion motion in gases and then used by numerous authors [3]: 'flight ' (1) where r&^t is a random number uniformly distributed in the [0,1] range and ytot is the total collision frequency including total electron-molecule collision frequency v and a null collision frequency unull chosen in order to have always utot constant: = v +v null = const (2) B. Trajectory between two successive collisions The trajectory between two successive collisions is obtained from the classical mechanic equations. In the framework of this paper, the electric field accelerating electrons (with charge —e, mass m, position r, and velocity v) is assumed to be antiparallel to the z axis. Under these conditions, the components vxl, vyl, and vzl in the laboratory frame of velocity VjUj) at time 11 at the end of the free flight can be written as a function of velocity v0(f0) (at initial time r0 and with components vx0, vy0, and vz0 in the laboratory frame). Then, new coordinates rj(x!,>>!,zx) of electron at time tl can be calculated from coordinates r0(x0,y0,z0) of an electron at time r0: Xi=x0 + vx0tRight , J'l^o + Vyo'flight . <3' 1 eEz 2 z\— Z0+Vz0{ flight + j m 'flight' with ?flight = f l ~ ro- So, starting from the electron parameters r0,v0,r0 at the beginning of the free flight, the new electron parameters ti,V[,r, at the end of the free flight are obtained, respectively, from relations (1) and classical mechanic equations. Then just after collision occurring at time 11, electron parameters become t\,v\,x\. However, it is Input: Choice of gas (cross sections) and simulation parameters: Vmax , l*max or t max , E/N , T and number of electron seeds Parameters of primary electron Vo.ro , to -1 Parameters of secondary electron vo =v5 , ro = rs, to= ts Electron parameters before collision: - Ifljgh, or ti from a random number - T1 , VI from classical mechanics equation Molecular velocity Vl from random numbers ==> Relative velocity Vr = VI - Vt_ Next Collision Storage of electron parameters for distribution function and swarm parameter calculations 4 Input parameters for next collision V0 = V1 ro = n to = ti Type of collision!*" . Real collisionJ** Yes Storage of electron parameters V s , r s , of ejected electron "^^Attachment ? ^0^>- Calculation of deviation and azimuth angles of scattered electron from random numbers Calculation of velocity v'i of scattered electron after the collision Input parameters for next collision VO - V'i ro - n to = ti Output: Electron distribution function and swarm parameters FIG. 1. Flow chart for Monte Carlo simulation of electron swarm motion in a gas under action of an accelerating electric field E (velocity umax, position rmax, and time tm&x correspond to the possible limits of the simulation domain: see Sec. II for an explanation of the other notations used). 3266 M. YOUSFI, A. HENNAD, AND A. ALKAA 49 necessary to calculate only electron velocity vj because the electron-molecule interaction is assumed to be instantaneous (t\ =ti) and local (rj =r,). In order to calculate the velocity vj, it is necessary to know the collision type. C. Type of collision The collision type necessitates knowledge of the likelihoods (pcoi,ei»PCoi>in»/'a)i,sup> or />col,null) of the occurrence of every collision kind (elastic, inelastic, superelastic, or null): P col, el" /> col.ii (4) Pcol,! J sup _ utui\\ • P col,null ~ with "col,el ^-Pcol.in "'"/'col.sup +/'col,null 1 • In fact, collision frequencies for the different processes (ud,Din,usup) depend on the relative velocity vr before the collision, which is denned as vr = |v,—Vj; V] is the known velocity of target molecule having a Maxwellian distribution. The collision type is then determined from a random number rcol uniformly distributed between 0 and 1. Several types of collision are possible: (i) if it is a null collision, velocities before and after the collision are the same; (ii) if it is, for instance, an attachment, the next primary electron is treated; and (iii) if it is another real collision, the velocity vj after interaction depends on the collision type; the components of electron velocity vj are given hereafter. D. Velocity vj after a real collision Let m and M be electron and molecule masses, vl and Vj their respective velocities in the laboratory frame before the collision, and vj and Vj after the collision. Unknown velocities v'j and V, can therefore be determined from the classical conservation equation of momentum transfer, yielding M m+M v, = - m m+M m+M m +M (5a) (5b) Of the right-hand term of relations (5), only the velocity (mv,+MVi)/(ffl +M) of the center-of-mass frame is already known. The unknown vector v'r (v'r = vj — Vj) represents the relative velocity after the collision. The modulus v'r and the direction of vector v'r are given hereafter. The relative speed v'r, which is obtained from the classical conservation equation of total energy, depends on collision type. For an elastic collision For an excitation of rotational, vibrational, or optical level from the lower level i with potential energy e, to the upper level j with energy e; 1 1/2 (6b) For a superelastic processes corresponding to deexcita-tion from the upper level j with potential energy e; to the lower level i with energy e, vl + — Ae,- Mr (6c) fir is the reduced mass: nr = mM/(m+M), and Then the knowledge of the scattering angle X and tne azimuthal angle 17 can give the direction of the vector v'r in the center-of-mass frame. The deflection angle x> between the relative velocity vr before and v). after a collision, varies between 0 and tt. It depends on the differential cross section a(v,x) and is determined from a uniform random number rx belonging to the [0,1] range: ' x fXa(v,x' hinx'dx' J n_ f"a(v,x')sinx'dx' (7a) When an isotropic scattering is assumed, this relation becomes cos#= \ — 2rx. The azimuthal angle 17 can be also assumed to be uniformly distributed in the [0,2ir] range and is calculated from a uniform random number r-. (7b) (6a) So, as the relative speed v'r is determined from relations (6) of the conservation of total energy and the x and f] angles from relations (7), the vector v'r is therefore completely defined in the center-of-mass frame. But in relations (5), the vector v'r is needed in the laboratory frame; such a transformation is undertaken using the classical Euler relations v'rx =v'r{ — sinx sinr/ sin^r +sin^ cost? cos0rcos0r + cosx sin8rcos4>r) , v'^ = v'r(sinx sinn cos<^>r + sin^ COS17 cos<9rsin0,. + cosjsin0rsin^r) , v'rz=v'r(— sinx C0SV sin0r + cosx cos#r) , where 6r is the polar angle and r the azimuthal angle in the laboratory frame of vector vr. Then, from relation (5a), it is possible to calculate electron velocity vj in the laboratory frame without neglecting, as usual the velocity V! of target molecule; the vector V! is of course completely defined because the distribution of background molecular gas is a known Maxwellian. Noting that if the target molecule is considered to be at rest (V1=0) and m negligible in comparison to M (i.e., m +M=M), relations (5a) and (6) reduce to the classical relations usually used in the literature for collision treatment in Monte Carlo method. 49 MONTE CARLO SIMULATION OF ELECTRON SWARMS . . . 3267 Concerning ionization processes, the previous relation (6b) used for excitation processes is still valid, but is not necessary. The reason is that ionization processes involve generally an energy amount much higher than the energy gas, so that the assumption of a background molecule at rest can be considered as a good approximation. Following this approximation, the residual energy, after a simple ionization process with threshold Ejon, is shared between scattered (e\) and ejected (cej) electrons following the relation E\ + Eej = jmv\ — eion. The energy sharing depends on the knowledge of the differential ionization cross section aion(eltz). As ejected and scattered electrons are not discernable, the energy of one of the two electrons (eej, for instance) can be obtained from a uniform random number rion from ion „ / „ \ y where 0\on(zx) is the integral ionization cross section and ex the incident energy (£j = -l-/nvf). The scattered electron, after an ionizing collision, is then deflected following an angle x> assumed to be isotropic, and the ejected electron is deflected following an angle x', assumed to be orthogonal to the X direction (#'=j+7r/2). Previous relations give the components of the velocity vectors of scattered v', (v\,x,rj) and ejected vey-(ve;-,^',ij) electrons in the center-of-mass frame. The corresponding components in the laboratory frame are then determined using the classical Euler transformation, which can be written in the case of, for example, the scattered electron with velocity v'(. v'xi—v\(— sin;^ sin»7 sin<£, + smx cost? cos0, cos0, +cos^'sin0,cos^l) , v'yi= v\{smx sin?7 cos^ + sinx cos-n cos^sin^ + cosxsin0,sin^1) , v'zl = v'j( — sin^ cost; sin0] +cos;f cos^) , where 6{ is the polar angle and ,(t,v,z), representing the isotropic part (/=0) and successive anisotropies (/>0) of the distribution function f(t,v,z), is directly obtained from a Monte Carlo simulation similarly to Penetrante, Bardsley, and Pitchford [4]. In other words, the energy e domain [0,emax], where e is the electron incident energy (E = \mv1) and emax its maximum value, is first divided into nv regular intervals: [0,emax] = [e0=:0,£1, . . . ,Ek-lEk,£k + l, . . . , e„v =emax] with constant energy step Ae = ek+ j — e^ . Then a discrete function Q>i(t,zk + l/2,z) is defined such as */(r,et + 1/2,z)= 2 2 6(Ei7)P;(cos0y; ;=i/=i (9) with 8(Ey)=l if zk 0(t,E,z) of the distribution function and anisotropy <^;(f,v,z) of order / can then be obtained from the following relation [4]: 4>i(t>Zk + \n,z)= .-~-r—;- . (10) n ne{t,z) where ne(t,z) is the electron number density used as normalization quantity for the distribution function E,t + l/2>z)- Concerning transport coefficients such as drift velocity W, longitudinal DL or transverse DT diffusion coefficients, mean energy (e), ionization (vion) or attachment frequency (vatt), etc., they are then calculated using a statistical mean based on conventional formulas. For instance, the time evolution of transport coefficients is obtained by first discretizing the time domain [0, tmax ] in nt regular intervals: [0,fmax] = [f0=0, tu ...,tm_x, tm, tm+l,. . . ,tnt=tmax] with constant time step to = tm + 1-tm. Then the relations for transport coefficients such as (e), W, (vion), or nj,5 + \n ; = 1 n/,s + l/2 , = 1 (12) where £yiJ + 1/2 is the kinetic energy for an electron number j undergoing collision number i in the space interval [zs,zs + l], while ny>s + i/2 represents the electron number counted in the interval [zs,zs + ]] and nts + l/2 the collision number undergone by the electron j in the same space interval. III. RESULTS AND DISCUSSIONS In Sees. Ill A and HIB, calculations of distribution functions and transport coefficients are carried out in molecular gases (H20) and atomic gas (He) chosen in order to check first the validity of the present Monte Carlo method and then to give insight into electron swarm data in the case of dominant low E/N collisional processes (i.e., thermal motion of target, inelastic, and superelastic collisions). Section III C is devoted to the case of strongly attaching gases. The set of electron-He cross sections is taken from the literature (see, e.g., [5]). Collision cross sections chosen for H20 are already partly fitted elsewhere [6] by comparing measured and calculated transport coefficients using a multiterm Boltzmann equation solution [8]. This set of electron-H20 cross sections includes elastic momentum transfer, excitations of ten rotational, two vibrational, and nine optical levels, and also ionization and attachment processes. Each electron-molecule collision cross section ay,(e) for deexcitation (superelastic processes) through the upper level j to the allowed lower level / is determined from the excitation cross section atj■{e) by using the well known principle of detailed balance (see, e.g., [9]): g, e + Ae,-- ah■ (e) =--er,,( e + Ae,-.-) , gj £ where g, and gj denote the statistical weight of the levels i and j, respectively [Ae;j is already defined in relations (6)]. A. Zero-field electron swarm results Probably one of the most convincing validity tests of the treatment of low energy electron-molecule collisions with the Monte Carlo method is to determine distribution functions and transport coefficients under zero-field conditions. Indeed, for an electron swarm or beam released (with known initial energetic and angular distributions) through a gas under zero-electric field conditions, it is well established that this electron swarm relaxes after a greater or lesser period of time (depending on initial conditions and background gas) towards an equilibrium distribution, whatever the initial distribution or the nature of the background gas [10]. Such an equilibrium is obviously characterized by a classical behavior: the electron distribution function becomes Maxwellian at the background gas temperature, there is no more electron drift, and diffusion becomes completely isotropic (i.e., longitudinal and transverse diffusion coefficients are identical). Such a behavior is perfectly illustrated in Figs. 2(a) and 2(b) showing electron mean energy, drift velocity, and also longitudinal and transverse diffusion coefficients in the case of an energetic electron swarm released in H20 background gas. Figure 2(a), showing an electron mean energy which relaxes towards gas energy whatever the initial electron energy (lower or higher than gas energy), corresponds to a long time scale. In Fig. 2(b), corresponding to shorter time scale, electrons emitted along the forward direction, after a relatively few collisions, lose their initial anisotropic angular distribution so that the initial directed velocity becomes rapidly negligible. In this short time scale, the longitudinal diffusion coefficient, after an overshoot effect due to the anisotropy of the initial distribution, tends towards transverse diffusion. Figure 3 shows another validity test corresponding to the zero-field electron mean energy in He calculated with and without including the thermal motion of background gas. Figure 3 clearly illustrates the consequence for assuming a target at rest since the electron mean energy tends towards zero energy instead of to gas energy (i.e., \kT). These results, completed with further validity tests not reported in this paper and undertaken in other gases, give us a good reliability concerning the Monte Carlo method described in Sec. II for the treatment of low energy collision processes (energy exchange during elastic, inelastic, and superelastic collisions). B. Low E /N electron swarm results Figures 4(a) and 4(b) show the electron mean energy ( e ) and drift velocity W as a function position z along the direction of the electric field acceleration in the case of two relatively low E/N values [31 and 93 Td OO"17 V cm2)] background gas. Such Monte Carlo calculations, which can correspond to the simulation of the classical steady-state Townsend experiment (see, e.g., [8]) are undertaken with and without including the effects of superelastic collisions in order to emphasize the influence of these low energy processes on electron transport. So, starting from the cathode (z =0) with an initial mono- 49 MONTE CARLO SIMULATION OF ELECTRON SWARMS . . . 3269 kinetic energy (0.038 eV) along the forward direction, electrons are accelerated by the electric field and undergo collisions with the background gas at the same time. Then, after a certain interelectrode distance, where an equilibrium between collisions and electric field can be (or not) reached depending on discharge conditions, electrons arrive at the anode (assumed to be completely absorbing), which breaks down the eventual equilibrium phase. For E/N =31 Td [Fig. 4(a)], after a short nonequilibri-um distance near the cathode, W and ( e } reach their 0.20 > 0.15 0.08 0.20 FIG. 3. Zero-field electron mean energy in He for p = 1 Torr, 7* = 300 K, and an initial electron beam energy £0 —0.1 eV distributed along forward direction: (..+..+...) including thermal motion of targets and (..A. . .A. . .) target assumed to be at rest. 3270 M. YOUSFI, A. HENNAD, AND A. ALKAA 49 DL//x and DT/fi data, calculated from two different sets of HzO cross sections (Yousfi and co-workers [6] and Hayashi [7]), are compared to measurements (Ref. [12] for DL/fi and Ref. [13] for DT/^i). It should be known that the main differences between the two cross section sets are situated in the low kinetic energy range: 0.08 0.06 0.04 h- ^ 0.02 - CO (y (a) _i_[_ 0.02 0.04 0.06 0.08 0.10 z (cm) 0.5 FIG. 4. (a) Variation between cathode (z=0) and anode (z=0.1 cm) of electron mean energy (- and • • • •) and drift velocity (---and -.-) in H20 for p =1 Torr, 7 = 300 K, V„ = \ V, and an initial low energy electron beam (e0=0.038 eV) emitted at the cathode along the forward direction with and without including superelastic collisions: with (--), (---) and without (• ■ • •), (-.-). (b) Variation between cathode (z =0) and anode (z =0.5 cm) of electron mean energy (- and ...) and drift velocity (---and -.-) in H20 for p = 1 Torr, T = 300 K, Va = 5V, and an initial low energy beam e0 = 0.038 eV) emitted at the cathode along the forward direction with and without including superelastic collisions: with (-), (---) and without (• • • •),(-.-). 10 10 10" :/n (id) 10 FIG. 5. Hydrodynamic values of longitudinal eDL/fi ( + +, —----; and -) and transverse eDr/fi (AA; O O; - ■ ■ ; and —--) characteristic energies in H20 for r = 294 K. Symbols: measurements ( + , Wilson et al. [12]; A and O, Elford and co-workers [13]). Lines: Monte Carlo simulations using Yousfi and co-workers [6] cross sections (• • • • and---) and Hayashi [7] cross sections (• • • - and---) (eDT/(i data are multiplied by 10). Hayashi's set does not include rotational excitation cross sections (i.e., for energy range from about 10"3 eV up to 10"1 eV), which are correctly considered in the set of Yousfi and co-workers (see also Ness and Robson [14]). In fact, at the low energy range, classical crossed beam experiments for cross sections measurements are not reliable enough due mainly to the difficulty of maintaining low energy beams. So it is usual in that case to confirm the calculated cross sections (using the Born approximation for H20 rotational cross sections) by comparing calculated and measured swarm parameters. It is thus easy to observe in Fig. 5 the good agreement between measurements and Monte Carlo calculations using H20 set of cross sections taken from Yousfi and co-workers [6] and the rather great sensitivity of electron swarm parameters (about 50% of deviation) due to the difference between cross sections at low kinetic energy dominated by rotational excitation, superelastic collisions, and also energy exchange and thermal motion of molecules which are properly considered in present Monte Carlo method. C. Low E /N swarm parameters in strongly electronegative gases It is known that the electron distribution function and swarm parameters are calculated from the Monte Carlo method by following seed electrons from initial conditions until they vanish either by collisions (e.g., attachment) or by passing beyond the limits of the simulation domain (arrival at anode or reaching maximum time allowed for simulation, etc.). 49 MONTE CARLO SIMULATION OF ELECTRON SWARMS ... 3271 At low E/N values in gases such as SF6 having high attachment cross sections at low energy (for reaction e+SF6—>SFJ~), most seed electrons, after a few free nights, can be attached. In this case, the classical Monte Carlo simulation is not able to calculate hydrodynamic electron swarm parameters with enough precision. This corresponds to the case where electron current measured, for example, using the classical time resolved experiment (see, e.g. [18]), which is too low to be accurate. At very E/N values, the Monte Carlo simulation can be completely stopped as all seed electrons are attached to SF6. Under the same conditions, there is no arrival of electrons at the anode in swarm experiments. In that case, it is known that swarm experiments are undertaken not in pure SF6 but in gas mixture including a buffer gas (e.g., N2) and a small admixture of SF6. In this section, an improved Monte Carlo method is proposed (see Yousfi and Hennad [15]) in order to obtain swarm parameters with a better accuracy in the first case (low E/N values). It is based on an additional fictitious ionization process with a constant ionization frequency ufic ion, which artificially increases the number of simulated electrons. Therefore, the total collision frequency utot defined by relation (2) becomes vtot = v +vnun + uflCiion = const. The collision probability of this new fictitious process is />coijflc,ion = ufiCiion/i;tot. Then, for every ficti- tious ionization collision a new electron is created. In that case, the primary (or scattered) and new (or ejected) electrons are not deflected during this fictitious ionization and both keep the velocity of the incident electron. As all electrons (primary or secondary created either by real ionization or fictitious ionization) are treated, the electron distribution, density, and swarm parameters thus obtained are those of the fictitious gas (including the additional ionization process). Therefore, the problem is to calculate the electron swarm data for the real gas. In the case of spatially infinite medium, it is easy to establish a simple relation between the electron distribution function ff(v,t) of the fictitious gas and f(v,t) of the real gas: f(v,t)=ff(v,t)e (13) A similar relation can be written between electron number density itf(t) of the fictitious gas and n(t) of the real gas: n(t) = nf(t)e (14) Therefore, Monte Carlo calculations will be performed in the fictitious gas including the additional ionization process in order to overcome the problem of the extensive vanishing of seed electrons due to the strong attachment processes. This leads to the calculation of the electron distribution function ff and the density rij of the fictitious gas. Then the distribution function / and the number density n of the real gas will be deduced using relations (13) and (14). In the following this improved Monte Carlo method is applied to SF6 at low E/N values. Collision cross sections of SF6 are taken from Phelps and Van Braunt [16] for excitation processes and from Yousfi [17] for the other collision processes (elastic momentum transfer, attach- ment, vibration, and ionization). Monte Carlo calculations are performed in the case of an unbounded system (i.e., 3/3r=0) for which relations (13) and (14) are valid and which corresponds (when the memory of initial conditions is lost) to the simulation of classical time of flight or time resolved swarm experiments (see, e.g., [8]). Figure 6 shows the electron drift velocity and the mean energy in SF6 for E/N = 10 Td [Fig. 6(a)] and 2 Td [Fig. 6(b)] with and without using fictitious ionization. The chosen initial energy distribution is monokinetic along > CO o J. > (D w \ E X 1 1 1 1 1 J ! 1 1 1 J 1 1 1 1 J 1 1 1 1 J 1 1 1 1 (a) W I I I I 1 I I I I -1 I I 1 I I 'l I 10 20 30 40 50 t Cns) id \ 1 £ co -1 I I I I I I I I [ I I I I I I I I I I I I \. I. 1 (ol I I ... I 10 20 30 t (ns) 40 50 FIG. 6. (a) Electron mean energy (---) and drift velocity (• • • ■) in SF6 with (-) and without (---and • • • •) fictitious ionization: E/N = \0 Td, 5000 seed electrons, initial electron beam with 2 eV emitted along the forward direction and ufic jon = 2X 10s s_1. (b) Electron mean energy (---) and drift velocity (• • • •) in SF6 with (-) and without (--- and - • • •) fictitious ionization: E/N =2 Td, 10000 seed electrons, initial electron beam with 1 eV emitted along the forward direction and UjCion = 3X 108 s-1. 3272 M. YOUSFI, A. HENNAD, AND A. ALKAA 49 the forward direction and the number of seed electrons considered for the simulation is 5000 for 10 Td and 10000 for 2 Td. Starting from their initial values, energy and drift velocity relax towards their equilibrium values in both cases (i.e., with and without using fictitious ionization). However, the statistical fluctuations are far from negligible when the fictitious ionization is not used to improve Monte Carlo results. It is to be noted that, due to the large fluctuations, the drift velocity can be- 10 o c \ +> £ C 10 10 10 10 10 10 10 10 10 10 10 - 1 ! 1 1 1 1 ! j 1 1 1 1 = r z V\ = - — E E T HIM 1 1 1 I I t 1 I i I I 1 I I I I 10 20 30 t (ns ) 40 50 10 10 10 10 (b) 10 20 30 t (ns) 40 50 FIG. 7. (a) Reduced electron number density n(t)/n(0) in SF6 with (-) and without (—) fictitious ionization and also nf(t)/n(0) (• • ■ •): E/N = \0 Td, 5000 seed electrons, initial electron beam with 2 eV emitted along the forward direction, and UfiC,i0n = 2X 10s s-1. (b) Reduced electron number density n(t)/n(0) in SF6 with (-) and without (—) fictitious ionization and also nf(t)/n(Q) (• ■ • •): E/N = 2 Td, 10000 seed electrons, initial electron beam with 1 eV emitted along the for- come negative for a lower E /N value (2 Td) while the improved Monte Carlo method avoids this problem [see Fig. 6(b)]. As expected, statistical fluctuations are more pronounced in the case of the drift velocity, which is a statistical mean only of component v2 of the velocity along the z axis [see relation (lib)], contrary to the mean energy value, which is an average of the square of the three components vx, vy, and vz [see relation (11(a)]. Thus, for the same number of collisions, the electron mean energy is necessarily more accurate than the drift velocity. Furthermore, the choice of the collision frequency Dflcion of the fictitious ionization depends on E/N values. As the E/N value decreases, the attachment efficiency in SF6 increases (see, e.g., the time resolved measurements of the attachment coefficients in SF6 of As-chwanden [18]), thus necessitating a larger value of v6cion in order to compensate for the vanishing of seed electrons by attachment. For example, the chosen values of v6c ion are 2X108 s"1 for E/N = 10 Td and 3X108 s_1 for E/N =2Td. Figures 7(a) and 7(b) show the electron number densities nit) and nfit) [see relation (14)] in SF6 including the additional fictitious ionization process for E/N = 10 and 2Td. The electron number density calculated without including fictitious ionization is also shown in these figures. First, the electron number densities calculated with and without using fictitious ionization are, as expected, in good agreement. Concerning the fluctuations, a similar behavior to the swarm parameters shown in Fig. 6 is observed. Indeed, the density calculated with the classical Monte Carlo method shows statistical fluctuations due to the decrease of electron number density. For a lower E/N \alue [Fig. 7(b)], the decrease of the number density n(t) is more pronounced due to the higher attachment efficiency. Figure 7 also shows the density «y(r) calculat- x ü c Q) 3 0~ Q) L 20 40 60 t (ns ) 80 100 ward direction, and = 3X108s" FIG. 8. Electron attachment frequency (vatt) (---) in SF6 with fictitious ionization wficio„ =2X 108 s~'(-): E/N = 10 Td, 5000 seed electrons, initial electron beam with 2 eV emitted along the forward direction. 49 MONTE CARLO SIMULATION OF ELECTRON SWARMS 3273 ed with the improved Monte Carlo method using fictitious ionization. In fact, the density rij-(t) varies following the relation «/(f) = n(0)e)', where is the macroscopic attachment frequency of real gas (SF6), already defined in relation (11c). Therefore, for a given E /N value, it is better to choose a uficion value higher than the attachment frequency (vatt) to be sure to compensate for electron attachment efficiency. Such a choice is illustrated in Fig. 8, showing the time evolution of the attachment frequency (vatt) compared to fictitious ionization frequency in SF6 for E/N = 10 Td. In that case [i.e., uflc,ion-(vatt>(r)>0], the density «y(f) increases as time evolves, as is shown in Fig. 9. The distribution function f(v,t) or the density n(t) [deduced from relation (13) or (14)] and other elec- tron swarm parameters are therefore obtained with better accuracy. In conclusion, it is to be noted that the present Monte Carlo method is well adapted for the low E/N cases since the dominant collision processes (elastic including energy exchange and thermal motion of gas, inelastic and su-perelastic) are properly taken into account. This method is also adapted for the cases of strongly electronegative gases such as SF6 by using an additional fictitious ionization process with constant collision frequency. ACKNOWLEDGMENTS The Laboratoire des Decharges Dans Les Gaz is "Unite de Recherche Associee du Centre National de la Recherche Scientifique No. 277." [1] Y. M. Li, L. C. Pitchford, and T. J. Moratz, Appl. Phys. Lett. 54, 1403 (1989). [2] H. R. Skullerud, J. Phys. D 1, 1567 (1968). [3] S. L. Lin and J. N. Bardsley, J. Chem. Phys. 66, 435 (1977); I. D. Reid, Aust. J. Phys. 32, 231 (1979); J. P. Boeuf and E. Marode, J. Phys. D 15, 2169 (1982); M. J. Brennan, IEEE Trans. Plasma Sei. 19, 256 (1991); M. Yousfi, A. Himoudi, and A. Gaouar, Phys. Rev. A 46, 7889 (1992). [4] B. M. Penetrante, J. N. Bardsley, and L. C. Pitchford, J. Phys. D 18, 1087 (1985). [5] M. Yousfi, G. Zissis, A. Alkaa, and J. J. Damelincourt, Phys. Rev. A 42, 978(1990). [6] G. Friedriech and M. Yousfi, in XXth International Conference on Phenomena in Ionized Gases, Pisa, 1991, edited by V. Palleschi and M. Vaselli (Institute of Atomic and Molecular Physics, CNR, Pisa, 1991); M. Yousfi, A. Poinsignon, and A. Hamani, in Non-Thermal Plasma Technique for Pollution Control, Vol. 34 of NATO Advanced Study Institute, Series G: Ecological Sciences, edited by B. M. Penetrante and S. E. Schultheis (Springer-Verlag, New York, 1993), pp. 299-329. [7] M. Hayashi, Swarm Studies an Inelastic Electron-Molecule Collisions, edited by L. C. Pitchford (Springer-Verlag, New York, 1985), pp. 167-187. [8] M. Yousfi, P. Segur, and T. Vassiliadis, J. Phys. D 18, 359 (1985). [9] L. G. H. Huxley and R. W. Crompton, The Diffusion and Drift of Electrons in Gases (Wiley Interscience, New York, 1974). [10] M. Yousfi and A. Chatwiti, J. Phys. D 20, 457 (1987). [11] E. E. Kunhardt, J. Wu, and B. Penetrante, Phys. Rev. A 37, 1654(1988). [12] J. F. Wilson, F. J. Davis, D. R. Nelson, R. N. Compton, and O. H. Crowford, J. Chem. Phys. 62,4204 (1975). [13] R. W. Crompton, M. T. Elford, and R. L. Jory, Aust. J. Phys. 20, 369 (1967); M. T. Elford, International Conference on Phenomena in Ionized Gases, edited by W. T. Williams (University College of Swansea, Swansea, 1987), Vol. 1, pp. 130-131. [14] K. F. Ness and R. E. Robson, Phys. Rev. A 38, 1446 (1988) . [15] M. Yousfi and A. Hennad, in Gaseous Electronic Conference, Montreal 1993, edited by M. Moisan (University of Montreal, Montreal, 1993). [16] A. V. Phelps and R. J. Van Bräunt, J. Appl. Phys. 64, 4269 (1989) . [17] M. Yousfi, these de Doctorat d'Etat, Universite Paul Saba- tier de Toulouse, 1986. [18] Th. Aschwanden, in Gaseous Dielectrics IV , edited by L. G. Christophorou (Pergamon, New York, 1984), pp. 24-33.