Denise Aregba-Driollet

Stéphane Brull

About viscous approximations of the bitemperature Euler system

published or not. The documents may come

Introduction

This paper is devoted to a viscous approximation of the bitemperature Euler system that has been studied in [START_REF] Aregba | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF]. This fluid model describes the interaction of a mixture of one species of ions and one species of electrons in thermal nonequilibrium, with applications in the field of Inertial Confinement Fusion where solutions with shocks occur. Quasineutrality being assumed, the electronic and ionic mass fractions are constant: subscripts e and i standing for electron and ions respectively,

⇢ e = m
e n e = c e ⇢, ⇢ i = m i n i = c i ⇢, c e + c i = 1 and the model consists of two conservation equations for mass and momentum and two nonconservative equations for each energy.

Moreover the pressure of each species is supposed to satisfy a gamma-law with its own constant:

p e = ( e 1)⇢ e " e = n e k B T e , p i = ( i 1)⇢ i " i = n i k B T i , (1.1)
where k B is the Boltzmann constant, " ↵ and T ↵ represent respectively the internal specific energy and the temperature of species ↵, ↵ 2 {e, i}.

The total energies are given by E

↵ = ⇢ ↵ " ↵ + 1 2 ⇢ ↵ u 2 , ↵ 2
{e, i}. We denote ⌫ ei 0 the interaction coe cient between electronic and ionic temperatures. The bitemperature Euler system is the following: 8 > > > > < > > > > :

@ t ⇢ + @
x (⇢u) = 0, @ t (⇢u) + @ x (⇢u 2 + p e + p i ) = 0, @ t E e + @ x (u(E e + p e )) u(c i @

x p e c e @

x p i ) = ⌫ ei (T i T e ), @ t E i + @ x (u(E i + p i )) + u(c i @ x p e c e @ x p i ) = ⌫ ei (T i T e ).
(1.2)

A first step in the comprehension of this system is to suppose that e = i = . In this case, one can define a global internal energy " = c

e " e + c i " i which satisfies p e + p i = ( 1)⇢". Denoting E = E e + E i the total energy one has

E = ⇢" + 1 2 ⇢u 2
and @ t E + @ x (u(E + p)) = 0 so that (⇢, ⇢u, E) satisfies the usual Euler 3 ⇥ 3 system with law. Nevertheless, even in this case one needs to solve a nonconservative equation in order to get T e and T i separately. In our context, the nonconservativity is not only due to source terms but especially to terms multiplying u by pressure gradients, making delicate the definition of admissible shocks. In order to define nonconservative products, Dal Maso, Le Floch and Murat proposed in [START_REF] Dal Maso | Definition and weak stability of nonconservative products[END_REF] a new theory based on the definition of family of paths. In [START_REF] Coquel | Numerical methods for weakly ionized gas[END_REF], the authors consider the bitemperature Euler system with di↵usive terms. By assuming that the electrons are isentropic, the system is transformed into a conservative model. In [START_REF] Sainsaulieu | Equilibrium velocity distribution functions for a kinetic model of twophase flows[END_REF], the authors consider a kinetic system for sprays and derive a nonconservative hyperbolic system that is studied in [START_REF] Raviart | A nonconservative hyperbolic system modeling spray dynamics. I. Solution of the Riemann problem[END_REF].

In [START_REF] Aregba | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF], the Euler bitemperature system has been derived by hydrodynamic limit of an underlying kinetic model which consists of a BGK model coupled with Poisson equation in the quasineutral regime. Moreover the obtained fluid system have been proved to be entropy dissipative by a direct approach and also by using the Boltzmann entropy. In particular, the nonconservative terms are obtained from the definition of the electric field according to a generalized Ohm's law. Numerical simulations have been performed with a finite volume scheme designed from this kinetic model, and have been proved to be physically realistic. Relaxation schemes based on the Suliciu approach have also been developped in [START_REF] Aregba | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF] for this sytem and in [START_REF] Brull | Modelling and entropy satisfying relaxation for the nonconservatioce Euler system with transverse magnetic[END_REF] in the case of a transverse magnetic configuration. In [START_REF] Wargniez | Numerical treatment of the nonconservative product in a multiscale fluid model for plasmas in thermal nonequilibrium: application to solar physics[END_REF], the authors perform a Chapman-Enskog expansion by introducing a small parameter representing the ratio between electronic and ionic molecular masses in the spirit of [START_REF] Graille | Kinetic theory of plasmas:. Maths models and methods in the[END_REF]. They obtain an hyperbolic system with a parabolic regularisation on the electrons.

One of the motivation of the derivation of viscous terms is the construction of traveling waves for the determination of shock profiles. In particular, in the context of a nonconservative system for two-phase flows, the authors construct in ( [START_REF] Raviart | A nonconservative hyperbolic system modeling spray dynamics. I. Solution of the Riemann problem[END_REF], [START_REF] Sainsaulieu | Ondes progressives solutions de systèmes convectifs-di↵usifs et systèmes hyperboliques non conservatifs[END_REF]) traveling waves solutions. The construction of traveling waves is also considered in [START_REF] Wargniez | Numerical treatment of the nonconservative product in a multiscale fluid model for plasmas in thermal nonequilibrium: application to solar physics[END_REF] by using a di↵usive equations on the electrons. Therefore this work is a first step before the construction of travelling waves for the bitemperature Euler system.

More precisely, in the present paper, we perform a Chapman-Enskog expansion of our kinetic model up to order one in order to get rigorously a viscous, Navier-Stokes type approximation of the bitemperature Euler system in the case e = i . As a result, we obtain conservative and nonconservative second order terms. To go into details, let us denote U = (⇢, ⇢u, E e , E i ). The Euler bitemperature system (1.2) being written in condensed form as

@ t U + A(U )@
x U = S(U ), for a fixed relaxation parameter ⌧ > 0 the obtained second order system can be written under the form

@ t U ⌧ + A(U ⌧ )@ x U ⌧ = S(U ⌧ ) + ⌧ (u ⌧ @ x (J(U ⌧ )@ x U ⌧ ) + @ x (D(U ⌧ )@ x U ⌧ )) . (1.3)
Here J(U ⌧ ) and D(U ⌧ ) are 4 ⇥ 4 matrices, while u ⌧ is the velocity. This result completes known models such as the one studied by C. Chalons and F. Coquel in [START_REF] Chalons | Navier-Stokes equations with several independent pressure laws and explicit predictor-corrector schemes[END_REF], by constructing rigorously some second order terms to their system. Next we prove the compatibility of the entropy of the bitemperature Euler system with the di↵usive terms. We recall that a dissipative entropy ⌘ exists for (1.2), namely

⌘(U ) = ⌘ e (⇢c e , " e ) + ⌘ i (⇢c i , " i ), ⌘ ↵ (⇢ ↵ , " ↵ ) = ⇢ ↵ m ↵ ( ↵ 1) ln ✓ p ↵ ⇢ ↵ ↵ ◆ , ↵ 2 {e, i} (1.4)
with flux Q(U ) = u⌘(U ) [START_REF] Aregba | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF]. We prove here that the solutions U ⌧ of (1.3) formally satisfy the following inequality:

@ t ⌘(U ⌧ ) + @ x (u⌘(U ⌧ ))  ⌫ ei k B T i T e (T i T e ) 2 ⌧ 5k B 2m ↵ X ↵=e,i @ x (n ↵ @ x T ↵ ) .
As a consequence, if (U ⌧ ) ⌧ is smooth and decreases fast enough at infinity one has

d dt Z R ⌘(U ⌧ (x, t))dx  ⌫ ei k B Z R (T ⌧ i T ⌧ e ) 2 T ⌧ i T ⌧ e dx .
This is a first step to prove that U ⌧ owns a limit U which is a weak entropy solution of the Euler bitemperature system, see [START_REF] Serre | Systems of conservation laws[END_REF] for detailed argumentation. In the limit we have

@ t ⌘(U ) + @ x Q(U )  ⌫ ei k B T i T e (T i T e ) 2 . (1.5)
This inequality is the one satisfied formally by a limit of the moments of the solution of the BGK system [START_REF] Aregba | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF].

In the case of BGK approximations of systems of conservation laws, the entropy dissipativity of the Chapman-Enskog expansion is wellknown, see [START_REF] Bouchut | Construction of BGK models with a family of kinetic entropies for a given system of conservation laws[END_REF], [START_REF] Cuesta | Kinetic shock profiles for nonlinear hyperbolic conservation laws[END_REF]. As a matter of fact if F is the flux function of such a system, the Chapman-Enskog expansion is

@ t U ⌧ + @ x F(U ⌧ ) = ⌧ @ x (D(U ⌧ )@ x U ⌧ ) (1.6) with D(U ) = 0 (U ) F 0 (u) F 0 (
u), the matrix 0 (U ) depending on the chosen BGK framework. Then one multiplies (1.6) by ⌘ 0 (U ⌧ ) and uses the equality

⌘ 0 (U )@ x (D(U )@ x U ) = @ x ⌘ 0 (U )D(U )@ x U ⌘ 00 (U ) (@ x U, D(U )@ x U ) . In this case, t D(U )⌘ 00 (U ) is symmetric nonnegative so that @ t ⌘(U ⌧ ) + @ x Q(U ⌧ )  ⌧ @ x ⌘ 0 (U ⌧ )D(U ⌧ )@ x U ⌧ .
In the case of the present paper, A(U ) is not a gradient, so we had to perform a direct calculation to prove the dissipation property. The paper is organized as follows. The Section 2 deals with the derivation of a Navier-Stokes system starting from the kinetic system proposed in [START_REF] Aregba | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF]. In section 3, the di↵usive terms are shown to be dissipative w.r.t. the entropy of the bitemperature Euler system. Finally, section 4 gives conclusions to this work.

2 Derivation of the Navier-Stokes system 2.1 Notations Kinetic models are described by the distribution function f ↵ of each species depending on the time variable t 2 R + , on the position x 2 R 3 and on the velocity v 2 R 3 . The macroscopic quantities can be obtained by extracting moments on these distribution functions w.r.t the velocity variable. Indeed density, velocity and total energy of the species ↵ can be defined as

n ↵ = Z R 3 f ↵ dv, u ↵ = 1 n ↵ Z R 3 v 1 f ↵ dv, E ↵ = 3 2 ⇢ ↵ k B m ↵ T ↵ + 1 2 ⇢ ↵ u 2 ↵ = Z R 3 m ↵ v 2 2 f ↵ dv. (2.7)
The present model is monoatomic ( = 5 3 ). Hence, the internal specific energy of species ↵ writes

" ↵ = 3 2m ↵ k B T ↵ .
In the following we shall use the moment operator P ↵ defined by

P ↵ (f ↵ ) = m ↵ Z R 3 0 @ 1 v 1 v 2 2 1 A f ↵ dv. (2.8) We denote P ↵ (f ↵ ) = U ↵ : U ↵ = 0 @ ⇢ ↵ ⇢ ↵ u ↵ E ↵ 1 A . (2.9)
Usually the velocity and the temperature of the mixture are defined by

u = ⇢ e u e + ⇢ i u i ⇢ e + ⇢ i , nk B T = X ↵ ( 1 2 ⇢ ↵ (u 2 ↵ u 2 )) + X ↵ (n ↵ k B T ↵ ), (2.10)
where n = n e + n i . Moreover, the current of the plasma j and the total charge ⇢ are defined by

⇢ = Z R 3 (q e f e + q i f i ) dv = n e q e + n i q i , j = Z R 3 v 1 (q e f e + q i f i ) dv = n e q e u e + n i q i u i . (2.11)

Chapman-Enskog expansion

Consider the following kinetic model in the quasi-neutral regime 8 < :

@ t f ↵ + v 1 @ x f ↵ + q↵ m↵ E@ v 1 f ↵ = 1 ⌧ (M ↵ f ↵ ) + 1 ⌧ ei (M ↵ f ↵ ), @ t E = j ⌧ 2 , @ x E = ⇢ ⌧ 2 , (2.12)
where ⌧ is a positive parameter proportional to the Knudsen number. M ↵ and M ↵ are the two Maxwellian distribution functions

M ↵ (f ↵ ) = n ↵ (2⇡k B T ↵ /m ↵ ) 3/2 exp( |v u ↵ | 2 2k B T ↵ /m ↵ ), ↵ = e, i, (2.13)
M ↵ (f e , f i ) = n ↵ (2⇡k B T /m ↵ ) 3/2 exp( |v u| 2 2k B T /m ↵ ), ↵ = e, i. (2.14)
Next we perform a first order Chapman-Engskog expansion up to order 1. Hence the solution of the system (2.12, 2.13, 2.14) f ↵ is researched as the expansion

f ↵ = M ↵ + ⌧ g ↵ , ↵ 2 {e, i}, (2.15)
with the constraints Z

R 3 f ↵ dv = Z R 3 M ↵ dv, Z R 3 v 1 f ↵ dv = Z R 3 v 1 M ↵ dv, Z R 3 v 2 f ↵ dv = Z R 3 v 2 M ↵ dv. (2.16) So Z R 3 g ↵ dv = 0, Z R 3 v 1 g ↵ dv = 0, Z R 3 v 2 g ↵ dv = 0.
(2.17)

Euler system

At this level, it is possible to compute the Euler system by arguing as in ( [START_REF] Aregba | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF]). Hence by plugging the expansion (2.15, 2.16, 2.17) into (2.12) and considering terms up to order 1, it holds that

@ t M ↵ + v 1 @ x M ↵ + q ↵ m ↵ E@ v 1 M ↵ = g ↵ + 1 ⌧ ↵ (M ↵ M ↵ ) + O(⌧ ), with j = O(⌧ 2 ) and ⇢ = O(⌧ 2 ). Hence by using that Z R 3 m ↵ M ↵ dv = ⇢ ↵ , Z R 3 m ↵ v 1 M ↵ dv = ⇢ ↵ u, Z R 3 m ↵ v 2 M ↵ dv = ⇢ ↵ " ↵ + 1 2 ⇢ ↵ u 2 ,
we get 8 > > < > > :

@ t ⇢ + @ x (⇢u) = O(⌧ ), @ t (⇢u) + @ x (⇢u 2 + p e + p i ) = O(⌧ ), @ t (⇢ e " e + 1 2 ⇢ e u 2 ) + @ x (u(⇢ e " e + 1 2 ⇢ e u 2 + p e )) u(c i @ x p e c e @ x p i ) = ⌫ ei (T i T e ) + O(⌧ ), @ t (⇢ i " i + 1 2 ⇢ i u 2 ) + @ x (u(⇢ i " i + 1 2 ⇢ i u 2 + p i )) + u(c i @ x p e c e @ x p i ) = ⌫ ei (T i T e ) + O(⌧ ). ( 2
g e = ⇣ (v 1 u) ✓ @ x n e n e 3 2
@ x T e T e ◆ + @ x u (v 1 u) 2 k B me T e 1 3 (v u) 2 k B me T e ! + ⌫ ei n e k B T e (T i T e ) (v u) 2 3 k B me T e 1 ! (v 1 u) (v u) 2 2 k B me @ x ( 1 T e ) (v 1 u) k B me T e ⇢ @ x (p e + p i ) ⌘ M e + q e m e E@ v 1 M e 1 ⌧ ei M e M e ! , (2.19)
g i = ⇣ (v 1 u) ✓ @ x n i n i 3 2 @ x T i T i ◆ + @ x u (v 1 u) 2 k B m i T i 1 3 (v u) 2 k B m i T i ! + ⌫ ei n i k B T i (T e T i ) (v u) 2 3 k B m i T i 1 ! (v 1 u) (v u) 2 2 k B m i @ x ( 1
T i ) (v 1 u) k B m i T i ⇢ @ x (p e + p i ) ⌘ M i + q i m i E@ v 1 M i 1 ⌧ ei M i M i ! . ( 2

.20)

Proof. g ↵ is given by the relation

g ↵ = ✓ @ t M ↵ + v 1 @ x M ↵ + q ↵ m ↵ E@ v 1 M ↵ 1 ⌧ ei (M ↵ M ↵ ◆ .
A direct computation gives

@ t M e = ( @ t n e n e 3 2
@ t T e T e ) + (v 1 u) @ t u k B me T e (v u) 2 2 k B me T e @ t ✓ 1
T e ◆ ! M e (2.

21)

and

v 1 @ x M e = ( v 1 @ x n e n e 3 2 v 1 @ x T e T e ) + (v 1 u) v 1 @ x u k B me T e (v u) 2 2 k B me T e v 1 @ x ✓ 1 T e ◆ ! M e .
By using the non-conservative Euler system (2.18), the time derivatives of (2.21) are computed in function of the space derivatives up to O(⌧ ) terms, as follows

@ t u = u@ x u 1 ⇢ @ x (p e + p i ) + O(⌧ ), @ t n i n i 3 2 @ t T i T i = u @ x n i n i + 3 2 u @ x T i T i ⌫ ei n i k B T i (T e T i ) + O(⌧ ), @ t T e T e = u @ x T e T e 2 3 @ x u + 2 3 ⌫ ei n e k B T e (T i T e ) + O(⌧ ) @ t T i T i = u @ x T i T i 2 3 @ x u + 2 3 ⌫ ei n i k B T i (T e T i ) + O(⌧ ).
Hence up to O(⌧ ) order terms, we get

@ t M e + v 1 @ x M e = u @ x n e n e + 3 2 u @ x T e T e ⌫ ei n e k B T e (T i T e ) + (v 1 u) k B me T e ✓ u@ x u 1 ⇢ @ x (p e + p i ) ◆ + (v u) 2 2 k B me T e ( u @ x T e T e 2 3 @ x u + 2 3 ⌫ ei n e k B T e (T i T e )) ! M e .
Therefore, we obtain

@ t M e + v 1 @ x M e = (v 1 u)( @ x n e n e 3 2

@

x T e T e ) + @ x u

(v 1 u) 2 k B me T e (v u) 2 3 k B me T e ! + ⌫ ei n e k B T e (T i T e ) (v u) 2 3 k B me T e 1 ! (v 1 u) (v u) 2 2 k B me @ x ( 1 T e ) (v 1 u) k B me ⇢T e @ x (p e + p i ) ! M e
and we recover (2.19). The same result holds for (2.20).

Proposition 2. The viscous approximation of the kinetic system (2.12, 2.13, 2.14) writes

@ t ⇢ + @ x (⇢u) = 0, (2.22) @ t (⇢u) + @ x (⇢u 2 + p e + p i ) 4 3 ⌧ @ x (p @ x u) = 0, (2.23) @ t (⇢ e " e + 1 2 ⇢ e u 2 ) + @ x (u(⇢ e " e + 1 2 ⇢ e u 2 + p e )) u(c i @ x p e c e @ x p i ) u ✓ 4 3 ⌧ c e @ x (p i @ x u) 4 3 ⌧ c i @ x (p e @ x u) ◆ 4 3 ⌧ @ x (up e @ x u) 5 2 ⌧ @ x ( k B m e p e @ x T e ) = ⌫ ei (T e T i ), (2.24) @ t (⇢ i " i + 1 2 ⇢ i u 2 ) + @ x (u(⇢ i " i + 1 2 ⇢ i u 2 + p i )) + u(c i @ x p e c e @ x p i ) +u ✓ 4 3 ⌧ c e @ x (p i @ x u) 4 3 ⌧ c i @ x (p e @ x u) ◆ 4 3 ⌧ @ x (up i @ x u) 5 2 ⌧ @ x ( k B m i p i @ x T i ) = ⌫ ei (T i T e ), (2.25)
where the electric field E is given by

✓ n e q e ⇢ e n i q i ⇢ i ◆ E = ⇢ ⇢ e ⇢ i n e q e E = ⇢ ⇢ e ⇢ i n i q i E = @ x p e ⇢ e @ x p i ⇢ i 4 3 ⌧ ⇢ e @
x (p e @ x u) + 4 3

⌧ ⇢ i @ x (p i @ x u) .
(2.26) Remark 1. In Proposition 2, the kinetic model is considered in the monoatomic case. This means that e = i = 5 3 . The polyatomic case can be managed by considering a kinetic model with internal energy variable like in [START_REF] Aregba | Construction and approximation of the polyatomic bitemperature Euler system[END_REF].

Remark 2. The relation (2.26) is an approximation at order ⌧ of the Ohm law given in [START_REF] Aregba | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF]. The nonconservative terms of the system (2.22, 2.23, 2.24, 2.25) are shown to appear from this relation defining E Proof. By expanding f ↵ like in (2.15), we get

@ t (M ↵ + ⌧ g ↵ ) + v 1 @ x (M ↵ + ⌧ g ↵ ) + q ↵ m ↵ E@ v 1 (M ↵ + ⌧ g ↵ ) = g ↵ + 1 ⌧ ei (M ↵ M ↵ ⌧ g ↵ ) (2.27)
and we extract the moments of (2.27) w.r.t 1, v, v 2 . From the orthogonality conditions (2.17) the viscous part of the expansion will be given by the moments of v 1 @ x g ↵ w.r.t 1, v, v 2 . Firstly, we obtain the same mass conservation equation as for the Euler system because from the relation (2.17).

In order to get the impulsion equation (2.23), we compute

Z R 3 m e v 2 1 g e dv = Z R 3 m e (v 1 u) 2 g e dv.
By using the expression of g e (2.19) and some orthogonality relations, it holds that Z

R 3 m e (v 1 u) 2 g e dv = @ x u Z R 3 (v 1 u) 4 k B me T e 1 3 (v 1 u) 2 (v u) 2 k B me T e ! M e m e dv.
Moreover, by using Z

R 3 (v u) 4 M e dv = n e ( k B T e m e ) 2 ,
we get finally Z

R 3 m e (v 1 u) 2 g e dv = 4 3 p e @
x u.

(2.28)

Hence we get for species

↵ @ t (⇢ ↵ u) + @ x (⇢ ↵ u + p ↵ ) q ↵ n ↵ E 4 3 ⌧ @ x (p ↵ @ x u) = 0. (2.29)
Then by adding the previous equations for each species (2.23). Next in order to establish the relation defining E, (2.29) is rewritten in the nonconservative form as

@ t u + u@ x u + 1 ⇢ ↵ @ x p ↵ q ↵ n ↵ ⇢ ↵ E 4 3 ⌧ ⇢ ↵ @ x (p ↵ @ x u) = 0. (2.30)
By substracting equation (2.30) for electrons and ions and by proceeding like in [START_REF] Aregba | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF], we get the generalized Ohm's law (2.26). By using the relation (2.17), we obtain Z

R 3 m e 2 v 1 v 2 g e dv = Z R 3 m e 2 (v 1 u)(v u) 2 g e dv + u Z R 3 m e v 2 1 g e dv.
By using orthogonality relations, it comes that

Z R 3 m e 2 (v 1 u)(v u) 2 g e dv = 5⇢ e ✓ k B m e T e ◆ 2 ( @ x n e n e 3 2
@ x T e T e ) + @ x ( 1 T e ) 35⇢ e 2 k B me ( k B m e T e ) 3 + Z R 3 (v 1 u) 2 (v u) 2 k B me ⇢T e @ x (p e + p i )M e dv + 5n
@ t (⇢ i " i ) + @ x (u(⇢ i " i )) + p i @ x u + 4 3 ⌧ @ x (up i @ x u) 5 2 ⌧ @ x ( k B m i p i @ x T i ) 4 3 ⌧ p i (@ x u) 2 = ⌫ ei (T i T e ).
(2.34)

Proof. Firstly, a straightforward computation gives

1 2 @ t (⇢ i u 2 ) = 1 2 @ x (⇢u 3 ) u@ x p i + q i n i uE + 4 3 ⌧ u@ x (p i @ x u) = @ x (u( 1 2 ⇢ i u 2 + p i ) + p i @ x u + q i n i uE + 4 3 ⌧ u@ x (p i @ x u).
Hence, by using the viscous Ohm's law (2.26), we get the relation

un i q i E = uc i @ x p e + uc e @ x p i + 4 3 c i ⌧ u@ x (p e @ x u) 4 3 c e ⌧ u@ x (p i @ x u). Therefore @ t (⇢ i " i + 1 2 ⇢ i u 2 ) + @ x (u(⇢ i " i + 1 2 ⇢ i u 2 + p i ) = @ t (⇢ i " i ) + @ x (u⇢ i " i ) + p i @ x u uc i @ x p e + uc e @ x p i + 4 3 c i u⌧ @ x (p@ x u).
Hence from the relation

@ x (up i @ x u) = u@ x (p i @ x u) + p i (@ x u) 2 ,
we obtain the equation (2.34).

We firstly remark that the viscous system of ([3], [START_REF] Chalons | Navier-Stokes equations with several independent pressure laws and explicit predictor-corrector schemes[END_REF]) contains the terms p ↵ @ x u and 4 3 ⌧ p ↵ (@ x u) 2 , for ↵ 2 {e; i}. However, the terms ⌧ 4 3 @

x (up ↵ @ x u) and 5

2 ⌧ @ x ( k B m i p i @ x T i ) do not appear in ([3], [ 6
]).

Dissipativity of the second order terms with respect to the entropy

This section is devoted to the proof of the entropy dissipativity of the viscous system (2.22-2.25).

Proposition 4. We assume that e = i = 5/3. Let U ⌧ be a solution of the second order system (2.22-2.25). Then U ⌧ satisfies the following entropy inequality:

@ t ⌘(U ⌧ ) + @ x (u ⌧ ⌘(U ⌧ ))  ⌫ ei k B T ⌧ i T ⌧ e (T ⌧ i T ⌧ e ) 2 ⌧ 5k B 2 X ↵=e,i 1 m ↵ @ x (n ⌧ ↵ @ x T ⌧ ↵ ) (3.35)
where ⌘ is defined by (1.4).

Proof. The result is obtained by multiplying (2.22)-(2.25) by ⌘ 0 (U ⌧ ). The system (2.22)-(2.25) being written in the synthetic form (1.3), we denote W the viscous terms

W = u ⌧ @ x (J(U ⌧ )@ x U ⌧ ) + @ x (D(U ⌧ )@ x U ⌧ ) .
In [START_REF] Aregba | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF] we have shown that

⌘ 0 (U ⌧ ) [@ t U ⌧ + A(U ⌧ )@ x U ⌧ S(U ⌧ )] = ⌫ ei k B T ⌧ i T ⌧ e (T ⌧ i T ⌧ e ) 2 .
It remains to prove that

⌘ 0 (U ⌧ )W  5k B 2 X ↵=e,i 1 m ↵ @ x (n ⌧ ↵ @ x T ⌧ ↵ ) .
Note that the first component of W is equal to zero, so that @ ⇢ ⌘ is not needed. We have

⌘ 0 (U ) = (@ 1 ⌘ e (⇢c e , " e ), @ 2 ⌘ e (⇢c e , " e )) 0 @ c e 0 0 0 E e c e ⇢ 2 + q 2 ⇢ 3 q ⇢ 2 1 ⇢c e 0 1 A +(@ 1 ⌘ i (⇢c i , " i ), @ 2 ⌘ i (⇢c i , " i )) 0 @ c i 0 0 0 E i c i ⇢ 2 + q 2 ⇢ 3 q ⇢ 2 0 1 ⇢c i 1 A .
Hence, denoting q = ⇢u, @ q ⌘(U ) =

T i @ x ✓ k B m i p i @ x T i ◆ = 4n 3 (@ x u) 2 X ↵=e,i ✓ 5 2m ↵ T ↵ @ x (n ↵ k B T ↵ @ x T ↵ ) ◆ .
Using the fact that

T 1 ↵ @ x (n ↵ T ↵ @ x T ↵ ) = @ x (n ↵ @ x T ↵ ) + n ↵ (@ x T ↵ ) 2 T ↵ ,
we thus have

⌘ 0 (U )W = 4n 3 (@ x u) 2 5k B 2 X ↵=e,i 1 m ↵ @ x (n ↵ @ x T ↵ ) + n ↵ (@ x T ↵ ) 2 T ↵ ! .

Conclusion

In this paper, starting from a kinetic model, we have derived a viscous approximation of the bitemperature Euler system from a Chapman-Enskog expansion. We have been able to compute explicitly all the viscous terms and we have obtained a generalization of the model proposed in [START_REF] Chalons | Navier-Stokes equations with several independent pressure laws and explicit predictor-corrector schemes[END_REF]. Then we have proved an entropy inequality. These results support the approach taken in [START_REF] Aregba | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF] where the same kinetic model is the basis of the numerical approximation of the system (1.2). The case e 6 = i can be handled by using a kinetic model with internal energy variable (see Remark 1). In a future work we plan to study the shocks obtained by limits of travelling waves constructed from this viscous model and to compare them to the ones numerically computed in our previous work [START_REF] Aregba | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF].

,

  ↵ = e, i. Now we multiply ⌘ 0 (U ) by W , that is :

  .18)2.4 Obtention of the viscous fluid systemWe expand f ↵ as in (2.15, 2.16, 2.17) and we extract the moments w.r.t. 1, v, v 2 . One important point to determine the viscous terms of the Navier-Stokes is to compute the term g ↵ of the expansion (2.15, 2.16, 2.17)). The calculus is performed in the following proposition.

	Proposition 1. The first order terms g	e and g	i of the expansion (2.15, 2.16, 2.17)) write

  Next, we use the generalized Ohm's law (2.26) that has been established previously. So neglecting the first order terms in ⌧ in the expression of E, it comes that Z

	R 3	m e 2	(v 1 u)(v u) 2	g	e dv =		5(	k m B e	T	e ) 2	@	x n	e +	15 2	n	e	k m B e	k	B @	x T	e	35 2	k	B	k m B e	T	e @	x T	e
													+ 5(	k	B	T	e )c	e @	x (p	e + p	i ) + 5T	e	k	B	(c	i @	x p	e	c	e @	x p	i ).
																			m e								m	e
	So
					Z	R 3	m 2 e	(v 1 u)(v u) 2	g	e dv =	5 2	k m B e	p	e @	x T	e
	and finally by using the relation (2.28), it holds that
					Z	R 3	m e 2	v 1 v 2	g	e dv =				4 3	up	e @	x u	5 2	k m B e	p	e @	x T	e .
	So we obtain (2.24). The same result holds for (2.25)
	2.5 Navier-Stokes viscous system written with internal energy
	Next, in order to compare with ([3], [6]), the system (2.22, 2.23, 2.24, 2.25) is rewritten by using
	internal energies.
	Proposition 3. The system (2.22, 2.23, 2.24, 2.25) is equivalent to the system
																													@	t ⇢ + @	x (⇢u) = 0,	(2.31)
													@	t (⇢u) + @	x (⇢u 2 + p	e + p	i )	4 3	⌧ @	x (p @	x u) = 0,	(2.32)
			@	t (⇢ e "	e ) + @	x (u(⇢ e "	e )) + p	e @	x u +	4 3	⌧ @	x (up	e @	x u)	5 2	⌧ @	x (	m k B e	p	e @	x T	e )
																													4 3	⌧ p	e (@	x u) 2 = ⌫	ei (T	e	T	i ),	(2.33)


Stéphane Brull
email: stephane.brull@math.u-bordeaux.fr

A KINETIC APPROACH OF THE BI-TEMPERATURE EULER MODEL

We are interested in the numerical approximation of the bi-temperature Euler equations, which is a non conservative hyperbolic system introduced in [3]. We consider a conservative underlying kinetic model, the Vlasov-BGK-Poisson system. We perform a scaling on this system in order to obtain its hydrodynamic limit. We present a deterministic numerical method to approximate this kinetic system. The method is shown to be Asymptotic-Preserving in the hydrodynamic limit, which means that any stability condition of the method is independant of any parameter Á, with Á ae 0. We prove that the method is, under appropriate choices, consistant with the solution for bi-temperature Euler. Finally, our method is compared to methods for the fluid model (HLL, Suliciu).

Introduction

Plasma physics is renowned for being a multiscale problem ( [START_REF] Chen | Introduction to Plasma Physics and Controlled Fusion[END_REF]). Phenomena range from plasma oscillations and charge separation at very small scales to fluid mechanics at macroscopic levels. Di erent physical models exist in order to describe a plasma according to the scale of interest. The kinetic model, constituted of the Vlasov-Boltzmann equations, coupled with Maxwell equations, are well suited for mesoscopic descriptions. In this model, each population of particles -(in the case of a plasma, these are electrons, ions and neutral atoms and molecules) is represented via the use of a distribution function f -. This function, defined over the phase space, denotes the number of particles that, at a given time t at a given point x, possess a velocity v. Via an hydrodynamic limit, equations on the macroscopic conserved quantities (mass, momentum and energy) can be derived. These equations form the so-called bi-temperature Euler system (see [START_REF] Aregba-Driollet | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF] for the electrostatic case, and [START_REF] Brull | Modelling and entropy satisfying relaxation scheme for the nonconservative bitemperature Euler system with transverse magnetic field[END_REF] for the transverse magnetic case). From a numerical point of view, these equations are less accurate than the underlying kinetic model, but they are cheaper to solve.

The bi-temperature Euler system is comprised of an equation over mass, an equation over momentum, and one equation over each species energy. These energy equations are non-conservative, accounting for the fact that, in a plasma, electrons and ions are exchanging energy through collisions and via electromagnetic interactions. Contrarily to the usual Euler system, for which the exact solution of any Riemann problem is known, these non-conservative terms prevent the knowledge of such solutions, especially when the solution contains a shock. This is a consequence that these non-conservative terms are written as products of distributions, which does not possess a clear definition. [START_REF] Maso | Definition and weak stability of nonconservative products[END_REF] presents a way of defining these products using a theory of paths. It is made clear that a certain amount of arbitrariness is needed in the choice of such paths. A major drawback of this property is that, when attempting numerical resolution of this system, every numerical scheme chooses its own path, and the Rankine-Hugoniot relation holding across shocks will be di erent for every scheme. This phenomenon is highlighted in [START_REF] Aregba-Driollet | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF]. More precisly, di erent schemes have di erent numerical viscosity, and the shock relations depend on such viscosities. Several solutions have been proposed in order to derive numerical methods that use the notion of paths to define jump relations ( [START_REF] Parès | Path-conservative numerical methods for nonconservative hyperbolic systems[END_REF], [START_REF] Chalons | A new comment on the computation of non-conservative products using Roe-type path conservative schemes[END_REF]). Other approaches have been proposed for di erent kinds of non-conservative systems (see [START_REF] Chalons | Large Time Step and Asymptotic Preserving Numerical Schemes for the Gas Dynamics Equations with Source Terms[END_REF], [START_REF] Dumbser | A new e cient formulation of the HLLEM Riemann solver for general conservative and non-conservative hyperbolic systems[END_REF], [START_REF] Lefloch | Numerical methods with controlled dissipation for small-scale dependent shocks[END_REF], [START_REF] Xu | Multiple temperature kinetic model and gas-kinetic method for hypersonic nonequilibrium flow computations[END_REF], [START_REF] Abgrall | A comment on the computation of non-conservative products[END_REF]). Finally, Bouchut's book [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF] presents an overview of numerical methods for non-conservative Shallow-water equations.

In this paper, the aim is to propose a reference numerical method for such solutions using the underlying kinetic model. This kinetic model is conservative, and hence does not exhibit the drawback of the macroscopic model. Hence, the idea dwells in solving the Vlasov-BGK-Ampère system in the hydrodynamic limit, which will be presented in section 2. [START_REF] Aregba-Driollet | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF]. In order to be able to compare the results with the scheme applied to the bi-temperature Euler system, it is necessary to describe identical scales. By performing a scaling on the kinetic system, dimensionless parameters (scaled Debye length, scaled plasma frequency) are introduced in the system. Taking the hydrodynamic limit then amounts to taking the limit when these parameters tend to zero. In the general case, if a naive numerical approach is used, extremely restrictive stability conditions will appear, rendering the scheme unusable in a decent amount of computation time. Hence, an Asymptotic-Preserving (AP) scheme needs to be derived. Such a scheme possess stability conditions independent of these small parameters and is then able to compute solutions for both cases when parameters are of the order of one and when these parameters are arbitrarily small. The adjective Asymptotic Preserving has been introduced by Shi Jin in [START_REF] Jin | E cient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF]. In the study of plasma simulation, numerous several Asymptotic Preserving methods have been proposed for various models, such as Euler-Maxwell, Euler-Poisson, Euler-Lorentz, Vlasov-Poisson ( [START_REF] Degond | Numerical approximation of the Euler-Poisson-Boltzmann model in the quasineutral limit[END_REF], [START_REF] Brull | Asymptotic-preserving scheme for a bi-fluid Euler-Lorentz model[END_REF], [START_REF] Degond | Numerical approximation of the Euler-Maxwell model in the quasineutral limit[END_REF], [START_REF] Crispel | An asymptotic preserving scheme for the two-fluid Euler-Poisson model in the quasineutral limit[END_REF], [START_REF] Crispel | An asymptotically stable discretization for the Euler-Poisson system in the quasi-neutral limit[END_REF], [START_REF] Degond | Asymptotic-preserving methods and multiscale models for plasma physics[END_REF], [START_REF] Belaouar | An asymptotically stable semi-Lagrangian scheme in the quasi-neutral limit[END_REF]).

The paper is organized as follows. Firstly, the macroscopic and kinetic models relative to our problem are introduced. A reduced 1D3V kinetic model is derived, and the relative properties of moment conservation and entropy dissipation are specified. Subsequently, a scaling of the kinetic model in order to exhibit the limiting behaviour of the system. Then, a Asymptotic-Preserving numerical method is presented. This method is able to compute solutions of the kinetic model in the asymptotic regime. Finally, test cases are presented in order to illustrate the derived method. In particular, the numerical results are compared with an HLL-inspired method and a Suliciu relaxation approach, both derived in [START_REF] Aregba-Driollet | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF].

Description of the models

2.1. Bi-temperature Euler model. The one-dimensional bi-temperature Euler model describes the behaviour of a two species fluid (constituted of electrons and ions) on a macroscopic scale (see [START_REF] Coquel | Numerical methods for weakly ionized gas[END_REF]):

Y _ _ _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ _ _ [ ˆtfl + ˆx(flu) = 0, (1a) ˆt(flu) + ˆx(flu 2 + p e + p i ) = 0, (1b) ˆt(fl i ' i + 1 2 fl i u 2 ) + ˆx(u(fl i ' i + 1 2 fl i u 2 + p i )) + u(c i ˆxp e ≠ c e ˆxp i ) = ≠‹ ei (T i ≠ T e ), (1c)
ˆt(fl e ' e + 1 2 fl e u 2 ) + ˆx(u(fl e ' e + 1 2

fl e u 2 + p e )) ≠ u(c i ˆxp e ≠ c e ˆxp i ) = ‹ ei (T i ≠ T e ). (1d)
Superscripts e and i denote quantities related to electrons and ions, respectively. fl = fl e + fl i = m e n e + m i n i is the total density, m -being the mass of a particle of species -, and n -the corresponding concentration. u is the macroscopic velocity of the fluid. p -, T -and ' -denote the partial pressures, temperatures and specific energies. c -= m - m -+m -denotes the mass fractions. ‹ ei is the temperature exchange rate between the two species, that can depend on space and time. 2.2. Macroscopic quantities. Define f -(t, x, v) : R + ◊ R 3 ◊ R 3 ' ae R + , the particle distribution function of species -, where t, x, v respectively denote the time, space and microscopic velocity variables. Integration with respect to v oe R 3 will be denoted as follows, for any function g that depends on v such that

(1 + v 2 )g oe L 1 (R 3 ): ÈÈgÍÍ = ⁄ R 3 gdv. Macroscopic quantities n -(t, x), u -(t, x)
and T -(t, x) are defined as moments of f -as follows:

ÈÈf -ÍÍ = n -, ÈÈvf -ÍÍ = n -u -, ÈÈm -v 2 2 f -ÍÍ = 1 2 m -n -(u -) 2 + 3 2 n -k B T -.
Similarly, mixture macroscopic quantities are defined as:

ÈÈm e f e + m i f i ÍÍ = fl, ÈÈv(m e f e + m i f i )ÍÍ = flu, ÈÈ v 2 2 (m e f e + m i f i )ÍÍ = 1 2 flu 2 + 3 2 (n e + n i )k B T.
Finally, we define electromagnetic quantities: E oe R is the electric field along x, q -is the electric charge of species -. Note that q e = ≠Zq i , with Z the ionization rate. j and fl are respectively the current and total electric charge, defined as:

ÈÈq e f e + q i f i ÍÍ = fl, ÈÈv(q e f e + q i f i )ÍÍ = j. (2)
2.3. Vlasov-BGK-Ampère system. Consider the Vlasov-BGK-Ampère model, for

t oe R, x oe R, v = (v 1 , v 2 , v 3 ) oe R 3 and -oe {e, i}: Y _ _ _ _ _ _ ] _ _ _ _ _ _ [ ˆtf -+ v 1 ˆxf -+ q -E m -ˆv1 f -= 1 • -(M (f -) ≠ f -) + 1 • --(M (f -) ≠ f -), (3a) ˆtE = ≠ j Á 0 , (3b) ˆxE = fl Á 0 , ( 3c
)
where • -and • --are, respectively, the relaxation rate towards intra-species and inter-species equilibria (• --= • --). Á 0 is the dielectric permittivity of vacuum.

The entropy-minimizing distribution function is the local Maxwellian distribution, denoted M (f -), defined by:

M (f -)(t, x, v) = n -(t, x) (2fik B T -(t,x) m -) 3 2 exp (≠ |v ≠ u -(t, x)| 2 2k B T -(t,x) m -
).

M (f -), the exchange Maxwellian distribution, is defined as:

M (f -)(t, x, v) = n -(t, x) (2fik B T (t,x) m -) 3 2 exp (≠ |v ≠ u(t, x)| 2 2k B T (t,x) m -
).

2.4. Reduced model.

Definition of the reduced distribution functions.

The problem is assumed to be one dimensional in space and 3 dimensional in velocity (1D3V). Hence, the distribution functions f -only depend on x = x 1 and v = (v 1 , v 2 , v 3 ). This means that the function f -is even in v 2 and v 3 , which can be expressed as:

Y _ _ _ _ ] _ _ _ _ [ ⁄ R 2 f -(t, x, v)v 2 dv 2 dv 3 = 0, (4a) ⁄ R 2 f -(t, x, v)v 3 dv 2 dv 3 = 0. (4b) Denote u -= (u - 1 , u - 2 , u - 3 ), this assumption gives u -= (u - 1 , 0, 0). Hence, the notation u -= u - 1 is used.
The number of variables is reduced, as in [START_REF] Andries | The Gaussian-BGK model of Boltzmann equation with small Prandtl number[END_REF]. Define the marginals, for -= e, i, by:

Y _ _ _ _ _ ] _ _ _ _ _ [ g - 0 (t, x, v 1 ) = ⁄ R 2 f -(t, x, v)dv 2 dv 3 , (5a) g - 2 (t, x, v 1 ) = ⁄ R 2 f -(t, x, v) m -(v 2 2 + v 2 3 ) 2 dv 2 dv 3 , (5b)
and denote integration with respect to v 1 , for any function g that depends on v 1 , as:

(6) ÈgÍ = ⁄ R gdv 1 .
2.4.2. Reduced system of equations and its properties. The idea is to describe the system of equations and the macroscopic quantities using these reduced functions. The system obtained after reduction reads as follows:

Proposition 1. Under the assumptions ( 4) and ( 5), the system (3) can be reduced to the following one, with -= e, i and p = 0, 2:

Y _ _ _ _ _ _ ] _ _ _ _ _ _ [ ˆtg - p + v 1 ˆxg - p + q -E m -ˆv1 g - p = 1 • -(M - p ≠ g - p ) + 1 • --(M - p ≠ g - p ), ( 7a
)
ˆxE| t=t0 = 1 Á 0 fl| t=t0 , (7b) ˆtE = ≠ 1 Á 0 j, ( 7c
)
where the di erent Maxwellian distributions are defined as:

Y _ _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ _ [ M - 0 (t, x, v 1 ) = n -(t, x) Ò 2fik B T -(t,x) m - exp (≠ (v 1 ≠ u -) 2 2k B T -(t,x) m - ), (8a) M - 2 (t, x, v 1 ) = k B T -(t, x)M - 0 (t, x, v 1 ), (8b) M - 0 (t, x, v 1 ) = n -(t, x) Ò 2fik B T (t,x) m - exp (≠ (v 1 ≠ u) 2 2k B T (t,x) m - ), (8c) M - 2 (t, x, v 1 ) = k B T (t, x)M - 0 (t, x, v 1 ). (8d)
Moreover, g - 0 and g - 2 satisfy the following relations:

Y _ _ ] _ _ [ n -= Èg - 0 Í, (9a) n -u -= Èg - 0 v 1 Í, (9b) 3 2 n -k B T -+ 1 2 m -n -(u -) 2 = Èg - 0 m -v 2 1 2 Í + Èg - 2 Í. (9c)
Proof. Firstly, by integrating equation (3a) with respect to v 2 and v 3 , it comes: [START_REF] Brull | Local discrete velocity grids for deterministic rarefied flow simulations[END_REF]. The definition of g - 0 given in (5a) allows us to rewrite [START_REF] Chalons | Large Time Step and Asymptotic Preserving Numerical Schemes for the Gas Dynamics Equations with Source Terms[END_REF] as (7a) with p = 0. We have obtained an equation on the function g - 0 . Let us express the moments using g - 0 . The 0 th -order moment can be expressed as:

ˆt ⁄ R 2 f -dv 2 dv 3 + v 1 ˆx ⁄ R 2 f -dv 2 dv 3 = 1 • -( ⁄ R 2 M (f -)dv 2 dv 3 ≠ ⁄ R 2 f -dv 2 dv 3 ) + 1 • --( ⁄ R 2 M (f -)dv 2 dv 3 ≠ ⁄ R 2 f -dv 2 dv 3 ). (10) Define M - 0 , M - 0 : R + t ◊ R x ◊ R v ' ae R using
ÈÈf -ÍÍ = Èg - 0 Í = n -.
Similarly, the first order moment writes:

n -u -= ÈÈf -vÍÍ = K ⁄ R 2 f - Q a v 1 v 2 v 3 R b dv 2 dv 3 L .
According to (4a) and (4b), the second and third components of the first order moment are equal to zero. Then, we get, for the first component:

K v 1 ⁄ R 2 f -dv 2 dv 3 L = Èv 1 g - 0 Í = n -u -.
The third moments reads as:

ÈÈf -m -|v| 2 2 ÍÍ = = g - 0 m -v 2 1 2 > + ÈÈf -m -(v 2 2 + v 2 3 ) 2 ÍÍ. (11)
Hence, another function g - 2 , defined by (5b), is introduced. Rewrite (11) as:

3 2 n -k B T -+ 1 2 m -n -(u -) 2 = = g - 0 m -v 2 1 2 > + Èg - 2 Í .
Finally, in order to close the system, an equation over g - 2 is needed. To obtain it, multiply (3a) by m -(v 2 2 +v 2 3 ) 2

, and integrate with respect to v 2 and v 3 :

⁄ R 2 ˆtf -m -(v 2 2 + v 2 3 ) 2 dv 2 dv 3 + ⁄ R 2 v 1 ˆxf -m -(v 2 2 + v 2 3 ) 2 dv 2 dv 3 + q -E m - ⁄ R 2 ˆv1 f -m -(v 2 2 + v 2 3 ) 2 dv 2 dv 3 = 1 • -( ⁄ R 2 M (f -) m -(v 2 2 + v 2 3 ) 2 dv 2 dv 3 ≠ ⁄ R 2 f m -(v 2 2 + v 2 3 ) 2 dv 2 dv 3 ).
Then we introduce the function M - 2 of (8b) corresponding to:

M - 2 = ⁄ R 2 M (f -) m -(v 2 2 + v 2 3 ) 2 dv 2 dv 3 .
We get equation (7a) with p = 2. ⇤ Remark 1. Equation (7b) means that, provided that Maxwell-Gauss equation is satisfied at t = 0 and if the Maxwell-Ampère equation is satisfied for all t, so will the Maxwell-Gauss equation.

Regarding the conservation properties of the BGK operators in equations (7a), the following proposition holds: Proposition 2. The intra-species BGK operators defined by the Maxwellians (8a) and (8b) conserve the partial masses, momenta and energies according to the following identities, for -= e, i:

(12) = (M - 0 ≠ g - 0 ) 3 1 v 1 4> = 0, = (M - 0 ≠ g - 0 ) m -(v 1 ) 2 2 > + È(M - 2 ≠ g - 2 )Í = 0.
The exchange Maxwellians (8c) and (8d) define inter-species BGK operators that enjoy the following conservation properties:

e M 0 -≠ g - 0 f = 0, ÿ - e (M 0 -≠ g - 0 )m -v 1 f = 0, ÿ - = (M 0 -≠ g - 0 ) m -(v 1 ) 2 2 > + ÿ - e M 2 -≠ g - 2 f = 0.
In order to obtain the entropy of the reduced system, consider the following definition of the Maxwellian distribution of the non-reduced system:

H -(M (f -)) = min f -Ø0 (H -(f -)s.t.g - 0 (v 1 ) = ⁄ R 2 f -(v)dv 2 dv 3 , g - 2 (v 1 ) = ⁄ R 2 f -(v) m -(v 2 2 + v 2 3 ) 2 dv 2 dv 3 ), with (13) H -(f -) = ÈÈf -log f -ÍÍ.
The resolution of this constrained minimization problem leads to the following expression of M (f -):

(14) M (f -)(g - 0 , g - 2 ) = (g - 0 ) 2 m - 2fig - 2 exp (≠ g - 0 m -(v 2 2 + v 2 3 ) 2g - 2 ).
Hence, by injecting ( 14) in [START_REF] Crispel | An asymptotically stable discretization for the Euler-Poisson system in the quasi-neutral limit[END_REF], it comes

H -(g - 0 , g - 2 ) = = g - 0 log (g - 0 ) 2 g - 2 ≠ 2g - 0 > .
Hence, define the following notations for the partial entropies of the reduced model:

(15) H -(g - 0 , g - 2 ) = È÷ -(g - 0 , g - 2 )Í , with ÷ -(g - 0 , g - 2 ) = g - 0 log (g - 0 ) 2 g - 2 ≠ 2g - 0 ,
for -= e, i, and the total entropy ( 16)

H(g e 0 , g e 2 , g i 0 , g i 2 ) = ÿ - H -(g - 0 , g - 2 ).
Firstly, this entropy is shown to be convex: Proposition 3. For -= e, i, H -is a strictly convex function. Consequently, H is strictly convex as well.

Proof. Compute the gradient of ÷ -:

(17)

Ò÷ -= 3 ˆg-0 ÷ - ˆg-2 ÷ - 4 = A log (g - 0 ) 2 g - 2 ≠ g - 0 g - 2 B
.

Compute the Hessian:

(18) Ò 2 ÷ -= A 2 g - 0 ≠ 1 g - 0 ≠ 1 g - 0 g - 0 (g - 2 ) 2 B .
By positivity of the distribution functions,

Tr(Ò 2 ÷ -) > 0 and det (Ò 2 ÷ -) = 1 (g - 2 ) 2 > 0.
Hence, the Hessian is positive definite and ÷ -is strictly convex. Consequently, H -is strictly convex. ⇤

The following property on entropy dissipation holds:

Proposition 4. The reduced system [START_REF] Brull | Modelling and entropy satisfying relaxation scheme for the nonconservative bitemperature Euler system with transverse magnetic field[END_REF] satisfies the H-theorem for the entropy defined by ( 16):

(19) ˆH ˆt + ˆG ˆx AE 0, with G the entropy flux associated to H: G(g e 0 , g e 2 , g i 0 , g i 2 ) = ÿ - G -(g - 0 , g - 2 ) = ÿ - Èv 1 ÷ -(g - 0 , g - 2 )Í . Proof. Multiply equation (7a) by ˆg- p ÷ -(g - 0 , g - 2 )
, sum over p and integrate over v 1 :

ˆH- ˆt + ˆG- ˆx = 1 • - + (M - 0 ≠ g - 0 )ˆg-0 ÷ -+ (M - 2 ≠ g - 2 )ˆg-2 ÷ -, + 1 • -- e (M - 0 ≠ g - 0 )ˆg-0 ÷ -+ (M - 2 ≠ g - 2 )ˆg-2 ÷ - f .
Then, summing overgives:

ˆH ˆt + ˆG ˆx = ÿ - 1 • - + (M - 0 ≠ g - 0 )ˆg-0 ÷ -+ (M - 2 ≠ g - 2 )ˆg-2 ÷ -, + ÿ - 1 • -- e (M - 0 ≠ g - 0 )ˆg-0 ÷ -+ (M - 2 ≠ g - 2 )ˆg-2 ÷ - f . ( 20
)
Consider every term of the right-hand side consecutively.

+ (M - 0 ≠ g - 0 )ˆg-0 ÷ -+ (M - 2 ≠ g - 2 )ˆg-2 ÷ -, = + (M - 0 ≠ g - 0 )(ˆg-0 ÷ -(g - 0 , g - 2 ) ≠ ˆg-0 ÷ -(M - 0 , M - 2 )) , + + (M - 2 ≠ g - 2 )(ˆg-2 ÷ -(g - 0 , g - 2 ) ≠ ˆg-2 ÷ -(M - 0 , M - 2 )) , + + (M - 0 ≠ g - 0 )ˆg-0 ÷ -(M - 0 , M - 2 ) , + + (M - 2 ≠ g - 2 )ˆg-2 ÷ -(M - 0 , M - 2 ) , . (21)
By convexity of ÷ -, the sum of the first two terms of the right-hand side of ( 21) is negative. The remaining terms are

+ (M - 0 ≠ g - 0 )ˆg-0 ÷ -(M - 0 , M - 2 ) , + + (M - 2 ≠ g - 2 )ˆg-2 ÷ -(M - 0 , M - 2 ) , = = (M - 0 ≠ g - 0 ) log ( (M - 0 ) 2 M - 2 ) > ≠ = (M - 2 ≠ g - 2 ) M - 0 M - 2 > = 1 k B T -È(M - 0 ≠ g - 0 ) log (M - 0 )Í ≠ 1 k B T -È(M - 2 ≠ g - 2 )Í = 1 k B T -log ( n - (2fik B T - m -) 1 2
)

È(M - 0 ≠ g - 0 )Í + K (M - 0 ≠ g - 0 )( vu - k B T - m - ≠ (u -) 2 2k B T - m - ) L ≠ 1 k B T - = (M - 0 ≠ g - 0 ) m -(v -) 2 2 + (M - 2 ≠ g - 2 ) > = 0,
according to [START_REF] Coquel | Numerical methods for weakly ionized gas[END_REF]. By an analogous computation on the remaining last term in [START_REF] Guisset | Asymptotic-preserving scheme for the M 1 -Maxwell system in the quasi-neutral regime[END_REF], ( 19) is obtained. ⇤ 2.5. Scaling of the equations. The regime of interest is the hydrodynamic limit of this model, which corresponds to the scale of the phenomena described by the bi-temperature Euler equations (see [START_REF] Aregba-Driollet | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF]). To do so, it is necessary to perform a scaling on this model, in order to exhibit the behaviour of the equations when the macroscopic characteristic scale is much larger than the ones involved by the kinetic equations (collisional scale, Debye scale). Denote t 0 , L, v th , q, m respectively the reference time, length, velocity, electric charge and particle mass. Denoting with a tilda the dimensionless quantities, it comes:

 t = t t 0 ,  x = x L ,  v 1 = v 1 v th , Ê q -= q - q , Á m -= m - m .
From this scaling, define also the following dimensionless quantities:

 g - p = v th n 0 g - p ,  E = qt 0 v th m E,  j = j qv th n 0 ,  f = f qn 0 ,  T = k B mv 2 th T.
where n 0 is the reference density.

Applying this scaling to our model, and removing the tildas for the sake of readability, it comes, for p = 0, 2 and -= e, i:

Y _ _ _ _ _ _ ] _ _ _ _ _ _ [ ˆtg - p + v 1 ˆxg - p + q - m -Eˆv 1 g - p = 1 Á e (M e p ≠ g e p ) + 1 Á ei (M e p ≠ g e p ), (22a) ˆtE = ≠ j -2 , (22b) ˆxE = fl -2 , ( 22c
)
where

Á -= • - t0 , Á --= • -- t0 and -= ⁄ De L , ⁄ De = ( Á0kB T0 q 2 n0 ) 1
2 , with T 0 = mv 2 th kB the reference temperature. ⁄ De is the so-called Debye length, which is the characteristic length of charge separation in a plasma. In the regime of interest, this quantity is extremely small compared to the reference length of the domain L. Similarly, the same-species mean free path related quantities Á -are small compared to the reference length L. Hence, reaching the hydrodynamic limit consists in taking the limit of the model when Á -andtend to 0. Notice that, however, Á --is finite. In the limit, the Maxwell equations (22b-22c) become j = 0 and fl = 0, which can be rewritten as:

u e = u i , n e = Zn i . ( 23
)
The equations ( 23) are called the quasi-neutrality constraints. This form of the system of equations is the singular quasi-neutral limit of our initial model, where two evolution equations degenerate into algebraic relations. The purpose of this work is to derive an Asymptotic-preserving (AP) numerical scheme, that is to say with a stability condition that is independant of every non-dimensional small parameters Á -, Á --and -.

2.6. Entropy inequality for hydrodynamic limits. In [START_REF] Aregba-Driollet | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF], the existence of an entropy-flux pair is shown for the bi-temperature Euler system [START_REF] Abgrall | A comment on the computation of non-conservative products[END_REF]. In this section, an entropy dissipation inequality in the limit is derived for the reduced model presented in Proposition 1.

Firstly, denote U = (fl, flu, ' e , ' i ) the unknowns of the bi-temperature model and (H macro (U), G macro (U)) the corresponding macroscopic entropy-flux pair. Note the following properties regarding the link between kinetic and macroscopic entropies: [START_REF] Xu | Multiple temperature kinetic model and gas-kinetic method for hypersonic nonequilibrium flow computations[END_REF] H(M e 0 , M e 2 , M i 0 , M i 2 ) = H macro (U), G(M e 0 , M e 2 , M i 0 , M i 2 ) = G macro (U). See [START_REF] Aregba-Driollet | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF] for further details on these properties. Proposition 5. Let U be a limit solution of the kinetic model from Proposition 1, i.e moments of solution distribution functions (g - 0 , g - 2 ) in the hydrodynamic limit. Then, the following inequality on entropy holds:

(25) ˆHmacro ˆt + ˆGmacro ˆx AE ≠ ‹ ei k B T i T e (T i ≠ T e ) 2 , with ‹ --= 3k B n -n - 2• --(n -+ n -) . Proof. Multiply equation (7a) by ˆg- p ÷ -(g - 0 , g - 2 )
, sum over p and integrate over v 1 . It comes, for -= e, i:

ˆH- ˆt + ˆG- ˆx = 1 Á - + (M - 0 ≠ g - 0 )ˆg-0 ÷ -+ (M - 2 ≠ g - 2 )ˆg-2 ÷ -, + 1 Á -- e (M - 0 ≠ g - 0 )ˆg-0 ÷ -+ (M - 2 ≠ g - 2 )ˆg-2 ÷ - f .
Using the entropy dissipation properties of the BGK operator shown in Section 2.4.2:

ˆH- ˆt + ˆG- ˆx AE 1 Á -- e (M - 0 ≠ g - 0 )ˆg-0 ÷ -+ (M - 2 ≠ g - 2 )ˆg-2 ÷ - f .

Vlasov-BGK equations: general method.

A time-splitting method is used on the Vlasov-BGK equations, in order to separate the Vlasov equations from the BGK operators. More precisely, the following equations are going to be consecutively solved:

(26)

ˆtg - p + v 1 ˆxg - p +
q - m -Eˆv 1 g - p = 0, the Vlasov-Maxwell equations for -= e, i, coupled with the Ampère equation:

(27) ˆtE = ≠ j • 2 .
Then, the e ects of the inter-species collisions are computed:

(28)

ˆtg - p = 1 Á --(M - p ≠ g - p )
, and finally, intra-species collisions:

(29) ˆtg - p = 1 • (M - p ≠ g - p )
, where • ae 0. A uniform discretization of the phase space is considered. Concerning the velocity discretization, in the present case, only static grids are considered, contrarily to [START_REF] Brull | Local discrete velocity grids for deterministic rarefied flow simulations[END_REF]. A suitable subset of R that contains the most part of the distribution function has to be chosen in order to perform correct numerical evaluations of the moments.

However, in this model, two particle distributions are considered, in order to follow the dynamics of both electrons and ions via a kinetic description. The major di erence between these two species is the value of the mass of a particle that can have a ratio of the order of one thousand. Consequently, the distribution functions have very di erent supports. Thus, if only one velocity grid is employed, then this grid has to be large enough to encompass the largest distribution and has to be fine enough to describe the smallest one, which is numerically very costly. Hence, two velocity grids, one for each species, are considered. The boundaries of the domains are chosen as follows: for a given initial condition, let us denote the initial moments by n - 0 , u - 0 and T - 0 . For -= e, i, we have:

V -= [min x (u - 0 (x)) ≠ l -max x Ú k B T - 0 (x) m -, max x (u - 0 (x)) + l -max x Ú k B T - 0 (x) m -], (30)
where l -is chosen to fit the initial maxwellian distribution. Note that this grid is the same for all points in space and is constant in time.

Notations for the discretization are chosen as follows: x denotes the discretization step of the domain [L min , L max ] µ R in N

x cells. For -= e, i, the velocity domains [V - min , V - max ] µ R are both discretized using N v cells and the discretization steps are denoted v e and v i . Then, we use indices k and j to denote quantities computed at the k-th cell in space and j-th cell in velocity. Moreover, we denote t the time step and we use a superscript n to denote the quantity computed at time t n = t 0 + n t.

Step 1: Computation of g and E

The following discretization is performed, for -= e, i:

(31) Â g -,n+1 0,k,j = g -,n 0,k,j ≠ t x (" -,n 0,k+ 1 2 ,j ≠ " -,n 0,k≠ 1 2 ,j ) ≠ t 2 v - q -E n+1 k m -(g -,n 0,k,j+1 ≠ g -,n 0,k,j≠1 ), coupled with (32) E n+1 k = E n k ≠ t • j n+1 k .
Note that the term ˆv1 g 0 has been discretized using a centred scheme. For the space discretization, " 0 represented the flux function of the spatial discretization which will be specified in the next section.

Then, the following procedure is applied to obtain  g -,n+1

2

: (33) Â g -,n+1 2,k,j = g -,n 2,k,j ≠ t x (" -,n 2,k+ 1 2 ,j ≠ " -,n 2,k≠ 1 2 ,j ) ≠ t v -(Â -,n 2,k,j+ 1 2 ≠ Â -,n 2,k,j≠ 1 2 ),
where

" -,n 2,k+ 1 2 ,j = v - j 2 (g -,n 2,k+1,j + g -,n 2,k,j ) ≠ |v - j | 2 (g -,n 2,k+1,j ≠ g -,n 2,k,j ), and
 -,n 2,k,j+ 1 2 = q -E n+1 k 2m -(g -,n 2,k,j+1 + g -,n 2,k,j ) ≠ |q -E n+1 k | 2m -(g -,n
2,k,j+1 ≠ g -,n 2,k,j ). Equations ( 31) are stable under the CFL condition :

(34) t AE min - ( 1
||V -||OE x + ||q - E n+1 ||OE m - v -
).

Moments of  g e,n+1

p and  g i,n+1

p are then computed using the midpoint quadrature formula. 'k oe {1, ..., N

x }:

 n -,n+1 k = Nv ÿ j=1  g -,n+1 0,k,j v -,  u -,n+1 k = 1  n -,n+1 k Nv ÿ j=1  g -,n+1 0,k,j v j v -,  T -,n+1 k = 1 3 n -,n+1 k k B S U Nv ÿ j=1 m - g -,n+1 0,k,j v 2 j v -+ 2 Nv ÿ j=1  g -,n+1 2,k,j v -≠ m - n -,n+1 k  u -,n+1 2 k T V .
In order to compute the mixture Maxwellian distribution in the next step, mixture quantities must also be computed:

 n n+1 k =  n e,n+1 k +  n i,n+1 k ,  u n+1 k = q -  fl -,n+ 1 2 k  u -,n+1 k q -  fl -,n+1 k ,  T n+1 k = 1  n n+1 k k B C ÿ - 1 3  fl -,n+1 k (( u -,n+1 k ) 2 ≠ ( u n+1 k ) 2 ) + ÿ -  n -,n+1 k k B  T -,n+1 k D .
Step 2: inter-species relaxation

We apply a backward Euler method to the equations (28):

ĝ-,n+1 p,k,j = 1 t + • ei 5 t Ê M -,n+1 p,k,j + • ei  g -,n+1 p,k,j 6 which is unconditionnaly stable. Ê M -,n+1
p,k,j is the mixture Maxwellian defined by the moments

 n -,n+1 k ,  u n+1 k and  T n+1 k
, which are conserved through this step. Hence, the mixture Maxwellian is known beforehand, which allows this implicitation. Partial velocities and temperatures (denoted in the next section with a hat) are then recomputed using the previous formulae.

Step 3: intra-species relaxation The next step consists in the relaxation towards Maxwellian equilibrium. Since the BGK operator preserves moments, we have the following properties: 'k oe {1, ..., N

x }, (35) n-,n+1

k = n -,n+1 k , û-,n+1 k = u -,n+1 k , T -,n+1 k = T -,n+1 k .
Similarly to the previous section, this property allows us to use an implicit scheme, the backward Euler method: 'k oe {1, ..., N

x }, 'j oe {1, ..., N v } and p = 0, 2,

(36) g -,n+1 p,k,j = ĝ-,n+1 p,k,j + t • (M -,n+1 p,k,j ≠ g -,n+1
p,k,j ), where

M -,n+1 0,k,j = n -,n+1 k Ò 2fik B T -,n+1 k m - exp (≠ (v j ≠ u -,n+1 k ) 2 2k B T -,n+1 k m - ), M -,n+1 2,k,j = k B T - k M -,n+1 0,k,j ,
which are expressed using quantities known at the end of step 2, according to the properties (35). Taking the limit • ae 0 in (36) leads to:

g -,n+1
p,k,j = M -,n+1 p,k,j . In the next section, enforcement of quasineutrality constraints (23) will be adressed. More precisely, specific spatial fluxes will be given and a reformulation of the Maxwell-Ampère equation will be performed in order to derive an Asymptotic-Preserving method, valid in the limit • ae 0. Note that if these constraints are enforced through the hyperbolic step of the scheme (step 1), then they will be trivially preserved by both step 2 and 3. Hence, quasineutrality only has to been enforced through step 1.

3.2. Modified-viscosity upwind scheme. The space divergence is discretized explicitly, via the following modified upwind scheme:

(37) " -,n 0,k+ 1 2 ,j = v - j 2 (g -,n 0,k+1,j + g -,n 0,k,j ) ≠ |V max | 2 (g -,n 0,k+1,j ≠ g -,n 0,k,j ),
where

V max = max (|V e min |, |V i min |, |V e max |, |V e min |
). This method is the classical upwind scheme, where the numerical viscosity |v j |, that depends directly of the microscopic velocity, has been replaced by the smallest velocity that ensures stability for all equations. The discrete system writes as follows:

(38) Â g -,n+1 0,k,j = g -,n 0,k,j ≠ t x (" -,n 0,k+ 1 2 ,j ≠ " -,n 0,k≠ 1 2 ,j ) ≠ t 2 v - q -E n+1 k m -(g -,n 0,k,j+1 ≠ g -,n 0,k,j≠1 ), and (39)
E n+1 k = E n k ≠ t • j n+1 k .
This level of implicitness allows for a reformulation of the equation on the electric field. More precisely, the equation used to compute the electric field in (39) is not valid in the limit • ae 0.

To solve this problem, we follow the idea developped in [START_REF] Guisset | Asymptotic-preserving scheme for the M 1 -Maxwell system in the quasi-neutral regime[END_REF]. By multiplying equation (38) by q -v - j v -and summing over j, it comes:

q - ÿ j v - j  g -,n+1 0,k,j v -= q - ÿ j v - j g -,n 0,k,j v -≠ q -t v - x ÿ j v - j (" -,n 0,k+ 1 2 ,j ≠ " -,n 0,k≠ 1 2 ,j ) ≠ t 2 ÿ j v - j ( (q -) 2 E n+1 k m - (g -,n 0,k,j+1 ≠ g -,n 0,k,j≠1 )). (40)
According to the definition of the current

j n k = q - q -q j v - j g -,n 0,k,j v -, taking (40) into account with -= e, i, it comes  j n+1 k = j n k ≠ t x ÿ - q - ÿ j v -v - j (" -,n 0,k+ 1 2 ,j ≠ " -,n 0,k≠ 1 2 ,j ) ≠ tE n+1 k 2 ÿ - (q -) 2 m - ÿ j v - j (g -,n 0,k,j+1 ≠ g -,n 0,k,j≠1 ).
By injecting this equation in (39), we get the following expression for the electric field:

E n+1 k = • E n k ≠ tj n k + t 2 x q - q -q j v -v - j (" -,n 0,k+ 1 2 ,j ≠ " -,n 0,k≠ 1 2 ,j ) • ≠ t 2 2 q - (q -) 2 m - q j v - j (g -,n 0,k,j+1 ≠ g -,n 0,k,j≠1 )
, which is valid for all values of • . In particular, we are interested in the limit • ae 0, which is:

E n+1 k = 1 x q - q -q j v -v - j (" -,n 0,k+ 1 2 ,j ≠ " -,n 0,k≠ 1 2 ,j ) q - (q -) 2 m -n -,n k ,
where it has been assumed that j n k = 0. Such a choice for E n+1 k ensures that  j n+1 k = 0, which is one of the two quasi-neutrality constraints [START_REF] Parès | Path-conservative numerical methods for nonconservative hyperbolic systems[END_REF]. The choice of " 0 (37) that has been made ensures that the other contraint  fn+1 k = 0 is also verified. More precisely, compute the discrete equation on  fn+1 k . 'k oe {1, ..., N

x }, provided that ÿ

j (Â -,n 0,k,j+ 1 2 ≠ Â -,n 0,k,j≠ 1
2 ) = 0, the electric charge is

ÿ j q -g -,n+1 0,k,j v -= ÿ j q -g -,n 0,k,j v -≠ t x ÿ j (" -,n 0,k+ 1 2 ,j ≠ " -,n 0,k≠ 1 2 ,j ) v -.
So, the following equation on the total electric charge is obtained:

fl n+1 k = fl n k ≠ t x ÿ - q - ÿ j (" -,n 0,k+ 1 2 ,j ≠ " -,n 0,k≠ 1 2 ,j ) v -.
According to (37), we get:

(41) fl n+1 k = fl n k ≠ t 2 x ÿ - q - ÿ j (v - j (g -,n
0,k+1,j ≠g -,n 0,k≠1,j )≠|V max |(g -,n 0,k+1,j ≠2g -,n 0,k,j +g -,n 0,k≠1,j )) v -.

By rearranging the terms, we obtain:

(42) fl n+1 k = fl n k ≠ t 2 x (j n k+1 ≠ j n k≠1 ) + t 2 x ÿ - q - ÿ j |V max |(g -,n
0,k+1,j ≠ 2g -,n 0,k,j + g -,n 0,k≠1,j ) v - By assuming that fl n k = j n k+1 = j n k≠1 = 0, it comes:

(43) fl n+1 k = t 2 x
V max (fl n k+1 ≠ 2fl n k + fl n k≠1 ) = 0. Hence, quasineutrality is conserved at all time. However, applying this modification to our scheme amounts to overevaluating the numerical viscosity by its highest value. This will create an important amount of dissipation in the numerical solutions obtained with this scheme. Hence, the discretization of the space domain needs to be very precise to compensate this viscosity, which will result in a longer computation time.

Numerical results

In this section, numerical resolution of Riemann problems with our method is presented. Such problems are characterized by initial conditions constituted of two constant states separated by a discontinuity. In the study of such problems, an exact solution can be provided for the conserved quantities (density, momentum and total energy). However, it is not true for non-conserved quantities, and no reference results exist.

Comparative results obtained by an HLL-type scheme and a Suliciu relaxation method on the bi-temperature Euler system are displayed along the solutions from the kinetic scheme that has been derived in this paper (see [START_REF] Aregba-Driollet | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF] for more details on these methods).

The spatial domain is chosen as the interval [0,1], supplemented with homogeneous Neumann boundary conditions. Physical constants (k B , Z, q i ) are chosen equal to one. All test cases are done in the hydrodynamic limit • ae 0. Concerning the choice of l -in (30), a convergence test has been performed on the initial conditions of the test cases in order to obtain a appropriate value that would optimize computational cost. Hence, in all test cases, l -= 8, foroe {e, i}.

A mass ratio of 10 (m e = 0.1 and m i = 1) is considered. The greater the mass ratio, the greater the value of V max in (37). Hence, a high mass ratio implies an important amount of numerical viscosity. Then, we are restricting ourselves to a mass ratio of 10, which is already demanding in terms of number of discretization points. However, in a proper physical context, mass ratios of the order of 1000 should be considered. Shock tube with identical initial temperatures. The Riemann problem of the shock tube with identical ionic and electronic temperature is given by the following initial condition:

I n -(t 0 , x) = 1, u -(t 0 , x) = 0, T -(t 0 , x) = 1 if x oe [0, 0.5], n -(t 0 , x) = 0.125, u -(t 0 , x) = 0, T -(t 0 , x) = 1 if x oe [0.5, 1].
Parameters are chosen as N x = 120000, N v = 40 and l = 8. The inter-species collision relaxation time is • ei = 0.1. Results are computed at time t = 0.1. For fluid methods, 4000 discretization points are used.

The results are provided in figure 1 and2.

The quantities described by conservative equations (densities and velocities) are consistant with the fluid methods HLL and Suliciu. Moreover, quasi-neutrality is shown to be achieved, since electronic and ionic densities are superposed, as well as electronic and ionic velocities.

For the partial temperatures presented in figure 2, the intermediate state reached between the rarefaction wave and the contact discontinuity is identical for all schemes, which is expected since this part of the solution is smooth. However, the states obtained between the contact discontinuity and the shock are di erent for every scheme. Figures 3 display a zoom on this part of the solution. The results can be seen to be di erent for all three methods. The behaviour of the di erent schemes should be investigated further. Shock tube with di erent initial temperatures. This test case is similar to the previous one, the only di erence being the initial left ionic and electronic temperatures. The following initial

I n -(t 0 , x) = 1, u -(t 0 , x) = 0, T -(t 0 , x) = 1 if x oe [0, 0.5], n -(t 0 , x) = 0.125, u -(t 0 , x) = 0, T e (t 0 , x) = 2 T i (t 0 , x) = 3 if x oe [0.5, 1].
Parameters are chosen as N x = 120000, N v = 40 and l = 8. The inter-species collision relaxation time is • ei = 0.1. Results are computed at time t = 0.1. The result are provided in figure 4 and5. The behaviour of the result if very similar to the previous one. Quasi-neutrality is achieved and temperatures exhibit di erent jump relations across discontinuities for each method. Rarefaction wave. This test case is constituted of two rarefaction waves going in opposite directions. It is given by the following data:

I n -(t 0 , x) = 1, u -(t 0 , x) = ≠1, T -(t 0 , x) = 1 if x oe [0, 0.5], n -(t 0 , x) = 1, u -(t 0 , x) = 1, T -(t 0 , x) = 1 if x oe [0.5, 1].
We have N x = 120000, N v = 40 and l = 8. The inter-species collision relaxation time is • ei = 0.1. Results are computed at time t = 0.1. The result are provided in figure 7 and8. In this test case, the solution is smooth. There are no shock waves, and then, all methods are expected to give the same results. Figure 9 shows that this is what happens. Shock wave. This test case is constituted of two shock waves. It is given by the following data:

I n -(t 0 , x) = 1, u -(t 0 , x) = 1, T -(t 0 , x) = 1 if x oe [0, 0.5], n -(t 0 , x) = 1, u -(t 0 , x) = ≠1, T -(t 0 , x) = 1 if x oe [0.5, 1].
We have N x = 120000, N v = 40 and l = 8. The inter-species collision relaxation time is • ei = 0.1. Results are computed at time t = 0.1.

The results are provided in Figure 10 and 11. Quasi-neutrality is achieved on the conserved quantities. The solution contains two shock waves and di erent jump relations can be observed for each numerical method. In Figure 12, a zoom is performed on the constant values reached by each method between the two shocks. This shows the di erent behaviour between the three methods.

Figure 5. Electronic and ionic temperatures of a shock tube test case with a mass ratio of 10 with 120000 space points, 40 velocity points and a domain length of 8 Figure 6. Di erent jump relations of electronic and ionic temperatures for a shock tube test case with a mass ratio of 10 with 120000 space points, 40 velocity points and a domain length of 8

Conclusion

In this article, a kinetic numerical method able to provide reference results for a non-conservative hyperbolic system, the bi-temperature Euler system, is proposed. This method is able to enforce all the desired properties that are relevant for comparisons with methods applied directly to the hyperbolic system. Such methods lack an unambiguous definition of the non-conservative products and exhibit di erent Rankine-Hugoniot relation when the solution contains shock waves. However, the proposed method introduces a important amount of numerical viscosity, which renders the method costly. A second-order method able to circumvent such a drawback is a current on-going work, which should be the subject of a subsequent article.

Figure 1 .Figure 2 .

 12 Figure 1. Density and velocity solutions of shock tube test case with a mass ratio of 10 with 120000 space points, 40 velocity points and a domain length of 8

Figure 3 .

 3 Figure3. Di erent jump relations through shocks of electronic and ionic temperatures for a shock tube test case with a mass ratio of 10 with 120000 space points, 40 velocity points and a domain length of 8

Figure 4 .

 4 Figure 4. Density and velocity solutions of shock tube test case with a mass ratio of 10 with 120000 space points, 40 velocity points and a domain length of 8

Figure 7 . 8 Figure 8 . 8 Figure 9 .

 78889 Figure 7. Density and velocity solution of a rarefaction wave test case with a mass ratio of 10 with 120000 space points, 40 velocity points and a domain length of 8

Figure 10 .Figure 11 . 8 Figure 12 .

 1011812 Figure 10. Density and veloctiy solutions of a shock wave test case with a mass ratio of 10 with 120000 space points, 40 velocity points and a domain length of 8

By formally passing to limit ' -ae 0, it comes g - p = Mg - p , n e = n i and u = u e = u i . Hence, summing overand using [START_REF] Xu | Multiple temperature kinetic model and gas-kinetic method for hypersonic nonequilibrium flow computations[END_REF], the inequality becomes

The right hand side reads

)

> .

The first and third term are equal to zero according to the conservation properties. The second term is equal, after computation, to :

> .

The total right hand side is then equal to

According to the definition of T , it comes:

and we get (25). ⇤

Numerical treatment

In this section, numerical treatment of system ( 22) is addressed. The method to compute the evolution of the particle distribution is presented. For the sake of simplicity, in the Ampère equation and BGK operator, every non-dimensional quantity Á e , Á i andare chosen equal to • . Reaching the hydrodynamic limit is then equivalent to the limit • ae 0.