Benoît Bruno Després

Benoît Perthame
email: benoit.perthame@upmc.fr

Uncertainty propagation; intrusive kinetic formulations of scalar conservation laws

Keywords: Uncertainty propagation, intrusive kinetic formulations of scalar conservation laws Uncertainty propagation, kinetic formulation of conservation laws, maximum principle, entropy dissipation, chaos polynomial. Mathematics Subject Classification (2010) 35L65, 35R60, 35A35

Introduction

We address the question of constructing intrusive kinetic methods for scalar conservation laws in view of uncertainty quantification (UQ) and propagation. The starting point is a scalar conservation law in dimension d

∂ t u + ∇ • F (u) = 0, ( 1
)
where u is the unknown, x ∈ R d is the space variable, t ≥ 0 is the time variable and F : R → R d is the flux (we use at some places that its second derivatives is locally bounded). We assume the solution u is a function of an additional variable ω, called the uncertainty variable. For simplicity of notations, we take in this work ω ∈ I ⊂ R, a bounded interval. It can be generalized to ω ∈ I p ⊂ R p , p > 1 by tensorisation.

The unknown u is therefore a function of (x, ω, t), meaning that we would like to solve an infinite series of conservation laws [START_REF] Bernardi | Polynomial interpolation results in Sobolev space[END_REF], where the dependency on ω stems from the initial data. Throughout all this paper, we assume it is nonnegative, upper bounded by some U M > 0 and belongs to L 1 (dxdµ(ω)) 0 ≤ u(x, ω, 0) = u init (x, ω) ≤ U M , u init (x, ω)dxdµ(ω) < ∞.

(

A central hypothesis is that the uncertainty obeys a certain given statistic, with probability density dµ(ω), that is dµ(ω) = 1. This information is used to derive reduced models. For this task, one determines the so-called chaos polynomials p i (ω), with d o p i = i, which are orthonormal for the measure dµ(ω) p i (ω)p j (ω)dµ(ω) = δ ij .

The terminology chaos polynomials has been coined in the seminal work of Wiener [START_REF] Wiener | The Homogeneous Chaos[END_REF]. Such methods are widely used in engineering, see for example [START_REF] Lin | The Stochastic Piston Problem[END_REF][START_REF] Bijl | Uncertainty Quantification in Computational Fluid Dynamics[END_REF] and references therein. Since the polynomials determine a Hilbert basis of L 2 µ , and using the vector space P N ω of polynomials of degree less or equal to N , any function v ∈ L 2 µ is expanded and approximated as

v = ∞ i=0 v i p i , v i = v(ω)p i (ω)dµ(ω), v (N ) = N i=0 v i p i ∈ P N ω . (3)
The approximation v (N ) has optimal spectral accuracy properties which have been proved in [7,[START_REF] Bernardi | Polynomial interpolation results in Sobolev space[END_REF]. We use the notation with parenthesis to distinguish it from other polynomial approximations to be introduced later.

In order to construct a closed system for the evolution of the moments for solutions of [START_REF] Bernardi | Polynomial interpolation results in Sobolev space[END_REF], one can use the so-called method of moments, and consider the model system of conservation laws,

∂ t u N i + ∇ • F (u N )p i (ω)dµ(ω) = 0, 0 ≤ i ≤ N. (4)
Following [8] and the references therein, such a method can be also rewritten as

           ∂ t u N + ∇ • F (N ) (u N ) = 0, F (N ) (u N ) = N i=0 F (u N (ω ))p i (ω )dµ(ω ) p i (ω), u N (x, ω, t) = N i=0
u N i (x, t)p i (ω).

(

) 5
In particular, in [8], the authors show how to model uncertainty in the coefficients F of the equations (1)-(2), using an additional variable. Even if the theory of existence, uniqueness and the theory of numerical approximation is well established for scalar conservation laws like [START_REF] Bernardi | Polynomial interpolation results in Sobolev space[END_REF], there is no such theory for systems of partial differential equations like [START_REF] Bourgault | Convergence rate for the method of moments with linear closure relations, 1-27[END_REF] or for the convergence of the solutions of (4) towards solutions of (1) parametrized by the parameter ω, see [8]. Some results may be found nevertheless in [START_REF] Gottlieb | Galerkin Method for Wave Equations with Uncertain Coefficients[END_REF] for the advection equation, in [8] where spectral convergence is proved with a weakstrong method but before any shock, and in [START_REF] Ch | Multi-level Monte Carlo finite volume methods for non-linear systems of conservation laws in multi-dimensions[END_REF] for Monte-Carlo methods applied to conservation laws. Note also the convergence proof in the recent work [START_REF] Bourgault | Convergence rate for the method of moments with linear closure relations, 1-27[END_REF], for the method of moments in the framework of kinetic equations. We also refer to [START_REF] Jin | A well balanced stochastic galerkin method for scalar conservation laws with random inputs[END_REF] for conservation laws with random right hand sides and to [START_REF] Jin | Asymptotic-preserving methods for hyperbolic and transport equations with random inputs and diffusive scalings[END_REF] for an asymptotic-preserving extension to transport equations with random inputs. One of the reason of these mathematical difficulties lies in the divorce between the standard L 1 -based theory for conservation laws, and the L 2 -based theory of spectral approximation. Similar concerns may be found in a PDE context in [START_REF] Brenier | Résolution d'équations d'évolution quasi-linéaires en dimension N d'espace à l'aide d'équations linéaires en dimension N + 1[END_REF][START_REF] Kand | Criteria on contraction for entropic discontinuities of the system of conservation laws via weighted relative entropy[END_REF], or in [START_REF] Jakeman | Enhancing l 1 minimization estimates of polynomial chaos expansions using basis selection[END_REF] for l 1 minimization algorithmic issues related to the use of polynomial chaos expansions.

To circumvent these difficulties, we propose to consider the kinetic formulation of conservation laws which has the advantage of reconciliating L 1 and L 2 theories. It writes as a Boltzmann equation for t ≥ 0, x ∈ R d and ξ ≥ 0, in a BGK (relaxation) form,

           ∂ t f ε + a(ξ).∇f ε + 1 ε f ε = 1 ε M (u ε ; ξ), a = ∇F, u ε (x, t) = f ε (x, ξ, t)dξ, f ε (t = 0) = M u init ; ξ , (6)
still assuming (2) for u init , and

M (u; ξ) = 1 I {0<ξ<u} (7)
is called a Maxwellian in the rest of this work. Notice that the non negativity u ≥ 0 is needed for this definition to make sense. That is why we assume the initial data is non negative u init ≥ 0 throughout this work. This assumption simplifies some non essential technicalities and allows us to disregard the negative part of M ; the reader can find in [START_REF] Perthame | A kinetic equation with kinetic entropy functions for scalar conservation laws[END_REF][START_REF] Lions | A kinetic formulation of multidimensional scalar conservation laws and related equations[END_REF][START_REF] Perthame | Kinetic formulation of conservation laws[END_REF] the adaptation for general sign of the initial data as well as convergence proofs of ( 6) to [START_REF] Bernardi | Polynomial interpolation results in Sobolev space[END_REF]. We recall, for later use, that M (u; ξ) is a universal minimizer for a family of entropy functionals [START_REF] Brenier | Résolution d'équations d'évolution quasi-linéaires en dimension N d'espace à l'aide d'équations linéaires en dimension N + 1[END_REF][START_REF] Perthame | A kinetic equation with kinetic entropy functions for scalar conservation laws[END_REF][START_REF] Lions | A kinetic formulation of multidimensional scalar conservation laws and related equations[END_REF][START_REF] Perthame | Kinetic formulation of conservation laws[END_REF]; for all convex functionals S(ξ),

M (u; •) = argmin u= gdξ, 0≤g≤1 S (ξ)gdξ. ( 8
)
The idea pursued in this work is firstly to write (6) for all ω, and secondly to modify it in a polynomial manner so as to consider the equation (we call it the intrusive kinetic formulation)

           ∂ t f N ε + a(ξ).∇f N ε + 1 ε f N ε = 1 ε M N u N ε ; ξ, ω , u N ε (x, ω, t) = f N ε (x, ξ, ω, t)dξ, f N ε (t = 0) = M N u init ; ξ, ω , (9)
where 0 ≤ M N (u N ε ; ξ, ω) ≤ 1 is a suitable polynomial modification of the Maxwellian M . Notice that f N ε (t = 0)dξdω = u init dω but the initial data needs not be at 'equilibrium' since u init usually does not belong to P N ω . The function f N ε is naturally a polynomial in ω of degree less or equal to N and we will show the intrusive kinetic formulation (9) of conservation laws is endowed with some convenient mathematical properties such as, for all t ≥ 0,

0 ≤ f N ε ≤ 1, f N ε (x, ξ, ω, t)dξ L ∞ xω ≤ f N ε (x, ξ, ω, 0)dξ L ∞ xω .
Notice that the solutions of ( 9) depend now on two parameters ε and N . We describe two families of methods to build the kinetic polynomial M N .

The first construction is based on convolution kernels such as Fejer or Jackson kernels and reads

M N (u N ε ; ξ, ω) = G N * ω M (u N ε ; ξ)
for some apropriate kernel G. Since this approach has natural comparison inequalities, it is possible to pass for all time t to the limit ε → 0, but with N → ∞. However εN 2 → ∞ is needed for the limit to satisfy the correct equation ( 1) for all ω, in the weak sense. Nevertheless, this result is a strong improvement with respect to [8] where the convergence is proved by means of a weak-strong method but only before the time of shock (for example if there is a discontinuity in the initial data, the result [8] simply does not apply). As shown below, a similar property holds for strong convergence, which is a corollary of Proposition 2.10. This unfortunate consequence -the limit ε → 0 independently of N is not possible-will be confirmed by considering the formal moment models which will be shown to be non consistent. Consequently this first family of results is probably of minor practical interest for numerical approximation even if endowed with rigorous approximation results.

The second family of construction intends to be less sensitive to N , so that it is possible to pass to the limit ε → 0 keeping N fixed. This will be performed with an original polynomial approximation of M , that we call a kinetic polynomial M N . It is described thanks to an algorithm based on sharp polynomial properties. The structure of M N is layered in ξ,

     M N (u N ; ξ, ω) = L l=0 h N l (ω)1 I {ξ l <ξ<ξ l+1 } , 0 = ξ 0 < ξ 1 < • • • < ξ L < ξ L+1 = u + = max I u N (ω), h N l ∈ P N ω , (10)
and is a natural option in the context of kinetic formulations. It defines an original approximation method which has the following form in dimension one for scalar conservation law

∂ t u N + ∂ x F N [u N ] = 0, F N [u N ] = a(ξ)M N (u N ; ξ, ω)dξ = l≥0 h N l (ω) (F (ξ l+1 ) -F (ξ l )) . (11)
An interesting asset of this kinetic polynomial method is the exact formula for F N [u N ], in this sense no quadrature rule is needed in opposition to the method of moments which requires to evaluate the flux F (u N )p i (ω)dµ(ω) in [START_REF] Bourgault | Convergence rate for the method of moments with linear closure relations, 1-27[END_REF]. Even if the theory of approximation is less advanced than for the convolution case, it seems to be much more adapted for the numerical approximation since this new method ( 11) is endowed with a natural maximum principle which is not the case for the general moment approximation (4). This property can be considered as one the main results of this work, see Theorem 4.1. For the Burgers flux f (u) = u 2 2 simple numerical results which satisfy the maximum principle illustrate the potential interest of this second family of methods.

The organization of this work is as follows. Section 2 is devoted to the convolution approach and the proof of Theorem 1.1. Section 2.3 deals with the design of the layered kinetic polynomial Maxwellian [START_REF] Gottlieb | Galerkin Method for Wave Equations with Uncertain Coefficients[END_REF] and the proof of Theorem 3.10 about the completion of the algorithm that constructs the layered Maxwellian. Theorem 4.1 about the stability in the maximum norm of the scheme [START_REF] Jakeman | Enhancing l 1 minimization estimates of polynomial chaos expansions using basis selection[END_REF] is proved in the last section 4. We complete this work with simple numerical illustrations of the new method [START_REF] Gottlieb | Galerkin Method for Wave Equations with Uncertain Coefficients[END_REF] for the case N = 2.

Convolution kernels

A simple idea to construct a polynomial approximation M N of M , is to use a convolution method under the form

M N (u N ε ; ξ, ω) = G N * ω M (u N ε ; ξ) := G N (ω, ω )M u N ε (ω ); ξ dµ(ω )
where the convolution kernel G N can be decomposed, using the orthonormal basis defined by (3), as

G N (ω, ω ) = N i=0 c i p i (ω)p i (ω ), (12)
where c i are appropriate coefficients, and where G N satisfies

G N ≥ 0, G N (ω, ω )dµ(ω ) = c 0 = 1 = G N (ω, ω )dµ(ω). ( 13
)
The theory of polynomial kernel approximation [START_REF] Devore | Constructive approximation[END_REF][START_REF] Weisse | The kernel polynomial method[END_REF] asserts that convolution kernels exist which satisfy the requirements ( 12)- [START_REF] Jin | A well balanced stochastic galerkin method for scalar conservation laws with random inputs[END_REF].

For example, considering the measure dµ(ω) = dω π √ 1-ω 2 , on the interval ω ∈ I = (-1, 1), they are built on the Tchebycheff orthonormal polynomials

T i (ω) = cos (i arcos ω) , -1 ≤ ω ≤ 1.
According to formula (3), a function v can be represented as the infinite series

v(ω) = µ 0 + 2 ∞ i=1 µ i T i (ω), µ i = I v(ω )T i (ω )dµ(ω ).
The following kernels are suitable modifications of the truncated Dirichlet series

v N D (ω) = µ 0 + 2 N i=1 µ i T i (ω) = G * ω v
for the choice c i ≡ 1, which generates oscillations [START_REF] Weisse | The kernel polynomial method[END_REF] and thus, is not convenient for our purposes since ( 13) is not satisfied.

Exemple 2.1 (Fejer Kernel) The Fejer kernel G N F is defined by the coefficients

c 0 = 1 and c i = 2 N + 1 -i N + 1 , 1 ≤ i ≤ N. The truncated Fejer series v N F (ω) = c 0 µ 0 + 2 N i=1 c i µ i T i (ω) is such that v ≥ 0 =⇒ v F N ≥ 0. It comes from the integral representation [9] v N F (ω) = 2π 0 v(cos(t -u))K N F (u)du = 2π 0 v(cos u)K N F (t -u)du, ω = cos t, (14)
with the kernel

K N F (u) = 1 2π(N +1) sin(N +1) u 2 sin u 2 2
Exemple 2.2 (Jackson Kernel) The Jackson integral representation writes

v N J (ω) = 2π 0 v(cos(t -u))K N J (u)du = 2π 0 v(cos u)K N J (t -u)du, ω = cos t. ( 15
)
The Jackson kernel

K N J (u) = λ N J sin(N +1) u 2 sin u 2 4
is a convenient renormalization of the square of the Fejer kernel [START_REF] Devore | Constructive approximation[END_REF], where the normalization coefficient λ N J is defined by

2π 0 K N J (u)du = 1.
By definition v N F is a polynomial of degree ≤ 2N in the variable ω. One has better approximation properties for the Jackson kernel than for the Fejer kernel, cf. [START_REF] Devore | Constructive approximation[END_REF][Theorem 2.2 page 205], since one can prove an optimal error estimate

v -v N J L 1 µ (I) ≤ Cmod 2 (v, 1 N ) ( 16
)
with mod 2 (v, α) := π 0 |v(cos(t + α)) -2v(cos(t)) + v(cos(t -α))| dt. ( 17
)
Notice that mod 2 (f, α) is an integral with respect to the trigonometric variable t ∈ (0, π), while the error is an integral with respect to the original variable ω = cos t. This is the reason of the weighted L 1 µ norm used for [START_REF] Lions | A kinetic formulation of multidimensional scalar conservation laws and related equations[END_REF]. One deduces that it also holds:

v -v N J L 1 µ (I) ≤ Cmod 1 (v, 1 N ) (18)
with

mod 1 (v, α) = 2π 0 |v(cos(t + α) -v(cos t)| dt for α > 0. ( 19
)
Exemple 2.3 (The modified Jackson Kernel) This kernel can be found in [START_REF] Meinardus | Approximation of functions: theory and numerical methods[END_REF] and is studied for modern physical calculations in [START_REF] Weisse | The kernel polynomial method[END_REF]. The polynomial is defined by

v N = c 0 µ 0 + 2 N i=1 c i µ i T i (ω), c i = (N + 2 -i) cos πi N +2 + sin πi N +2 cot π N +2 N + 2 . ( 20
)
This is also a positive kernel [START_REF] Weisse | The kernel polynomial method[END_REF] and it admits a convolution representation like [START_REF] Lin | The Stochastic Piston Problem[END_REF]. The approximation properties in L 1 µ are similar to those of the Jackson kernel due to the estimates (61)-(69) in [START_REF] Weisse | The kernel polynomial method[END_REF]. For practical computations, this kernel has two major interests. The coefficients c n are known by a simple analytical formula, and its degree is arbitrary, unlike the Jackson kernel which has even order. The corresponding kernel is denoted as K N mod,J .

Entropy and a priori bounds

We review some a priori bounds and properties satisfied by the equation ( 9) recalling that we use well-prepared initial data

M N (u init ; ξ, ω) = G N * ω M (u init ; ξ, ω). (21)
The stability estimates are the same for any of the three kernels mentioned above. However, we detail the approximation estimates only for the Jackson kernel, for which the theory of approximation is well established. The generalization to the modified Jackson kernel is left to the reader together with the study of the non optimal approximation properties of the Fejer kernel.

Proposition 2.4 For any of the three kernels, there is a unique solution of ( 9) and we have the properties 1.

0 ≤ f N ε ≤ 1. 2. 0 ≤ u N ε ≤ U M := sup x,ω u init (x, ω), f N ε (x, ξ, ω, t) ≡ 0 for ξ ≥ U M . 3.
For all smooth convex functions S(•), we have

∂ t S (ξ)f N ε (x, ξ, ω, t)dξdµ(ω) + div a(ξ)S (ξ)f N ε (x, ξ, ω, t)dξdµ(ω) ≤ 0.
4. (Contraction principle) Consider two solutions f N ε and g N ε of (9), then

f N ε (t) -g N ε (t) dxdξdµ(ω) ≤ f N ε (0) -g N ε (0) dxdξdµ(ω). 5. (Comparison principle) Also, if f N ε (0) ≤ g N ε (0),then f N ε (t) ≤ g N ε (t)
for all t ≥ 0. 6. Finally, the BV bounds in space are propagated, for all t ≥ 0,

|∇ x f N ε |dxdξdµ(ω) ≤ C init , |∇ x u N ε |dxdµ(ω) ≤ C init . (22)
Proof. 1. Indeed, from [START_REF] Jin | A well balanced stochastic galerkin method for scalar conservation laws with random inputs[END_REF] we have

G N * ω M (u N ε ; ξ) ≥ 0 and thus f N ε ≥ 0. 2. Also, we have G N * ω M (u N ε ; ξ) ≤ G N (ω, ω )dµ(ω ) = 1 and thus f N ε ≤ 1.
3. Similarly, for ξ ≥ U M , we have M (u init ; ξ) ≡ 0 and thus f N ε (t = 0) = 0. This property is propagated by the equation because

f (ξ) = 0 for ξ ≥ U M implies u = f dξ ≤ U M , which itself implies M = 0 for ξ ≥ U M and thus G * ω M = 0 for ξ ≥ U M .
4. To prove the entropy inequality, we can always assume S(0) = 0. Then, after multiplying the equation by S (ξ), integrating in ξ and ω, we notice that the right hand side is 1 ε times

S (ξ)G N * ω M (u N ε ; ξ)dξdµ((ω) -S (ξ)f N ε (x, ξ, ω, t)dξdµ(ω) = G N * ω S(u N ε )dµ(ω) -S (ξ)f N ε (x, ξ, ω, t)dξdµ(ω) = S(u N ε )dµ(ω) -S (ξ)f N ε (x, ξ, ω, t)dξdµ(ω) ≤ 0
because this is true ω by ω, thanks to the universal entropy minimisation principle (8). Here we have also used G N (ω, ω )dµ(ω) = 1.

One has

∂ t |f N ε -g N ε | + a(ξ).∇ x |f N ε -g N ε | + 1 ε |f N ε -g N ε | = sgn(f N ε -g N ε ) G * ω M (u N ε ; ξ) -M (v N ε ; ξ) ≤ G * ω M (u N ε ; ξ) -M (v N ε ; ξ) . So d dt |f N ε -g N ε |dxdξdµ(ω) + 1 ε |f N ε -g N ε |dxdξdµ(ω) ≤ 1 ε G * ω M (u N ε ; ξ) -M (v N ε ; ξ) dxdξdµ(ω) ≤ 1 ε M (u N ε ; ξ) -M (v N ε ; ξ) dxdξdµ(ω) = 1 ε u N ε -v N ε dxdµ(ω) = 1 ε |f N ε -g N ε |dxdξdµ(ω).
6. The comparison principle is a simple variant using (...) + in place of the absolute value. 7. The BV bound, is an immediate consequence of the contraction principle combined with translational invariance.

The proof is complete.

The entropy inequality of Proposition 2.4 is integrated over ω. It is possible to get a sharper estimate by using a test function ϕ ≥ 0 for the ω variable. The price is to work with weak topology. Proposition 2.5 (Estimates in weak norms in ω) Consider solutions of ( 9) with the Jackson kernel. For all smooth convex functions ξ → S(ξ) and all non negative function ω → ϕ(ω), there exists a constant C S,U M such that

∂ t ϕ(ω)S (ξ)f N ε (x, ξ, ω, t)dξdµ(ω) + div ϕ(ω)a(ξ)S (ξ)f N ε (x, ξ, ω, t)dξdµ(ω) ≤ C S,U M mod 2 (ϕ, 1 N ) ε .
Proof. We only have to evaluate the right hand side, that is 1/ε times A with

A = ϕ(ω)S (ξ)G N * ω M (u N ε ; ξ)dξdµ(ω) -ϕ(ω)S (ξ)f N ε (x, ξ, ω, t)dξdµ(ω).
The kernel being symmetric one has (normalizing S with S(0) = 0)

ϕ(ω)S (ξ)G N * ω M (u N ε ; ξ)dξdµ(ω) = S (ξ)M (u N ε ; ξ) G N * ω ϕ(ω) dξdµ(ω) ≤ S (ξ)M (u N ε ; ξ)ϕ(ω)dξdµ(ω) + S (ξ)M (u N ε ; ξ)r N (ω)dξdµ(ω)
where the convolution error r N = G N * ω ϕ -ϕ is controled as, recalling ( 16)-( 17),

r N (ω) L 1 µ ≤ mod 2 (ϕ, 1 N ).
Therefore, using the L ∞ bounds of Proposition 2.4, the last term is bounded as

S (ξ)M (u N ε ; ξ)r N (ω)dξdµ(ω) ≤ sup ω |S(u N ε )| r N (ω) L 1 µ ≤ max 0≤ξ≤U M S(ξ) mod 2 (ϕ, 1 N ).
Finally, using the universal entropy minimisation principle (8), A ≤ C S,U M mod 2 (ϕ, 1 N ) which ends the proof.

Next question is to obtain a BV bound with respect to ω. In view of the comparison principle of Proposition 2.4, it is appealing to rely on transformations which commute with the operators of the system. Fortunately, the kernels based on Tchebycheff polynomials are endowed with such a natural transformation. This is immediately visible in [START_REF] Lin | The Stochastic Piston Problem[END_REF] which is written as convolutions. This is also true for the modified Jackson kernel [START_REF] Ch | Multi-level Monte Carlo finite volume methods for non-linear systems of conservation laws in multi-dimensions[END_REF] which can also be written as a convolution [START_REF] Weisse | The kernel polynomial method[END_REF].

For any f and 0 < α, we define the linear operator (which mimicks translations)

g = Q α f ⇐⇒ g(ω) = f (cos(t + α)) -f (cos(t)) α , ω = cos t ∈ [-1, 1]. ( 23
)
The aforementioned convolution formulas have the consequence that the operators commute

Q α G N * ω = G N * ω Q α . ( 24
)
Proposition 2.6 (BV bound in ω) Consider solutions of ( 9) with any of the three kernels. As a consequence of the commutation property (24), one has that

Q α f N ε (t) dxdξdµ(ω) ≤ Q α f N ε (0) dxdξdµ(ω). ( 25
)
Proof. The proof is a slight modification of the comparison principle in Proposition 2.4. We set 9) with any of the three kernels. The time derivative is bounded as follows:

g N εα = Q α f N ε which satisfies ∂ t g N εα + a(ξ).∇g N εα + 1 ε g N εα = 1 ε Q α G N * ω M (u N ε ; ξ) = 1 ε G N * ω Q α M (u N ε ; ξ). Therefore ∂ t |g N εα | + a(ξ).∇ x |g N εα | + 1 ε |g N εα | = sgn(g N εα ) G * ω Q α M (u N ε ; ξ) ≤ G * ω Q α M (u N ε ; ξ) and d dt |g N εα | dxdξdµ(ω) + 1 ε |g N εα |dxdξdµ(ω) ≤ 1 ε G * ω Q α M (u N ε ; ξ) dxdξdµ(ω) ≤ 1 ε Q α M (u N ε ; ξ) dxdξdµ(ω) ≤ 1 ε Q α u N ε dxdµ(ω) ≤ 1 ε Q α f N ε dxdξdµ(ω) = 1 ε g N εα dxdξdµ(ω
∂ t f N ε (t) dxdξdµ(ω) ≤ ∂ t f N ε (0) dxdξdµ(ω) ≤ C U M |∇ x u init |dxdµ(ω).
Proof. This is once again a consequence of the comparison principle. Set

g N ε (t) = f N ε (t + α), with α > 0. The comparison principle of Proposition 2.4 yields that f N ε (t + α) -f N ε (t) α dxdξdµ(ω) ≤ f N ε (α) -f N ε (0) α dxdξdµ(ω).
We can pass to the limit α = 0 + using that

∂ t f N ε (0) = -a(ξ) • ∇ x f N ε (0)
for the well-prepared data under considerations. Because we have

|∇ x M N u init ; ξ, ω |dxdξdµ(ω) ≤ |∇ x M u init ; ξ |dxdξdµ(ω) ≤ |∇ x u init |dxdµ(ω),
the proof is complete.

For our next statement, we use truncation functions in x with the properties (with R a parameter allowing localization on sets as large as we wish)

χ ∈ C 2 (R d ), χ ∈ L 1 (R d ), χ > 0, χ = 1 in B R , |∇χ| ≤ χ. (26)
In other words we choose truncation functions which decay as e -|x| at infinity.

Proposition 2.8 (Estimate of derivative in ξ) Consider any of the three kernels and assume u init ∈ BV x and consider a nonnegative truncation function satisfying (26). Then, solutions of ( 9) have bounded ξ derivatives

χ(x) ∂ ξ f N ε (t) dxdξdµ(ω) ≤ e -t ε χ(x) ∂ ξ f N ε (0) dxdξdµ(ω) + C R
where C R depends on the radius in (26) and on the x -BV estimates [START_REF] Wiener | The Homogeneous Chaos[END_REF].

Proof. Indeed one can differentiate the equation with respect to ξ

∂ t ∂ ξ f N ε + a(ξ).∇∂ ξ f N ε + 1 ε ∂ ξ f N ε = 1 ε G N * ω ∂ ξ M (u N ε ; ξ) -a (ξ).∇f N ε where ∂ ξ M (u N ε ; ξ) = δ(ξ -u N ε ) is a measure. It yields ∂ t ∂ ξ f N ε + a(ξ).∇ ∂ ξ f N ε + 1 ε ∂ ξ f N ε = 1 ε G N * ω ∂ ξ M (u N ε ; ξ) + a (ξ).∇f N ε ,
Here the right hand side can be bounded. Firstly one has that

G N * ω ∂ ξ M (u N ε ; ξ) dξdµ(ω) ≤ ∂ ξ M (u N ε ; ξ) dξdµ(ω) = dµ(ω) = C 1 ,
Consequently, we find

∂ t χ(x) ∂ ξ f N ε dξdµ(ω) + a(ξ).∇[χ(x) ∂ ξ f N ε ] dξdµ(ω) + χ(x) 1 ε ∂ ξ f N ε dξdµ(ω) ≤ χ(x) C 1 ε + χ(x) a (ξ).∇f N ε dξdµ(ω) + ∂ ξ f N ε a(ξ).∇χ(x)dξdµ(ω)
Secondly one has

a (ξ).∇f N ε dxdξdµ(ω) ≤ C 2 ∇f N ε dxdξdµ(ω) ≤ C 3
where we used Proposition 2.4, and |a | ≤ C 3 because ξ ≤ U M in these integrals. Therefore, using the truncation properties of (26),

d dt χ(x) ∂ ξ f N ε dxdξdµ(ω)+ 1 ε χ(x) ∂ ξ f N ε dxdξdµ(ω) ≤ C 1 ε +C 3 χdx+C 4 χ(x) ∂ ξ f N ε dxdξdµ(ω).
This yields the result using Gronwall's lemma.

Convergence

We establish some conditions for the convergence of f N ε towards a correct limit f as N → ∞ and ε → 0. To begin with, we notice that the solutions are BV in the domain

(x, ξ, ω, t) ∈ D α = R d × [0, U M ] × [1 -α, 1 + α] × [0, T ], 0 < α < 1,
due to the previous propositions.

The restriction in direction ω comes from the transformation (23), since we recognize the standard BV criterion for the variable t but mapping it back to the variable ω, we loose the BV bound because of the degeneracy at the end points due to the vanishing derivative of the cosine function. This is why the BV property (in particular for the measure with respect to ω which writes

dω π √ 1-ω 2 ) is obtained readily only in D α for α > 0. Nevertheless the functions f N ε also belong to L ∞ . Therefore, 1. as N → ∞, ε → 0, from any sequence f N ε ∈ L ∞ ∩ D α for all 0 < α, one can extract a subsequence that converges strongly in L 1,loc xξµ(ω) (D 0 ) to a limit f which satisfies 0 ≤ f ≤ 1 and f dxdξdµ(ω) ≤ u init (x, ω)dxdµ(ω). 2. for the same subsequence, u N ε = f N ε dξ converges strongly in L 1,loc xµ(ω) (D 0 ) to the limit u = f dξ, 3. and M (u N ε ; ξ) converges strongly in L 1,loc xξµ(ω) (D 0 ) to the limit M (u; ξ), 4. For the Jackson kernel, M N (u N ε ; ξ, ω) converges strongly in L 1,loc xξµ(ω) (D 0 ) to M (u; ξ). This is because M N (u N ε ; ξ, ω) -M (u N ε ; ξ) = (G N -I) * ω [M (u N ε ; ξ) -M (u; ξ)] + (G N -I) * ω M (u; ξ).
The first term converges to 0 thanks to (13) and 3., while the second term converges to 0 because

|(G N -I) * ω M (u; ξ)|dxdξdµ(ω) ≤ C mod 1 M (u; ξ), 1 N ≤ C mod 1 M (u init ; ξ), 1
N
where we use the standard L 1

x contraction inequality for solutions of the scalar conservation law (1). A BV bound of the initial data with respect to the variable ω shows the convergence to zero of this term (with the same argument as before for the degeneracy at the endpoints), 5. and, still for the Jackson kernel, f = M (u; ξ, ω) because the equation [START_REF] Devore | Constructive approximation[END_REF] gives

M N (u N ε ; ξ, ω) -f N ε = ε[∂ t f N ε + a(ξ).∇f N ε ] → 0 in D .
Consequently, the questions of interest are to establish some conditions which imply that f = M (u; ξ) and u satisfy the correct limit equations and to obtain error estimates.

Proposition 2.9 (The limit is a weak solution) Consider the Jackson kernel. Assume that N 2 ε → ∞. Then, the full sequence f N ε (see the construction before) converges to f = M (u; ξ) and it is a weak solution of

∂ t f + a(ξ).∇f = ∂ ξ m, f (t = 0) = M (u init ; ξ), ( 27
)
where m is a non negative measure, and thus u = f dξ is the unique entropy solution of (1), that is for a.e. ω.

Proof. The theory in [START_REF] Perthame | A kinetic equation with kinetic entropy functions for scalar conservation laws[END_REF][START_REF] Lions | A kinetic formulation of multidimensional scalar conservation laws and related equations[END_REF][START_REF] Perthame | Kinetic formulation of conservation laws[END_REF] immediately shows that u is the entropy solution of the conservation law, a.e. with respect to ω, and thus uniqueness as soon as ( 27) is established (we recall the L 1 ∩ L ∞ assumption for u init following ( 2)).

To prove that (27), holds, we write

∂ t f N ε + a(ξ).∇f N ε = 1 ε M (u N ε ; ξ) -f N ε + 1 ε G N * ω -I ω M (u N ε ; ξ). (28)
The first three terms pass to the limit from the observations above. It remains to prove that

1 ε G N * ω -I ω M (u N ε ;
ξ) tends to zero in the weak sense. For this purpose, it is sufficient to use a smooth test function ϕ(x, ξ, ω, t) with compact support C, and observe that

1 ε C ϕ(x, ξ, ω, t) (G N * ω -I ω )M (u N ε ; ξ) dxdξdµ(ω)dt = 1 ε C (G N * ω -I ω )ϕ M (u N ε ; ξ)dxdξdµ(ω)dt.
To conclude the proof, we use the same argument as in Proposition 2.5, and notice that for the test function we can assume it has compact support in ω ∈ (-1, 1)

(G N * ω -I ω )ϕ L 1 (dµ) ≤ Cmod 2 (ϕ, 1 N ) ≤ C N 2
and the equation is established.

Our next result quantifies strong convergence, based on the previous comparison estimates. The compatibility relation between ε and N is nevertheless more stringent than for weak convergence since we need N ε → 0. This is the lowest rate for N that we can use to handle the right hand side of (29).

Proposition 2.10 (Strong error bounds) Consider the Jackson kernel. One has the inequalities

f N ε (t) -G N * ω f ε (t) dxdξdµ(ω) ≤ C t ε mod 1 (u init , 1 N )dxdξ, (29)
f N ε (t) -f ε (t) dxdξdµ(ω) ≤ C(1 + t ε ) mod 1 (u init , 1 N )dxdξ, (30)
f N ε (t) -M (u; ξ) dxdξdµ(ω) ≤ c √ ε + C(1 + t ε ) mod 1 (u init , 1 N )dxdξ. (31)
Proof. 1. The convolution of (6) yields the identity

∂ t g N ε + a(ξ).∇g N ε + 1 ε g N ε = 1 ε G N * ω M (u ε ; ξ), g N ε = G N * ω f ε . Set v N ε = g N ε dξ = G N * ω u ε and r N ε = G N * ω M (u ε ; ξ) -M (v N ε ; ξ) so that ∂ t g N ε + a(ξ).∇g N ε + 1 ε g N ε = 1 ε G N * ω M (v N ε ; ξ) + 1 ε r N ε .
By definition, see [START_REF] Weisse | The kernel polynomial method[END_REF], the initial conditions are the same for f N ε and g N ε :

f N ε (0) = g N ε (0) = G N * ω M (u init ; ξ) = M N (u init ; ξ, ω).
Since by definition

∂ t f N ε + a(ξ).∇f N ε + 1 ε f N ε = 1 ε G N * ω M (u N ε ; ξ)
, one more use of the comparison principle between f N ε and g N ε yields

f N ε (t) -g N ε (t) dxdξdµ(ω) ≤ 1 ε t 0 r N ε (s) dxdξdµ(ω) ds.
Using the definition of r N ε and the property ( 13), one has

r N ε (x, ξ, ω, s) dxdξdµ(ω) ≤ M (u ε ; ξ) -M (v N ε ; ξ) dxdξdµ(ω) ≤ u ε (x, ω, s) -v N ε (x, ω, s) dxdξdµ(ω) ≤ u ε (x, ω, s) -G N * ω u ε (x, ω, s) dxdµ(ω) ≤ C mod 1 (u ε (s), 1 N )dx.
The last inequality comes from the approximation property [START_REF] Lions | A kinetic formulation of multidimensional scalar conservation laws and related equations[END_REF] in the L 1 µ norm of the Jackson kernel. It remains to use the standard L 1 xξ contraction estimate for solutions of the kinetic equation ( 6). One gets mod

1 (u ε (s), 1 N )dx = 2π 0 u ε x, cos(τ + 1 N ), s) -u ε (x, cos τ, s) dx dτ ≤ 2π 0 u ε x, cos(τ + 1 N ), 0 -u ε (x, cos τ, 0) dx dτ = mod 1 (u ε (0), 1 N )dx.
Since the initial data is independent of ε, that is u ε (0) = u init , it ends the proof of the first estimate.

2. We have, using the contraction principle for the BGK equation, cf [START_REF] Perthame | A kinetic equation with kinetic entropy functions for scalar conservation laws[END_REF][START_REF] Lions | A kinetic formulation of multidimensional scalar conservation laws and related equations[END_REF][START_REF] Perthame | Kinetic formulation of conservation laws[END_REF],

G N * ω f ε (t) -f ε (t) dxdξdµ(ω) ≤ mod 1 (f ε (t), 1 N )dxdξ ≤ 2π 0 |f ε(x, ξ, ω + 1 N , t) -f ε(x, ξ, ω, t)|dxdξdω ≤ 2π 0 |f ε(x, ξ, ω + 1 N , 0) -f ε(x, ξ, ω, 0)|dxdξdω = mod 1 (u init , 1 N )dx,
which proves the second inequality 3. The third estimate is just the consequence of the second and of the standard convergence rate from the BGK model to the scalar conservation law, cf [START_REF] Perthame | A kinetic equation with kinetic entropy functions for scalar conservation laws[END_REF][START_REF] Lions | A kinetic formulation of multidimensional scalar conservation laws and related equations[END_REF][START_REF] Perthame | Kinetic formulation of conservation laws[END_REF].

Remark 2.11

The estimate (29) can be slightly enhanced under the form

f ε (t) -G N * ω f ε (t) dxdξdµ(ω) ≤ C ε mod 2 (f ε (t), 1 N )dxdξ.
However, propagation bounds hold true for mod 1 (f ε , 1 N ), but not for mod 2 (f ε , 1 N ).

Moment equations

Even if the weak or strong convergence estimates do not allow to pass to the limit ε independently of N , it is instructive in view of practical numerical computations to write the formal limit in the regime εN = O(1). The unknown of the resulting moment system are the quantities

u N ε,i (x, t) = f N ε,i (x, ω, t)dξ, f N ε,i (x, ω, t) = f N ε (x, ξ, ω, t)T i (ω)dµ(ω).
We now explain why an artificial damping phenomenon arises.

For convenience we set N + 1 = 1 ε . The projected equation for the modified Jackson kernel are

∂ t u N ε,i + div a(ξ)f N ε,i dξ = 1 ε c modJ i u N ε,i -u N ε,i = (N + 1) (N + 1 -i) cos πi N +1 + sin πi N +1 cot π N +1 N + 1 -1 u N ε,i = (N + 1 -i) cos πi N + 1 + sin πi N + 1 cot π N + 1 -N -1 u N ε,i = -h N (i)u N ε,i .
Elementary calculations show that h N (0) = 0, and that h N (x) > 0 for 0 < x < N with h N (x) → 0 for all x as N → 0. One also has that 0 < h N (i) < i for 0 < i ≤ N . It implies after integration in x

∂ t u N ε,i dx = -h N (i) u N ε,i dx =⇒ u N ε,i dx(t) = e -h N (i)t u N ε,i dx(0) =⇒ lim t→∞ u N ε,i dx(t) = 0.
This damping phenomenon of the moments i = 0 also shows up if one uses the Jackson kernel, and is even stronger starting from the Fejer kernel. This is the price to pay for the good theoretical properties. However, with this regard this formulation is less satisfactory than moment methods like (4) which do not damp. This spurious damping is not satisfactory for practical purposes and motivates the intrinsic method studied in the next section.

Kinetic polynomials

In order to define another polynomial system without damping, we introduce directly a different polynomial approximation of the Maxwellian M (u; ξ) = 1 I {0<ξ<u} , which we call 'Kinetic Polynomials'. We use a minimization procedure instead of the convolution in ω according to a construction depicted in Figures 1 and2.

We first present the minimization principles underlying the construction of kinetic polynomials. Then, we present a layered algorithm that allow to build practically and efficiently these polynomials.

M=1 ω ξ

M(ξ,ω) M=0

Figure 1: Plot of the upper limit of the support of M (ω, ξ) in the (ω, ξ) plane, identical to the graph of the function ω → u(ω). The convolution method is based on some averages on horizontal lines, and does not preserve the integral on the vertical lines which is preserved, by construction, for the minimization method of section 3.1.

Minimization principle for kinetic polynomials

The purpose in this section is to generalize the universal entropy principle (8) and to construct, for any N ≥ 0, an equilibrium M N (u N ; ξ, ω), independent of S, which satisfies

M N (u N ) = argmin g N ∈K N (u N ) S (ξ)g N dξdµ(ω), for all admissible functions S, (32)
where K N (u N ) is the set of states defined by

K N (u N ) = g N (•, •) ∈ P N ω , u N (ω) = g N (ξ, ω)dξ, 0 ≤ g N ≤ 1 , for u N (•) ∈ P N ω .
In order to reformulate this problem, we remark that

S (ξ) = ∞ 0 S (s)a s (ξ)ds, a s (ξ) = 1 I {0<s<ξ} (33)
meaning that any function S such that S ≥ 0 and S (0) = 0 is a non-negative integral of functions a s (ξ) which also satisfy a s ≥ 0 and a s (0) = 0. So we replace (32) with a family of similar problems

M N (u N ) = argmin g N ∈K N (u N ) ∞ ξ g N dsdµ(ω), ∀ξ. (34)
Since the mass must be preserved, that is g N (s, ω)dsdµ(ω) = u N (ω)dsdµ(ω), this problem can be rewritten with the alternative formulation

M N (u N ) = argmax g N ∈K N (u N ) ξ 0 g N dsdµ(ω), ∀ξ. (35)
Remark 3.1 Even if the formulation (35) is simpler than (32), it is still involved. Indeed M N is required to be the single maximizer of an infinite number of maximization problem, i.e. for all ξ.

Considering the theory of polynomial approximation, it is possible to imagine that a clever use of (35) combined with the polynomial structure could generate one single maximization problem that encompasses all the properties of M N . But it is still an open problem to find this unique maximization formulation.

In the case N = 0, the formulation is simply (8). As shown in the literature [START_REF] Brenier | Résolution d'équations d'évolution quasi-linéaires en dimension N d'espace à l'aide d'équations linéaires en dimension N + 1[END_REF][START_REF] Perthame | A kinetic equation with kinetic entropy functions for scalar conservation laws[END_REF][START_REF] Lions | A kinetic formulation of multidimensional scalar conservation laws and related equations[END_REF], the degeneracy question is untied by using a strictly convex entropy S > 0, and noticing that the unique minimum is the same for all strictly convex entropy S > 0. But in our general case N > 0, the formulations (32), ( 34) and ( 35) are definitely linearly degenerate in the variable ω. This is the reason of the difficulties encountered for the moment with these maximization formulations. This is why we try less to analyze the abstract properties of this problem in terms of existence and uniqueness of the solutions, but more to reformulate it in a way adapted to design a constructive algorithm. Hence we will be able to define another formulation, denoted below as the second maximization problem, which reveals to be constructive and is the one that we use in the sequel.

Let u + = max ω∈I u N (ω) and u -= min ω∈I u N (ω). We notice that K N (u N ) is non empty since

g N = 1 I {0<ξ<u -} + u N (ω) -u - u + -u - 1 I {u -<ξ<u + } ∈ K N (u N ).
This function vanishes for u + < ξ and is identically equal to 1 for 0 < ξ < u -. Considering (34) and (35), we deduce that if M N (u N ) exists, it also vanishes identically for u + < ξ and must be identically equal to 1 for ξ < u -. Therefore we introduce these constraints in the set of possible solutions

S N (u N ) = g N ∈ K N (u N ) : g N = 1 for 0 < ξ < u -and g N = 0 for u + < ξ where u N (•) ∈ P N ω and 0 ≤ u N . Since g N ∈ S N (u N ), one has S N (u N ) = ∅.
We obtain the constraint maximization problem.

Problem 3.2 (First maximization problem) Find a unique Maxwellian polynomial M N (u N ) ∈ S N (u N ) such that ξ 0 M N dsdµ(ω) ≥ ξ 0 g N dsdµ(ω), ∀ξ, ∀ g N ∈ S N (u N ).
This is a family of maximization problems of a linear functional, over a convex set S N (u N ). An interpretation is that if a solution exists for (34) or (35), then the mass at each level ξ → M N (u N ; ξ, ω)dµ(ω) is maximized for smaller ξ: it is already a property of the solution of the initial problem (8). The whole problem is to construct the maximizer (if it exists) for values of ξ between u -and u + . As stressed in remark 3.1 the existence and uniqueness of a solution to this problem is nevertheless an open problem. We will build on the idea that the problem must be analyzed level after level (i.e. layer after layer) to construct layered feasible solutions under the form:

M N (u N ; ξ, ω) = L l=0 h N l (ω)1 I {ξ l <ξ<ξ l+1 } with 0 = ξ 0 < • • • < ξ L+1 = u + and h N l ∈ P N ω polynomial
inside each layer. See Figure 2.

In what follows, we assume that the problem 3.2 admits a unique solution, and detail some immediate consequences which explains the interest of this problem. Proposition 3.3 Under the assumption that a solution exists to the maximization problem 3.2, then it is a minimizer of (32) (that is for any S convex).

Proof. This is a consequence of the identity (33). Proposition 3.4 Cconsider a nonnegative polynomial approximation of the initial data u init,N ≥ 0. Under the assumption that a solution exists to the maximization problem 3.2, then the solution of the kinetic equation

           ∂ t f N ε + a(ξ).∇f N ε + 1 ε f N ε = 1 ε M N (u N ε ; ξ, ω), u N ε (x, ω, t) = f N ε (x, ξ, ω, t)dξ, f N ε (t = 0) = M N (u init,N ; ξ), (36)
satisfies the maximum principle. For all smooth convex functions S(•), we have the entropy inequality

∂ t S (ξ)f N ε (x, ξ, ω, t)dξdµ(ω) + div a(ξ)S (ξ)f N ε (x, ξ, ω, t)dξdµ(ω) ≤ 0.
Proof. Immediate. Notice the assumption on the non negativity of the initial data can easily be realized by using the previous method with convolution kernels:

u init,N = G N * ω u init .
Proposition 3.5 (Derivation of the polynomial system) Under the assumption that a solution exists to the maximization problem 3.2, and if u N ε converges strongly to some u N , then we can pass to the limit ε → 0 in (36) and obtain the system of conservation laws

∂ t u N i + div F N i [u N ] = 0, 0 ≤ i ≤ N, F N i [u N ] := a(ξ)M N (u N ; ξ, ω)T i (ω)dξdµ(ω), ( 37
)
with the entropy inequalities, for all smooth convex function S(•),

∂ t S N [u N ] + div G N [u N ] ≤ 0,
where the entropy and entropy fluxes are defined by

S N [u N ] := S (ξ)M N (u N ; ξ, ω)dξdµ(ω), G N [u N ] := S (ξ)a(ξ)M N (u N ; ξ, ω)dξdµ(ω).
Proof. Because of the bounds 0 ≤ f N ε (x, ξ, ω, t) ≤ 1, we may extract subsequences such that f N ε

f N ∈ P N ω in L ∞ -w * . Because, u N ε converges strongly, M N (u N ε
) also converges strongly to M N (u N ). From the equation (36) (once multiplied by ε), we conclude that f N = M N (u N ). To find the system of conservation laws, it remains to integrate in ξ the equation (36) and pass to the limit. Similarly, we pass to the limit in the entropy inequality.

This system, which is equivalent to [START_REF] Jakeman | Enhancing l 1 minimization estimates of polynomial chaos expansions using basis selection[END_REF] once projected on the orthonormal basis, has a much better structure than those in section 2.3 because they do not contain relaxation terms in the right hand side. Also, the family of entropy inequalities explains that it satisfies the maximum principle.

The rest of this work is devoted to construct a reasonable, or feasible, solution to the maximization problem. It is possible to prove the following results. For N = 0, the solution of the first maximization problem 3.2 exists, is unique and is of course equal to M (u) (single layer). For N = 1 the solution of the first maximization problem 3.2 exists, is simple to construct in two layers and is also unique. This is explained in section 3.2.

For N > 1, we reformulate the problem. The idea is that if one determines the function M N (u N ; ξ, ω) under a certain threshold ξ < ξ * , then the maximization formulation 3.2 implies that

ξ * +α 0 M N (u N ; ξ, ω)dξdµ(ω) ≥ ξ * 0 M N (u N ; ξ, ω)dξdµ(ω) + ξ * +α ξ * g N (ξ, ω)dξdµ(ω), α > 0, that is ξ * +α xi * M N (u N ; ξ, ω)dξdµ(ω) ≥ ξ * +α ξ * g N (ξ, ω)dξdµ(ω), α > 0.
The limit α → 0 + gives another characterization as a maximization problem one layer after the other:

Problem 3.6 (Second maximization problem) Assume M N (u N ) ∈ S N (u N
) is a solution of the first maximization problem 3.2. Then it is solution of another maximization problem

lim α→0 + 1 α ξ+α ξ M N (ξ, ω)dξdµ(ω) ≥ lim α→0 + 1 α ξ+α ξ g N (ξ, ω)dµ(ω), ∀ξ (38)
and for all g N such that g

N (ξ , ω) = M N (ξ , ω)1 I {0<ξ <ξ} + g N (ξ , ω)1 I {ξ<ξ <ξ+α} ∈ S N (u N ).
Remark 3.7 Like the first maximization problem, this is an infinite family (that is for all ξ) of maximization problems of a linear functional, over a convex set S N (u N ). But the difference if that it is now ordered, in the sense that we can try to maximize for small ξ, and after that to maximize on a small layer ξ + ε, and so on. We will use this principle to construct the unique solution of the second maximization problem 3.6 with a layered technique: that is we will construct M N (u N ) step by step (layer by layer).

A natural interpretation is that M N (u N ) is a function that tries to have as much mass as possible for small ξ. If we see ξ = 0 as where a ground state M N (0, ω) = 1 is, and the function M N (u N ) as some generalized "particles" distribution like in statistical physics, then the algorithm just tries to pill up "particles" above the ground state.

N = 1

A preliminary result is the following, which is the consequence that polynomials of degree 1 can be characterized by their values at the two end points. Proposition 3.8 Consider the case N = 1 and I = [-1, 1]. Denote the constraint as u 1 (ω) = αω + β ≥ 0 and assume for instance that α ≥ 0. Then, the solution of the first maximization problem 3.2 exists, is unique with the the layered form

M 1 (u 1 ; ξ, ω) = 1 I {0<ξ<min(u(-1),(u(1))} + 1 + ω 2 1 I {u(-1)<ξ<u(1)} + 1 -ω 2 1 I {u(1)<ξ<u(-1)}
where one of the two last terms vanishes.

Proof. For a first order polynomial in ω, g 1 (ξ, ω), we set c(ξ) = g 1 (ξ, -1) and d(ξ) = g 1 (ξ, 1), so that

   g 1 (ξ, ω) = c(ξ) 1-ω 2 + d(ξ) 1+ω 2 , g 1 ∈ S 1 (u 1 ) ⇐⇒ c(ξ)dξ = u(-1), d(ξ)dξ = u(1), 0 ≤ c(ξ) ≤ 1 and 0 ≤ d(ξ) ≤ 1.
The bounds 0 ≤ c, d ≤ 1 guarantee that the linear polynomial is in bounds 0 ≤ g 1 ≤ 1. We define the positive coefficients γ = 1-ω 2 dµ(ω) and δ = 1+ω 2 dµ(ω), with

S (ξ)g 1 (ξ, ω)dω = γ S (ξ)c(ξ)dξ + δ S (ξ)d(ξ)dξ.
Therefore the maximization problem N = 1 is equivalent to two separate maximization problems N = 0 like (8). The unknowns are c and d with the constraints 0 ≤ c, d ≤ 1. The solution is given by c(ξ) = 1 I {0<ξ<u(-1)} and d(ξ) = 1 I {0<ξ<u(1)} which ends the proof.

The general case N > 1

We design the solution of the second maximization problem 3.6 hereafter by a constructive method (an algorithm) under the form

M N (u N ; ξ, ω) = l≥0 h N l (ω)1 I {ξ l <ξ<ξ l+1 } , 0 = ξ 0 < ξ 1 < • • • < ξ L < ξ L+1 = u + = max I u N (ω). ( 39
)
The construction shows the uniqueness of the solution. The layer structure of this function is illustrated in Figure2. The integral identity

u + 0 M N (u N ; ξ, ω)dξ = u N (ω) writes l≥0 (ξ l+1 -ξ l )h N l (ω) = u N (ω), ω ∈ I. ( 40
)
This function is constructed step by step departing from the bottom, the first step being trivial. The second step is the critical one where all the ideas of the method are carefully explained, in particular the role of the Bojavic-Devore theorem for one sided approximation. The other steps are designed with the same method, and Theorem 3.10 guarantees the completion of the algorithm.

with the layer size

ξ 1 = u -. (42)
After this stage, it remains to construct the next layers. The integral relation (40) written for the next layers l = 1, 2, . . . becomes

l≥1 (ξ l+1 -ξ l )h N l (ω) = v N (ω) := u N (ω) -u -, ω ∈ I. (43)
By construction v N reaches its minimum 0 at some points denoted as ω i , 1 ≤ i ≤ p. It also reaches its maximum D = u + -ξ 1 at some points denoted as µ j , 1 ≤ j ≤ q, and one has that

0 ≤ 1 D v N ≤ 1. ( 44
)
In the sequel the points ω i and µ j are called points of contact, as illustrated in Figure 3. We naturally define an integer referred to as the local order of the contact. It is an even number if ω j (resp. µ i ) is inside the interval I. It can be any non zero natural number if ω j = ±1 (resp. µ i = ±1) is on the extremities. The local order of contact at ω i (rest. µ j ) will be denoted as r i + 1 (resp. s j + 1). The total order of contact of the polynomial 1 D v N is defined by i (r i + 1) + j (s j + 1). These notions play critical role in the construction of h N 1 which is detailed below.

N ω µ ω ω Here ω 1 has order 1, ω 2 has order 2 and µ 1 has order 2. The total order is 5 = 1 + 2 + 2. Also is represented at ω 2 a function which has locally a contact order greater or equal to the one of the main function.

Construction of h N

1 in the first layer The two unknowns that will be defined at the end of this section are the polynomial h N 1 (ω) and ξ 2 which gives the length ξ 2 -ξ 1 of the layer. Due to the points of contact, this construction, detailed in 6 steps, is largely constrained.

Step of the constraints: Let us examine the function

g N (ξ, ω) = l≥1 h N l (ω)1 I {ξ l <ξ<ξ l+1 } .
Since we desire to impose

D ξ 1 g N (ξ, ω)dξ = l≥1 (ξ l+1 -ξ l )h N l (ω) = v N (ω) for all ω ∈ I, and 0 ≤ h N l ≤ 1 for all l, (45)
it means that at all points of contact and for all l ≥ 1 h N l (ξ, ω i ) = 0 for 1 ≤ i ≤ p, and h N l (ξ, µ j ) = 1 for 1 ≤ j ≤ q.

(46)

The ω j and µ i , defined at the end of the previous step, are the points of contact of the graph of 1 D v N at values 0 and 1. In particular the points of contact of 1 D v N are also points of contact for h N 1 . Actually it can be proved that the local contact orders for h N 1 are greater or equal to the contact orders for 1 D v N . The example depicted in figure 3 illustrates that property: Let us assume that ω 2 ∈ (-1, 1) is a contact order r 1 = 4 for v N , and the contact order of h N 1 is less than 4, let us take 2. In this case we can write h

N 1 (ω) = c(ω -ω 2 ) 2 + O(ω -ω 2 ) 3 , c > 0. On the other hand v N (ω) = c(ω -ω 2 ) 4 + O(ω -ω 2 ) 5 , c > 0. In this case it is sure that ∀τ > 0, ∃ ω close to ω 2 with v N (ω) < τ h N 1 (ω).
So it will not be possible to enforce (45). This is a contradiction. It means the local expansion of h N 1 at the point of contact ω 2 is of the form

h N 1 (ω) = c(ω -ω 2 ) 2p + O(ω -ω 2 ) 2p+1 , c > 0, p ≥ 2.
This situation, when the contact order of h N 1 is greater or equal to the one of v N , is represented in dashed in the figure 3. This example is easily generalized to all other points of contact: we obtain that all points of contact of v N are also points of contact for h N 1 with a contact order greater or equal to r i (resp. s j ).

Let us now consider the polynomial H =1 D v N . By construction one has that

(h N 1 -H) (r) (ω i ) = 0 0 ≤ r ≤ r i 1 ≤ i ≤ q, and (h N 1 -H) (s) (µ j ) = 0 0 ≤ s ≤ s j 1 ≤ j ≤ q.
It implies that h N 1 can be decomposed as

h N 1 (ω) = H(ω) + W (ω)r(ω) ( 47
)
where W is a polynomial which can be written as

W (ω) = Π p 1 |ω -ω i | r i Π q 1 |ω -µ j | s j , deg(W ) ≤ N. (48)
W is indeed a polynomial over I due since r i and s j are even for interior points of contact.

Step 2-Maximization under constraints: We define the admissible set

S 1 = r ∈ P N -deg(W ) , 0 ≤ H + W r ≤ 1 .
Since 0 ≤ H (see 44), one has that 0 ∈ S 1 which is non empty. Therefore the maximization problem

r = argmax r∈S 1 (H(ω) + W (ω)r(ω)) dµ(ω) (49)
makes sense. It can be recast also as

r = argmax r∈S 1 W (ω)r(ω)dµ(ω). ( 50
)
Step 3-Design of h N 1 : The integral (49) defines a linear functional. The set S 1 is non empty, compact and convex. So there exists at least one solution to the maximization problem (49). The Bojavic-Devore theorem [START_REF] Bojanovic | On polynomials of best one side approximation[END_REF][Theorem 3] yields the uniqueness of the solution. Actually the original Bojavic-Devore theorem considers only one bound, but it is immediate to extend 1 the technical part [START_REF] Bojanovic | On polynomials of best one side approximation[END_REF][Lemma 3] of its proof to the case with two bounds, a lower bound and an upper bound, as in g ≤ P ≤ f assuming g < 0 < f .

Step 5-Proof that ξ 2 ≤ u + : It stems directly from (56). Indeed (54) yields

(ξ 2 -ξ 1 ) -D + v N ≤ (ξ 2 -ξ 1 )h N 1 ≤ v N =⇒ ξ 2 ≤ D + ξ 1 = u + .
In other words if ξ = u + the algorithm just completed and there is nothing more to say. In the other case ξ 2 < ξ 1 + D = u + . It also means that one can forget the u + in the definition (56) of ξ 2 .

Step 6-Increase of the contact order. We want to show that the total contact order of v N is strictly greater than the total contact order of v

N . Since v N = v N -(ξ 2 -ξ 1 )h N
1 and h N 1 has all the points of contact of v N (with greater or equal multiplicity order), then the total contact order cannot decrease. It remains to show that is increases strictly.

We analyze this property for the case ξ 2 < D, that is either min

ω∈I w N (ω) < D or min ω∈I z N (ω) < D.
So let us assume min ω∈I w N (ω) ≤ min ω∈I z N (ω) and that min ω∈I w N (ω) < D. The other case is symmetric

. Set d = ξ 2 -ξ 1 < D so that 0 ≤ v N (ω) -dh N 1 (ω) ∀ω ∈ I. One also has for all sufficiently small ε > 0, ∃ω ε ∈ I, v N (ω ε ) -(d + ε)h N 1 (ω ε ) < 0. (57)
Notice that ω ε = ω i , if not it is in contradiction with (57) since both polynomials vanish at ω i . Passing to the limit after extraction of a subsequence in I which is a closed bounded interval, there exists ω * ∈ I such that

ω ε → ω * and v N (ω * ) -dh N 1 (ω * ) = 0. ( 58
)
The discussion considers two cases.

• If ω * = ω i for all i, where ω i is a point of contact of v N , then it adds one more point of contact in the list of the ω i 's. In this case the property is proved.

• The other case when ω * = ω i for a given i needs a little more work. Assume the local expansion of v N can be written

v N (ω) = c(ω -ω i ) r i +2 + O(ω -ω i ) r i +3 , 0 < c, 0 < r i ,
where -1 < ω i < 1 for the simplicity of notations, so that r i is necessarily even. The local expansion of h N 1 can be written under a similar form

h N 1 (ω) = c(ω -ω i ) r i +2 + O(ω -ω i ) r i +3 , 0 < c, 0 < r i ≤ r i . But in view of (57-58), it is necessary that r i = r i . Then v N (ω ε ) -(d + ε)h N 1 (ω ε ) = (c -(d + ε) c)(ω ε -ω i ) r i +2 + O(ω ε -ω i ) r i +3 Using (57) and ω ε -ω i = 0, we see that (c -(d + ε) c) + O(ω ε -ω i ) < 0.
Passing to the limit we obtain that c -d c ≤ 0. But at the same time 0

≤ v N -dh N 1 . In view of the local expansion v N (ω ε ) -dh N 1 (ω ε ) = (c -d c)(ω ε -ω i ) r i + O(ω ε -ω i ) r i +1 it implies 0 ≤ c -d c. Therefore c -d c = 0, meaning that v N (ω ε ) -dh N 1 (ω ε ) = O(ω ε -ω i ) r i +3
. It shows that the local contact number is increased at least by 1.

In all cases there is strict increase of the total contact number.

detailed formula is (37), accordingly the construction presented in the previous section. The value at next time step t n+1 = (n + 1)∆t in cell j is denoted with a bar u N j ∈ P N (ω). We assume the initial data is a positive and bounded polynomial. With our notations it can be written as 0 ≤ U m ≤ u N j (ω) ≤ U M < ∞, ∀j and ∀ω ∈ I.

We consider the archetype of a convex flux which is the Burgers flux F (ξ) = ξ 2 2 . This is compatible with the upwinding of the discrete spatial derivative visible in (59) which corresponds to the characteristic lines ω per ω of a non negative initial data.

The following result states that the explicit Euler scheme satisfies the maximum principle (this is a minimal stability requirement) under a CFL condition which is independent of N . The property is here checked directly on the scheme (59) but can also be derived as a consequence of the underlying kinetic formulation.

Proof. Using v N j = F N [u N j ] to simplify the notations, one notices that by construction

u N j (ω) 2 2 ≤ v N j (ω) ≤ U 2 M 2
, where the lower bound is, by the standard Brenier inequality, a consequence of (8). One rewrites the Euler scheme as The condition 0 ≤ α N j (ω) ≤ 1 yields the propagation of the propagation of the upper bound u N j ≤ U M . This inequality is guaranteed under the CFL condition.

To prove the lower bound in (61), we need a sharper upper bound on v j (ω). Let us start with A Brenier type inequality [START_REF] Brenier | Résolution d'équations d'évolution quasi-linéaires en dimension N d'espace à l'aide d'équations linéaires en dimension N + 1[END_REF] yields the needed upper bound

v j (ω) ≤ U 2 m 2 + U M U M -(u N (ω)-Um) ξdξ = U 2 m 2 + U 2 M 2 - U M -u N (ω) -U m 2 2 ≤ U 2 m 2 + U M (u N (ω) -U m ).
Plugging in (62), we get evident at time t = 0.1 > T . The moment method generates oscillations. On the contrary, the new method based on kinetic polynomials satisfies the maximum principle at all times.

u N j (ω) ≥ u N j (ω) - ∆t ∆x U 2 m 2 - ∆t ∆x U M (u N j (ω) -U m ) + ∆t ∆x v N j-1 (ω).
It is interesting also to compare with the non intrusive method which is fierce competitor due to its simplicity and versatility in engineering. This last method amounts to decide quadrature points denoted as ω i with 1 ≤ N + 1, to solve in parallel N + 1 standard conservation laws, and after that to reconstruct the function in the ω direction in order to get an approximation for all ω. This method also has its difficulties. The first one is the approximation of the initial data if it is discontinuous in ω. So if a quadrature point if at the discontinuity, the initialization is ambiguous. The second one is the oscillations may show up after reconstruction with Lagrange interpolation techniques.

In the test, we use the general prescription for quadratures. We take the roots of the Tchebycheff polynomial T 3 (ω) = 4ω 4 -3ω: ω 1 = - 3 4 , ω 2 = 0 and ω 3 = 3 4 . We do not perform a test for the first initial data (69) due to the ambiguity of the initialization at ω 2 . The results for the second problem (70) are displayed in figure 6. Without any surprise, we observe that the maximum principle is not satisfied, an artefact also observed with the moment method.

Figure 2 :

 2 Figure2: Layered structure of the kinetic polynomial M N , a polynomial of degree less or equal to N in ω which is constant in each layer (ξ i , ξ i+1 ). The vertical line indicates that the integral on the line is equal to the integral on the same line in Figure1.

3. 3 . 1

 31 Construction of h N 0 in the lower layer The solution in the lower layer is by definition of S N h N 0 (ω) = 1 ∀ω,

1 2 1 vFigure 3 :

 13 Figure 3: Example of contact points.Here ω 1 has order 1, ω 2 has order 2 and µ 1 has order 2. The total order is 5 = 1 + 2 + 2. Also is represented at ω 2 a function which has locally a contact order greater or equal to the one of the main function.

Theorem 4 . 1

 41 Assume the CFL condition U M ∆t ≤ ∆x. Then U m ≤ u N j (ω) ≤ U M , ∀j and ∀ω ∈ I.

  α N j (ω))u N j (ω)+α N j (ω)U M , α N j (ω) = ∆t ∆x × u N j (ω) + U M 2 .

  M N j (u N ; ξ, ω) = 1 for ξ < U m and U M ξ=Um M N j (u N ; s, ω)ds = u N j (ω) -U m . N ; s, ω)ds.

5 ' 4 Figure 4 :

 544 Figure 4: Burgers flux. Illustration of the maximum principle satisfied by the method based on kinetic polynomials for the initial data (69). The usual moment solution does not satisfy the maximum principle which is a salient property of the exact solution.

1 Figure 5 :

 15 Figure 5: Burgers flux. Illustration of the maximum principle satisfied by the new method, for the initial data (70). The usual moment solution does not satisfy the maximum principle which is a salient property of the exact solution.

  Theorem 1.1 (Convolved Maxwellian) Consider the Jackson kernel. Assume the initial data is BV in all variables. One has convergence in L 1 loc of f N ε towards f with a rate of strong convergence O 1 εN , and a rate of weak convergence O 1 εN 2 .

  ).

	So that we conclude d dt	|g N εα |dxdξdµ(ω) ≤ 0, and thus |g N εα (t)|dxdξdµ(ω) ≤ |g N εα (0)|dxdξdµ(ω).
	which is the announced result.
	Proposition 2.7 (BV bound in time) Consider solutions of (

A sketch of the proof is as follows, using the notation of the seminal Bojavic-Devore paper. The key part concerns the polynomial P + ηQ defined bottom of page 146, where by definition g ≤ P ≤ f . To take into account the lower bound which is needed in our formulation, it is needed to check that g ≤ P + ηQ for all x. The verification goes as follows. Case a) If x ∈ [a, b]/I, then Qε(x) ≥ 0 by construction. So g ≤ P (x)+ηQ (x). Case b). If x ∈ I, then Qε(x) ≤ 0 by construction so it is the dangerous case. But at the same time P (xν ) = f (xν ) where xν the point of contacts of P on f . Since f ≥ α > 0 and for all x ∈ I, there exists xν ∈ I such that |x -xν | ≤ ε, it is sufficient to take ε > 0 as small as needed to guarantee that P (x) ≥ β > g for all x ∈ I. Taking η > 0 small enough ends the proof of the technical Lemma 3. The same method works for the Lemma 4. The rest of the proof is the same.

Theorem 3.9 (Bojavic-Devore [START_REF] Bojanovic | On polynomials of best one side approximation[END_REF]) Let w be a non negative Lebesgue integrable weight function:

w(ω)dω > 0. Let g < 0 < f be two functions both differentiable on (a, b), except possibly at a finite number of points τ c (c = 1, . . . , C) where f (τ c ) = +∞ or g(τ c ) = -∞. Let n ≥ 0 and Q n (z) = {r n ∈ P n , 0 ≤ r n ≤ f }. Then, there exists a unique maximizer r n = argmax s n ∈Q n (z) b a w(ω)s n (ω)dω. Considering the maximization problem (50) with the definition (48) of W , we can apply the Bojanic-Devore theorem with the weight w(ω) = W (ω) dµ(ω) dω , the function f = 1-H W > 0, the function g = -H

W < 0 and the degree is n = N -degW . Notice that since W has all the points of contact of H with orders, then g < 0 < f . The value of the function f is infinite at a root of W which is not a root of 1 -H (with multiplicity). Similarly the function -g is infinite at a root of W which is not a root of -H (with multiplicity). So that all the assumptions of the theorem are fulfilled. It defines a unique r solution of ( 49)-( 50) and we set h N 1 = H + W r. By definition, the solution satisfies the comparison inequality

Step 4-Design of ξ 2 : Now that h N 1 has been determined, it remains to decide of a value for

Whatever the value of ξ 2 , we will have to design the next layers so that

This problem has the same structure than the original problem over a smaller length D = D -(ξ 2 -ξ 1 ) and for a function

It can be solved with the method we used to construct

Indeed this inequality is the same as (44) which is at the starting point of the construction of the solution in the this layer. So if we manage to guarantee (53) it will be able to continue the construction. Therefore we define ξ 2 as the largest value so that (53) is true. It writes

By construction the polynomial h N 1 divides the polynomial v N , and has the same points of contact (with equal of greater order of multiplicity) of v N . We define two rational fractions

In general, w N and z N are polynomials, but they may be rational fractions if the contact order of h N 1 at a given point is strictly greater than the contact order of v N at the same point. We take the largest as possible ξ 2 which satisfies (54), that is

(56)

Next steps

The problems in the second layer, and in the third layer, fourth layer, . . . , have exactly the same structure as in the first layer. It constructs by iteration a sequence (v N l , ξ l , h N l ). The first steps are

. The next steps are constructed with the method used to construct (v N l , ξ l , h N l ) from (v N 0 , ξ 0 , h N 0 ). Some properties are guaranteed by construction

• ii) one has the bound 0

• iii) the contact number of h N l (and v N l ) increases strictly at each step,

• iv) the generalization of (51) holds

Theorem 3.10 (Completion of the algorithm) The method ends in a finite number of steps L ≤ N , and the layers fill the interval [0,

By construction the function M N (u N ; ξ, ω) = l≥0 h N l (ω)1 I {ξ l <ξ<ξ l+1 } is a solution of the second maximization problem 3.6.

Proof. The proof is by contradiction. Assume the algorithm never completes. It construct a infinite series of layers, all of them [ξ l , ξ l+1 ] ⊂ [0, u + ] with ξ l < ξ l+1 . Since the total contact number increases at least by 1 at each level, it reaches N + 1 or more at a certain step L of the algorithm. N + 1 contacts imply N + 1 equality constraints for h N L . There is only one possibility

The rest of the proof is evident. The completion property of the kinetic polynomials is also related to the identity u N (ω) = M N u N (ξ, ω)dξ. Applying this identity for any µ j which is an upper point of contact at every layer, one obtains

4 A simple numerical scheme for the projected equations [START_REF] Jakeman | Enhancing l 1 minimization estimates of polynomial chaos expansions using basis selection[END_REF] We discretize in time and space and implement the method under the form

where u N j ∈ P N (ω) is a polynomial in ω of degree N (fixed), in cell j and at the current time step t n = n∆t and the generic flux F N [u N j ] is constructed with the kinetic polynomial formula [START_REF] Jakeman | Enhancing l 1 minimization estimates of polynomial chaos expansions using basis selection[END_REF], a more It yields

where the last inequality is a consequence of the CFL condition. By iteration, the maximum principle (61) holds true one step after the other and the proof is complete.

We now show an explicit construction of the solution of the Theorem 3.10 for the case N = 2, and we use it to obtain preliminary numerical results for the corresponding projected equations.

Construction of the kinetic polynomial for N = 2

Let u 2 ∈ P 2 be a non negative polynomial u 2 (ω) = aω 2 + bω + c ≥ 0. We assume a = 0 since the solution for a = 0 can be determined by the method for N = 1. We first determine u + and u -. Let ω 0 = -b 2a be the solution of u (ω 0 ) = 0.

• If ω 0 ≥ 1 and a > 0:

• If ω 0 ≤ 1 and a < 0: u + = u 2 (-1) = a -b + c and u -= u 2 (1) = a + b + c.

• If -1 < ω 0 < 1 and a > 0: u + = max(u 2 (-1), u 2 (1)) = a+|b|+c and and u -= u 2 (ω 0 ) = c 2 -b 2 4a .

• If -1 < ω 0 < 1 and a < 0:

and and u -= min(u 2 (-1), u 2 (1)) = a -|b| + c.

We introduce the notations:

D . In case the absolute extrema of u 2 strictly belongs to I, it is the same for w 2 . The only possibility to reconstruct a second order polynomial when prescribing one value at ±1, one value at ω 0 and the zero derivative at ω 0 is to keep the same polynomial. In this case we set h 2 1 = w 2 and

In the other case, w 2 has necessarily its extrema at the boundary. Assume w 2 (-1) = 0 and w 2 (1) = 1. We factorize

In view of the representation formula, the maximum of the integral is for r = 1 4 , and the minimum is for r = -1 4 . Let ξ 1 be the solution of

), which is given by ξ 1 = 2D(r + 1 4 ) where D = u + -u -. One obtains

The case where w 2 (-1) = 1 and w 2 (1) = 0 is deduced by symmetry

With these formulas one can compute R + ξM 2 (ξ, ω)dξ. One finds, in the first case (63),

In the second case (64) the flux is

The third case (65) is symmetric and yields

Numerical examples

We compare numerical results obtained for N = 2 with the usual moment method, with the new method based on kinetic polynomials and also with a non intrusive method. We use the measure

1-ω 2 and consider first the Burgers equation. Also we restrict the numerics to N = 2. The moment method, that is solving (4), is easy to compute since the orthonormal basis is formed of the Tchebycheff polynomials and the Burgers flux is also a polynomial. It writes

The kinetic polynomials method use the fluxes defined in (66-68). We consider the initial data

The exact solution is a shock at velocity 2 for ω < 0, and another shock at velocity 3 for 0 < ω. The results plotted in Figure 4 show the gain in term of stability of the new method with respect to the moment method. We consider still the Burgers equation, but with another initial data, which is continuous,

The exact solution is a compressive ramp on all lines, and a shock at time T = 1 11 . Therefore the exact solution is continuous in x and ω directions for t < T , and is discontinuous in the ω direction for T < t.

Notice that we preprocess the initial data with the modified Jackson kernel in order to make sure the initial data is in bounds. As one can see it in Figure 5, the compressive nature of the solution is


Thierry Menand
email: t.menand@bristol.ac.uk

Michel De Saint-Blanquat
email: michel@lmtg.obs-mip.fr

Catherine Annen

Thierry Menand1

Michel De Saint-Blanquat2

Catherine Annen1

Emplacement of magma pulses and growth of magma bodies

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Introduction

Magmatism is responsible for the formation of continental and oceanic crusts. It is the main agent of mass and heat transfers from the mantle towards the crust, the hydrosphere (oceanic and continental hydrothermalism) and the atmosphere (emissions of volcanic gas and ashes). Its expressions are the crystallization of intrusive rocks and the eruption of volcanic products. However, with a ratio between the volumes of extrusive and intrusive magmatic rocks estimated to be of the order of 1:5 for most magmatic systems [START_REF] White | Long-term volumetric eruption rates and magma budgets[END_REF], one of the main characteristics of magmatism is that magmas rarely reach the Earth's surface. This suggests that the prevailing conditions in the crust are not favourable to the arrival of magma on the Earth's surface, but instead lean heavily towards the formation of intrusive bodies in the crust. These intrusive bodies, referred to here as plutons, constitute the elementary building brick of the continental crust. An increasing number of geophyscial and geochronological data as well as geological observations is currently modifying our understanding of pluton construction. When plutons were initially envisaged as quasi-spherical bodies growing slowly and essentially by an overall inflation, they are now recognized as growing incrementally by the accretion of successive and relatively small magma pulses, over variable periods of time, from hundred to millions of years, depending on geodynamic setting and source fertility. This brings new challenges and has far-reaching implications. According to this new model of plutonism, the evolution of magma bodies is related to the processes that control the timescale and the spatial distribution of the successive pulses. Depending on their emplacement rate and on their ability to amalgamate, repeated magma pulses can either rapidly solidify or ultimately build up an active magma chamber. Thus understanding how magma bodies grow has fundamental implications for the link between volcanism and plutonism as well as for magma differentiation and ultimately for our understanding of the growth and evolution of the Earth crusts. The concept of pluton incremental growth challenges our understanding as well as our field interpretations of the processes involved during pluton construction. Indeed, these processes and how they operate in governing the emplacement and growth of plutons, both in space and time, are still debated.

A state-of-the-art session on these very issues was held at the 2008 General Assembly of the European Geosciences Union. As a follow-up of this session, this special volume Emplacement of magma pulses and growth of magma bodies brings together both theoretical models and field studies that cover most aspects of the emplacement and growth of plutons.

This volume

The volume starts with three reviews of the thermal, mechanical and structural aspects of the emplacement and construction of crustal magma bodies. Annen (this issue) reviews the thermal evolution of magma bodies growing by slow amalgamation of discrete pulses, and suggests that incremental growth could explain the bi-modal character of magmatic provinces with intermediate magmas resulting from the open nature of the systems through the mixing of both mafic ans silicic end-members. This review also supports a model of magma bodies construction as a multi-time-scale process with the development of large magma chambers corresponding to the

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

highest magma fluxes, several orders of magnitude higher than the average pluton emplacement rate. Menand (this issue) then reviews the different mechanisms controlling the formation and emplacement depths of sills, and how these can amalgamate to form larger igneous bodies. The main conclusion of this review is that the dominant control on sill emplacement is exerted by crustal heterogeneities via their rigidity and rheological anisotropy, and their interaction with the local and tectonic stress fields. In turn, this determines whether and how successive sills can amalgamate to form larger magma bodies. Finally, Saint Blanquat et al. (this issue) review the mechanisms and duration of pluton construction in continental magmatic arcs. Field studies on plutons of various sizes (from 1 to more than 1000 km3) show that each pluton size is related to a pulsed magmatic activity with a characteristic time scale, and that each of these coupled time-space scales is related to a specific process.

High magma fluxes relative to crustal tectonic strain rates ensure that magmatic processes control the system from below. One of the main consequences is that for high fluxes or over short time scales the dynamics of the pulsed magmatism observed in continental arc plutonic systems are a direct proxy for deep lithospheric and magmatic processes. These reviews are complemented by two thermal and mechanical studies. First, Michaut and Jaupart (this issue) discuss two thermal models for the formation of magma reservoirs by small increments. The authors show that if a magma body grows by amalgamation of relatively thick individual igneous sheets, crystallization in that body proceeds at equilibrium as described by Annen (this issue). If the body grows instead by amalgamation of relatively thin sheets, the thermal evolution of the body is kinetically controlled such that the initial sheets do not crystallize completely but instead preserve a glassy residue. As the temperature of the system subsequently increases, devitrification of this glass occurs which leads crystallization to proceeds catastrophically in a positive feedback loop involving latent heat release and temperature rise. These two thermal models have different implications for the formation of magma chambers as well as the amount of melt generated. Then, Gudmundsson (this issue) discusses the mechanisms for the arrest of dykes and their deflection into sills at discontinuities in the upper crust in relation to the formation of shallow plutons and magma chambers. The emphasis is put on the role of hostrock layering and three associated mechanisms for the arrest of dykes and their deflection into sills. The Cook-Gordon mechanism corresponds to a weak contact that opens up as a result of tensile stresses induced by a nearby dyke, a mechanism that is likely to operate primarily at shallow crustal depths. A second mechanism is the development of stress barriers, which occurs within rock layers where the minimum compressive stress becomes vertical and thus opposes the advance of dykes and favours sill emplacement. In the third mechanism, adjacent layers with different elastic properties can lead to material-toughness ratios that are favourable to sill formation. This situation is most likely to arise when a rock layer is stiffer than the one underneath, such as in the case of a relatively stiff lava flow overlying comparatively a softer pyroclastic layer; a propagating dyke would then have a greater tendency to be deflected into the contact between the adjacent rock layers. The volume closes with four field-based studies that cover a range of pluton emplacement settings, sizes and localities. They also illustrate the challenges faced when interpreting field observations in the light of our improved understanding of pluton emplacement and growth. Numerous evidence suggest and support an incremental growth of plutons. Yet well-defined contacts within plutons remain commonly elusive. Miller et al. (this issue) compare and contrast the structural and zircon features of the Spirit Mountain batholith and the smaller Aztec Wash pluton in southern Nevada, USA, to explain this discrepancy. Both are proposed to have grown by magma replenishments and sheet-stacking into low-strength mush zones, but the contacts between successive sheets would have disappeared gradually owing to the continuous increase in magma content and heat during their growth. According to Miller et al. (this issue), the smaller Aztec Wash pluton retained a more diverse composition and texture because of its smaller volume and history, whereas a much more protracted history lead to a much larger and more homogenized

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

Spirit Mountain batholith except for the latest stage of its construction. This interpretation explains not only the apparent contradiction between incremental and protracted pluton growth and a seeming absence of intrusive contacts, but is also consistent with the recent thermal models exposed in this special volume on the thermal evolution of plutons growing by amalgamation of sills and the conditions for the existence of long-lived shallow crustal magma chambers (Annen this issue; Michaut and Jaupart, this issue). Allibon et al. (this issue) document a rare example of a pluton located in the oceanic crust which could be interpreted as a feeder of the Fuerteventura volcano, in Canary Islands. It is a vertically layered mafic dyke-like pluton constructed in a transtensive tectonic environment. The structural features of the pluton attest of an incremental growth by periodic injections and amalgamations of dykes similar to a sheeted-dyke complex, and the authors interpret this emplacement geometry as controlled by the extensional regional tectonic setting that prevailed during the Miocene pluton growth. The authors attribute also the observed sequence of internal differentiation of individual dykes, their subsequent compaction and the migration of interstitial melt that followed, as reflecting the interplay between the regional tectonic setting and the rate and volume of the successive dyke injections that were controlled instead by source-related processes. Paterson et al. (this issue) review field, geochronological and geochemical data from North America and Mongolia to constraint the size, geometry and duration of magma chambers resulting form incremental construction. This study stresses the importance of both intrusive rates and pulse sizes on the spatial and temporal characteristics of magma chambers, and so on the relationship between plutonism and volcanism, and thus echoes other theoretical studies presented in this volume (Annen, this issue; Michaud and Jaupart; this issue) that large, long-lived magma chambers can exist but their volumes and lifetime are strongly dependent on intrusive rates, their temporal variation and their spatial distribution. Finally, Marko and Yoshinobu (this issue) document the structural and thermal aureoles o the White Horse Pluton in eastern Nevada, USA. The authors use the structures of the contact aureole to show both spatial and temporal variations in the rheology of the host rocks associated with the repeated injections of magma that lead to the formation of the pluton, and thus to constrain the temporal and spatial construction of this pluton. The authors show evidences of both ductile and brittle deformation within the contact aureole, and that brittle deformation occurred in the late stage of the pluton construction following previous episode of ductile deformation. Additionally, the authors place geometric constraints on the shape and growth of the White Horse Pluton, and their estimates of the strain associated with the pluton growth suggest that shortening of the contact aureole could accounts for half the space created during the pluton formation. Other deformation mechanisms are thus required to account for the remaining space created during pluton emplacement. Our understanding of pluton construction is sharpening, and we hope this special volume Emplacement of magma pulses and growth of magma bodies provides a balanced and critical overview of our current knowledge of the incremental growth of plutons. With a set of both theoretical models and field studies that, we think, covers most aspects of the emplacement and growth of plutons, we also hope that this volume will be of interest for both the volcanic and plutonic communities, including theoreticians, field geologists and petrologists.

Acknowledgments

We would like to acknowledge the effort and expertise of all the reviewers who helped us for this volume: Pierre Barbey, Keith Benn, George Bergantz, Barrie Clarke, Sandy Cruden, Jo Dufek, Adelina Geyer, Agust Gudmundsson, Olivier Merle, Calvin F Miller, Philippe Olivier, Alessandro Tibaldi, Jean-Louis Vigneresse, Roberto Weinberg, Robert Wiebe. Furthermore, we would like to thank the editorial support from Mike Sandiford, Mian Liu, Andrés Villavicencio, Tim Horscroft, Frank Wang, Sheilagh Douma and Yvonne Philippo.

A C C E P T E D

M A N U S C R I P T