C Courtes

D Coulette

E Franck

L Navoret

Vectorial kinetic relaxation model with central velocity. Application to implicit relaxations schemes

We apply flux vector splitting (FVS) strategy to the implicit kinetic schemes for hyperbolic systems. It enables to increase the accuracy of the method compared to classical kinetic schemes. The method also allows to tackle multi-scale problems, such as the low Mach number limit, for which waves with different speeds are involved. We present several possible kinetic relaxation schemes based on FVS and compare them on one-dimensional test-cases. We discuss stability issues for this kind of methods.

Introduction

Hyperbolic systems often involve multi-scale phenomena, resulting from different characteristic speeds. A typical example is given by material and acoustic waves in the (compressible) Euler system. In oder to avoid overly stringent stability condition due to the fast waves, implicit methods have been developed. However, they often lead to hard-to-solve nonlinear problems. As introduced in [START_REF] Coulette | Palindromic discontinuous Galerkin method for kinetic equations with stiff relaxation[END_REF], an alternative is to consider implicit methods based on vectorial kinetic relaxation schemes, which decouple linear transport from the non-linear dynamics. In this paper, we study and compare several variants of this method, which are based on the Flux Vector Splitting strategy and show that they allow to capture better the various waves.

Kinetic BGK (Bhatnagar, Gross and Krook) relaxation methods have been introduced in [START_REF] Bhatnagar | A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems[END_REF] to solve hyperbolic systems. The unknown of the hyperbolic system is expressed as the macroscopic moment of a kinetic distribution fonction, whose components correspond to various velocities. The kinetic distribution function satisfies a transport equation at constant velocities with a BGK source term, that makes the kinetic distribution relax to a Maxwellian distribution depend only on the macroscopic quantities. Consistency relations link the Maxwellian dsitribution and the fluxes so that the kinetic equation tends to the hyperbolic system in the limit of large relaxation frequency.

Vectorial kinetic BGK relaxation methods, introduced in [START_REF] Aregba-Driollet | Discrete Kinetic Schemes for Systems of Conservation Laws[END_REF][START_REF] Natalini | A Discrete Kinetic Approximation of Entropy Solutions to Multidimensional Scalar Conservation Laws[END_REF], are a systematic way to devise relaxation methods for any hyperbolic system: each component of the hyperbolic system is represented with the same number of components in the kinetic vector, which are associated to the same set of velocities. This is the strategy we consider in this work. Note that the Jin-Xin relaxation method [START_REF] Jin | The Relaxation Schemes for Systems of Conservation Laws in Arbitrary Space Dimensions[END_REF] enters in this class of methods. Other kinetic BGK relaxation methods take advantage of the structure of the hyperbolic problem to build less generic but more compact kinetic representations are used in lattice Boltzmann community.

The kinetic BGK model has the advantage of concentrating the nonlinearity of the system in the local source term: the transport term is linear with constant advection velocities. Therefore, implicit schemes can be easily implemented using splitting techniques. The transport equations can be solved independently using implicit schemes, but with explicit cost, while the local source term is solved using for instance a θ-scheme. In [START_REF] Coulette | Palindromic discontinuous Galerkin method for kinetic equations with stiff relaxation[END_REF], the transport equations were solved with a Discontinuous Galerkin solver. We propose here to use a semi-Lagrangian method and we restrict ourselves to uniform meshes. In order to be able to use large time steps while still preserving accuracy, we consider the second order scheme propose in [START_REF] Coulette | Palindromic discontinuous Galerkin method for kinetic equations with stiff relaxation[END_REF].

The simplest vectorial kinetic BGK model is the so-called [D1Q2] N , where N is the size of the hyperbolic system and [D1Q2] stands for dimension 1 with two velocities. It is equivalent to the Jin-Xin relaxation system [START_REF] Jin | The Relaxation Schemes for Systems of Conservation Laws in Arbitrary Space Dimensions[END_REF]. The associated scheme has been already studied in [START_REF] Graille | Approximation of mono-dimensional hyperbolic systems: A lattice Boltzmann scheme as a relaxation method[END_REF] and has good stability properties. However, in return, the scheme suffers from a too large diffusion error. More precisely, the diffusion coefficient in the equivalent equation involves the difference between the kinetic speed and the maximal characteristic speed. When the last ones exhibit large variations in space or when they have very large ratio with each other (in the low Mach number Euler system for example), this [D1Q2] N scheme induces large error and is consequently difficult to use in practice. When increasing the order of the method, the results are improved but it does not enable to consider larger time steps.

In the present work, we consider [D1Q3] N schemes that involve a zero/central velocity in the set of velocities in order to better capture the slow scale of the hyperbolic problem. The design of the scheme will be based on a Flux Vector Splitting (FVS) strategy, which consists in decomposing the flux between an upwind and a downwind part. We note that the equivalence between FVS finite-volume scheme and BGK method has been pointed out in [START_REF] Bouchut | Construction of BGK models with a family of kinetic entropies for a given system of conservation laws[END_REF]. We here apply the strategy for vectorial kinetic BGK model. Several methods have already been proposed in this framework: an upwind splitting for the advection equation [START_REF] Dubois | Stable lattice Boltzmann schemes with a dual entropy approach for monodimensional nonlinear waves[END_REF], an Osher-Salomon type splitting [START_REF] Natalini | A Discrete Kinetic Approximation of Entropy Solutions to Multidimensional Scalar Conservation Laws[END_REF] or a specific splitting for the isothermal Euler system in [START_REF] Bechereau | Élaboration de méthodes Lattice Boltzmann pour les écoulements bifluides à ratio de densité arbitraire[END_REF]. We also introduced new splittings: a Lax-Wendroff type splitting, a splitting for the low Mach number Euler system and two other based on standard FVS methods (Van Leer and AUSM).

The stability of these methods is still a difficult issue and we only give here partial results. Following the works [START_REF] Bouchut | Entropy satisfying flux vector splittings and kinetic BGK models[END_REF][START_REF] Dubois | Stable lattice Boltzmann schemes with a dual entropy approach for monodimensional nonlinear waves[END_REF], we are able to determine sufficient conditions ensuring entropic stability in the case where the transport step is performed exactly. Upwind and Osher-Salomon splitting satisfy these conditions, while the Lax-Wendroff splitting does not. When stability is not guaranteed, numerical simulations may be sensitive to the numerical parameters. Finding out the correct stability region of these schemes is an open problem.

The outline of this article is as follows. In Section 2, we introduce the implicit kinetic method and the standard [D1Q2] N scheme with consistency equation and its equivalent equation and stability properties. Section 3 presents the kinetic scheme including a central velocity, based on FVS and its equivalent equation. In the following section (Section 4), we introduce examples for different one-dimensional models, based on specific FVS. We also present a method adapted to capture the low Mach number regime of the Euler system. In Section 5, the methods are compared on several test-cases using the implicit kinetic method. Finally, in the last section (Section 6), we briefly discuss stability issues for this kind of methods. An Appendix details some proofs of properties for the classical [D1Q2] N scheme.

Implicit kinetic relaxation schemes (IKR schemes)

We are interested in the approximation of one-dimensional hyperbolic systems:

∂ t U + ∂ x F(U) = 0, (1)
with the unknown U(t, x) ∈ R N and the hyperbolic flux F : R N → R N . We follow the ideas already developed in [START_REF] Aregba-Driollet | Discrete Kinetic Schemes for Systems of Conservation Laws[END_REF][START_REF] Natalini | A Discrete Kinetic Approximation of Entropy Solutions to Multidimensional Scalar Conservation Laws[END_REF] and choose to approximate the hyperbolic system (1) by a kinetic BGK model with a small set of velocities. We introduce a kinetic vector f (t, x) ∈ R M , whose components are associated to M velocities {λ 1 , λ 2 , λ 3 , ..., λ M } ∈ R M . The macroscopic unknown U, associated to f , is defined through the linear relation U = P f , where P ∈ M N,M (R). The kinetic BGK model is then given by:

∂ t f + Λ ∂ x f = 1 ε ( f eq (U) -f ), ( 2
)
where Λ is the diagonal matrix:

Λ =                  λ 1 0 λ 2 . . . 0 λ M                 
, and f eq ∈ R M denotes an equilibrium vector, which is a function of U = P f , and satisfies the following consistency conditions:

       P f eq = U, PΛ f eq = F(U). (3)
Thus, the equilibrium function is chosen such that, as ε tends to zero, this kinetic BGK model (2) tends to the hyperbolic system (1). The generic case is detailed in [START_REF] Aregba-Driollet | Discrete Kinetic Schemes for Systems of Conservation Laws[END_REF].

Thereafter, we will focus on a specific kinetic BGK representation called "vectorial" representation, which is detailed in [START_REF] Graille | Approximation of mono-dimensional hyperbolic systems: A lattice Boltzmann scheme as a relaxation method[END_REF]. This vectorial representation is strongly related to the classical relaxation [START_REF] Jin | The Relaxation Schemes for Systems of Conservation Laws in Arbitrary Space Dimensions[END_REF] and consists in performing what was explained previously to each component of U. More precisely, we pick out q velocities {λ 1 , ..., λ q } and associate them to each component of U. The set of velocities is given by {λ 1 , ..., λ q } N and we thus obtain a total of M = qN velocities. To the i-th velocity will be associated N components of the kinetic vector: we write them

f i = ( f i,k ) 1 k N ∈ R N .
We order the components of vector f as follows

f = ( f 1 , . . . , f q ) T = ( f 1,1 , . . . , f 1,N , f 2,1 , . . . , f 2,N , . . . , f q,1 , . . . , f q,N ) T ∈ R qN ,
According to the naming convention [START_REF] Graille | Approximation of mono-dimensional hyperbolic systems: A lattice Boltzmann scheme as a relaxation method[END_REF], we will refer to this kinetic BGK vectorial representation as D1Qq N .

Vectorial implicit kinetic relaxation schemes

To solve the kinetic representation (2) and capture the limit system as ε → 0, we consider an implicit scheme similar to the one introduced in [START_REF] Coulette | Palindromic Discontinuous Galerkin Method[END_REF][START_REF] Coulette | Palindromic discontinuous Galerkin method for kinetic equations with stiff relaxation[END_REF].

Let ∆t > 0 be the time step and h > 0 be the space step. We denote f n = ( f n i,k ) the kinetic vector at time t n = n∆t. The kinetic BGK representation (2) being composed of a transport part (left-hand side) and a relaxation part (right-hand side), we consider a splitting method.

Since the matrix Λ is diagonal, the transport part consists in M independent linear transport equations. To solve these M equations, we here consider the Backward Semi-Lagrangian scheme which allows to obtain a high-order method in space and exact method in time. We obtain the following transport step: Definition 1 (Transport step). The transport step T (∆t) is given by

f * i,k (x) = I h ( f n i,k )(x -λ i ∆t), ∀ 1 ≤ i ≤ q and ∀ 1 ≤ k ≤ N, (4)
where, for any g : R → R, I h (g) is a piece-wise polynomial interpolation of the values taken by g on the mesh points.

For the relaxation step, we consider the following θ-scheme, with θ ∈ [1/2, 1],

f n+1 -f n ∆t = θ f eq (U n+1 ) -f n+1 ε + (1 -θ) f eq (U n ) -f n ε .
Remark 2.1. Applying matrix P to the previous relation and using the consistency condition P f = P f eq = U shows that the macroscopic vector is unchanged during this relaxation step:

U n+1 = U n .
This enables us to rewrite the relaxation step into

f n+1 = ε -(1 -θ)∆t ε + θ∆t f n + ∆t ε + θ∆t f eq (U n ).
Eventually, we can rewrite the relaxation step on the following canonical form.

Definition 2 (Relaxation step). The relaxation step R ε (∆t, θ) is given by

f n+1 = f n + ω f eq (U n ) -f n , (5)
with ω = ∆t ε + θ∆t . Supposing that ε ∆t, we have ω ∈ ]0, 2].
In the following we consider ε ≈ 0. In this case the choice ω = 1 corresponds to the projection on f eq (U n ), whereas ω = 2 corresponds to a symmetry with respect to f eq (U n ). They are respectively obtained for θ = 1 (fully implicit relaxation) and θ = 1/2 (Crank-Nicolson scheme) when ε = 0. In numerical results of Section 5, we will refer to parameter ω rather than θ. We thus consider the following scheme.

Definition 3 (First order splitting scheme). The scheme is defined by:

• Transport step: f * i,k (x) = I h ( f n i,k )(x -λ i ∆t), ∀ 1 ≤ i ≤ q and ∀ 1 ≤ k ≤ N. (6)
• Relaxation step:

f n+1 = f * + ω ( f eq (U n ) -f * ), with ω ∈]0, 2]. (7)
A second order in time scheme can be obtained using a Strang splitting combined with a Crank-Nicolson relaxation step

φ(∆t) = T ∆t 2 • R ε ∆t, θ = 1 2 • T ∆t 2 .
Since the first transport step of one time stage can be merged with the last transport step of the previous time stage, the second order splitting scheme is equivalent to the first order Crank-Nicolson scheme with a first and last half step of transport. Splitting schemes with higher order accuracy in time have been proposed in [START_REF] Coulette | Palindromic Discontinuous Galerkin Method[END_REF][START_REF] Coulette | Palindromic discontinuous Galerkin method for kinetic equations with stiff relaxation[END_REF].

Remark 2.2. If we choose the velocities {λ 1 , • • • , λ q }, the time step ∆t and the space step h such λ i ∆t/h is an integer for any i ∈ {1, . . . , q}, then the transport is exact (we recall that the interpolation is exact if the characteristic foot is exactly a point of the mesh) and the scheme becomes a Lattice Boltzmann Method (LBM) [START_REF] Graille | Approximation of mono-dimensional hyperbolic systems: A lattice Boltzmann scheme as a relaxation method[END_REF].

Consistency conditions (3) write in that particular case

                 i∈{-,+} f i,k = U k , ∀ 1 ≤ k ≤ N, i∈{-,+} λ i f i,k = F k (U), ∀ 1 ≤ k ≤ N, .
Hence, the equilibrium necessary writes f eq = ( f eq -, f

eq

+ ) T with f eq ± (U) = 1 2 U ± F(U) 2λ . (8)
In the limit ε → 0, asymptotic developments [START_REF] Graille | Approximation of mono-dimensional hyperbolic systems: A lattice Boltzmann scheme as a relaxation method[END_REF] show that this [D1Q2] N kinetic model is consistent with the following equation

∂ t U + ∂ x F(U) = ε ∂ x λ 2 Id-| ∂F(U) | 2 ∂ x U + O(ε 2 ),
where ∂F(U) denotes the Jacobian matrix of the hyperbolic flux F(U) and for any diagonalizable matrix A with real eigenvalues, |A| stands for the absolute value of the matrix1 . Therefore, the kinetic representation [D1Q2] N is an approximation to the initial hyperbolic problem when ε is small.

One of the advantages of the previous [D1Q2] N scheme is its high stability properties. For general hyperbolic systems, we may state the following entropic stability property, whose main steps of the proof are mentioned in [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 Systems of Conservation Laws[END_REF][START_REF] Dubois | Stable lattice Boltzmann schemes with a dual entropy approach for monodimensional nonlinear waves[END_REF] and are detailed in Appendix ??. We introduce the following notation: for any matrix A, λ max (A) denotes its spectral radius (its largest eigenvalue in modulus). •

ω ∈ [1/2, 1] • λ > | λ max (∂F(U)) |,
• the transport is exact (LBM scheme) the splitting scheme (6)- [START_REF] Charles | Enhanced Convergence Estimates for Semi-Lagrangian Schemes Application to the Vlasov-Poisson Equation[END_REF] with equilibrium (8) and velocities {-λ, λ} N is entropy stable in the sense that there exists an entropy function H( f ) such that:

H( f n+1 )(x) dx H( f n )(x) dx.
When the limit hyperbolic model is linear, the following L 2 -stability property holds.

Proposition 2.4 (L 2 stability). The splitting scheme (6)- [START_REF] Charles | Enhanced Convergence Estimates for Semi-Lagrangian Schemes Application to the Vlasov-Poisson Equation[END_REF] with equilibrium (8) and velocities {-λ, λ} N for linear model F(U) = AU is L 2 -stable on the following conditions:

• ω ∈ [0, 2], • λ > | λ max (A) |.
The guidelines of the proof can be found in [START_REF] Graille | Approximation of mono-dimensional hyperbolic systems: A lattice Boltzmann scheme as a relaxation method[END_REF].

Proposition 2.5 (Equivalent equation). The splitting scheme (6)- [START_REF] Charles | Enhanced Convergence Estimates for Semi-Lagrangian Schemes Application to the Vlasov-Poisson Equation[END_REF] with equilibrium (8) and velocities {-λ, λ} N is consistent in time with the following PDE

∂ t U + ∂ x F(U) = 1 ω - 1 2 ∆t ∂ x λ 2 Id-| ∂F(U) | 2 ∂ x U + O(∆t 2 ). (9)
We don't detail the proof since the [D1Q2] N case is a particular case to the following method.

Considering the mentioned properties, this [D1Q2] N scheme has many advantages. If an exact transport (Lattice Boltzmann Method) is considered, the scheme has first or second order accuracy in time and space and is proved to be entropy stable (Prop.2.3). However, a CFL condition, ∆t = h/λ should be satisfied (as mentioned in Remark 2.2). If a semi-Lagrangian method is used for the transport step, the scheme has high order accuracy in space and first or second order accuracy in time. Moreover, since the semi-Lagrangian is unconditionally stable [START_REF] Charles | Enhanced Convergence Estimates for Semi-Lagrangian Schemes Application to the Vlasov-Poisson Equation[END_REF], the method inherits the same property: the stability is proved for linear fluxes and is numerically observed for the nonlinear case. Consequently, with this method, nonlinear hyperbolic equations can be solved without CFL condition and without nonlinear inversion of large matrices, whereas the classical implicit method are CPU time and memory demanding.

However, as explained in the introduction, the main drawback comes from the large diffusion error appearing in the right-hand side of equation [START_REF] Coulette | Palindromic discontinuous Galerkin method for kinetic equations with stiff relaxation[END_REF]. In the [D1Q2] N scheme, all the waves are diffused using an upper bound of the fastest wave. This approximation gives rise to very large errors in two cases : in presence of large spatial gradient of the wave speeds or large ratio of the characteristic speeds. In both situations, the slow waves involved in the dynamics are not correctly captured by the scheme. We here propose a smarter choice of the velocities and equilibrium which suits best multi-scale problems.

     λ - . . . 0 λ - λ 0 . . . λ 0 λ + 0 . . . λ +                                                  (10)
and the equilibrium becomes f eq = ( f eq -, f eq 0 , f eq + ) T , defined hereafter. Since an additional degree of freedom with λ 0 and f eq 0 has been introduced, consistency conditions (3) do not determine uniquely the equilibrium f eq . As explained in the following subsection, we will based the definition of this equilibrium on a flux vector splitting (FVS). Remark 3.1. As will be discussed later, in practice, we often choose min(|λ + |, |λ -|) > λ max (∂F).

Decentered flux vector splitting

A flux splitting is a decomposition of the flux into two parts, abusively called "positive" and "negative" parts, as follows

F(U) = F + (U) + F -(U). ( 11
)
We could define the equilibrium f eq directly with this splitting method but a more judicious choice is to take into account the central velocity λ 0 in this splitting method in order to better approach the physical steady (or quasi steady) state of the hyperbolic system (1). The key point is then to apply this splitting method on a translation of F(U):

F(U) -λ 0 U = (F(U) -λ 0 U) + + (F(U) -λ 0 U) -.
For simplicity, we denote F + 0 (U) = (F(U) -λ 0 U) + and F - 0 = (F(U) -λ 0 U) -, which gives the following flux vector splitting

F(U) = F + 0 (U) + F - 0 (U) + λ 0 U. (12)
By derivating the previous relation, we obtain the same splitting on the Jacobian matrices

∂F(U) = ∂F + 0 (U) + ∂F - 0 (U) + λ 0 Id.
Let us introduce the notion of an entropic flux vector splitting for an entropy function η: Definition 4. A flux vector splitting [START_REF] Dubois | Simulation of strong nonlinear waves with vectorial lattice Boltzmann schemes[END_REF] is entropic for the entropy η, if ∂F + 0 (U) is a positivedefinite matrix (respectively ∂F - 0 (U) is a negative-definite matrix) and if there exists

ζ ± (U) de- fined by ∂ U ζ ± (U) = ∂F ± 0 (U)∂ U η(U) such that ∂ t η(U) + ∂ x ζ ± (U) ≤ 0.
Remark 3.2 (Linear flux). For a linear hyperbolic flux F(U) = AU, a possible method to split F(U) into two parts is to decompose the spectrum of A (denoted by σ(A)) between the non negative eigenvalues and the negative ones. Hence,

F + (U) = A + U and F -(U) = A -U with σ(A + ) = σ(A)χ σ(A)≥0 and σ(A -) = σ(A)χ σ(A)<0
, where χ is the indicator function. The flux vector splitting [START_REF] Dubois | Simulation of strong nonlinear waves with vectorial lattice Boltzmann schemes[END_REF] becomes in this particular case

F(U) = A + 0 U + A - 0 U + λ 0 U with σ(A + 0 ) = (σ(A) -λ 0 )χ σ(A)≥λ 0 and σ(A - 0 ) = (σ(A) -λ 0 )χ σ(A)<λ 0 .
3.2 Determination of the equilibrium f eq from a flux vector splitting

The flux vector splitting [START_REF] Dubois | Simulation of strong nonlinear waves with vectorial lattice Boltzmann schemes[END_REF] enables us to construct the equilibrium f eq for the [D1Q3] N scheme, as follows. The consistency conditions (3) give for [D1Q3] N scheme

       f eq -+ f eq 0 + f eq + = U, λ -f eq -+ λ 0 f eq 0 + λ + f eq + = F(U), which is equivalent to        f eq -+ f eq 0 + f eq + = U, (λ --λ 0 ) f eq -+ (λ + -λ 0 ) f eq + = F(U) -λ 0 U.
Thanks to the flux vector splitting [START_REF] Dubois | Simulation of strong nonlinear waves with vectorial lattice Boltzmann schemes[END_REF], we obtain the following equation for the equilibrium

f eq        f eq -+ f eq 0 + f eq + = U, (λ --λ 0 ) f eq -+ (λ + -λ 0 ) f eq + = F + 0 (U) + F - 0 (U).
This sytem is still under-determined. By analogy with the kinetic flux splitting scheme [START_REF] Mandal | Kinetic flux vector splitting for Euler equations[END_REF], we propose the following decomposition

             f eq -+ f eq 0 + f eq + = U, (λ --λ 0 ) f eq -= F - 0 (U), (λ + -λ 0 ) f eq + = F + 0 (U),
which gives the following expression of the equilibrium

                           f eq -(U) = - 1 (λ 0 -λ -) F - 0 (U), f eq 0 (U) = U - F + 0 (U) (λ + -λ 0 ) - F - 0 (U) (λ 0 -λ -) , f eq + (U) = 1 (λ + -λ 0 ) F + 0 (U). ( 13
)
Remark 3.3. There are various ways to extend the previous model with more velocities. We propose here an easy extension named the vectorial IKR [D1Q(2m + 1)] N scheme. First, we consider the following velocities set {λ -m , ..., λ -1 , λ 0 , λ +1 , ...., λ +m } N . Suppose we have a splitting of the following form

F(U) = m i=1 F -,i 0 (U) + F +,i 0 (U) + λ 0 U.
Using the consistency conditions (3) and generalizing the previous decomposition, the vectorial IKR [D1Q(2m + 1)] N scheme for the system of conservation laws is defined by the equilibrium f eq defined by:

                 f eq ±i (U) = 1 (λ ±i -λ 0 ) F ±,i 0 (U), ∀i ∈ {1, . . . , m} f eq 0 (U) = U - m i=1        F +,i 0 (U) (λ +i -λ 0 ) - F -,i 0 (U) λ 0 -λ -i )        , (14)

Equivalent equation

In this subsection, we propose studying the consistency of the vectorial IKR [D1Q3] N scheme. The stability property is more delicate to establish in a general setting. An overview of the stability issue will be detailed in Section 6.

Proposition 3.4. The equivalent equation of the splitting scheme (6)-( 7) with velocities {λ -, λ 0 , λ + } N and equilibrium (13) is given by

∂ t U + ∂ x F(U) = 1 ω - 1 2 ∆t ∂ x D(U) ∂ x U + O(∆t 2 ) ( 15
)
with the following diffusion matrix

D(U) = λ + ∂F + 0 (U) + λ -∂F - 0 (U) + λ 0 ∂F(U) -| ∂F(U) | 2 . ( 16
)
Proof. We plug the exact solution in the general [D1Q3] N scheme to obtain

f * i,k (t + ∆t, x) = f i,k (t, x -λ i ∆t), ∀i ∈ {-, 0, +} and ∀k ∈ {1, ..., N}, (17)
and

f (t + ∆t, x) = f * (t + ∆t, x) + ω f eq (U * (t + ∆t, x)) -f * (t + ∆t, x) . (18)
To obtain Relation [START_REF] Mandal | Kinetic flux vector splitting for Euler equations[END_REF], we recall that the relaxation step may be written as

f n+1 = f * + ω( f eq (U n+1 ) -f * ) since U n+1 = U n in the relaxation step (see Remark 2.1)
. Thanks to a Taylor expansion in time, we expand both sides of Eq. ( 17), ∀i ∈ {-, 0, +} and ∀k ∈ {1, ..., N},

f * i,k (t + ∆t, x) = f * i,k (t, x) + ∆t ∂ t f * i,k (t, x) + ∆t 2 2 ∂ tt f * i,k (t, x) + O(∆t 3 ), ( 19
)
f i,k (t, x -λ i ∆t) = f i,k (t, x) -λ i ∆t ∂ x f i,k (t, x) + λ 2 i ∆t 2 2 ∂ xx f i,k (t, x) + O(∆t 3 ). ( 20
)
We have

U = P f , we denote V = PΛ f = λ -f -+ λ 0 f 0 + λ + f + and W = PΛ 2 f = λ 2 -f -+ λ 2 0 f 0 + λ 2 + f + ,
where Λ is the diagonal matrix defined by [START_REF] Dubois | Third order equivalent equation of lattice Boltzmann scheme[END_REF]. Using the consistency conditions (3) and relations ( 13), we have

U = P f eq (U), V eq = PΛ f eq (U) = F(U), W eq = PΛ 2 f eq (U) = λ -F - 0 (U) + λ 0 F(U) + λ + F + 0 (U).
From ( 17) and ( 20), we obtain for those macroscopic variables the following system

                 U * (t + ∆t, x) = U(t, x) -∆t ∂ x V(t, x) + ∆t 2 2 ∂ xx W(t, x) + O(∆t 3 ), V * (t + ∆t, x) = V(t, x) -∆t ∂ x W(t, x) + O(∆t 2 ), W * (t + ∆t, x) = W(t, x) + O(∆t), (21)
where the precision of the expansion in ∆t is function of the moment order. From [START_REF] Mandal | Kinetic flux vector splitting for Euler equations[END_REF], we obtain

               U(t + ∆t, x) = U * (t + ∆t, x), V(t + ∆t, x) = V * (t + ∆t, x) + ω V eq (U(t + ∆t, x)) -V * (t + ∆t, x) , W(t + ∆t, x) = W * (t + ∆t, x) + ω W eq (U(t + ∆t, x)) -W * (t + ∆t, x) . (22)
Plugging expansions ( 21) into ( 22) results in

                         U(t + ∆t, x) =U(t, x) -∆t∂ x V(t, x) + ∆t 2 2 ∂ xx W(t, x) + O(∆t 3 ), V(t + ∆t, x) =V(t, x) -∆t∂ x W(t, x) + ω V eq U(t, x) -∆t∂ x V(t, x) -V(t, x) + ∆t∂ x W(t, x) + O(∆t 2 ), W(t + ∆t, x) =W(t, x) + ω W eq (U(t, x)) -W(t, x) + O(∆t).
We then make a Taylor expansion in time of the left-hand sides. For variable U, we have:

U(t + ∆t, x) = U(t, x) + ∆t ∂ t U(t, x) + ∆t 2 2 ∂ tt U(t, x) + O(∆t 3 ),
We perform the same expansion for V and W.

We then consider a Chapman-Enskog expansion in ∆t by considering the following expansions:

U = U 0 + ∆tU 1 + ∆t 2 U 2 + O(∆t 3 ), V = V 0 + ∆tV 1 + O(∆t 2 ) and W = W 0 + O(∆t).
In the following, we gather all terms of the same order in ∆t.

Terms in O(1). We have:

               U 0 = U 0 , V 0 = V 0 + ω V eq (U 0 ) -V 0 , W 0 = W 0 + ω W eq (U 0 ) -W 0 ,
which gives V 0 = V eq (U 0 ) and W 0 = W eq (U 0 ).

Terms in O(∆t).

We have the following system:

       U 1 + ∂ t U 0 = U 1 -∂ x V 0 , V 1 + ∂ t V 0 = V 1 -∂ x W 0 + ω ∂V eq (U 0 )U 1 -∂V eq (U 0 )∂ x V 0 -V 1 + ∂ x W 0 .
Using the relations obtained previously V 0 = V eq (U 0 ), W 0 = W eq (U 0 ), the system becomes

       ∂ t U 0 + ∂ x V eq (U 0 ) = 0, ∂ t V eq (U 0 ) + ∂ x W eq (U 0 ) = ω ∂V eq (U 0 )U 1 -∂V eq (U 0 )∂ x V eq (U 0 ) -V 1 + ∂ x W eq (U 0 ) . ( 23
)
From the second equation, we get the following expression of

V 1 V 1 = ∂V eq (U 0 )U 1 -∂V eq (U 0 )∂ x V eq (U 0 ) + ∂ x W eq (U 0 ) - 1 ω (∂ t V eq (U 0 ) + ∂ x W eq (U 0 )) = ∂V eq (U 0 )U 1 -(∂V eq (U 0 )) 2 ∂ x U 0 + ∂W eq (U 0 )∂ x U 0 - 1 ω (∂V eq (U 0 )∂ t U 0 + ∂W eq (U 0 )∂ x U 0 ).
Thanks to the first equation, that gives the relation

∂ t U 0 = -∂ x V eq (U 0 ), we obtain V 1 = ∂V eq (U 0 )U 1 - 1 ω -1 ∂W eq (U 0 )∂ x U 0 -∂V eq (U 0 ) 2 ∂ x U 0 .
Terms in O(∆t 2 ). We have the following equation

U 2 + ∂ t U 1 + 1 2 ∂ tt U 0 = U 2 -∂ x V 1 + 1 2 ∂ xx W 0 .
It becomes:

∂ t U 1 + ∂ x V 1 = - 1 2 ∂ tt U 0 + 1 2 ∂ xx W eq (U 0 ). ( 24
)
Using Equation (23), we get

∂ tt U 0 = -∂ tx V eq (U 0 ) = -∂ t ∂V eq (U 0 )∂ x U 0 = -∂ 2 V eq (U 0 )∂ t U 0 ∂ x U 0 -∂V eq (U 0 )∂ x ∂ t U 0 = ∂ 2 V eq (U 0 )∂V eq (U 0 ) [∂ x U 0 ] 2 + ∂V eq (U 0 ) 2 ∂ xx U 0 + ∂V eq (U 0 )∂ 2 V eq (U 0 ) [∂ x U 0 ] 2 = ∂ x ∂V eq (U 0 ) 2 ∂ x U 0 .
By replacing V 1 in (24), we obtain

∂ t U 1 + ∂ x ∂V eq (U 0 )U 1 - 1 ω -1 ∂W eq (U 0 ) -∂V eq (U 0 ) 2 ∂ x U 0 = 1 2 ∂ x -∂V eq (U 0 ) 2 ∂ x U 0 + 1 2 ∂ x ∂W eq (U 0 )∂ x U 0 ,
and then

∂ t U 1 + ∂ x ∂V eq (U 0 )U 1 = ∂ x 1 ω - 1 2 ∂W eq (U 0 ) -∂V eq (U 0 ) 2 ∂ x U 0 . (25)
Conclusion Adding together Equations ( 23) and (25), we finally have

∂ t (U 0 + ∆tU 1 )+∂ x V eq (U 0 ) + ∆t∂V eq (U 0 )U 1 = ∆t∂ x 1 ω - 1 2 ∂W eq (U 0 ) -∂V eq (U 0 ) 2 ∂ x U 0 .
We use then V eq (U) = F(U) and

W eq (U) = λ -F - 0 (U) + λ 0 F(U) + λ + F + 0 (U).
Remark 3.5. The diffusion matrix (16) associated to the [D1Q3] N scheme is different from the one of the [D1Q2] N scheme (see [START_REF] Coulette | Palindromic discontinuous Galerkin method for kinetic equations with stiff relaxation[END_REF]). For the latter, the "positive" term (the term λ 2 Id) depends only on the speed λ (= λ + = -λ -) and is constant in space. For the [D1Q3] N scheme, the "positive" term (the term ∂W eq (U 0 )) depends on the Jacobian of the flux: it thus varies in space.

In the next section, we will provide its expression for several flux splittings and we will see that this is a way to reduce the error when the solution exhibits large variations in space.

Applications in 1D

In this section, we present different flux splittings for the [D1Q3] N scheme for various particular equations and discuss their relative performance. Some of them were already proposed in a relaxation/Lattice Boltzmann formalism but, we also investigate other flux splittings coming from classical ideas in finite volume schemes literature.

Flux splittings for the Burgers equation

We consider the Burgers equation

∂ t ρ + ∂ x ρ 2 2 = 0, (26)
and we will consider three examples of flux splittings F ± 0 which allow to define the vectorial IKR [D1Q3] 1 scheme.

Scheme 1 (D1Q3 Rusanov). For this [D1Q3] 1 scheme, the velocities are λ -≤ λ 0 = 0 ≤ λ + and the flux splitting is of Rusanov type, given by

F - 0 (ρ) = - λ - λ + -λ - 1 2 ρ 2 -λ + ρ and F + 0 (ρ) = λ + λ + -λ - 1 2 ρ 2 -λ -ρ . ( 27
)
Using Equation ( 15), the associated scheme has the following diffusion coefficient

D(ρ) = (λ + + λ -) ρ -λ -λ + -ρ 2 .
Note that if λ -= -λ + , the scheme reduces to the classical [D1Q2] 1 scheme and then the diffusion coefficient is the same as the one of the [D1Q2] 1 scheme. Considering non-symmetric velocities, λ --λ + , can enable to reduce the diffusion. For instance, if ρ is positive, λ -can be chosen equal to 0. Then the diffusion coefficient of the [D1Q3] 1 Rusanov scheme equals ρ(λ + -ρ) and is thus weaker than the one of the classical [D1Q2] 1 scheme, equal to λ 2 + -ρ 2 (since λ + > ρ(t, x)) . Scheme 2 (D1Q3 upwind). For this [D1Q3] 1 scheme, the velocities are λ -≤ λ 0 ≤ λ + and the flux splitting is of Osher (or upwind) type, given by

F - 0 (ρ) = χ ρ<λ 0 1 2 ρ 2 -λ 0 ρ and F + 0 (ρ) = χ ρ>λ 0 1 2 ρ 2 -λ 0 ρ . ( 28
)
where, for any set E, χ E denotes the indicator function.

This flux splitting has already been proposed in [START_REF] Dubois | Stable lattice Boltzmann schemes with a dual entropy approach for monodimensional nonlinear waves[END_REF]. Thanks to [START_REF] Jin | The Relaxation Schemes for Systems of Conservation Laws in Arbitrary Space Dimensions[END_REF], this scheme has the following diffusion coefficient

D(ρ) = (ρ -λ 0 )(λ --ρ)χ ρ<λ 0 + (ρ -λ 0 )(λ + -ρ)χ ρ>λ 0 .
This [D1Q3] 1 upwind enables to reduce diffusion. For instance, in the case λ 0 = 0, the diffusion coefficient equals ρ 1 Rusanov scheme, there is no need to adapt the values λ -and λ + with respect to the sign of ρ and we can thus consider solution with variable sign with lower diffusion coefficient. Moreover, if λ 0 is a good approximation of an average of ρ, the diffusion is even lowered.

(λ + -ρ) if ρ > 0 and ρ (λ --ρ) if ρ < 0. Compared to the [D1Q3]
Scheme 3 (D1Q3 Lax-Wendroff). For this [D1Q3] 1 scheme, the velocities are λ 0 = 0 and -λ -= λ + = λ > 0 and the flux splitting is of Lax-Wendroff type, given by

F - 0 (ρ) = 1 4 ρ 2 -α ρ 3 6λ and F + 0 (ρ) = 1 4 ρ 2 + α ρ 3 6λ , with α ∈ [1, 2]. ( 29
)
Using this splitting, the scheme has the following diffusion coefficient

D(ρ) = (α -1) ρ 2 .
Hence, the diffusion vanishes for α = 1: this leads to a second-order in time scheme.

Flux splittings for the advection equation with variable velocity

The advection and advection-diffusion cases are already treated in [START_REF] Dubois | Third order equivalent equation of lattice Boltzmann scheme[END_REF] and [START_REF] Dubois | Stable lattice Boltzmann schemes with a dual entropy approach for monodimensional nonlinear waves[END_REF]. Here we consider an advection equation with a variable velocity a(x) ∈ R

∂ t ρ + ∂ x a(x)ρ = 0. ( 30
)
We further suppose that a is increasing: a (x) > 0, for all x ∈ R. This condition ensures the L 2 dissipation and consequently stabilty for smooth solutions of [?]. Indeed, this model has the following entropy-entropy flux pair (η(ρ), ζ(ρ)) = ( ρ 2 2 , a(x)ρ 2 2 ):

∂ t η(ρ) + ∂ x ζ(ρ) = -a (x) ρ 2 2 ≤ 0.
Because of the space dependent flux, this model does not directly enter the framework presented so far. However, the method can be adapted to this case. The three flux splittings introduced for the Burgers equation can be transposed for this equation. Scheme 4 (D1Q3 Rusanov). Given λ -≤ λ 0 = 0 ≤ λ + , the flux splitting is set to

F - 0 (ρ, x) = -λ - (a(x)ρ -λ + ρ) λ + -λ - and F + 0 (ρ, x) = λ + (a(x)ρ -λ -ρ) λ + -λ - . (31)
Scheme 5 (D1Q3 upwind). Given λ -≤ λ 0 ≤ λ + , the flux splitting is set to

F - 0 (ρ, x) = χ a(x)<λ 0 (a(x) -λ 0 )ρ and F + 0 (ρ, x) = χ a(x)>λ 0 (a(x) -λ 0 )ρ. ( 32
)
with χ the indicator function.

Scheme 6 (D1Q3 Lax-Wendroff). Given λ -≤ λ 0 = 0 ≤ λ + , with -λ -= λ + = λ > 0, the flux splitting is set to

F - 0 (ρ, x) = 1 2 a(x)ρ -α a(x) 2 ρ 2λ , F + 0 (ρ, x) = 1 2 a(x)ρ + α a(x) 2 ρ 2λ with α ∈ [1, 2]. ( 33
)
The diffusion coefficients associated to the corresponding schemes are given by

D Rusanov (ρ, x) = (λ + + λ -) a(x) -λ -λ + -a(x) 2 , D upwind (ρ, x) = (λ --a(x)) (a(x) -λ 0 ) χ a(x)<λ 0 + (λ + -a(x)) (a(x) -λ 0 ) χ a(x)>λ 0 , D Lax-Wendroff (ρ, x) = (α -1) a(x) 2 .
As in the previous case (Burgers equation), the same remarks still hold as regards the magnitude of these diffusion coefficients.

Flux splittings for the isothermal Euler system

In this subsection, we consider the 1D isothermal Euler system

     ∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρu) + ∂ x ρu 2 + c 2 ρ = 0, (34)
where ρ(t, x) ∈ R is the density of the fluid, u(t, x) ∈ R its velocity and c > 0 the sound speed.The different flux vector splittings that we propose are listed below. They all involve a zero central velocity: λ 0 = 0.

Scheme 7 (D1Q3 Rusanov). For this [D1Q3] 2 scheme, the velocities are λ -< λ 0 = 0 < λ + and the flux splitting is of Rusanov type, given by

F - 0 (U) = - λ - λ + -λ - ρu -λ + ρ ρu 2 + c 2 ρ -λ + ρu and F + 0 (U) = λ + λ + -λ - ρu -λ -ρ ρu 2 + c 2 ρ -λ -ρu .
As mentioned previously, for symmetric velocities λ -= -λ + , we recover exactly the classical [D1Q2] 2 scheme. Let us mention additionally that this Rusanov type flux splitting allows to have an entropy stable [D1Q3] 2 scheme (see Proposition 6.1 below).

Scheme 8 (D1Q3 HLL). For this [D1Q3] 2 scheme, the velocities are λ -< λ 0 = 0 < λ + and the flux splitting is of HLL type, given by

F - 0 (U) = - (u -c) 2c ρu -(u + c)ρ ρu 2 + c 2 ρ -(u + c)ρu and F + 0 (U) = (u + c) 2c ρu -(u -c)ρ ρu 2 + c 2 ρ -(u -c)ρu .
This splitting has been obtained from the Rusanov splitting by replacing the diffusion velocities (equal to the maximal and minimal numerical velocities λ ± ) by local approximations (uc and u + c).

Scheme 9 (D1Q3 Van-Leer). For this [D1Q3] 2 scheme, the velocities are λ -< λ 0 = 0 < λ + and the flux splitting is of Van-Leer type, given by

F ± 0 (U) =                    0, if ∓ M > 1, ± 1 4 ρc(M ± 1) 2        1 ±2c        , if |M| < 1, F(U), if ± M > 1,
where M = u/c > 0 is the Mach number.

This scheme is devised such that F + 0 (U) = F(U) and F - 0 (U) = 0 (resp. F + 0 (U) = 0 and F - 0 (U) = F(U)) in the case of hypersonic flows M > 1 (resp. subsonic flows).

Scheme 10 (D1Q3 Bechereau). For this [D1Q3] 2 scheme, the velocities are λ -< λ 0 = 0 < λ + and the flux splitting is given by

F ± 0 (U) = ± 1 λ + -λ - ρu(u -λ ∓ ) + c 2 ρ ρu 2 (u -λ ∓ ) + c 2 ρ(3u -λ ∓ )
.

This flux splitting has been proposed in [START_REF] Bechereau | Élaboration de méthodes Lattice Boltzmann pour les écoulements bifluides à ratio de densité arbitraire[END_REF] to tackle test-cases with large ratio of density. The system is obtained with the consistency relations and two other constraints on the second discrete moment by mimicking the moment theory for the Euler equation.

Scheme 11 (D1Q3 Osher-Salomon). For this [D1Q3] 2 scheme, the velocities are λ -< λ 0 = 0 < λ + and the flux splitting is given by

F ± 0 (U) = F(U) ± 1 0 | ∂F(U 0 + s (U -U 0 )) | ds (U -U 0 ), U 0 = 1 |Ω| Ω U(x) dx.
In practice, the integral is computed using a quadrature method.

The splitting is based on the Osher-Salomon approximated scheme proposed by Dumsber and Toro in [START_REF] Dumbser | On Universal Osher-Type Schemes for General Nonlinear Hyperbolic Conservation Laws[END_REF]. This scheme ensures that ∂F + 0 (U) -∂F - 0 (U) is a good approximation of | ∂F(U) | and consequently, in the case where λ + = -λ -= λ, the diffusion matrix satisfies

λ + ∂F + 0 (U) + λ -∂F - 0 (U) -| ∂F(U) | 2 ≈ λ | ∂F(U) | -| ∂F(U) | 2 .
This can be viewed as a generalization of the diffusion coefficient of the upwind scheme in the scalar case (Burgers or varying velocity advection equations). This scheme is actually not specific to the Euler system: like the Rusanov scheme and contrary to the other ones of this subsection, it can be used for any hyperbolic system.

Flux splittings for the Euler system

In this subsection, we consider the 1D full Euler system

             ∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρu) + ∂ x ρu 2 + p = 0, ∂ t E + ∂ x (Eu + pu) = 0, (35)
where ρ(t, x) ∈ R is the density of the fluid, u(t, x) ∈ R its velocity, p(t, x) ∈ R the pressure. The total energy E(t, x) ∈ R is given by E = p γ-1 + ρ u 2 2 . As before, we will propose different flux splittings to construct [D1Q3] 3 schemes. We first define the D1Q3 Rusanov and D1Q3 Osher-Salomon schemes, which are generic schemes.

Scheme 12 (D1Q3 Rusanov). For this [D1Q3] 3 scheme, the velocities are λ -< λ 0 = 0 < λ + and the flux splitting is of Rusanov type, given by

F - 0 (U) = - λ - λ + -λ -           ρu -λ + ρ ρu 2 + p -λ + ρu Eu + pu -λ + E           and F + 0 (U) = λ + λ + -λ -           ρu -λ -ρ ρu 2 + p -λ -ρu Eu + pu -λ -E           .
Scheme 13 (D1Q3 Osher-Salomon). For this [D1Q3] 3 scheme, the velocities are λ -< λ 0 = 0 < λ + and the flux splitting is of Osher-Salomon type, given by

F ± 0 (U) = F(U) ± 1 0 | ∂F(U 0 + s (U -U 0 )) | ds (U -U 0 ).
U 0 is chosen as the average quantities at the time n. We use an approximation of the absolute value for the Jacobian based on a Halley approximation [6].

Now, we present two schemes which are known in the finite volume community to be efficient in the supersonic case, that is when all the waves move in the same direction. By using ideas of [START_REF] Van Leer | Flux-vector splitting for the Euler equations[END_REF], we first propose a Van-Leer type flux splitting. We also consider the Advection Upstream Splitting Method (AUSM) flux splitting introduced in [START_REF] Liou | A New Flux Splitting Scheme[END_REF] and [START_REF] Liou | A Sequel to AUSM: AUSM+[END_REF]. Scheme 14 (D1Q3 Van Leer). For this [D1Q3] 3 scheme, the velocities are λ -< λ 0 = 0 < λ + and the flux splitting is of Van-Leer type, given by

F ± 0 (U) =                                  0, if ∓ M > 1, ± 1 4 ρc(M ± 1) 2                     1 (γ -1)u ± 2c γ ((γ -1)u ± 2c) 2 2(γ + 1)(γ -1)                     , if |M| < 1, F(U), if ± M > 1,
where c = γp/ρ denotes the acoustic speed and M = u/c.

Scheme 15 (D1Q3 AUSM).

For this [D1Q3] 3 scheme, the velocities are λ -< λ 0 = 0 < λ + and the flux splitting is of AUSM type, given by

F ± 0 (U) = M ±           ρc ρuc ρHc           +           0 p ± 0          
, with c = γp/ρ the acoustic speed and H = E+p ρ the enthalpy and where M ± and p ± are given by

M ± =            M ± |M| 2 , if |M| > 1, ± 1 4 (M ± 1) 2 ± β (M 2 -1) 2 , if |M| ≤ 1, p ± =            1 ± sign(M) 2 p, if |M| > 1, 1 2 (1 ± M)p, if |M| ≤ 1.
with M = u/c and β ∈ R is a parameter, whose choice is discussed in [START_REF] Liou | A Sequel to AUSM: AUSM+[END_REF][START_REF] Liou | A New Flux Splitting Scheme[END_REF].

In the incompressible low-Mach regime, the fluid velocity has magnitude much lower than the acoustic speed, the pressure and the fluid velocity are almost constant and the density is advected at the velocity u (material wave). This asymptotic dynamics is not correctly captured by the D1Q2 model since the diffusion error is of order of the acoustic velocity instead of the velocity u. For the four previous flux splitting, this material wave condition is still not preserved.

Therefore, we devise a specific flux splitting which allows to preserve the material wave. We split the flux into an advection part and an acoustic part as follows

F(U) =           ρu ρu 2 + p Eu + pu           =           ρu (ρu)u Eu           +           0 p pu           .
The same splitting has been considered in [START_REF] Zha | Numerical solutions of Euler equations by using a new flux vector splitting scheme[END_REF][START_REF] Toro | Flux splitting schemes for the Euler equations[END_REF]. We then consider flux splittings for each part separately.

For the advection part of the flux, we can consider flux splittings proposed for the advection equation in Subsection 4.2. We here consider the Lax-Wendroff flux vector splitting, with velocities -λ, 0, λ. For the acoustic part of the flux, the splitting is devised using a reconstruction of the pressure similar to the one appearing in AUSM but where the acoustic speed c is replaced by λ + . This modification enables to preserve steady states with constant fluid velocity and pressure.

The proposed splitting thus writes Scheme 16 (D1Q3 Low-Mach). For this [D1Q3] 3 scheme, the velocities are λ 0 = 0 and -λ -= λ + = λ. The flux splitting is adapted for the low Mach number flows

F ± 0 (U) = 1 2             (ρu ± α u 2 λ ρ) ± p λ (ρu 2 ± α u 3 λ ρ) + p (1 ± γ u λ ) (Eu ± α u 2 λ E) + pu ± 1 λ γ(u 2 + λ 2 )p             .
A pressure term has been added in the first equation to stabilize the scheme. Like the D1Q3 Lax-Wendroff for transport, the dissipation depend on u and not λ, we have a minimal dissipation of the convection part for α = 1 and increasing α enables to stabilize the method. However, there is no guarantee that the dispersive terms error in the equivalent equation do not depend on λ.

To avoid this possible problem, we propose a [D1Q5] 3 vectorial scheme (see construction ( 14)) with the following velocity set -λ f , -λ s , 0, λ s , λ f

3

. The slow velocities ±λ s are used for the convection part and the fast velocities ±λ f for the acoustic one. We consider the following flux splitting for slow and fast velocities:

F ±,s 0 (U) = 1 2              (ρu ± α u 2 λ s ρ) ± p λ s ρu 2 ± α u 3 λ s ρ Eu ± α u 2 λ s E              F ±, f 0 (U) = 1 2             0 p(1 ± γ u λ f ) pu ± 1 λ f γ(u 2 + λ 2 f )p             . ( 36
)
5 Numerical results for IKR schemes

In this section, we numerically compare all the flux splittings introduced in Section 4 for the different equations. The IKR method enables to use large time steps. Consequently, the time accuracy of the different splittings will be examined in details. We also focus on the ability of the schemes to capture the involved slow waves.

The Burgers equation

Here, we compare the accuracy in time for the different flux splittings presented in Section 4.1.

Test case 1 : smooth wave. We consider the Burgers equation on the domain [0, 1] with homogeneous Dirichlet boundary conditions. The initial datum is given by

ρ(t = 0, x) = sin(2πx),
and the final time is set to T f = 0.12 since the analytical solution can be computed up to that time. The transport step is solved using a Semi-Lagrangian scheme of order 17 in space (see ( 4)) with 10000 cells. In Table 1, we compare the different splittings with velocities λ 0 = 0 and -λ -= λ + = 1.5. We note that we obtain the expected order of accuracy: the second order for the Lax-Wendroff (LW) scheme with α = 1 and the first order for the others. We also observe that the D1Q3 Rusanov method, which is equivalent to the [D1Q2] 1 scheme in that case, is the least accurate.

Rusanov

Upwind LW α = 1.0 LW α = 2.0 Error Order Error Order Error Order Error Order ∆t = 0.02

7.4E -2 - 1.9E -2 - 8.0E -3 - 1.7E -2 - ∆t = 0.01 3.9E -2 0.88 1.1E -2 0.79 2.3E -3 1.80 1.0E -2 0.77 ∆t = 0.005
2.1E -2 0.89 6.4E -3 0.78 6.0E -4 1.94 5.7E -3 0.81 ∆t = 0.0025 1.1E -2 0.93 3.5E -3 0.87 1.5E -4 2.00 3.0E -3 0.93 ∆t = 0.00125 5.4E -3 1.03 1.8E -3 0.96 3.9E -5 1.95 1.5E -3 2.00 Table 1: (Burgers equation, test 1) Convergence in time different flux splittings using the first order in time scheme (ω = 1). LW stands for Law-Wendroff.

Test case 2: rarefaction wave. Still in the domain [0, 1], the initial condition is now given by

ρ(t = 0, x) = ρ L 1 [0,0.5[ (x) + ρ R 1 [0.5,1] (x),
with 1 E is the indicator function of the set E and ρ L = 2 > ρ R = 1. The final time is set to T f = 0.5. The transport step is solved using a Semi-Lagrangian scheme of order 11 in space (see [START_REF] Bouchut | Construction of BGK models with a family of kinetic entropies for a given system of conservation laws[END_REF]) with 2000 cells. The velocity set is chosen equal -λ -= λ + = 2.5 and λ 0 = 0. We perform the numerical simulations for two different time steps: ∆t = 0.002 and ∆t = 0.01. They respectively correspond to a finite volume CFL about 10 and 50. The results for the first order in time scheme are gathered in Figure 1. Here again, we observe that the [D1Q3] 1 scheme with a Lax-Wendroff splitting with α = 1 leads to the most accurate result. However, the solution oscillates due to the numerical dispersion. As previously, the D1Q3 Rusanov splitting is the most diffusive. Comparing the Lax-Wendroff with α = 1.5 and the upwind splittings, the results are quite similar except that the maximum principle is only preserved by the second one.

The variable advection equation

Here, we compare the different flux splittings proposed in Subsection 4.2 for the advection equation (30) with constraints a(x) > 0 and a (x) > 0. Like in the previous subsection, we focus on the study of their time accuracy. For the transport step, we thus use a Semi-Lagrangian scheme with 17-th order of accuracy with 10000 cells.

Test case 1. We consider the advection equation with velocity given by a(x) = x on the domain Ω = [0, 2.5] and with final time T f = 0.3. The initial density is taken equal to

ρ(t = 0, x) = e -x 2 σ .
with σ = 0.005. We compare the different schemes using the same velocities set, λ 0 = 0 and -λ -= λ + = 3, in Tables 2 and3. As expected, we observe that the three schemes (Rusanov, upwind and Lax-Wendroff for α > 1) have convergence close to the first order accuracy for ω = 1 (Table 2) and larger than the second order accuracy for ω = 2 (Table 3). This super-convergence results for all the schemes may result from the choice of the simple velocity. We also remark that the Rusanov scheme is less accurate than the others. For the Lax-Wendroff scheme with α = 2, we obtain second order accuracy with ω = 1 since the dissipation error term cancels out exactly. When ω = 2, this splititng does not lead to a better order of convergence than the other since the dispersive term is not suppressed. Note also that the error when using large time steps is quite large since the numerical parameters may be closed to their stability limit.

Test case 2. We consider the advection velocity a(x) = 1.0 + 0.01x 2 on the domain Ω = [0, 4] and with final time T f = 1 The initial density is the same as previously. For this test case, we do not have a closed-form analytical solution. The aim is to compare the numerical diffusion/dispersion in time when using large time steps for the different flux splittings.

On Figure 2, we observe that the Lax-Wendroff flux splitting is the most accurate when using the first order in time scheme and the upwind splitting is the second most accurate. For the former, the dispersive effects are visible when using ∆t = 0.1 (right figure) which indicates that in that case the scheme actually has second order accuracy in time. When

= 0.1 9.6E -2 - 1.2E -3 - 2.5E -1 - 9.2E -3 - ∆t = 0.05 3.8E -2 1.33 1.2E -4 3.32 1.2E -0 - 2.5E -2 - ∆t = 0.025
5.3E -3 2.84 8.1E -6 3.90 4.1E -1 1.60 1.7E -3 3.87 ∆t = 0.0125 3.7E -4 3.86 5.3E -7 3.82 1.1E -4 11.8 1.2E -5 6.97 ∆t = 0.00625 2.3E -5 4.00 3.3E -8 4.00 6.2E -6 4.19 6.8E -7 4.14 Table 3: (Variable advection, test 1) Convergence in time for the second order time scheme (ω = 2) [D1Q3] 1 and different flux splittings. LW stands for Law-Wendroff. order in time scheme, in Figure 3, the Lax-Wendroff splitting is unstable and the upwind splitting seems the most accurate. Hence, the Lax-Wendroff splitting seems the best choice for the first order in time scheme whereas, for the second order in time scheme, a more judicious choice would be the upwind splitting. In Figure 4, we study the effect of the choice of the velocities λ ± . For the Rusanov and the upwind splitting, we compare the use of symmetric maximal velocities (-λ -= λ + = 1.5) with the one of unsymmetric velocities closer to the characteristic advection velocity (-λ -= 0.5 and λ + = 1.5). For the two splittings, the second set of velocities leads a significant increase of accuracy.

The isothermal Euler system

In this subsection, we present results for the splittings introduced for the isothermal Euler equations in Subsection 4.3. For the two test cases below, the domain is taken equal to Ω = [-1, 1] and is discretized with 4000 cells. The transport equation are solved using a Semi-Lagrangian scheme with 11-th order of accuracy interpolation.

Test case 1: wave propagation. We have consider the following initial data with σ = 0.005 and the initial constant fluid velocity either equals U 0 = 0 or U 0 = 0.98. The final time is set to T f = 0.5 in the first case (U 0 = 0) and to T f = 0.3 in the second one (U 0 = 0.98). We choose c = 1 and the kinetic velocities are taken equal to λ ± = ±1.2 and λ 0 = 0.

       ρ(t = 0, x) = 1.0 + 0.1e -x 2 σ , u(t = 0, x) = U 0 ,
In Figures 5 and6, we compare the different splittings for three values of the relaxation parameter ω introduced in ( 5): ω = {1, 1.5, 2}. Figure 5 corresponds to the initial condition U 0 = 0 and Figure 6 to U 0 = 0.98. As expected, the [D1Q3] 2 scheme with a Rusanov splitting (equivalent to [D1Q2] 2 scheme) leads to the most diffusive solution. For the low Mach number test case (M = u/c = 0) when U 0 = 0, the Bechereau type splitting and the Osher one are the most accurate when using the first order accuracy in time scheme (ω < 2). With the second order accuracy in time scheme (ω = 2), the Bechereau type splitting is very dispersive contrary to the Osher or the Van-Leer splittings, which provide good results. Contrary to the low Mach number case, we will see that the Van-Leer [D1Q3] 2 scheme gives the best results (see Figure 6).

       ρ(t = 0, x) = 1.0 1 [-1,0] (x) + 0.125 1 [0,1] (x), u(t = 0, x) = 0.
The kinetic velocities are set to λ ± = ±1.7 and λ 0 = 0. The final time equals T f = 0.4 and the As for the previous test cases, the Rusanov splitting is the least accurate and the Bechereau, Osher and Van-Leer splittings are more accurate. We also note that the Osher splitting produces oscillations around the shock wave. A possible explanation would be that the dissipation term is small and the dispersive term is dominant. Indeed, with the [D1Q3] 2 scheme, the wave structure of the BGK system better approximates the original wave structure and the first order error is thus smaller than the one of the [D1Q2] 2 classical scheme.

The Euler system

In this subsection, we give numerical results for the Euler system, using splittings introduced in Subsection 4.4.

Test case 1: Acoustic wave. We consider the domain Ω = [-2, 2] with 4000 cells and a Semi-Lagrangian scheme with 11-th order accuracy in space for the transport step. The initial data is given by

             ρ(t = 0, x) = 1 + 0.1e -x 2 σ , u(t = 0, x) = 0, p(t = 0, x) = ρ(t = 0, x).
with σ = 0.005. The kinetic velocities equal λ ± = ±2 and λ 0 = 0. In Figures 8,9 and 10, we compare the different splittings for various ∆t and various values of ω (ω = 1 in Figure 8, ω = 1.5 in Figure 9 and ω = 2 in Figure 10 left).

From Figures 8 and9, we still observe that the Rusanov splitting is far less accurate than the three others splittings. By comparing the Van-Leer and the Osher splittings, the second one is actually a little bit more accurate.

In Figure 10 right, we compare the results given by the [D1Q3] 3 scheme with the Low-Mach splitting for different time steps and ω = 1. Comparing with Figure 8, the Low-Mach splitting is more accurate for all the waves than the other splittings. On the one hand, for the acoustic waves, the Low-Mach splitting for ∆t = 0.01 is closed to the results given by the other splittings with ∆t = 0.002. On the other hand, the Low-Mach splitting does not dissipate the material wave as done by the other splittings. However, for large time steps, the material wave is a little overestimated. These results confirm that the Low-Mach splitting is accurate for the low Mach number flow. However, the Low-Mach splitting suffers from from a lack of stability. Indeed, if ω is increased, λ (= λ + = -λ -) have to be increased too in order to preserve the stability of the scheme. A similar behavior has already been observed for the Lax-Wendroff splitting for scalar equations. Test case 2: smooth contact. We consider the domain Ω = [-1, 1] with 4000 cells. The smooth contact solution is given by 3 scheme is that the solution with M = 0 is not preserved and the numerical error does not depend on ∆t and u but on ∆t and λ(= λ + = -λ -) or c. The numerical error is therefore disproportionate for the (steady) material wave. A more adapted scheme would be such that the numerical error actually depends on u. Thus, a more accurate scheme would be a scheme for which the choice ∆t = O(1/u) would be possible and would lead to a constant numerical error.

             ρ(t, x) = ρ L w(x, t) + ρ R (1 -w(x, t)), u(t, x) = M, p(t, x) = 1.
First, we compare the different splittings for a given time step ∆t = 0.02 and for various values of u = M. The results are given in Tables 4 and5 where VL stands for the Van-Leer splitting and LM for the Low-Mach splitting. We remark an increase in accuracy when using the Van-Leer splitting and the Osher splitting instead of the Rusanov one. However, all these three splittings (Rusanov, Van-Leer and Osher) generates numerical errors on u, p and ρ, whose magnitude is independent of the value of u. The Low-Mach splitting behaves differently: the numerical errors decrease as u descreases and the material wave is preserved for u = 0. Finally, we observe that considering larger ω (Table 5), results are not significantly better for the Low-Mach splitting: the error is reduced by a factor 1.8 instead of a factor 5 for the other splittings. This may be due to the fact that increase ω allows to decrease diffusion term (not the higher terms) error which is already small in for the low-mach regimes. Let us further study the low Mach number regime. Appropriate scheme for this regime are expected to filter the acoustic wave and enable the use of time steps related to the convection of the density. In order to know whether the Low-Mach splitting is adapted to this aim, we make the time step and the final time vary as follows: ∆t = 0.001/u β , T f = 0.2/u, with β ∈ {1, 1/2, 1/4} and u ∈ 10 -1 , 10 2 , 10 -3 , 10 -4 . Table 6 summarizes the results on the numerical error on the density when using the first order in time scheme (ω = 1) and the coefficient α = 1 (the results are similar with other values of α). We note that the error increases with β. These results seem to indicate that the numerical error (on the density) of the Low-Mach splitting would probably

Rusanov VL Osher LM LM LM α = 1 α = 1.2 α = 1.5 u = 10 -2 ρ(t, x) 0.26 1.0E -1 8.4E -2 1.0E -3 1.0E -3 1.0E -3 u(t, x) 0 3.4E -3
LM for α = 1 u = 10 -1 u = 10 -2 u = 10 -3 u = 10 -4 Error 2.5E -3 2.5E -3 2.5E -3 2.5E -3 λ s 2 0.2 0.02 0.002 λ f 2 20 200 2000
Table 7: Errors for different u and ∆t = 0.001/u for the [D1Q5] 3 scheme with Low-Mach splitting (α = 1) and ω = 1, at final time T f = 0.2/u and we consider the following time step: ∆t = 0.002. On Figure 11 left, we compare the different splittings for ω = 1. As expected, the Rusanov splitting is less accurate than the others on the three waves. Moreover, the Low-Mach splitting oscillates around the contact. These oscillations are coherent with the choice of α = 1 (which leads to a diffusion term equal to zero for the advection part). These oscillations are thus due to dispersive effects which seems depend on λ.

On Figure 11 right, we study the effect of the parameter α of the [D1Q3] 3 scheme with Low-Mach splitting. We note that increasing α enables to add numerical diffusion on the contact wave and to avoid oscillations while keeping the good property since the error also depends on u and not only on λ. 6 Stability issue for the [D1Q3] N scheme

The stability of the schemes presented before is a harder issue than their consistency: it would be difficult to state a general stability analysis. Here, we only provide sufficient conditions on the flux splitting ensuring entropic stability. • ω ∈ [0, 1],

• The matrix Id -

∂F + 0 (λ + -λ 0 ) - ∂F - 0 (λ 0 -λ -)
, and ∂F ± 0 are η -dissipative for the dual of the entropy η call η .

• the transport step is exact,

• the flux vector splitting F ± 0 (U) is entropic, in the sense of Definition 4. Then, the splitting scheme (6)-( 7) with equilibrium (13) and velocities {λ -, λ 0 , λ

+ } N is entropy stable in the sense of x∈Ω H( f n+1 )(x) dx ≤ x∈Ω H( f n )(x) dx, ∀n ∈ N,
where H( f ) denotes a kinetic entropy for the kinetic representation (2).

Proof. First of all, the splitting scheme ( 6)-( 7) with equilibrium ( 13) and velocities {λ -, λ 0 , λ

+ } N rewrites                        ∂ t f + + λ + ∂ x f + = 1 ε ( f eq + -f + ), ∂ t f 0 + λ 0 ∂ x f 0 = 1 ε ( f eq 0 -f 0 ), ∂ t f -+ λ -∂ x f -= 1 ε ( f eq --f -).
In [START_REF] Dubois | Stable lattice Boltzmann schemes with a dual entropy approach for monodimensional nonlinear waves[END_REF] and [START_REF] Dubois | Simulation of strong nonlinear waves with vectorial lattice Boltzmann schemes[END_REF], the authors propose an entropy H for the scheme [D1Q3] N given by

H( f ) = h + ( f + ) + h 0 ( f 0 ) + h -( f -)
where h -, h 0 and h + are three convex functions. We define the dual variables φ

= ∂ U η(U), η (φ) = sup U (φ • U -η(U)),
h + (φ) + h 0 (φ) + h -(φ) = η (φ), λ -h -(φ) + λ 0 h 0 (φ) + λ + h + (φ) = ζ (φ), (37)
and

f eq + = ∂ φ h + , f eq 0 = ∂ φ h 0 , f eq -= ∂ φ h -, where h * j (y) = sup f (y • f -h j ( f ))
denotes the Legendre transform of the h j , we obtain that the equilibrium distribution minimizes the entropy H (see [START_REF] Dubois | Stable lattice Boltzmann schemes with a dual entropy approach for monodimensional nonlinear waves[END_REF]).

At the discrete level, the transport is assumed to be exact. Consequently, the transport step preserves the kinetic entropy

x∈Ω H( f * )(x) dx = x∈Ω H( f n )(x) dx. (38)
Then, the relaxation step gives

x∈Ω H( f n+1 )(x) dx = x∈Ω H (1 -ω) f * + ω f eq (x)dx ≤ (1 -ω) x∈Ω H( f * )(x)dx + ω x∈Ω H( f eq )(x)dx ≤ (1 -ω) x∈Ω H( f n )(x)dx + ω x∈Ω H( f eq )(x)dx ≤ x∈Ω H( f n )(x)dx.
where we first use the convexity of H and the hypothesis ω ∈ [0, 1] and then (38) and the fact that f eq minimizes the entropy. The property is thus established if we are able to build the convex functions h i .

To define functions h i , we replace the decomposition (37) by

       h + (φ) + h 0 (φ) + h -(φ) = η (φ), λ -h -(φ) + λ 0 h 0 (φ) + λ + h + (φ) = ζ -, 0 (φ) + ζ +, 0 (φ) + λ 0 η (φ)
, where ζ ±, 0 (φ) are the dual entropy fluxes of F ± 0 . Thus, ∂ φ ζ ±, 0 (φ) = F ± 0 (U) and we thus have

∂ η (ζ -, 0 (φ) + ζ +, 0 (φ) + λ 0 η (φ)) = F(U) = ∂ φ ζ * (φ) as expected. We can write the relation as follows        h + (φ) + h 0 (φ) + h -(φ) = η (φ), (λ --λ 0 ) h -(φ) + (λ + -λ 0 ) h + (φ) = ζ -, 0 (φ) + ζ +, 0 (φ). A particular choice of h + (φ), h 0 (φ), h -(φ) can be              h + (φ) + h 0 (φ) + h -(φ) = η (φ), (λ --λ 0 )h -(φ) = ζ -, 0 (φ), (λ + -λ 0 )h + (φ) = ζ +, 0 (
φ), We then have the following definition:

                             h -(φ) = - 1 (λ 0 -λ -) ζ -, 0 (φ) , h 0 (φ) = η (φ) -        ζ +, 0 (φ) (λ + -λ 0 ) - ζ -, 0 (φ) (λ 0 -λ -)        , h + (φ) = 1 (λ + -λ 0 ) ζ +, 0 (φ) .
We then look for condition ensuring their convexity. We recall that we have the relation

∂ φφ ζ ±, 0 (φ) = ∂ U F ± 0 (U) ∂ φ η (φ). We thus get ∂ φφ h -(φ) = - 1 (λ 0 -λ -) ∂ U F - 0 (U) ∂ φφ η (φ), (39)
∂ φφ h 0 (φ) = Id - ∂ U F + 0 (U) (λ + -λ 0 ) - ∂ U F - 0 (U) (λ 0 -λ -) ∂ φφ η (φ), (40)
∂ φφ h + (φ) = 1 (λ + -λ 0 ) ∂ U F + 0 (U) ∂ φφ η (φ). ( 41
)
Satisfying the condition of the proposition we obtain the result.

Corollary 6.2. If F ± 0 (U) are diagonalizable in the same basis than F(U) the condition of η dissipation introduced in the theorem (6.1) are equivalent to: the three matrices Id-∂F + 0 (λ + -λ 0 ) -∂F - 0 (λ 0 -λ -) , and ∂F ± 0 are nonnegatives. Proof. To obtain this we take the last line of the previous prove and use the property that the hessian of the dual entropy is diagonal in the eigenvector basis of F(U). Consequently we can diagonalise and obtain the result.

The second condition becomes more classical in some particular cases. For example, let us study the stability of the Rusanov and the upwind splitting for the variable advection equation. Corollary 6.3. For the advection equation (30), scheme (6)-( 7) with exact transport and the Rusanov (31) or the upwind (32) flux splitting with λ -≤ λ 0 ≤ λ + and λ -≤ 0 ≤ λ + is entropy stable on the conditions λ -< a(x) < λ + and ω ∈ [0, 1] .

Proof. Firstly we verify that the flux splitting is entropy for the original system. For the Rusanov splitting (31), the two equations linked to the flux splitting are

∂ t ρ ± λ ± ∂ x a(x) -λ ∓ λ + -λ - ρ = 0.ρ ou f ± ,
and then, we obtain the following entropy equations

∂ t ρ 2 2 ± λ ± ∂ x (a(x) -λ ∓ ) λ + -λ - ρ 2 2 = ∓ λ ± λ + -λ - ρ 2 2
∂ x a(x) ≤ 0 since ∂ x a(x) ≥ 0. Since λ -< a(x) < λ + , the following two conditions hold ∂ ρ F - 0 (ρ, x) = -λ - a(x) -λ + λ + -λ - < 0, and ∂ ρ F + 0 (ρ, x) = λ + a(x) -λ - λ + -λ - > 0.

Moreover,

1 - ∂ ρ F + 0 (ρ, x) λ + - ∂ ρ F - 0 (ρ, x) -λ - = 1 - a(x) -λ - λ + -λ - - a(x) -λ + λ + -λ - = 0.
So the flux splitting is an entropic one and, according to the previous theorem, the scheme is entropy stable. For the upwind splitting (32), we obtain the following entropy equation

∂ t ρ 2 2 + ∂ x χ a(x)<λ 0 (a(x) -λ 0 ) ρ 2 2 = -χ a(x)≶λ 0 ρ 2 2 ∂ x a(x) ≤ 0.
since ∂ x a(x) ≥ 0. We also have ∂ ρ F - 0 (ρ, x) = χ a(x)<λ 0 (a(x) -λ 0 ) < 0, ∂ ρ F + 0 (ρ, x) = χ a(x)>λ 0 (a(x) -λ 0 ) > 0.

Moreover, 1 -∂ ρ F + 0 (ρ, x) λ + -λ 0 -∂ ρ F - 0 (ρ, x) λ 0 -λ - = 1 -a(x) -λ 0 λ + -λ 0 χ a(x)>λ 0 -a(x) -λ 0 λ 0 -λ - χ a(x)<λ 0 = λ +a(x) λ + -λ 0 χ a(x)>λ 0 + a(x) -λ - λ 0 -λ - χ a(x)<λ 0 ≥ 0.

So the flux splitting is entropic and, according to the previous theorem, the scheme is entropy stable.

Let us now study the entropic stability in the particular Burgers equation case, with splittings detailed in Subsection 4.1. Corollary 6.4. For the Burgers equation (26), the scheme (6)-( 7) with the Rusanov splitting (27) or the Osher splitting (28), with λ -≤ λ 0 ≤ λ + and λ -≤ 0 ≤ λ + is entropy stable on the conditions λ -< ρ < λ + and ω ∈ [0, 1] .

Proof. Firstly we verify that the flux splitting is entropy for the original system. The Rusanov splitting (27) on the Burgers equation (26) leads to

∂ t ρ ± λ ± ∂ x ρ 2 2(λ + -λ -) - λ ∓ ρ (λ + -λ -) = 0.
Thus, the following entropy equations is obtained

∂ t ρ 2 2 ± λ ± ∂ x ρ 3 6(λ + -λ -) - λ ∓ ρ 2 2(λ + -λ -) ≤ 0 Moreover, ∂F - 0 = -λ - ρ -λ + (λ + -λ -) < 0, ∂F + 0 = λ + ρ -λ - (λ + -λ -)
> 0, and

1 - ∂F + 0 λ + - ∂F - 0 -λ - = 1 - ρ -λ - λ + -λ - - ρ -λ + λ + -λ - = 0.
Hence, the flux splitting is entropic and, according to the previous theorem, the scheme is entropy stable. For the Osher splitting (28), we obtain the following entropy equations

∂ t ρ 2 2 + ∂ x χ ρ≶λ 0 ρ 3 6 -λ 0 ρ 2 2 ≤ 0.
We have also ∂F - 0 = χ ρ<λ 0 (ρλ 0 ) < 0, ∂F + 0 = χ ρ>λ 0 (ρλ 0 ) > 0. and

1 - ∂F + 0 λ + -λ 0 - ∂F - 0 λ 0 -λ - = 1 - ρ -λ 0 λ + -λ 0 χ ρ>λ 0 - ρ -λ 0 λ 0 -λ - χ ρ<λ 0 = λ + -ρ λ + -λ 0 χ ρ>λ 0 + ρ -λ - λ 0 -λ - χ ρ<λ 0 ≥ 0.
So the flux splitting is entropic and, according to the previous theorem, the scheme is entropy stable.

Remark 6.5. The Lax-Wendroff splittings (33) for the variable advection equation and (29) for the Burgers equation do not satisfy the conditions of the entropic-stability Proposition 6.1, since actually these splittings are not entropic: the sign conditions of the Jacobian flux ∂F ± 0 is not fullfilled.

Conclusion

In [START_REF] Coulette | Palindromic discontinuous Galerkin method for kinetic equations with stiff relaxation[END_REF], an implicit solver for hyperbolic system without matrix inversion has been proposed. It is mainly based on [D1Q2] N vectorial kinetic BGK approximation. The method is competitive since it is easy to parallelize and very low storage, but it suffers from large diffusion/dispersion error. It was partially corrected by a high order in time scheme. However, it still lacks the ability to capture accurately different scales in space or time.

In the present paper, we present a [D1Q3] N kinetic representation, which gives a supplementary degree of freedom to treat this drawback. This representation has been studied theoritically [START_REF] Natalini | A Discrete Kinetic Approximation of Entropy Solutions to Multidimensional Scalar Conservation Laws[END_REF][START_REF] Bouchut | Construction of BGK models with a family of kinetic entropies for a given system of conservation laws[END_REF] and numerically only for the Burgers equation [?]. We here present these different kinetic representations and we compare them with new ones in the context of our implicit kinetic scheme. We show that it enables to reduce errors with a similar computational cost, when the characteristic speeds have large spatial variation. We have also designed one scheme which correctly captures contact waves and low Mach number regime for the full Euler equations.

The main drawback is that most of the schemes are not entropy stable, the diffusion analysis is not sufficient as soon as ω > 1 to ensure stability. Reduced stability conditions should be investigated especially for second-order in time schemes (ω = 2).

The extension to two or three-dimensions must also be investigated, especially for the low Mach number regime of the Euler equations.

Proposition 2 . 3 (

 23 Entropic stability). Let us consider the hyperbolic system (1) with an entropy / entropy-flux pair (η, ζ), satisfying the entropy inequality ∂ t η(U) + ∂ x ζ(U) ≤ 0. Under the following hypotheses:

Figure 1 :

 1 Figure 1: (Burgers equation, test 2) Comparison between different splittings with the first order time scheme. Left ∆t = 0.002. Right ∆t = 0.01. Reference (black), Rusanov (yellow), Upwind (violet), Lax-Wendroff α = 1.0 (green), Lax-Wendroff α = 1.5 (blue).

Figure 2 :

 2 Figure 2: (Variable advection, test 2) Comparison between different flux splittings with the first order in time scheme (ω = 1). Left ∆t = 0.01. Right ∆t = 0.1. Reference (black), Rusanov (violet), Upwind (green), Lax-Wendroff α = 1 (blue), Lax-Wendroff α = 2 (Yellow).

Figure 3 :

 3 Figure 3: (Variable advection, test 2) Comparison between different flux splittings with the second order in time scheme (ω = 2). Left ∆t = 0.01. Right ∆t = 0.1. Reference (black), Rusanov (violet), Upwind (green), Law-Wendroff α = 2 (Yellow).

Figure 4 :

 4 Figure 4: (Variable advection, test 2) Left: Rusanov splitting with second order in time scheme (ω = 2) and ∆t = 0.1. Reference (black), Rusanov for {-1.5, 0.0, 1.5} (violet), Rusanov for {-0.5, 0.0, 1.5} (green). Right: Upwind splitting with second order in time scheme (ω = 2) and and ∆t = 0.1. Reference (black), Upwind for {-1.5, 0.0, 1.5} (violet), Upwind for {-0.5, 1.0, 1.5} (green)

Figure 5 :

 5 Figure 5: (isothermal Euler, test 1) Comparison between different splittings: Reference (black), Rusanov (yellow), HLL (violet), Van Leer (green), Bechereau (blue), Osher (red) for the first option U 0 = 0. Top left: ω = 1 and ∆t = 0.01, Top right: ω = 1.5 and ∆t = 0.01. Bottom: ω = 2 and ∆t = 0.02

Figure 6 :

 6 Figure 6: (isothermal Euler, test 1) Comparison between different splittings: Reference (black), Rusanov (yellow), HLL (violet), Van Leer (green), Bechereau (blue), Osher (red) fot the second option U 0 = 0.98. Left: ω = 1.5 and ∆t = 0.0075. Right: ω = 2 and ∆t = 0.015

Figure 7 :

 7 Figure 7: (isothermal Euler, test 2) Comparison between different splittings : Reference (black), Rusanov (yellow), HLL (violet), Van Leer (green), Bechereau (blue), Osher (red). ω = 1 and ∆t = 0.0075.

Figure 8 :

 8 Figure 8: (Euler, test 1) Results for ω = 1. Left ∆t = 0.002. Right ∆t = 0.005. Reference (black), Rusanov (orange), Van-Leer (green), Osher (violet), AUSM (red).

Figure 9 :

 9 Figure 9: (Euler, test 1) Results for ω = 1.5. Left ∆t = 0.002. Right ∆t = 0.005. Reference (black), Rusanov (orange), Van-Leer (green), Osher (violet), AUSM (red).

Figure 10 :

 10 Figure 10: (Euler, test 1) Left: Results for ω = 2 and ∆t = 0.01. Reference (black), Rusanov (orange), Van-Leer (green), Osher (violet), AUSM (red). Right: Results for [D1Q3] 3 Low-Mach with ω = 1. Reference (black), ∆t = 0.002 (orange), ∆t = 0.005 (green), ∆t = 0.01 (violet)

Figure 11 :

 11 Figure 11: (Euler, test 3) Left: Results for ω = 1.0. Reference (black), Rusanov (orange), Van-Leer (green), Osher (violet), Low-Mach with α = 1 (red). Right: Results for ω = 1.0, ∆t = 0.001. Reference (black), α = 1 (orange), α = 1.2 (green), α = 1.5 (violet), α = 2 (red), α = 1 + u (blue).

Definition 5 .Proposition 6 . 1 .

 561 A matrix A(U) with U a vector of N fucntion is called η-dissipative for the entropy η if the the product A(U)∂ UU η(U) is a nonnegative matrix. Consider the hyperbolic system (1) with an entropy/entropy-flux pair satisfying the following entropy inequality ∂ t η(U) + ∂ x ζ(U) ≤ 0. Let us suppose that the following conditions holds true:

  the Legendre transform of the entropy and ζ (φ) = sup U (φ • F(U)ζ(U)) the dual entropy flux. We thus have ∂ φ ζ (φ) = F(U). On the following conditions

Table 2 :

 2 9E -2 1.00 4.2E -3 1.50 5.5E -4 3.40 1.1E -2 1.30 ∆t = 0.00625 7.9E -3 1.26 1.3E -3 1.70 5.3E -5 3.40 3.9E -3 1.50 (Variable advection, test 1) Convergence in time for the first order time scheme (ω = 1) [D1Q3] 1 and different flux splittings. LW stands for Law-Wendroff.

	considering the second

Table 4 :

 4 (Euler, test 2) Errors for the different splittings for ω = 1 and λ = 2. VL: Van Leer,

					6.0E -7	0	0	0
		p(t, x)	0	5.0E -4 4.3E -8	0	0	0
		ρ(t, x)	0.26	1.0E -1 8.4E -2 1.0E -5 1.0E -5 1.0E -5
	u = 10 -4	u(t, x)	0	3.4E -3 6.0E -7	0	0	0
		p(t, x)	0	5.0E -4 4.3E -8	0	0	0
		ρ(t, x)	0.26	1.0E -1 4.8E -2	0.0	0.0	0.0
	u = 0	u(t, x)	0	3.4E -3 6.0E -7	0	0	0
		p(t, x)	0	5.0E -4 4.3E -8	0	0	0
	LM: Low-Mach.
				.
			Rusanov	VL	Osher	LM	LM	LM
						α = 1 α = 1.2 α = 1.5
		ρ(t, x)	0.13	4.8E -2 3.8E -2 5.5E -4 5.5E -4 5.5E -4
	u = 10 -2	u(t, x)	0	1.1E -3 1.2E -7	0	0	0
		p(t, x)	0	1.0E -4 1.7E -8	0	0	0
		ρ(t, x)	0.13	4.8E -2 3.8E -2 5.5E -6 5.5E -6 5.5E -6
	u = 10 -4	u(t, x)	0	1.1E -3 1.2E -7	0	0	0
		p(t, x)	0	1.0E -4 1.7E -8	0	0	0
		ρ(t, x)	0.13	4.8E -2 3.8E -2	0.0	0.0	0.0
	u = 0	u(t, x)	0	1.1E -3 1.2E -7	0	0	0
		p(t, x)	0	1.0E -4 1.7E -8	0	0	0

Table 5 :

 5 (Euler, test 2) Errors for the different splittings for ω = 1.5 and λ = 2.5. VL: Van Leer, LM: Low-Mach.

If A = PDP -1 ∈ M N (R) where P ∈ M N (R) is invertible and D = diag(d 1 , . . . , d N ) ∈ M N (R) is diagonal, then |A| = P|D|P -1 with |D| = diag(|d 1 |, . . . , |d N |).

be homogeneous to

The first term is the diffusion error term coming from the Lax-Wendroff flux splitting of the advection part. The second error term is coherent with the observations: for a given time step, the error decreases with u. Consequently, with this Low-Mach splitting, we can choose ∆t proportional to 1/ √ u for a given final time T f and keep the same error independently of the value of u. However, we can not have a constant error if we increase at the same time the final T f as 1 u (as defined previously).

Error u = 10 -1 u = 10 -2 u = 10 -3 u = 10 -4 β = 1 2.5E -3 0.12 0.25 0.20

The main problem raised by the previous remark comes from the choice of λ. Indeed, λ is determined such as to be larger than the acoustic speed. In the density equation, this choice leads to a disproportional dispersive term. To overcome this difficulty, we consider the [D1Q5] 3 scheme introduced in (36) that enables to have a kinetic velocity λ s , adapted to the slow scale, and an other one λ f , adapted to the acoustic scale. In Table 7, we have performed the same study for the [D1Q5] 3 scheme with the Low-Mach splitting with the same ∆t and T f as previously with β = 1 and, for each value of u, we consider well-adapted λ s and λ f . These results confirm that [D1Q5] 3 scheme with the Low-Mach splitting enables to choose the time step comparing to u while obtaining a constant error (independent of u). Thus, this scheme allows to filter the acoustic part. These results confirm also the form of the error

since ∆t ∝ 1/u, T f ∝ 1/u and λ s ∝ u enable to obtain a constant error. The fact that the velocity λ f have to be increased is dictated by stability condition. The numerical results suggest a stability condition for the [D1Q5] 3 scheme with Low-Mach splitting like

Test case 3: Sod problem. The last test case proposed in this subsection is the classical Sod problem with contact discontinuity, shock and rarefaction waves. The domain is given by Ω = [-1, 1] with 4000 cells and a first order in space Semi-Lagrangian scheme for the transport step. The initial condition is given by: