A hierarchy of diffusion models for partially ionized plasmas. ∗

Partially ionized plasmas corresponding to different ionization degrees are derived and connected one with each other by the diffusion approximation methodology. These plasmas are the following electrical discharges: a thermal arc discharge, glow discharges in local thermodynamic equilibrium -LTE- and in non-LTE, and a non-LTE glow discharge interacting with an electron beam (or flow).


Introduction
Quasi-neutral gases of charged and possibly neutral particles which exhibit collective behaviour can either be natural or manmade. Well known examples of natural plasmas are solar corona, solar wind, nebula, lightning strokes, and aurora Borealis. Manmade plasmas have existed for almost two centuries and developments useful for practical applications are still in progress.
In this study we investigate cold plasmas, that means partially ionized with electron energies less than a few hundred of electron volts. The most common way to form and maintain a cold plasma is, up to now, the electric discharge in a gas, or gas discharge. Various electric discharges exist and produce plasma characterized by different parameters such as luminosity, electron energy, and ionization degree. Most manmade cold plasmas have electron densities lying within the range of 10 6 to 10 18 electrons/cm 3 , and electron energies from 0.1 to 20 eV. They can thus be used for a broad variety of applications. A simple way of making discharges (developed by Townsend) is to introduce a gas into a glass tube ended by two planar electrodes distant by l and apply a voltage V [18]. Consider for instance a DC discharge and a low pressure neon gas (1 Torr) with l = 50 cm [12], [19]. Varying the applied voltage V , three types of discharges may be observed and classified as function of the measured current I flowing through it: the dark discharge (where I < 10 −5 A), the glow discharge (10 −5 ≤ I ≤ 1 A), and the arc discharge ( I ≥ 0.1 A). In addition, each type of discharge is subdivided into regimes. The dark discharge (invisible to the eye except for corona) covers the background ionization regime (10 −10 ≤ I < 10 −9 A) and the saturation regime (I ≈ 10 −9 A) where ionization is only produced by surrounding sources such as cosmic rays. It also covers the Townsend regime (10 −9 < I ≤ 10 −5 A). There, electrons in the discharge start being enough accelerated by the applied electric field for initiating ionization reactions and producing secondary electrons, but not enough to self-sustain the discharge. As the applied electric field further increases, the secondary electrons may also in turn initiate ionization, leading to an avalanche. The Townsend regime includes the unipolar corona discharge (where visible radiations start being emitted) and sparking or electrical breakdown. Corona discharges can indeed be either initiated from electrodes with sharp geometries or from asperities creating a localized electric field larger than the breakdown voltage. Discharges get self-sustained when their voltage reaches the breakdown voltage. Then, the cathode layer generates enough electrons to balance the plasma current in the positive column. This delimits the transition from dark to glow (luminous plasma) discharge. The glow discharge includes the normal glow regime (10 −5 ≤ I ≤ 10 −2 A) characterized by almost independent current-voltage and the abnormal glow regime (10 −2 ≤ I ≤ 1 A) where the current increases with the voltage up to the glow-arc transition zone (0.1 ≤ I ≤ 1 A ). Dark and glow discharges are both non-thermal with gas (that is ions and neutral) energy much lower than the electron energy. Moreover electrons do not always sufficiently interact among themselves to achieve kinetic equilibrium, so that their energy distribution function is often non-Maxwellian. The arc discharge also has a non-thermal regime, which is associated with a falling voltage (1 ≤ I ≤ 10 A) and, as in this example, a low pressure (10 −3 -100 Torr). In addition, and contrary to the previous discharges, it posseses a thermal regime which is associated with an increasing voltage (I ≥ 10 A) and a high pressure (0.1-100 atm). There, all the species (including electrons) are usually in local or partial local thermodynamical equilibrium, and temperatures are of same order.
These three types of electric discharges are involved in various domains such as: wave absorption for stealthness [17], energy transport [14], industrial products and manufacturing [19], [20]. As an illustration, corona is the operating principle of electrostatic precipitators, xerography, surface modification of polymer films (wettability, adhesivity), antistatic applications for photographic films and plastic sheets, or plasma chemistry (for producing ozone for instance). Many applications are also based on high voltage plasma discharges, that is above the electrical breakdown. Glow discharge plasma is the operating principle of lightning devices such as neon and fluorescent lights, plasma screens, plasma chemical reactors. Glow may also be used for depositing thin films, as active media for gas lasers, and for forming charged particle beams allowing doping micro-electronics components for instance. Arc discharges are often coupled with a gas flow to form plasma jets with temperatures above the melting (and/or vaporization) point of many metals and ceramics; they are thus commonly used for material processing [2]. An important application field is in metallurgy to cut, melt or weld. Another concerns plasma torches used to spray protective coatings (for a better resistance to wear, corrosion, oxidation, thermal fluxes, also for electric or electronic purposes and bio-compatibility) or to vitrify toxic wastes for instance. In addition, these types of discharges may be associated to take advantage of their respective properties. A typical example is the gliding (or auto-oscillating periodic) arc that combines non-thermal arc properties (promoting selective chemical reactions) to thermal arc properties (high power level), [15]. Gliding arcs were first applied in the beginning of the 1900s for producing nitrogen-based fertilizers. Nowadays they are used for many chemical applications that cannot be made with conventional means (combustion): environmental control to reduce emission of indutrial exhausts such as chlorofluorocarbons for instance.
The cold plasmas investigated in this study are the typical non-thermal and thermal discharges, that is the glow and arc discharges. The ionization mechanism of a plasma glow is mostly provided by direct electron impact with non-excited atoms and molecules. Conventional values of glow discharge parameters are: high voltage between electrodes (100-1000 V) to promote enough secondary electron emissions, low electrode current (10 −4 -0.5 A), a power level of around 100 W. Concerning the plasma core (or positive column), the electron energy ranges from 1 to 3 eV, the gas temperature is close to room temperature, the ionization degree is weak with an electron density ranging from 10 9 to 10 11 electrons per cm 3 , and the gas pressure is low (0.03 -30 Torr). In the vicinity of atmospheric pressure, the ionization mechanism of a thermal plasma arc is mostly provided by direct electron impact with preliminary excited atoms and molecules. At larger pressure radiative ionization also gets important. Conventional values of thermal arc discharge parameters are: a lower voltage between electrodes (10-100V) leading to Joule heating of the gas, a larger arc current (30 A -30k A), and a power level per unit length larger than 1kW/cm. In the positive column electron and gas energies are of same order (1-10 eV), the ionization degree is moderate with an electron density within the range of 10 15 to 10 19 electrons per cm 3 , while the gas pressure is often atmospheric or even larger (0.1-100 atm).
We here focus on plasmas defined according to Langmuirs description: the region of a discharge not influenced by walls and/or electrodes. In other words, we study the plasma core, or positive column. The transition zone (sheath) between the plasma and its boundaries is not considered here. It screens electrically the plasma from the influence of its surrounding and has properties that differ from the plasma properties. Its modelling still raises difficulties. There exist various models developed for low potential fall and very narrow sheath as in arc discharges, also for larger potential fall and thick multi-layer sheath as in glow discharges. Further details are available in [12], [19].
As mentioned earlier, a gas discharge provides a mean to produce either ion or electron beams. The beam particles may initially have a rather low drift velocity and broad thermal velocity distribution. External electromagnetic fields are thus often used for accelerating them, narrowing their thermal velocity distribution, also for deflecting, focussing, or keeping parallel the beam in order to transport and utilize it effectively. Charged particle beams were first used for atomic and nuclear physics, and are now also applied to plasma diagnostics, space propulsion applications, film deposition [21] and ion implantation for microelectronics, to perform precision electron beam welding, rapid cutting of thermosetting plastics, cross-linking of thermoplastics to improve their physical properties, promote or increase plasma chemical activity, to invert the population of a gas laser and give rise to light amplification (from the soft X-ray region to the far infrared) [13], to control thermonuclear reactions via plasma heating [5], to process of surface treatment and depollution of high-volume exhaust streams [13], to support externally nonself-sustained discharges, or to study stellar plasma. Some of these listed applications combine an electron beam interacting with a plasma. We also investigate that case, considering a glow plasma that interacts with an electron beam characterized by a low drift velocity and a broad thermal velocity distribution. In the sequel, we will use the term electron flow rather than beam, to avoid any confusion with focussed and mono-energetic beams.
In this study, we start from the kinetic scale to derive macroscopic hydrodynamic/diffusion limits suited to the modelling of the plasma column of: an atmospheric thermal arc discharge (denoted by case 1), a glow discharge (case 2), and a glow discharge interacting with an electron flow (case 3). We account for impact ionization and recombination, and neglect radiative ionization and recombination (that get significant for high pressure thermal discharges). We thus investigate partially ionized plasma whose electrons, ions and neutral molecules are subject to elastic binary collisions as well as impact ionization and its reverse recombination reaction. The activation energy ∆ of ionization reactions is supposed to be constant and given by the impacting electron. We will see that the coupling between electrons and heavy species plays a major role.
Let us recall that the derivation of hydrodynamic/diffusion limits for a binary plasma gas mixture can be found in [9] for instance. The ternary gas mixture corresponding to a very weakly ionized plasma, such as a glow discharge where ionization occurs very seldom, is studied in [10]. A problem with dominant impact ionization and its reverse recombination is investigated in [7]- [8] within the frame of semiconductors (where all charged particles have masses of same order) and in [6] in the arc discharge context (where electron and ion masses differ by orders of magnitude). This paper is a followup of reference [6], where we give a more precise description of the model. We start this study introducing the kinetic model: a system of Boltzmann type transport equations governing the distribution functions of electrons, ions and neutral molecules. This system is coupled through collision operators that involve three collisional processes: i) elastic binary collisions where at least one particle is neutral (Boltzmann), ii) elastic binary collisions between charged particles (Fokker-Planck), and iii) inelastic collisions with impact ionization and its reverse recombination. This system is scaled in section 3, based on its two small parameters. The first parameter measures the relative smallness of the electron mass with respect to the neutral particles. The second parameter δ measures the ionization level of the plasma. The three different cases: themal arc, glow discharge, and glow interacting with an electron flow are then introduced. Sections 4 to 6 are devoted to some preparatory results: the conservation relations of the collision terms, the moment method, and the diffusion scaling, that will be used when investigating each of the three different cases. A model of thermal arc discharge is then derived in section 7. For readability purpose, the proofs of sections 6 and 7 are detailed in section 8. A model of glow discharge is then developed in section 9, and its interaction with an electron flow is detailled in section 10. The hierarchy between these various models, as well as the macroscopic limit linking two successive steps of this hierarchy are schematicaly summarized in section 11.

The kinetic model
Let us consider a mixture made of three species: electrons (e), ions (i) and neutrals (n), which interact all together, through various collisional processes which can be of elastic or inelastic type (taking thus into account ionization or recombination processes). Denoting by f α (α = e, i, n) the distribution function of the α species, the kinetic system modelling this mixture is given by [6]: where m α is the mass of the α species and F α denotes a force term; moreover, the notation (∂ t f n ) c stands for the collision terms which are given by: where the superscript ir stands for ionization-recombination, and α, β, γ = e, i, n with α = β = γ = α. Let us now describe more precisely these different collision terms, starting with the elastic ones.
Let us first consider binary elastic collisions between the two particles α and β. When one of these particles (or both) is neutral, collisions are governed by short range forces, significant only when the particles are in close proximity to each other. Then, the binary collisions are described by Boltzmann operators of the form: where α = n and β = e, i, n or α = e, i, n and β = n. In this expression, v α [resp. v β ] is the velocity of particle α [resp. β] before collision, and f α [resp. f β ] denotes: The post-collisional velocities v α and v β are defined from the pre-collisional velocities v α and v β by: where µ αβ = m α m β /(m α + m β ) is the reduced mass, and Ω ∈ S 2 + denotes a unit vector of part of the unit sphere S 2 of IR 3 defined by: S 2 + := Ω ∈ S 2 ; (v α − v β ) · Ω > 0 . The notations f α and f β stand for f α (t, x, v α ) and f β (t, x, v β ), respectively. The scattering cross section σ B αβ is a function of two variables: where E = µ αβ |v α − v β | 2 is the reduced kinetic energy and χ denotes the angle vα−v β |vα−v β | , Ω . While the former belongs to IR + , the latter lies within the range [0, 1].
Elastic collisions between two charged particles α and β are conversely governed by long range Coulomb interactions which act between each and every charged particle in the plasma. They are modelled by Fokker-Planck-Landau operators: where α, β = e, i and ∇ vα f α = (∇f α )(v α ), while S(w) denotes the matrix S(w) = Id − w⊗w |w| 2 , Id being the identity matrix. Here, the scattering cross section for grazing collisions σ F αβ only depends on the reduced kinetic energy: Radiative ionization and recombination are supposed to be negligible; the ionization process we consider is thus impact ionization. Its mechanism can be schematized by the following direct and reverse reactions: where e represents an electron, A + a single charged ion, and A the related neutral atom or molecule. σ d and σ r stand for the direct and reverse reaction cross sections. They are supposed to be positive. Applying the principle of detailed balance, we assume in the sequel that these cross sections are linked through where F 0 is a positive constant, which represents the efficiency of the dissociation with respect to the recombination. The ionization-recombination operators are then given by: We suppose that the activation energy of impact ionization reactions is given by the electron, and not by a heavy particle, so that the reverse reaction cross section writes and σ r = σ r (v e , v e ; v e ). The notations δ E and δ v hold for the energy and momentum conservation during the ionization-recombination process; more precisely, we have: where δ denotes the Dirac measure, and ∆ the ionization energy (which is a constant). In the same way, the notations δ E and δ v stand for: Notice that the factor 2 in Eq. (7-a) is a consequence of the indistinguishability of the two electrons in the right hand side of equations (5). This indistinguishability and the principle of detailed balance imply that Remark. In this study we do not consider the internal energy of atoms or molecules and the related ions. This is justified as far as glow discharges are concerned since then heavy particles are "cold". For arc discharges this assumption is a simplification since, as mentioned in introduction, the heavy species are pre-excited by Joule effect. The internal energy of atoms or molecules may thus contribute, jointly with the electron kinetic energy, to promote the ionization reactions.
The reference values of the problem are now introduced in order to scale the kinetic system (1).

Different scalings of the kinetic system
Let us first introduce the different small parameters involved in this study. First, ε denotes the parameter measuring the relative smallness of the electron mass with respect to the neutral particle: Case 1: For a thermal arc discharge the plasma column is in local thermal equilibrium (generally partial) so electrons, ions and neutral species have temperatures of same order of magnitude T 0 . Case 2 and 3: For a glow discharge, local thermal equilibrium is not satisfied. Electrons have energies of the order of the electronVolt (where 1 eV may be asociated to 11 065 K), while ions and neutral species have temperatures close to room temperature. So heavy and electron temperatures do not have the same order of magnitude; the difference is however rather small and gets almost negligible as far as thermal velocities are concerned (the square root of the characteristic temperatures ratio being less than one order of magnitude).

Case 3:
The electron flow retainned in this study is also supposed to have a temperature of order T 0 .
The characteristic velocities of all species (v α ) 0 are thus defined in all cases as the respective thermal velocities, i.e.
(v α ) 0 = kT 0 m α , with α = e, i, n, k being the Boltzmann constant. Consequently, these velocities only depend on the masses, and more precisely we have: Besides, we will choose x 0 = t 0 (v e ) 0 as reference length. The reference time t 0 is specified latter on, for different physical situtations. We assume that the mean densities of the charged particles (denoted by (ρ e ) 0 and (ρ i ) 0 ) are smaller than the typical density (ρ n ) 0 of the neutral particles; but, a priori, they can have different orders of magnitude. We thus denote by δ e and δ i the two small paremeters defined by: In the particular case (ρ e ) 0 = (ρ i ) 0 , we simply denote by δ the ratio and call it the ionization level. We get the following general orderings (compare with ([6]): and: The dimensionless kinetic equations then write: , , where we have set for simplicity: The scaled collision operators are now detailed. In the Boltzmann case, we have (note that the factor 1/ε just below is due to the fact that the integral term in the expression of Q ε ni is of order ε; we refer to [9] for details, and to Lemma A.1): and Concerning the Fokker-Planck-Landau case, the scaled collision operators read: We refer to Lemma A.1 and A.2 for a precise development, in terms of the small parameter ε, of these elastic inter species collision operators. Let us finally consider the inelastic collision operators. The scaled conservation equations (9) are written: where ∆ holds for the ionization energy scaled by the thermal energy kT 0 , so that these operators also depend on ε. Moreover, the factor F 0 has to be rescaled according to the relation: where k is the Boltzmann constant, and F 0 (which will be later simply denoted by F 0 ) is of order one. We refer to Lemma A.5 for a development of these operators in terms of ε and also to some properties (weak formulation, entropy ...) of their leading order terms.
All along the present study, we assume that the smallest time scale unit is the one related to the collisions with the neutrals, and more precisely that we have: We also investigate different orderings of the ionization processes (i.e. of the parameter τ ir ) and various values for the parameters δ e and δ i corresponding to different physical situations. More precisely, we consider three cases, which share in common that: Case 1: For a thermal arc discharge, the plasma column is free of space charge, i.e. quasineutral, thus δ i = δ e = δ, [6]. The ionization level δ lies within the range 10 −3 to 10 −1 , so that δ ε. This case thus corresponds to: Case 2: For a glow problem the plasma column is also free of space charge, i.e. quasi-neutral, thus δ i = δ e = δ. But the ionization level is lower than for an arc discharge; it lies within the range 10 −8 to 10 −5 , so that δ ≤ ε 2 . Here we will retain δ ε 2 which corresponds to (see [10]): Case 3: In the last problem, a glow discharge (as in case 2) interacts with an electron flow of relative numerical density δ e,beam = ε compared to the neutrals of the plasma glow. The electron and ion densities of the system (glow discharge + electron flow) have thus different orders of magnitude. More precisely we suppose that: So in all the cases under consideration here, we have: From now on, we suppose assumptions (15), (16) and (20) fulfilled (we will particularize each specific case later). So, the collision terms (13) have the following orderings (all the parameters involved in the equations below are small parameters): Within this framework, we remark that impact ionization can be a leading order collisional process. More precisely, in the first and third cases (i.e. under assumptions (17) or (19)), it is actually a dominant collision term for the ions. As a direct consequence, this will lead to a generalized Saha law (55) for the corresponding equilibrium states (see Proposition 6.1 below). In order to derive a macroscopic model, we use, if possible, a classical moment method, that we briefly recall in paragraph 5. But for this method to work, we first need some conservation relations concerning the collision operators that we now state.

Conservation relations of the collision terms
Concerning the elastic collision operators, we have the following classical conservations [3] , [4] : (i) The intra species elastic collision operators conserve mass, momentum and energy, which writes: (ii) The inter species elastic collision operators conserve mass, i.e.: while they globaly conserve the impulse and energy of the considered coupled system of particles, which means that we have: (iii) The inelastic collision operators conserve the electric charge and the number of heavy particles, i.e.: (iv) Last, these ionization-recombination operators conserve the global momentum and energy, which gives:

The moment method
The aim in this paper is to derive (when it is possible) a fluid model for the mixture. This one involves the following scaled macroscopic quantities associated with each species α: the mean density ρ α , mean velocity u α and total mean energy W α , which are defined by: One can also split the total energy W α into the kinetic energy ( 1 2 n α |u α | 2 ) and the internal energy according to: A way to derive an evolution system for the above macroscopic quantities consists in multiplying the kinetic equation for the α species by 1, v α , |v α | 2 , and integrating with respect to the velocity variable v α . If we do this manipulation directly on the kinetic system (12), we get the following system (setting for simplicity ε α = 1, for α = e, ε α = ε/ √ 1 − ε 2 for α = i and ε α = ε for α = n): where we have introduced the pressure tensor P α and the heat flux vector Q α defined by: we have also set, for the source terms: Let us note that, as all the elastic collision operators conserve mass, the source term R α only depends on the ionization-recombination term Q ε α,ir (f e , f i , f n ); moreover, the source terms S α and U α do not depend on the intra-species collision operator Q α,α (f α , f α ) because this one also preserves momentum and energy. But system (36) is not closed, because the space derivatives cannot be a priori expressed in terms of the macroscopic variables (ρ α , u α , W α ). One thus need a "closure relation", in order to compute the pressure tensor and the heat flux, which consists in supposing that the distribution function f α has a particular form which can be expressed only with the macroscopic variables (ρ α , u α , W α ). This is the usually called the "Local Thermodynamical Equilibrium" (LTE) assumption: the distribution function is supposed to be close to a Maxwellian, i.e. a function of the following form (still in scaled variables): For such functions we have in fact: with here: and also: A way to justify this LTE assumption consists in viewing the system at large time and space scales, far larger than the typical time and space units related to the collisional processes. So, from now on, we write the kinetic system at far larger time and space scales, and more precisely at the "diffusion" scale (t → ε 2 t, x → εx) (this terminology will get more clear later).

The diffusion scaling
We start from the scaled system of kinetic equations (12) with collision terms given by (21).
Introducing the diffusion scaling of small parameter ε: t → ε 2 t , x → εx, we obtain: Let us however point out that the sources terms in (44) and (45) are in fact of order 1 ε because, on account of (20), Q ε i and Q ε n are both of order ε. Now, taking the moments of these equations with respect to the velocity variables leads to the following system (apart from the index ε, the notations are those defined in the previous paragraph): and, for the heavy species (i.e. α=i,n) by: but with R ε α , S ε α and U ε α of order ε, still on account of (20). We now look for the limit ε → 0 of this macroscopic system (46) -(47). To this aim, we have to perform the limit ε → 0 in the microscopic system (43)-(45), with collision terms given by (21). We first replace in (21) the elastic (and inelastic for neutrals) inter species collision operators by their expansions in terms of ε (according to Lemma A.1, A.2 and A.5). Then we expand the solutions in powers of ε according to: and insert these expansions in the system (43)-(45). Finally, we identify terms of equal powers of ε, starting with the lowest order terms, which are of order ε −2 for the electrons (in equation (43)) and ε −1 for the other species (i.e. in equations (44) (45)). We get: Concerning the ions, if δ e = ε (which corresponds to cases 1 (17) and 3 (19)), we get: and, if δ e = ε 2 (i.e. in the second case (18)), then: This allows to derive the equilibrium distribution functions f 0 n , f 0 i and f 0 e . More precisely, we first obtain (the proof is postponed to section 8): Proposition 6.1: The equilibrium distribution function of electrons f 0 e is an isotropic function. The equilibrium states for neutral particles (f 0 n ) and ions (f 0 i ) are Maxwellians characterized by the same mean velocity u and temperature T : Moreover, if δ e = ε (which corresponds to case 1 (17) or to case 3 (19)), we have: and, if f 0 e = 0, then the density ρ i is given in terms of the two other species by: where δ 0 E stands for |v e | 2 = |v e | 2 + |v e | 2 + 2∆. Remark. Let us point out that relation (55) is a generalization of the Saha law (62) exhibited in [6] in the first case (i.e. under assumption (17)).
As the heavy species share in common the same mean velocity u and temperature T , we can drop the two last equations in the fluid system (47) for one species, for example for α = i, and keep the complete system only for α = n. Moreover, if δ e = ε (i.e. in case 1 (17) or 3 (19)), ρ i is completely determined, in terms of the other species, by (55); so, in that case, we can drop the whole fluid system (47) for α = i. We also remark that, at this level, the electrons have not necessarely reached their local thermodynamical equilibrium: the order zero term f 0 e is only for the moment an isotropic function. Moreover, on account of the assumption (20) (and the resulting scaling (21) for the collision terms), it will not be possible to go further (i.e. to obtain a Maxwellian distribution) if we suppose for example δ i = ε 2 , which corresponds to the two last cases under consideration here. But, concerning the first one, it will be (see section 7 below). As a consequence, the macrosopic models will be completely different, especially concerning the electrons: we obtain an "energytransport" (ET) model in the first case (i.e. a diffusion system on the electronic density ρ e and temperature T e ), and a model, that we will denote by "SHE-FP" model ("SHE" for "Spherical Harmonic Expansion" -using a common terminology of semi-conductors theory-and "FP" for "Fokker-Planck" model -referring to plasma context-), in the last two cases (i.e. a diffusion equation on the isotropic distribution function f 0 e ). But in both cases, the electronic models depend on the heavy species through transport coefficients and source terms. Notice that these results are in agreement with the physical properties of glow and arc discharges (recalled in introduction): in glow discharges, and contrary to arc discharges, electrons do not always sufficiently interact among themselves to achieve kinetic equilibrium, so that their energy distribution function is often non-Maxwellian.
So, we now propose to examine separately each situation, combining a classical moment method with a Hilbert method, when necessary, for electrons. We start with the first case, which corresponds to the arc discharge problem.

Case 1: a thermal arc plasma
We recall that this case corresponds to a partially ionized plasma with assumptions (15), (16) and: δ i = δ e = δ = ε [6]. The scaled collision terms are here given by: Replacing in (56) the elastic (and inelastic for neutrals) inter species collision operators by their expansions in terms of ε (according to Lemma A.1, A.2 and A.5), we get the following kinetic system: Next, we expand the solutions in powers of ε according to (48) and identify terms of order ε −1 in the kinetic equation (57). We get, using the notations of Lemma A.5 for the development of the ionization-recombination operator Q e,ir : The solvability condition for this equation of unknown f 1 e specifies the equilibrium state for the electrons. In fact, thanks to some entropy properties satisfied by the operator Q 0 e,ir (see Lemma A.5 below), we get the following result: The distribution function f 0 e is a centered Maxwellian, i.e. of the following form: Remark. According to (55), the densities ρ e , ρ i and ρ n are now linked to the electron temperature T e by the following Saha law (see [6]): We now see how the informations concerning these equilibrium states translate into the macroscopic system. In order to do this, we use a classical Hilbert expansion of all the macroscopic quantities, in the following way: At the lowest order, we get: ρ 0 e = ρ e , u 0 e = 0, W 0 e = w 0 e = 3 2 ρ e T e , ρ 0 n = ρ n , u 0 n = u, W 0 n = 1 2 ρ n |u| 2 + w 0 n , w 0 Concerning the pressure tensor and the heat flux, the order zero terms are given by: with p α = ρ α T α and α = e, i, n. This gives the leading order terms in the left hand side of the fluid equations satisfied by the heavy species, i.e. in system (47). But concerning the electrons, the corresponding leading order terms (i.e. in the left hand side of system (46)) are also connected to the order one corrections u 1 e and Q 1 e , which are given in terms of f 1 e by: where φ 1,o e is given by: In this expression, α is the isotropic function defined by: This allows deriving the expected terms in (66).

Corollary 7.3:
We have: with: where the diffusion matrix is given by: As a consequence, the leading divergence term in the electronic energy equation of system (46) is given by: First, as the leading order velocity term for electrons is perfectly known (in terms of the other macroscopic parameters), we can drop the second equation in the system (46). It now remains to look at the limits, when ε → 0, of the right hand sides in systems (46),(47). For that, let us do an asymptotic expansion, with respect to ε, in these source terms, according to: Concerning the neutrals, we get, first on account of the scaling of Q ε n and secondly thanks to (50), that: Moreover, as Q 0 so that the source terms in the fluid system (47) are all equal to zero for α = n. As a consequence, we get the following fluid model for neutrals, which is completely independant of the two other species: Proposition 7.4: The density ρ n of neutral particles, their velocity u and their temperature T are governed by the following fluid system (t > 0, x ∈ IR 3 ) ∂ t ρ n + div(ρ n u) = 0, where W 0 n = 1 2 ρ n |u| 2 + 3 2 ρ n T . Let us now consider the electrons. On account of the scaling (56) of the ionization-recombination collision term in Q ε e , and the fact that Q 0 e,ir (f 0 e , f 0 i , f 0 n ) = 0, we first get: R 0 e = R 1 e = 0, so that: R ε e = ε 2 R 2 e + O(ε 3 ). In the same way, thanks to (49) and (60), we have: U ε e = ε 2 U 2 e + O(ε 3 ). It thus remains to compute R 2 e and U 2 e . Lemma 7.5: We have: and where the source term R e is defined by: with the notations φ 1 α = f 1 α /f 0 α , for α ∈ {i, n}, φ 1,e e = f 1,e e /f 0 e and δ 0 E = δ(|v e | 2 −|v e | 2 −|v e | 2 −2∆). The relaxation coefficient λ is defined by: Gathering all the above results and notations, we obtain: Lemma 7.6: The electronic density ρ e and the temperature T e satisfy the following ET model: Remark. The first equation in (83) looks like a "mass conservation equation"; in fact, let us point out that the right hand side R e is here unknown, because it depends (in particular) on the arbitrary isotropic part f 1,e e of f 1 e . Now, thanks to (140), the ion density ρ i is such that In fact, we can remark that the equality R 2 e = R 2 i also results from the charge conservation (30) (at the kinetic level), which gives here: R ε e = R ε i . Taking into account (84), we can write system (83) in an equivalent form, where the source term R e completely disappears. More precisely, we obtain the following global two temperatures fluid model: Theorem 7.7: The equilibrium states for the heavy species, defined by (53), are characterized by the same mean velocity u and temperature T , while the electronic distribution function f 0 e is the centered Maxwellian given by (61). Moreover, the ion density ρ i is given in terms of the two other densities ρ n and ρ e , and the electronic temperature T e , by the following Saha law: The neutral particles satisfy an Euler system (setting W 0 n = 1 2 ρ n |u| 2 + 3 2 ρ n T ): which is totally independant on the other species. Finally, the electronic macroscopic quantities (ρ e , T e ) satisfy the following modified ET model: where we have introduced, for simplicity: with u J and v J defined by (72) and: Remark. The first equation in (86) corresponds, at the macroscopic level, to conservation of the number of heavy particles; the two last ones to the global momentum and energy conservation laws. The first equation in (87) is linked to the charge conservation, while the last one reflects the balance-sheet of the electronic energy.

Proofs of parts 6 and 7
Proof of Proposition 6.1: The equilibrium distribution function f 0 n is given by equation (50); applying the classical Boltzmann theory, we obtain (53) for α = n. Referring to Lemma A.1 and A.3, equation (49) means that f 0 e is a (non negative) isotropic function, i.e. f 0 e (v e ) = f 0 e (|v e |). Let us now determine f 0 i , first under the assumption δ = ε. In that case, equation (51) reads: We look for f 0 i in the form f 0 i = M u,T φ 0 i , with unknown φ 0 i . Let us set for simplicity: Using Lemma A.5 (point (i)), we obtain: where a 1 and a denote the following non negative constants (thanks to (8), a 1 and a are in fact independent of v i ): with the notation δ 0 E = δ |v e | 2 − |v e | 2 − |v * e | 2 − 2∆ . Referring to the definition (131) of L in , the determination of f 0 i solution of (90) reduces then to the derivation of the positive function φ 0 i solution of: There are two cases. First, if f 0 e = 0, then a 1 = a = 0, and the above equation admits all the constant functions as solutions; we then deduce that: f 0 i = ρ i M u,T , with arbitrary density ρ i . In the opposite case, i.e. if f 0 e = 0, then a > 0 and the operator L in − aId (Id denotes the identity operator) is invertible, so that equation (93) admits a unique solution. Let us identify this solution. We first observe that, integrating this equation against the Maxwellian M u,T , we get: aρ i = a 1 ρ n , which means that ρ i is uniquely determined in terms of the two other species by relation (55). Secondly, equation (93) With the notations introduced in Lemma A.3, we get, by integration on any arbitrary sphere S W (W > 0): Now, as f 0 e is an even function of the velocity, we know from Lemma A.1 iii) that Q 1 en (f 0 e , f 0 n ) is an odd function of v e and, remarking that v e · ∇ x f 0 e + F e · ∇ ve f 0 e is also odd, the above relation then reduces to: Let us now consider the isotropic function H f e = log F −1 0 ρ −1 n ρ i f e ; we get in particular: But, thanks to the classical H-theorem applied to the operator Q ee , and to Lemma A.5 iii), we have separately: so that (96) implies that each one of these two terms has to be zero. Now from the first one, we classicaly deduce that f 0 e is a Maxwellian; and this Maxwellian is necessarely centered, i.e. of the form (61), because it is isotropic. Conversely, let us suppose (61); then, we have Q ee (f 0 e , f 0 e ) = 0 and also Q 0 e,ir (f 0 e , f 0 i , f 0 n ) = 0, on account of (55), which here becomes (62) (see Lemma A.5, point (iv)), so that (95) is fulfilled.
Proof of Lemma 7.5: Let us introduce for simplicity Q 2 e defined by (Q 2 e is the term of order ε 0 in the right hand side of equation (57)): We have: where we have introduced: From the definition of Q 0 en (see Lemma A.1), and the properties of q B e (self adjointness, kernel made of isotropic functions, see [9]), we have, for any f and g (and setting C g = ( R 3 g(v)dv)): In the same way, as Q 1 en (f 0 e , f 1 n ) = q B e (∇f 0 e ) R 3 v n f 1 n dv n , we also have: From [11], we get: Concerning the linearized Fokker Planck operator, we classicaly have: It remains to compute the contribution due to the inelastic collision term. Referring to the definition of LQ e,ir in Lemma A.6, and to (134), simple computations show that we have: where we have set for simplicity: dv n dv e dv e dv e . But thanks to (138), we have C = R e , with R e given by (82). Now, taking into account the conservation relations (140) for the remainder term, we get which concludes the proof.
The scaled collision terms are here given by: At the diffusion scale, we get: We remark in particular that, concerning the ions, Q i,ir is no longer a leading order term, and, for the electrons, Q e,ir is ε 2 smaller than Q e,n . Expanding each distribution function in terms of ε according to (48), we still have (49), (50) and (52), so that the equilibrium states f 0 e , f 0 i and f 0 n are always given by Proposition 6.1; but we do not have the generalized Saha law (55) anymore. This means in particular that, at the macroscopic level, we willl have to keep the first equation in the fluid system (47) for α = i. Moreover, Proposition 7.1 is not valid anymore, because f 0 e is only an isotropic function (and not a priori a Maxwellian). In fact, f 1 e solves the following equation: which right hand side is an odd function of the velocity variable, so that its integral over any arbitrary sphere S W , of fixed energy W = |v| 2 /2, is zero. Moreover, we have: where the operator∇ x is defined, for any isotropic function by: Thanks to (97), we can precisely compute f 1 e (like in [16]), and we obtain: Lemma 9.1: The order one correction f e 1 is given by: where the operator L is defined, for any isotropic function f , by: Concerning the neutral particles, the macroscopic system of Proposition 7.4 is still valid. Concerning the ions, we have, with simple computations: Lemma 9.2: The density ρ i of ions satisfies the following equation: where R i is given by: using the notations: a 1 = A 1 (f 0 e ) and a = A(f 0 e ), with A 1 and A defined in Lemma A.5. Remark. We first remark that the source term R i only depends on the equilibrium states f 0 e , f 0 i and f 0 n and not on the corrections of order one, like in the previous scaling (i.e. the arc discharge problem). Moreover, this term, which is not a priori equal to zero here, is thus a relaxation term on the "general" Saha law (55), which is not satisfied here anymore.
If we go on for electrons and identify terms of constant order (with respect to ε) in (124), we get, thanks to the fact that Q 0 en (f 0 e , f 2 n ) = 0: With the notations of Lemma A.3, we integrate this equation on any sphere S W (W > 0) and obtain: The source term S e is given by (using the oddness of Q 1 en (f 1,e e , f 0 n ), Q 1 en (f 0 e , f 1 n ) and property (102) ): and does not depend on the arbitrary isotropic part f 1,e e . We then get the following result: The isotropic function f 0 e satisfies the following SHE-FP model: where ,∇ x is defined by (110) and D(W ) is the following diffusion matrix: Id.
The electron source term S e is given by (116). Gathering all the above results and notations, we have in summary: Proposition 9.4: The equilibrium states are given by: and we obtain the following hydrodynamic/SHE-FP system: Remark 1. If we multiply the SHE-FP equation satisfied by the isotropic function f 0 e by 1 and W , and integrate it with respect to W , we get the following energy-transport (ET) model for the electrons: where the electronic temperature T e is defined by: 3ρ e T e = R 3 f 0 e (v)|v| 2 dv and where we have set:ū Moreover, we have: The source term R e is thus also a relaxation term on the generalized Saha law (55); the relation R e = R i still reflects the charge conservation. If moreover, we suppose that f 0 e is Maxwellian, i.e. f 0 e = ρ e M 0,Te (this can be obtained, for example, by rescaling Q ee by a factor 1/ε), then we get for the electrons the above energytransport model with in addition: so that the second term in S E e is a relaxation term on the temperatures. The closed fluid model we obtain is then the following hydrodynamic/ET system: • ∂ t ρ e + div(ρ e (u + u J )) = R e , where R e = R i , which is given by (114), is a relaxation term of the Saha law (85); moreover, S E e = U e − ∆R e , with U e given by (89). So, up to the Saha relaxation term R e , we recover for the electrons the same ET model than in Lemma 7.6.

Remark 2.
We can recover the complete hydrodynamic/ET system of Theorem 7.7, either by rescaling the source term R e in the above system (121) by a factor 1/ε, or directly from the hydrodynamic/SHE-FP system of Proposition 9.4, by rescaling Q e,ir (in the source term S e ) by the same factor 1/ε (thanks to point (iv) of Lemma A.5, f 0 e is then a Maxwellian).

Case 3: a glow plasma interacting with an electron flow
In the last example under consideration, we assume that the densities of the charged particles do not have the same order of magnitude and more precisely that we have: The collision terms have here the following orderings: We observe in particular that, concerning the ions, Q i,ir is now again a leading order term (like for the arc discharge problem). At the diffusion scale, the kinetic system (12) writes: Next, we expand the distribution functions in powers of ε. At the lowest order, we get exactly the same system (49),(50), (51) as for the arc discharge problem; as a consequence, the equilibrium states are still given by Proposition 6.1, and we have: f 0 e is isotropic, while f 0 i and f 0 n are both Maxwellians, with same mean velocity u and temperature T ; moreover, the ion density is given in terms of the two other species by the general Saha law (55), and we also have: and we obtain the following hydrodynamic/SHE-FP system: with W 0 n = 1 2 ρ n |u| 2 + 3 2 ρ n T • ρ i given by (55) ( i.e. generalized Saha law ) Let us remark that this new system can be recovered from the previous one (i.e. the hydrodynamic/SHE-FP system of Proposition 9.4) by rescaling the ion source term R i by a factor 1/ε. Remark 1. If we multiply the SHE-FP equation satisfied by the isotropic function f 0 e by 1 and W , and integrate it with respect to W , we get the following energy-transport (ET) model for the electrons: • ∂ t ρ e + div(ρ e (u +ū J )) = 0, withū J ,v J defined by (118) and where denotes the electronic temperature T e (i.e. 3ρ e T e = R 3 f 0 e (v)|v| 2 dv). Moreover, the mass source term is here equal to zero (because it is linked to the equilibrium states, and not to the order one corrections): this is linked to the charge conservation, at the macroscopic level. Finally, S E e is given by: − ρ e F e ] + 3ρ e ρ n λ +λT n , withλ andλ defined by (120). If moreover, we suppose that f 0 e is Maxwellian, i.e. f 0 e = ρ e M 0,Te (rescaling Q ee by a factor 1/ε), we obtain the following closed hydrodynamic/ET system: n ) + div[u W 0 n + ρ n T ] − ρ n u · F n = 0, with W 0 n = 1 2 ρ n |u| 2 + 3 2 ρ n T, • ρ i given by (85) (i.e. Saha law), • ∂ t ρ e + div(ρ e (u + u J )) = 0, with U e given by (89). So, up to the Saha relaxation term R e (which is here zero), we recover for the electrons the same ET model than in Lemma 7.6. Last, this fluid system (127) can also be recovered from (121) by rescaling the ion source term by a factor 1/ε.

Remark 2.
We can also recover the first hydrodynamic/ET system of Theorem 7.7 from the hydrodynamic/SHE-FP system of Proposition 10.1, by rescaling Q e,ir (in the source term S e ) by the factor 1/ε.

Conclusion
The hierarchy between the various models previously derived, as well as the macroscopic limit linking two successive steps of this hierarchy are now summarised Fig. 1. Each of the considered models appears in a box. For simplicity, we have not specified the macroscopic Euler system for neutrals. An arrow between two boxes indicates the macroscopic limit connecting two models. Along the arrow is written the specific collision mechanism (i.e. the leading order collisional process associated with a specific space and time scale) used in the corresponding macroscopic limit. For simplicity, we have simply denoted by "neutrals" the elastic collisions against neutrals for each species. Note that we can also replace Q e,ir by Q ee each time in this figure.
A similar hierarchy has been established in [1] within the frame of semiconductors using the diffusion approximation methodology. An important difference lies in the fact that charged particles have masses of same order of magnitude when addressing semiconductor problems, while they differ by at least two orders of magnitude concerning plasma problems. In addition, impact ionization was not included in the set of collisions accounted for in [1].
.  i) Let α, β = i, n and α = β. Then ii) Let α, β = e, n and α = β. Then The superscript s indicates that a tensor is symmetrized, and the notation A : B, where A and B are two matrices with respective entries A ij , B ij , denotes the contracted product: i,j A ij B ij . Moreover, the linear operator q B e is defined by: iii) For any f e , f n we have: iv) Mass conservation imply that: Q j en f e , f n (v e ) dv e = 0, ∀j ∈ IN.
In the Fokker-Planck-Landau case, we obtain: i) Let f α with α = e, i, be sufficiently regular functions. Then ii) Mass conservation imply that: We now examine some properties of the linear operators involved in the different steps of the Hilbert expansion. In the sequel, we denote by M uα,Tα the normalized (i.e. with mean density equal to 1) Maxwellian of mean velocity u α and temperature T α defined by: We also denote by L 2 M u,T the weighted Hilbert space defined by: First, concerning the electrons, the linear operator involved is the operator L en defined by: where ρ n = R 3 f 0 n (v n )dv n is the density of neutral particles. Let us first recall a result of Ref. [1]. Lemma A.3: (i) The operator L en is self-adjoint on the weighted Hilbert space L 2 M 0,T . The kernel of the operator L en is made of isotropic functions, i.e. functions φ = φ(v) such that φ(v) =φ(|v|). In particular, if ϕ is an odd function of the velocity variable, then the equation L en φ = ϕ has a unique odd solution φ 0 and any other solution φ writes φ = φ 0 +φ, whereφ is isotropic.
(ii) More generally, let us introduce the energy variable W (v) = |v| 2 /2, and the sphere S W = {v ∈ R 3 , W (v) = W }; we recall the co-area formula (dS W is the euclidian surface element on S W ). Then the equation L en φ = ϕ has a solution if and only if the right hand side satisfies the following orthogonality relation: For the ions, we define the linear operator L in by with Q 0 in given in Lemma A.1. We have Proof: The proof partly results from the following weak formulation We now turn investigating some properties of the ionization-recombination collision operator: in Lemma A.5 below, we first state a weak formulation, and an entropy inequality, for the dominating part of this operator (in terms of ε), while Lemma A.6 is devoted to the computation of its linearization.
For any α ∈ {e, i, n}, we simply denote by Q 0 α,ir the limit, when ε goes to zero, of the operator Q α,ir . Taking into account the scaled conservation equations (14), we first get: so that we deduce the following result, which is a generalization of Lemma 4.5 of [6]: Lemma A.5: Let us set: ρ n = R 3 f n (v n )dv n and ρ i = R 3 f i (v i )dv i . i) We have, for any α ∈ {e, i, n}: where the leading order operators Q 0 α,ir are still given by expression (7-a)-(7-c), but where the conservation equations δ v , δ E (14) have been replaced by their limit when ε goes to zero, i.e. by: Moreover, introducing the two following contants ii) For electrons, Q 0 e,ir f e , f i , f n only depends on the heavy species through their densities ρ n , ρ i , and we have the following weak formulation: iv) In particular, if f e is isotropic (i.e. f e (v e ) = f e (|v e |)) and such that Q 0 e,ir f e , f i , f n = 0, then we have: Conversely, if f e is given by (135), then we have, for any f n , f i with fixed mean ratio of the densities ρ n /ρ i : Q 0 e,ir f e , f i , f n = 0.
Let us now precise the expansion, in terms of ε, of each ionization-recombination operator around the equilibrium state. Lemma A.6: We suppose that: f i 0 = ρ i M u,T , f n 0 = ρ n M u,T and f e 0 = ρ e M 0,Te . Next, we expand the distribution functions in terms of ε by setting: Then, if (135) is satisfied, we have Q 0 α,ir f 0 e , f 0 i , f 0 n = 0, for α = e, n, i and: with: where δ 0 E still stands for |v e | 2 = |v e | 2 + |v e | 2 + 2∆. The remainder terms R α,ir 1 are linked to the order one corrections in the asymptotic expansions (in terms of ε) of the conservation equations (132); thanks to the following property: we have: and R 1 e,ir is an odd function of the velocity variable. With (135), we deduce the following conservations: e,ir (v e ) dv e = R 3 R 1 e,ir (v e ) |ve| 2 2 dv e = 0, Proof: Under the assumption (135), we deduce from the expansion (132) of the conservation equations that: Now, inserting the expansion (136) in each ionization-recombination collision term, we obtain, after some easy computations, the expected results. We just precise that the order one remainder terms R 1 i,ir and R 1 n,ir , which are given by: are in fact zero, thanks to relation (138).