F Berthelin
email: florent.berthelin@unice.fr

T Goudon

S Minjeaud
email: minjeaud@unice.fr

Kinetic schemes on staggered grids for barotropic Euler models: entropy-stability analysis

We introduce, in the one-dimensional framework, a new scheme of finite volume type for barotropic Euler equations. The numerical unknowns, namely densities and velocities, are defined on staggered grids. The numerical fluxes are defined by using the framework of kinetic schemes. We can consider general (convex) pressure laws. We justify that the density remains non negative and the total physical entropy does not increase, under suitable stability conditions. Performances of the scheme are illustrated through a set of numerical experiments.

Introduction

This work is concerned with the numerical simulation of the following system of conservation laws

∂ t ρ + ∂ x (ρV ) = 0, (1.1) ∂ t (ρV ) + ∂ x (ρV 2 + p(ρ)) = 0. (1.2)
This is the Euler model for compressible fluids (in the absence of external forces) with a barotropic equation of state: the pressure p(ρ) is a function of the density only. The unknowns are the density ρ and the velocity V of the fluid. We restrict the discussion to the one-dimension framework, but the ideas can be extended to higher dimensions. As a relevant 1 specific case we can deal with isentropic flows for polytropic ideal gases where p(ρ) = kρ γ , for some constants γ 1 and k > 0.

We set U = (ρ, J = ρV ) and the system can be recast as

∂ t U + ∂ x F (U ) = 0, with F (U ) = J , J 2 ρ + p(ρ) .
Roughly speaking, given fixed time and (homogeneous) space steps δt, δx respectively, a numerical scheme for this system reads

U k+1 j -U k j + δt δx (F k j+1/2 -F k j-1/2 ) = 0
and the cornerstone of the method relies on a suitable definition of the numerical flux

F k j+1/2
as a function of the numerical unknowns U k ℓ for certain values of ℓ, neighbouring the considered index j (say for a three points approximation flux U k j-1 , U k j , U k j+1 ). We wish to discuss numerical schemes based on the framework of the so-called Boltzmann schemes [START_REF] Coron | Numerical passage from kinetic to fluid equations[END_REF][START_REF] Deshpande | Kinetic theory based new upwind methods for inviscid compressible flows[END_REF][START_REF] Deshpande | On the Maxwellian distribution, symmetric form and entropy conservation for the Euler equations[END_REF][START_REF] Perthame | Second order Boltzmann schemes for compressible Euler equations in one and two space dimension[END_REF][START_REF] Pullin | Direct simulation methods for compressible gas flow[END_REF][START_REF] Sanders | On the origin of the 3 kiloparsec arm[END_REF]. More precisely the system is seen as the limit ǫ → 0 of the following BGK-like system

∂ t f + ξ∂ x f = 1 ǫ (M [f ] -f )
in the spirit of hydrodynamic limits which allow to derive the Euler equations from the Boltzmann (or BGK) equation, see e. g. [START_REF] Raymond | Hydrodynamic limits of the Boltzmann equation[END_REF]. Here, (t, x, ξ) → f (t, x, ξ) is a vector valued function, and the "Maxwellian" state ξ ∈ R → M [f ](ξ) ∈ R 2 is a function of the auxiliary variable ξ, which is parametrized by the zeroth moments of f in such a way that, denoting (ρ, ρV ) = f dξ, we have M [f ] dξ = (ρ, ρV ), ξM [f ] dξ = (ρV, ρV 2 + p(ρ)).

Basically, the scheme works in two steps for solving the system of conservation laws: knowing ρ k , V k , approximation of density, velocity at the discrete time kδt

• First, solve 1 δt (f ⋆ -f k ) + ξ∂ x f k = 0.
This is a mere linear transport equation.

• Second, project the solution to the equilibrium state

f k+1 = M [f ⋆ ]
the Maxwellian having the same zeroth moments as f ⋆ .

In practice, we get rid of the extra velocity variable ξ by integrating the formula with respect to ξ: it provides a scheme of Finite Volume type for updating ρ k+1 , V k+1 . Precisely, by using the basic upwind discretization in the convection step, the reasoning leads to the following definition of numerical fluxes

F k j+1/2 = ξ 0 ξM k j dξ + ξ 0 ξM k j+1 dξ, 2
where ξ → M k j (ξ) stands for the Maxwellian state associated to the moments

U k j = (ρ k j , ρ k j V k j ) = M k j dξ. It is convenient to set F k j+1/2 = F + (U k j ) + F -(U k j+1 ) with F + (U ) = ξ 0 ξM (ξ) dξ, F -(U ) = ξ 0 ξM (ξ) dξ. (1.3)
We refer the reader to [START_REF] Coron | Numerical passage from kinetic to fluid equations[END_REF][START_REF] Deshpande | Kinetic theory based new upwind methods for inviscid compressible flows[END_REF][START_REF] Deshpande | On the Maxwellian distribution, symmetric form and entropy conservation for the Euler equations[END_REF][START_REF] Perthame | Second order Boltzmann schemes for compressible Euler equations in one and two space dimension[END_REF][START_REF] Pullin | Direct simulation methods for compressible gas flow[END_REF][START_REF] Sanders | On the origin of the 3 kiloparsec arm[END_REF] for the design of such kinetic schemes in the context of gas dynamics. The analysis of such schemes is thoroughly detailed in the textbooks [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF][START_REF] Perthame | Kinetic formulation of conservation laws[END_REF]. In particular, the schemes can be reinterpreted by means of approximate Riemann solvers, see [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF]Section 2.5]. The consistency of the scheme is simply embodied into the property F + (U ) + F -(U ) = F (U ). The numerical analysis aims at exhibiting stability conditions which, at least, preserves the natural positivity of the density ρ. Another crucial issue is related to the behaviour of certain nonlinear functionals of the unknowns, the socalled entropies: admissible solutions should dissipate these quantities. For kinetic schemes, these properties are intimately connected to the design of the Maxwellian functions M . In particular, preserving the positivity of the density makes appealing the choice of equilibrium with compact support (with respect to the ξ variable). Dealing with isentropic flows for a polytropic ideal gas (p(ρ) = kρ γ ), it is possible to identify a convenient Maxwellian in order to dissipate the natural "physical" entropy of the problem. It leads to solve a minimization problem under constraints, but the effective computation relies strongly on the homogeneity property of the pressure law [START_REF] Bouchut | Construction of BGK models with a family of kinetic entropies for a given system of conservation laws[END_REF]; and it is not clear how to apply the method when dealing with intricate pressure laws, like the laws described e. g. in [START_REF] Houzé | Dissipative issue of high-order shock capturing schemes with non-convex equations of state[END_REF]. Furthermore, the resulting numerical fluxes do not have in general an explicit expression by means of the numerical unknowns, which might lead to practical difficulties (the case γ = 2 being a remarkable exception).

Here, we revisit the design of kinetic schemes for (1.1)-(1.2). Our approach differs from the standard one in the following three directions.

• Firstly, we propose a non-classical definition of the Maxwellian M . The motivation is to consider quite general pressure laws, for which it is not obvious to find dissipative equilibrium states just by solving minimization problems. However, our computation remains reminiscent of ideas in [START_REF] Kaniel | Approximation of the hydrodynamic equations by a transport process[END_REF] considering compactly supported equilibria, with a support driven by the propagation speeds of (1.1)-(1.2).

• Secondly, our version of the kinetic scheme works on staggered grids where the discrete density ρ and the material velocity V are not stored on the same location. While the approaches are completely different in spirit, the idea dates back to [START_REF] Neumann | A method for the numerical calculation of hydrodynamic shocks[END_REF] and [START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF]: it is used in industrial contexts in the framework of Lagrangian methods, see e.g. [START_REF] Houzé | Dissipative issue of high-order shock capturing schemes with non-convex equations of state[END_REF]. We also refer the reader to [START_REF] Arakawa | A potential enstrophy and energy conserving scheme for the shallow-water equations[END_REF] for the conception of such a scheme on staggered grids for the Shallow-Water system, based on a Finite Difference reasoning and the preservation of certain physical quantities. More recently, staggered strategies have been developed for gas dynamics systems in [START_REF] Herbin | Consistent explicit staggered schemes for compressible flows; part I: the barotropic Euler equations[END_REF][START_REF] Herbin | Explicit staggered schemes for the compressible Euler equations[END_REF].

• Thirdly, by contrast to the most common strategy adopted for hyperbolic systems, our scheme upwinds with the material velocity as a privileged speed, instead of using the full wave structure of the system. Although the derivation of the scheme is based on different principles, this idea appears in the so-called AUSM schemes [START_REF] Liou | A new flux splitting scheme[END_REF][START_REF] Liou | A sequel to AUSM: AUSM+[END_REF], see also [START_REF] Gallouët | Short notes: Kinetic energy control in explicit finite volume discretizations of the incompressible and compressible Navier-Stokes equations[END_REF][START_REF] Gastaldo | Staggered discretizations, pressure correction schemes and all speed barotropic flows[END_REF][START_REF] Herbin | Staggered schemes for all speed flows[END_REF][START_REF] Herbin | Consistent explicit staggered schemes for compressible flows; part I: the barotropic Euler equations[END_REF][START_REF] Herbin | Explicit staggered schemes for the compressible Euler equations[END_REF] for recent analysis of schemes in the same vein.

Our motivation is two fold. On the one hand staggered discretizations can be expected to fulfil better stability properties in low-Mach regimes, because they naturally avoid odd-even decoupling of the pressure, and the possible occurrence of spurious modes. We refer the reader to [START_REF] Gastaldo | Staggered discretizations, pressure correction schemes and all speed barotropic flows[END_REF][START_REF] Herbin | Staggered schemes for all speed flows[END_REF] and the references therein on this aspect. On the other hand, the method is well adapted to treat coupled models describing mixtures, which involve an intricate constraint on the material velocity. The staggered framework allows to design a scheme for such complex flows in order to conserve exactly the total mass of the mixture. This issue is detailed in [START_REF] Berthelin | Multifluid flows: a kinetic approach[END_REF]. This work is organized as follows. We start by recalling briefly a few basic facts about (1.1)-(1.2) in Section 2 where we set up the notation. In Section 3, we introduce the Maxwellian states on which the scheme is based. Then, we detail the space discretization. We also identify the stability condition which preserves the positivity of the density ρ. Section 4 is devoted to the analysis of the behaviour of the discrete (physical) entropy. It turns out that working on staggered grids helps in proving a dissipation property. We establish both local and global inequalities, the latter being not direct consequences of the former. Finally, in Section 5 we discuss a few numerical simulations, dealing either with the simple equation of state p(ρ) = kρ γ , or with more intricate pressure laws. It shows that the method is efficient and reliable.

Basic facts on the system of conservation laws

In all what follows, we assume that the pressure function p satisfies

p ∈ C 2 [0, ∞) , p(ρ) > 0, p ′ (ρ) > 0, p ′′ (ρ) > 0, ∀ρ > 0.
(2.4)

Consequently, the sound speed c : ρ -→ p ′ (ρ) is well defined; furthermore c : ρ → c(ρ) is strictly increasing. This property plays a central role in the analysis of the proposed numerical fluxes (see Lemma 3.2 below). Assumption (2.4) also implies that the function τ → p( 1 τ ) is strictly convex. This property enters in the analysis of the invariant regions of the system (1.1)-(1.2), but it will not appear in the analysis of the scheme. Note that non-convex pressure laws arise in real-life applications, and it leads to specific difficulties, see [START_REF] Houzé | Dissipative issue of high-order shock capturing schemes with non-convex equations of state[END_REF].

At least formally, the system (1.1)-(1.2) can be rewritten in the non-conservative form

∂ t U + A(U )∂ x U = 0, A(U ) = ∇ U F (U ),
where A is the jacobian matrix of the flux function

F : (ρ, J = ρV ) = U → (J , J 2 /ρ + p(ρ)), namely A(U ) = 0 1 p ′ (ρ) -J 2 /ρ 2 2J /ρ .
The eigenvalues of this matrix define the characteristic speeds of the system (1.1)-(1.2):

λ -(U ) = V -c(ρ), λ + (U ) = V + c(ρ).
Therefore for any ρ > 0, V ∈ R, A(U ) admits two distinct real eigenvalues, and the system (1.1)-(1.2) is hyperbolic.

In view of its physical meaning the density ρ is expected to remain non negative. Next,

let us set W ± (U ) = V ± G(ρ) where G ′ (ρ) = √ p ′ (ρ) ρ = c(ρ) ρ .
As far as the solution of (1.1)-(1.2) is smooth, these quantities, the so-called Riemann invariants, are simply advected at the speed λ ± (U ): we can check that (∂

t + λ ± ∂ x )W ± = 0. Owing to (2.4), ρ → ρG(ρ) is convex, which is equivalent to d dρ (ρ 2 G ′ (ρ)) = d dρ (ρ p ′ (ρ)) 0.
Then, for any κ ∈ R, the sets {U = (ρ, ρV ), V + G(ρ) κ} and {U = (ρ, ρV ), V -G(ρ) κ} are convex and they are left invariant by the dynamics, see for instance [START_REF] Serre | Systèmes de lois de conservation[END_REF]Th. 8.3.8]. This observation provides uniform estimates on the solutions of (1.1)-(1.2). Accordingly, we can deduce L ∞ estimates on the density ρ and the velocity V , by means of the initial data. We refer for instance to [START_REF] Chueh | Positively invariant regions for systems of nonlinear diffusion equations[END_REF] for the analysis of such invariant sets through viscous approximations. However showing that a numerical scheme preserves these natural estimates is far from obvious: it can be justified for Godunov's scheme, which is based on the exact resolution of Riemann problems, or Lax-Friedrichs' schemes, see e. g. [START_REF] Roux | Numerical stability for some equations of gas dynamics[END_REF], and these estimates are then the first step towards the analysis of the existence of solutions to (1.1)-(1.2), see [START_REF] Perna | Convergence of approximate solutions to conservation laws[END_REF][START_REF] Ding | Convergence of the fractional step Lax-Friedrichs scheme and Godunov scheme for the isentropic system of gas dynamics[END_REF].

We set Φ : ρ -→ Φ(ρ) such that ρΦ ′ (ρ) -Φ(ρ) = p(ρ).
Note that Φ ′′ (ρ) = p ′ (ρ) ρ , so that Φ is convex. Smooth solutions of (1.1)-(1.2) can be shown to satisfy the local balance law

∂ t ρ V 2 2 + Φ(ρ) + ∂ x ρ V 2 2 + Φ(ρ) + p(ρ) V = 0.
It motivates an admissibility criterion to select among weak solutions. In turn, they satisfy the following global entropy inequality (for suitable boundary conditions)

d dt ρ V 2 2 + Φ(ρ) dx 0 (2.5)
indicating that entropy is dissipated in the admissible discontinuities of the solutions of (1.1)-(1.2). It is therefore an issue to determine whether or not a numerical scheme produces solutions that satisfy the entropy criterion. In order to proceed with the analysis of the scheme we propose, it is worth having in mind how the computation works at the continuous level. On the one hand, we check that

∂ t Φ(ρ) + ∂ x Φ(ρ)V = Φ(ρ) -ρΦ ′ (ρ) ∂ x V = -p(ρ)∂ x V. (2.6)
On the other hand, for the kinetic energy we get

∂ t ρ V 2 2 + ∂ x ρ V 2 2 V = ∂ t (ρV ) + ∂ x (ρV 2 ) V = -∂ x p(ρ) V. (2.7)
Adding the two relations yields the local relation

∂ t ρ V 2 2 + Φ(ρ) + ∂ x ρ V 2 2 + Φ(ρ) + p(ρ) V = 0 (2.8)
and the conclusion follows by integrating and using the boundary conditions. Of course time and space discretizations break this structure and the challenge consists in identifying discrete version of the derivatives in (2.6) and (2.7).

3 Definition of the kinetic scheme on staggered grids

Maxwellian states

The kinetic scheme we shall study is based on the following definition

M 0 (ρ, V, ξ) = ρ 2c(ρ) 1I |ξ-V | c(ρ) , (3.9) M 1 (ρ, V, ξ) = V M 0 (ρ, V, ξ) + M (ρ, V, ξ) with M (ρ, V, ξ) j = ξL(ρ, V )1I |ξ| |V |+c(ρ) , L(ρ, V ) = 3 2 (|V | + c(ρ)) -3 p(ρ).
We remind that c(ρ) stands for the sound speed. In particular the definition of M 0 is reminiscent of Kaniel's strategy for gas dynamics [START_REF] Kaniel | Approximation of the hydrodynamic equations by a transport process[END_REF] with a support of the equilibrium that exactly contains all the velocities in the interval [V -c(ρ), V + c(ρ)] defined by the fundamental wave speeds of the system (1.1)-(1.2). From now on we adopt a slight abuse of notation for λ ± compared to Section 2, denoting λ ± (ρ, V ) = V ± c(ρ). With U = (ρ, ρV ), we set

F ± (U ) = ξ≷0 ξ M 0 (ρ, V, ξ) M 1 (ρ, V, ξ) dξ.
Consistency of the corresponding flux-splitting method is a consequence of the following claim.

Proposition 3. [START_REF] Arakawa | A potential enstrophy and energy conserving scheme for the shallow-water equations[END_REF] We have

(M 0 , M 1 )(ρ, V, ξ) dξ = (ρ, ρV ), ξ(M 0 , M 1 )(ρ, V, ξ) dξ = (ρV, ρV 2 + p(ρ)).
The flux decomposition is consistent since

F + (U ) + F -(U ) = F (U ) holds.
The definition of the Maxwellian (M 0 , M 1 ) is very specific. In the momentum fluxes, the convection terms are dictated by the mass fluxes, while the pressure term has a very simple expression. This remark, which is crucial in the analysis of the scheme, is made clear through the following statement.

Lemma 3.1 For any ρ > 0, V ∈ R, we have ξ≷0 ξ M (ρ, V, ξ) dξ = p(ρ) 2 and ξ≷0 ξM 1 (ρ, V, ξ) dξ = V F ± (ρ, V ) + p(ρ) 2 .
For further purposes, it is convenient to introduce the following mappings

F ± : [0, ∞) × R -→ R (ρ, V ) -→ ξ≷0 ξM 0 (ρ, V, ξ) dξ. (3.10)
Next, we set

F (ρ, V ) = F + (ρ, V ) + F -(ρ, V ) = ξM 0 (ρ, V, ξ) dξ = ρV, F |•| (ρ, V ) = F + (ρ, V ) -F -(ρ, V ) = |ξ|M 0 (ρ, V, ξ) dξ 0. (3.11)
Figure 1 provides the graphs of V → F ± (ρ, V ) for a fixed ρ > 0. In particular, we observe the difference with the simple UpWind definition of the mass flux F UpW,± (ρ, V ) = ±ρ[V ] ± . The analysis of the scheme relies on the properties of the functions F ± , that we collect in

V c(ρ) -c(ρ) F + (ρ, V ) upwind F -(ρ, V ) upwind Figure 1: Graphs of V → F + (ρ, V ) and V → F -(ρ, V ) for a fixed ρ > 0.
the following claim. Lemma 3.2 Assume (2.4). Then, the functions F + , F -, F |•| satisfy the following properties:

(i) 0 F + (ρ, V ) ρ λ + (ρ, V ) + , ∀V ∈ R, ∀ρ 0. (ii) -ρ λ -(ρ, V ) - F -(ρ, V ) 0, ∀V ∈ R, ∀ρ 0. (iii) (ρ, V ) ∈ (0, ∞) × R → F ± (ρ, V ) are C 1 functions and F + (ρ, -V ) = -F -(ρ, V ). (iv) ρ ∈ [0, ∞) → F + (ρ, V ) is increasing for V ∈ R, and strictly increasing for V 0. (v) ρ ∈ [0, ∞) → F -(ρ, V ) is decreasing for V ∈ R, and strictly decreasing for V 0. (vi) ρ ∈ [0, ∞) → F |•| (ρ, V ) 0 is strictly increasing for V ∈ R,
Proof. The conclusion follows by direct inspection of the following formula, where the expression of F + (ρ, V ) changes depending on the Mach number V /c(ρ):

F + (ρ, V ) =            0 if V + c(ρ) 0, ρ 4c(ρ) (V + c(ρ)) 2 = ρ 4c(ρ) λ + (ρ, V ) 2 if V -c(ρ) < 0 < V + c(ρ), ρV = ρ 4c(ρ) (λ + (ρ, V ) 2 -λ -(ρ, V ) 2 ) if 0 < V -c(ρ).
The estimates can be obtained by comparing λ + (ρ, V ) and 4c(ρ) when we consider separately each case. Similarly, we have

F -(ρ, V ) =            ρV = ρ 4c(ρ) (λ + (ρ, V ) 2 -λ -(ρ, V ) 2 ) if V + c(ρ) 0, - ρ 4c(ρ) (V -c(ρ)) 2 = - ρ 4c(ρ) λ -(ρ, V ) 2 if V -c(ρ) < 0 < V + c(ρ), 0 if 0 < V -c(ρ).

Remark 3.3

The explicit expression of the functions F ± are no longer used in the sequel. We only use the properties stated in Lemma 3.2 so that the result could be directly applied to another flux splitting satifying all these properties. More precisely, properties (i) and (ii) are used to derive CFL conditions ensuring the positivity of the discrete density ρ (or more generally, guaranteeing non uniform bound from above and below for the density ρ). Properties (iii)-(vi) are used in the analysis of the entropy-stability.

Staggered grids

From now on, we consider that the problem (1.1)-(1.2) holds on the bounded domain (0, L).

It is completed by the boundary condition

V (t, 0) = 0 = V (t, L), ∀t > 0. (3.12)
Note that it belongs to the framework designed in [START_REF] Dubois | Boundary conditions for nonlinear hyperbolic systems of conservation laws[END_REF] for the analysis of initial boundary value problems for systems of conservation laws. Roughly speaking, only one field is incoming and we do not need further boundary data. Notice that with this boundary condition the entropy balance (2.5) holds.

We now wish to discuss the adaptation of the kinetic scheme on staggered grids. While we present the framework in the one dimension case, the method can be adapted to higher dimensions. We consider a set of J + 1 points of the computational domain [0, L] such that x 1 = 0 < x 2 < ... < x J < x J+1 = L. These points define a subdivision of [0, L] which is called the primal mesh. The x j 's are referred to as the edges (or vertices) of the primal mesh. We set δx j+1/2 = x j+1 -x j and x j+1/2 = 1 2 (x j + x j+1 ) for j ∈ {1, ..., J}. The J points x 3/2 < ... < x J+1/2 are the centers of the primal cells and realize the dual mesh. We set δx 1 = δx 3/2 /2, δx J = δx J+1/2 /2 and δx j = 1 2 (δx j-1/2 + δx j+1/2 ) for j ∈ {2, ..., J -1}. For the numerical unknowns:

• Densities are evaluated at centers of the primal mesh : ρ j+1/2 , with j ∈ {1, ..., J},

• Velocities are evaluated at edges of the primal mesh : V j with j ∈ {1, ..., J + 1}.

The density is updated with a Finite Volume approximation on the primal mesh, which thus requires an approximation of the fluxes ρV at the interfaces x = x j . Namely, we have

δx j+1/2 δt (ρ k+1 j+1/2 -ρ k j+1/2 ) + F k j+1 -F k j = 0, ∀j ∈ {1, .., J}. (3.13)
We adopt the kinetic scheme. For internal edges, quite naturally, we use the value of the velocity at the edge x = x j and we upwind the density:

F k j = F + (ρ k j-1/2 , V k j ) + F -(ρ k j+1/2 , V k j ), ∀j ∈ {2, . . . , J}. (3.14)
For external edges, since we use homogeneous boundary conditions (3.12), the fluxes are set to zero:

F k 1 = 0, F k J+1 = 0. (3.15)
We proceed similarly to define a Finite Volume approximation of the momentum equation on the cells (x j-1/2 , x j+1/2 ) of the dual mesh. We introduce an approximation of ρ at the edges of the primal mesh

ρ k j = δx j+1/2 ρ k j+1/2 + δx j-1/2 ρ k j-1/2 2δx j , ∀j ∈ {2, . . . , J}. (3.16)
Since at the kth time iteration, the approximate density is seen as the piecewise constant function J j=1 ρ k j+1/2 1 [x j ,x j+1 [ , ρ k j is nothing but the mean value of the density on the cell (x j-1/2 , x j+1/2 ). We update the velocity with

δx j δt (ρ k+1 j V k+1 j -ρ k j V k j ) + F k j+1/2 -F k j-1/2 = 0, ∀j ∈ {2, . . . , J}. (3.17)
Since we want to impose boundary conditions (3.12), we set

V k+1 1 = V k+1 J+1 = 0. (3.18)
The fluxes F k j+1/2 are also defined by the kinetic scheme. As remarked in Lemma 3.1, M only contributes to the pressure. More precisely, it yields a centered difference of the pressure term since for any ρ, V, V ′ we have

ξ>0 ξ M (ρ, V ) dξ + ξ<0 ξ M (ρ ′ , V ′ ) dξ = 1 2 (p(ρ) + p(ρ ′ )).
This term is evaluated at the center x = x j+1/2 , it is thus natural to make use of the available value of the density at this point: the corresponding contribution to the flux F k j+1/2 is therefore p(ρ k j+1/2 ). The convection flux is given by an approximation of ξV M 0 (ρ, V ) dξ at x = x j+1/2 . We mimic the formula obtained with the mass flux: we upwind the quantity ρV , which is advected by the velocity V . The latter is evaluated at the interface. We are thus led to

F k j+1/2 = ξ>0 ξV k j M 0 (ρ k j , V k j+1/2 , ξ) dξ + ξ<0 ξV k j+1 M 0 (ρ k j+1 , V k j+1/2 , ξ) dξ + p(ρ k j+1/2 )
where we need to make ρk j , ρk j+1 and V k j+1/2 precise. Instead of using the basic interpolation ρ k j , we bear in mind that ξ>0 ξM 0 (ρ k j , V k j+1/2 , ξ) dξ represents the mass flux going from left to right through the interface located at x = x j+1/2 . We evaluate it as the average of the (already known) mass fluxes from left to right at the interfaces x = x j and x = x j+1 . Reasoning the same way with the mass flux going from right to left (ξ < 0) yields

F k j+1/2 = V k j 2 F + (ρ k j-1/2 , V k j ) + F + (ρ k j+1/2 , V k j+1 ) + V k j+1 2 F -(ρ k j+1/2 , V k j ) + F -(ρ k j+3/2 , V k j+1 ) + p(ρ k j+1/2 ), ∀j ∈ {2, ..., J -1}. ( 3.19)
For j = 1, (resp. j = J) we remind that V k 1 = 0 (resp. V k J+1 = 0) so that the contribution associated to the positive (resp. negative) ξ's vanishes. Hence only the mass flux from right to left at the interface x 3/2 (resp. the mass flux from left to right at the interface x J-/2 ) has to be considered. Since there is no mass flux at x = 0 (resp. x = L), we arrive at

F k 3/2 = V k 2 2 F -(ρ k 5/2 , V k 2 ) + p(ρ k 5/2 ), F k J+1/2 = V k J 2 F + (ρ k J-1/2 , V k J ) + p(ρ k J-1/2 ).
Due to the very specific form of the Maxwellian M 1 , the scheme treats differently inertia and pressure, in the spirit of AUSM schemes [START_REF] Liou | A new flux splitting scheme[END_REF][START_REF] Liou | A sequel to AUSM: AUSM+[END_REF]. As pointed out in the Introduction, the definition of the numerical mass and momentum fluxes does not involve the resolution of Riemann problems, nor the computation of intricate integrals, that could be quite costly.

Remark 3.4 A simple variant of the scheme is obtained by replacing

M 0 (ρ, V, ξ) by ρδ(ξ = V )
in the definition of the mass fluxes, in the spirit of [START_REF] Sanders | On the origin of the 3 kiloparsec arm[END_REF]. With such a definition of the equilibrium state, convection terms are UpWinded, and the pressure is approached by a centered difference. For the sake of simplicity, let us assume that we are working with a uniform mesh with constant mesh size δx > 0.

With [V ] ± = 1 2 (|V | ± V ) = max(±V, 0), the scheme reads ρ k+1 j+1/2 = ρ k j+1/2 - δt δx F k,UpW j+1 -F k,UpW j , with F k,UpW j = ρ k j-1/2 V k j + -ρ k j+1/2 V k j -. ρ k+1 j V k+1 j = ρ k j V k j - δt δx G k,UpW j+1/2 -G k,UpW j-1/2 - δt δx p(ρ j+1/2 ) -p(ρ j-1/2 ) ,
where

ρ k j = 1 2 (ρ k j+1/2 + ρ k j-1/2 ) and G k,UpW j+1/2 = V k j 2 ρ k j-1/2 V k j + + ρ k j+1/2 V k j+1 + - V k j+1 2 ρ k j+1/2 V k j -+ ρ k j+3/2 V k j+1 -.
This scheme is very close to the ones proposed and analysed in [START_REF] Herbin | Consistent explicit staggered schemes for compressible flows; part I: the barotropic Euler equations[END_REF][START_REF] Herbin | Explicit staggered schemes for the compressible Euler equations[END_REF] which are also constructed on staggered grids with UpWinding strategies based on the material velocity. It produces consistent results but with oscillations, that remain of controlled amplitude, when the velocity vanishes. It will be further discussed in Appendix B.

Remark 3.5 Note that the numerical fluxes naturally incorporates the boundary condition (3.12). In particular, since the density is evaluated only on "interior points", we do not need any ghost cells to treat the pressure gradient.

We conclude this Section with the following stability statement, which gives conditions in order to produce non negative densities. Proposition 3.2 Assume that the initial data satisfies ρ 0 j+1/2 0 for any j ∈ {1, . . . , J}. We assume the CFL-like condition

δt δx j+1/2 λ + (ρ k j+1/2 , V k j+1 ) + + λ -(ρ k j+1/2 , V k j ) - 1, ∀j = 1, . . . , J, (3.20)
at every time step. Then the scheme preserves the positivity of ρ.

Proof. We rewrite the evolution of the discrete density as follows

ρ k+1 j+1/2 = ρ k j+1/2 + δt δx j+1/2 F -(ρ k j+1/2 , V k j ) -F + (ρ k j+1/2 , V k j+1 ) + δt δx j+1/2 F + (ρ k j-1/2 , V k j ) -F -(ρ k j+3/2 , V k j+1 ) .
By assumption the components ρ k ℓ+1/2 are non negative for any ℓ, and the contribution of the last two terms is non negative. Next, Lemma 3.2 tells us that

ρ k+1 j+1/2 ρ k j+1/2 1 - δt δx j+1/2 λ + (ρ k j+1/2 , V k j+1 ) + + λ -(ρ k j+1/2 , V k j ) - .
Hence ρ k+1 j+1/2 0 when (3.20) is fulfilled.

Remark 3.6

Remark that the stability condition involves the characteristic speeds V ± c(ρ), and not the material velocity V only. This is by contrast to the schemes which use UpWinding strategies based on the material velocity, see Remark 3.4 and Appendix B for further comments.

Stability analysis: entropy dissipation

In this Section, we wish to establish a discrete analog of (2.5). Let us set up a few notation. We set (for j ∈ {2, ..., J}):

ρ k j,min := min ρ k j-1/2 , ρ k j+1/2 , ρ k j,Max := max ρ k j-1/2 ,
ρ k j+1/2 .
In order to avoid technical difficulties due to the presence of vacuum, we shall assume that ρ k min > 0. Note however that the scheme performs well in vacuum formation, see Test 5 in Section 5. In fact, we are going to analyze a slightly modified version of the scheme (3.13)- (3.19). According to (3.19), we split the momentum flux as follows:

F k j+1/2 = G k j+1/2 + π k+1/2 j+1/2 , j ∈ {1, . . . , J}, (4.21)
with

G k 3/2 = V k 2 2 F -(ρ k 5/2 , V k 2 ), G k j+1/2 = V k j 2 F + (ρ k j-1/2 , V k j ) + F + (ρ k j+1/2 , V k j+1 ) + V k j+1 2 F -(ρ k j+1/2 , V k j ) + F -(ρ k j+3/2 , V k j+1 ) , j ∈ {2, . . . , J -1}, G k J+1/2 = V k J 2 F + (ρ k J-1/2 , V k J ),
and

π k+1/2 j+1/2 = ρ k j+1/2 Φ ′ (ρ k+1 j+1/2 ) -Φ(ρ k j+1/2 ), j ∈ {1, . . . , J}. (4.22)
The flux G k j+1/2 represents the contribution of the inertial terms whereas the term π k+1/2 j+1/2

discretizes the pressure forces. The expression of the latter relies on the relation p(ρ) = ρΦ ′ (ρ) -Φ(ρ). We have thus replaced p(ρ k j+1/2 ) in (3.19) by π k+1/2 j+1/2 . It has the flavor of an implicit relation; however, we should bear in mind that the density is updated by (3.13) before computing the velocity and we thus have ρ k+1 j+1/2 at hand without the need of an intricate fixed point method. The motivation of this modification will appear clearly in Section 4.1. Throughout this section, for some k ∈ N, the state (ρ k j+1/2 ) j∈{1,...,J} and (V k j ) j∈{1,...,J+1} is given and we suppose that the following set of assumptions is fulfilled: We shall prove the following discrete global entropy inequality, which can be seen as a stability statement.

(h1) The pressure function p satisfy (2.4), (h2) ρ k min > 0, (h3) V k 1 = V k J+1 = 0, (

Theorem 4.1 Assume (h1)-(h4).

Then, there exists τ ⋆ > 0 such that for any 0 < δt < τ ⋆ , the updated state verifies

1 δt J j=2 δx j 1 2 ρ k+1 j V k+1 j 2 + J j=1 δx j+1/2 Φ(ρ k+1 j+1/2 ) 1 δt J j=2 δx j 1 2 ρ k j V k j 2 + J j=1 δx j+1/2 Φ(ρ k j+1/2 ) .
The time step τ ⋆ depends only on the state ρ k , V k , on the parameters of the space discretization, and on the properties of the pressure law; its identification relies on the combination of quite intricate but explicit stability conditions (see Corollary 4.12 for a precise statement). By contrast to standard proofs in hyperbolic theory, we consider separately, in the two next sections, the evolution of the internal energy Φ(ρ) (Section 4.1) and of the kinetic energy ρV 2 /2 (Section 4.2). Roughly speaking, the evolution of the kinetic energy is obtained by multiplying the momentum equation by V . It can be split into two contributions: the work of the pressure forces (see (2.6)) and the contribution of the inertial terms (see (2.7)). We shall adopt the same splitting at the discrete level and the modification of the scheme will be useful in order to compensate the work of the pressure forces with a similar contribution coming from the evolution of the internal energy. We proceed into two steps. We start by establishing local estimates, namely discrete versions of (2.6) and (2.7) in Sections 4.1 and 4.2 respectively. The former holds on the primal mesh, the latter on the dual cells. Unfortunately, neither a local nor a global entropy inequality (that is the discrete analog of (2.8) or (2.5) as stated in Theorem 4.1) can be obtained directly by summing the two local inequalities. It requires a further, quite intricate, analysis presented in Section 4.3. Note that local estimates are interesting in itself, in view of a full consistency analysis. Indeed, it is likely that we can adapt arguments from [START_REF] Herbin | Consistent explicit staggered schemes for compressible flows; part I: the barotropic Euler equations[END_REF][START_REF] Herbin | Explicit staggered schemes for the compressible Euler equations[END_REF] in order to establish a Lax-Wendrofflike statement for this scheme, which would prove that the limit of a converging sequence of stepwise constant functions defined from the scheme is a weak solution of the system of conservation laws that satisfies the entropy inequality. This question is however beyond the scope of the present paper.

Evolution of the internal energy

This Section is devoted to the proof of the discrete analog of equality (2.6). The following notation will be useful:

d ± (ρ 1 , ρ 2 , V ) =      F ± (ρ 1 , V ) -F ± (ρ 2 , V ) ρ 1 -ρ 2 if ρ 1 = ρ 2 , ∂ ρ F ± (ρ 1 , V ) if ρ 1 = ρ 2 , ( 4.23)
and

d |•| (ρ 1 , ρ 2 , V ) =      F |•| (ρ 1 , V ) -F |•| (ρ 2 , V ) ρ 1 -ρ 2 if ρ 1 = ρ 2 , ∂ ρ F |•| (ρ 1 , V ) if ρ 1 = ρ 2 . (4.24) Owing to Lemma 3.2, d + (ρ 1 , ρ 2 , V ) is non negative for any V ∈ R, positive when V 0, while d -(ρ 1 , ρ 2 , V ) is non positive for any V ∈ R, negative when V 0. Consequently d |•| = d + -d -is always positive.
Let us also introduce the two following quantities

ρ k+1 j+1/2 = ρ k j+1/2 - 2δt δx j+1/2 F + (ρ k j+1/2 , V k j ) -F + (ρ k j-1/2 , V k j ) , ∀j ∈ {2, . . . , J}, ρ k+1 j-1/2 = ρ k j-1/2 - 2δt δx j-1/2 F -(ρ k j+1/2 , V k j ) -F -(ρ k j-1/2 , V k j ) , ∀j ∈ {2, . . . , J}. (4.25)
and, at the boundary,

ρ k+1 J+1/2 = ρ k J+1/2 , ρ k+1 3/2 = ρ k 3/2
. By combining the equality F + (ρ, V ) + F -(ρ, V ) = ρV and the discrete mass balance (3.13), we readily obtain the following equality

ρ k+1 j+1/2 + ρ k+1 j+1/2 2 = ρ k+1 j+1/2 + δt δx j+1/2 ρ k j+1/2 V k j+1 -V k j , ∀j ∈ {1, . . . , J}. (4.26)
We can now state the main result of this section.

Proposition 4.2 Assume (h1)-(h4). Then, for any

δt > 0 such that δt < min ( ρ k j,min , ρ k j,Max ) Φ ′′ max ( ρ k j,min , ρ k j,Max ) Φ ′′ min δx j-1/2 , δx j+1/2 4d |•| (ρ k j-1/2 , ρ k j+1/2 , V k j ) , ∀j ∈ {2, . . . , J}, (4.27)
the following inequality holds

δx j+1/2 δt Φ(ρ k+1 j+1/2 ) -Φ(ρ k j+1/2 ) + G k j+1 -G k j + π k+1/2 j+1/2 V k j+1 -V k j 0, ∀j ∈ {1, . . . , J},
where

G k 1 = 0 and G k j = Φ(ρ k j-1/2 )V k j - δx j-1/2 2δt Φ ρ k+1 j-1/2 -Φ(ρ k j-1/2
) , ∀j ∈ {2, . . . , J + 1}.

Remark 4.3

The flux G k j is defined as a function of

ρ k j-1/2 , ρ k j-1/2 and V k j G k j = G(ρ k j-1/2 , ρ k j-1/2 , V k j ).
The flux is consistent in the sense that G(ρ, ρ, V ) = Φ(ρ)V . Suprisingly, the function G depends on the ratio δt δx j-1/2 . It is harmless if this ratio remains bounded away from zero but it could become an issue in the semi-discrete limit when the time step δt tends to zero (independently of the mesh size). Nevertheless, formally, assuming that the density and the velocity converge, say

ρ k j±1/2 -→ δt→0 ρ j±1/2 and V k j -→ δt→0 V j , we obtain G k j -→ δt→0 Φ(ρ j-1/2 )V j + F -(ρ j+1/2 , V j ) -F -(ρ j-1/2 , V j ) Φ ′ ρ j-1/2 ,
which is still a consistent discretization of the internal energy flux.

The proof of Proposition 4.2 makes use of the following technical result.

Lemma 4.4 Let Φ be a strictly convex function of class

C 2 . Let ρ 1 , ρ 2 > 0, V ∈ R. We denote ρ = min(ρ 1 , ρ 2 ), ρ = max(ρ 1 , ρ 2 ). Let λ, µ ∈ R verify λ, µ > 2 max ( ρ, ρ ) Φ′′ min ( ρ, ρ ) Φ′′ d |•| (ρ 1 , ρ 2 , V ). ( 4

.28)

Let us set

ρ 1 = ρ 1 - 1 µ F -(ρ 2 , V ) -F -(ρ 1 , V ) , ρ 2 = ρ 2 - 1 λ F + (ρ 2 , V ) -F + (ρ 1 , V ) .
Then, there holds

Φ(ρ 2 )-Φ(ρ 1 ) V +λ Φ ρ 2 -Φ(ρ 2 ) +µ Φ ρ 1 -Φ(ρ 1 ) + 1 4 |V | + d |•| 2 min ( ρ, ρ ) Φ′′ (ρ 1 -ρ 2 ) 2 0. (4.29)
We postpone to Appendix A the details of the proof; the arguments rely on the convexity of the function Φ and on the properties of the fluxes F + and F -stated in Lemma 3.2.

= V /λ if V > 0 and α = -V /µ if V < 0) Φ(ρ 2 ) -Φ(ρ 1 ) + 1 α Φ ρ 2 -α(ρ 2 -ρ 1 ) -Φ(ρ 2 ) + 1 4 min ( ρ, ρ ) Φ′′ (ρ 1 -ρ 2 ) 2 0.
It holds when 0 < α < min

( ρ, ρ ) Φ′′ /(2max ( ρ, ρ ) Φ′′ ).
We now go back to the proof of Proposition 4.2.

Proof of Proposition 4.2. For j ∈ {2, . . . , J}, we apply Lemma 4.4 with

ρ 1 = ρ k j-1/2 , ρ 2 = ρ k j+1/2 , V = V k j , λ = δx j+1/2 2δt , µ = δx j-1/2
2δt . Note that assumption (4.27) on the time step ensures that (4.28) is satisfied. Hence, for all j ∈ {2, . . . , J}, we have,

Φ(ρ k j+1/2 ) -Φ(ρ k j-1/2 ) V k j + δx j+1/2 2δt Φ ρ k+1 j+1/2 -Φ(ρ k j+1/2 ) + δx j-1/2 2δt Φ ρ k+1 j-1/2 -Φ(ρ k j-1/2 ) + T 1,j 0, (4.30) with T 1,j = 1 4 |V k j | + d |•| (ρ k j-1/2 , ρ k j+1/2 , V k j ) 2 min ( ρ k j,min , ρ k j,Max ) Φ ′′ (ρ k j+1/2 -ρ k j-1/2 ) 2 .
The term T 1,j is non negative. Inequality (4.30) ensures that, for all j ∈ {2, . . . , J},

G k j+1 -G k j -Φ(ρ k j+1/2 )(V k j+1 -V k j )+ δx j+1/2 δt   Φ ρ k+1 j+1/2 + Φ ρ k+1 j+1/2 2 -Φ(ρ k j+1/2 )   0. (4.31)
Moreover, it is easy to see that the left hand side vanishes when j = 1. We now end the proof by using the convexity of Φ again and (4.26) which yields, for all j ∈ {1, . . . , J},

Φ ρ j+1/2 + Φ ρ j+1/2 2 Φ ρ k+1 j+1/2 + δt δx j+1/2 ρ k j+1/2 V k j+1 -V k j Φ ρ k+1 j+1/2 + δt δx j+1/2 ρ k j+1/2 Φ ′ ρ k+1 j+1/2 V k j+1 -V k j . (4.32)
We conclude by combining (4.31) (which holds for any j ∈ {1, . . . , J}) and (4.32).

Evolution of the kinetic energy

We turn to the discrete analog of the evolution equation (2.7) of the kinetic energy. We start by observing that the mean density ρ k j (defined by (3.16)) still satisfies a discrete conservation law.

Lemma 4.6 Let us set

F k j+1/2 = F k j + F k j+1 2
, ∀j ∈ {1, . . . , J}.

Then, we have

δx j δt (ρ k+1 j -ρ k j ) + F k j+1/2 -F k j-1/2 = 0, ∀j ∈ {2, . . . , J -1}. (4.33)
Furthermore, the relation

G k j+1/2 = V k j F k,+ j+1/2 + V k j+1 F k,- j+1/2 = 1 2 (V k j + V k j+1 )F k j+1/2 + 1 2 (V k j -V k j+1 )F k,|•| j+1/2 (4.34)
holds for all j ∈ {1, . . . , J} with

F k,+ 3/2 = 1 2 F + (ρ k 3/2 , V k 2 ) 0, F k,- 3/2 = 1 2 F -(ρ k 5/2 , V k 2 ) 0, F k,+ j+1/2 = 1 2 F + (ρ k j+1/2 , V k j+1 ) + F + (ρ k j-1/2 , V k j ) 0, ∀j ∈ {2, . . . , J -1}, F k,- j+1/2 = 1 2 F -(ρ k j+3/2 , V k j+1 ) + F -(ρ k j+1/2 , V k j ) 0, ∀j ∈ {2, . . . , J -1}, F k,+ J+1/2 = 1 2 F + (ρ k J-1/2 , V k J ) 0, F k,- J+1/2 = 1 2 F -(ρ k J+1/2 , V k J ) 0, F k,|•| j+1/2 = F k,+ j+1/2 -F k,- j+1/2
0, ∀j ∈ {1, . . . , J}.

Proof. We just use the definitions (3.16) and (3.13) of ρ k j and ρ k+1 j+1/2 to obtain

δx j δt (ρ k+1 j -ρ k j ) = 1 δt δx j+1/2 2 (ρ k+1 j+1/2 -ρ k j+1/2 ) + δx j-1/2 2 (ρ k+1 j-1/2 -ρ k j-1/2 ) = - 1 2 (F k j+1 -F k j ) + (F k j -F k j-1 ) = -(F k j+1/2 -F k j-1/2 ).
We split the expression of F k j+1/2 into positive and negative contributions

F k j+1/2 = F k,+ j+1/2 + F k,- j+1/2 , ∀j ∈ {1, . . . , J}.
Finally, the right hand side in (4.34) reads

V k j F k,+ j+1/2 + V k j+1 F k,- j+1/2 = V k j 2 F + (ρ k j+1/2 , V k j+1 ) + F + (ρ k j-1/2 , V k j ) + V k j+1 2 F -(ρ k j+3/2 , V k j+1 ) + F -(ρ k j+1/2 , V k j )
and we recognize the convection terms given by the momentum flux G k j+1/2 , as defined in (4.21). The statement makes a clear connection appear between the mass fluxes F k j+1/2 and the momentum fluxes G k j+1/2 , and it brings out the role of upwinding. With this Lemma at hand, we can establish the evolution equation for the discrete kinetic energy.

Proposition 4.7 Assume (h1)-(h4). Then, for any δt > 0 such that

δx j δt ρ k j F k,|•| j+1/2 + F k,|•| j-1/2 , (4.35)
the following inequality holds, for all j ∈ {2, . . . , J},

δx j 2δt ρ k+1 j V k+1 j 2 -ρ k j V k j 2 + Γ k j+1/2 -Γ k j-1/2 + π k+1/2 j+1/2 -π k+1/2 j-1/2 V k+1 j + δx j 4δt ρ k+1 j V k+1 j -V k j 2 0,
where the quantities

F k,|•| j-1/2 , F k,+ j-1/
Γ k j+1/2 = 1 2 V k j V k j+1 F k j+1/2 + 1 2 (V k j -V k j+1 ) 2 F k,|•| j+1/2 , ∀j ∈ {1, . . . , J}.
Proof. While the upwinding strategy is quite different (see Appendix B), the proof of Proposition 4.7 is inspired from [START_REF] Gallouët | Short notes: Kinetic energy control in explicit finite volume discretizations of the incompressible and compressible Navier-Stokes equations[END_REF]. Let j ∈ {2, . . . , J}. At the discrete level, the inertial terms are defined as follows:

C j = ρ k+1 j V k+1 j -ρ k j V k j δt + G k j+1/2 -G k j-1/2 δx j .
Multiplying the discrete momentum equation by V k+1 j we are led to

δx j C j + π k+1/2 j+1/2 -π k+1/2 j-1/2 V k+1 j = 0. (4.36)
Let us split the convection term to estimate as follows

δx j C j V k+1 j = T 1 + T 2 (4.37) with T 1 = δx j δt ρ k+1 j V k+1 j -ρ k j V k j V k+1 j and T 2 = G k j+1/2 -G k j-1/2 V k+1 j .
We begin our study with the flux term T 2 , rewritten as

T 2 = G k j+1/2 -G k j-1/2 V k j T 2,1 + G k j+1/2 -G k j-1/2 V k+1 j -V k j T 2,2
.

In T 2,1 we make use of (4.34) and we obtain

T 2,1 = 1 2 (V k j + V k j+1 )F k j+1/2 -(V k j-1 + V k j )F k j-1/2 V k j + 1 2 (V k j -V k j+1 )F k,|•| j+1/2 -(V k j-1 -V k j )F k,|•| j-1/2 V k j . (4.38)
The first term of the right hand side of (4.38) is equal to

1 2 (V k j + V k j+1 )F k j+1/2 -(V k j-1 + V k j )F k j-1/2 V k j = 1 2 F k j+1/2 -F k j-1/2 V k j 2 + 1 2 V k j V k j+1 F k j+1/2 -V k j-1 V k j F k j-1/2 .
A similar reasoning allows us to rewrite the second term of the right hand side of (4.38) as follows

1 2 (V k j -V k j+1 )F k,|•| j+1/2 -(V k j-1 -V k j )F k,|•| j-1/2 V k j = 1 2 (V k j -V k j+1 ) 2 F k,|•| j+1/2 + 1 2 (V k j -V k j+1 )F k,|•| j+1/2 V k j+1 -(V k j-1 -V k j )F k,|•| j-1/2 V k j . We conclude that T 2,1 = 1 2 F k j+1/2 -F k j-1/2 V k j 2 + 1 2 (V k j -V k j+1 ) 2 F k,|•| j+1/2 + Γ j+1/2 -Γ j-1/2 .
We turn to T 2,2 . To this end, we rewrite G k j-1/2 and G k j+1/2 as follows

G k j-1/2 = V k j F k j-1/2 + (V k j-1 -V k j )F k,+ j-1/2 , G k j+1/2 = V k j F k j+1/2 + (V k j+1 -V k j )F k,- j+1/2 . We obtain the following expression of T 2,2 T 2,2 = (V k+1 j -V k j )V k j (F k j+1/2 -F k j-1/2 ) + T 3 where T 3 is defined by T 3 = T 3,1 + T 3,2 with T 3,1 = (V k+1 j -V k j )(V k j+1 -V k j )F k,- j+1/2 , and T 3,2 = -(V k+1 j -V k j )(V k j-1 -V k j )F k,+ j-1/2 . The Young inequality yields T 3,1 1 2 (V k+1 j -V k j ) 2 F k,- j+1/2 + 1 2 (V k j+1 -V k j ) 2 F k,- j+1/2 , and T 3,2 - 1 2 (V k+1 j -V k j ) 2 F k,+ j-1/2 - 1 2 (V k j -V k j+1 ) 2 F k,+ j+1/2
. By summing these two inequalities, we find

T 3 - 1 2 (V k+1 j -V k j ) 2 F k,+ j-1/2 -F k,- j+1/2 - 1 2 (V k j+1 -V k j ) 2 F k,|•| j+1/2 .
The last two terms are non negative. Finally, we obtain the following bound for T 2

T 2 (V k+1 j - 1 2 V k j )V k j (F k j+1/2 -F k j-1/2 ) - 1 2 (V k+1 j -V k j ) 2 F k,+ j-1/2 -F k,- j+1/2 + Γ j+1/2 -Γ j-1/2 .
We can now use the mass balance on edges (4.33) to find

T 2 - δx j 2δt (2V k+1 j -V k j )V k j (ρ k+1 j -ρ k j ) - 1 2 (V k+1 j -V k j ) 2 F k,+ j-1/2 -F k,- j+1/2 + Γ j+1/2 -Γ j-1/2 .
It remains to treat the term T 1 which recasts as

T 1 = δx j δt ρ k+1 j (V k+1 j -V k j )V k+1 j + δx j δt (ρ k+1 j -ρ k j )V k j V k+1 j = δx j 2δt ρ k+1 j V k+1 j 2 -V k j 2 + V k+1 j -V k j 2 + δx j δt (ρ k+1 j -ρ k j )V k j V k+1 j .
Gathering the last two estimates on T 1 and T 2 , and using (4.36) and (4.37) we obtain

-π k+1/2 j+1/2 -π k+1/2 j-1/2 V k+1 j δx j 2δt ρ k+1 j V k+1 j 2 -ρ k j V k j 2 + Γ j+1/2 -Γ j-1/2 + δx j 2δt ρ k+1 j - 1 2 F k,+ j-1/2 -F k,- j+1/2 V k+1 j -V k j 2 .
To conclude the proof, it remains to prove that under the assumption (4.35) we have

δx j 4δt ρ k+1 j 1 2 F k,+ j-1/2 -F k,- j+1/2
which holds since, by virtue of Lemma 4.6, we have

δx j δt ρ k+1 j + 2(F k,- j+1/2 -F k,+ j-1/2 ) = δx j δt ρ k j -F k,|•| j+1/2 -F k,|•| j-1/2 .

Global entropy inequality

We shall prove the decay of the discrete global entropy which is defined, at time t k (k 0), by the following formula

E k = J j=2 δx j 2 ρ k j V k j 2 + J j=1 δx j+1/2 Φ(ρ k j+1/2 ).
Roughly speaking, the decay of the discrete global entropy E k should be proved by summing the pointwise inequalities for internal and kinetic energy obtained in Proposition 4.2 and Proposition 4.7 respectively. When summing these inequalities, the flux terms vanish (thanks to boundary conditions) and we obtain the two following inequalities:

1 δt J j=2 δx j 1 2 ρ k+1 j V k+1 j 2 - 1 2 ρ k j V k j 2 + J j=2 π k+1/2 j+1/2 -π k+1/2 j-1/2 V k+1 j 0, 1 δt J j=1 δx j+1/2 Φ(ρ k+1 j+1/2 ) -Φ(ρ k j+1/2 ) + J j=1 π k+1/2 j+1/2 V k j+1 -V k j 0. (4.39)
Unfortunately, these two global inequalities are not sufficient to conclude to the global entropy decay since the two terms

J j=2 π k+1/2 j+1/2 -π k+1/2 j-1/2 V k+1 j and J j=1 π k+1/2 j+1/2 V k j+1 -V k j
do not compensate when summing the two relations (even after a discrete integration by part, as they do in the continuous case). Indeed, the first one involves the velocities (V k+1 j ) j∈{1,...,J+1} at time t k+1 whereas the second one involves the velocities (V k j ) j∈{1,...,J+1} at time t k . In other words, the difficulty relies on the fact that Proposition 4.2 and Proposition 4.7 do not lead to a discrete version of (2.8). To cope with this technical difficulty, we shall prove a global internal energy inequality that differs from (4.39). To this end, we need to introduce some technical ingredients.

Two additional auxiliary quantities. We introduce the two following auxiliary quantities

ρ k+1 j-1/2 = ρ k+1 j-1/2 + 2δt δx j-1/2 ρ k j-1/2 V k+1 j -V k j , ∀j ∈ {2, . . . , J + 1}, ρ k+1 j+1/2 = ρ k+1 j+1/2 - 2δt δx j+1/2 ρ k j+1/2 V k+1 j -V k j , ∀j ∈ {1, . . . , J}, (4.40)
which are defined so that

ρ k+1 j+1/2 + ρ k+1 j+1/2 2 = ρ k+1 j+1/2 + δt δx j+1/2 ρ k j+1/2 V k+1 j+1 -V k+1 j , j
(ii) ρ k+1 j+1/2 -ρ k+1 j-1/2 ρ k j+1/2 -ρ k j-1/2 , (4.43) (iii) ρ k+1 j+1/2 -ρ k+1 j-1/2 ρ k j+1/2 -ρ k j-1/2 + 4δt min(δx j-1/2 , δx j+1/2 ) ρ k j,Max V k+1 j -V k j . (4.44)
Proof.

We shall use the shorthand notation

d k,± j = d ± (ρ k j-1/2 , ρ k j+1/2 , V k j ). Conditions (4.27) imply 0 2δt δx j+1/2 d k,+ j 1 2
, and 0 -

2δt δx j-1/2 d k,- j 1 2 .
(i) The definition of ρ k+1 j-1/2 yields

ρ k+1 j-1/2 = ρ k j-1/2 - 2δt δx j-1/2 d k,- j ρ k j+1/2 -ρ k j-1/2 .
We are thus led to

ρ k+1 j-1/2 ρ k j-1/2 + 1 2 ρ k j+1/2 -ρ k j-1/2 3 2 ρ k j,Max .
(ii) Similarly, the definition of ρ k+1 j+1/2 and ρ k+1 j-1/2 yields

ρ k+1 j+1/2 -ρ k+1 j-1/2 = 1 - 2δt δx j+1/2 d k,+ j + 2δt δx j-1/2 d k,- j ρ k j+1/2 -ρ k j-1/2 ,
which lead to (4.43).

(iii) By definition of ρ k+1 j+1/2 and ρ k+1 j-1/2 , we have:

ρ k+1 j+1/2 -ρ k+1 j-1/2 = ρ k+1 j+1/2 -ρ k+1 j-1/2 - 2δt δx j ρk j V k+1 j -V k j ,
where ρk j stands for the following mean of ρ k j-1/2 and ρ k j+1/2 : ρk

j δx j = ρ k j-1/2 δx j-1/2 + ρ k j+1/2 δx j+1/2 , with 1 δx j = 1 δx j-1/2 + 1 δx j+1/2 .
We have ρk Non uniform bound for the density. As mentioned in Section 2, when assuming (2.4), uniform a priori estimates can be established for the continuous problem from the properties of the invariant regions, see [START_REF] Chueh | Positively invariant regions for systems of nonlinear diffusion equations[END_REF]. In particular the density remains bounded when the initial density lies in L ∞ and is non negative. It is not obvious that a numerical scheme preserves such a strong property, see [START_REF] Roux | Numerical stability for some equations of gas dynamics[END_REF][START_REF] Ding | Convergence of the fractional step Lax-Friedrichs scheme and Godunov scheme for the isentropic system of gas dynamics[END_REF] for discussions on Lax-Friedrichs and Godunov schemes or [5, Section 2.2] for general conditions. For this reason, we provide stepwise estimates on the updated solution. To this end, we can adapt the proof of Proposition 3.2, which itself relies on Lemma 3.2, to show that the discrete density at time t k+1 remains bounded from above and below at the price of a slightly strengthened CFL condition compared to (3.20) (which ensures only the positivity of the discrete density). Note that the bounds for the density at time t k+1 are defined from the bounds on the density at time t k and, hence, they are not uniform in time. In particular, the strengthened CFL conditions ensure that the discrete density is always positive but do not prevent from small (near vacuum) or large densities.

Lemma 4.9 Assume (h1)-(h4). Then, we have

0 < ρ k min 2 ρ k j+1/2 2 ρ k+1 j+1/2 2 max(ρ k j-1/2 , ρ k j+1/2 , ρ k j+3/2 ) 2ρ k Max , ∀j ∈ {1, . . . , J}
provided the following CFL-like conditions are satisfied:

2δt δx j+1/2 λ + (ρ k j+1/2 , V k j+1 ) + + λ -(ρ k j+1/2 , V k j ) - 1, ∀j ∈ {1, . . . , J}, δt δx j+1/2 λ + (ρ k j-1/2 , V k j ) + + λ -(ρ k j+3/2 , V k j+1 ) - 1, ∀j ∈ {2, . . . , J -1}, δt δx 3/2 λ -(ρ k 5/2 , V k 2 ) - 1, δt δx J+1/2 λ + (ρ k J-1/2 , V k J ) + 1.
(4.45)

A bluff extension of the internal energy. We shall work with auxiliary densities ρ k+1 j+1/2 ,

ρ k+1 j+1/2 which might be out the interval [ρ k min /2, 2ρ k Max ].
For this reason, we introduce a bluff extension of the internal energy, which is a mere quadratic polynomial function for small and large values of ρ:

Φ(ρ) =                L Φ ρ k min 2 (ρ) if ρ ρ k min 2 , Φ(ρ) if ρ k min 2 ρ 2ρ k Max , L Φ 2ρ k Max (ρ) if ρ 2ρ k Max ,
where

L Φ [ρ 0 ](ρ) = 1 2 Φ ′′ (ρ 0 )(ρ -ρ 0 ) 2 + Φ ′ (ρ 0 )(ρ -ρ 0 ) + Φ(ρ 0 )
is the second order Taylor expansion of Φ. In particular, we will use in the sequel the following properties:

Φ′ (ρ) a M + b M ρ, ∀ρ > 0, (4.46) Φ′ (ρ 1 ) -Φ′ (ρ 2 ) b M ρ 1 -ρ 2 , ∀ρ 1 , ρ 2 > 0, (4.47) with a M = max (ρ k min /2, 2ρ k Max ) |Φ ′ | + 2ρ k Max b M , b M = max (ρ k min /2, 2ρ k Max ) Φ ′′ . ( 4
Φ ′′ |V k j | + d |•| (ρ k j-1/2 , ρ k j+1/2 , V k j ) , (4.50)
where a M is defined by (4.48). Then, for all j ∈ {1, . . . , J}, the following inequality holds

δx j+1/2 δt Φ(ρ k+1 j+1/2 ) -Φ(ρ k j+1/2 ) + G k j+1 -G k j + π k+1/2 j+1/2 V k+1 j+1 -V k+1 j 1 4δt δx j ρ k+1 j V k+1 j -V k j 2 , (4.51) with G k 1 = G k J+1 = 0 and G k j = Φ(ρ k j-1/2 )V k+1 j - δx j-1/2 2δt Φ ρ k+1 j-1/2 -Φ(ρ k j-1/2 ) ,
G j since it involves ρ k j-3/2 , ρ k j-1/2 , ρ k j+1/2 , ρ k j+3/2 , V k j-1 , V k j and V k j+1 .
We state now the local entropy inequality and the decay of the global entropy which are the equivalents of the continuous inequality (2.8) and (2.5).

Corollary 4.12 Assume (h1)-(h4)

. The time step δt > 0 is assumed to satisfy the constraints (4.27), (4.35), (4.45), (4.49) and (4.50) which thus define the time τ ⋆ in Theorem 4.1. Then, the updated state verifies

E k+1 E k .
Proof. The assumptions allows to apply both Proposition 4.7 and Proposition 4.10. The result is obtained by summing all the obtained inequalities.

Note that, by contrast to the scheme studied in [START_REF] Herbin | Consistent explicit staggered schemes for compressible flows; part I: the barotropic Euler equations[END_REF], we prove here the decay of the global entropy, under suitable stability constraints: in [START_REF] Herbin | Consistent explicit staggered schemes for compressible flows; part I: the barotropic Euler equations[END_REF], the pointwise entropy production does not have a definite sign, but it is shown to be the sum of non positive contributions and vanishing (as the discretization parameters tend to 0) ones.

Proof of Proposition 4.10.

The assumption allow to apply Proposition 4.2. Hence, we can (re)start from (4.30) applied to Φ. Since Φ is convex, we have, for all j ∈ {2, . . . , J},

Φ ρ k+1 j-1/2 Φ ρ k+1 j-1/2 - 2δt δx j-1/2 ρ k j-1/2 V k+1 j -V k j Φ′ ρ k+1 j-1/2 , Φ ρ k+1 j+1/2 Φ ρ k+1 j+1/2 + 2δt δx j+1/2 ρ k j+1/2 V k+1 j -V k j Φ′ ρ k+1 j+1/2 .
Going back to (4.30), we find, for all j ∈ {2, . . . , J},

Φ(ρ k j+1/2 ) -Φ(ρ k j-1/2 ) V k j + δx j+1/2 2δt Φ ρ k+1 j+1/2 -Φ(ρ k j+1/2 ) + δx j-1/2 2δt Φ ρ k+1 j-1/2 -Φ(ρ k j-1/2 ) + T k 1,j -T k 21,j -T k 22,j 0, (4.52)
with

T k 1,j = 1 4 |V k j | + d |•| (ρ k j-1/2 , ρ k j+1/2 , V k j ) 2 min ( ρ k j,min , ρ k j,Max ) Φ′′ (ρ k j+1/2 -ρ k j-1/2 ) 2 , T k 21,j = -ρ k j+1/2 Φ′ (ρ k+1 j+1/2 ) -Φ′ (ρ k+1 j-1/2 ) V k+1 j -V k j , T k 22,j = -ρ k j+1/2 -ρ k j-1/2 Φ′ (ρ k+1 j-1/2 ) V k+1 j -V k j .
Inequality (4.52) readily leads to

G k j+1 -G k j -Φ(ρ k j+1/2 )(V k+1 j+1 -V k+1 j ) + δx j+1/2 2δt Φ ρ k+1 j+1/2 + Φ ρ k+1 j+1/2 -2 Φ(ρ k j+1/2 ) + T k 1,j -T k 21,j -T k 22,j -T k 3,j 0 for j ∈ {1, ..., J}, if we set T k 1,j = T k 2,j = T k 3,j = 0 for j = 1 and
T k 3,j = Φ(ρ k j+1/2 ) -Φ(ρ k j-1/2 ) (V k+1 j -V k j ), ∀j ∈ {2, . . . , J}.
Using the convexity of Φ again and (4.41) yields, for any j ∈ {1, ..., J},

Φ ρ k+1 j+1/2 + Φ ρ k+1 j+1/2 2 Φ ρ k+1 j+1/2 + δt δx j+1/2 ρ k j+1/2 V k+1 j+1 -V k+1 j Φ ρ k+1 j+1/2 + δt δx j+1/2 ρ k j+1/2 Φ′ ρ k+1 j+1/2 V k+1 j+1 -V k+1 j .
Finally, we can conclude (using the definition (4.22) of π k+1/2 j+1/2 ) that, for any j ∈ {1, ..., J},

G k j+1 -G k j + δx j+1/2 δt Φ(ρ k+1 j+1/2 ) -Φ(ρ k j+1/2 ) + π k+1/2 j+1/2 V k+1 j+1 -V k+1 j + T k 1,j -T k 21,j -T k 22,j -T k 3,j 0
T k 21,j b M ρ k j,Max ρ k+1 j+1/2 -ρ k+1 j-1/2 V k+1 j -V k j .
Owing to (4.44), we deduce the following inequality

T k 21,j b M ρ k j,Max ρ k j+1/2 -ρ k j-1/2 V k+1 j -V k j + 4 δt b M ρ k j,Max 2 min(δx j-1/2 , δx j+1/2 ) V k+1 j -V k j 2 .
By using the Young inequality, with condition (4.50), we arrive at

T k 21,j 12δt(b M ) 2 (ρ k j,Max ) 2 δx j ρ k j ρ k j+1/2 -ρ k j-1/2 2 + 1 48 + 1 48 1 δt δx j ρ k j V k+1 j -V k j 2 .
Estimate for T k 22,j . For T k 22,j , we use (4.46) and the definition of ρ k+1 j-1/2 to find

T k 22,j a M + b M |ρ k+1 j-1/2 | ρ k j+1/2 -ρ k j-1/2 V k+1 j -V k j + 2b M δt δx j-1/2 ρ k j-1/2 ρ k j+1/2 -ρ k j-1/2 V k+1 j -V k j 2 .
Using (4.42), we find

T k 22,j a M + 3 2 b M ρ k j,Max ρ k j+1/2 -ρ k j-1/2 V k+1 j -V k j + 3b M δt δx j-1/2 ρ k j,Max 2 V k+1 j -V k j 2 .
Applying the Young inequality with condition (4.49) then leads to

T k 22,j 3 2 δt a M + 3 2 b M ρ k j,Max 2 δx j ρ k j ρ k j+1/2 -ρ k j-1/2 2 + 1 6 + 1 48 1 δt δx j ρ k j V k+1 j -V k j 2 .
Estimate for T k 3,j . Using the Young inequality and since Φ ′ is a non decreasing function, we find

T k 3,j 12δt δx j ρ k j Φ ′ (ρ k j,Max ) 2 ρ k j+1/2 -ρ k j-1/2 2 + δx j 48δt ρ k j (V k+1 j -V k j ) 2 .

Conclusion. Finally

, gathering together all these information, we deduce an estimate from below of T k 1,j -T k 21,j -T k 22,j -T k 3,j with two terms containing either the factor

|V k+1 j -V k j | 2 , or |ρ k j+1/2 -ρ k j-1/2 | 2 .
The former precisely reads

- 1 4δt δx j ρ k j V k+1 j -V k j 2 (4.54)
and the latter brings a non negative contribution since δt satisfy (4.50). We finally remark that 1 2 ρ k j ρ k+1 j (as a consequence of Lemma 4.9) to obtain the conclusion from (4.53) and (4.54).

With a slight and obvious modification of the proof, we can strengthen Corollary 4.12 by including an estimate on J j=1 V k+1 j -V k j 2 . This refinement can be useful for the consistency analysis.

Numerical simulation

In this section, we present several numerical simulations to illustrate the behaviour of the scheme. We numerically solve Riemann problems: the initial data is made of two constant states (ρ l , V l ) and (ρ r , V r ) with a discontinuity located at x = 0. For such initial data, the structure of the solution is well known: it is made of three constant states (the two initial states (ρ l , V l ) and (ρ r , V r ), and an additional intermediate state (ρ m , V m )); these constant states are linked by two propagating waves, each being associated with an eigenvalue of the system. Each wave can be either a rarefaction wave or a shock wave depending on the particular values of the initial left and right states. For the simulation, the computational domain [a, b] is fitted to the region of interest, depending on the considered case. In order to keep the structure described above, we use Neumann like boundary conditions: as far as the waves do not reach the boundary, the solution coincides with the solution of the problem set on the whole line. The numerical parameters δt and δx are defined consistently with the stability assumption (3.20). For each test case, we precise the value of the following quantity:

cfl = δt δx 1 min λ ± (ρ l , V l ), λ ± (ρ m , V m ), λ ± (ρ r , V r ) .
We should bear in mind that the entropy-stability analysis of the scheme requires further restrictions which can be, in some circumstances, significantly more constrained than (3.20). This is verified in the numerical experiments.

Polytropic ideal gases

We first present simulations using the state law of polytropic ideal gases:

p(ρ) = kρ γ , ( 5.55)
where k > 0 and γ > 1 are two real constants. We begin with three test cases coming from [START_REF] Chalons | Relaxation approximation of the Euler equations[END_REF]. The state law is given by (5.55) with k = (γ-1) 2 4γ and γ = 1.6. The computational domain [a, b] and the initial data are defined as follows Test 1 (shock-shock) a=-0.2 b=0.8

ρ l = 1 ρ r = 2 V l = 1 V r = 0.5 Test 2 (rarefaction-rarefaction) a=-0.7 b=0.3 ρ l = 0.5 ρ r = 1 V l = -0.5 V r = -0.2 Test 3 (rarefaction-shock) a=-0.7 b=0.3 ρ l = 1 ρ r = 0.5 V l = -0.5 V r = -0.5
The corresponding Riemann solutions develop two shocks, two rarefaction waves and a rarefaction wave followed by a shock wave, respectively. For the simulation, we make the number J of cells within the grid vary; the time step is fixed according to the relation Jδt = 0.25 (that is cfl = 0.3 for Test 1, cfl = 0.2 for Test 2 and 3). For each test case, we plot the approximate density and velocity obtained for J = 100, 400, 3200 compared to the exact solution at time T = 0.5. We also plot the evolution of the discrete L 1 norm of the error e between the approximate solution and the exact solution (ρ ex , V ex ) at the final time T = 0.5 = nδt: The next examples are more difficult. The state law is given by (5.55) with k = 1 and γ = 1.4. Test 4 is inspired from [START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF]: the Riemann solution present two strong shocks. Test 5 is inspired from [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF]: the Riemann solution is made of two symmetric rarefaction waves and the difficulty relies on the formation of near-vacuum in the intermediate region. The computational domain and the initial data for these test cases are given by Test 4 (shock-shock) a=-0.1 b=0.15

J j=1 δx ρ n j+1/2 -ρ ex (T, x j ) and J+1 j=1 δx V n j -V ex (T, x j ) , ( 5
ρ l = 10 ρ r = 20 V l = 50 V r = 0 Test 5 (rarefaction-rarefaction) a=-0.5 b=0.5 ρ l = 1 ρ r = 1 V l = -5 V r = 5
The results are presented in Figures 5 and6 respectively. For Test 4, as previously, we plot the approximate density and velocity obtained for J = 400, 800, 3200 compared to the exact solution at T = 0.005. For these runs, we impose δtJ = 0.0004 (that is cfl = 0.08). We also plot the evolution of the discrete L 1 norm of the error e between the approximate solution and the exact solution at the final time T = 0.005 as a function of the mesh size and we provide the associated convergence rates. We obtain a convergence rate close to 1. The exact solution is well approximated, again. Nevertheless, near the first shock for the velocity, we observe a small overshoot the amplitude of which decreases with the time step. For Test 5, we plot the approximate density and velocity but also the momentum obtained for J = 200, 400, 3200 compared to the exact solution at T = 0.07 (here we impose δtJ = 0.01, that is cfl = 0.06). The velocity is poorly approximated in the near vacuum region but the evaluation of the momentum, which is the quantity of interest, is fair. As previously, we also plot the evolution of the discrete L 1 norm of the error e between the approximate solution and the exact solution at T = nδt = 0.07 as a function of the mesh size and we provide the associated convergence rates. The L 1 error norm for the momentum is defined as follows J j=1 δx ρ n j V n j -ρ ex (T, x j )V ex (t 0 , x j ) .

(5.57)

We observe a convergence rate close to 1 for the density and the momentum.

Other examples

We complete the numerical illustration with examples that departs from the standard polytropic ideal gas law. We start with the following state law:

p(ρ) = k ρ ρ * -ρ γ with k = (γ -1) 2 4γ
, γ = 0.6, ρ * = 3.

(5.58)

This constitutive law appears as a particular case of the Van der Waals state law; it is used in the modeling of dusty gases, see [START_REF] Benzoni-Gavage | Multidimensional Hyperbolic Partial Differential Equations. First-order Systems and Applications[END_REF][START_REF] Jena | Self-similar shocks in a dusty gas[END_REF] and the references therein. The interaction forces between gas molecules are ignored here but the constant ρ * is intended to introduce a correction accounting for the finite size of the molecules. In particular, we note that ρ → p(ρ) is not an homogeneous function. (In particular it is not clear how to define a co-localized kinetic scheme that makes the physical entropy decay.) Of course it is important for this problem to preserve the natural bound ρ < ρ ⋆ . Such a discrete maximum principle can be incorporated in the stability condition, as discussed in [START_REF] Berthelin | Multifluid flows: a kinetic approach[END_REF] for close-packing models in fluid-particles flows. The computational domain [a, b] and the initial data used in our simulation are defined as follows Test 6 (shock-shock) a=-0.2 b=0.8

ρ l = 1 ρ r = 2 V l = 1 V r = 0.5
The corresponding Riemann solutions develop two shocks. We perform simulations for several number J of cells in the grid; the time step being imposed by the relation Jδt = 0.25 (that is cfl = 0.3). We plot the approximate density and velocity obtained for J = 100, 200, 1600 compared to the exact solution at T = 0.5. The results are given in the Figure 7. The exact solution is well approximated: the intermediate constant state (ρ m , V m ) and the propagation speed of shocks are correctly computed. We also plot the evolution of the discrete L 1 norm of the error e between the approximate solution and the exact solution at T = 0.5 as a function of the mesh size. The convergence rate is close to 1.

Finally, we investigate the performances of the scheme with the following complicated state law

p(ρ) = -C v 0 T 0 Γ 0 ρ 0 + K 0 2 (χ + 1) 2 χ 2f 0 (χ) + χf ′ 0 (χ) + exp Γ 0 1 - ρ 0 ρ , ( 5.59)
where χ = ρ ρ 0 -1, and the function f 0 is defined by:

f 0 (χ) = 1 + s 3 -2 χ + qχ 2 + rχ 3 1 -sχ , with s = 1.5, q = - 42080895 14941154 , r = 727668333 149411540 .
The parameters K 0 , ρ 0 , T 0 , C v 0 and Γ 0 are the following constants: Proof. Let us introduce the following shorthand notation

K 0 =
p 1 = PΦ ρ 1 , ρ 1 and p 2 = PΦ ρ 2 , ρ 2 .
We observe that

ρ 1 = ρ 1 - 1 µ d -(ρ 1 , ρ 2 , V )(ρ 2 -ρ 1 ), ρ 2 = ρ 2 - 1 λ d + (ρ 1 , ρ 2 , V )(ρ 2 -ρ 1 ).
Owing to Lemma 3.2, since V 0, we know that

d + (ρ 1 , ρ 2 , V ) > 0, d -(ρ 1 , ρ 2 , V ) 0 and d |•| = d + -d -> 0.
We assume that ρ 1 < ρ 2 ; the case ρ 2 < ρ 1 can be treated in a similar way.

By (A.61), λ, µ > d |•| (ρ 1 , ρ 2 , V
), and we obtain

ρ 1 ρ 1 + d -(ρ 1 , ρ 2 , V )(ρ 1 -ρ 2 ) d |•| (ρ 1 , ρ 2 , V ) = ρ 2 + d + (ρ 1 , ρ 2 , V )(ρ 1 -ρ 2 ) d |•| (ρ 1 , ρ 2 , V ) < ρ 2 .
Thus, we have

ρ 1 ρ 1 < ρ 2 < ρ 2 .
Since Φ is strictly convex, we deduce from the three chord lemma that p 1 < p 2 .

We can now define ρ 1/2 as the solution of the following linear equation

Φ(ρ 2 ) + p 2 (ρ 1/2 -ρ 2 ) + 1 4 min ( ρ, ρ ) Φ′′ (ρ 1 -ρ 2 ) 2 = Φ(ρ 1 ) + p 1 (ρ 1/2 -ρ 1 ).
We are going to prove that ρ 1 < ρ 1/2 < ρ 2 . We start with the following equality Φ′′ .

(p 2 -p 1 )(ρ 1/2 -ρ 2 ) = Φ(ρ 1 ) -Φ(ρ 2 ) + p 1 (ρ 2 -ρ 1 ) - 1 4 min ( ρ, ρ ) Φ′′ (ρ 1 -ρ 2 ) 2 . (A.
However, owing to (A.61), we have: min

( ρ, ρ ) Φ′′ - d + (ρ 1 , ρ 2 , V ) λ max ( ρ, ρ ) Φ′′ > 1 2 min ( ρ, ρ )
Φ′′ .

Going back to (A.63), it proves that (p 1 -p 2 )(ρ 1/2 -ρ 1 ) < 0, and finally that ρ 1/2 > ρ 1 .

Proof Lemma 4.4. The result is trivial when ρ 1 = ρ 2 . Hence, we assume ρ 1 = ρ 2 . We first focus on the case V 0. We adopt the shorthand notations d ± and d |•| instead of d ± (ρ 1 , ρ 2 , V ) and d |•| (ρ 1 , ρ 2 , V ), respectively. We bear in mind the equality V = d + + d -.

B Comments on numerical diffusion and schemes comparison

It is worth discussing the expression of the numerical diffusion produced by our scheme. For the sake of simplicity, we drop the superscript k. The mass flux is given by We remind that ρV = F + (ρ, V ) + F -(ρ, V ). It allows to rewrite

F j = F + (ρ j-1/2 , V j ) + F -(ρ j+1/2 , V j ), F ± (ρ, V ) =
F j = ρ j+1/2 + ρ j-1/2 2 V j - 1 2 F + (ρ j+1/2 , V j ) -F -(ρ j+1/2 , V j ) +F -(ρ j-1/2 , V j ) -F + (ρ j-1/2 , V j ) = ρ j+1/2 + ρ j-1/2 2 V j - 1 2 F |•| (ρ j+1/2 , V j ) -F |•| (ρ j-1/2 , V j ) .
It makes the numerical mass diffusion appear, since where we remind that

F j+1 -F j = ρ j+3/2 + ρ j+1/2 2 V j+1 - ρ j+1/2 + ρ j-1/2 2 V j -∆ mass
d |•| (ρ, ρ ′ , V ) = F |•| (ρ,V )-F |•| (ρ ′ ,V ) (ρ-ρ ′ )
for ρ = ρ ′ , see (4.24). Next, let us evaluate the numerical diffusion for the momentum equation. We remind that the pressure gradient is already approximated by a centered difference. Then, we study

G j+1/2 = V j F + j + F + j+1 2 + V j+1 F - j + F - j+1 2 . (B.67)
which recasts as

G j+1/2 = F j + F j+1 2 V j + V j+1 2 + V j 4 (F + j -F - j + F + j+1 -F - j+1 ) + V j+1 4 (F - j -F + j + F - j+1 -F + j+1 ) = F j + F j+1 2 V j + V j+1 2 + 1 4 (F |•| j + F |•| j+1 )(V j -V j+1 ).
42 Thus, we identify the numerical diffusion for the momentum as

∆ q j+1/2 = 1 4 (F |•| j + F |•| j+1 )(V j+1 -V j ) - 1 4 (F |•| j-1 + F |•| j )(V j -V j-1 ) = F |•| j + F |•| j+1 4 V j+1 - 2F |•| j + F |•| j+1 + F |•| j-1 4 V j + F |•| j-1 + F |•| j 4 V j-1 .
Indeed, we have

G j+1/2 -G j-1/2 = F j + F j+1 2 V j + V j+1 2 - F j-1 + F j 2 V j-1 + V j 2 -∆ q j+1/2 .
As noticed above, both coefficients F |•| (ρ, V ) and d |•| (ρ, ρ ′ , V ) are always strictly positive for any ρ, ρ ′ > 0, V ∈ R with our choice of Maxwellian state M 0 .

A tempting simplification of the scheme consists in replacing M 0 (ρ, V )(ξ) by ρδ(ξ = V ), see Remark 3.4. Then, the transport flux reduces to the UpWind flux:

F UpW j = ρ j-1/2 V j + -ρ j+1/2 V j -. (B.68)
Note that the CFL condition guaranteeing the positivity of the density has to be modified accordingly, and now it involves the material velocity V only instead of the characteristic speeds V ± c(ρ). The stability analysis (see also Remark We illustrate this discussion with two additional numerical tests (coming from [START_REF] Herbin | Explicit staggered schemes for the compressible Euler equations[END_REF]). We use the state law (5.55) with k = 1 and γ = 2, which corresponds to the case where (1.1)-(1.2) is the so-called shallow water system. The computational domain For both test, we compare the numerical results obtained for our scheme (using mass fluxes F j defined by (B.66)) and the variant with the mass fluxes F UpW j defined by (B.68). The former is denoted by 'Kinetic' and the latter by 'UpWind' in the legend of the pictures.

The results for Test 10 are plotted in Figure 13 (approximate densities on the left hand side and approximate velocities on the right hand side). We also plot the sound speed (c j = max(c(ρ j-1/2 ), c(ρ j+1/2 )) j∈{2,...,J} obtained with our scheme which allows to easily identify the subsonic regions (defined by |V j | < c j ) where the two fluxes F j and F UpW j differ. We observe that the two schemes perfectly match out of the subsonic region but are slightly different inside. According to the discussion above, the UpWind scheme seems to be less

Remark 4 . 5

 45 In the case of the variant of the scheme presented in Remark 3.4, we haved |•| (ρ 1 , ρ 2 , V ) = |V |and inequality (4.29) reduces to the following, very simple, convexity inequality (with α

  j-1/2 , δx j+1/2 ) , and consequently, (4.44).

  10 11 , ρ 0 = 10 4 , T 0 = 300, C v 0 = 10 3 , Γ 0 = 1.5.

( a )

 a Density at T = 0.5 (b) Velocity at T = 0.5 (c) L 1 -error norm at T = 0.5 as a function of δx (with δt = 0.25δx).

Figure 2 :

 2 Figure 2: Results for Test 1 p(ρ) = k (ρ) γ .

Figure 3 :

 3 Figure 3: Results for Test 2 p(ρ) = k (ρ) γ .

Figure 4 :

 4 Figure 4: Results for Test 3 p(ρ) = k (ρ) γ .

Figure 5 :

 5 Figure 5: Results for Test 4 p(ρ) = k (ρ) γ .

Figure 6 :

 6 Figure 6: Result for Test 5 p(ρ) = k (ρ) γ .

Figure 7 :

 7 Figure 7: Results Test 6 p(ρ) = k (ρ) γ (ρ * -ρ) -γ .

Figure 8 :

 8 Figure 8: Results for Test 7 (p defined by (5.59), δt = 10 -8 , J = 1600).

( a )Figure 9 :Figure 10 :Figure 11 :Figure 12 :

 a9101112 Figure 9: Results for Test 8 (p defined by (5.59), δt = 10 -8 , J = 1600). Density at T = 10 -5 .

2 - 1 2 d

 212 | (ρ j+3/2 , ρ j+1/2 , V j+1 )(ρ j+3/2 -ρ j+1/2 ) -1 2 d |•| (ρ j+1/2 , ρ j-1/2 , V j )(ρ j+1/2 -ρ j-1/2 ) = 1 2 d |•| (ρ j+3/2 , ρ j+1/2 , V j+1 ) ρ j+3/2 + 1 2 d |•| (ρ j+1/2 , ρ j-1/2 , V j ) ρ j-1/|•| (ρ j+3/2 , ρ j+1/2 , V j+1 ) + d |•| (ρ j+1/2 , ρ j-1/2 , V j ) ρ j+1/2

4 . 5 )

 45 and the computation of the numerical diffusion proceeds exactly the same way, with nowF UpW,|•| (ρ, V ) = ρ|V | and d UpW,|•| (ρ, ρ ′ , V ) = |V |.Therefore, the numerical mass diffusion and momentum diffusion vanish when the material velocity vanishes. Note also that the right hand side in the CFL condition (4.50) vanishes when V = 0 for the flux (B.68). It certainly explains the difficulties observed in the simulations.

  [a, b] and the initial data are defined as follows Test 10 (shock-rarefaction) a=-0.5 b=0.5 ρ l = 1 ρ r = 10 V l = 3 V r = 5.5 Test 11 (shock-rarefaction) a=-0.5 b=0.5 ρ l = 1 ρ r = 10 V l = 1.89 V r = 4.39

  (3.17)-(3.18) with the momentum fluxes modified as in (4.21)-(4.22).

h4) The updated state (ρ k+1 j+1/2 ) j∈{1,...,J} and (V k+1 j ) j∈{2,...,J} is defined by (3.13)-(3.15) and

  ∀j ∈ {2, . . . , J}.

Remark 4.11 Both Propositions 4.2 and 4.10 state a discrete version of (2.6), with the identification of discrete energy fluxes. However, in Proposition 4.10, the discrete analog of p(ρ)∂ x V involves V k+1 whereas V k appears in Proposition 4.2. The price to be paid is a more stringent CFL-like condition (see (4.49) and (4.50)) and a non negative term that appears in the right hand side of the inequality. It can be compensated by the dissipation term in the kinetic energy balance stated in Proposition 4.7. Note also that the stencil of the flux G j is more extended than for the flux

  .56) as a function of the mesh size. It provides an evaluation of the convergence rate. The results are given in Figures2, 3 and 4, respectively. The exact solution is well approximated: the intermediate constant state (ρ m , V m ) and the propagation speed of the waves are correctly computed. For each test case, we obtain a convergence rate close to 1.

  62)Since Φ is strictly convex and ρ 1 < ρ 1 < ρ 2 , we arrive at(p 2 -p 1 )(ρ 1/2 -ρ 2 ) Φ(ρ 1 ) -Φ(ρ 2 ) + p 1 (ρ 2 -ρ 1 ) (ρ 2 -ρ 1 ) -Φ(ρ 1 ) -Φ(ρ 2 )This proves ρ 1/2 < ρ 2 . Similarly, we have(p 1 -p 2 )(ρ 1/2 -ρ 1 ) = Φ(ρ 2 ) -Φ(ρ 1 ) + p 2 (ρ 1 -ρ 2 ) + Φ′′ (ρ 1 -ρ 2 ) 2 . (A.63)Using the Taylor expansion of the function Φ, we can prove that the right hand side is non positive. Indeed, let h = λ -1 d + (ρ 1 , ρ 2 , V )(ρ 2 -ρ 1 ) > 0. We havep 2 = Φ(ρ 2 -h) -Φ(ρ 2 )

							ρ 1 -ρ 2	+	Φ(ρ 1 ) -Φ(ρ 1 ) ρ 1 -ρ 1	< 0.
							1 4	min ( ρ, ρ )
		-h		= Φ′ (ρ 2 ) -	1 h	ρ 2 ρ 2 -h
		ρ 2 ρ 1	Φ′′ (u)(u -ρ 1 ) du +	ρ 2 -ρ 1 h	ρ 2 ρ 2 -h	Φ′′ (u) u -(ρ 2 -h) du
	-	(ρ 1 -ρ 2 ) 2 2	min ( ρ, ρ )	Φ′′ -	h ρ 2 -ρ 1	max

Φ′′ (u)(u -(ρ 2 -h)) du, while Φ(ρ 2 ) -Φ(ρ 1 ) = Φ′ (ρ 2 )(ρ 2 -ρ 1 ) -ρ 2 ρ 1 Φ′′ (u)(u -ρ 1 ) du. It yields Φ(ρ 2 ) -Φ(ρ 1 ) + p 2 (ρ 1 -ρ 2 ) = -( ρ, ρ ) Φ′′ .

Bearing in mind the definition of h, we obtain

Φ(ρ 2 ) -Φ(ρ 1 ) + p 2 (ρ 1 -ρ 2 ) -(ρ 1 -ρ 2 ) 2 2 min ( ρ, ρ ) Φ′′ -d + (ρ 1 , ρ 2 , V ) λ max ( ρ, ρ )

(a) Density at T = 10 -5 (b) Velocity at T = 10 -5

Acknowledgements

We thank the referees for their interesting comments, that include most of the material in Appendix B.

This example is an isentropic version of a model introduced in [START_REF] Houzé | Dissipative issue of high-order shock capturing schemes with non-convex equations of state[END_REF], referred to with the nickname "Bizarrium". The equation of state is non-convex. This toy-model has been proposed to serve as a benchmark that reproduces the main features of "real-life" applications, in order to evaluate how numerical schemes select the solution when the convexity of the state law might vary. We refer the reader to [START_REF] Houzé | Dissipative issue of high-order shock capturing schemes with non-convex equations of state[END_REF] for detailed motivations and comparisons of several numerical methods, and to [START_REF] Benzoni-Gavage | Multidimensional Hyperbolic Partial Differential Equations. First-order Systems and Applications[END_REF] for further details on the mathematical theory for such general equations of state.

For the simulation, the computational domain [a, b] and the initial data are defined as follows:

Test 7 a=-0.2 b=0.8 ρ l = 11000 ρ r = 10000

Test 8 a=-0.2 b=0.8 ρ l = 14285 ρ r = 10000

We plot the approximated density and velocity obtained for J = 1600, δt = 10 -8 at time T = 10 -5 . The results for Test 7 are given in Figure 8. For this test, the density remains bounded between ρ r = 10000 and ρ l = 11000; on this range the pressure ρ → p(ρ) is a convex function so that the assumption (2.4) is satisfied. The scheme behaves very well in this case. However, this is not the case for Test 8 and 9. In the range [10000, 14285] (Test 8), the pressure p(ρ) has two convexity changes and in the range [12000, 13000] (Test 9) the pressure p(ρ) is concave. In these cases where assumption (2.4) is not satisfied, the structure of the solutions of the Riemann problems is more complex than the structure described above (see for instance [START_REF] Wendroff | The Riemann problem for materials with nonconvex equations of state. I. Isentropic flow[END_REF] and [START_REF] Houzé | Dissipative issue of high-order shock capturing schemes with non-convex equations of state[END_REF]) and the numerical results exhibit oscillations. These results illustrate the role of (2.4). Nevertheless, the scheme can be improved to treat such cases with loss of convexity. Roughly speaking the idea consists in extending the support of the Maxwellian M 0 . Designing and analyzing a scheme for a general state law is beyond the scope of this work but to illustrate the capability of the scheme we perform simulations with a simple adaptation of M 0 directly inspired from [START_REF] Kurganov | Adaptive semidiscrete central-upwind schemes for nonconvex hyperbolic conservation laws[END_REF]. In order to replace M 0 , we define the following function

where

and H is the unique cubic polynomial function that satisfies the following interpolation conditions:

For internal edges (for j ∈ {2, . . . , J}), the discrete mass fluxes are then defined as follows:

Recall that ρ k j is defined by (3.16). The momentum fluxes are then deduced from these mass fluxes as explained in Section 3.2. Figures 9 and 10 present the results for Test 8. The modification of the support of the Maxwellian allows to reduce the amplitude of the oscillations near the discontinuity: it becomes of the order of 2% (resp. 8%) of the height discontinuity instead of 30% (resp. 20%) for the density (resp. velocity). Figures 11 and12 present the results for Test 9. Simulations with the Maxwellian M 0 show very strong oscillations whereas there is only one oscillation located near each discontinuity when the Maxwellian M 0 is used. The amplitude of this oscillation is less than 25% (resp. 15%) of the height of the discontinuity for the density (resp. velocity).

A Proof of Lemma 4.4

We begin with some notation which are useful in the sequel. Let a, b ∈ R and let Φ : R → R be a function of class C 1 . We denote by P Φ (a, b) the following fraction:

We remind that d ± and d

We denote ρ = min(ρ 1 , ρ 2 ) and ρ = max(ρ 1 , ρ 2 ) and we set

Then, there exists ρ 1/2 ∈ ( ρ, ρ ) such that:

holds.

The equality we wish to prove is equivalent to

where

With the definition of ρ 2 , ρ 1 , we find that this is again equivalent to

However, both λ and µ satisfy the condition (A.61). Thus, we can apply Lemma A.1: there exist

We multiply this equality by d + + d -and substract it to (A.64). It turns out that the inequality we have to prove is equivalent to

We can rearrange terms so that we are led to prove

Since ρ 1/2 ∈ ( ρ, ρ ) and sgn(p 2 -p 1 ) = sgn(ρ 2 -ρ 1 ) (see the proof of Lemma A.1), we know that

And since d -0, it is sufficient to prove

However, equality (A.62) shows that the left hand side is exactly equal to

which is non positive since Φ is a convex function.

The result for V 0 is obtained by applying the obtained equality with -V , inverting the role of ρ 1 , λ and ρ 2 , µ, and using Lemma 3.2-(iii). diffusive but the discrepancies between the two schemes are small whereas the mesh size is quite large (δx = 10 -2 ). Results for Test 11 are plotted in Figure 13 (approximate densities on the left hand side and approximate velocities on the right hand side, Kinetic scheme at the top and UpWind scheme at the bottom). The values of initial right and left velocities V r and V l are chosen so that the velocity V m in the intermediate state nearly vanishes. We observe spurious oscillations on the solution for the Upwind scheme whereas the results obtained with the Kinetic scheme are satisfactory.

The scheme proposed and analyzed in [START_REF] Herbin | Consistent explicit staggered schemes for compressible flows; part I: the barotropic Euler equations[END_REF][START_REF] Herbin | Explicit staggered schemes for the compressible Euler equations[END_REF] has many similarities with the method devised in this paper. Indeed, [START_REF] Herbin | Consistent explicit staggered schemes for compressible flows; part I: the barotropic Euler equations[END_REF][START_REF] Herbin | Explicit staggered schemes for the compressible Euler equations[END_REF] works with staggered discretization:

• the mass flux is defined by the simple UpWinding (B.68),

• in the momentum equation, the pressure is still treated with the centered approximation, as it is quite natural in the staggered framework,

• and the convection term is discretized with the following "UpWind" flux

The convection flux can be cast as follows

Hence the numerical momentum diffusion becomes ∆ HLN,q j+1/2 = 1 2

). Again, the numerical diffusion might vanish, and spurious oscillations have been reported in the simulations with this scheme. For this reason, in practice, the authors of [START_REF] Herbin | Consistent explicit staggered schemes for compressible flows; part I: the barotropic Euler equations[END_REF] introduce an artificial viscosity (proportional to ρ∆x) in the momentum equation. The differences between our scheme and the scheme in [START_REF] Herbin | Consistent explicit staggered schemes for compressible flows; part I: the barotropic Euler equations[END_REF][START_REF] Herbin | Explicit staggered schemes for the compressible Euler equations[END_REF] can be summarized as follows:

• Definition of the mass flux: the Kinetic flux (B.66) and the UpWind flux (B.68) differ only in subsonic regions, as it is clear from Figure 1. Accordingly, stability conditions for our scheme involve the full characteristic speeds, not only the material velocity. Furthermore, our scheme induces naturally a numerical mass diffusion which is always positive.

• Definition of the momentum flux: to define the UpWinding of the convection flux (B.69) at x j+1/2 , [START_REF] Herbin | Consistent explicit staggered schemes for compressible flows; part I: the barotropic Euler equations[END_REF][START_REF] Herbin | Explicit staggered schemes for the compressible Euler equations[END_REF] uses the positive and negative parts of the average of the mass fluxes at x j and x j+1 ; instead, to define our convection flux (B.67), we use the averages of the positive mass fluxes F + and, respectively, of the negative mass fluxes F -at x j and x j+1 . Note that we define the entire mass flux by F = F + + F -but F ± in general do not coincide with its positive and negative parts [F ] ± . Furthermore, the two operations, averaging and taking positive/negative parts, do not commute. Again, due to the definition of the mass fluxes, our scheme induces a non degenerate numerical diffusion in the momentum balance. Hence, we do not need to introduce an artificial viscosity and the scheme naturally performs well when the material velocity vanishes.