NON EQUILIBRIUM IONIZATION IN MAGNETIZED TWO-TEMPERATURE THERMAL PLASMA

. A thermal plasma is studied accounting for both impact ionization, and an electromagnetic ﬁeld. This plasma problem is modeled based on a system of Boltzmann type transport equations. Electron-neutral collisions are assumed to be much more frequently elastic than inelastic, to com- plete previous investigations of thermal plasma [4]-[6]. A viscous hydrody-namic/diﬀusion limit is derived in two stages doing an Hilbert expansion and using the Chapman-Enskog method. The resultant viscous ﬂuid model is characterized by two temperatures, and non equilibrium ionization. Its diﬀusion coeﬃcients depend on the magnetic ﬁeld, and can be computed explicitely.

(Communicated by Pierre Degond) Abstract. A thermal plasma is studied accounting for both impact ionization, and an electromagnetic field. This plasma problem is modeled based on a system of Boltzmann type transport equations. Electron-neutral collisions are assumed to be much more frequently elastic than inelastic, to complete previous investigations of thermal plasma [4]- [6]. A viscous hydrodynamic/diffusion limit is derived in two stages doing an Hilbert expansion and using the Chapman-Enskog method. The resultant viscous fluid model is characterized by two temperatures, and non equilibrium ionization. Its diffusion coefficients depend on the magnetic field, and can be computed explicitely. 1. Introduction. In this study we pursue the investigation of arc discharges in thermal regime treated in the references [4]- [6], and extend this investigation of thermal plasma to induction thermal plasma. A thermal plasma is characterized by dominant thermal effects coupled with fluid dynamic. When made by an arc discharge or an imposed magnetic field, its thermal and kinetic energies result from a transfer of electromagnetic energy. Examples of applications of thermal plasma concern i) material processing based on thermal energy transfer [10]: cutting, melting, welding, plasma spraying and vitrification, and ii) propulsion with electrothermal thrusters [15]. Material processing generally makes a direct use of the thermal energy, while electrothermal thrusters further convert it into directed kinetic energy.
Plasma with a lower ionization degree than a thermal plasma, or glow plasma, was studied by Degond et al. in [9], and by Choquet et al. in [6] accounting for impact ionization. On the other hand, the case of a fully ionized plasma was investigated by Degond et al. in [7], [8], and by Lucquin-Desreux in [18] for multicharged ions and in [17] in the presence of a magnetic field.

ISABELLE CHOQUET AND BRIGITTE LUCQUIN-DESREUX
The remaining part of this introduction is organized as follows: we first recall some properties of thermal plasma usefull for setting the problem, the scaling, and discussing the results. Next we do summarize specificities of previous studies on thermal plasma modelling [4]- [6], to underline the new aspects of the present study. Then we do present some of the main results of this paper. Based on reference models available in the literature, we discuss the inviscid fluid model obtained in this study. Finaly we outline the steps of the model derivation, which are further developed in the rest of this paper.
1.1. Physical properties of thermal plasma. A thermal plasma can be made in different ways. A traditional way consists in forming an arc discharge by applying a potential difference between electrodes separated by a neutral gas, as in [11], [14], [27], and [30]. More details about arc discharges and comparisons with other electric discharges can be found in [5]. The structure of a thermal plasma arc is not homogeneous. For this reason it is usually modeled by doing a splitting into different sub-regions: the main core or plasma column, and the thin electrode layers (made of an ionization zone or presheath and a space charge layer or sheath). A further detailed description of these sub-regions can be found in [4]. A thermal plasma can also be made applying instead a magnetic field to a neutral gas, to form a so-called induction thermal plasma as in [31], [25]. It can as well result from the interaction of a glow discharge (cold, luminous and self-sustained electric arc) with an electron beam for instance, as modeled in [6]. Or it can result from a very strong compression wave, as in hypersonic flow. In the sequel we focus on the first cases, cases in which an external electric or magnetic field is imposed. From now on we will simply call them thermal plasma.
Two important physical characteristics specific to thermal plasma (except in an arc sheath), and useful for the present study are [22], [28] : i) a weak ionization degree δ of the order of 10 −3 to 10 −1 , and ii) a moderate electromagnetic field. The ionization degree δ is defined as the ratio of the electron number density ρ e to the neutral number density ρ n : δ = ρ e /ρ n . The electromagnetic field comprises the imposed electric or magnetic field and the resultant induced field. This last one can be obtained based on the Maxwell equations. A moderate electromagnetic field means here a field large enough to promote impact ionization but small enough to avoid runaway electrons.
In this context neutrals, ions and electrons coexist. The numerical density of neutrals is much larger than the numerical density of any charged particle, so that collisions with neutral particles are the most frequent. They govern the velocity of neutral particles, which is random.
Charged particles are also subject to collisions (mostly with neutrals) that govern the random component of their velocity. They do have an additional velocity component: the drift velocity induced by the electromagnetic field. The proportionality factor between drift velocity and electromagnetic field, or mobility, is inversely proportional to the mass of the charged particle. Because of a much lower mobility, the drift velocity of an ion is at least two orders of magnitude lower than the drift velocity of an electron. In the context of a moderate field, the drift velocity is then negligible for ions while significant for electrons. Furthermore, the characteristic drift velocity of electrons does not exceed their characteristic random velocity (no runaway). This last property dictates in the sequel the scaling of the electromagnetic forces.
If the thermal plasma is locally neutral, the electric current only results from the drift motion of charged particles. It is then mainly due to electron drift. The electron drift is damped by the electron collisions resulting in a loss of electron momentum along the electromagnetic field. Electron-electron and electron-ion collisions do not produce any significant resistance to electron drift. The former because it conserves the total current of the pair, and the later because it happens too rarely in a weakly ionized plasma. Electron-neutral collision is here the most frequent type of collision for electrons. The characteristic time for electron-neutral collision per electron, τ en , set to is thus less than the charateristic time for electron-electron collision per electron, τ ee . When elastic, electron-neutral collisions can allow transferring to a neutral the momentum and energy gained by an electron from the electromagnetic field along its free flight. Electron-neutral collisions can thereby significantly impede the electron drift along the field. This mechanism is also the basic principle of operation of thermal plasma used for propulsion and material processing: the plasma ability to transfer electromagnetic energy into thermal energy (via Joule heating). The thermal energy is distributed among the particles via various elastic and inelastic collision processes. The elastic ones include: i) collisions between charged particles, governed by long range Coulomb interaction (modeled by Fokker-Planck-Landau operator), and ii) binary collisions involving at least one neutral particle (modeled by Boltzmann operator). The inelastic collision retained here is the impact ionization and its reverse recombination, where the ionization energy is provided by an electron. We thus consider a plasma pressure of the order of the atmospheric pressure, so that radiative ionization and recombination are negligible compared with impact ionization and recombination. In addition we do the usual assumption of optically thin plasma, implying that the emitted radiation is not reabsorbed and radiative phenomena can be neglected compared with collision phenomena. Notice that in practice thermal plasmas are often optically thick to some wavelength [28].
The efficiency of energy transfer by elastic collision is proportional to the mass ratio of the colliding particles. This may lead to differences in thermal equilibrium for the various types of species present in the plasma. The energy transfer between heavy particles (neutrals, ions) is indeed much more efficient than the energy transfer between heavy particle and electron. As a result, the heavy particles can tend towards the same thermal equilibrium. In contrast, electrons may have a larger temperature (assuming that they do reach a thermal equilibrium) in regions where their frequency of elastic collision with heavy particles is low, such as in the presheath of an arc. In a similar way, in regions where the frequency of inelastic collision between electron and neutral is low, the ionization balance (also called Saha balance) is not reached. These two types of deviations from equilibrium (thermal and Saha) can decisively determine the thermal plasma behavior. They do increase when the electron density decreases, such as in the outer regions of the plasma. The density threshold below which thermal and Saha deviation appears is not strictly the same. The characteristic frequency of both elastic and inelastic collisions between electrons and neutrals are thus important parameters that dictate the order of magnitude of collision operators and in turn govern the behavior of the macroscopic limit.

1.2.
Context of this study. Thermal plasmas were previously studied in [4]- [6] in the presence of an electric field (but no magnetic field) and accounting for impact ionization. These studies share a common point with the present work: the small parameters defined by the ionisation degree δ (see i) above) and the mass ratio of electron to neutral ε = m e /m n are in all cases of same order, δ = ε. The developments done below thus involve only one small parameter, namely ε. The basic difference lies in the ordering of the elastic and inelastic collision frequencies between electrons and neutrals, ν el ∝ τ −1 en and ν inel ∝ τ −1 ir , respectively. • Thermal plasma with strong ionisation (ν inel = ν el ) is treated in [6]. The impact ionization is then a leading collisional process for both ions and electrons. As a main consequence, the equilibrium distribution function of electrons is a Maxwellian, ionization equilibrium is satisfied (Saha balance), and the diffusion coefficients of the fluid model derived for electrons depend on the ionization process. These coefficients are defined implicitely as functions of the first order correction of the electron distribution function. • Thermal plasma with moderate ionisation (ν inel = εν el ) is investigated in [4]- [5]. The impact ionization is then a leading collisional process for ions, but no longer for electrons. The equilibrium distribution function of electrons is still Maxwellian, ionization equilibrium is also satisfied (Saha balance), but the diffusion coefficients of the fluid model derived for electrons do not depend on the ionization process. These coefficients are again defined implicitely as functions of the first order correction of the electron distribution function. • The case of thermal plasma with weak ionization (ν inel = ε 2 ν el ) is investigated in the present study. In this frame, and using (1), we have thus where τ ir denotes the characteristic collision time per electron for inelastic electron-neutral collision.
1.3. Results of this study. Some of the main results of this paper are now summarized and discussed. In the frame of a thermal plasma with ν inel = ε 2 ν el , the impact ionization is no longer a leading order collisional process. The equilibrium states obtained for neutral particles and ions are Maxwellians characterized by the same mean velocity u and temperature T . As expected at the order ε 0 and for low concentrations of charged particles, electrons do not affect the heavy particle distributions, ions do not affect the neutrals, and the energy gained by ion drift along the electric field is negligible (a higher order development is needed to observe it). The resultant Euler system governing the density, momentum and energy of the neutral species is thus totally independant of the other species.
The equilibrium distribution obtained for the electrons is a centered Maxwellian at temperature T e . The electrons mean velocity thus differs from the mean velocity of the whole fluid. The partial local thermal equilibrium (T e = T ) expresses the poor efficiency of energy transfer in collisions between particles of disparate masses. During elastic light-heavy collisions, the velocity of the heavy particle is indeed almost unchanged, the velocity norm of the electron is also almost unchanged, only its velocity direction is changed. The energy relaxation occurs at the higher order in the developments. Ionization is not at equilibrium (no Saha balance) so that the density of electrons is governed by a continuity equation, instead of being defined from the electron temperature and the densities of the heavy particles via a generalized Saha law as in [4]- [6]. The electron density and temperature satisfy an energy-transport model in which electrons gain energy by drift along the electric field, and the diffusion matrix depends on the magnetic field B. Each submatrix of the diffusion matrix is anisotropic due to the presence of the magnetic field. It is moreover perfectly computable in terms of B, the neutral density, and the Boltzmann kernel.
The first order correction derived for the electron distribution function can be computed explicitely. This allows pushing further the developments to account for viscous effects, and deriving a macroscopic (hydrodynamic/diffusion) limit at a higher order than previously done in [4]- [6], using here the rigorous Chapman-Enskog theory. Notice that a particular attention to the scaling of the forces applied on the charge particles is then needed. Computing the corrective terms of higher order, we do obtain a viscous fluid model which is globally conservative, up to the second order. This system governs the number density of neutrals ρ n , ions ρ i , electrons ρ e , the momentum ρ n u and energy W n of neutrals, and the electron energy. It is given in scaled form in Theorem 5.1, section 5. The variables are now changed to the total density , the total charge χ, the electron density, the heavy species momentum u and energy W , and the electron energy to facilitate the discussion. The scaled system then reads, up to the second order ε 2 , The density of electrons does not enter the total density , as the fluid system (3)- (4) results from developments at the order ε. = α m α ρ α is thus the mass density of the mixture made of heavy species of mass m α and number density ρ α , and u is the plasma (also heavy particle) velocity. χ = Zρ i − ρ e is the charge number density of the plasma where Z denotes the ion charge number. From now on we will assume that Z = 1, to simplify the computations. j tot = χu − ρ e u J is the total current density (per unit charge), j ρe = ρ e (u + u J ) is the electron current density (per unit charge) with u J given in (27)- (28) and (32)-(34), and j ρ T is given in (40). The first term χu of the total current density, j tot , represents charge conduction by the plasma flow. The second term ρ e u J contains the following contributions: the electron drift induced by the electric field, and the diffusion current. The diffusion current is induced by both the ordinary diffusion of electrons, and the thermal diffusion (Soret effect) resulting from electron temperature gradients. Notice that induction current and Hall current are not observed at this order (ε). The diffusion coefficients, obtained from the explicit computation of the first order correction of the electron distribution function do depend on the magnetic field. They are given by (27)- (28) and (32)-(36). In regions where local electroneutrality is satisfied (such as in the plasma core of an arc) the charge conservation equation reduces to the usual current continuity law (div j tot = 0) and the diffusion coefficients to ambipolar diffusion coefficients. The source term in the continuity equation governing the electron number density is the rate of change of number density due to impact ionization and its reverse recombination reaction. This source termρ e = R e , is given in (25), (29). It represents a relaxation term towards a generalized Saha law that depends on the electron temperature and not on the heavy particle temperature. This Saha law turns out to be a dimensionless version of Eindhoven's formulation [12], as previously underlined in [4]. This result thus supports previous argumentations, such as in [12], in favour of Eindhoven's generalization of the Saha law. Other generalized formulations, given by Potapov [21] and Van de Sanden [29], can depend on both the electron and the heavy particle temperature. In that case they do not seem to be consistent with the entropy inequality, which is needed to derive the macrosopic limit.
The heavy species energy is W = 1 2 |u| 2 + 3 2 T , and σ(u) is the traceless rate of strain tensor where (σ ij = ∂ui ∂xj + ∂uj ∂xi − 2 3 divuδ ij ). µ n (T ) and κ n (T ) are the usual viscosity and thermal conductivity of neutrals ( [2], [3]). At this order (ε) the contribution of ions (and not only of electrons) to the viscosity and thermal conductivity of the plasma indeed turns out to be negligible.
The source terms of energy for the heavy species are at the order ε, and include the rate of drift energy gained by ions while moving along the electric field (first term) and the energy gained per unit time by neutrals during collisions with electrons (inẆ ch ). Notice that in regions where local electroneutrality is verified the energy gained by electrons and ions from the electric field do cancel each other. As expected,Ẇ ch also appears as a sink term in the equation governing the thermal energy of electrons.Ẇ ch consists of electron thermal diffusion (in the first term), energy flux due to electron concentration gradients (or Dufour effect, also in the first term), electron drift energy built up along the electric field (second term), induced magnetic energy (third term), and thermal relaxation due to cooling of electrons by neutrals (last term).

1.4.
Comparison with models used for simulating thermal plasma. Models used for simulating thermal plasma often assume thermodynamic and chemical equilibrium. A model in chemical non-equilibrium and thermodynamic equilibrium proposed by Tanaka et al. for a pulsed arc can be found in [26]. The discussion will focus here on models accounting for partial thermodynamic equilibrium. Examples of such models have been proposed by Selezneva et al. for induction plasma [25], by Ghorui et al. [11] to study oxygen-plasma arc cuting torches, by Gonzales et al. [14] to simulate a SF 6 arc plasma, by Trelles et al. [27] to investigate arc plasma torches used in plasma spraying for instance, and by Wendelstorf [30] to study gas tugnsten arc welding. An important difference compared with the present study is that most of these models (except in [30]) are obtained assuming local electroneutrality, implying that χ = 0. So the arc models above listed do apply to the plasma core and do not adress the presheath.
All these models generally assume that each species satisfy the ideal gas law. However, the closure equation is not yet specified in the fluid system (3)-(4). An ideal gas law or a closure equation for thermaly expansible plasma could either be used. For consistency, the comparison with existing models is now done considering also the ideal gas model, and in the particular case χ = 0 (except for [30]). Then the main differences are the following: -In the fluid system (3)-(4) derived here, ions do not contribute to the heavy species viscosity and thermal conductivity (simply as a result of the low ionization degree of a thermal plasma). This is a difference compared with models used in practice for simulating thermal plasma, where the viscosity and thermal conductivity are instead computed for the mixture of heavy species (so accounting for ions too, as clearly detailed in [30] for instance).
-A term of thermal diffusion of electrons is present in both the momentum and the energy equation for the heavy particles, in the system (3)-(4) derived here using the Chapman-Enskog method. This term is not accounted for in [25], [11], [14], [27], and [30].
-The magnetic induction is handled differently depending on the author. It is neglected in [14], and [30]. It is accounted for into the current in [11], and [27]. In the fluid model (3)-(4) derived here, the magnetic induction comes out at an intermediate order, through the term of energy exchangeẆ ch between electrons and heavy particles.
-The diffusion matrix derived in this study depends on the magnetic field. On the contrary, the diffusion matrix used in [25] for an induction thermal plasma, for the thermal arc plasma of [11] and [27] based on Murphy's model [19], and [14] based on Gleize et al. model [12], and in [30], do not account for the magnetic field. The magnetic field should be significant when applied as external field, such as for induction thermal plasma [31], [25]. When the external field is electric, the magnetic field is induced. In regions of the plasma where local electroneutrality is not verified (such as in the presheath of an arc) the non-zero conductive current induces a magnetic field that should also be significant and could in turns influence the diffusion coefficients. In regions where local electroneutrality is satisfied, such as in most of the plasma core of a thermal arc, the induced magnetic field is generally rather small, and it could thus be neglected in the diffusion matrix. However, this simplification may not hold if the electric potential presents large gradients.
Notice that there exist studies focussing on the modeling of transport properties for a two-temperature thermal plasma, such as [24] and [23] for instance. Compared with the present work, the form of the first order perturbation function is presumed in [24] and [23], rather than derived. The expressions proposed in [24] and [23] differ from each other: the electron temperature is scaled by the heavy species temperature in [24] and not in [23]. On the other hand the diffusion coefficients account for the magnetic field in [23] and not in [24].

1.5.
Outline. To obtain the fluid system (3)-(4) we start this study at the kinetic level, introducing in section 2 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 collision 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. The system involves two small parameters: the ionization level of the plasma, δ, and the relative smallness of the electron mass with respect to the neutral particles, ε 2 . Within the frame of thermal plasma these parameters are of same order [5]: δ ≈ ε, so that the scaling is done as function of ε alone. An inviscid fluid limit is then derived rigorously in section 5, doing a first order Hilbert expansion at the electron diffusion scale. For readability purpose, preparatory results concerning some properties of the collision operators, and conservation relations of the collision terms, as well as proofs of the inviscid fluid limit, are postponed to appendix A (in section 7), appendix B (in section 8), and section 6, respectively. The higher order (or viscous) fluid limit discussed above is derived rigorously in section 6, doing a Chapman Enskog expansion. The entropy of this system is also determined in view of a forthcoming numerical study.
2. The kinetic model. We consider a magnetized plasma made of electrons (e), ions (i) and neutrals (n), which interact all together, through elastic and inelastic collisions. The inelastic collision processes is related to ionization and its reverse recombination. Denoting by f α (α = e, i, n) the distribution function of the α species (of mass m α and charge q α ), the kinetic system modeling this mixture, [4], is given by where F α denotes the Lorentz force given in terms of the electric field E and the magnetic field B by F α = q α (E + v α ∧ B). We suppose that the charges are of the same order of magnitude, more precisely that q i = −Zq e = Ze, where e is the elementary charge of an electron and Z = 0(1). For simplificity, we will suppose from now on that Z = 1. 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 α = β = γ = α. We now recall the expressions of these different collision terms, starting with the elastic ones. When at least one of the two particles α or β is neutral, binary collisions are governed by short range forces and described by Boltzmann operators of the form In this expression, v α [resp. v β ] is the velocity of the particle α [resp. β] before collision, and 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 respectively. The scattering cross section σ B αβ is a function of the reduced kinetic energy E = µ αβ |v α − v β | 2 and of the angle Elastic collisions between two charged particles α and β are conversely governed by long range Coulomb interactions; they are modelled by Fokker-Planck-Landau operators, [16], Here, the scattering cross section for grazing collisions, σ F αβ , only depends on the reduced kinetic energy. S(w) denotes the matrix S(w) = I − w⊗w |w| 2 , I being the identity matrix. Radiative ionization and recombination are supposed to be negligible; the ionization process we consider is thus impact ionization. When Z = 1, 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, while σ d and σ r (which are positive) stand for the direct and reverse reaction cross sections. Applying the principle of detailed balance, we assume in the sequel that these cross sections are linked through σ d = F 0 σ r , 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 assume that the activation energy of impact ionization reactions is given by the electron, so that σ r = σ r (v e , v e ; v e ), and σ r = σ r (v e , v e ; v e ). The notations δ E and δ v hold for the energy and momentum conservation during the ionizationrecombination 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 indistinguishability of the two electrons in the right hand side of equations (6) and the principle of detailed balance imply that σ r = σ r . We now introduce the characteristic physical units of the problem, in order to scale the kinetic system (5).
3. Scaling of the kinetic system. We first introduce the small parameter that represents the mass ratio between the electron and the neutral particles. We assume that electrons, ions and neutral species have temperatures of same order of magnitude T 0 . The characteristic velocities are the respective thermal velocities given by k B being the Boltzmann constant. Consequently, these velocities only depend on the masses, and more precisely we have We choose x 0 = t 0 (v e ) 0 as reference length; the reference time t 0 is specified latter on. Concerning the force fields, we choose the electric field unit E 0 and the magnetic one B 0 such that the electron drift energy gained during a mean free path does not exceed the thermal energy: We assume that the mean numerical densities of the charged particles (ρ e ) 0 and (ρ i ) 0 are both ε smaller than the typical numerical density (ρ n ) 0 of the neutral particles, i.e.
Denoting by τ αβ [resp. Q αβ 0 ] the characteristic collision time [resp. characteristic elastic collision operator] of a particle of β species against a particle of α species, we have: Q αβ 0 = (f α ) 0 /τ αβ . We then get the following orderings (see [5] for further details): The dimensionless kinetic equations then write , where As mentioned in the introduction, we assume that the smallest time scale unit is the one related to the collisions with the neutrals. Using (2) we thus set As underlined in the introduction (in the sub-section "Context of this study"), these orderings are completely new compared with previous ones [4], [5], [6]. The collision terms (12) of the dimensionless kinetic equations (11) have then the following orderings: The scaled collision operators are now detailed. Note that the factor 1/ε just below is due to the fact that the integral term in the expression of Q ε ne is of order ε; we refer to [7] for details, and to Lemma 7.1. In the Boltzmann case, we thus have and , Ω) f n f n − f n f n dv n dΩ.

ISABELLE CHOQUET AND BRIGITTE LUCQUIN-DESREUX
The scaled Fokker-Planck-Landau operators read (α = e, i) We refer to Lemma 7.1 and 7.2 for a precise expansion, 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 (10) are written where ∆ holds for the ionization energy scaled by the thermal energy k B T 0 , so that these operators also depend on ε. Moreover, the factor F 0 has to be rescaled according to the relation where F 0 (which will be later simply denoted by F 0 ) is of order one. We refer to Lemma 7.5 for a development of these operators in terms of ε and also to some properties (weak formulation, entropy, ...) of their leading order terms. Within this framework, we remark that impact ionization is no longer a leading order collisional process, compared with references [4], [5], [6]. As a consequence, the fluid model we will derive at the macroscopic scale will no longer contain any Saha law, like in previous works. The "Saha" law will just appear here as a relaxation term in the right hand side of the equations (see Proposition 2 below).
In order to derive a fluid model, we first consider the above kinetic system at the macroscopic level, more precisely at the electron diffusion scale, and then use a classical Hilbert method; this is the object of the following paragraph. 4. Hilbert expansion at the electron diffusion scaling. We start from the scaled system of kinetic equations (11) with the collision terms given by (13), and rewrite it at the most relevant macroscopic scale, which is the electron diffusion scale (t 0 /ε 2 , x 0 /ε, εE 0 , B 0 ). At this scale, the equations write showing in particular that the magnetic field is a leading order term for the dynamic of the electrons, like in [17], [18]. We point out that the source terms in (16) and (17) are in fact of order 1 ε because Q ε i and Q ε n are both of order ε. The aim is now to look for the limit ε → 0 in this microscopic system, in order to derive a macroscopic model. We first replace in (13) the collision operators by their expansions in terms of ε (see Lemma 7.1, 7.2 and 7.5). Then we expand the solutions in powers of ε according to , α = e, i, n, and insert these expansions in the system (15)- (17). Finally, we identify terms of equal powers of ε, starting with the lowest order terms, which are of order ε −2 for the electrons (in equation (15)) and ε −1 for the other species (i.e. in equations (16) and (17)). We get Concerning the first equation, (18), we deduce thanks to Lemma 7.3 that f 0 e is an isotropic function (with respect to the velocity variable). Then, using the classical mono-species theory [2], we get from (19) while f 0 e is an isotropic function. We remark that, at this stage, the electrons do not seem to have reached their local thermodynamic equilibrium; they will however reach it, as now shown by looking at higher order terms. Identifying terms of order ε −1 in equation (15), and recalling that Q 0 The solvability condition for this equation of unknown f 1 e specifies the equilibrium state for the electrons. We find that the isotropic function f 0 e is actually at local thermodynamic equilibrium, so that it is a centered Maxwellian (see Proposition 1 below, whose proof is postponed to section 6). Moreover, it is possible to compute explicitely the first order correction f 1 e (up to the addition of an arbitrary isotropic function), which is specific to Lorentz operators (see also [18] in the Fokker-Planck context).

Proposition 1.
The distribution function f 0 e is a centered Maxwellian, i.e. of the following form: Introducing the entropic variables ( µe Te , − 1 Te ), where µ e is the chemical potential defined by µ e T e = Log ρ e In this expression, α is the isotropic function defined by Remark 1. As expected, we recover that the velocity vector of an electron after collision is of average zero, so completely randomized. During its free flight, the electron builds up the drift momentum and the related drift kinetic energy along the electromagnetic field. In the next collision, this gained momentum and energy is transferred from drift to random part.
Identifying terms of order ε 0 in equation (17), and remarking that Q 0 ni (f 0 n , f 0 i ) = 0 and Q 0 ne (., f 0 e ) = 0 because f 0 e is an even function, we get the following equation of unknown f 1 n : which does not depend on the other species. Thanks to the classical mono-species Boltzmann theory [2], the solvability conditions for this equation give a classical Euler system for the density, mean velocity and temperature of the neutral species, which is totally independent of the other species (see Proposition 2 below). Let us now look for the ion order one correction f 1 i , by identifying terms of order ε 0 in equation (16). Noting that Q 0 ii (f 0 i , f 0 i ) = 0 and Q 0 ie (., f 0 e ) = 0 still because f 0 e is even, we get the following equation : Again, we look for f 1 i = M u,T φ 1 i , with φ 1 i unknown. Thanks to Lemma 7.4, the solvability condition reduces to one equation on the ion density, with a source term which is directly linked to the non-conservative property of the ionizationrecombination operator Q 0 i,ir . On the contrary, this source term does not depend on the operator Q 0 in , as this last one conserves mass. We then deduce the forthcoming result.
Proposition 2. 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 ): where W n = 1 2 ρ n |u| 2 + 3 2 ρ n T . The ion density ρ i satisfies the following equation: A 1 and A are defined in Lemma 7.5, and the constant a 0 is given by Remark 2. The relation R i = 0, or equivalently is called "Saha law" (see [4], [5]). In the general case, the source term R i is not equal to zero; it can thus be interpreted as a relaxation term to this Saha law.
In order to derive a fluid model for the electrons, we have to go further and identify terms of order ε 0 in the kinetic equation (15); this gives the following equation of unknown f 2 e : as Q 0 en (f 0 e , .) = Q 0 ei (f 0 e , .) = 0. Establishing the solvability condition, we get (see section 6 for the proof): Proposition 3. The electron density ρ e and the temperature T e satisfy the following energy-transport model : where D denotes the diffusion matrix D is a positive definite matrix which sub-blocks are given by Moreover, we have and where the relaxation coefficient λ is defined by Remark 3. The matrix D = D (B) actually depends on the magnetic field B; it satisfies the relation D (B) * = D (−B), which represents the so called "Onsager symmetry" [20]. Moreover, thanks to the explicit computation of Φ 1 e done in Proposition 1, it is also possible to write each sub-block D ij in the following form (with W = |v| 2 /2): where D(W ) is an isotropic matrix. This matrix is obtained as in [18] but with the notations of Lemma 7.3, and setting N (W ) = S W dN (v) = 4π √ 2W . It is defined by (B 1 , B 2 , B 3 ) being the coordinates of the magnetic field B. We remark that each sub-matrix D ij is anisotropic, due to the presence of the magnetic field. It is moreover perfectly computable in terms of the magnetic field B and of the Boltzmann kernel B B (via the isotropic function α). If we introduce for simplicity the three scalar isotropic functions defined by where the different scalars are given by Note that these scalars are all computable in terms of the kernel B B ; they depend on the magnitude of the magnetic field B and also on the neutral density ρ n . In the particular case B = 0 we have D ij = d 0 ij I, i.e. each sub-matrix D ij is scalar, and To summarise, we have shown that the discharge satisfies the following coupled fluid model: Theorem 4.2. The equilibrium states for the heavy species, defined by (21), are characterized by the same mean velocity u and temperature T , while the electron distribution function f 0 e is the centered Maxwellian given by (23). The neutral particles satisfy an Euler system which is totally independant of the other species, ∂ t ρ n + div ρ n u = 0, where W n = 1 2 ρ n |u| 2 + 3 2 ρ n T . The ion density ρ i satisfies the following equation: with R i given by (25). Finally, the electron macroscopic quantities (ρ e , T e ) satisfy the following modified energy-transport model: with R e , U e given by (29), (30), and with u J and v J defined by (27).

5.
A fully conservative fluid model. Starting from the inviscid fluid model derived in Theorem 4.2, the aim is now to take into account corrective terms of order ε in order to obtain a viscous fluid model which would be globally conservative up to the order ε 2 . This can be achieved doing a classical Chapman Enskog expansion.
As the numerical densities of the two charged species are here ε smaller than the numerical density of neutrals, it is in fact sufficient to find corrective terms of order ε only for the neutral fluid system. We obtain the following globally conservative fluid system: Theorem 5.1. The neutral particles satisfy the following Navier-Stokes system: where W n = 1 2 ρ n |u| 2 + 3 2 ρ n T , and π = µ(T )σ(u) with σ(u) the traceless rate of strain tensor (σ(u) ij = ∂ui ∂xj + ∂uj ∂xi − 2 3 divuδ ij ). q is the classical heat flux defined by q = κ(T )∇ x T , and µ(T ) and κ(T ) are the usual viscosity and thermal conductivity of neutrals ( [2], [3]). The source terms are given by with λ defined by (31). The two other species satisfy with j ρe , j Te given by (40), and the source terms R i , and R e , U e by (25), (29) and (30). Moreover, up to the order ε 2 , this coupled fluid model globally conserves mass, momentum and energy.
Proof of Theorem 5.1. The ε viscosity and heat flux terms in the system (41) are quite classical (see [2], [3] and also [17], [18]). It thus remains to compute the corrective source terms, which are nothing but the three first moments of the term of order ε 3 in the asymptotic expansion of Q ε n . As all the elastic collision terms MAGNETIZED AND NON EQUILIBRIUM THERMAL PLASMA 687 conserve mass (see (65) and (66) below), we first have according to (63). Let us now compute S n . Using again (65), it remains According to (63), we get for the first term on the right hand side Then, easy computations show that which gives recalling that The two last terms in S n are finally computed using the conservation relation (68). It gives in particular, for the 0(ε 0 ) term where the second equality results from equation (24). The last equality arises from (38) and the second equation in (37). The following conservation for the ion momentum (up to the order 0(ε)) is then deduced from (47): with a source term which writes

ISABELLE CHOQUET AND BRIGITTE LUCQUIN-DESREUX
Gathering (44), (46) and (47), we get the expression given in (42) for S n . Let us now compute U n . Following the same steps, we first have From (45), we get recalling that Next, using the conservation relation (69), we successively get The last equality results from the following conservation for the ion energy W i = 1 2 ρ i |u| 2 + 3 2 ρ i T (still up to the order 0(ε)): where we have set for simplicity The final expression of U n then results from (48), (49) and (50). To conclude, let us introduce the source term S e given by This is nothing but the momentum of the collisional term of order 1/ε in the kinetic equation (15). The first equality results from the conservation relation (67) and the last one from (46). We then get from the above expressions of the different source terms the following relations: This shows that, up to the order ε 2 , our fluid model globally conserves mass (first equation), momentum (second one) and energy (last one), as the electron energy conserved at the microscopic level is |v| 2 /2 + ∆.
Let us now determine an entropy for this system. We have: with Then, if ε < 2 3 ρn ρi , H is a strictly convex function with respect to (ρ n , ρ n u, W n , ρ i , ρ e , W e ), which satisfies, up to the order ε 2 , the following relation: where P is the non negative term given by The function H is thus an entropy for the whole system (41)-(43) which is compatible with both the diffusion and the source terms.
Remark 5. Let us remark that the unusual term ρeu J Te · (u ∧ B) is a function, so that it can be considered as a source term; it does not load the shocks. The other terms in equation (55) are either non negative or conservative.
Proof of Theorem 5.2. Let us first show (55), starting from the system (41) which can be equivalently written ∂ t ρ n + u · ∇ x ρ n + ρ n divu = εR n , remarking that We deduce from the first and the last equation that Multiplying by ρ n , we get from which we deduce that the function H n given by (52) satisfies Let us now find a similar equation for the ions, starting from the first equation in (43) which gives Following the same approach, we successively get and ∂ t ρ i Log ρ i T 3/2 +div u ρ i Log

MAGNETIZED AND NON EQUILIBRIUM THERMAL PLASMA 691
The function H i defined in (53) thus satisfies Let us finally consider the electrons. The two last equations in (43) give As before, we deduce that Last, multiplying by ρ e , we get from which we deduce that the function H e given by (54) satisfies Gathering (51), (57), (58) and (59), and neglecting the terms of order ε 2 , we get We now consider each term on the right hand side. The last term is clearly non positive; the first one is also. This indeed results from the definition (25) of R e = R i , and the monotonicity of the Logarithm. Now, thanks to (27), the sum of the three terms in the middle writes, up to the factor ε, So that, up to the first term which is conservative, the sum of these three terms is also non positive. We then get (55)-(56), using the definition q = κ∇ x T given in Theorem 5.1. Let us end with the convexity of H. Using classical arguments [13], we can do the proof in the case u = 0. Moreover, as H e does not depend on the neutrals, it is sufficient to show that the two following mappings: (ρ n , W n , ρ i ) → H n + εH i and (ρ e , W e ) → H e are strictly convex. Concerning the second one, easy computations show that the Hessian matrix writes Its two diagonal minors, respectively given by D 1 = 5 2ρe , and D 2 = 3 2W 2 e , are both positive. We thus deduce that this Hessian matrix is positive definite, so that the mapping (ρ e , W e ) → H e is strictly convex. Let us now consider the mapping (ρ n , W n , ρ i ) → H n + εH i , which Hessian matrix is given by The three diagonal minors, respectively denoted by D 1 , D 2 , D 3 , are now computed. We first have D 1 = 1 2ρ 2 n [5ρ n − 3ερ i ], which shows that this minor is positive for ε < 5 3 ρn ρi . The second one writes so that this minor is positive if ε ≤ 2 3 ρn ρi . The third minor is given by The roots of the second polynomial inside the brackets are 3/2 ε and −ε. We deduce that, if ε < 2 3 ρn ρi , this polynomial, and thus D 3 , is positive; which ends the proof.
Proof of Proposition 1. Let us set f 1 e = f 0 e φ 1 e , with φ 1 e unknown. Equation (22) writes L en φ 1 e = S 1 e , where L en is defined by (61) and With the notations of Lemma 7.3, the solvability condition (62) for this equation writes As f 0 e is an even function of the velocity variable v e , we know from Lemma 7.1 iii) that Q 1 en (f 0 e , f 0 n ) is an odd function of v e . Remarking that v e · ∇ x f 0 e − E · ∇ ve f 0 e is also odd, the above relation then reduces to Next, multiplying by the isotropic function Log(f 0 e ) and integrating with respect to the energy variable W , we get Thanks to the classical H-theorem, we deduce that f 0 e is a Maxwellian. This Maxwellian is necessarely centered, i.e. of the form (23), because it is isotropic. The converse is straightforward.
The precise expression of the order one correction f 1 e is a simple adaptation, in the Boltzmann case, of the computation done in reference [18] (see expression (2.58)) for the Fokker-Planck case. The only difference lies in the definition of the isotropic function α.
Proof of Proposition 3. Let us look for f 2 e = f 0 e φ 2 e , with φ 2 e unknown. Equation (26) writes L en φ 2 e = S 2 e , with With the notations of Lemma 7.3, the solvability condition (62) for this equation writes As We remark that the arbitrary isotropic part f 1,e e only appears through the collision term 2Q ee f 1,e e , f 0 e = dQee dfe (f 0 e ) (f 1,e e ). We thus deduce that the condition (60) entirely determines f 1,e e (W ), for any W . Next, multiplying (60) by f 0 e (W ), f 0 e (W )W , and integrating with respect to W , we get in particular We now use the definition of Q 0 en (see Lemma 7.1), and the properties of q B e : self adjointness, kernel made of isotropic functions (see [7]). Setting C g = I R 3 g(v)dv, we then have, for any f and g, In the same way, as Q 1 en (f 0 e , f 1 n ) = q B e (∇f 0 e ) I R 3 v n f 1 n dv n , we also have From [18], 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 (64), we have ii) Let α, β = e, n and α = β. Then The superscript s indicates that a tensor is symmetrized. 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 implies that
Q j en f e , f n (v e ) dv e = 0, ∀j ∈ IN.
In the Fokker-Planck-Landau case, we obtain: ii) Mass conservation implies that Q j ei f e , f i (v e ) dv e = 0, ∀j ∈ IN.
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 with where ρ n = I R 3 f 0 n (v n )dv n is the density of neutral particles. Thanks to a result of reference [1], we know that q B e is a self-adjoint operator on the weighted Hilbert space L 2 M 0,T , and its kernel is made of isotropic functions (with respect to the velocity variable v). Remarking that Γ B φ = 0 for any isotropic function φ, and that Γ B (M 0,T φ) = M 0,T Γ B φ, we deduce that Γ B is an antisymmetric operator on the same space ( [18]), and that the kernel of the operator L en is still made of isotropic functions. We then get the following Lemma which is just a generalization of the result of reference [1] (see also [18] for a similar result in the Fokker-Planck case): Lemma 7.3. (i) 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 ∈ IR 3 , W (v) = W }; we recall the co-area formula , and 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 7.1. We have Proof of Lemma 7.4. The proof partly results from the following weak formulation: We now turn investigating some properties of the ionization-recombination collision operator. In Lemma 7.5 below, we first state a weak formulation, and an entropy inequality, for the dominating part of this operator (in terms of ε). 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 and |v e | 2 − |v e | 2 + |v e | 2 + 2∆ = |v i | 2 − |v n | 2 = 2 ε v n · (v e − v e − v e ) + O(ε 2 ) , so that we deduce the following result, which is a generalization of the Lemma 4.5 of [4]: Lemma 7.5. Let us set ρ n = I R 3 f n (v n )dv n , and ρ i = I R 3 f i (v i )dv i . i) We have, for any α ∈ {e, i, n}, Q α,ir = Q 0 α,ir + O(ε), where the leading order operators Q 0 α,ir are still given by expression (7)- (9). But the conservation equations δ v , δ E in (14) have now been replaced by their limit when ε goes to zero, i.e. by δ 0 E = δ |v e | 2 − [|v e | 2 + |v e | 2 + 2∆] .

ISABELLE CHOQUET AND BRIGITTE LUCQUIN-DESREUX
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: so that we get The function ψ defined by ψ(v) = 1 + |v| 2 2∆ is the only collisional invariant.
8. Annex B : Conservation relations of the collision terms. Concerning the elastic collision operators, we have the following classical conservations [2], [3]: (i) The intra species elastic collision operators conserve mass, momentum and energy, which writes