D Aregba-Driollet

J Breil

S Brull

B Dubroca

E Estibals

Mathematical Modelling and Numerical Analysis Will be set by the publisher Modélisation Mathématique et Analyse Numérique MODELLING AND NUMERICAL APPROXIMATION FOR THE NONCONSERVATIVE BITEMPERATURE EULER MODEL

Keywords: 1991 Mathematics Subject Classification. 65M08, 35L60, 35L65. Secondary: 82D10, 76X05 Relaxation method, nonconservative hyperbolic system, kinetic schemes, BGK models, hydrodynamic limit, entropy dissipation

This paper is devoted to the study of the nonconservative bitemperature Euler system.

We firstly introduce an underlying two species kinetic model coupled with the Poisson equation. The bitemperature Euler system is then established from this kinetic model according to an hydrodynamic limit. A dissipative entropy is proved to exist and a solution is defined to be admissible if it satisfies the related dissipation property. Next, four different numerical methods are presented. Firstly, the kinetic model gives rise to kinetic schemes for the fluid system. The second approach belongs to the family of the discrete BGK schemes introduced by Aregba-Driollet and Natalini. Finally, a quasi-linear relaxation approach and a Lagrange-remap scheme are considered.

Introduction

The present paper is devoted to the approximation of a nonconservative bitemperature compressible Euler system by different numerical schemes.

This fluid model consists of two conservation equations for mass and momentum and two nonconservative equations, that is to say, one for each energy. Physically, this model describes the interaction of a mixture of one species of ions and one species of electrons in thermal nonequilibrium. The pressure of each species is supposed to satisfy a gamma-law with its own γ constant.

Solving nonconservative hyperbolic systems is a delicate problem because the definition of weak admissible solutions remains unclear. In order to define nonconservative products, Dal Maso, Le Floch and Murat proposed in [START_REF] Maso | Definition and weak stability of nonconservative products[END_REF] a new theory based on the definition of family of paths. However, it is shown, in [START_REF] Abgrall | A comment on the computation of non-conservative products[END_REF], that even if the correct path is known, the numerical solution can be far from the expected solution. In [START_REF] Berthon | Why many theories of shock waves are necessary : Kinetic relations for nonconservative systems[END_REF], the authors consider an hyperbolic system having np equations in a conservative form, the remaining p equations being nonconservative. In [START_REF] Pares | Numerical methods for nonconservative hyperbolic systems: a theoretical framework[END_REF], path-conservative schemes are defined by using the concepts developped in [START_REF] Maso | Definition and weak stability of nonconservative products[END_REF]. In [START_REF] Pelanti | A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves[END_REF], the authors use a Roe solver, and an HLLC solver which neglects the nonconservative part of the system. They validate the approach by comparing their results to theoretical temperatures/pressures curves ( [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF], [START_REF] Saurel | Modelling phase transition in metastable liquids: application to cavitating and flshing flows[END_REF]).

In [START_REF] Coquel | Numerical methods for weakly ionized gas[END_REF], the authors assume that the electronic entropy is conserved by all weak solutions including shocks and the system is approached by a system of conservation laws.

The bitemperature Euler system can be constructed from an underlying kinetic model which consists of a BGK model coupled with Poisson equation in the quasi-neutral regime. This BGK model possesses different interspecies collision frequencies in order to take into account discrepancies in the particle masses. This point has been in particular mentionned in [START_REF] Greene | Improved bhatnagar-gross-krook model of electron-ion collisions[END_REF], where a conservative formulation is proposed. Hence, following this idea we obtain such a model satisfying a H theorem. Next, performing an hydrodynamic limit, we get the nonconservative Euler system. To take into account the different γ constants, we should consider polyatomic kinetic models. We do not make so in order to avoid heavy technical proofs. However, the ideas are analogous to the monoatomic case. We show that the fluid model owns a dissipative strictly convex entropy which can be recovered from the Boltzmann entropy. In this paper, we consider the entropic solutions of the nonconservative Euler system which are the hydrodynamic limits of solutions of the kinetic system. Using the Boltzmann entropy, we prove that the hydrodynamic limits are entropy-dissipative. Moreover this approach leads at the discrete point of view to numerical schemes which are consistent with the physics.

In another way, the present bitemperature Euler system is approached by using numerical methods inspired from the approximation of the Euler system. The first one is a kinetic scheme ( [START_REF] Perthame | Boltzmann type schemes for gas dynamics and the entropy property[END_REF]) which is based on the underlying kinetic model. The second method is based on discrete BGK schemes introduced in [START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF]. One important point of this paper is their generalization to the present nonconservative setting and the fact that those models provide entropy-dissipative solutions. We also prove discrete entropy inequalities for the related numerical schemes. The third method is based on a Suliciu relaxation approach for the pressure variables ( [START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF]). The interest of such an approach is that all the characteristic fields of the related homogeneous system are linearly degenerate and hence the Riemann solver is easier to handle. Finally, the last method of the paper is a Lagrangian method ( [START_REF] Maire | A cell-centered lagrangian scheme for two dimensional compressible flow problems[END_REF]) based on Lagrangian and projection steps. This method is used in the multiphysics CHIC code ( [START_REF] Breil | Multi-material ALE computation in inertial confinement fusion code CHIC[END_REF]) to solve a bitemperature model in the field of inertial confinement fusion (ICF) and high energy density physics (HEDP). Lagrangian formalism is well suited to deal with the multi-material flow encountered in the field of ICF and numerous papers devoted to laser plasma interaction use this code ( [START_REF] Olazabal-Loumé | Linear and non-linear high mode perturbations amplification at ablation front in HIPER targets[END_REF]).

This paper is organised as follows. The second section is dedicated to the physical models that are involved in this paper. Firstly Euler bitemperature macroscopic model is given and we consider the Vlasov-BGK model from which the Euler system is derived. Starting from an ad-hoc scaling the construction of the Euler system is performed. Next, we study the entropy properties of this system. In section 3, we design the different numerical schemes and in section 4, some numerical results are presented. The goal of those first tests is to compare the schemes on 1D examples. It is out of the scope of this paper to propose higher order improvments or multi-dimensionnal computations. Finally, section 5 deals with some conclusions and perspectives to this paper.

The physical models

This part is devoted to the presentation of the different physical models that are involved in this paper. Firstly, we give the Euler bitemperature system, the underlying kinetic model and next we derive the fluid model from an hydrodynamic limit. For the sake of simplicity we consider a new BGK model to represent the interaction between ions and electrons coupled with Poisson equation considered at the quasi-neutral regime. Then we prove the existence of a strictly convex dissipative entropy for the fluid model.

There are many BGK models devoted to gas mixtures ( [START_REF] Andries | A consistent BGK-type model for gas mixtures[END_REF], [START_REF] Brull | Derivation of BGK models for mixtures[END_REF], [START_REF] Brull | An ellipsoidal statistical model for gas mixtures[END_REF], [START_REF] Garzo | A kinetic model for a multicomponent gas[END_REF], [START_REF] Hamel | Kinetic model for binary gas mixtures[END_REF], [START_REF] Morse | Kinetic model equation for a gas mixture[END_REF]). However except ( [START_REF] Andries | A consistent BGK-type model for gas mixtures[END_REF], [START_REF] Brull | Derivation of BGK models for mixtures[END_REF], [START_REF] Brull | An ellipsoidal statistical model for gas mixtures[END_REF]), there are few BGK models for gas mixtures enjoying fundamental properties (H theorem, equilibrium states, nonegativity of the distribution function, . . . ) inherited from the Boltzmann operator ( [START_REF] Aoki | Knudsen layer for gas mixtures[END_REF]). In order to take into account disparate masses, we consider in this paper the situation where the interspecies collision frequencies are different. This point has been in particular mentionned in ( [START_REF] Greene | Improved bhatnagar-gross-krook model of electron-ion collisions[END_REF]) where a conservative formulation is proposed. Here we are able to prove an H-theorem and to make the link with the macroscopic entropy. Moreover, one originality of this work is to take into account the electric field in the rescaled kinetic model.

The Euler bitemperature model

The nonconservative two species Euler equations are given by the hyperbolic system:

       ∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρu) + ∂ x (ρu 2 + p e + p i ) = 0, ∂ 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 ), ∂ 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 ), (1)
where ρ = ρ e + ρ i ≥ 0 is the total density of the plasma, u is the average velocity of the plasma. T e and T i represent the temperatures of electrons and ions. ρ e = n e m e , ρ i = n i m i are the density of the electrons and ions, where the concentrations of electrons n e and of ions n i are related by the average ionization number Z = n e /n i ≥ 1. Z will be consider here constant and as a physical property of the ions. m e and m i are the mass of the electrons and ions particles. Hence the mass fractions

c α = ρ α ρ , α = e, i (2)
are also constant and c e and c i write

c e = Zm e m i + Zm e , c i = 1 -c e . (3)
The electronic and ionic pressures and temperatures are related by

p e = n e k B T e , p i = n i k B T i ,
where k B is the Boltzmann constant. The internal energies are given by

ε e = k B T e m e (γ e -1)
,

ε i = k B T i m i (γ i -1) , (4)
where γ e , γ i are constant numbers belonging to the interval [START_REF] Abgrall | A comment on the computation of non-conservative products[END_REF][START_REF] Aoki | Knudsen layer for gas mixtures[END_REF] and ν ei ≥ 0 is the frequency exchange between temperatures.

Along this article we denote E α = ρ α ε α + 1 2 ρ α u 2 for α = e, i, and

U = (ρ, ρu, E e , E i ). ( 5
)
We also use in this article the expression of the bitemperature model in Lagrangian formalism, which writes as

       ρd t τ -∂ x u = 0, ρd t u + ∂ x (p e + p i ) = 0, ρ e d t ε e + p e ρd t τ = ν ei (T i -T e ), ρ i d t ε i + p i ρd t τ = -ν ei (T i -T e ), (6)
where d t = ∂ t + u∂ x represents the material derivative and τ = 1 ρ is the specific volume.

The kinetic model

Notations

Kinetic models are described by the distribution function f α of each species depending on the time variable t ∈ R + , on the position x ∈ R 3 and on the velocity v ∈ 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 α = R 3 f α dv, u α = 1 n α R 3 vf α dv, E α = 3 2 ρ α k B m α T α + 1 2 ρ α u 2 α = R 3 m α v 2 2 f α dv, (7)
where m α is the mass particle, ρ α = m α n α , and T α the temperature of species α.

The internal specific energy of species α can be defined as

ε α = 3 2m α k B T α .
Those definitions are consistent with (4) for γ α = 5 3 . As already pointed out in introduction, this value corresponds to the monoatomic case. This is not a conceptual restriction because the general case could be obtained with the same approach by considering polyatomic models. This generalization is postponed to a forthcoming paper. In the following we shall use the moment operator P α defined by

P α (f α ) = m α R 3   1 v v 2 2   f α dv. ( 8
)
We denote P α (f α ) = U α :

U α =   ρ α ρ α u α E α   . (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 = α ( 1 3 ρ α (u 2 α -u 2 )) + α (n α k B T α ), (10)
where n = n e + n i . Finally we define the entropy of the mixture by

H(f e , f i ) = H s (f e ) + H s (f i ), with H s (f ) = R 3 (f ln(f ) -f )dv. ( 11
)
The related entropy flux is then

Φ(f e , f i ) = Φ s (f e ) + Φ s (f i ), with Φ s (f ) = R 3 v(f ln(f ) -f )dv. ( 12
)

Description of the BGK model

In this section, we present the kinetic model and we show the fundamental properties of the BGK model describing the plasma interacting with an electric field E ∈ R 3 . This model writes

∂ t f α + v∇ x f α + q α m α E.∇ v f α = 1 τ α (M α (f α ) -f α ) + 1 τ αβ (M α (f e , f i ) -f α ), {α, β} ∈ {e, i}, α = β, (13)
with τ α > 0, τ αβ > 0. q α is the charge of the species α.

The relaxation term 1 τα represents the collision frequency for the interaction between α particles and 1 τ αβ corresponds to the collision frequency for the interaction between ions and electrons.

In this paper, contrarily to classical BGK models, we consider as in ( [START_REF] Greene | Improved bhatnagar-gross-krook model of electron-ion collisions[END_REF]) that τ ei = τ ie , due to the discrepancy of ion and electron masses. In particular, τie τei is of the same order as ci ce and can reach important values. Moreover this quantity is kept fixed in the physical model presented in this paper. The model ( 13) is coupled to the Maxwell-Ampère and the Poisson equations through the electric field E as follows

∂ t E = - j ε 0 , (14)
∇ x • E = ρ ε 0 . ( 15
)
j represents the current in the plasma, ρ the total charge, and ε 0 is the vacuum permitivity. j and ρ are defined by

ρ = R 3
(q e f e + q i f i ) dv = n e q e + n i q i , j = R 3

v(q e f e + q i f i ) dv = n e q e u e + n i q i u i .

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, (17)
M α (f e , f i ) = n α (2πk B T # /m α ) 3/2 exp(- |v -u # | 2 2k B T # /m α ), α = e, i, (18)
where

u # = 1 τei m e n e u e + 1 τie m i n i u i 1 τei m e n e + 1 τie m i n i , (19)
T # = 3 2 1 τei k B n e T e + 3 2 1 τie k B n i T i + 1 τei ( 1 2 ρ e (u e ) 2 -1 2 ρ e (u # ) 2 ) + 1 τie ( 1 2 ρ i (u i ) 2 -1 2 ρ i (u # ) 2 ) 3 2 k B 1 τei n e + 3 2 k B 1 τie n i . (20)
In [START_REF] Degond | The asymptotics of collision operators for two species of particles of disparate masses[END_REF][START_REF]Transport coefficients of plasmas and disparate mass binary gases[END_REF], u # and T # are fictitious quantities and are defined in such a way that the BGK model is conservative.

In the equilibrium case, the relations u e = u i = u and T e = T i = T correspond to u # = u and T # = T . Moreover, u # being defined in [START_REF] Degond | The asymptotics of collision operators for two species of particles of disparate masses[END_REF] according to a convex combinaison of u e and u i , it holds that

( 1 τ ei ρ e + 1 τ ie ρ i )(u # ) 2 ≤ 1 τ ei ρ e u 2 e + 1 τ ie ρ i u 2 i .
Therefore T # defined by relation [START_REF]Transport coefficients of plasmas and disparate mass binary gases[END_REF] is positive. In the situation where τ ei = τ ie , we recover u # = u and T # = T , where u and T are given by [START_REF] Breil | Multi-material ALE computation in inertial confinement fusion code CHIC[END_REF]. Remark that in [START_REF] Andries | A consistent BGK-type model for gas mixtures[END_REF], [START_REF] Brull | Derivation of BGK models for mixtures[END_REF], [START_REF] Brull | An ellipsoidal statistical model for gas mixtures[END_REF], the authors consider only one global BGK operator per species in order to reproduce the interaction between the given species with the other species.

The following properties are well known:

Property 2.1. For α = e, i, let F α be the flux function of Euler equations and η α (U α ) = η α (ρ α , ε α ), Q α = uη α the usual entropy-flux pair for the Euler system with the γ-law (4):

η α (ρ α , ε α ) = - ρ α m α (γ α -1) ln (γ α -1)ρ α ε α ρ γα α + C . ( 21
)
Here C is a nonnegative constant. If γ e = γ i = 5/3, then for α = e, i:

P α (M α (f α )) = U α , P α (vM α (f α )) = F α (U α ), (22)
and

H s (M α (f α )) = η α (U α ), Φ s (M α (f α )) = Q α (U α ). ( 23
)
Proposition 2.1. The model [START_REF] Brull | Derivation of BGK models for mixtures[END_REF][START_REF] Coquel | Numerical methods for weakly ionized gas[END_REF][START_REF] Maso | Definition and weak stability of nonconservative products[END_REF][START_REF] Degond | The asymptotics of collision operators for two species of particles of disparate masses[END_REF][START_REF]Transport coefficients of plasmas and disparate mass binary gases[END_REF] conserves the mass per species, the total momentum and the total energy.

Proof. The two Maxwellian distributions [START_REF] Coquel | Numerical methods for weakly ionized gas[END_REF][START_REF] Maso | Definition and weak stability of nonconservative products[END_REF] satisfy the constraints

R 3 (M α -f α )   m α m α v m α v 2 2   dv = 0, R 3 (M α -f α )dv = 0, α = e, i. (24)
Straightforward computations give

R 3 1 τ ei (M e -f e )m e vdv + R 3 1 τ ie (M i -f i )m i vdv = 1 τ ei m e n e (u # -u e ) + 1 τ ie m i n i (u # -u i ) R 3 1 τ ei (M e -f e ) 1 2 m e v 2 dv + R 3 1 τ ie (M i -f i ) 1 2 m i v 2 dv = 3 2 k B ( 1 τ ei n e (T # -T e ) + 1 τ ie n i (T # -T i )) + 1 2τ ei (ρ e (u # ) 2 -ρ e u 2 e ) + 1 2τ ie (ρ i (u # ) 2 -ρ i u 2 i ).
Then the definition of u # and T # given by ( 19) and ( 20) leads to

1 τ ei R 3 ((M e -f e ) m e v m e v 2 2 + 1 τ ie (M i -f i ) m i v m i v 2 2
)dv = 0.

Therefore the conservation properties follow.

Theorem 2.1. The model [START_REF] Brull | Derivation of BGK models for mixtures[END_REF][START_REF] Coquel | Numerical methods for weakly ionized gas[END_REF][START_REF] Maso | Definition and weak stability of nonconservative products[END_REF] satisfies the H theorem.

(1) The model satisfies the entropy inequality

1 τ e R 3 (M e -f e )ln(f e )dv + 1 τ i R 3 (M i -f i )ln(f i )dv + 1 τ ei R 3 (M e -f e )ln(f e )dv + 1 τ ie R 3 (M i -f i )ln(f i )dv ≤ 0. (26)
(2) The equality holds in the above equation if and only if there exists (n α , u, T

) ∈ R + × R 3 × R + such that f α = n α (2πk B T /m α ) 3/2 exp(- |v -u| 2 2k B T /m α ), α = e, i. (27)
Proof. For each α = e, i, by using the conservation properties [START_REF] Hamel | Kinetic model for binary gas mixtures[END_REF][START_REF] Liu | Hyperbolic conservation laws with relaxation[END_REF], it holds that

R 3 (M α -f α )ln(f α )dv = R 3 (M α -f α )ln( f α M α )dv (28)
and

R 3 ( 1 τ ei (M e -f e )ln(f e ) + 1 τ ie (M i -f i )ln(f i ))dv = R 3 ( 1 τ ei (M e -f e )ln( f e M e ) + 1 τ ie (M i -f i )ln( f i M i ))dv.
Therefore the convexity of the function x → x ln x yields the first part of the theorem. Moreover if the equality holds in [START_REF] Maire | A cell-centered lagrangian scheme for two dimensional compressible flow problems[END_REF], it comes that f α = M α , M α = M α , α = e, i. Next, by using the expressions [START_REF] Coquel | Numerical methods for weakly ionized gas[END_REF][START_REF] Maso | Definition and weak stability of nonconservative products[END_REF][START_REF] Degond | The asymptotics of collision operators for two species of particles of disparate masses[END_REF][START_REF]Transport coefficients of plasmas and disparate mass binary gases[END_REF] of M α and M α , it comes that u e = u i = u = u # , T i = T e = T = T # and we get [START_REF] Morse | Kinetic model equation for a gas mixture[END_REF]. Therefore we recover the same equilibrium states as for the Boltzmann operator for a two component gas ( [3]).

Hydrodynamic limit

In this section we derive as an hydrodynamic limit, a nonconservative Euler system. We firstly rescale the model and next perform the hydrodynamic limit.

Scaling on the one dimensionnal BGK model

For the sake of clarity, we assume that the system [START_REF] Brull | Derivation of BGK models for mixtures[END_REF][START_REF] Chalons | Computing material fronts with a Lagrange-projection approach[END_REF][START_REF] Chalons | An all-regime lagrange-projection like scheme for the gas dynamics equations on unstructured meshes[END_REF][START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF][START_REF] Coquel | Numerical methods for weakly ionized gas[END_REF][START_REF] Maso | Definition and weak stability of nonconservative products[END_REF][START_REF] Degond | The asymptotics of collision operators for two species of particles of disparate masses[END_REF][START_REF]Transport coefficients of plasmas and disparate mass binary gases[END_REF] is even in (v 2 , v 3 ). Such property being preserved by the kinetic equation [START_REF] Brull | Derivation of BGK models for mixtures[END_REF][START_REF] Chalons | Computing material fronts with a Lagrange-projection approach[END_REF][START_REF] Chalons | An all-regime lagrange-projection like scheme for the gas dynamics equations on unstructured meshes[END_REF][START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF][START_REF] Coquel | Numerical methods for weakly ionized gas[END_REF][START_REF] Maso | Definition and weak stability of nonconservative products[END_REF][START_REF] Degond | The asymptotics of collision operators for two species of particles of disparate masses[END_REF][START_REF]Transport coefficients of plasmas and disparate mass binary gases[END_REF], the distribution function f α of the species α depends on the time variable t ∈ R + , the space variable x ∈ R and the velocity variable v (v 1 , v 2 , v 3 ) ∈ R 3 . In the following, the macroscopic velocity u(t, x), the current j(t, x) belong to R and [START_REF] Berthon | A local entropy minimum principle for deriving entropy preserving schemes[END_REF][START_REF] Brull | Local discrete velocity grids for deterministic rarefied flow simulations[END_REF][START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF] are to be read with

u α = 1 n α R 3 v 1 f α dv, Φ s (f ) = R 3 v 1 (f ln(f ) -f )dv, j = R 3 v 1 (q e f e + q i f i ) dv.
The electric field belongs to R, and P α (f α ) = U α belongs to R 3 :

P α (f α ) = m α R 3   1 v 1 v 2 2   f α dv. ( 29
)
The property [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] reads then as

P α (M α (f α )) = U α , P α (v 1 M α (f α )) = F α (U α ). (30)
In order to obtain the quasi-neutral limit, the system (13-15) is rescaled in the following way

   ∂ t f α + v 1 ∂ x f α + qα mα E∂ v1 f α = 1 ε (M α -f α ) + 1 τ αβ (M α -f α ), α = β ∂ t E = -j ε 2 , ∂ x E = ρ ε 2 , ( 31
)
where ε is a positive parameter proportional to the Knudsen number. In particular the Maxwellian distributions given in [START_REF] Coquel | Numerical methods for weakly ionized gas[END_REF][START_REF] Maso | Definition and weak stability of nonconservative products[END_REF] write

M α (f α ) = n α (2πk B T α /m α ) 3/2 exp(- (v 1 -u α ) 2 + v 2 2 + v 2 3 2k B T α /m α ), α = e, i (32)
and

M α (f e , f i ) = n α (2πk B T # m α ) 3/2 exp(- (v 1 -u # ) 2 + v 2 2 + v 2 3 2k B T # /m α ), α = e, i, (33)
where u # and T # are defined by [START_REF] Degond | The asymptotics of collision operators for two species of particles of disparate masses[END_REF] and [START_REF]Transport coefficients of plasmas and disparate mass binary gases[END_REF].

Remark 2.1. In ( [START_REF] Degond | The asymptotics of collision operators for two species of particles of disparate masses[END_REF], [START_REF]Transport coefficients of plasmas and disparate mass binary gases[END_REF]), the authors perform formally hydrodynamic limits toward nonconservative Euler systems in the context of mixtures of ions and electrons. In particular they consider for their scaling the mass ratio between ions and electrons as a small parameter. In the present case, the scaling proposed in (31) is different. Our aim is to provide a kinetic approximation of the system [START_REF] Pelanti | A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves[END_REF] in the spirit of ( [START_REF] Perthame | A kinetic scheme for the saint-venant system with a source term[END_REF]). The comparison between the different scalings at the physical point of view has to be investigated and is postponed to a future paper.

Derivation of the Euler equations

The two species Euler system is obtained starting from [START_REF] Pelanti | A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves[END_REF][START_REF] Perthame | Boltzmann type schemes for gas dynamics and the entropy property[END_REF][START_REF] Perthame | A kinetic scheme for the saint-venant system with a source term[END_REF] and performing an hydrodamic limit.

Proposition 2.2. The system [START_REF] Pelanti | A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves[END_REF][START_REF] Perthame | Boltzmann type schemes for gas dynamics and the entropy property[END_REF][START_REF] Perthame | A kinetic scheme for the saint-venant system with a source term[END_REF] formally converges to the nonconservative two species Euler equations (1) where E is given according to the Ohm's law

1 ρ e ∂ x p e - 1 ρ i ∂ x p i = ( 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 (34)
and

ν ei = 3k B n e n i 2(τ ie n e + τ ei n i ) . ( 35
)
Remark that the three dimensional version of ( 1) can be obtained by considering the system (13-15) with the same scaling as in [START_REF] Pelanti | A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves[END_REF].

Proof. For this proof, for any g belonging to

L 1 2 = {g ∈ L 1 /(1 + v 2 )f ∈ L 1 }, we use the notation ∀g ∈ L 1 2 , g = R 3
g dv.

To obtain the Euler system (1), we perform a Chapman-Enskog expansion for each distribution function

f α = f 0 α + εf 1 α , α ∈ {e, i},
with the constraints

R 3 f α dv = R 3 f 0 α dv, R 3 f 1 α dv = 0, (36)
R 3 v 1 (m e f e + m i f i )dv = R 3 v 1 (m e f 0 e + m i f 0 i )dv, R 3 v 1 (m e f 1 e + m i f 1 i )dv = 0, (37)
R 3 v 2 m α f α dv = R 3 v 2 m α f 0 α dv, R 3 v 2 m α f 1 α dv = 0. ( 38
)
By identifying the terms of (31) of order ε -1 , we get

f 0 α = M α , α ∈ {e, i},
where M α is given by [START_REF] Perthame | Boltzmann type schemes for gas dynamics and the entropy property[END_REF]. Letting ε tend to 0, it comes that ρ = 0 and j = 0. Next, we use the relations

q e = -e, q i = Ze ( 39
)
where Z represents the ionisation rate and e the elementary charge. Then we deduce

n e = Zn i , u e = u i = u. (40)
Therefore,

f 0 α = M α , with u α = u.
(41) Moreover, since u e = u i , the relation [START_REF] Degond | The asymptotics of collision operators for two species of particles of disparate masses[END_REF] implies u # = u = u e = u i . The terms of order ε 0 of (31) are given by

∂ t f 0 α + v 1 ∂ x f 0 α + q α m α E∂ v1 f 0 α = 1 τ αβ (M α -f 0 α ) -f 1 α , α, β ∈ {e, i}, α = β. (42)
Multiply (42) by m α , integrate w.r.t v, use that ∂ v1 f 0 α = 0 and f 1 = 0 leads to

∂ t m α f 0 α + ∂ x m α v 1 f 0 α = 0.
Hence we obtain the mass conservation equation for species α

∂ t ρ α + ∂ x (ρ α u) = 0, α ∈ {e, i}.
By summing on α this equation, we obtain the conservation of total mass. In order to obtain the conservation of momentum, multiply (42) by m α v 1 and integrate w.r.t v. Hence we get

∂ t m α v 1 f 0 α + ∂ x m α v 2 1 f 0 α + q α E v 1 ∂ v1 f 0 α = 1 τ αβ m α v 1 (M α -f 0 α ) -m α v 1 f 1 α , α = β. ( 43
)
The scaling of the second equation of [START_REF] Pelanti | A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves[END_REF] gives

R 3 v 1 (q e f 1 e + q i f 1 i ) dv = 0. ( 44
)
Hence combining (44) with (37) leads to

m α v 1 f 1 α = 0. Moreover, as m α v 1 f 0 α = ρ α u, it comes that m α v 2 f 0 α = ρ α u 2 + p α , where p α is the pressure of species α = e, i. As v 1 ∂ v1 f 0 α = ∂ v1 (v 1 f 0 α ) -f 0 α , we get v 1 ∂ v1 f 0 α = -n α . Moreover, using that m α v 1 (M α -f 0 α ) = 0, we obtain ∂ t (ρ α u) + ∂ x (ρ α u 2 + p α ) -q α n α E = 0. (45)
By summing term by term (45) and using [START_REF] Liu | Hyperbolic conservation laws with relaxation[END_REF], we get the second equation of the system (1). Moreover (45) can be rewritten on the form

∂ t u + u∂ x u + 1 ρ α ∂ x p α - q α n α ρ α E = 0. (46)
By substracting equation ( 46) for electrons and ions, we get Ohm's law [START_REF] Saurel | Modelling phase transition in metastable liquids: application to cavitating and flshing flows[END_REF]. In order to obtain one equation on each internal energy for each species, we multiply equation ( 31) by m α v 2 /2 and integrate w.r.t v

∂ t m α v 2 2 f 0 α + ∂ x m α v 2 2 f 0 α v 1 + q α E v 2 2 ∂ v1 f 0 α = 1 τ αβ m α v 2 2 (M α -f 0 α ) -m α v 2 2 f 1 α , α, β ∈ {e, i}, α = β. (47) Firstly, we compute m α v 2 2 f 0 α = ρ α ε α + 1 2 ρ α u 2 , m α v 2 2 f 0 α v 1 = u(ρ α ε α + 1 2 ρ α u 2 + p α ).
For the last term of the left-hand side of (47), the relation

v 2 2 ∂ v1 f 0 α = ∂ v1 ( v 2 2 f 0 α ) -v 1 f 0 α leads to q α E v 2 2 ∂ v1 f 0 α = -n α q α Eu. Let α = β ∈ {e, i}. Then (34) reads 1 ρ α ∂ x p α - 1 ρ β ∂ x p β = ρ ρ α ρ β . ( 48
)
Next for the right-hand side of (47), a direct computation yields

m α v 2 2 (M α -f 0 α ) = 3 2 n α k B (T # -T α ).
Moreover, the relation ( 20) implies that

T # -T α = 1 τ βα n β 1 τei n e + 1 τie n i (T β -T α ), β = α, {α, β} ∈ {e, i}, (49)
because

u # = u = u e = u i . So m α v 2 2 (M α -f 0 α ) = 3 2 n e n i τ ie n e + τ ei n i k B (T β -T α ), β = α, α ∈ {e, i}. (50)
Hence, using that m α v 2

2 f 1 α = 0 and (50), (47) leads to the 2 last equations of (1). Remark 2.2. The relation (48) shows that the force term of the kinetic equations leads as in ( [START_REF] Perthame | A kinetic scheme for the saint-venant system with a source term[END_REF]) to the nonconservative terms of the bitemperature Euler system (1).

Entropy dissipation

This section is devoted to the entropy study for the system (56) and is one of the crucial part of the paper. We firstly prove the existence of an entropy-entropy flux pair directly on the fluid system. Next, the entropy dissipation property is obtained through the kinetic point of view. It is shown that the two obtained entropy dissipation properties are compatible.

Existence of a dissipative entropy for the Euler bitemperature model

In this part we aim to identify an entropy for the system (1). More precisely, we are looking for a dissipative entropy-entropy flux pair (η, Q), defined as follows: Definition 2.1. Consider a hyperbolic quasilinear system

∂ t U + A(U)∂ x U = G(U). (51)
Let η be a strictly convex real valued function defined on an open domain Ω of R n . Let Q be a real valued function defined on Ω. (η, Q) is a dissipative entropy-entropy flux pair for the system (51) if for all U ∈ Ω:

η ′ (U)A(U) = Q ′ (U) (52) and η ′ (U)G(U) ≤ 0. ( 53
)
Let us perform the change of variable U = φ(V). We denote

Ã(V) = (φ ′ (V)) -1 A(φ(V))φ ′ (V), η(V) = η(φ(V)), Q(V) = Q(φ(V)).
It is an easy matter to prove that U is a solution of system (51) if and only if V is a solution of

∂ t V + Ã(V)∂ x V = (φ ′ (V)) -1 G(φ(V)), (54)
and that condition (52) is satisfied if and only if

η′ (V) Ã(V) = Q′ (V). ( 55
)
In our case, the primal variable U = (ρ, ρu, E e , E i ) can be changed into V = (ρ, u, ε e , ε i ). A straightforward calculation leads to the following system for V:

           ∂ t ρ + u∂ x ρ + ρ∂ x u = 0, ∂ t u + u∂ x u + ρ -1 ∂ x (p e + p i ) = 0, ∂ t ε e + u∂ x ε e + ρ -1 e p e ∂ x u = ρ -1 e ν ei (T i -T e ), ∂ t ε i + u∂ x ε i + ρ -1 i p i ∂ x u = ρ -1 i ν ei (T e -T i ). ( 56
)
Lemma 2.2. For α = e, i, let η α (ρ α , ρ α u, E α ) = η α (ρ α , ε α ) be as defined in [START_REF] Garzo | A kinetic model for a multicomponent gas[END_REF], and Q α = uη α . We define (η, Q) for system (56) as η(ρ, u, ε e , ε i ) = µ e η e (ρc e , ε e )

+ µ i η i (ρc i , ε i ), Q(ρ, u, ε e , ε i ) = µ e Q e (ρc e , u, ε e ) + µ i Q i (ρc i , u, ε i ) (57)
where µ e and µ i are constant. Hence relation (55) is satisfied.

We skip the proof of this lemma: it is a direct calculation where one uses the fact that the entropy η α does not depend on the variable u.

As a consequence, if we define η(U) = η(ρ, u, ε e , ε i ) and the fonction U α as

U α (U) = (c α ρ, c α ρu, E α ), α = e, i, (58)
equality ( 52) is satisfied with

η(U) = µ e η e (U e (U)) + µ i η i (U i (U)), Q(U) = uη(U). (59)
We can now prove the existence of a strictly convex entropy for system (1).

Theorem 2.3. Let (η, Q) be defined by (59), with η e , η i defined in lemma 2.2 and U e , U i defined in (58). Let us set µ e = µ i = 1.

(60) Then η is a strictly convex dissipative entropy for system (1) and Q is the related entropy flux. More precisely, any smooth solution of the system satisfies the following equality:

∂ t η(U) + ∂ x Q(U) = - ν ei k B T i T e (T i -T e ) 2 . (61)
Proof. It remains to prove convexity and dissipation by the source-term. Consider U = W. Then U α (U) = U α (W) for at least one α. As η e and η i are strictly convex and µ e and µ i are positive:

η(W) > η(U) + α=e,i µ α η ′ α (U α (U))(U α (W) -U α (U)).
As the U α are linear,

U α (W) -U α (U) = U ′ α (U)(W -U), and therefore η(W) > η(U) + η ′ (U)(W -U),
which proves that η is strictly convex.

Then by using (4) one has

µ e ∂(η e • U e ) ∂E e (U) -µ i ∂(η i • U i ) ∂E i (U) (T i -T e ) = - 1 k B (T i -T e ) µ e T e - µ i T i
and (61) follows.

Entropy inequality for hydrodynamic limits

Theorem 2.3 is valid for any value of γ e , γ i . In the present case, the kinetic model [START_REF] Pelanti | A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves[END_REF] being monoatomic only the values of γ e and γ i equal to 5/3 can be considered. However, recall that the problem can be overcome by introducing a polyatomic kinetic model. The existence of a strictly convex dissipative entropy allows us to precise which solutions are admissible. In that goal, we study the entropy dissipation induced by the kinetic model.

Proposition 2.3. Let U be a solution of [START_REF] Abgrall | A comment on the computation of non-conservative products[END_REF] with γ e = γ i = 5/3. If U is a limit of a solution to the kinetic model [START_REF] Pelanti | A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves[END_REF][START_REF] Perthame | Boltzmann type schemes for gas dynamics and the entropy property[END_REF][START_REF] Perthame | A kinetic scheme for the saint-venant system with a source term[END_REF] then the following inequality holds:

∂ t η(U) + ∂ x Q(U) ≤ - ν ei k B T i T e (T i -T e ) 2 . (62)
Proof. Multiply equation for f α of (31) by ln(f α ) and integrate w.r.t v gives

∂ t R 3 f α ln(f α ) -f α dv + ∂ x R 3 v 1 (f α ln(f α ) -f α )dv + q α m α R 3 E ∂ v1 (f α ln(f α ) -f α )dv = 1 ε R 3 (M α -f α ) ln(f α )dv + 1 τ ei R 3 (M α -f α ) ln(f α ) dv.
The third term of the left-hand side is zero. By using [START_REF] Natalini | A discrete kinetic approximation of entropy solutions to multidimensional scalar conservation laws[END_REF] and notations [START_REF] Brull | An ellipsoidal statistical model for gas mixtures[END_REF][START_REF] Brull | Local discrete velocity grids for deterministic rarefied flow simulations[END_REF]:

∂ t H s (f α ) + ∂ x Φ s (f α ) ≤ 1 τ αβ R 3 (M α -f α ) ln(f α ) dv, α ∈ {e; i}, β = α. ( 63
)
Now we sum over α = e, i and formally pass to the limit ε = 0. We obtain f α = M α , with ρ = 0 and j = 0. So P α (f α ) = U α (U). The left-hand-side tends to the left-hand-side of inequality (62) as a direct consequence of ( 23). Moreover we get for α ∈ {e, i}:

R 3 (M α -M α ) ln(M α ) dv = R 3 (M α -M α ) ln( n α (2πk B Tα mα ) 3 2
) dv -

R 3 (M α -M α ) (v -u) 2 2k B Tα mα dv.
The first term of the right-hand-side is equal to zero by [START_REF] Hamel | Kinetic model for binary gas mixtures[END_REF]. For the second, a straightforward computation gives:

R 3 (M α -M α ) (v -u) 2 2k B Tα mα dv = 3 2 n α T α (T # -T α ).
Therefore the right-hand side of (63) reads

A = - 3 2 1 τ ei n e T e (T # -T e ) + 1 τ ie n i T i (T # -T i ) .
We conclude by using (49) and ( 35).

This result is our motivation to define the admissibility of a solution of the bitemperature Euler model even in the general case γ e = γ i . Definition 2.2. A solution U of system (1) is said to be admissible if it satisfies inequality (62).

Numerical Approximation

The spatial discretization is defined by a step ∆x and discretization cells C j =]x j-1 2 , x j+ 1 2 [. We consider that ∆x is constant, except in the Lagrangian formalism below. The time step can be variable in the applications and is denoted ∆t: t 0 = 0, t n+1 = t n + ∆t.

We adopt the finite volume viewpoint: for an unknown V (x, t), we look for approximations V n j of the average of V at time t n on cells C j .

A kinetic scheme

In this section we design a numerical approximation of the bitemperature Euler model (1) by using the kinetic approximation [START_REF] Pelanti | A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves[END_REF][START_REF] Perthame | Boltzmann type schemes for gas dynamics and the entropy property[END_REF][START_REF] Perthame | A kinetic scheme for the saint-venant system with a source term[END_REF]. We use the well-known transport-projection method that can be described as follows.

An initial data U 0 being given, for α = e, i, we set ρ 0 α = ρ 0 c α , n 0 α = ρ 0 α /m α , and u 0 α = u 0 . Suppose that at time t n an approximate solution U n = (ρ n , (ρu) n , E n e , E n i ) is known and that we have been able to define U n e and U n i such that

ρ n α = ρ n c α , u n α = u n , α = e, i. (64)
As a consequence, n n α = ρ n α /m α is also well defined for α = e, i and we are able to compute the microscopic quantities f n e (v),

f n i (v) by f n α (v) = M α (U n α ), α = e, i.
(65) We recall that the moment operator P α is defined by [START_REF] Olazabal-Loumé | Linear and non-linear high mode perturbations amplification at ablation front in HIPER targets[END_REF] and satisfies [START_REF] Pares | Numerical methods for nonconservative hyperbolic systems: a theoretical framework[END_REF]. We obtain:

P α (f n α ) = U n α , P α (v 1 M α (U n α )) = F α (U n α ), ( 66
)
where F α is the flux function of Euler equations with γ = 5/3. Now, at the microscopic level, we approximate the equations

∂ t f α + v 1 ∂ x f α + q α m α E∂ v1 f α = 1 τ αβ (M α -f α ), α = e, i, β = α. ( 67
) Definition 3.1. For each value of v we define a numerical flux h α,j+ 1 2 (v) = h α (f α,j (v), f α,j+1 (v), v)
such that for all v, h α (., ., v) is Lipschitz continuous and for all f :

h α (f, f, v) = v 1 f. ( 68
)
For α ∈ {e, i} with α = β, we define f

n+ 1 2 α,j and U n+1 α as f n+ 1 2 α,j = f n α,j - ∆t ∆x (h n α,j+ 1 2 (v) -h n α,j-1 2 (v)) -∆t q α m α E n+1 j ∂ v1 f n+ 1 2 α,j + ∆t τ αβ (M α (f n+ 1 2 e,j , f n+ 1 2 i,j ) -f n+ 1 2 α,j ), (69)
and U n+1 α,j = P α (f

n+ 1 2 α,j ). ( 70
) E n+1 j
not being defined, we impose the two following relations, which are the relaxed limit of the last two equations of [START_REF] Pelanti | A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves[END_REF]:

     q e m e ρ n+1 e,j + q i m i ρ n+1 i,j = 0, q e m e
ρ n+1 e,j u n+1 e,j +

q i m i ρ n+1 i,j u n+1 i,j = 0. ( 71
)
Therefore u n+1 i,j = u n+1 e,j and we define U n+1 j by setting

ρ n+1 j = ρ n+1 e,j + ρ n+1 i,j , u n+1 j = u n+1 i,j = u n+1 e,j .
Consequently, relations (64) are satisfied at level n + 1. For α = e, i, we set

F α,j+ 1 2 = F α (U α,j , U α,j+1 ) (72)
with

F α (U α , V α ) = P α (h α (M α (U α ), M α (V α ), •)). (73)
This is a consistent numerical flux for Euler system of species α:

F α (U α , U α ) = P α (v 1 M α (U α )) = F α (U α ).
Now we proceed as in section 2.3.2. First we note again that

P α (∂ v1 f α ) = (0, -ρ α , -ρ α u).
We thus obtain the following relations:

c e ρ n+1 j u n+1 j = c e ρ n j u n j - ∆t ∆x F n e,j+ 1 2 ,2 -F n e,j-1 2 ,2 +
∆t q e m e E n+1 j ρ n+1 e,j ,

c i ρ n+1 j u n+1 j = c i ρ n j u n j - ∆t ∆x F n i,j+ 1 2 ,2 -F n i,j-1 2 ,2 + ∆t q i m i E n+1 j ρ n+1 i,j .
Hence, denoting

δ n j+ 1 2 = -c i F n e,j+ 1 2 ,2 + c e F n i,j+ 1 2 ,2 , (74)
we have

c i c e ρ n+1 j E n+1 j - q e m e + q i m i = 1 ∆x δ n j+ 1 2 -δ n j-1 2 .
Therefore by ( 3) and (39):

-

q e m e E n+1 j ρ n+1 e,j = q i m i E n+1 j ρ n+1 i,j = 1 ∆x δ n j+ 1 2 -δ n j-1 2 .
It is clear that δ n j+ 1 2 is consistent with -c i p e + c e p i :

δ n j+ 1 2 = δ(U n j , U n j+1 ), δ(U, U) = -c i p e + c e p i .
Finally, if we define F j+ 1 2 as

F j+ 1 2 =      F e,j+ 1 2 ,1 + F i,j+ 1 2 ,1 F e,j+ 1 2 ,2 + F i,j+ 1 2 ,2 F e,j+ 1 2 ,3 F i,j+ 1 2 ,3      (75)
we obtain the following consistent scheme.

Proposition 3.1. A consistent scheme for Euler bitemperature system (1) is constructed as follows. For all n ≥ 0 if U n = {U n j } j∈Z is the approximate solution of system (1) at time t n , we set

U n α,j = (c α ρ n j , c α ρ n j u n j , E α ), α = e, i. (76)
A kinetic flux h α is chosen as in definition 3.1. We then define the numerical fluxes F α,j+ 1 2 , F j+ 1 2 and δ j+ 1 2 by (72), (73), (74), (75). The approximate solution at time t n+1 is defined by the implicit system:

                     ρ n+1 j = ρ n j - ∆t ∆x F n j+ 1 2 ,1 -F n j-1 2 ,1 , ρ n+1 j u n+1 j = ρ n j u n j - ∆t ∆x F n j+ 1 2 ,2 -F n j-1 2 ,2 , E n+1 e,j = E n e,j - ∆t ∆x F n e,j+ 1 2 ,3 -F n e,j-1 2 ,3 -u n+1 j ∆t ∆x δ n j+ 1 2 -δ n j-1 2 + ∆tν ei (T n+1 i,j
-T n+1 e,j ),

E n+1 i,j = E n i,j - ∆t ∆x F n i,j+ 1 2 ,3 -F n i,j-1 2 ,3 + u n+1 j ∆t ∆x δ n j+ 1 2 -δ n j-1 2 -∆tν ei (T n+1 i,j -T n+1 e,j ). ( 77
)
At the implementation level, the values of ρ and u are computed explicitly. It just remains to handle the source terms, and this happens to be a linear problem.

Remark 3.1. The relations (77) lead to a conservative numerical scheme for the conservative system satisfied by ρ, ρu and E = E e + E i . Remark 3.2. (77) gives rise to a numerical scheme as soon as the fluxes coming from a discretisation of the conservative Euler system. In particular, starting from (77), the numerical fluxes do not need to come from a kinetic process. For example, HLLC fluxes for Euler can be chosen. This point needs further investigations.

As an example, suppose that

h α (f, g, v) = v 1 2 (f + g) - |v max | 2 (g -f ). ( 78
) Then F α (U α , V α ) = 1 2 P α (v 1 M α (U α ) + v 1 M α (V α ) -|v max |(M α (V α ) -M α (U α ))),
that is

F α (U α , V α ) = 1 2 (F α (U α ) + F α (V α )) - |v max | 2 (V α -U α ) . (79)
Consequently:

F j+ 1 2 = 1 2           ρ j u j + ρ j+1 u j+1 -|v max |(ρ j+1 -ρ j ) ρ j u 2 j + ρ j+1 u 2 j+1 + p e,j + p e,j+1 + p i,j + p i,j+1 -|v max |(ρ j+1 u j+1 -ρ j u j )
u j+1 (E e,j+1 + p e,j+1 ) + u j (E e,j + p e,j ) -|v max |(E e,j+1 -E e,j )

u j+1 (E i,j+1 + p i,j+1 ) + u j (E i,j + p i,j ) -|v max |(E i,j+1 -E i,j )          
.

The nonconservative part is approximated with

δ j+ 1 2 = - c i 2
(p e,j + p e,j+1 ) + c e 2 (p i,j + p i,j+1 ).

Finally, denoting η e = 1, η i = -1 the scheme reads as:

ρ n+1 j = ρ n j - ∆t 2∆x ρ n j+1 u n j+1 -ρ n j-1 u n j-1 + |v max |∆t 2∆x ρ n j+1 -2ρ n j + ρ n j-1 , (80)
ρ n+1 j u n+1 j = ρ n j u n j - ∆t 2∆x ρ n j+1 (u n j+1 ) 2 + p n e,j+1 + p n i,j+1 -ρ n j-1 (u n j-1 ) 2 -p n e,j-1 -p n i,j-1 (81)
+ |v max |∆t 2∆x ρ n j+1 u n j+1 -2ρ n j u n j + ρ n j-1 u n j-1 , E n+1 α,j = E n α,j - ∆t 2∆x u n j+1 (E n α,j+1 + p n α,j+1 ) -u n j-1 (E n α,j-1 + p n α,j-1 ) (82)
+ |v max |∆t 2∆x E n α,j+1 -2E n α,j + E n α,j-1 + η α ∆t ν ei (T n+1 i,j -T n+1 e,j ) -η α u n+1 j ∆t 2∆x -c i (p n e,j+1 -p n e,j-1 ) + c e (p n i,j+1 -p n i,j-1 ) , α = e, i.
For the numerical tests, we used the more general HLL flux, which is obtained by the choice

h α (f, g, v) = v 1 λ 3 λ 3 -λ 1 f - λ 1 λ 3 -λ 1 g - λ 1 λ 3 λ 3 -λ 1 (g -f ). ( 83
)
Here λ 1 and λ 3 are constants to be fixed. The fact that this scheme satisfies discrete entropy inequalities is not clear at this point because the transport scheme with flux (78) or ( 83) is not monotone for all v, but it is actually true, as it will be proved in section 3.2. In the situation where for any v 1 , one has |v 1 | ≤ |v max |, a discrete entropy inequality can be proved. This condition is fullfilled, when the velocity grid that is used is compactly supported. This situation is classical in kinetic theory for gas dynamics, where the authors consider velocity grids equal to [u -4 √ RT , u + 4 √ RT ] (see [START_REF] Brull | Local discrete velocity grids for deterministic rarefied flow simulations[END_REF]). In particular, the Maxwellian distributions are truncated beyond this domain. There is another way to proceed which consists in replacing the Maxwellian by compactly supported distributions shearing the same macroscopic quantities as the associated Maxwellian.

Discrete BGK schemes

In this section we start from the discrete BGK approach of [START_REF] Natalini | A discrete kinetic approximation of entropy solutions to multidimensional scalar conservation laws[END_REF], [START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF] to construct another family of schemes for system [START_REF] Abgrall | A comment on the computation of non-conservative products[END_REF]. An important feature of those schemes is that they satisfy discrete entropy inequalities. As a particular case, we recover the scheme (80-82), which is therefore entropic.

The model developped in section 2 deals with a monoatomic setting. That is why the values γ i = γ e = 5 3 are imposed. Here the situation γ i = γ e is allowed.

Models

We take a discrete BGK model with L velocities ( [START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF]) for Euler equations with γ = γ α , α = e, i: denoting F α (U ) the Euler flux function, there exists L functions M α,l : R 3 → R 3 and λ 1,α , . . . , λ L,α ∈ R such that

∀U ∈ R 3 , L l=1 M α,l (U ) = U, L l=1 λ α,l M α,l (U ) = F α (U ). ( 84
)
We denote M α = (M α,l ) 1≤l≤L , Λ α = diag(λ α,1 I 3 , . . . , λ α,L I 3 ) where I 3 is the 3 × 3 identity matrix, P = (I 3 , . . . , I 3 ) ∈ M 3,3L (R), so that relations (84) may be written as follows:

∀U ∈ R 3 , P M α (U ) = U, P Λ α M α (U ) = F α (U ). (85)
Moreover we denote

N α = q α m α   0 0 0 1 0 0 0 1 0   , N α = diag(N α , . . . , N α ) ∈ M 3L (R).
For all f = (f 1 , . . . , f L ) ∈ R 3L , the macroscopic variables ρ, u, E and ε are defined by:

P f = U =   ρ ρu E = ρε + 1 2 ρu 2   .
For α, β ∈ {e, i}, α = β, we define

τ α,l = ν ei m α k B (γ α -1) 1 ρ α f α,3,l - ρ α u α 2ρ 2 α f α,2,l , B αβ,l (f e , f i ) = (0, 0, τ β,l -τ α,l ) , B αβ = (B αβ,l ) 1≤l≤L
and we set

                       ∂ t f ε e + Λ e ∂ x f ε e -E ε (x, t)N e f ε e = 1 ε (M e (U ε e ) -f ε e ) + B ei (f ε e , f ε i ), ∂ t f ε i + Λ i ∂ x f ε i -E ε (x, t)N i f ε i = 1 ε (M i (U ε i ) -f ε i ) + B ie (f ε e , f ε i ), ∂ t E ε = - 1 ε 2 q e m e ρ ε e u ε e + q i m i ρ ε i u ε i , ∂ x E ε = 1 ε 2 q e m e ρ ε e + q i m i ρ ε i , (86)
with U ε α defined by U ε α = P f ε α , α = e, i. If f ε α → f α when ε tends to 0 (α = e, i), then passing to the limit in the last two equations of (86) we obtain that u e = u i = u, q e m e ρ e + q i m i ρ i = 0 .

Moreover q e = -e and q i = Ze, so that m i ρ e = Zm e ρ i . As c e and c i satisfy (3), we obtain that ρ α = c α ρ, α = e, i.

We have also for α = e, i: M α (U α ) = f α . We apply P on the two first equations of (86) and we obtain the relaxed system:

                   ∂ t ρ α + ∂ x (ρ α u) = 0, α = e, i, ∂ t (ρ α u) + ∂ x (ρ α u 2 + p α ) - q α m α Eρ α = 0, α = e, i, ∂ t (ρ e ε e + 1 2
ρ e u 2 ) + ∂ x (u(ρ e ε e + 1 2 ρ e u 2 + p e )) -q e m e Eρ e u = ν ei (T i -T e ),

∂ t (ρ i ε i + 1 2 ρ i u 2 ) + ∂ x (u(ρ i ε i + 1 2 ρ i u 2 + p i )) - q i m i Eρ i u = -ν ei (T i -T e ). ( 88
)
We then proceed as in the proof of Proposition 2.2: we retrieve Ohm's law and U is a solution of the bitemperature Euler system.

Compatibility of microscopic entropies

In this paragraph we suppose that the Maxwellian functions are of the following form:

M α,l (U α ) = ξ α,l U α + ζ α,l F α (U α ), 1 ≤ l ≤ L, α = e, i. (89)
Here ξ α,l , ζ α,l are real constants. For instance, this is the case for the model (100-101) below. We also suppose that those functions are one-to-one and that for all U α under consideration:

σ(M ′ α,l (U α )) ⊂]0, +∞[, α = e, i. (90)
This condition is related to the well-known Liu's subcharacteristic condition, see [START_REF] Liu | Hyperbolic conservation laws with relaxation[END_REF], [START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF]. The Euler entropy being fixed as in property 2.1, we define

G α,l (U ) = ξ α,l η α (U ) + ζ α,l Q α (U ), H α,l (f l ) = G α,l (M -1 α,l (f α,l )), 1 ≤ l ≤ L, α = e, i.
Then F. Bouchut proved that the following properties hold for α = e, i ( [START_REF] Bouchut | Construction of BGK models with a family of kinetic entropy for a given system of conservation laws[END_REF]):

• for l = 1, . . . , L, H α,l is convex, (E0)

• L l=1 H α,l (M α,l (U )) = η α (U ), (E1).
• for all f , if U f = P f , one has

L l=1 H α,l (M α,l (U f )) ≤ L l=1 H α,l (f l ). (E2)
Those properties are similar to the ones of the Boltzmann entropy for the physical Maxwellian distributions of paragraph 2 and the same formalism could be used for both of them. Using moreover the fact that

η ′ (U α )N α U α = 0, α ∈ {e, i} (91)
we can prove straightforwardly that our relaxation model preserves the entropy properties of the solutions:

Proposition 3.2.
If U is a solution of system (1) obtained as a limit of the discrete BGK model (86) then U is an admissible solution.

The counterpart of this result is that the related numerical scheme is entropic, as it is stated in the following paragraph.

An entropic numerical scheme

We proceed as in section 3.1. Here h α is a numerical flux for the linear system of transport equations

∂ t f ε α + Λ α ∂ x f ε α = 0,
that is Definition 3.2. For α = e, i we define a numerical flux h α = (h α,l ) 1≤l≤L :

∀l ∈ {1, . . . , L}, h α,j+ 1 2 ,l = h α,l (f α,j,l , f α,j+1,l ) : h α is Lipschitz continuous and satisfies the following consistency property

(∀ f ∈ R 3L ) h α (f, f ) = Λ α f. ( 92
)
An approximate solution U n being known and U n e , U n i having been defined, we set

f n α,l = M α,l (U n α ). (93)
For α, β ∈ e, i, with β = α, we define (f

n+ 1 2 α,j ) α=e,i as f n+ 1 2 α,j = f n α,j - ∆t ∆x (h n α,j+ 1 2 -h n α,j-1 2 ) + ∆tE n+1 j N α f n+ 1 2 α,j + ∆tB αβ (f n+ 1 2 α,j , f n+ 1 2 β,j ). ( 94
)
We define U n+1 α , α = e, i, as

U n+1 α,j = P (f n+ 1 2 α,j ) (95)
and we impose (71). The remaining of the method follows the lines of paragraph 3.1, simply replacing the numerical flux (73) by

F α (U α , V α ) = P (h α (M α (U α ), M α (V α ))).
(96) It is easy to see that this is a consistent numerical flux for Euler system of the species α. As a final result, we obtain the following consistent scheme. Proposition 3.3. A consistent scheme for Euler bitemperature system (1) is constructed as follows. For all n ≥ 0 if U n = {U n j } j∈Z is the approximate solution of system (1) at time t n , U n e,j and U n i,j are defined by (76). A kinetic flux h α is chosen as in definition 3.2. We then define the numerical fluxes F α,j+ 1 2 , F j+ 1 2 and δ j+ 1 2 by ( 72), ( 96), ( 74), (75). The approximate solution at time t n+1 is defined by the implicit system (77).

In all the sequel, we choose h α as the upwind flux: denoting for λ ∈ R

λ -= max(0, -λ), λ + = max(0, λ), ( 97
) we define ∀f, g ∈ R 3L , h α,l (f l , g l ) = λ + α,l f l -λ - α,l g l . ( 98
)
We then find

F α (U α , V α ) = L l=1 -λ - α,l M α,l (V α ) + λ + α,l M α,l (U α ) . (99)
A first choice of parameters is as follows, for both α = e, α = i. As the parameters are the same, we omit to mention α in the notations. The characteristic velocities are:

λ 1 < λ 3 , λ 2 = 1 2 (λ 1 + λ 3 ) (100)
while the Maxwellian functions depend on a parameter κ ≥ 0:

M 1 (U ) = (λ 3 -κ)U -F (U ) λ 3 -λ 1 , M 2 (U ) = 2κU λ 3 -λ 1 , M 3 (U ) = -(λ 1 + κ)U + F (U ) λ 3 -λ 1 . (101)
The subcharacteristic condition ∀λ ∈ σ(F ′ (U )), λ 1 + κ < λ < λ 3κ ensures that if κ > 0 then condition (90) is satisfied for l = 1, 2, 3. If κ = 0, then M 2 0, the model reduces to a 2 × 2 one and the condition (90) is satisfied.

The related numerical flux reads as an upwind flux if λ 1 λ 3 ≥ 0. Otherwise we obtain, setting again α:

F α (U α , V α ) = - λ 1 λ 3 -λ 1 F α (V α ) + λ 3 λ 3 -λ 1 F α (U α ) (102) + λ 1 λ 3 λ 3 -λ 1 (V α -U α ) + κ λ 3 -λ 1 ((-λ 1 -(λ 1 + λ 3 ) -)V α + (-λ 3 + (λ 1 + λ 3 ) + )U α ).
If κ = 0, one retrieves the HLL flux given by (83). If κ = 0 and λ 1 = -|v max | = -λ 3 , one retrieves (79). The entropic behaviour of our schemes is proved by the same method as in [START_REF] Aregba-Driollet | Explicit diffusive kinetic schemes for nonlinear degenerate parabolic systems[END_REF] because (91) holds. We denote

f n+ 1 2 α,j,l = f n α,j,l - ∆t ∆x h α,l (f n α,j,l , f n α,j+1,l ) -h α,l (f n α,j-1,l , f n α,j,l ) , 1 ≤ l ≤ L. ( 103
)
Theorem 3.1. Suppose that the Maxwellian functions are of form (89) and that condition (90) is satisfied. Suppose moreover that for α = e, i and for all l ∈ {1, . . . , L} there exists a numerical entropy flux

G n α,j+ 1 2 ,l = G α,l (f n α,j,l , f n α,j+1,l ) such that H α,l ( f n+ 1 2 α,j ) -H α,l (f n α,j ) ∆t + G n α,j+ 1 2 ,l -G n α,j-1 2 ,l ∆x ≤ 0. ( 104
)
Then the following discrete entropy inequality holds for the numerical scheme defined in proposition 3.3:

η(U n+1 j ) -η(U n j ) ∆t + Q n j+ 1 2 -Q n j-1 2 ∆x ≤ - ν ei k B T n+1 i,j T n+1 e,j (T n+1 i,j -T n+1 e,j ) 2 , ( 105
)
where the numerical entropy flux

Q n j+ 1 2 = Q(U n j , U n j+1 ) is defined by Q(U, V) = α=e,i L l=1 G α,l (M α,l (U α ), M α,l (V α )) .
In the applications, as we choose the upwind flux for the transport equations, the numerical entropy fluxes G α,l exist under the CFL condition

max{|λ α,l | ∆t ∆x , 1 ≤ l ≤ L, α = e, i} ≤ 1. ( 106
)
They read as

G α,l (f, g) = λ + α,l H α,l (f ) -λ - α,l H α,l (g)
. This is different from the kinetic scheme of section 3.1 where the CFL condition cannot be satisfied for all v, so that the scheme with flux (78) is not entropic.

In the present case we have

G α,l (M α,l (U α ), M α,l (V α )) = ξ α,l λ + α,l η α (U α ) -λ - α,l η α (V α ) + ζ α,l λ + α,l Q α (U α ) -λ - α,l Q α (V α ) .

3.3.

Relaxation approach (Suliciu)

Relaxed system and relations of Suliciu

We consider a relaxation approach to solve the bitemperature Euler equations [START_REF] Abgrall | A comment on the computation of non-conservative products[END_REF]. As for the Euler system described in [START_REF]Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF] for one species, the partial pressures p e and p i are replaced by passive scalar variables π e and π i . These new variables are relaxed towards partial pressures p e and p i . We consider the system under the form (56). The Suliciu relaxation system can be written as :

                                     ∂ t ρ + u∂ x ρ + ρ∂ x u = 0, ∂ t u + u∂ x u + 1 ρ ∂ x (π e + π i ) = 0, ∂ t ε e + u∂ x ε e + π e ρ e ∂ x u = ν ei ρ e (T i -T e ), ∂ t ε i + u∂ x ε i + π i ρ i ∂ x u = - ν ei ρ i (T i -T e ), ∂ t π e + u∂ x π e + a 2 c e ρ ∂ x u = 1 τ (p e -π e ), ∂ t π i + u∂ x π i + a 2 c i ρ ∂ x u = 1 τ (p i -π i ). ( 107
)
When the positive relaxation parameter τ tends to zero one has

π α = p α + δτ + O(τ 2 ), α = e, i. (108)
From equation (108) and system (107), we can deduce that :

δ = c α ρ (ρ 2 a 2 α -a 2 )∂ x u + O(τ ), a α = γ α p α ρ α .
Hence the stability condition needed by the parameter a is :

a 2 ≥ ρ 2 max(a 2 e , a 2 i ) . ( 109
)
It is necessary to diminish the numerical diffusion induced by the increasing of parameter a to insure stability. Hence, a new equation about a is introduced:

∂ t a + u∂ x a = 0. ( 110
)
Finally, with respect to the primal variable, system (107, 110) reads as follows:

                               ∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρu) + ∂ x (ρu 2 + π e + π i ) = 0, ∂ t E e + ∂ x (u(E e + π e )) -u∂ x (c i π e -c e π i ) = ν ei (T i -T e ), ∂ t E i + ∂ x (u(E i + π i )) + u∂ x (c i π e -c e π i ) = ν ei (T e -T i ), ∂ t ( ρ e π e (ac e ) 2 ) + ∂ x (u( ρ e π e (ac e ) 2 + 1)) = ρ e (ac e ) 2 τ
(p eπ e ),

∂ t ( ρ i π i (ac i ) 2 ) + ∂ x (u( ρ i π i (ac i ) 2 + 1)) = ρ i (ac i ) 2 τ (p i -π i ), ∂ t (ρa) + ∂ x (ρau) = 0. (111)
The Suliciu relaxation scheme for system (1) consists in three steps. The approximate solution U n of (1) being known, we define a n such that condition (109) is satisfied and we project the pressures onto equilibrium:

π n α = p n α , α = e, i.
Then we use an exact Riemann solver for the homogeneous system related to (111):

                               ∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρu) + ∂ x (ρu 2 + π e + π i ) = 0, ∂ t E e + ∂ x (u(E e + π e )) -u∂ x (c i π e -c e π i ) = 0, ∂ t E i + ∂ x (u(E i + π i )) + u∂ x (c i π e -c e π i ) = 0, ∂ t ( ρ e π e (ac e ) 2 ) + ∂ x (u( ρ e π e (ac e ) 2 + 1)) = 0, ∂ t ( ρ i π i (ac i ) 2 ) + ∂ x (u( ρ i π i (ac i ) 2 + 1)) = 0, ∂ t (ρa) + ∂ x (ρau) = 0. (112)
Remark that the nonconservative product u∂ x (c i π ec e π i ) makes difficult the definition of a weak solution for this system. However, we show here below that the solution of the Riemann problem is well defined and this is the only property we need. Let U n+ 1 2 be the obtained solution. The third step takes the temperatures interaction into account implicitly: the approximate solution of system (1) at time t n+1 is defined by

ρ n+1 = ρ n+ 1 2 , u n+1 = u n+ 1 2 and    E n+1 e = E n+ 1 2 e + ∆t ν ei (T n+1 i -T n+1 e
),

E n+1 i = E n+ 1 2 i + ∆t ν ei (T n+1 e -T n+1 i ).
This system is linear and owns an explicit solution. 3.3.2. Riemann solver for the homogeneous system (112)

ρ L , u L , ε e,L , ε i,L π e,L , π i,L , a L ρ * L , u * , ε * e,L , ε * i,L π * e,L , π * i,L , a L ρ * R , u * , ε * e,R , ε * i,R π * e,R , π * i,R , a R π e,R , π i,R , a R ρ R , u R , ε e,R , ε i,R
By using the variable U = (ρ, u, ε e , ε i , π e , π i , a), one can easily compute the eigenvalues of system (112). They read as {ua/ρ, u, u + a/ρ} where u is an eigenvalue of order 5. All the fields are linearly degenerated.

Let U L , U R be a Riemann data for system (112). We look for a weak solution as a superposition of three contact discontinuities propagating with velocities u L -aL ρL , u * (to be determined), u R + aR ρR with

u L - a L ρ L ≤ u * ≤ u R + a R ρ R , (113)
Riemann invariants can be given by :

Field (u -a ρ ) : u -a ρ , a, π e + a 2 ce ρ , π i + a 2 ci ρ , ε e - π 2 e 2(cea) 2 , ε i - π 2 i 2(cia) 2 . Field (u) : u, π e + π i . Field (u + a ρ ) : u + a ρ , a, π e + a 2 ce ρ , π i + a 2 ci ρ , ε e - π 2 e 2(cea) 2 , ε i - π 2 i 2(cia) 2 .
As a consequence, c i π ec e π i is a Riemann invariant for both extreme eigenvalues. This means that this quantity remains constant through the related contact discontinuities, so that u∂ x (c i π ec e π i ) = 0 there. For the central discontinuity, u is constant so that u∂ x (c i π ec e π i ) = ∂ x (u(c i π ec e π i )). Hence, this product is also well defined in the usual weak sense. Hence the superposition of the three discontinuities is actually a weak solution of system (112). The Riemann invariants lead to the following system:

                                         u L - a L ρ L = u * - a L ρ * L , u R + a R ρ R = u * + a R ρ * R , π e,L + a 2 L c e ρ L = π * e,L + a 2 L c e ρ * L , π e,R + a 2 R c e ρ R = π * e,R + a 2 R c e ρ * R , π i,L + a 2 L c i ρ L = π * i,L + a 2 L c i ρ * L , π i,R + a 2 R c i ρ R = π * i,R + a 2 R c i ρ * R , ε e,L - π 2 e,L 2(c e a L ) 2 = ε * e,L - (π * e,L ) 2 2(c e a L ) 2 , ε e,R - π 2 e,R 2(c e a R ) 2 = ε * e,R - (π * e,R ) 2 2(c e a R ) 2 , ε i,L - π 2 i,L 2(c i a L ) 2 = ε * i,L - (π * i,L ) 2 2(c i a L ) 2 , ε i,R - π 2 i,R 2(c i a R ) 2 = ε * i,R - (π * i,R ) 2 2(c i a R ) 2 , π * e,L + π * i,L = π * e,R + π * i,R . Denoting π = π e + π i :        1 ρ * L = 1 ρ L + a R (u R -u L ) + π L -π R a L (a L + a R ) , 1 ρ * R = 1 ρ R + a L (u R -u L ) + π R -π L a R (a L + a R ) , u * = a R u R + a L u L + π L -π R a L + a R , (114)
                         π * e,L = π e,L + a L c e π R -π L -a R (u R -u L ) a L + a R , π * i,L = π i,L + a L c i π R -π L -a R (u R -u L ) a L + a R , π * e,R = π e,R + a R c e π L -π R -a L (u R -u L ) a L + a R , π * i,R = π i,R + a R c i π L -π R -a L (u R -u L ) a L + a R ,                          ε * e,L = ε e,L + 1 2(c e a L ) 2 ((π * e,L ) 2 -π 2 e,L ), ε * i,L = ε i,L + 1 2(c i a L ) 2 ((π * i,L ) 2 -π 2 i,L ), ε * e,R = ε e,R + 1 2(c e a R ) 2 ((π * e,R ) 2 -π 2 e,R ), ε * i,R = ε i,R + 1 2(c i a R ) 2 ((π * i,R ) 2 -π 2 i,R ).
We must now provide some sufficient conditions on a L and a R in order to satisfy (113) and the realisability of the intermediate states, that is positivity of

ρ * L , ρ * R , ε * e,L , ε * i,L , ε * e,
       a L (a L + a R ) ≥ ρ L (p e,R + p i,R -p e,L -p i,L ), a R (a L + a R ) ≥ -ρ R (p e,R + p i,R -p e,L -p i,L ), a L ≥ ρ L c L , a R ≥ ρ R c R ,
where the two first conditions are set to satisfy (113).

The so defined solver is entropic for the case of Euler one temperature conservative system (see [START_REF] Berthon | A local entropy minimum principle for deriving entropy preserving schemes[END_REF] paragraph 4.2). By analogy, one can expect that this property remains true for the bitemperature system, but this needs to be proved.

Numerical Scheme

In this paragraph we detail the second step of the scheme, namely the use of the Riemann solver to define the approximate solution U

n+ 1 2 = {U n+ 1 2 j } j∈Z . At each interface {x j+ 1 2 } × [t n , t n+1 ] we use the Riemann solver with U L = U n j , U R = U n j+1
. The speed of the different waves are

λ 1,j+ 1 2 = u n j - a n j ρ n j , λ 2,j+ 1 2 = u * j+ 1 2 , λ 3,j+ 1 2 = u n j+1 + a n j+1 ρ n j+1 .
We define U(x, t) on [x j-1 2 , x j+ 1 2 ] × [t n , t n+1 ] as the superposition of the Riemann solutions (see fig 2). In order to avoid overlapping Riemann problems on [x j-1 2 , x j+ 1 2 ], the following CFL condition is used :

λ max ∆t ∆x ≤ 1 2 , λ max = max j∈Z ( max 1≤k≤3 (λ + k,j-1 2 , λ - k,j+ 1 2
)).

We then define U ∆x ∆t

x j-1 2 x j x j+ 1 2 U n j-1 U L j-1 2 U R j-1 2 U n j U L j+ 1 2 U R j+ 1 2 U n j+1 Figure 2. Structure of Riemann solutions on [x j-1 2 , x j+ 1 2 ] × [t n , t n+1 ].
We set :

U - j = 1 ∆x xj x j-1 2 U(x, ∆t)dx, U + j = 1 ∆x x j+ 1 2 xj U (x, ∆t)dx.
Finally, using notation (97) we obtain

             U - j = 1 2 U n j - ∆t ∆x λ + 3,j-1 2 (U n j -U R j-1 2 ) + λ + 2,j-1 2 (U R j-1 2 -U L j-1 2 ) + λ + 1,j-1 2 (U L j-1 2 -U n j-1 ) , U + j = 1 2 U n j + ∆t ∆x λ - 3,j+ 1 2 (U n j+1 -U R j+ 1 2 ) + λ - 2,j+ 1 2 (U R j+ 1 2 -U L j+ 1 2 ) + λ - 1,j+ 1 2 (U L j+ 1 2 -U n j ) , U n+ 1 2 j = U - j + U + j .

Lagrangian approximation

The system ( 6) is solved by using a splitting. The gas dynamic is first solved using Lagrangian formalism. Then the computation of the new ion and electron internal energies is done in two steps. First, after the Lagrangian step we have a new total internal energy which is the sum of ion and electron internal energy. To compute the electron internal energy we make the assumption that the entropy deposition due to the shock is done on the ion internal energy [START_REF] Zel | Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena[END_REF]. Thus, electrons follow an isentropic evolution. Then the exchange of energy between ion and electron is computed. Finally, the Lagrangian solution is conservatively interpolated on the initial grid. This approach is also called Lagrange-projection scheme ( [START_REF] Breil | Multi-material ALE computation in inertial confinement fusion code CHIC[END_REF], [START_REF] Chalons | Computing material fronts with a Lagrange-projection approach[END_REF], [START_REF] Chalons | An all-regime lagrange-projection like scheme for the gas dynamics equations on unstructured meshes[END_REF]).

Lagrangian phase

In Lagrangian hydrodynamic methods, a computational cell moves with the flow velocity its mass being timeinvariant. Thus, Lagrangian methods can capture contact discontinuities sharply in multi-material fluid flows.

Here the Lagrangian scheme is based on a cell-centered discretization of Lagrangian hydrodynamics equations. This scheme is written in total energy form. Here, e is the total specific energy.

   ρd t τ -∂ x u = 0, ρd t u + ∂ x p = 0, ρd t e + ∂ x pu = 0, (115)
where p = p i + p e and e are respectively the total pressure and the total energy. From total energy we can get internal energy from ε = e -1 2 u 2 . As the mesh is moving we have to add the local kinetic equation : d t x = u.

Numerical scheme

In order to find an approximation of the flux of (115) we use the Godunov Riemann solver and we denote p * , u * the intermediate quantities. By using the Riemann invariants (u * , p * ) are given by

       u * = p l -p r + a r u r + a l u l a r + a l , p * = a r p l + a l p r + a l a r (u l -u r ) a r + a l . ( 116
)
The problem is solved on the interval [0; L] which is discretized in N elements C 1 , . . . , C N . Furthermore, we can define the j th cell center as x j = 1 2 (x j+ 1 2 + x j-1 2 ). We perform first the integration in space on ∆x j = [x j-1/2 , x j+1/2 ] :

                     ρ d dt x j+ 1 2 x j-1 2 τ (x, t) dx -u(x j+ 1 2 , t) -u(x j-1 2 , t) = 0, ρ d dt
x j+ 1 2

x j-1 2 u(x, t) dx + p(x j+ 1 2 , t)p(x j-1 2 , t) = 0, ρ d dt

x j+ 1 2

x j-1 2 e(x, t) dx + pu(x j+ 1 2 , t)pu(x j-1 2 , t) = 0.

(117)

Let be ϕ j = 1 ∆x j

x j+ 1 2

x j-1 2 ϕ(x, t) dx, now we integrate in time :

                   ρ∆x j tn+1 tn dτ j dt dt - tn+1 tn u(x j+ 1 2 , t) -u(x j-
       m j (τ n+1 j -τ n j ) -∆t(u * j+ 1 2 -u * j-1 2 ) = 0, m j (u n+1 j -u n j ) + ∆t(p * j+ 1 2 -p * j-1 2 ) = 0, m j (e n+1 j -e n j ) + ∆t(p * j+ 1 2 u * j+ 1 2 -p * j-1 2 u * j-1 2 ) = 0, (119)
with m j = ρ j ∆x j , (pu) * = p * u * and where (p * , u * ) are given by (116). Thanks to (119) we can compute a discrete form for internal energy using the difference between the total energy and the kinetic energy.

Next remark that the scheme (119) is entropic. More precisely, considering a macroscopic expression s of the entropy (see [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]), the second law of thermodynamics writes

T ds dt = dε dt + p dτ dt .
In that case, it is shown in ( [START_REF] Maire | A cell-centered lagrangian scheme for two dimensional compressible flow problems[END_REF]), that the numerical scheme defined by (119) satisfies

m j T j ds j dt = ρ j c j u * j+ 1 2 -u j 2 + ρ j+1 c j+1 u * j-1 2 -u j 2 ≥ 0, (120)
corresponding to a semi-discrete entropy inequality.

Bitemperature model

Based on total energy conservation, the Lagrangian scheme must deal with the exchange between the kinetic energy and the internal energy of ion and electron. The system we aim to solve for the internal energy of ion and electron is

ρ dε e dt + p e dτ dt = ν ei (T i -T e ), ρ dε i dt + p i dτ dt = ν ei (T e -T i ).
As the scheme is written in total energy formulation after the Lagrangian hydrodynamic step, we need to solve the exchange of energy in the system. The new ion and electron internal energy are computed in two steps. First, after the Lagrangian step we have a new internal energy which is the sum of ion and electron internal energy. To compute the electron internal energy we make the assumption that the entropy deposition is done on the ion internal energy [START_REF] Zel | Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena[END_REF]. Thus, electrons follow the isentropic evolution

m d ε e dt + p e d V dt = 0,
where m = ρ V is the mass, V the volume and d V dt is the change of volume during the Lagrangian step. We can deduce ion internal energy from ε i = εε e , where ε = e -1 2 u 2 is the mean total internal energy. Then, to solve the energy transfer between ion and electron we solve

ρ dε e dt = ν ei (T i -T e ), ρ dε i dt = ν ei (T e -T i ).

Remap step

As we deal here with a Lagrange-projection like scheme the remapping phase corresponds here to a conservative interpolation of the physical variables from the Lagrangian grid at time t n+1 onto the initial grid. Since we are using a Lagrangian scheme where in the placement of the variables is cell-centered, we are developing a cell-centered remapping phase. In what follows, all quantities related to the initial grid are denoted here with the tilde accent. Let ψ be a physical variable of the flow defined on a cell by its piecewise constant representation ψj . Being given ψ(x, t), on the Lagrangian grid we want to compute ψj = 1 Ṽj

xj+ 1 2 xj-1 2 ψ(x, t) dx. (121)
That is, knowing the mean value of ψ(x, t) = ψ over each cell of the Lagrangian grid, we want to compute its mean value over each cell of the initial grid. After the remap of all the conservative quantities we can reconstruct our solution.

Numerical results

In this part, some comparisons of the numerical schemes designed in Section 3 are performed. When γ e = γ i , the total density, total momentum and total energy E = E e + E i are solutions to the conservative 3 × 3 Euler system. Therefore exact solutions of the Riemann problem are known for those quantities.

Firstly, we consider a robustness test case dedicated to a double rarefaction wave whose parameters are

k B = 1, γ e = γ i = γ = 5 3 , Z = 1, m e = 1, m i = 1.
Next we present two other test cases in order to check the validity of the source term approximation and to study the behaviour of shocks. In order to compare our results to more realistic physical situations ( [START_REF] Zel | Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena[END_REF]), the physical parameters are set as follows:

k B = 1, γ e = γ i = γ = 5 3 , Z = 1, m e = 10 -3 , m i = 1.
This choice allows us to isolate the small m e /m i mass ratio feature. The characteristic velocities σ i of system (1) are u (double), u ± a, with a = γ e p e + γ i p i ρ .

In all the test cases that are presented, the time step ∆t n satisfies the CFL condition over all the cells:

max 1≤i≤4 |σ i | ∆t n ∆x ≤ 0.25.

Double rarefaction wave

The double rarefaction waves test case does not involve shocks. The goal of this test case is to check the robustness of the different schemes.

The left and the right states of the Rieman problem are the following

ρ L = 1, u L = -10, T e,L = 1000, T i,L = 1, ρ R = 1, u R = 10, T e,R = 1000, T i,R = 1.
The solutions are computed for x ∈ [0, 1] at time t = 0.05 for 200 cells. On the pictures, the term AN3V is for the discrete kinetic scheme of subsection 3.2 with the choice κ = (λ 3λ 1 )/5 in (101), while HLL is for the kinetic scheme of subsection 3.1 with the flux (83). In both cases, λ 1 and λ 3 are computed locally, at each cell interface. The term relaxation refers to the scheme presented in subsection 3.3.

The numerical results presented in figures 3 and 4 are in good agreement with the exact solutions for the conserved quantities. The results of the Lagrange remaping method show important descrepancies on the ionic temperature. This fact is well known and is due to the isentropic assumptions on the electrons. Then the entropy of the ions is overestimated in a problem without any shock. We remark also that AN3V and HLL give very close results. Indeed if κ = 0, HLL and AN3V coincide. A different value of this parameter does not seem to improve the results. For these reasons, in the following we chose to focus on HLL and relaxation schemes. Others schemes based on the approach of subsection 3.2 are to be tested in a forthcoming work.

Source-term approximation: an analytical solution

We take initial data such that ∀x ∈ R ρ(x, 0) = 1, u(x, 0) = 10, T i (x, 0) + T e (x, 0) = T = 2.

As Then the solution of the 3 × 3 Euler system is a stationary shock. The solutions are computed for x ∈ [0, 1] at time t = 0.05 with 1000 cells. The results are depicted figures 6 and 7.

ρ e ε e + ρ i ε i = p e + p i γ -1 = ρk B (γ -1)(m i + m e ) (T e + T i ), 0
For ν ei = 0, we observe a contact discontinuity propagating at velocity u R = 2.749750250, see figure 7 left. For the 3 × 3 Euler equations, u is a single eigenvalue and the left and right values of a contact discontinuity lie on an integral curve of the eigenvector (-ρ, 0, ε), in (ρ, u, ε) variables. Thus, a contact discontinuity must involve a jump of ρ. For the Euler bi-temperature system, the eigenvalue u is double and the eigenvectors are r 2 = (0, 0, -(γ i -1)c i , (γ e -1)c e ) and r 3 = (-ρ, 0, ε e , ε i ). The observed contact discontinuity lies on an integral curve of r 2 : only electronic and ionic temperatures jump. We also observe that the intermediate values of those temperatures differ with respect with the scheme. This is due to the fact that each scheme has a different viscosity, and therefore converges to a different interpretation of the nonconservative products in the equations.

For ν ei = 100, both schemes give similar results. For all x, density, velocity and total energy remain the same and are not depicted again. Electronic and ionic temperatures are represented on figure 7 right. For x < 0.5, as T i = T e , the value of ν ei does not influence the solution. For x > tu R , by finite propagation speed, one can compute the value of T e and T i as in subsection 4.2. For x ∈ [0, tu R ] the result is qualitatively in coherence with the physical behavior of the plasma, as predicted in Zeldovitch ( [START_REF] Zel | Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena[END_REF]). In particular, we observe a high ionic temperature at the shock, and then a decrease, while the electronic temperature increases.

Conclusion

In this paper we have investigated at the modelling and numerical point of views the bitemperature Euler system.

At the modelling point of view, we have proposed a multicomponent BGK kinetic system enjoying fundamental mathematical properties (positivity of the distribution function, conservation properties, H theorem) coupled with Ampère and Poisson equations. Next, be perfoming an hydrodynamic limit the bitemperature Euler model has been recovered. Moreover, we have shown for this system the existence of an entropy that is consistent with Boltzmann entropy.

At the numerical level, we have investigated and compared four numerical schemes for the bitemperature Euler system. We have firstly presented a kinetic scheme based on a discrete hydrodynamic limit. Next, the numerical method of Aregba-Natalini based on discrete BGK schemes has been generalized to a nonconservative setting. Doing this we show that the previous kinetic scheme is entropic and we recover a HLL scheme for the non-conservative system. A Suliciu relaxation approximation has also been designed. Finally all these schemes are compared to a Lagrangian scheme. Moreover, the kinetic and discrete BGK approaches are consistent with We shall address to forthcoming papers the following different points. Firstly, we plan to propose high order and multi-dimensionnal versions of the schemes developped in this paper. Moreover, the derivation of a polyatomic BGK model based on an additional continous internal energy variable has to be performed. Indeed, this model will be devoted to lead to general γ e and γ i for the fluid system. The case of Navier-Stokes asymptotics has also to be considered. In another direction, we aim to take into account the magnetic field in the equations with a coupling through the Maxwell system.

Figure 1 .

 1 Figure 1. Structure of the Riemann problem.

x j+ 1 2 x j-1 2 U

 22 (x, ∆t)dx.

Figure 3 .

 3 Figure 3. Density, velocity and pressure computed with AN3V, HLL and Lagrangian remaping schemes compared with exact solution for the double rarefaction waves problem.

Figure 4 .

 4 Figure 4. Electronic and ionic temperatures computed with AN3V, HLL, Lagrangian remaping and relaxation schemes for the double rarefaction waves problem.

Figure 5 .

 5 Figure 5. Electronic and ionic temperatures computed with HLL and relaxation schemes for the analytical test case. Left : ν ei = 0, right : ν ei = 100.

Figure 6 .

 6 Figure 6. Density, velocity and pressure computed with HLL and relaxation schemes compared with exact solution for the stationary shock problem.

Figure 7 .

 7 Figure 7. Electronic and ionic temperatures computed with HLL and relaxation schemes for the stationary shock problem. Left : ν ei = 0, right : ν ei = 100.

  R and ε * i,R . First we can note that (113) implies the positivity of ρ

* L and ρ * R . On the other hand, from a straightforward calculation, ε * e,L , ε * i,L , ε * e,R and ε * i,R are positive if condition (109) is fulfilled. Finally we can deduce that a must satisfy