Christophe Chalons
email: chalons@math.jussieu.fr

Damien Kah
email: damienk@stanford.edu

Marc Massot
email: marc.massot@ecp.fr-

BEYOND PRESSURELESS GAS DYNAMICS : QUADRATURE-BASED VELOCITY MOMENT MODELS

Keywords: Quadrature-based moment methods, gas-particle flows, kinetic theory, particle trajectory crossing, entropic measure solution, frontier of the moment space subject classifications 76N15 (35L65 65M99 76M25 82C40)

Following the seminal work of F. Bouchut on zero pressure gas dynamics, extensively used for gas particle-flows, the present contribution investigates quadrature-based velocity moments models for kinetic equations in the framework of the infinite Knudsen number limit, that is, for dilute clouds of small particles where the collision or coalescence probability asymptotically approaches zero. Such models define a hierarchy based on the number of moments and associated quadrature nodes, the first level of which leads to pressureless gas dynamics. We focus in particular on the four moment model where the flux closure is provided by a two-node quadrature in the velocity phase space and provides the right framework for studying both smooth and singular solutions. The link with both the kinetic underlying equation as well as with zero pressure gas dynamics, i.e. the dynamics at the frontier of the moment space of order four, is provided. We define the notion of measure solutions and characterize the mathematical structure of the resulting system of four PDEs. We exhibit a family of entropies and entropy fluxes and define the notion of entropic solution. We study the Riemann problem and provide entropic solutions in particular cases. This leads to a rigorous link with the possibility of the system of macroscopic PDEs to allow particle trajectory crossing (PTC) in the framework of smooth solutions. Generalized δ-shock solutions resulting from Riemann problem are also investigated. Finally, using a kinetic scheme proposed in the literature in several areas, we validate such a numerical approach and propose a dedicated extension at the frontier of the moment space in the framework of both regular and singular solutions. This is a key issue for application fields where such an approach is extensively used.

Introduction

The physics of particles and droplets in a carrier gaseous flow field are described in many applications (fluidized beds, spray dynamics, alumina particles in rocket boosters, ...) by a number density function (NDF) satisfying a kinetic equation introduced by [Williams (1958)]. Solving such a kinetic equation relies either on a sample of discrete numerical parcels of particles through a Lagrangian-Monte-Carlo approach or on a moment approach resulting in a Eulerian system of conservation laws on velocity moments eventually conditioned on size. In the latter case investigated in the present contribution, the main difficulty for particle flows with high Knudsen numbers (i.e. weakly collisional flows), where the velocity distribution can be very far from equilib-rium, is the closure of the convective transport at the macroscopic level. One way to proceed is to use quadrature-based moment methods where the higher-order moments required for closure are evaluated from the lower-order transported moments using quadratures in the form of a sum of Dirac delta functions in velocity phase space (see [Yuan & Fox (2010), Kah (2010)] and the references therein).

Such an approach also allows for a well-behaved kinetic numerical scheme in the spirit of Bouchut [START_REF] Bouchut | [END_REF]] (See references from [Laurent (2002), Fréret & al (2008), de Chaisemartin (2009), [START_REF] Massot | [END_REF], de Chaisemartin & al (2008)] to [START_REF] Kah | Eulerian quadrature-based moment models for dilute polydisperse evaporating sprays[END_REF], [START_REF] Fréret | On the role of preferential segregation in flame dynamics in polydisperse evaporating sprays[END_REF], Yuan & Fox (2010)]) where the fluxes in a cell-centered finite-volume formulation are directly evaluated from the knowledge of the quadrature abscissas and weights with guaranteed realizability conditions and singularity treatment. Such a quadrature approach and the related numerical methods have been shown to be able to capture particle trajectory crossing (PTC) in a Direct Numerical Simulation (DNS) context, where the distribution in the exact kinetic equation remains at all times in the form of a sum of Dirac delta functions. Such methods can be extended to partially high-order numerical schemes ( [Vikas et al. (2010)]).

In another component of the literature devoted to multiphase semiclassical limits of the Schrödinger equation [Jin & Li (2003), Gosse et al. (2003), [START_REF] Gosse | [END_REF]], the series of Wigner measures obtained from the Wigner transform for studying the semiclassical limits can be shown to converge towards a measure solution of the Liouville equation. Such an equation naturally unfolds the caustics and can generate the proper multiphase solutions globally in time. Two approaches have been used to solve this equation with a moment approach, either the Heaviside closure [Brenier & Corrias (1998)] as it is called in [Jin & Li (2003), Gosse et al. (2003)], or, the one which is related to the present work, the delta closure (see [Jin & Li (2003), Gosse et al. (2003)] and references therein). It leads to a weakly hyperbolic systems of conservation laws by taking moments of a Liouville equation exactly identical to the Williams-Boltzmann equation studies in gas-particle flows previously mentioned. Such approaches naturally degenerate towards the pressureless gas system of equation in the context of monokinetic velocity distributions [START_REF] Massot | [END_REF], [START_REF] Kah | Eulerian quadrature-based moment models for dilute polydisperse evaporating sprays[END_REF], Kah (2010), Runborg (2000)].

Numerical algorithms in order to simulate such systems of conservations laws with the related delta closure or quadrature-based closure have been proposed in [Jin & Li (2003), Gosse et al. (2003)] and [START_REF] Desjardins | A quadrature-based moment method for dilute fluid-particle flows[END_REF]] independently, from the work for [START_REF] Bouchut | [END_REF], Laurent (2002)] using naturally kinetic scheme with finite volume methods. However, many issues are still to be tackled in order to reach fully high order numerical schemes that preserve the vector of moments inside or at the frontier of the moment space, thus leading to several possibilities of degeneration from a given number of abscissas to a lower one. In fact, such models are meant to capture a given level of complexity in the phase space which is fixed in advance by the number of moments and related quadrature nodes. In some particular situations, for perfectly controlled dynamics, it can be guaranteed that the solutions will remain smooth and consist in free boundary value (contact discontinuities) problem associated with switches between various numbers of quadrature abscissas. However, in most cases the numerical schemes have to tackle the possibility of singular solutions when the dynamics complexity goes beyond the one allowed by the model. In such cases the solution of the resulting system of PDEs is the viscosity solution and does not reproduce the exact dynamics in phase space and measure so-lutions are expected, for which we need a precise framework. More specifically, even if for the pressureless gas system, [Bouchut (1994)] had set the correct mathematical background to define general entropic solutions, such a work had not yet been performed for higher order moment methods in the cited publications and no rigorous link has been provided between these and zero pressure gas dynamics. This is the purpose of the present contribution for both theoretical and numerical points of view.

The paper is organized as follows. First we introduce the quadrature-based or delta closure velocity moment models for kinetic equations and focus on the four moment model in order to generalize to what can be done to higher orders but which would be difficult to expose due to algebra complications. The behavior at the frontier of the moment space is characterized as well as the mathematical structure of the system of conservation laws. We then define entropy conditions and provide, for smooth solution, the one-to-one kinetic-macroscopic relation. We then tackle the Riemann problem and define entropy measure solutions. Three examples of piecewise linear and singular solutions are then provided for which we rigorously identify the entropic character of the solution and which are then reproduced numerically.

Quadrature-based velocity moment models for kinetic equations

Consider the solution f = f (t,x,v) of the free transport kinetic equation

∂ t f + v∂ x f = 0, t > 0, x ∈ R, v ∈ R, f (0,x,v) = f 0 (x,v) (2.1)
The exact solution is given by f

(t,x,v) = f (0,x -vt,v) = f 0 (x -vt,v). Defining the i-order moment M i = v f (t,x,v)v i dv, i = 1,.
..,N, N ∈ N, the associated governing equations are easily obtained from (2.1) after multiplication by v i and integration over v, and write

∂ t M i + ∂ x M i+1 = 0, i ≥ 0.
For the sake of simplicity, but without any restriction, we will focus our attention hereafter on the four-moment model

       ∂ t M 0 + ∂ x M 1 = 0, ∂ t M 1 + ∂ x M 2 = 0, ∂ t M 2 + ∂ x M 3 = 0, ∂ t M 3 + ∂ x M 4 = 0.
(2.

2)

It will be convenient to write (2.2) under the following abstract form

∂ t M + ∂ x F(M) = 0, (2.3) with M = (M 0 ,M 1 ,M 2 ,M 3 ) t and F(M) = (M 1 ,M 2 ,M 3 ,M 4 ) t .

Quadrature inside the moment space

This model is closed provided that M 4 is defined as a function of M. In quadrature-based moment methods, the starting point to define this closure relation consists in representing the velocity distribution of f (t,x,v) by a set of two Dirac delta functions, that is a two-node quadrature :

f (t,x,v) = ρ 1 (t,x)δ(v -v 1 (x,t)) + ρ 2 (t,x)δ(v -v 2 (x,t)),
(2.4)

where the weights ρ 1 (t,x) > 0, ρ 2 (t,x) > 0 and the velocity abscissas v 1 (t,x), v 2 (t,x) are expected to be uniquely determined from the knowledge of M(x,t). Dropping the (x,t)-dependance to avoid cumbersome notations, such a function f has exact moments of order i = 0,...,4 given by ρ 1 v i 1 + ρ 2 v i 2 . The next step then naturally consists in setting

M 4 = ρ 1 v 4 1 + ρ 2 v 4 2 (2.5)
where ρ 1 , ρ 2 and v 1 , v 2 are defined from M by the following nonlinear system :

       M 0 = ρ 1 + ρ 2 , M 1 = ρ 1 v 1 + ρ 2 v 2 , M 2 = ρ 1 v 2 1 + ρ 2 v 2 2 , M 3 = ρ 1 v 3 1 + ρ 2 v 3 2 .
(2.6)

At last, it remains to prove that this system is well-posed, which is the matter of the next proposition. We refer to [Jin & Li (2003), Gosse et al. (2003), [START_REF] Desjardins | A quadrature-based moment method for dilute fluid-particle flows[END_REF]] for the proof.

Proposition 2.1. System (2.3)-(2.5)-(2.6) is well-defined on the convex phase space Ω, also called the moment space, given by

Ω = {M = (M 0 ,M 1 ,M 2 ,M 3 ) t ,M 0 > 0,M 0 M 2 -M 2 1 > 0}. Moreover, setting U = (ρ 1 ,ρ 2 ,ρ 1 v 1 ,ρ 2 v 2 ) t , the function U = U(M)
is one-to-one and onto as soon as we set for instance v 1 > v 2 . Moreover we have 0 < ρ 1 < M 0 and 0 < ρ 2 < M 0 .

Proposition 2.1 can be extended to the more general case of a 2k-moment models, k > 1. The velocity distribution is represented in this situation by a set of k Dirac delta functions, leading to

M i = k j=1 ρ j v i j , i = 0,...,2k -1, and M 2k = k j=1 ρ j v 2k j .

Hyperbolic structure inside the moment space

The two-moment model, corresponding to k = 1 (one-node quadrature) writes

∂ t ρ + ∂ x ρv = 0, ∂ t ρv + ∂ x ρv 2 = 0,
which is the well-known pressureless gas dynamics system. Recall that this model is weakly hyperbolic (the jacobian matrix is not diagonalizable) with v as unique eigenvalue, the characteristic field being linearly degenerate. Since there can be areas in the solution where a single quadrature node is sufficient at the frontier of the moment space in order to describe the dynamics, the solution in such zones will satisfy the previous system of two conservation laws. However we will first work inside the moment space and leave the behavior at the frontier for the next subsection.

Actually, we will observe in the course of the next section that the four-moment model (2.3) is equivalent for smooth solutions (only) to two decoupled pressureless gas dynamics systems associated with (ρ 1 ,ρ 1 v 1 ) and (ρ 2 ,ρ 2 v 2 ) respectively. Then (2.3) is expected to admit two eigenvalues v 1 and v 2 and to be weakly hyperbolic with linearly degenerate characteristic fields, as stated in the following proposition.

Proposition 2.2. ( [Jin & Li (2003), Gosse et al. (2003)]) System (2.3)-(2.5)-(2.6) is weakly hyperbolic on Ω and admits the two eigenvalues v 1 and v 2 , v 1 = v 2 . The associated characteristic fields are linearly degenerate.

Proof. For the sake of completeness, we propose here a direct proof of the eigenvalues v 1 and v 2 of (2.3)-(2.5)-(2.6). By (2.6), we first easily get

       M 0 = ρ 1 + ρ 2 , M 1 -v 1 M 0 = ρ 2 (v 2 -v 1 ), M 2 -v 1 M 1 = ρ 2 v 2 (v 2 -v 1 ), M 3 -v 1 M 2 = ρ 2 v 2 2 (v 2 -v 1 ),
and then, setting

σ 0 = v 1 v 2 and σ 1 = -(v 1 + v 2 ), M 0 M 1 M 1 M 2 σ 0 σ 1 = - M 2 M 3 .
This system is invertible in the phase space Ω (M 0 M 2 -M 2 1 = 0) and uniquely defines σ 0 and σ 1 with respect to M :

σ 0 σ 1 = 1 M 0 M 2 -M 2 1 M 1 M 3 -M 2 2 M 1 M 2 -M 0 M 3 . (2.7)
Then, we have

M 4 = ρ 1 v 4 1 + ρ 2 v 4 2 = ρ 1 v 3 1 v 1 + ρ 2 v 3 2 v 2 = (ρ 1 v 3 1 + ρ 2 v 3 2 )(v 1 + v 2 ) -(ρ 1 v 2 1 + ρ 2 v 2 2 )v 1 v 2 = -M 2 σ 0 -M 3 σ 1 , (2.8)
which finally gives M 4 with respect to M. The Jacobian matrix J = ∇ M F is given by

J =     0 1 0 0 0 0 1 0 0 0 0 1 a b c d     with        a = ∂ M0 M 4 , b = ∂ M1 M 4 , c = ∂ M2 M 4 , d = ∂ M3 M 4 .
Using (2.7) and (2.8), the calculations of the last row coefficients eventually lead to

J =     0 1 0 0 0 0 1 0 0 0 0 1 -σ 2 0 -2σ 0 σ 1 -2σ 0 -σ 2 1 -2σ 1     .
Finally, the characteristic polynomial p(λ) of J is easily shown to equal

p(λ) = (λ -v 1 ) 2 (λ -v 2 ) 2 .
This concludes the proof. Propositions 2.1 and 2.2 show that System (2.3)-(2.5)-(2.6) is well-defined and weakly hyperbolic only on Ω, which gives in particular v 1 = v 2 in the interior of the moment space. At a first sight, this might appear to be restrictive in the sense that one of the main objectives of the model is to allow particle trajectory crossing, that is in particular to deal with initial data consisting of two colliding particle packets such that v 1 = v 2 at each point initially (see for instance Section 7). Thus, in the last part of the present section, we characterize the behavior at the frontier Γ of the moment space when M 0 > 0 :

Γ = {M = (M 0 ,M 1 ,M 2 ,M 3 ) t ,M 0 > 0,M 0 M 2 -M 2 1 = 0}.

Behavior at the frontier of the moment space

As mentioned previously, it is rather natural to envision the coexistence, in a single smooth moment solution, of zones where the number of quadrature nodes are different. More specifically, we will examine the coexistence of zones where only one quadrature node is needed (v 1 = v 2 ), that is where M = (M 0 ,M 1 ,M 2 ,M 3 ) t with M 0 > 0 and M 0 M 2 -M 2 1 = 0, and zones where M 0 M 2 -M 2 1 > 0 which are inside the moment space Ω, whereas the vector of moments are smooth everywhere.

After easy calculations in terms of ρ 1 , ρ 2 , v 1 and v 2 , the latter equality M 0 M 2 -M 2 1 = 0 writes ρ 1 ρ 2 (v 1 -v 2 ) 2 = 0, so that if the vector (ρ 1 ,ρ 2 ,ρ 1 v 1 ,ρ 2 v 2 ) t exists, this actually corresponds to the case v := v 1 = v 2 (still under the assumption ρ 1 = 0 and ρ 2 = 0) or to the case where one of the weights is zero. We also note that in both cases M k = M 0 v k , whatever k in this case, so that the whole set of moments should be provided once M 0 and M 1 are given, in close connection to the case of pressureless gas dynamics. There are in fact two possibilities :

-either M = (M 0 ,M 1 ,M 2 ,M 3 ) t is such that M 0 M 3 -M 1 M 2 = 0 : in this case (2.6) cannot be solved and the vector (ρ 1 ,ρ 2 ,ρ 1 v 1 ,ρ 2 v 2 ) t does not exist, -or M = (M 0 ,M 1 ,M 2 ,M 3 ) t is also such that M 0 M 3 -M 1 M 2 = 0 : in this case (2.6
) can be solved and we have v = v 1 = v 2 = M 1 /M 0 , together with ρ 1 and ρ 2 defined by the one-parameter equation ρ 1 + ρ 2 = M 0 . As we will see just below, the choice ρ 1 = ρ 2 = M 0 /2 is the most natural one when we have to deal with an isolated point at the frontier of the moment space.

In order to justify the choice

ρ 1 = ρ 2 = M 0 /2, we first observe that both condi- tions M 0 M 2 -M 2 1 = 0, M 0 M 3 -M 1 M 2 = 0, are equivalent to conditions e = 0, q = 0,
, where

e = M 0 M 2 -M 2 1 , q = (M 3 M 2 0 -M 3 1 ) -3M 1 (M 0 M 2 -M 2 1 ),
and we consider ρ 1 , ρ 2 , v 1 and v 2 as functions of (M 0 ,M 1 ,q,e) with M 0 > 0, and e > 0 1 . We then propose to study the asymptotic behavior of these functions when e → 0 + , considering that M 0 > 0, M 1 and q are fixed. Note that Γ = {(M 0 ,M 1 ,e,q) t ,M 0 > 0,e = 0}. We get the following result.

Lemma 2.3. Let be given M 0 > 0, M 1 and q. Then we have

lim e→0 + ρ 2 =    M 0 if q > 0, 0 if q < 0, M0 2 if q = 0, lim e→0 + ρ 1 =    0 if q > 0, M 0 if q < 0, M0 2 if q = 0,
1 The definitions of e and q naturally comes out after noticing that setting

ρ 1 = ρ 1 M 0 , ρ 2 = ρ 2 M 0 , v 1 = v 1 -M 1 M 0 , v 2 = v 2 -M 1 M 0 , solving (2.6) is equivalent to solving        1 = ρ 1 + ρ 2 , 0 = ρ 1 v 1 + ρ 2 v 2 , e = ρ 1 v 2 1 + ρ 2 v 2 2 , q = ρ 1 v 3 1 + ρ 2 v 3 2 , with e = (M 0 M 2 -M 2 1 )/M 2 0 and q = ((M 3 M 2 0 -M 3 1 ) -3M 1 (M 0 M 2 -M 2 1 ))/M 3 0 . lim e→0 + v 2 =    M1 M0 if q > 0, -∞ if q < 0, M1 M0 if q = 0, lim e→0 + v 1 =    +∞ if q > 0, M1 M0 if q < 0, M1 M0 if q = 0, lim e→0 + ρ 2 v 2 =    M 1 if q > 0, 0 if q < 0, M1 2 if q = 0, lim e→0 + ρ 1 v 1 =    0 if q > 0, M 1 if q < 0, M1 2 if q = 0.
Proof. The admissible change of variables (M 0 ,M 1 ,M 2 ,M 3 ) → (M 0 ,M 1 ,e,q) allows to write after easy calculations

v 1 = M 1 M 0 + q + q 2 + 4e 3 2M 0 e , ρ 2 = M 0 e(v 1 M 0 -M 1 )
q2 + 4e3 ,

and if q = 0 v 1 = M 1 M 0 + q M 0 e 1 + sign(q) 1 + 4e 3 /q 2 2 , ρ 2 = M 0 e q v 1 M 0 -M 1 sign(q) 1 + 4e 3 /q 2 ,
where we have set

sign(q) = 1 if q > 0, -1 if q < 0.
It is then an easy matter to get the expected results distinguishing between the three cases q < 0, q > 0 and q = 0. It is then clear by a continuity argument that the proposed choice ρ 1 = ρ 2 = M 0 /2 when e = q = 0 is actually natural.

An important consequence of this lemma is that in the half plane e > 0, the region close to the frontier Γ for a non-zero q corresponds to abscissas going to infinity with arbitrary small weights. Moreover, when the velocity distributions at the kinetic level have compact support in the initial distribution, such a property will be preserved in the dynamics of the system and we want to be able to switch continuously from twonode to one-node quadrature without pathological behavior on abscissas and weights.

Let us provide a first example where such a behavior is present. We consider a path in the moment space parametrized by the variable x, such that ρ 1 = x 3 , ρ 2 = 1, v 1 = 1/x and v 2 = 0. As x approaches zero, the smooth moment vector M = (1 + x 3 ,x 2 ,x,1) t has a very regular limit at the frontier of the moment space along the lines presented before with an unbounded abscissa. Indeed we have here e = x and q = 1 -x 3 approaches the fixed non-zero value of 1. Note that if we replace the first weight by ρ 1 = x 4 , we still have an unbounded abscissa even if we converge toward the point (0,0) in the (e,q) plane (e = x 2 ,q = x(1 -x 4 )).

We thus have to find a framework in a subset of the plane (e,q) such that the limits are better behaved. A natural choice presented above is the line q = 0 but it is too restrictive. In order to naturally introduce the relevant subset of Ω, let us consider the other example with ρ 1 = αx β , ρ 2 = 1, v 1 = γx δ and v 2 = 0, with α > 0, γ > 0, β ≥ 0 and δ ≥ 0. As x approaches zero, the moment vector M = (1 + αx β ,αγx β+δ ,αγ 2 x β+2δ ,αγ 3 x β+3δ ) t reaches the frontier Γ of the moment space 2 . Two cases are interesting; firstly when β = 0, we reach the point (0,0) in the (e,q) plane asymptotically along the line q/e 3/2 = (1 -α)/α 1/2 and no weight is approaching zero, whereas the two abscissas are joining (see formulas in the proof above). Secondly, when δ = 0, one of the weights is reaching zero, whereas the two abscissas remain different at a distance of γ at the limit and we reach the point (0,0) in the (e,q) plane asymptotically along the line q/(M 0 e) = γ at the limit x → 0 (see again formulas in the proof above). We will prove in the following proposition that the proper framework is a symmetric cone in the (e,q) plane centered at the point (0,0) corresponding the |q/(M 0 e)| ≤ η, where η is a measure of the maximal distance allowed between the two abscissas.

Definition 2.4 (Regular path).

We define a regular path parametrized by x in the moment space, M x = (M 0 x ,M 1 x ,M 2 x M 3 x ) which admits a limit as x goes to zero and is at least C 1 up to the limit M 0 . Moreover, we define its reduced second and third order moments

e x = M 0 x M 2 x -M 2 1 x and q x = (M 3 x M 0 x -M 3 1 x ) -3M 1 x (M 0 x M 2 x -M 2 1 x ).
Its limit further satisfies e 0 = e x=0 = 0 and we assume e x>0 > 0, M 0 x > ν > 0, |q x /(M 0 x e x )| ≤ η, where η > 0.

Proposition 2.5. We then have the following properties :

• lim x→0 q x = q 0 = 0.

• the weights and abscissas admit limits

ρ i0 = lim x→0 ρ ix , v i0 = lim x→0 v ix . If
we assume that ρ i0 > 0 for both i, then v 1 0 = v 2 0 , or, if one weight approaches zero, such as ρ 1 0 = 0 then we have |v 1 0 -v 2 0 | ≤ η and η is then a bound on the distance between the two abscissas.

• M 0 = (M 0 0 ,M 1 0 ,M 2 1 0 /M 0 0 ,M 3 1 0 /M 2 0 0 ) t . Proof. It is first clear that lim x→0 q x = q 0 = 0 since |q x | ≤ ηM 0 x e
x and e 0 = 0. Then, easy calculations give

v 1 -v 2 = q 2 M 2 0 e 2 + 4e
M 2 0 so that denoting l = lim x→0 |q x |/(M 0x e x ) ≥ 0, we clearly have v 1x -v 2x → l when x → 0 and η represents an upper bound for v 1 0 -v 2 0 . Let us now distinguish between the cases l > 0 and l = 0. We first note the following expression for q/(M 0 e) :

q M 0 e = (ω 1 -ω 2 )(v 2 -v 1 ), ω i = ρ i /M 0 , ω 1 + ω 2 = 1. If l > 0, one can write |ρ 1 x -ρ 2 x | = | q x M 0 x e x | × 1 |v 1 x -v 2 x | × M 0 x ,
and this quantity clearly tends to M 0 0 as x goes to zero. Which means that one weight approaches zero, ρ 1 0 = 0 or ρ 2 0 = 0, and the other one M 0 0 (with

ρ 1 x + ρ 2 x = M 0 x ).
If l = 0, it is clear by the following formula for v 1 (see the proof above)

v 1 = M 1 M 0 + q + q 2 + 4e 3 2M 0 e that v 1 0 = v 2 0 = M 1 0 /M 0 0 .
Using now the definition of ρ 1 and ρ 2 (see again the proof above), one easily get

ω 2 = ρ 2 M 0 = 1 + 4e 3 /q 2 + sign(q) 2 1 + 4e 3 /q 2 , ω 1 = ρ 1 M 0 = 1 + 4e 3 /q 2 -sign(q)
2 1 + 4e3 /q 2 , so that both weights have limits depending on the limit of q/e 3/2 . Clearly, M 0 = (M 0 0 ,M 1 0 ,M 2 1 0 /M 0 0 ,M 3 1 0 /M 2 0 0 ) t , which completes the proof. A very important consequence of the previous proposition is the fact that along smooth paths inside the proposed cone which reach the point (0,0) in the (e,q) plane, the flux introduced in equation 2.3 is regular up to the frontier of the moment space, even if the mapping of M onto U is not 3 Proposition 2.6. For any regular path in the moment space satisfying the assumptions of the previous proposition, that is living in the proper cone in the (e,q) plane and reaching smoothly the point (0,0), the flux F(M) is continuous up to the frontier of the moment space and C 1 at (0,0) in any direction inside the proposed cone.

Proof. The proof is rather straightforward when one has noticed the two equations, the first of which is the expression (2.8) of M 4 = -σ 0 M 2 -σ 1 M 3 as a function of M 2 , M 3 , σ 0 and σ 1 , and the second is the expression of σ 0 and σ 1 :

    σ 0 σ 1     =      q M 0 e M 1 M 0 + M 1 M 0 2 - e M 2 0 - q M 0 e -2 M 1 M 0      .
Finally, following the same lines for the evaluation of the Jacobian matrix of the flux (see matrix J in section 2.2), it becomes clear that the flux is continuous and continuously differentiable in any direction inside the proposed cone. Besides, it can be easily seen that the expression of the flux as a function of (M 0 ,M 1 ,e,q) t becomes

M 4 M 0 = - q M 0 e M 1 M 0 + M 1 M 0 2 - e M 2 0 e M 2 0 + M 1 M 0 2 + q M 0 e + 2 M 1 M 0 q M 3 0 + M 1 M 0 3 + 3 M 1 M 0 e M 2 0 .
M 4 /M 0 tends to (M 10 /M 00 ) 4 when e goes to 0 + which concludes the proof. Remark 2.1. Let us emphasize that in the various configurations we have proposed when the convergence toward the frontier of the moment space does not lie inside the cone |q x /(M 0 x e x )| ≤ η in the (e,q) plane, the flux can dramatically loose regularity. It can either have a limit without being differentiable or even not have a limit at all. The impact of the previous proposition thus becomes clear and sets the proper framework for solutions which will reach the frontier of the moment space. Remark 2.2. In the case M 0 > 0, at the frontier of the moment space within the previous proposed framework, we have e = 0, q = 0; the model is made of the two unknowns M 0 and M 1 and then degenerates to the usual pressureless gas dynamics which is weakly hyperbolic with a single eigenvalue v = M 1 /M 0 . It should be noticed that for smooth solutions, the last two equations of system 2.2 on

M 2 = M 2 1 /M 0 and M 3 = M 3 1 /M 2 0 are still satisfied with M 4 = M 4
1 /M 3 0 coherent with the previous limit obtained for the flux. As a consequence, we can notice, that at least for smooth solutions, the system of partial differential equations (2.2) can describe the dynamics inside and at the frontier of the moment space.

Entropy conditions

In this section, we will work in the interior of the moment space and we exhibit natural entropy inequalities for the following small viscosity system associated with (2.2) :

       ∂ t M 0 + ∂ x M 1 = ε∂ xx M 0 , ∂ t M 1 + ∂ x M 2 = ε∂ xx M 1 , ∂ t M 2 + ∂ x M 3 = ε∂ xx M 2 , ∂ t M 3 + ∂ x M 4 = ε∂ xx M 3 , (3.1)
which gives in condensed form

∂ t M + ∂ x F(M) = ε∂ xx M. (3.2)
Throughout this section, we will consider smooth solutions only. We thus have

∂ t M + J∂ x M = ε∂ xx M with J = ∇ M F.
Setting A = ∇ U M, we then get

∂ t U + A -1 JA∂ x U = εA -1 ∂ x (A∂ x U). (3.3)
Our objective is to prove that

∂ t η + ∂ x q ≤ ε∂ xx η, (3.4)
for a natural choice of couple (η,q) given by

η = ρ 1 S(v 1 ) + ρ 2 S(v 2 ), q = ρ 1 v 1 S(v 1 ) + ρ 2 v 2 S(v 2 ). (3.5)
Here S denotes a convex function from R to R, and we will especially consider the case where S(v) = v 2α , α ≥ 0. Of course, the densities ρ 1 , ρ 2 and velocities v 1 , v 2 involved in (3.5) are naturally defined by means of the one-to-one and onto function U = U(M). In the following and with a little abuse in the notations, we will consider without distinction η and q as functions of M or U.

We first observe

∂ t η + ∂ x q = ∂ t η(U) + ∂ x q(U) = ∇ U η∂ t U + ∇ U q∂ x U = ∇ U η{-A -1 JA∂ x U + εA -1 ∂ x (A∂ x U)} + ∇ U q∂ x U.
The following two lemmas, the proofs of which are left to the reader, will be useful in order to estimate the entropy dissipation rate D defined by

D = ∇ U η{-A -1 JA∂ x U + εA -1 ∂ x (A∂ x U)} + ∇ U q∂ x U.
Lemma 3.1. The matrices J and A and A -1 JA are given by

J =     0 1 0 0 0 0 1 0 0 0 0 1 -v 2 1 v 2 2 2v 1 v 2 (v 1 + v 2 ) -2v 1 v 2 -(v 1 + v 2 ) 2 2(v 1 + v 2 )     , A =     1 1 0 0 0 0 1 1 -v 2 1 -v 2 2 2v 1 2v 2 -2v 3 1 -2v 3 2 3v 2 1 3v 2 2     , A -1 JA =     0 0 1 0 0 0 0 1 -v 2 1 0 2v 1 0 0 -v 2 2 0 2v 2     .
Lemma 3.2. The gradients ∇ U η and ∇ U q are given by

∇ U η =     S(v 1 ) -v 1 S ′ (v 1 ) S(v 2 ) -v 2 S ′ (v 2 ) S ′ (v 1 ) S ′ (v 2 )     t , ∇ U q =     -v 2 1 S ′ (v 1 ) -v 2 2 S ′ (v 2 ) S(v 1 ) + v 1 S ′ (v 1 ) S(v 1 ) + v 1 S ′ (v 1 )     t ,
and we have

∇ U q = ∇ U ηA -1 JA.
Before going on, let us make the following two remarks. We first note that thanks to the first lemma, (3.3) with ε = 0 inside the moment space gives :

       ∂ t ρ 1 + ∂ x ρ 1 v 1 = 0, ∂ t ρ 2 + ∂ x ρ 2 v 2 = 0, ∂ t ρ 1 v 1 -v 2 1 ∂ x ρ 1 + 2v 1 ∂ x ρ 1 v 1 = 0, ∂ t ρ 2 v 2 -v 2 2 ∂ x ρ 2 + 2v 2 ∂ x ρ 2 v 2 = 0, which is equivalent to        ∂ t ρ 1 + ∂ x ρ 1 v 1 = 0, ∂ t ρ 1 v 1 + ∂ x ρ 1 v 2 1 = 0, ∂ t ρ 2 + ∂ x ρ 2 v 2 = 0, ∂ t ρ 2 v 2 + ∂ x ρ 2 v 2 2 = 0. (3.6)
Besides, at the frontier of the moment space, we obtain the pressureless gas dynamics on a single quadrature node. We then observe that for smooth solutions, thanks to the remark at the end of the previous section, the system (2.3) is nothing but either two decoupled or one single pressureless gas dynamics system of equations. We then observe that still with ε = 0, D = 0, by lemma 3.2, in both cases :

∂ t η + ∂ x q = 0, (3.7)
that is for smooth solution both inside and at the frontier of the moment space.

Let us go back to the case ε > 0. We thus have the following equality,

D = ε∇ U ηA -1 ∂ x (A∂ x U),
from which it is natural to isolate ε∂ xx η :

D = ε∇ U ηA -1 ∂ x (A∂ x U) = ε∇ U η∂ xx U + ε∇ U ηA -1 ∂ x A∂ x U = ε∂ x (∇ U η∂ x U) -ε∂ x (∇ U η)∂ x U + ε∇ U ηA -1 ∂ x A∂ x U = ε∂ xx η + ε ∇ U ηA -1 ∂ x A∂ x U -∂ x (∇ U η)∂ x U .
By lemma 3.2 giving ∇ U η, we easily get

∂ x (∇ U η)∂ x U = ρ 1 (∂ x v 1 ) 2 S ′′ (v 1 ) + ρ 2 (∂ x v 2 ) 2 S ′′ (v 2 ).
It is now a matter to calculate

∇ U ηA -1 ∂ x A∂ x U = A -t (∇ U η) t ∂ x A∂ x U.
We first observe that

∂ x A∂ x U =     0 0 2ρ 1 (∂ x v 1 ) 2 + 2ρ 2 (∂ x v 2 ) 2 6ρ 1 v 1 (∂ x v 1 ) 2 + 6ρ 2 v 2 (∂ x v 2 ) 2    
so that only the last two components of A -t (∇ U η) t are actually needed. Finally, easy calculations not reported here lead to

A -t (∇ U η) t ∂ x A∂ x U = - 1 (v 1 -v 2 ) 4 (2ρ 1 (∂ x v 1 ) 2 X 1 + 2ρ 2 (∂ x v 2 ) 2 X 2 ),
where we have set

   X 1 = (v 1 -v 2 ) 2 3 S(v 1 ) -S(v 2 ) -2(v 1 -v 2 )S ′ (v 1 ) -(v 1 -v 2 )S ′ (v 2 ) , X 2 = -(v 1 -v 2 ) 2 3 S(v 1 ) -S(v 2 ) -(v 1 -v 2 )S ′ (v 1 ) -2(v 1 -v 2 )S ′ (v 2 ) .
The entropy inequality (3.4) is then valid if and only if

∇ U ηA -1 ∂ x A∂ x U -∂ x (∇ U η)∂ x U ≤ 0, that is, setting S i = S(v i ), S ′ i = S ′ (v i ) and S ′′ i = S ′′ (v i ), i = 1,2, (v 1 -v 2 ) 2 ρ 1 (∂ x v 1 ) 2 6 S 1 -S 2 -4(v 1 -v 2 )S ′ 1 -2(v 1 -v 2 )S ′ 2 + (v 1 -v 2 ) 2 S ′′ 1 + (v 1 -v 2 ) 2 ρ 2 (∂ x v 2 ) 2 -6 S 1 -S 2 + 2(v 1 -v 2 )S ′ 1 + 4(v 1 -v 2 )S ′ 2 + (v 1 -v 2 ) 2 S ′′ 2 ≥ 0 . (3.8) A sufficient condition is given by 6 S 1 -S 2 -4(v 1 -v 2 )S ′ 1 -2(v 1 -v 2 )S ′ 2 + (v 1 -v 2 ) 2 S ′′ 1 ≥ 0, -6 S 1 -S 2 + 2(v 1 -v 2 )S ′ 1 + 4(v 1 -v 2 )S ′ 2 + (v 1 -v 2 ) 2 S ′′ 2 ≥ 0.
(3.9)

Let us focus for instance on the first inequality (the second one is treated in a similar way), and let us consider the left-hand side as a function of v 2 , for any given v 1 :

F 1 (v 2 ) = 6(S 1 -S 2 ) -4(v 1 -v 2 )S ′ 1 -2(v 1 -v 2 )S ′ 2 + (v 1 -v 2 ) 2 S ′′ 1 .

Differentiation yields

F ′ 1 (v 2 ) = 4(S ′ 1 -S ′ 2 ) -2(v 1 -v 2 )(S ′′ 1 + S ′′ 2 ), F ′′ 1 (v 2 ) = 2(S ′′ 1 -S ′′ 2 ) -2(v 1 -v 2 )S ′′′ 2 , F ′′′ 1 (v 2 ) = 2(v 2 -v 1 )S ′′′′ 2 . It is then clear that F 1 (v 1 ) = F ′ 1 (v 1 ) = F ′′ 1 (v 1 ) = F ′′′ 1 (v 1 ) = 0. Then, provided that S ′′′′ (
v) ≥ 0, ∀v, we easily get by a chain argument based on the sign of the derivative and the monotonicity property that F 1 (v 2 ) ≥ 0, ∀v 1 ,v 2 . We have thus proved the following proposition : Proposition 3.3. Smooth solutions of (3.2) satisfy the entropy inequality (3.4) for any entropy entropy-flux pair (η,q) defined by (3.5) provided that v → S(v) is a smooth function from R to R with nonnegative fourth-order derivative. In particular, the natural choice S(v) = v 2α with α ≥ 2 is suitable.

Remark 3.1. Any third-order polynomial may of course be added to the leading term of S, without changing the sign of the fourth-order derivative. However, if we focus on (strictly) convex functions v → S(v) in order to get a (strictly) convex entropy η = η(U), only first-order polynomials may be added without changing the convexity property. Remark 3.2. If we consider S(v) = 1,v,v 2 ,v 3 , it is easily checked that (3.9) holds true with two equalities. In agreement with (3.1), these choices that lead to the pairs (η,q) = (M i ,M i+1 ), i = 0,...,3, are admissible. Remark 3.3. In the case M 0 > 0, e = 0, q = 0 with the additional conditions associated with the connection between the interior and the frontier of the moment space presented in subsection 2.3, we clearly have

M k = M k 1 /M k-1 0 .
In this case, the entropy pairs clearly work also in such a case, for smooth solutions, and admits a smooth behavior at the frontier of the moment space in the cone we have defined previously.

Kinetic-macroscopic relation for smooth solutions

For smooth solutions, we established in the previous section that the four-moment model (2.2) is equivalent to the following two decoupled pressureless gas dynamics model

       ∂ t ρ 1 + ∂ x ρ 1 v 1 = 0, ∂ t ρ 1 v 1 + ∂ x ρ 1 v 2 1 = 0, ∂ t ρ 2 + ∂ x ρ 2 v 2 = 0, ∂ t ρ 2 v 2 + ∂ x ρ 2 v 2 2 = 0, (4.1)
where ρ 1 , ρ 2 , v 1 and v 2 are defined by the nonlinear system (2.6). The aim of this section is to prove a rigorous equivalence result, still for smooth solutions, between this macroscopic model and the free transport kinetic formulation (2.1) when the velocity distribution is given by a set of two Dirac delta functions. This result is nothing but a generalization of the one given in [Bouchut (1994)] for the usual pressureless gas dynamics model.

Proposition 4.1. Let T > 0 and ρ i (t,x), v i (t,x) in C 1 (]0,T [×R) for i = 1,2. Let us define f (t,x,v) = 2 i=1 ρ i (t,x)δ(v -v i (t,x)).
Then, ρ i and v i solve (2.2), or equivalently (4.1), in ]0,T [×R if and only if

∂ t f + v∂ x f = 0, in ]0,T [×R × R (4.2) in the distributional sense, i.e. if and only if ∀φ ∈ C ∞ c (]0,T [×R) and χ ∈ C ∞ c (R) T 0 R 2 i=1 ρ i (t,x) ∂ t φ(t,x) + v i (t,x)∂ x φ(t,x) χ v i (x,t) = 0. (4.3)
Proof. Let us first assume that (4.3) holds true. Since the velocity functions v i are in particular locally bounded, one can successively choose χ ∈ C ∞ c (R) such that χ(v) = v k , k = 0,...,3 for all v = v i (t,x) and then get

T 0 R 2 i=1 ρ i (t,x) ∂ t φ(t,x) + v i (t,x)∂ x φ(t,x) v i (t,x) k = 0 for all φ ∈ C ∞ c (]0,T [×R).
Invoking the closure relation (2.6), this gives the fourmoment model (2.2), or equivalently (4.1), as ρ i (t,x) and v i (t,x) are smooth functions. Conversely, let us assume that the partial differential equations of (4.1) are satisfied. Using the mass conservation equations, it is then usual to show that for i = 1,2

ρ i (∂ t v i + v i ∂ x v i ) = 0,
and then multiplying by χ ′ for any smooth function χ,

∂ t ρ i χ(v i ) + ∂ x ρ i χ(v i )v i = 0.
Summing over i = 1,2 and integrating past a test function φ ∈ C ∞ c (]0,T [×R) gives the expected result (4.3). This concludes the proof.

This Proposition allows, as a corollary, to introduce a particular type of solution which will be denoted piecewise free boundary C 1 solutions. Such solutions correspond to a discontinuous connection from the interior of the moment space to the frontier through a contact discontinuity for which the Rankine Hugoniot solutions are trivially satisfied as well as the entropy conservation equation (3.7). Corollary 4.2. We consider the following distribution at the kinetic level

f (t,x,v) = 2 i=1 ρ i (t,x)δ(v -v i (t,x))
, where ρ 1 (t,x) > 0 and v 1 (t,x) are taken as constants (or sufficiently smooth in some time interval), whereas ρ 2 (t,x) is zero except in a compact connected subset K 0 at time t = 0 of R, where ρ 2 (0,x) > 0 and v 2 (0,x) are two constants (or sufficiently smooth in some time interval) such that v 2 = v 1 . The resulting solution at the moment level exhibits two discontinuities at the frontier of the compact set K t which is the translation of set K 0 at velocity v 2 ; however the system 4.1 as well as the entropy conservation equation (3.7) are satisfied in the weak sense, that is the equations are satisfied in the usual sense where the solution is smooth and Rankine-Hugoniot conditions are satified at discontinuity points.

Proof. Clearly, the moment solution will satisfy the system of conservation equations everywhere except at the frontier of the K t set. The Rankine-Hugoniot jump conditions are trivially satisfied at the discontinuities where the mass flux associated to the first abscissa is ρ 1 (v 1 -v 2 ) in the referential of the discontinuity and leads to zero jump conditions for the part of the flux associated to the first abscissa by continuity, whereas the mass flux associated to the second abscissa is zero, as usual in contact discontinuities, which also allows to conclude. The same path allows to conclude that for any entropy-flux pair, the conservation equation is satisfied in the weak sense.

Let us underline the fact that in the region where the second weight is zero outside the compact set K t , we have used so far the convention that in such a region where a single quadrature node is to be found, the two weights are equal and the two abscissas are equal. The results proposed in the previous corollary are of course independent of such a choice since the point at the frontier of the moment space is isolated. Besides, such a corollary can be extended to as many quadrature nodes as needed as long as the number of nodes allows to describe the dynamics at the kinetic level at any point and time. Finally, the collision of two particle packets presented in subsection 6.1 satisfies the assumptions of Corollary 4.2 and will be an entropic solution.

Riemann problems and entropic measure solutions

In this section, we focus on the Riemann problem, which is associated with the inital condition for two constant states M L and M R in Ω.

M(x,0) = M L if x < 0, M R if x > 0, (5.1)
The solution of (2.3)-(5.1) is sought in the form

M(x,t) =            M L if x < σ L t, M δ L (t)δ(x -σ L t) if x = σ L t, M ⋆ if σ L t < x < σ R t, M δ R (t)δ(x -σ R t) if x = σ R t, M R if x > σ R t,
(5.2) which corresponds to the juxtaposition of two Dirac delta functions with mass M δ β and position x -σ β t, β = L,R, and separated by a constant state M ⋆ in Ω. Here, σ β denotes a real number, we have m β (t) ≥ 0, m β (0) = 0 and β = L,R, and M δ β (t) is defined by

M δ β (t) =     m β (t) m β (t)σ β m β (t)σ 2 β m β (t)σ 3 β     (5.3)
We introduce the following natural definitions of (entropy) measure solutions.

Definition 5.1. Let η = ρ 1 S(v 1 ) + ρ 2 S(v 2 ) and q = ρ 1 v 1 S(v 1 ) + ρ 2 v 2 S(v 2 ) with S(v) = v 2α , α ∈ N. We say that (5.2) is a measure solution of (2.3)-(5.1) if and only if ∂ t η + ∂ x q = 0 in D ′ (]0,∞[×R) for S(v) = v 2α with α = 0, 1 2 ,1, 3 2 , that is for (η,q) = (M i ,M i+1 ), i = 0,1,2,3. We say that (5.2) is an entropic measure solution of (2.3)-(5.1) if and only if

∂ t η + ∂ x q ≤ 0 in D ′ (]0,∞[×R) for S(v) = v 2α with α = 0, 1 2 ,1, 3 2 and α ≥ 2.
We can prove the following equivalence result.

Theorem 5.2. The solution given by (5.2) is a measure solution of (2.3)-(5.1) iff

   σ L (M L -M ⋆ )t -F(M L ) -F(M ⋆ ) t + M δ L (t) = 0 σ R (M ⋆ -M R )t -F(M ⋆ ) -F(M R ) t + M δ R (t) = 0.
(5.4)

The solution given by (5.2) is a entropic measure solution of (2.3)-(5.1) if and only if in addition, for all S(v) = v 2α with α ≥ 2 :

   σ L η(M L ) -η(M ⋆ ) t -q(M L ) -q(M ⋆ ) t + η(M δ L (t)) ≤ 0 σ R η(M ⋆ ) -η(M R ) t -q(M ⋆ ) -q(M R ) t + η(M δ R (t)) ≤ 0,
(5.5)

Proof. We first recall that ∂ t η + ∂ x q = 0 in D ′ (]0,∞[×R) means that for any smooth function with compact support ϕ ∈ C ∞ c (]0,∞[×R), we have

< ∂ t η + ∂ x q,ϕ >= -< η,∂ t ϕ > -< q,∂ x ϕ >= - ∞ 0 +∞ -∞ η(x,t)∂ t ϕ(x,t)dxdt - ∞ 0 +∞ -∞ q(x,t)∂ x ϕ(x,t)dxdt = 0,
where by definition (5.2)

η(x,t) =            η(M L ), x < σ L t, η(M δ L (t))δ(x -σ L t), x = σ L t, η(M ⋆ ), σ L < x t < σ R , η(M δ R (t))δ(x -σ R t), x = σ R t, η(M R ), x > σ R t, q(x,t) =            q(M L ), x < σ L t, q(M δ L (t))δ(x -σ L t), x = σ L t, q(M ⋆ ), σ L < x t < σ R , q(M δ R (t))δ(x -σ R t), x = σ R t, q(M R ),
x > σ R t.

Here we note that q(M δ β (t

)) = σ β η(M δ β (t)) with η(M δ β (t)) = m β (t)S(σ β ), β = L,R. For all ϕ ∈ C ∞ c (]0,∞[×R) we thus have < ∂ t η + ∂ x q,ϕ >= - ∞ 0 dt σLt -∞ η(M L )∂ t ϕ(x,t)dx - ∞ 0 η(M δ L (t))∂ t ϕ(σ L t,t)dt - ∞ 0 dt σRt σLt η(M ⋆ )∂ t ϕ(x,t)dx - ∞ 0 η(M δ R (t))∂ t ϕ(σ R t,t)dt - ∞ 0 dt ∞ σRt η(M R )∂ t ϕ(x,t)dx - ∞ 0 dt σLt -∞ q(M L )∂ x ϕ(x,t)dx - ∞ 0 q(M δ L (t))∂ x ϕ(σ L t,t)dt - ∞ 0 dt σRt σLt q(M ⋆ )∂ x ϕ(x,t)dx - ∞ 0 q(M δ R (t))∂ x ϕ(σ R t,t)dt - ∞ 0 dt ∞ σRt q(M R )∂ x ϕ(x,t)dx that is, using in particular q(M δ β (t)) = σ β η(M δ β (t)), < ∂ t η + ∂ x q,ϕ >= = - ∞ 0 dtη(M L ) d dt σLt -∞ ϕ(x,t)dx -σ L ϕ(t,σ L t) - ∞ 0 q(M L )ϕ(σ L t,t)dt - ∞ 0 dtη(M ⋆ ) d dt σRt σLt ϕ(x,t)dx -σ R ϕ(t,σ R t) + σ L ϕ(t,σ L t) - ∞ 0 q(M ⋆ ) ϕ(σ R t,t) -ϕ(σ L t,t) dt - ∞ 0 dtη(M R ) d dt ∞ σRt ϕ(x,t)dx + σ R ϕ(t,σ R t) + ∞ 0 q(M R )ϕ(σ R t,t)dt - ∞ 0 η(M δ L (t)) d dt [ϕ(σ L t,t)]dt - ∞ 0 η(M δ R (t)) d dt [ϕ(σ R t,t)]dt = ∞ 0 ϕ(σ L t,t) σ L η(M L ) -q(M L ) dt + ∞ 0 ϕ(σ R t,t) σ R η(M ⋆ ) -q(M ⋆ ) dt - ∞ 0 ϕ(σ L t,t) σ L η(M ⋆ ) -q(M ⋆ ) dt - ∞ 0 ϕ(σ R t,t) σ R η(M R ) -q(M R ) dt - ∞ 0 η(M δ L (t)) d dt [ϕ(σ L t,t)]dt - ∞ 0 η(M δ R (t)) d dt [ϕ(σ R t,t)]dt = ∞ 0 ϕ(σ L t,t) σ L (η(M L ) -η(M ⋆ )) -(q(M L ) -q(M ⋆ )) dt + ∞ 0 ϕ(σ R t,t) σ R (η(M ⋆ ) -η(M R )) -(q(M ⋆ ) -q(M R )) dt - ∞ 0 η(M δ L (t)) d dt [ϕ(σ L t,t)]dt - ∞ 0 η(M δ R (t)) d dt [ϕ(σ R t,t)]dt = ∞ 0 ϕ(σ L t,t) d dt σ L (η(M L ) -η(M ⋆ ))t -(q(M L ) -q(M ⋆ ))t + η(M δ L (t)) dt + ∞ 0 ϕ(σ R t,t) d dt σ R (η(M ⋆ ) -η(M R ))t -(q(M ⋆ ) -q(M R ))t + η(M δ R (t)) dt.
Then, since m β (0) = 0 and then η(M δ β (0)) = 0, β = L,R, it is clear that (5.2) is a measure solution of (2.3)-(5.1) if and only if (5.4) is valid for all t ≥ 0, and an entropic measure solution if and only if in addition (5.5) holds true for all t ≥ 0 and all η = ρ 1 S(v 1 ) + ρ 2 S(v 2 ) and q

= ρ 1 v 1 S(v 1 ) + ρ 2 v 2 S(v 2 ) with S(v) = v 2α , α ≥ 2. Remark 5.1. Let us recall that η(M δ β (t)) = m β (t)S(σ β ). Since (5.5
) is made of equalities when S(v) = 1 and S(v) = v (we get in these cases the first two components of (5.4)), the validity of (5.5) for all S(v) = v 2α , α ≥ 2, is equivalent to the validity of

   σ L η(M L ) -η(M ⋆ ) -q(M L ) -q(M ⋆ ) ≤ 0 σ R η(M ⋆ ) -η(M R ) -q(M ⋆ ) -q(M R ) ≤ 0, (5.6) for all S(v) = v 2α -σ 2α L v-σR σL-σR -σ 2α R v-σL σR-σL , α ≥ 2.

Examples of entropic solutions

In this section, we propose three particular entropic solutions. The first one models the collision of two particles packets with free boundary C 1 smooth solution, that is a solution for which an exact link with the kinetic level is preserved and for which the entropy equation is exactly satisfied. In such a situation the fourmoment model does not develop δ-shock Dirac delta functions and is actually able to properly represent the crossing of the two packets which correspond to the dynamics at the kinetic level. The second one models the collision of four particles packets. In this case and as expected since the number of moments is set to four, the entropic solutions involves two δ-shock Dirac delta functions singularities. Whereas the first case corresponds to a connection from the interior of the moment space to the frontier through a contact discontinuity, or free boundary solution, resulting in an isolated point at the frontier, we consider in a third example a smooth connection to the frontier of the moment space, such that the point at the frontier is an accumulation point of a trajectory inside the moment space.

Collision of two particles packets

We consider a Riemann initial data (5.1) where

M L = M(U L ) and M R = M(U R ) are such that U L = 1 2     ρ L ρ L ρ L v L ρ L v L     and U R = 1 2     ρ R ρ R ρ R v R ρ R v R    
for two given densities ρ L > 0 and ρ R > 0 and velocities v L > 0 and v R < 0. We recall that the function M = M(U) is defined by (2.6). We define

M(x,t) =    M L if x < v R t, M ⋆ if v R t < x < v L t, M R if x > v L t, (6.1) with M ⋆ = M(U ⋆ ) given by U ⋆ = (ρ L ,ρ R ,ρ L v L ,ρ R v R ) t .
Our objective here is to prove that the following Dirac delta functions free solution is an entropy solution of (2.3)-(5.1). Conditions (5.4) and (5.5) write here

v R (M L -M ⋆ ) -F(M L ) -F(M ⋆ ) = 0 v L (M ⋆ -M R ) -F(M ⋆ ) -F(M R ) = 0. (6.2) and v R η(M L ) -η(M ⋆ ) -q(M L ) -q(M ⋆ ) ≤ 0 v L η(M ⋆ ) -η(M R ) -q(M ⋆ ) -q(M R ) ≤ 0, (6.3) with η(U) = ρ 1 S(v 1 ) + ρ 2 S(v 2 ) and q(U) = ρ 1 v 1 S(v 1 ) + ρ 2 v 2 S(v 2 ) for all S(v) v 2α
with α ≥ 2. We will focus only on the first equality of (6.2) and the first inequality of (6.3), the second ones being treated similarly. We clearly have

M L -M ⋆ =     ρ L ρ L v L ρ L v 2 L ρ L v 3 L     -     ρ L + ρ R ρ L v L + ρ R v R ρ L v 2 L + ρ R v 2 R ρ L v 3 L + ρ R v 3 R     = -     ρ R ρ R v R ρ R v 2 R ρ R v 3 R     , while F(M L ) -F(M ⋆ ) =     ρ L v L ρ L v 2 L ρ L v 3 L ρ L v 4 L     -     ρ L v L + ρ R v R ρ L v 2 L + ρ R v 2 R ρ L v 3 L + ρ R v 3 R ρ L v 4 L + ρ R v 4 R     = -     ρ R v R ρ R v 2 R ρ R v 3 R ρ R v 4 R     .
It is then clear that the first equality of (6.2) holds true. Let us now check that the proposed Riemann solution fulfills the entropy condition. We clearly have

v R η(M L ) -η(M ⋆ ) -q(M L ) -q(M ⋆ ) = v R ρ L S(v L ) -(ρ L S(v L ) + ρ R S(v R )) -ρ L v L S(v L ) -(ρ L v L S(v L ) + ρ R v R S(v R ))
= 0 which allows to prove that the proposed solution is an entropic smooth solution.

Collision of four particles packets

We consider a Riemann initial data (5.1) where

M L = M(U L ) and M R = M(U R ) are such that U L = 1 2     ρ ρ ρv 1 ρv 2     and U R = 1 2     ρ ρ -ρv 2 -ρv 1    
for a given density ρ > 0 and two velocities v 2 > v 1 > 0. We have

M L = 1 2     2ρ ρ(v 1 + v 2 ) ρ(v 2 1 + v 2 2 ) ρ(v 3 1 + v 3 2 )     and M R = 1 2     2ρ -ρ(v 1 + v 2 ) ρ(v 2 1 + v 2 2 ) -ρ(v 3 1 + v 3 2 )     F(M L ) = 1 2     ρ(v 1 + v 2 ) ρ(v 2 1 + v 2 2 ) ρ(v 3 1 + v 3 2 ) ρ(v 4 1 + v 4 2 )     and F(M R ) = 1 2     -ρ(v 1 + v 2 ) ρ(v 2 1 + v 2 2 ) -ρ(v 3 1 + v 3 2 ) ρ(v 4 1 + v 4 2 )     .
We define

M(x,t) =            M L if x < -σt, M δ L (t)δ(x + σt) if x = -σt, M ⋆ if -σt < x < σt, M δ R (t)δ(x -σt) if x = σt, M R if x > σt, (6.4) with σ > 0, M ⋆ = M(U ⋆ ) given by U ⋆ = 1 2     ρ ⋆ ρ ⋆ -ρ ⋆ v ⋆ ρ ⋆ v ⋆     , M ⋆ =     ρ ⋆ 0 ρ ⋆ v 2 ⋆ 0     , ρ ⋆ > 0, v ⋆ > 0, and M δ L (t), M δ R (t) given by M δ L (t) = m(t)     1 -σ σ 2 -σ 3     , M δ R (t) = m(t)     1 σ σ 2 σ 3     , F(M ⋆ ) =     0 ρ ⋆ v 2 ⋆ 0 ρ ⋆ v 4 ⋆     ,
with m(t) ≥ 0. The generalized Rankine-Hugoniot jump conditions (5.4) write here

   -σ(M L -M ⋆ )t -F(M L ) -F(M ⋆ ) t + M δ L (t) = 0 σ(M ⋆ -M R )t -F(M ⋆ ) -F(M R ) t + M δ R (t) = 0.
that is, equivalently

                           2σM ⋆ -σ(M L + M R ) + (F(M R ) -F(M L )) + m(t) t     2 0 2σ 2 0     = 0, 2F(M ⋆ ) + σ(M R -M L ) -(F(M L ) + F(M L )) - m(t) t     0 2σ 0 2σ 3     = 0. (6.5)
This is made of eight equalities, four are trivial (zero equals zero), so that four are left to determine the four unknowns ρ ⋆ , v ⋆ , σ and m(t)/t. We propose below to numerically solve this nonlinear system for a specific set of values for ρ, v 1 and v 2 .

Remark 6.1. We conjecture existence and uniqueness of a solution to this nonlinear system. Indeed, the initial condition involves four different velocities while the model is able to represent two different velocities only. More precisely and by analogy with the usual pressureless gas dynamics model (see [START_REF] Bouchut | [END_REF]], [Bouchut (1994)]), velocity v 1 is to "bump" into velocities -v 1 and -v 2 to create a first Dirac delta function. By symmetry, another Dirac delta function is expected.

Regarding the entropy inequalities (5.5), we first remark that for

S(v) = v 2α , η(M L ) = η(M R ) = ρ 2 (S(v 1 ) + S(v 2 )) η(M ⋆ ) = ρ ⋆ S(v ⋆ ), η(M δ L (t)) = η(M δ R (t)) = m(t)S(σ), while q(M L ) = -q(M R ) = ρ 2 (v 1 S(v 1 ) + v 2 S(v 2 )), q(M ⋆ ) = 0.
As an immediate consequence, both inequalities in (5.5) are equivalent and the entropy condition is

σ ρ ⋆ S(v ⋆ ) - ρ 2 (S(v 1 ) + S(v 2 )) - ρ 2 (v 1 S(v 1 ) + v 2 S(v 2 )) + m(t) t S(σ) ≤ 0. (6.6)
A concrete example. We propose to take ρ = 1, v 1 = 0.8 and v 2 = 1.2. Numerically solving (6.5) gives ρ ⋆ = 1.88265, v ⋆ = 1.06026, σ = 0.87983 and m(t)/t = 0.22342. If the left-hand side of (6.6), which represents the entropy dissipation rate associated with S(v) = v 2α , is denoted D(α), a simple calculation gives for instance D(2) = -0.27324, D(3) = -0.86854, D(4) = -1.88678,... The proposed solution (6.4) is then actually an entropy measure solution of the four-moment model. Such an exact entropic solution will be used in the following to prove the relevance of the numerical scheme proposed hereafter with respect to the exact solution when singularities are present.

Piecewise linear solution connected with the frontier of the moment space

As a last example, we introduce a piecewise linear solution which allows to connect zones inside the moment space with zone at the frontier within the proper framework introduced in subsection 2.3. Particles are initiated in the domain [0, 0.5]. In the domain [0, 0.1], a monomodal velocity distribution is reconstructed, with v 1 = v 2 = 1, and ρ 1 = ρ 2 = 0.5. On the contrary, a bimodal velocity distribution is reconstructed in the domain [0.1, 0.4], with ρ 1 = ρ 2 = 0.5 and for the abscissas, v 1 = 1 + x-0.1 0.3 and v 2 = 1. The initial conditions are represented in Fig. 6.1.

The ground difference with the first test case is that the transition between the two zones is smooth, and so the numerical strategy to account for this transition is important. The analytical solution of this problem, in smooth areas, consists of a decoupled transport of each of the quadrature nodes as two independant pressureless gas as showed in system (3.6). This comes from the fact that the number of Dirac delta functions reconstructed from the moments, two in this case, is always sufficient to capture the problem dynamics. The equivalence between the kinetic and macroscopic equations is preserved, and the solution in terms of moments satisfies the entropy equation. Therefore, the solution is the superposition of a translation at constant velocity v 2 = 1, and a transport with the following velocity field:

v 1 =    1 x ∈ [0, 0.1] 1 + x-0.1 0.3 x ∈ [0.1, 0.4] 2 x ∈ [0.4, 0 .5] (6.7)
The analytical solution, displayed in Fig. 6.2, has two fronts at x = 0.2 and x = 0.9 moving with a velocity v = 1 and v = 2 respectively. The first front corresponds to particles initiated with velocity v = 1 between x = 0 and x = 0.1. The second front corresponds to particles initiated with velocity v = 2 between x = 0.4 and x = 0.5. The square wave between x = 0.8 and x = 0.9 is the final location of the particles initiated with velocity v = 2 between x = 0.4 and x = 0.5. The value of ρ 1 between x = 0.3 and x = 0.8 corresponds to the expansion of the density field due to transport with a linear velocity field with positive slope. The value of the density in that area is the solution of the equation

∂ t ρ 1 + v 1 ∂ x ρ 1 = -ρ 1 ∂ x v 1 , which yields ρ 1 = 0.3 at time t = 0.2 4 .
According to the initial conditions, both weights have the same profile, the quantities q/(M 0 e) and q/e 3/2 are null. At time t = 0.2, the weight profiles are different in the interval [0.3,0.8] corresponding to the expansion of ρ 1 . Therefore, q/e 3/2 is non null, as well as q M0e , since the velocities have also different values, and can be exactly calculated (see Fig. 7.10).

Numerical simulations via kinetic schemes

This section is devoted to the discretization of (2.3)-(2.5)-(2.6). As already stated, we use as a building block a natural first-order kinetic scheme already proposed in the literature [Jin & Li (2003), Gosse et al. (2003), [START_REF] Desjardins | A quadrature-based moment method for dilute fluid-particle flows[END_REF]] and briefly recalled here for the sake of completeness. Let us first introduce a time step ∆t > 0 and a space step ∆x > 0 that we assume to be constant for simplicity. We set λ = ∆t/∆x and define the mesh interfaces x j+1/2 = j∆x for j ∈ Z, and the intermediate times t n = n∆t for n ∈ N. In the sequel, M n j denotes the approximate value of M at time t n and on the cell C j = [x j-1/2 ,x j+1/2 ). For n = 0, we set M 0 j = 1 ∆x x j+1/2

x j-1/2 M 0 (x)dx, j ∈ Z, where M 0 (x) is the initial condition. Let us now assume as given (M n j ) j∈Z the sequence in Ω of approximate values at time t n . In order to advance it to the next time level t n+1 , the kinetic scheme is decomposed into two steps.

First step : transport (t n → t n+1-) We first set U n j = U(M n j ) and define the function (x,v) → f n (x,v) by

f n (x,v) = (ρ 1 ) n j δ v -(v 1 ) n j + (ρ 2 ) n j δ v -(v 2 ) n j , ∀ (x,v) ∈ C j × R, j ∈ Z.
We then solve the transport equation

∂ t f + v∂ x f = 0, (x,v) ∈ R × R, f (t = 0,x,v) = f n (x,v), (7.1)
the solution of which is given by f

(t,x,v) = f n (x -vt,v). At last, we set f n+1-(x,v) = f n (x -v∆t,v).
Second step : collapse (t n+1-→ t n+1 ) The first four moments at time t n+1 are now naturally defined by setting

(M i ) n+1 j = 1 ∆x x j+1/2 x j-1/2 +∞ -∞ v i f n+1-(x,v)dvdx.
Then, we have M n+1 j = (M 0 ) n+1 j ,(M 1 ) n+1 j ,(M 2 ) n+1 j ,(M 3 ) n+1 j t for all j ∈ Z, which completes the algorithm description.

Remark 7.1. It is easy to see that this scheme preserves the moment space Ω, see for instance [START_REF] Desjardins | A quadrature-based moment method for dilute fluid-particle flows[END_REF]].

Remark 7.2. Under the natural CFL condition ∆t max j∈Z ((v 1 ) n j ,(v 2 ) n j ) ≤ CF L∆x, with CF L ≤ 1, integrating (7.1) over (t,x,v) ∈ (0,∆t) × C j × R and against v i , i = 0,...,3 easily leads to the equivalent update formula

M n+1 j = M n j - ∆t ∆x F n j+1/2 -F n j-1/2 , j ∈ Z,
where we have set

F n j+1/2 = (M 1 ) n j+1/2 ,(M 2 ) n j+1/2 ,(M 3 ) n j+1/2 ,(M 4 ) n j+1/2
t and

(M i ) n j+1/2 = (M i ) n+ j+1/2 + (M i ) n- j+1/2
, and with

(M i ) n- j+1/2 = (ρ 1 ) n j+1 min(0,(v 1 ) n j+1 ) (v 1 ) n j+1 i-1 + min(0,(v 2 ) n j+1 )(ρ 2 ) n j+1 (v 2 ) n j+1 i-1 , (M i ) n+ j+1/2 = (ρ 1 ) n j max(0,(v 1 ) n j+1 ) (v 1 ) n j i-1 + (ρ 2 ) n j max(0,(v 2 ) n j+1 ) (v 2 ) n j i-1 .
7.1. Numerical quadrature strategy at the frontier of the moment space

This paragraph addresses the issue of how to numerically handle the transition between a vector in the interior of the moment space, and a vector lying at its frontier. For an isolated point at the frontier of the moment space Γ, there is no specific problem since we use a single quadrature node and the quadrature is not a problem. The two problems we have to face are related to the preservation of the cone in which we envision to work in section 2 as well as to deal with finite precision algebra in the neighborhood of the point (0,0) in the (e,q) plane.

Consequently we introduce two constants in the numerical quadrature we use. First, for finite precision algebra and in order to avoid numerical errors, we only evaluate the two quadrature nodes when e/M 2 0 > ǫ 1 , where ǫ 1 is a small number related to machine precision. Under this threshold, the velocity dispersion is considered null (e = 0), and the set of Dirac delta functions are reconstructed as suggested in section 2:

ρ 1 = ρ 2 = M 0 /2, v 1 = v 2 = M 1 /M 0 .
Secondly, since we want to deal with compactly supported velocity distribution at the kinetic level, we will introduce another constant, η which is a bound for the distance between the two abscissas. In fact we would like to impose to the solution to remain inside the cone in the (e,q) plane : |q| M0e ≤ η. It will be shown in the following examples that the cone is in fact automatically preserved by the proposed algorithm and that such a bound does not have to be imposed but is satisfied by the numerical solution.

If |q| M0e > η, then we set q/(M 0 e) = sign(q)η, so that : q M 0 e = sign(q)η, q = sign(q)η M 0 e, M 3 = q + M 1 M 2 /M 0 .

As illustrated in Fig 7 .1 and explained in Section 2.3, the variable q naturally lives in a cone delimited by the straight lines of slope ±η and becomes null when e ≤ ǫ 1 .

ε1 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 11111111111111111 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 00000000000000000 Let us remark that limiting q M0e does not lead to a limitation on the quantity q e 3/2 in such a way as to allow the possibility to reach large ratio in terms of weights (one weight can approach a zero value while the quantity q e 3/2 grows indefinitely) but without allowing the abscissas distance to grow beyond of fixed limit naturally inherited from the initial solutions at the kinetic level. We will come back to this point in the following result section. For the simulations we present, the parameters are such that: ǫ 1 = 10 -9 and η = 2.

Numerical results

This section is devoted to numerical illustrations of the two Riemann problems and one dedicated to the case of a free boundary connecting zone where e = 0 and a zone where e > 0, discussed in Section 6.

In all the figures, we choose to represent ρ 1 and v 1 by solid lines, and ρ 2 and v 2 by lines with circle markers. In the representation of weights and abscissas, values have to be assigned for v 1 and v 2 : we thus decide that v 1 is the maximum of the relative values of velocity.

Two packet collision. Figure 7.2 represents the initial conditions for the two particle packet case. Figures 7.3 and 7.4 presents the numerical and analytical solutions for the two particle packet case with ρ L = ρ R = 1, and v L = 1, v R = -1. The computation is run with a 1000 cell grid on the spatial domain [0,1], with CF L = 1. The length of each packet is 0.4 (ρ 1 = ρ 2 = v 1 = v 2 = 0 for x ≤ 0.1 and x ≥ 0.9) and the two packets start to collide exactly at time t = 0. Moving in opposite direction one across the other, with the same opposite speed, they then overlay and we note in particular that ρ 1 = ρ 2 = 1 and v 1 = v 2 = 0 in the mixing zone (see for instance the plots at time t = 0.1). As expected, they finally get separated again and we note that a perfect agreement is obtained with the exact entropic solution.

Four packet collision. Figures 7.5 presents the initial conditions. Figures 7.6 and 7.8 present the numerical and analytical solutions respectively for the moments and the weights and abscissas. The computation is run with a 1000 cell grid on the spatial domain [0,1], at CF L = 1. Here, we observe the presence of two Dirac delta functions as already discussed in Section 6. The agreement between exact and numerical solutions for the moments is very good, showing that the numerical solution converges to the analytical one. The disparities encountered between the analytical and numerical solution in the case of the weights and abscissas are due to the fact that the mapping U(M) is discontinuous at the moment space border. In area of numerical diffusion (x ≈ 0.22 and x ≈ 0.88), when the dispersion e is under the threshold ǫ 1 , weights and abscissas are reconstructed as explained in section 2.

The wave propagating at velocity 1.2 and separating the constant states (v 1 = 0.8,v 2 = 0) and (v 1 = 1.2, v 2 = 0.8) is steep and coincide with the analytical wave whereas the wave propagating at velocity 0.8 is actually smooth since the CFL number is based on the highest value of velocity, which is 1.2 in this studied case. The same explanation holds for the symmetric jump at location x = 0.78. Meanwhile, because of the conservation of the velocity moments, the numerical velocity jump (at x = 0.154) happens before the analytical velocity jump (at x = 0.18). One can here notice that the quadrature method provides the expected value of velocity in the numerical diffusion zones. The same explanation holds for the different velocity jump locations between the analytical and numerical solution at x = 0.82 and x = 0.845. The same phenomenon is responsible for the disparities between the analytical and numerical solutions at the δ-shocks locations, at x = 0.4 and x = 0.6.

Figure 7.7 displays the final profile of the quantities q M0e and q e 3/2 . Thus, q

M0e

has significative values in areas where the abscissa distance as well as the weight difference are important, whereas q e 3/2 reaches high value in areas where the weight ratio is important. Since in the domain, except for the singularities, the weights have the same value, both quantities are equal to zero. At the singularities, q e 3/2 is roughly proportional to the square root of the weight ratio, and q M0e is bounded, accounting for the fact that the velocity field is bounded.

We have thus provided numerical simulations in the two cases for which we have at our disposal an analytical entropic solution, either in the piecewise contant case, or in the singular case where δ-shock mesure solutions are present. In the first case, the crossing of the two droplets monokinetic packets is very properly reproduced without numerical diffusion since we work at CFL one, even if this infinite Knudsen number limit, where the various droplet packets cross without interacting, differs from the entropic solution of the system of partial differential equations (2.2) obtained through the quadrature-based closure. This is the same type of behavior as seen in the case of pressureless gas dynamics at a lower level.

Free boundary case. The last test case, explained in subsection 6.3, assesses the ability of the method to solve free boundary cases connecting in a continuous manner states lying in the interior of the moment space and at its frontier. The chosen grid contains 400 cells, CFL number is set to 0.98. This value of the CFL number is chosen is order to prevent high frequency instabilities to occur. The computation is run until t = 0.2. The analytical solution has been provided in subsection 6.3.

Results are displayed in Fig. 7.9. Let us first focus on the fronts present at of PTC where the solution remains smooth and where the scheme allows to describe the phase space dynamics properly as well as for cases where the complexity of the dynamics in phase space leads to generalized δ-shocks, as observed for pressureles gas dynamics due to the weakly hyperbolic structure of the system of conservation laws. Two stumbling blocks still remains to be treated. First, we would need a uniqueness theory and a convergence analysis in a general framework in order to fully justify the use of the kinetic schemes for the simulation of such models. However, as explained already in [Bouchut (1994)], the framework of entropic solution is not sufficient in order to provide uniqueness since one can exhibit multiple entropic solutions for measure solutions. Let us underline that it is easy to construct the same type of measure solutions for system (2.3), which is the exact same collision case used by Bouchut, but with a motionless Dirac delta function in density localized at the collision point of the other two incoming "particles". An infinite set of entropic solution can then be exhibited depending on the nature of the collision. As a consequence, it would be first useful to investigate such a point on the pressureless gas dynamics and then to extend it to the present system of higher order quadrature-based moment models.

Besides, the construction of fully high order methods is still an open question and requires further developments. At last, let us mention that most of the results of the present paper do naturally extend to higher order moment systems, but at the price of algebra complications. In fact, the key point lies in the extension of the proposed study of the behavior of the quadrature at the frontier of the moment space, namely when the two velocities v 1 and v 2 become equal. If we consider for instance the 6-moment model and assume that one of the three velocities v 1 , v 2 and v 3 is smooth while the other two become equal, we are in the same framework as in the present paper. But the case when the three velocities tend to be equal needs to be precised in a future work.

  Fig. 6.1. Moment dynamics for a free boundary connecting areas where e = 0 and e > 0. Initial conditions. Top-left: M 0 . Top-right: M 1 , Bottom-left: weights. Bottom-right: abscissas. The solid line corresponds to (ρ 1 ,v 1 ), the dashed line with circles to (ρ 2 ,v 2 ).

  Fig. 6.2. Moment dynamics for a free boundary connecting areas where e = 0 and e > 0. Solution at t = 0.2. Top-left: M 0 . Top-right: M 1 . Bottom-left: weights, Bottom-right:abscissas. The solid line corresponds to (ρ 1 ,v 1 ), the dashed line with circles to (ρ 2 ,v 2 ).

  Fig. 7.1. Handling of the moment space border with the admissible cone.

  Fig. 7.2. Initial fields of weights (left) and abscissas (right) for the two particle packet case. The first quadrature node is represented by solid lines whereas the second node is represented by circles.

  Fig. 7.7. Four particle packet case, at time t = 0.1. Left: ratio q/(M 0 e). Right: ratio q/e 3/2 .

The corresponding values of e and q are e = αγ 2 x β+2δ and q = αγ

x β+3δ (1 -αx β ).

The mapping of M onto U will never be smooth at point (0,0) in the (e,q) plane since the limit of U will depend on the limit of q/(M 0 e), whereas the limit value of M in the proposed cone is always fixed.

In the interval [0.8,0.9], although ρ 2 should be null (the square wave at velocity v 2 = 1 should be bounded between the front x = 0.2 and x = 0.7) we computed ρ 2 = 0.5 and v 2 = 2. This is consistent with the conditions in Section 2.3 for moment vectors at the frontier of the moment space.

Acknowledgement. We would like to thank the referees for a very careful reading of the paper and for useful hints to improve and make complete the present work.

This research was partially supported by an ANR Young Investigator Award (French Agence Nationale de la Recherche), contract ANR-08-JCJC-0132-01 -INTOCS -2009-2013 This research was supported by an ANR Young Inverstigator Award for M. Massot, ANR-05-JCJC-0013 -jéDYS -2005-2009. We would like to thank the Center for Turbulence Research at Stanford University and its director, Parviz Moin, for their kind invitation during the Summer Program 2010 where the present work was completed and Laboratory EM2C for supporting the visit of C. Chalons at Stanford as well as Damien Kah's Master Thesis during which the study was initiated.

is not symptomatic of the numerical diffusion such methods will encounter with a first order method in realistic configurations [START_REF] Kah | Eulerian quadrature-based moment models for dilute polydisperse evaporating sprays[END_REF]]. In the second case, the numerical method is able to capture the creation of the measure singular solutions associated to the fact the we have limited the number of quadrature node to two. With this node number limitation, the proper physical solution, in the x = 0.2 and x = 0.9 for the analytical solution. Since the CFL number is based on the highest velocity value (2 in this case) and is taken as 0.98, the corresponding wave is less diffused at x = 0.9, contrary to the front wave located at x = 0.2 moving at velocity v 2 = 1. Note that in these areas, ρ 1 = ρ 2 , since e = 0 or e < ǫ 1 . The borders of these areas are clearly seen at x ≈ 0.28 and x = 0.8. The constant profiles for ρ 2 and ρ 1 observed respectively between x = 0.5 and x = 0.7 and between x = 0.5 and x = 0.8 correspond to the one observed in the analytical solution. In the area between x = 0.28 and x = 0.5, the density profiles would be expected to be constant, with the same value as before. Instead of that, one observes a peak value for ρ 2 and a low value for ρ 1 , the sum ρ 1 + ρ 2 being constant. This results from a coupling between the behaviour of the quadrature method when e tends to zero and the numerical diffusion. We are here in the situation q > 0 and e small but e > ǫ 1 . See Lemma 2.3. Further away from the discontinuities the density values tend to their analytical value. M0e and q e 3/2 at time t = 0.2. First, it is interesting to note that q M0e is naturally bounded by 1 in the numerical simulation, as shown in Fig 7 .10 such that the bound imposed with η = 2 in never effective. Let us underline that in the present case, the numerical scheme allows to preserve the proposed cone associated with a maximal distance between the abscissas which is also invariant. Comparing the analytical value of |q/(M 0 e)| to the numerical resolution in Fig. 7.10, it is instructive to observe that the numerical diffusion is creating zones where it can be far from zero, whereas it should be zero in the analytical solution. However, this is done is such a way as to preserve the maximal value foreseen as the maximal distance between the abscissa at time t = 0 which is one. Let us also underline that we have rerun this case with various values of ǫ 1 varying from 10 -9 up to 10 -7 without any effect on the solution. It can also be noticed that the quantity q/e 3/2 has Results at t = 0.2. Left: ratio q/(M 0 e). The dashed line with crosses represents the analytical solution, whereas the solid line represents the numerical solution, bounded by its extremal initial values (dotted-dashed curve). Right: ratio q/e 3/2 . large values in the regions of connection between the interior and the frontier of the moment space, but has no reason to be naturally limited as opposed to the previous quantity as it can be seen in Fig. 7.9-bottom left.

Concerning the convergence behavior of the numerical solution, we display the error in L 1 norm relative to the analytical solution to assess convergence quantitatively, Tab 7.1, for 400, 800, 1600, and 3200 cell grids. As expected, we get an experimental order of convergence of 0.5. These data clearly show the convergence towards the analytical solution for each moment.

As a consequence, it can be seen that we have designed the proper theoretical setting for the transition from the interior of the moment space towards its frontier since the cone we have defined seems to be automatically preserved by the kinetic scheme we have used. Such a point would be worth a detailed study which is beyond the scope of the present paper.

Conclusion

In this paper, we have extended the notion of entropic measure solution of quadrature-based moment method for kinetic equations. Such kinetic equations are frequently encountered in many application fields where a complex dynamics in phase space is involved. Following the contribution of [START_REF] Bouchut | [END_REF]] for the pressureless gas dynamics which is the one-node quadrature version of a more general system of conservation laws for quadrature-based moment models, we have been able to provide a few problem test-cases showing that the numerical solution of the resulting system of conservation laws through kinetic schemes reproduces the defined entropic solution as well as the proper theoretical setting for the transition from the interior to the frontier of the moment space. It is an important point for the case