STUDY OF AN ASYMPTOTIC PRESERVING SCHEME FOR THE QUASI NEUTRAL EULER–BOLTZMANN MODEL IN THE DRIFT REGIME

. We deal with the numerical approximation of a simpliﬁed quasi neutral plasma model in the drift regime. Speciﬁcally, we analyze a ﬁnite volume scheme for the quasi neutral Euler–Boltzmann equations. We prove the unconditional stability of the scheme and give some bounds on the numerical approximation that are uniform in the asymptotic parameter. The proof relies on the control of the positivity and the decay of a discrete energy. The severe non linearity of the scheme being the price to pay to get the unconditional stability, to solve it, we propose an iterative linear implicit scheme that reduces to an elliptic system. The elliptic system enjoys a maximum principle that enables to prove the conservation of the positivity under a CFL condition that does not involve the asymptotic parameter. The linear L 2 stability analysis of the iterative scheme shows that it does not request the mesh size and time step to be smaller than the asymptotic parameter. Numerical illustrations are given to illustrate the stability and consistency of the scheme in the drift regime as well as its ability to compute correct shock speeds.


A fluid model for a quasi-neutral plasma dynamic: the scaled Euler-Boltzmann system
We consider a confined plasma that is made of one species of ions. To model this plasma, we use a fluid approach where the plasma position is denoted x and belongs to the unit periodic segment 3) Here, the constant parameters ε > 0 and γ > 1 denote respectively the thermal energy relatively to the ions kinetic energy and the adiabatic constant. As it is, the system (2.1)-(2.5) is not a closed set of equations, we therefore use the physical hypothesis of quasi-neutrality with adiabatic electrons [24]. Namely, it mathematically writes ∀ (t, x) ∈ [0, T ) × [0, 1], n (t, x) = e φ(t,x) Te with T e > 0 a fixed parameter. (2.6) The set of equations (2.1)-(2.6) now provides a closed system for (n, u, w, p, φ). The hypothesis of quasineutrality with adiabatic electrons enables to reduce the number of unknowns of the system and to re-write it into a set of conservation laws. Indeed, the algebraic relation (2.6) enables to eliminate the electrostatic potential thanks to the following formal computations: from (2.6) we get, φ (t, x) = T e ln (n (t, x)) and from (2.4) we have E (t, x) = −T e ∂xn(t,x) n(t,x) . Therefore the momentum equation writes in the conservative form: Some algebra also enables to write the energy equation (2.3) into the conservative form: ∂ t (w + T e n (ln (n) − 1)) + ∂ x (u (w + p + T e n ln (n))) = 0. (2.8) It stems from writing the right hand side of (2.3) into the sum of a temporal and spatial derivate: nuE = −T e nu ∂ x n n = −T e nu∂ x ln (n) = −T e ∂ x (nu ln (n)) + T e ∂ x (nu) ln (n) , then using the equation (2.1) we obtain: nuE = −T e ∂ x (nu ln (n)) − T e ∂ t (n) ln (n) = −∂ t (T e n (ln (n) − 1)) − T e ∂ x (nu ln (n)) which is enough to get the conservative energy equation (2.8). Finally, (2.1)-(2.6) is formally equivalent to the following system of conservation laws of unknown (n, u, w, p, φ): (P ε ) : ε ∂ x (p + T e n) = 0, ∂ t (w + T e n (ln (n) − 1)) + ∂ x (u (w + p + T e n ln (n))) = 0, w = εnu 2 2 + 1 γ−1 p, φ = T e ln (n) .
The three first equations of (P ε ) stands for the local conservation of mass, momentum and energy while the two last algebraic equations are constitutive relations relating pressure and electrostatic potential to energy and density. The system (P ε ) is supplemented with a periodic initial data n 0 , u 0 , w 0 , p 0 , φ 0 where p 0 and φ 0 obey the constitutive relations.

Conservation properties and admissible set
Smooth periodic solutions to (P ε ) system conserve mass, momentum and energy as stated precisely in the following proposition.
It is also mandatory that solutions to (P ε ) belong to the admissible set Ω = {(n, u, w, p, φ) ∈ R 5 : n > 0, w− εnu 2 2 ≥ 0}. When designing a numerical scheme, it is important to check that the numerical approximation shares the same conservation properties as solutions to the continuous problem. It is a necessary step to assess the robustness of the numerical scheme and often a necessary step towards theoretical convergence result. Significant efforts have been made in this direction in the context of fluid dynamics, for example with the work of Herbin [12,27].

Existence theory for (P ε )
The system (P ε ) can be re-cast into a system of conservation laws for the variables (n, nu, w + T e n (ln (n) − 1)). It has the form: . Therefore (P ε ) is strictly hyperbolic since the Jacobian of the flux f ε has three real distincts eigen values: is the acoustic speed. In the case T e = 0, the system is nothing but the compressible Euler equations for which existence theory is available [9,22,29]. The case T e > 0 should not rise specific difficulties since considering the specific entropy S (n, q, e) = ln p n γ the system (P ε ) admits the additional conservation laws for any h ∈ C 1 (R) and smooth solutions: ∂ t (nh (S)) + ∂ x (nuh (S)) = 0. (2.10) Particularly, if h ∈ C 2 (R) satisfies the inequality h (S) h (S) < 1 γ the couple (−nh (S (n, q, e)) , −nuh (S (n, q, e))) is an entropy-flux pair and therefore the system (P ε ) is symmetrizable which enables to apply the Kato-Friedrichs theory [18].
2.3. The drift limit ε → 0 + The asymptotic ε → 0 + is called the drift limit. In this limit, the pressure force balances the electric-field. To be more precise, let us assume that smooth solutions to (P ε ) admits the formal asymptotic Hilbert expansion: Plugging this expansion into (2.1)-(2.7) and balancing the equal order terms in ε, we get the following system: ∂ t (w 0 + T e n 0 (ln (n 0 ) − 1)) + ∂ x (u 0 (w 0 + p 0 + T e n 0 ln (n 0 ))) = 0, w 0 = 1 γ−1 p 0 , φ 0 = T e ln (n 0 ) , Substituting w 0 in the in the fourth equation and expanding the computation leads to an equation for the pressure that writes: Therefore the formal limit system reads: where π 1 = p 1 +T e n 1 is the first order correction of the total pressure. The system (P 0 ) of unknown (n 0 , u 0 , p 0 , π) is closed. The study of the convergence of the solutions to (P ε ) toward the solutions of (P 0 ) when ε → 0 + is a classical problem arising in fluid mechanics and more generally in the framework of the so called singular limits of hyperbolic systems [1,19,23].

Non conservative formulation P ε
It is sometimes easier to work with a non conservative formulation, that is with a set of variables that is not conserved. As far as the system (P ε ) is concerned, the energy equation (2.3) can be replaced by an equation on the pressure. Indeed, if we assume a smooth solution (n, u, w, p, φ) and then multiply the momentum equation (2.2) by u, we get: Then thanks to the momentum equation (2.7) and the continuity equation (2.1) one has for the right hand side of the second equation: We therefore deduce: Substracting the last equation to the conservative form of the energy equation (2.8) gives the following equation for the pressure: where the right hand side of the last equality vanishes thanks to the continuity equation (2.1). Under the smoothness hypothesis, the system (P ε ) is equivalent to the following non conservative system: with a periodic initial data n 0 , u 0 , w 0 , p 0 , φ 0 where w 0 and φ 0 obey the constitutive relations.
3. Numerical issues related to the use of either (P ε ) or (P ε ) In this section, we briefly discuss the motivation and the numerical issues we care about in this work concerning the numerical simulation of (P ε ).

The three dimension Euler-Lorentz model as a motivation
First and foremost, we mention though it is one dimensional in space, our model contains the essential of the numerical difficulties encountered in the numerical simulation of the three dimensional Euler-Lorentz model for strongly magnetized plasmas. Specifically, in the strong magnetic field regime, the parallel (to the magneticfield) momentum equation degenerates into a balance of force between the parallel gradient of pressure and the parallel electric field. In [7,28] asymptotic preserving schemes are proposed for the three dimensional Euler-Lorentz model. A special care is dedicated to the reformulation of the model so as to ensure that the spatial discretization would yield formally an estimate of the parallel gradient of total pressure that should scale like the ion Larmor radius. This estimate would ensure the asymptotic consistency of the scheme. One can make an analogy with our one dimensional model, to this effect if one assumes the spatial variable x as being the coordinate aligned to the magnetic field. The essential of the discussion in [7,28] is to ensure an estimate of the form ∂ x (p + T e n) = O (ε) in some norm. Linear L 2 stability issues in the asymptotic ε → 0 are also discussed.

Positivity and asymptotic consistency: use of P ε
To derive a numerical scheme that is asymptotically consistent, one has to ensure two things at the discrete level: (a) The positivity of the total pressure p + T e n so as to avoid unphysical gradient.
The use of an equation on the pressure at the continuous level has the advantage to verify a maximum principle, it then gives a hope to ensure the positivity of at the discrete level. As far as the second point is concerned, we mention that it is not mandatory to employ a non conservative form to ensure the formal asymptotic estimate ∂ x (p + T e n) = O (ε), it is also possible to ensure it with the conservative system (P ε ) (see [15] for a proof). Here the non conservative formulation enables to derive an equation on the total pressure p + T e n that is uniformly elliptic in ε [5].
3.3. Asymptotic stability: on the necessity to implicit the gradient of total pressure and the divergence of velocity It is known that the numerical discretization of hyperbolic systems involving a small parameter may suffer from severe stiffness on the numerical discretization parameter [13,26]. A classical approach to overcome this difficulty is to use the so called IMEX splitting schemes [3,31] that is the discrete counterpart of splitting operator techniques: stiff operators are implicit while the a priori non stiff operators are explicit. Here we propose to explain what is the minimal level of implicitness our numerical scheme shall implement to avoid stability problem in the drift regime ε 1. Let us consider a constant solution n 0 , u 0 , w 0 , p 0 , φ 0 to (P ε ). Let us consider a solution (n, u, w, p, φ) to (P ε ) associated with an initial data that is a small perturbation of the constant state. Let then write the solution as: Neglecting the second order terms in (P ε ), it yields the following linearized equations: where the second term of the momentum equation vanishes thanks to the continuity equation, we therefore obtain the linearized equations for ñ,ũ,p,w,φ : To now identify which terms are needed to be implicit, we consider a semi-discretization in time of the form: where ∆t > 0 andñ k ,ũ k ,p k stand for an approximation at time t k = k∆t > 0 ofñ t k , . ,ũ t k , . ,p t k , . . If we take k * = k, and use the linearized momentum equation to expressũ k * as a function of ∂ x p k + T eñ k , therefore we obtain by inserting u k * in both the linearized density and pressure equations, an elliptic equation for the total pressure that takes the form: wherer k is a residual term. Such an explicit discretization would yield a stability constrain that is dependent on ε. From a computational point of view, it is unacceptable. However, if we chose k * = k + 1 the stiff term becomes implicit and the discretization would yield a stability constrain that is dependent on u 0 but not on ε. In this scope, the minimal level of implicitness is: Let us mention however that the implicitness of the gradient of total pressure is not the only way to ensure stability independently on ε. Another way to do it [15,26], is to split the gradient of pressure 1 ε ∂ x p k+1 + T eñ k+1 , and to substitute it by α∂ x p k + T eñ k + 1 ε − α ∂ x p k+1 + T eñ k+1 where α > 0 is a numerical parameter that must be correctly chosen so as to ensure the stability of the scheme.

Shock speed computation for the non conservative formulation P ε and energy conservation
Though it is convenient to work with (P ε ) to ensure the positivity of the pressure (because it is not in this case the difference of two non negative terms) and the formal AP consistency, it brings another difficulty related to the hyperbolic structure of the Euler-Boltzmann model. It is well-known that for hyperbolic systems [20], solutions are rarely classical even if the initial data and the boundary conditions are smooth. Solutions develop discontinuities owing to the fact that characteristics cross in finite time. In this case the pressure equation and the energy equation are not equivalent. Therefore (P ε ) and (P ε ) are no longer equivalent. To palliate this eventuality, we use the strategy developed in [27]. We shall add a corrective source term in the discretization of (2.11) that accounts for the numerical representation of measures appearing when the solution is no longer smooth. It is proven in [27] that under suitable a priori estimates, such an approach allows to recover the Rankine-Hugoniot relations when the mesh size and time step tend to zero. To clarify the discussion, let us briefly sketch what information is lost in the computation of Section 2.4 when solutions are no longer smooth. Consider the simple inviscid Burgers equation posed in (0, +∞) × R: In the case s is smooth, we can multiply by s (3.1) and perform differentiation in time and space to obtain: Assuming now a solution of the form s (t, x) = s l 1 x≤σt (t, x) + s r 1 x>σt (t, x) yields the Rankine-Hugoniot relation for the first equation σ = s l +sr 2 while for the second σ = 2 3 (s l + s r ). It therefore gives different shock speeds. One may ask what is missing in (3.2) and why do we obtain different shock speeds. To understand, let us write the weak formulation of the Burgers equation (3.1). We shall say that s ∈ L ∞ ((0, +∞) × R) satisfies the Burgers equation ( To mimick the formal computation of Section 2.4, or simply to try to obtain (3.2), we are tempted to take ϕ = sψ for some ψ ∈ C 1 c ((0, +∞) × R). However it is not possible since a priori It is easy to see that the left hand side tends to R +∞ 0 s 2 ∂ t ψ + f (s) s∂ x ψdtdx as n → +∞. As for the right hand side, if s∂ t s n +f (s) ∂ x s n is only bounded in L 1 loc ((0, +∞) × R) then up to a subsequence, s∂ t s n +f (s) ∂ x s n tends to a measure. Thus in the case solutions are discontinuous, weak solutions to (3.1) are different from weak solutions to (3.2), they have both different weak formulations. To ensure the consistency between the two formulations at the discrete level, one has to take into account this measure. Our last numerical concern in this work, is the conservation of a discrete analogue of the energy equation (2.3). To obtain the conservation of the energy, one shall perform the formal computation of Section 2.4 in the opposite sense. To do so one has to ensure the discrete analogue of the following properties: • The cancellation at the discrete level: • Make a crucial use of the the continuity equation (2.1): ∂ t n + ∂ x (nu) = 0.
A convenient tool to ensure these two properties is to use a staggered discretization where density and pressure are discretized on a primal mesh, while the velocity is discretized on a dual mesh. It enables to build discrete differential operator such that the duality between the "divergence" of the velocity and the gradient of total pressure holds. Another important property is that the continuity equation needs to be valid both on the primal and the dual mesh.

An unconditionally stable non linear implicit scheme for (P ε )
In this section, we introduce the discretization and propose a non linear implicit scheme to approximate solutions to (P ε ) . Let us begin with the mesh notation. Let a time step ∆t > 0 and define for all k ∈ {0, 1, . . . , T ∆t } the discrete times t k := k∆t. Let N ∈ N * and ∆x := 1 N +1 be the mesh size. The the periodic unit segment line [0, 1] per is discretized with the points defined for all i ∈ Z by x i := i∆x where for all (i, j) ∈ Z 2 , the two points x i and x j denote the same grid point and we note it In the spirit of finite volume framework, we introduce the control cell We also define the cell centered on We use a staggered discretization, namely the density and pressure are discretized on the primal mesh In all the sequel, we shall omit to precise the dependence on T e > 0 and γ > 1 of the solutions to (P ε ) for any ε > 0. For any ε > 0, the discretization of (P ε ) consists in approaching for each k ∈ {0, 1, . . . , T ∆t }, the unknown functions n ε t k , . , p ε t k , . , u ε t k , . by: ⊂ R are assumed to be a solution to the following non linear scheme P ε,∆t,∆x : with the periodic boundary conditions: where all terms will be defined here under. The non linear scheme is supplemented with a periodic initial data For consiceness, we shall use the notation n k ε , p k ε , u k ε ∈ R Z × R Z × R Z to denote the sequences n k ε,i i∈Z ⊂ The staggered discretization (4.1)-(4.11) features two important properties: (a) For all k ∈ {0, . . . , T ∆t }, the duality formula The validity of the discrete continuity equation on the dual mesh T * h : and F k ε,i are defined here under by (4.15) and (4.16).
Now we shall defined precisely, the flux and the density at the interfaces. Our scheme follows the standard finite volumes approach, the discrete continuity (4.5) equation is obtained by integrating the continuity equation (2.1) on C i and considering the discrete values n k ε as an approximation of the mean values of the density at times t k in the cell C i . The flux of mass at the interface is upwind on the density with respect to the sign of the velocity: where (.) + := max (., 0) and (.) − := max (−., 0) denote respectively the positive and the negative part functions.
The discrete momentum equation (4.6) is obtained by integrating momentum equation (2.7) on C i− 1 2 . The density is here defined in the cell C i− 1 2 as the mean value over the two neighbouring cells: The flux of mass and the velocity at the interfaces of the cell C i− 1 2 are defined by: The velocity is here upwind with respect to the sign of the flux of mass: In (4.6), δ i− 1 2 denotes the two points finite difference operator centered at x i− 1 2 , namely, Finally, the discrete pressure equation is obtained by integrating the pressure equation (2.11) in C i . The flux of pressure at the interface is upwind on the pressure with respect to the sign of the velocity: In the discrete pressure equation (4.7), δ i denotes the two point finite difference operator centered at x i , namely (4.20) Lastly, the corrective source term is non negative and here defined for all i ∈ {0, . . . , N } by: The source term S k+1 ε,i is a numerical representation of measures appearing when solutions develop discontinuities. It is designed to compensate in the cell C i the contribution of the residual terms R k+1 ε,i− 1 2 and R k+1 ε,i+ 1 2 that originate from the kinetic energy balance (4.23).
Let us now make some comments on different aspects of the scheme: • The two properties (4.12) and (4.13) together with the presence of the source term (4.21) are crucial to solve the problem of having a conservative scheme that computes the correct shock speeds (see Sect. 3.4). Specifically, the purpose of the duality formula (4.12) is to ensure the cancellation of the sum of the integral of non conservative products. The validity of the continuity equation on the dual mesh (4.13) aims at recovering a discrete total energy balance: namely, we want to mimic the formal computation of Section 2.4 so as to derive an equation for the kinetic energy (2.8). Altogether, these properties ensure the consistency of the scheme P ε,∆t,∆x with (P ε ) when ∆t and ∆x tend to zero. It is necessary, if one wants to recover the Rankine-Hugoniot relations (for further details we refer to [12]). In this respect, we performed discrete time and space differentiation on the discrete momentum equation and continuity equation (4.6), (4.5). Using (4.12) and (4.13), we obtained in Proposition 4.2 a discrete kinetic energy balance with a residual term and in Proposition 4.3 a discrete potential energy balance with another residual term. On the contrary to our formal computation of Section 3.4 where the solution was assumed to be smooth, in Propositions 4.2 and 4.3 the discrete numerical approximation is a priori discontinuous. It thus yields the presence of Dirac measures. The residual terms are therefore the numerical manifestation of this measures. For instance, one can check that if the velocity u ε has a discontinuity then the residual term in Proposition 4.2 does not tend to zero when ∆x and ∆t do. The role of S k+1 ε,i is therefore to compensate the contribution of this residual term on each cell.
• It is proven in Lemmas 2.1 and 3.1 of [11] that such a discretization of the continuity equation and the pressure equation yields the unconditional positivity of the density (provided it is at initial time) and the unconditional non negativity of the pressure (provided it is at initial time). Essentially, the proof relies on the algebraic structure of the discretization of the continuity equation that plays a central role in the derivation of stability estimate. The discrete continuity equation (4.5) is used to prove the non-negativity of the pressure. Without going too much in the details, with the choice of an upwind flux the discrete continuity equation (4.5) can be re-written into the form: is an M -matrix. As far as the nonnegativity of the pressure is concerned, the authors in [11] work with the internal energy. One can re-write the discrete pressure equation (4.7) in terms of the internal energy by decomposing the pressure as p k ε,i = n k ε,i s k ε,i and obtain that the pressure equation re-writes Keeping this in mind, to recover some energy balance, it also requires to implicit the flux of mass also in the momentum equation (2.2). This is why our scheme is fully implicit. It may possible that another approach, as in the spirit of pressure-correction scheme [12] may avoid the implicitness of the flux to obtain stability under a CFL number that does not depend on ε. • The flux of mass and pressure (4.14)- (4.19) at the interfaces can be written under the form: which can be seen as a Rusanov flux where the viscosity only depends on the velocity. In the standard Rusanov flux [21], the viscosity depends on the maximal characteristics speeds of the system which in our case would include the acoustic speed c ε (2.9).
We have the following.

Unconditional stability and uniform in ε bounds
In this section, we are concerned with the well posedness of the scheme P ε,∆t,∆x . More precisely, we are going to prove that: • The scheme is well-defined. Namely it has a solution for any positive numbers ε, ∆t, ∆x.
• Solutions to the scheme are unconditionnaly stable. Namely, for any positive numbers ε, ∆t, ∆x, the scheme preserves the positivity of the density, the non negativity of the pressure, the total mass, the total momentum. It dissipates some discrete energy functional and provides L ∞ bounds on the pressure and density that are uniform in ε as soon as the initial data are O (1) as ε → 0 + .
By periodicity, summing the discrete continuity equation (4.5) and the discrete momentum equation (4.6) yields the conservation of the total mass and momentum, namely for all k ∈ {0, . . . , T ∆t }: Moreover, because the scheme also provides the density to be positive, we control the L 1 -norm of the discrete density n k ε,h and thus any norm since we are working in finite dimensional functional spaces. This therefore gives the L ∞ boundedness of the density. The cornerstone to get the existence of a solution to (P ε,∆t,∆x ) and an L ∞ bound on the pressure, is to obtain an estimate for the discrete energy. More precisely, let us define for a solution n k ε , u k ε , p k ε to the scheme P ε,∆t,∆x the discrete energy functional: To prove the decay of this functional, we need a discrete analogue of the conservative energy equation (2.8). To do so, we crucially uses the duality property (4.12) and the validity of the discrete continuity equation (4.13) on the dual mesh. One has the following.
∆T } and n k ε ∈ R Z , p k ε ∈ R Z , u k ε ∈ R Z a solution to P ε,∆t,∆x . Then the following relation holds for all k ∈ {0, . . . , T ∆t }, i ∈ {0, . . . , N }: where the residual term is given by

(4.24)
Moreover, because of the definition of velocity at the interface, one has R k+1 Proof. It is a consequence of Lemma A.1 by taking the convex function s ∈ R → ψ (s) = s 2 2 . The non negativity of the residual term, follows from the definition of the velocity at the interfaces of each cell C i− 1 2 . Indeed one has by definition We therefore deduce that We obtained a discrete balance for the kinetic energy, but we need a discrete balance for the total energy. Note that in the kinetic energy balance (4.23), we have two non conservative terms: 1 . When carrying the discrete summation of the kinetic energy balance, the sum corresponding to this term is compensated by the sum of the non conservative term of the discrete pressure equation (4.7) thanks to the duality formula (4.12). In order to get rid of the second term, we shall mimick at the discrete level the formal computation of Section 2.4. We have the following.  where the dissipative term is non negative and given by: Proof. Let us consider the function ψ : s ∈ R * + → s (ln (s) − 1) . Let k ∈ {0, . . . , T ∆t }, and i ∈ {0, . . . , N }. Let us multiply the discrete continuity equation (4.5) by ψ n k+1 ε,i , therefore we have: Since the function ψ is C 2 R * + and that the density is everywhere positive, a Taylor-Lagrange expansion around n k+1 ε,i yields the existence of α k+1 i ∈] min n k ε,i , n k+1 Therefore we can re-write the first term: For the second term, we simply rewrite it as: .
Using the definition of the flux, the expansion of the second term of the right hand side yields since the logarithm function is C 2 R * + , a Taylor-Lagrange expansion around both n k+1 ε,i and n k+1 ε,i+1 yields the existence of β k+1 i , γ k+1 i ∈] min n k+1 ε,i+1 , n k+1 ε,i , max n k+1 ε,i+1 , n k+1 ε,i+1 [ such that: We therefore obtain that Collecting all the terms together yields Gathering all the term together yields the result.
Remark 4.4. Note that the dissipation term D k+1 ε,i does not vanish when ∆t and ∆x tend to zero if the density has a discontinuity.
Remark 4.5. The source term S k+1 ε,i is by construction so as We can now prove the following Theorem.
Theorem 4.6 (Existence and unconditional stability). Let ε > 0 and let n 0 Besides, for all k ∈ {0, . . . , T ∆t } the following holds: (e) Provided the initial data n 0 ε , u 0 ε , p 0 ε is O (1) as ε → 0, we have the L ∞ -boundedness uniformly in ε of the density and pressure: sup Proof. We begin with proving the estimates: (a) it suffices to sum the discrete continuity and momentum equation. (b) and (c) We begin with summing the kinetic energy balance (A.1) together with the pressure equation (4.7) (divided by γ −1) and the discrete potential balance (4.25). Then thanks to the periodic boundary conditions and the discrete duality between the divergence of the velocity and the gradient of pressure (4.12), the non conservative terms cancel. Moreover, since by construction where the residual D n+1 ε,i is defined in (4.3). Then it suffices to sum the previous equality from n = 0 to n = k − 1 to get the desired result. (d) Note that for all x ∈ R + * , T e x (ln (x) − 1) ≥ −T e therefore we have, The stability estimate (d) is obtained thanks to the discrete continuity equation (4.5) combined with the conservation of the L 1 -norm and the positivity of the density, one has The last estimate (e) is a consequence of (d). For the sake of simplicity, we now omit to precise the dependence on ε of the discrete density, pressure and velocity. Let us now prove that the scheme is well-defined, we shall apply the Brouwer fixed point theorem. By periodicity, it suffices to prove the existence of the scheme over one period. To do so, let us consider k ∈ {0, . . . , T ∆t } and a solution n k , u k , p k . We look for (n, u, p) ∈ R N +1 × R N +1 × R N +1 solution to: We shall build (n, u, p) as the fixed point of a certain map. Let us introduce for M > 0, which is a convex and compact subset of R N +1 . Let us now denote by S : u ∈ B M → (n (u) , p (u)) ∈ R N +1 * + × R N +1 + → u * (n (u) , p (u)) the mapping that associates with a given velocity u, the density n (u) > 0 solution to (4.28), the pressure p (u) ≥ 0 solution to (4.30) and u * (n (u) , p (u)) solution to (4.29). Note that this map is well-defined since for fixed u ∈ B M , the system of equation on n (u) is linear, invertible and yields n (u) ∈ (R + * ) N +1 , the system of equation on p (u) is also linear, invertible and yields p (u) ∈ (R + ) N +1 and the equation on u * (n (u) , p (u)) is explicit. Moreover, the map S is continuous on B M for every M > 0. It remains to prove that there exists M > 0 (that depends on the mesh size and the previous state at time t k ) such that S (B M ) ⊂ B M . The discrete total energy balance (that is also valid for the system, given the velocity field u ∈ R N +1 ) yields the following where E k is the discrete total energy associated with the solution n k , u k , p k . To conclude, it suffices to prove that we can find M > 0 such that To do so, let us estimate the lower bound inf i=0,...,N n i (u) . One has that n (u) satisfies for all i ∈ {0, . . . , N }: Since n (u) is positive we infer that Because u ∈ B M , one has u + Then we claim that there exists M > 0 (large enough) such that In virtue of the Brouwer fixed-point theorem there exists at least one (n, u, p) solution to (4.28)-(4.31). By induction we deduce that the scheme is well-defined for all k ∈ {0, . . . , T ∆t }. Remark 4.7. Provided the initial data are O (1) as ε → 0 + , the scheme provides bounds that are uniform in ε for the pressure and the density. Up to subsequences, the sequences n k ε , p k ε converges to some limit n k 0 , p k 0 . Nevertheless, it is not enough to conclude the asymptotic consistency of the scheme. We need a uniform estimate on the velocity u k ε . Unfortunately, we have not been able to get more informations on the velocity. In the same spirit, it worth noticing that in the case one has a uniform in ε estimate for the velocity, the kinetic energy must tend to zero as ε → 0 + and the source term S k+1 ε,i and the residual term R k+1 ε,i as well.
5. An iterative linear scheme to solve the non linear one

Definition of the linear scheme
Here we propose an iterative procedure to solve the non linear scheme P ε,∆t,∆x based on a linear implicit scheme. The reason for introducing such an iterative process lies in the fact that, the non linear scheme requires to solve a genuinely non linear system. Standard Newton-like methods may be very cumbersome and delicate to implement. Our motivation here, is therefore to avoid this specific difficulty by offering an iterative scheme where each iteration requires to determine the solution of a linear scheme that is proven to preserve positivity and is linearly L 2 stable independently of the small parameter ε. Let us now define our iterative scheme: given ε > 0 and k ∈ {0, . . . , T ∆t } and n k ε ∈ (R + * ) Z , p k ε ∈ (R + ) Z , u k ε ∈ R Z a solution to P ε,∆t,∆x , we define at time t k the iterative scheme I k P ε,∆t∆x : where the integer r here refers to the internal number of iterations. To lighten the notation we shall omit in this section the dependence on ε of the internal iteration n r,k , p r,k , u r,k . Provided the density remains positive, at each iteration r ∈ N, the velocity u r+1,k is defined from equation (5.4) depends implicitly on p r+1,k and n r+1,k but explicitly on the flux. Here the mass flux is defined at each iteration r ∈ N by: and the source term is here given by Other quantities are defined similarly as for the non linear scheme (4). Note that if the sequence n r,k , p r,k , u r,k ∈ R Z × R Z × R Z converges to some (n, u, p) ∈ R Z × R Z × R Z then (n, u, p) solves the non linear scheme P ε,∆t,∆x at the iteration k + 1.

Reduction to an elliptic system and positivity conservation
Let k ∈ {0, . . . , T ∆t }, provided the density is positive, each internal iteration r ∈ N of I k P ε,∆t∆x is equivalent to solve a linear elliptic system of unknown n r+1,k , p r+1,k . Then the update of the velocity u r+1,k becomes an explicit step. For clarity of exposure, we shall defined some discrete finite difference operators. Namely, for a vector q r ∈ R Z , and for i ∈ Q the discrete centered Laplacian operator centered at x i is defined by: We also define the the diffusion operators: We can now re-write the iterative scheme. Thanks to the momentum equation (5.4) one has for all r ∈ N and for all i ∈ {0, . . . , N } : If n r,k i± 1 2 > 0 one has: . Plugging this expression into the pressure equation (5.5) yields ε (∆x) 2 ∆ i p r+1,k + T e n r+1,k = n k ε,i +n r,k i , It thus yields a linear elliptic system of unknown n r+1,k , p r+1,k that writes given by (5.19),p r,k i given by (5.17), and, ∀i ∈ Z, n r+1,k i = n r+1,k i+(N +1) , p r+1,k i = p r+1,k i+(N +1) .
Provided n r is positive, there exists a unique n r+1,k , p r+1,k ∈ R Z × R Z solution to L r+1,k ε .

M. BADSI
Remark 5.1. In the limit ε → 0 + , the system L r+1,k ε degenerates into which is an ill-posed problem since any constant vector is solution. This difficulty is overcomed using a Duality-Based decomposition [25] which here consists in the decomposition of the solution into L r+1,k ε into its mean and fluctuation.
We are now going to prove that the solution to L r+1,k ε is such that p r+1,k + T e n r+1,k ∈ (R + * ) Z . To do so, we shall need the following Lemma.
There exits δ r,k > 0 that depends only on n k ε , p k ε , u k ε , n r,k , p r,k , u r,k and T e > 0 such that if 0 < ∆t ∆x < δ r,k then p k ε,i + T e n k ε,i +p r,k i + T en r,k i > 0 for all i ∈ Z.
Remark 5.7. Note that the veracity of the condition: there is 0 < δ < 1 such p0 Ten0+(γ−1)p0 < 1 − δ depends on the kind of physics we consider. Typically, if we write the ionic pressure as p 0 = T i n 0 where T i is a ionic reference temperature then the inequality writes: It therefore depends both on the temperature ratio Te Ti and the value of γ. In the plasma physics context, fluid models are derived from kinetic ones and γ takes the form γ = 2 dv + 1 where d v is the dimension of the velocity space. When d v ≥ 2, the stability essentially depends on the temperature ratio.
Remark 5.8. If one replaces in the previous analysis u 0 = 0, one can show that the scheme is unconditionally stable.

Numerical results
The aim of this section is to illustrate the main features of our non linear scheme P ε,∆t,∆x and to test its ability to be stable and consistent in the drift regime ε 1. We shall explain how we proceed numerically: we fix N ∈ N, T > 0 and take ∆x = 1 N +1 . For k ∈ {0, .., T ∆t } and r ∈ N we define the relative l 2 error at the r-th iteration of the solver I k P ε,∆t,∆x by e r+1,k k := W r+1,k −W r,k 2 W r,k 2 where we denote W r,k = n r,k , u r,k , p r,k the solution of I k P ε,∆t,∆x at iteration r ∈ N and where . 2 denotes the euclidian norm on R 3(N +1) . For each iteration k ∈ {0, . . . , T ∆t }, the iterative solver I k P ε,∆t,∆x computes at each internal iteration r ∈ N the time step according to the CFL condition of Lemma 5.2, that is ∆t = ∆x min i∈{0,...,N } sup{λ > 0 | p k ε,i + T e n k ε,i − p r,k −,i (λ) − T en r,k −,i (λ) > 0}. The iterative scheme stops iterating when there is r ∈ N such that e r+1 k < 10 −10 . In all the numerical test cases presented under, the electron temperature is T e = 1.

Shock speed computation in ε = O (1) regime
Here we test the ability of our scheme to capture the correct shock speed. For this test case, homogeneous Neumann boundary conditions are prescribed on each boundary. We consider the Riemann problem with the initial condition defined for all x ∈ [0, 1] by n 0 (x) = 1.0, p 0 (x) = 10 3 − 1 1 x≤0.5 (x) + 10 −2 − 1 1 x>0.5 (x) , u 0 (x) = 0. (6.1) We discretize the problem with ∆x = 1 512 and we choose ε = 1 and γ = 1.4. We represent in Figure 1 the exact solution and the approximation at time T = 0.012 and the evolution of the discrete energy (4.22). We see that the scheme computes the correct shock speed and is slightly diffusive. Refining the mesh size diminishes the diffusion. The energy as expected theoretically is non increasing.

The drift regime ε 1
For this test case, we consider an initial condition of the form: ∀x ∈ [0, 1], n 0 ε (x) = 1 + ε cos (2πx) , u 0 ε (x) = 1.0, p 0 ε (x) = 1 − ε cos (2πx) , (6.2) with γ = 3. We discretize the problem with a spatial step ∆x = 1/64 and a temporal step that is fixed ∆t = ∆x 100 , then the CFL is fixed at 0.01. We look at ε values of 10 −4 and 10 −6 . We compare the numerical solution to P ε,∆t,∆x with the analytical solution to (P ε ) which is given for all t ≥ 0 and x ∈ [0, 1] per by n ε (t, x) = n 0 ε (x − t), u ε (t, x) = 1.0, p ε (t, x) = p 0 ε (x − t) . In Figure 2, we represent the error between the analytical solution and the numerical approximation given by P ε,∆t,∆x at time T = 1.0. We observe that there is a little error between the limit solution and the numerical approximation. Note that in the case ε = 10 −6 , the CFL ∆t ∆x is ten times greater than √ ε and the scheme is still able to compute a solution which evidences its asymptotic stability.

Spatial convergence of the scheme
Here we measure the spatial convergence of P ε,∆t,∆x towards an analytical solution. The test case we pay attention to, is a pure transport problem where the velocity remains constant. For ε > 0 given, the functions defined for t ≥ 0 and x ∈ [0, 1] per by: n ε (t, x) = n 0 ε (x − t) , p ε (t, x) = p 0 ε (x − t) , u ε (t, x) = 1.0, where p 0 ε (x) = 1.0 + ε cos (2πx), n 0 ε (x) = −p 0 ε (x)+2.0 Te , and u 0 ε (x) = 1 are solutions to (P ε ) . We perform a spatial convergence study where the CFL is 0.1. We represent in Figure 3, the total L 1 -error between the analytical solution and its approximation as a function of ∆x at the final time T = 0.01 for ε ∈ {1, 10 −2 , 10 −4 }. We measure a spatial convergence rate equal to one for all the curves, which is the best we could expect from the scheme since the exact solution is smooth and the flux are first order consistent.

Conclusion
We considered the numerical discretization of the quasi-neutral Euler-Boltzmann equations. We proposed a non linear implicit scheme based on staggered grids that was proven to be unconditionally stable and to provide some uniform bounds with respect to ε. Because of its non linearity, we proposed an iterative linear implicit scheme to solve it. The iterative scheme was proven to preserve the positivity and to be L 2 linearly stable under a CFL condition that does not involve ε. We test the ability of the scheme to compute the correct shock speed through a Riemann problem and we showed that the scheme is stable when the CFL number is larger than √ ε. Let us mention some perspectives of amelioration of this work. Due to the lack of uniform estimate for the velocity, the theoretical convergence of the scheme in the drift regime is not proven. This problem is postponed to a future work. Extension of the scheme to higher order in space flux is also a possible perspective. Eventually, the extension of the scheme to a more realistic magnetized plasma model in a three dimensional geometry is an on going work.
where the residual term is given by In particular, since ψ is convex and by definition of the flux, the residual term is non negative, R k+1 i− 1 2 ≥ 0.