Francis Filbet

Shi Jin

A CLASS OF ASYMPTOTIC PRESERVING SCHEMES FOR KINETIC EQUATIONS AND RELATED PROBLEMS WITH STIFF SOURCES

In this paper, we propose a general time discrete framework to design asymptotic preserving schemes for initial value problem of the Boltzmann kinetic and related equations. Numerically solving these equations are challenging due to the nonlinear stiff collision (source) terms induced by small mean free or relaxation time. We propose to penalize the nonlinear collision term by a BGK-type relaxation term, which can be solved explicitly even if discretized implicitly in time. Moreover, the BGK-type relaxation operator helps to drive the density distribution toward the local Maxwellian, thus naturally imposes an asymptotic-preserving scheme in the Euler limit. The scheme so designed does not need any nonlinear iterative solver or the use of Wild Sum. It is uniformly stable in terms of the (possibly small) Knudsen number, and can capture the macroscopic fluid dynamic (Euler) limit even if the small scale determined by the Knudsen number is not numerically resolved. It is also consistent to the compressible Navier-Stokes equations if the viscosity and heat conductivity are numerically resolved. The method is applicable to many other related problems, such as hyperbolic systems with stiff relaxation, and high order parabolic equations. Contents 1. Introduction 1 2. Asymptotic Preserving (AP) stiff ODE solvers 5 3. Hyperbolic systems with relaxations 9 4. The Boltzmann equation 11 5. Numerical tests 15 6. Other applications: a nonlinear diffusion equation 20 7. Conclusion 23 References 25

Introduction

The Boltzmann equation describes the time evolution of the density distribution of a dilute gas of particles when the only interactions taken into account are binary elastic collisions. For space variable x ∈ Ω ⊂ R dx , particle velocity v ∈ R dv (d v ≥ 2), the Boltzmann equation reads:

(1.1) ∂f ∂t + v • ∇ x f = 1 ε Q(f ),
where f := f (t, x, v) is the time-dependent particles distribution function in the phase space. Here for simplicity, we do not study the case of Maxwell diffusion boundary condition for which boundary layers may be generated, but only consider specular or periodic boundary condition in space. The parameter ε > 0 is the dimensionless Knudsen number defined as the ratio of the mean free path over a typical length scale such as the size of the spatial domain, which measures the rarefiedness of the gas. The Boltzmann collision operator Q is a quadratic operator,

(1.2) Q(f )(v) = R dv S dv -1 B(|v -v ⋆ |, cos θ) (f ′ ⋆ f ′ -f ⋆ f ) dσ dv ⋆ ,
where we used the shorthanded notation

f = f (v), f ⋆ = f (v ⋆ ), f ′ = f (v ′ ), f ′ ⋆ = f (v ′ ⋆ )
. The velocities of the colliding pairs (v, v ⋆ ) and (v ′ , v ′ ⋆ ) are related by

       v ′ = v - 1 2 (v -v ⋆ ) -|v -v ⋆ | σ , v ′ ⋆ = v - 1 2 (v -v ⋆ ) + |v -v ⋆ | σ ,
with σ ∈ S dv -1 . The collision kernel B is a non-negative function which by physical arguments of invariance only depends on |vv ⋆ | and cos θ = u • σ (where u = (vv ⋆ )/|vv ⋆ | is the normalized relative velocity). In this work we assume that B is locally integrable and we will simply take

(1.3) B(|u|, cos θ) = C α |u| α ,
for some α ∈ (0, 1] and a constant C α > 0.

Boltzmann's collision operator has the fundamental properties of conserving mass, momentum and energy: at the formal level (1.4)

R dv Q(f ) φ(v) dv = 0, for φ(v) = 1, v, |v| 2 ,
and it satisfies the well-known Boltzmann's H theorem

- d dt R dv f log f dv = - R dv Q(f ) log(f ) dv ≥ 0.
The functionalf log f is the entropy of the solution. Boltzmann's H theorem implies that any equilibrium distribution function, i.e., any function which is a maximum of the entropy, has the form of a local Maxwellian distribution

M ρ,u,T (v) = ρ (2πT ) dv/2 exp - |u -v| 2 2T ,
where ρ, u, T are the density, macroscopic velocity and temperature of the gas, defined by

ρ = R dv f (v) dv = R dv M ρ,u,T (v), u = 1 ρ R dv v f (v) dv = 1 ρ R dv v M ρ,u,T (v) dv, (1.5) T = 1 d v ρ R dv |u -v| 2 f (v) dv = 1 d v ρ R dv |u -v| 2 M ρ,u,T (v) dv. (1.6)
Therefore, when the Knudsen number ε > 0 becomes very small, the macroscopic model, which describes the evolution of averaged quantities such as the density ρ, momentum ρ u and temperature T of the gas, by fluid dynamics equations, namely, the compressible Euler or Navier-Stokes equations, become adequate [START_REF] Bardos | Fluid dynamic limits of kinetic equations. I. Formal derivations[END_REF][START_REF] Bouchut | Kinetic equations and asymptotic theory[END_REF]. More specifically, as ε → 0, the distribution function will converge to a local Maxwellian M, and system (1.2) becomes a closed system for the 2 + d v moments. The conserved quantities satisfy the classical Euler equations of gas dynamics for a mono-atomic gas:

(1.7)

                   ∂ρ ∂t + ∇ x • ρ u = 0, ∂ρ u ∂t + ∇ x • (ρ u ⊗ u + p I) = 0, ∂E ∂t + ∇ x • ((E + p) u) = 0,
where p is the pressure, E represents the total energy

E = 1 2 ρ u 2 + d v 2 ρ T,
and I is the identity matrix. These equations constitute a system of 2 + d v equations in 3 + d v unknowns. The pressure is related to the internal energy by the constitutive relation for a polytropic gas

p = (γ -1) E - 1 2 ρ |u| 2 ,
where the polytropic constant γ = (d v + 2)/d v represents the ratio between specific heat at constant pressure and at constant volume, thus yielding p = ρ T . For small but non zero values of the Knudsen number ε, the evolution equation for the moments can be derived by the so-called Chapman-Enskog expansion [START_REF] Cercignani | The Boltzmann equation and its applications[END_REF], applied to the Boltzmann equation. This approach gives the Navier-Stokes equations as a second order approximation with respect to ε to the solution of the Boltzmann equation:

(1.8)

                   ∂ρ ε ∂t + ∇ x • ρ ε u ε = 0, ∂ρ ε u ε ∂t + ∇ x • (ρ ε u ε ⊗ u ε + p ε I) = ε ∇ x • [µ ε σ(u ε )], ∂E ε ∂t + ∇ x • (E ε + p ε ) u ε ) = ε ∇ x • (µ ε σ(u ε ) u + κ ε ∇ x T ε ) .
In these equations σ(u) denotes the strain-rate tensor given by

σ(u) = ∇ x u + (∇ x u) T - 2 d v ∇ x • u I
while the viscosity µ ε = µ(T ε ) and the thermal conductivity κ ε = κ(T ε ) are defined according to the linearized Boltzmann operator with respect to the local Maxwellian [START_REF] Bardos | Fluid dynamic limits of kinetic equations. I. Formal derivations[END_REF].

The connection between kinetic and macroscopic fluid dynamics results from two properties of the collision operator [START_REF] Bardos | Fluid dynamic limits of kinetic equations. I. Formal derivations[END_REF][START_REF] Bouchut | Kinetic equations and asymptotic theory[END_REF]:

(i) conservation properties and an entropy relation that imply that the equilibria are Maxwellian distributions for the zeroth order limit; (ii) the derivative of Q(f ) satisfies a formal Fredholm alternative with a kernel related to the conservation properties of (i). Past progress on developing robust numerical schemes for kinetic equations that also work in the fluid regimes has been guided by the fluid dynamic limit, in the framework of asymptotic-preserving (AP) scheme. As summarized by Jin [START_REF] Jin | Effcient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF], a scheme for the kinetic equation is AP if

• it preserves the discrete analogy of the Chapman-Enskog expansion, namely, it is a suitable scheme for the kinetic equation, yet, when holding the mesh size and time step fixed and letting the Knudsen number go to zero, the scheme becomes a suitable scheme for the limiting Euler equations • implicit collision terms can be implemented explicitly, or at least more efficiently than using the Newton type solvers for nonlinear algebraic systems. Comparing with a multi-physics domain decomposition type method [START_REF] Bourgat | Coupling Boltzmann and Euler equations without overlapping[END_REF][START_REF] Degond | A smooth transition model between kinetic and diffusion equations[END_REF][START_REF] Degond | A smooth transition model between kinetic and hydrodynamic equations[END_REF][START_REF] Günther | Numerical modeling of gas flows in the transition between rarefied and continuum regimes[END_REF][START_REF] Klar | Transition from kinetic theory to macroscopic fluid equations: a problem for domain decomposition and a source for new algorithm[END_REF][START_REF] Tallec | Coupling Boltzmann and Navier-Stoke s equations by half fluxes[END_REF], the AP schemes avoid the coupling of physical equations of different scales where the coupling conditions are difficult to obtain, and interface locations hard to determine. The AP schemes are based on solving one equation-the kinetic equation, and they become robust macroscopic (fluid) solvers automatically when the Knudsen number goes to zero. A generic way to prove that an AP scheme implies a numerical convergence uniformly in the Knudsen number was given by Golse-Jin-Levermore for the linear discrete-ordinate transport equation in the diffusion regime [START_REF] Golse | The Convergence of Numerical Transfer Schemes in Diffusive Regimes I: The Discrete-Ordinate Method[END_REF]. This result can be extended to essentially all AP schemes, although the specific proof is problem dependent. We refer to AP schemes for kinetic equations in the fluid dynamic or diffusive regimes [START_REF] Larsen | Asymptotic solutions of numerical transport problems in optically thick, diffusive regimes[END_REF][START_REF] Larsen | Asymptotic solutions of numerical transport problems in optically thick, diffusive regimes II[END_REF][START_REF] Coron | Numerical passage from kinetic to fluid equations[END_REF][START_REF] Caflish | Uniformly accurate schemes for hyperbolic systems with relaxation[END_REF][START_REF] Jin | Diffusive Relaxation Schemes for Discrete-Velocity Kinetic Equations[END_REF][START_REF] Jin | Uniformly accurate diffusive relaxation schemes for multiscale transport equations[END_REF][START_REF] Jin | Discretization of the multiscale semiconductor Boltzmann equation by diffusive relaxation schemes[END_REF][START_REF] Klar | An asymptotic-induced scheme for nonstationary transport equations in the diffusive limit[END_REF][START_REF] Klar | An asymptotic preserving numerical scheme for kinetic equations in the low Mach number limit[END_REF][START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF][START_REF] Bennoune | Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible NavierStokes asymptotics[END_REF][START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF]. The AP framework has also been extended in [START_REF] Crispel | An asymptotically preserving scheme for the two-fluid Euler-Poisson model in the quasi-neutral limit[END_REF][START_REF] Degond | An asymptotically stable Particle-in-Cell (PIC) scheme for collisionless plasma simulations near quasineutrality[END_REF] for the study of the quasi-neutral limit of Euler-Poisson and Vlasov-Poisson systems, and in [START_REF] Degond | Mach-number uniform asymptotic-preserving gauge schemes for compressible flows[END_REF][START_REF] Degond | All speed scheme for the low mach number limit of the Isentropic Euler equation[END_REF][START_REF] Haack | An all-speed asymptotic-preserving schemes for compressible flows[END_REF] for all-speed (Mach number) fluid equations bridging the passage from compressible flows to the incompressible flows. One should note that underresolved computation may not yield accurate or even physically correct approximations in areas with sharp transitions, such as shock and boundary layers. In these areas one may want to use resolved calculations. The AP schemes allow one to use suitable mesh size and time step at needed domains with one first-principle equation, thus is especially suitable for problems with localized sharp transitions where macroscopic simulation is necessary.

To satisfy the first condition for AP, the scheme must be driven to the local Maxwellian when ε → 0. Let t n (n = 0, 1, 2, • • • ) be the discrete time, and U n = U (t n ) for a general quantity U . Then an AP scheme requires that, for ∆t ≫ ε,

(1.9) f n -M n = O(ε) , n ≥ 1
for any initial data f 0 . Namely, the numerical solution projects any data into the local Maxwellian, with an accuracy of O(ε), in one step. This can usually be achieved by a backward Euler or any L-stable ODE solvers for the collision term [START_REF] Jin | Runge-Kutta Methods for Hyperbolic Conservation Laws with Stiff Relaxation Terms[END_REF]. Such a scheme requires an implicit collision term to guarantee a uniform stability in time. However, how to invert such an implicit, yet nonlocal and nonlinear, collision operator is a delicate numerical issue. Namely, it is hard to realize the second condition for AP schemes. One solution was offered by Gabetta, Pareschi and Toscani [START_REF] Gabetta | Relaxation schemes for nonlinear kinetic equations[END_REF]. They first penalize Q by a linear function λf , and then absorb the linearly stiff part into the time variable to remove the stiffness. The remaining implicit nonlinear collision term is approximated by finite terms in the Wild Sum, with the infinite sum replaced by the local Maxwellian. This yields a uniformly stable AP scheme for the collision term, capturing the Euler limit when ε → 0. Such a time-relaxed method was also used to develop AP Monte Carlo method, see [START_REF] Caflisch | An implicit Monte Carlo method for rarefied gas dynamics I: The space homogeneous case[END_REF][START_REF] Pareschi | Time relaxed Monte Carlo methods for the Boltzmann equation[END_REF]. When the collision operator Q is the BGK collision operator (1.10)

Q BGK = M -f ,
it is well-known that even an implicit collision term can be solved explicitly, using the property that Q preserves mass, momentum and energy [START_REF] Coron | Numerical passage from kinetic to fluid equations[END_REF]. Our new idea in this paper is to utilize this property, and penalize the Boltzmann collision operator Q by the BGK operator:

(1.11) Q = [Q -λ(M -f )] + λ[M -f ]
where λ is the spectral radius of the linearized collision operator of Q around the local Maxwellian M. Now the first term on the right hand side of (1.11) is either not stiff, or less stiff and less dissipative compared to the second term, thus it can be discretized explicitly, so as to avoid inverting the nonlinear operator Q. The second term on the right hand side of (1.11) is stiff or dissipative, thus will be treated implicitly. As mentioned earlier, the implicit BGK operator can be inverted explicitly. Therefore we arrive at a scheme which is uniformly stable in ε, with an implicit source term that can be solved explicitly. In other words, in terms of handling the stiffness, the general Boltzmann collision operator can be handled as easily as the much simpler BGK operator, thus we significantly simplify an implicit Boltzmann solver! A related problem is hyperbolic systems with relaxations. Such systems arise in reacting gases, shallow water equations, discrete-velocity kinetic models, etc. [START_REF] Whitham | Linear and Nonlinear Waves[END_REF], and have been mathematically studied extensively in recent years (see for example [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF][START_REF] Jin | The Relaxation Schemes for Systems of Conservation Laws in Arbitrary Space Dimensions[END_REF][START_REF] Bianchini | Hyperbolic limit of the Jin-Xin relaxation model[END_REF][START_REF] Natalini | Recent results on hyperbolic relaxation problems[END_REF]). A prototype example is the following 2 × 2 nonlinear hyperbolic system with relaxation:

(1.12)

     ∂u ∂t + f 1 (u, v) x = 0 , ∂v ∂t + f 2 (u, v) x = 1 ε R(u, v) .
The relaxation term R : R 2 → R is dissipative in the sense of [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF]:

(1.13) ∂ v R ≤ 0 .
It possesses a unique local equilibrium, namely, R(u, v) = 0 implies v = g(u). At the local equilibrium, one has the macroscopic system u t + f 1 (u, g(u)) x = 0 . This system can be derived by sending ε → 0 in (1.12), the so-called zero relaxation limit ( [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF]). This limit is analogous to the passage from kinetic equations to their fluid limit, and in the last decade the development in these two areas -both analytic studies and numerical approximations-have strongly intervened. The numerical methods for such systems are similar to those developed for the Boltzmann equations, especially for discrete-velocity kinetic models [START_REF] Caflish | Uniformly accurate schemes for hyperbolic systems with relaxation[END_REF][START_REF] Jin | Runge-Kutta Methods for Hyperbolic Conservation Laws with Stiff Relaxation Terms[END_REF]. The guiding principle for the AP schemes is the same for both classes of problems, and in this paper we will study both applications whenever appropriate.

Let V n and U n be the time-discrete approximations to v and u respectively in (1.12). A classical AP scheme requires that, for ∆t ≫ ε,

(1.14) V n -g(U n ) = O(ε) , n ≥ 1
for any initial data V 0 . Namely, the numerical solution projects any data V into the local equilibrium V = f (U ), with an accuracy of O(ε), in one step. This is the analogy of (1.9). Our new method is not necessarily AP in the classical sense of (1.14). Nevertheless, we can show that, for any ε, and ∆t ≫ ε, there exists an

N ε ≥ 1 , such that (1.15) V n -g(U n ) = O(ε) , n ≥ N ε
for any initial data V 0 . Namely, the numerical solution projects the solution into the local equilibrium after the initial transient time, for any initial data. This is a slightly weaker condition than (1.14), but is enough to guarantee the desired numerical performance as good as the classical AP schemes. Although a linear penalty (by removing M on the right hand side of (1.11)) can also remove the stiffness, we can show that, when applied to the relaxation system (1.12), it only has the following property:

(1.16) V n -g(U n ) = O(∆t) , n ≥ N
for any initial data V 0 , when ∆t ≫ ε. Since in the fluid regime, we really want to take ∆t ≫ ε, schemes with a weak AP property (1.16) is much less accurate than our scheme which has the property (1.15).

The BGK operator that we use in (1.11) helps to drive f into M (or V into g(u)) more effectively than a linear damping -λf , thus preserves the Euler limit more accurately. Moreover, if v ng(U n ) = O(ε) (well-prepared initial data), then our method implies that v n+1g(U n+1 ) = O(ε), while the linear penalty method always yields v n+1g(U n+1 ) = O(∆t) even for well-prepared initial data.

For the Boltzmann equation, although we cannot analytically prove an analogy of (1.15) for f -M, our numerical examples show that this is true. We can prove, however, that if the initial data are well-prepared,

f n -M n = O(ε) for some n = N ≥ 0 .
then the scheme captures the correct Euler limit for later time n > N . Moreover, for suitably small time-step, our method is also consistent to the Navier-Stokes equations (1.8) for ε ≪ 1.

Our method is partly motivated by the work of Haack, Jin and Liu [START_REF] Haack | An all-speed asymptotic-preserving schemes for compressible flows[END_REF], where by subtracting the leading linear part of the pressure in the compressible Euler equations with a low Mach number, the nonlinear stiffness in the pressure term due to the low Mach number is removed and an AP scheme was proposed for the compressible Euler or Navier-Stokes equations that capture the incompressible Euler or Navier-Stokes limit when the Mach number goes to zero. In terms of removing the stiffness of nonlinear parabolic equations Smereka used the idea of adding and subtracting a linear elliptic operator. However his approach was not aimed at achieving the AP property.

Our method is not restricted to the Boltzmann equation. It applies to general nonlinear hyperbolic systems with stiff nonlinear relaxation terms [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF][START_REF] Jin | Numerical schemes for hyperbolic conservation laws with stiff relaxation terms[END_REF][START_REF] Jin | Runge-Kutta Methods for Hyperbolic Conservation Laws with Stiff Relaxation Terms[END_REF][START_REF] Coquel | Relaxation of energy and approximate Riemann solvers for general pressure laws in fluid dynamics[END_REF]], as will be shown in Section 3, and higherorder parabolic equations (see Section 6). Indeed, it applies to any stiff source term that admits a stable local equilibrium.

We will present and study this framework for stiff ODEs (Section 2), nonlinear hyperbolic system with relaxation (Section 3), and the Boltzmann equation (Section 4). We present different numerical tests on the Boltzmann equation in Section 5 to illustrate the efficiency of the present method. In particular, we will include a multi-scale problem where the Knudsen number ε depends on the space variable and takes different values ranging from 10 -4 (hydrodynamic regime) to 1 (kinetic regime). Finally, in Section 6, we design a scheme for the nonlinear Fokker-Planck equations for which the asymptotic preserving scheme can be used to remove the CFL constraint of a parabolic equation. We conclude the paper in Section 7.

Asymptotic Preserving (AP) stiff ODE solvers

We first present our method for stiff ordinary differential equations. Let us consider a Hilbert space H and the following nonlinear autonomous ordinary differential system (2.1)

     df ε dt (t) = Q(f ε ) ε , t ≥ 0, f ε (0) = f 0 ∈ H,
where the source term Q(f ) satisfies the following properties:

• there exists a unique stationary solution M to (2.1), namely, Q(M) = 0;

• the solution to (2.1) converges to the steady state M when time goes to infinity, and the spectrum of ∇Q(f

) ⊂ C -= {z | z ∈ C -, Re(z) < 0}, (2.2) 0 < α m ≤ ∇Q(f ) ≤ α M , ∀ f ∈ H \ {0}.
where ∇Q(f ) denotes the Frechet derivative of Q.

Remark 2.1. The second hypothesis above is certainly not the most general, but is convenient for our purpose. The lower bound implies that the solution converges to the steady state M, while the upper bound is a sufficient condition for existence and uniqueness of a global solution.

When ε becomes small, the differential equation (2.1) becomes stiff and explicit schemes are subject to severe stability constraints. Of course, implicit schemes allow larger time step, but new difficulty arises in seeking the numerical solution of a fully nonlinear problem at each time step. Here we want to combine both advantages of implicit and explicit schemes : large time step for stiff problems and low computational complexity of the numerical solution at each time step.

Two classical procedures handle the aforementioned difficulties well. One is to linearize the unknown Q(f n+1 ) at time step t n+1 around f at the previous time step f n :

(2.3) Q(f n+1 ) ≈ Q(f n ) + ∇Q(f n )(f n+1 -f n ) .
This yields a problem that only needs to solve a linear system with coefficient matrices depending on ∇Q(f n ) [START_REF] Yee | A Class of High-Resolution Explicit and Implicit Shock-Capturing Methods[END_REF]. This approach gives a uniformly stable time discretization without nonlinear solvers.

The second approach, introduced in [29], takes

(2.4) Q(f ) = [Q(f ) -µf ] + µf .
In [START_REF] Jin | Uniformly accurate diffusive relaxation schemes for multiscale transport equations[END_REF], the second µf term in absorbed into the time derivative, which removes the stiffness, and then Q(f ) is approximated by the Wild sum which is truncated at finite terms with the remaining infinite series replaced by the local Maxwellian in order to become AP. If one is just interested in removing the stiffness, one can just approximate the right hand side of (2.4) by

[Q(f n ) -µf n ] + µf n+1 .
For sufficiently large µ, this yields a scheme with stability independent of ε, yet can be solved explicitly. However, a disadvantage of the linear penalty method, as well as method (2.3), is that the operators on the right hand size do not preserve exactly the mass, momentum and total energy as the BGK operator does.

As will be shown in Section 3, these two classical approaches project the data into the local equilibrium in the sense of (1.16).

We propose to split the source term of (2.1) as the sum of a stiff-dissipative part and a non-(or less) stiff and non-dissipative part as

(2.5) Q(f ) ε = Q(f ) -P (f ) ε less stiff part + P (f ) ε stiff, dissipative part
, where P (f ) is a well balanced, i.e. preserving the steady state, P (M) = 0, linear operator and is asymptotically close to the source term Q(f ). For instance, performing a simple Taylor expansion, we get

Q(f ) = Q(M) + ∇Q(M) (f -M) + O( f -M 2
H ) and we may choose

P (f ) := ∇Q(M) (f -M
). Since it is not always possible to compute exactly ∇Q(M), we may simply choose

P (f ) := β (M -f ) ,
where β is an upper bound of ∇Q(M) or some approximation of it such as

[Q(f )-Q(M)]/(f -M).
In the following we propose a discretization to (2.5) based on IMEX schemes. We simply apply a first order implicit-explicit (IMEX) scheme for the time discretization of (2.1):

(2.6)

f n+1 -f n ∆t = Q(f n ) -P (f n ) ε + P (f n+1 ) ε , or f n+1 = [ε I -∆t ∇Q(M)] -1 [ε f n + ∆t (Q(f n ) -P (f n )) -∆t ∇Q(M) M] .
This method is easy to implement, since f n+1 is linear in the right hand side of (2.6). For linear problems, we have the following result: Theorem 2.2. Consider the differential system (2.1) with Q(f ) = -λ f , where Re(λ) > 0. Set P (f ) := -ν λ f with ν ≥ 0. Then, the scheme (2.6) is A-stable and L-stable for ν > 1/2.

Proof. For linear systems with Q(f ) = -λf , the scheme simple reads

f n+1 = ε + (ν -1) λ ∆t ε + νλ ∆t f n = 1 - λ∆t ε + ν λ∆t f n .
Observe that ν = 0 gives the explicit Euler scheme, which is stable only for ∆t ≤ ε/λ, whereas for 0 ≤ ν ≤ 1, it yields the so-called θ-scheme, which is A-stable for ν > 1/2. For ν = 1 it corresponds to the A-stable implicit Euler scheme. Moreover,

f n+1 H ≤ 1 - λ∆t ε + ν λ∆t f n H ∼ 1 - 1 ν f n H for ε ∼ 0 or λ∆t ≫ 1 , where |1 -1 ν | < 1 for ν > 1/2
. This is also the condition for the L-stability [START_REF] Gear | Numerical Initial Value Problems in Ordinary Differential Equations[END_REF]. Clearly λ ∼ 1 gives the fastest convergence to the equilibrium.

Concerning nonlinear problems, we observe that the scheme (2.6) is not AP in the sense of (1.9). However, we can prove that it is AP in the sense of (1.15). Theorem 2.3. Assume that the operator Q satisfies (2.2) and

(2.7) Q(f n ) -Q(M) f n -M < 0.
Assume that ∆t ≫ ε. Then, for β sufficiently large, there exists an 0 < r < 1, independent of ∆t and ε, such that

|f n -M| ≤ r n |f 0 -M| . Consequently scheme (2.6) is AP in the sense of (1.15). Proof. We choose β > 0 such that β > 1 2 sup f ∈H Q(f ) -Q(M) f -M = 1 2 α M .
Scheme (2.6) can be written as

[f n+1 -M] -[f n -M] ∆t = 1 ε Q(f n ) -Q(M) f n -M + β (f n -M) - β [f n+1 -M] ε .
This gives

1 + β ∆t ε [f n+1 -M] = 1 + ∆t ε D n [f n -M],
where D n is given by

D n = Q(f n ) -Q(M) f n -M + β.
Clearly, under the assumption (2.7),

|D n | ≤ β - 1 2 α M .
Thus,

f n+1 -M = ε + ∆t D n ε + ∆t β [f n -M n ]. For ∆t ≫ ε, r = sup ε,n,∆t ε + ∆t D n ε + ∆t β ∼ sup n β -1 2 α M β < 1 hence (2.8) implies |f n+1 -M| ≤ r |f n -M|. From here it is simple to see that |f n -M| ≤ r n |f 0 -M|.
So for any ε > 0, and any initial data, there exists an N ε ≥ 1 such that when n ≥ N ε , f n -M = O(ε). This is the AP property defined in (1.15).

To improve the numerical accuracy, second order schemes are sometimes more desirable. Thus, we propose the following second order IMEX extension. Assume that an approximate solution f n is known at time t n , we compute a first approximation at time t * using a first order IMEX scheme and next apply the trapezoidal rule and the mid-point formula. The scheme reads (2.8)

         2 f ⋆ -f n ∆t = Q(f n ) -P (f n ) ε + P (f ⋆ ) ε . f n+1 -f n ∆t = Q(f ⋆ ) -P (f ⋆ ) ε + P (f n ) + P (f n+1 ) 2 ε .
For Q = -λf, P = -νλf , (2.8) gives

f n+1 = 1 + ∆t ε λ(ν -1) + 1 4 ∆t ε 2 λ 2 (ν 2 -4ν + 2) 1 + ∆t 2ε νλ 2 f n .
For ∆t ≫ ε this gives

f n+1 ∼ ν 2 -4ν + 2 ν 2 f n Note that r = ν 2 -4ν + 2 ν 2 < 1 if ν > 1 2 ,
thus the second order IMEX scheme has the same AP property as the first order scheme (2.6). Moreover, we can prove a theorem similar to Theorem 2.3 for (2.8) but the details are omitted here.

To illustrate the efficiency of (2.6) and (2.8) in various situations, we consider a simple linear problem with different scales for which only some components rapidly converge to a steady state whereas the remaining part oscillates. We solve (2.9)

Q(f ) = A f,
where

A =   -1000 1 0 -1 -1000 0 0 0 i   (2.10)
for which the eigenvalues are Sp(A) = {-1000 + i, -1000i, i}. The first block represents the fast scales whereas the last one is the oscillating part. Indeed, the first components go to zero exponentially fast whereas the third one oscillates with respect to time with a period of 2π. We want to solve accurately the oscillating part with a large time step without resolving the small scales. Then, we apply the first order (2.6) and second order (2.8) schemes by choosing

P (f ) = ν A f,
with ν ≥ 0. Here we take a large time step ∆t = 0.3 and ν = 2, which means that P (f ) has the same structure of Q(f ) but the eigenvalues are over estimated. Thus, fast scales are under-resolved whereas this times step is a good discretization of the third oscillating component. Therefore, an efficient AP scheme would give an accurate behavior of the slow oscillating scale with large time step with respect to the fast scale.

In Figure 1, we present the real part of the numerical solution to the differential system (2.9)-(2.10) corresponding to the initial datum f (0) = (2, 1, 1) on the time interval [0, 15]. We compare the numerical solution obtained with our first (2.6) and second (2.8) order AP schemes using a large time step (∆t = 0.3) and the one obtained with a first and second order explicit Runge-Kutta scheme using a small time step (∆t = 0.0001) for which the numerical solution is stable.

It clearly appears in Figure 1 (1) that the time step is too large to give accurate results for the first order scheme (2.6): the solution is stable but the oscillation of the third component is damped for this time step which is too large. This approximation is compared with the one obtained with a first order explicit Euler using a times step 300 times smaller. Thus, the first order AP scheme gives a numerical solution which is stable for large time step but the accuracy is not satisfying.

Therefore, we also compare the numerical solution of the second order scheme (2.8) with the one obtained using a second order explicit Runge-Kutta scheme corresponding to ν = 0 with a time step three hundred times smaller. In Figure 1, we observe the stability and good accuracy of the second order scheme (2.8). Let us emphasize that for the same time step, the numerical solution given by an explicit Runge-Kutta scheme blows-up (hence the result is not reported here)! t=0.001 (1)

∆ t f(t) X f(t) --- ∆ =0.3 1st
Figure 1. Comparison of the time evolution of the numerical approximation to the differential system (2.9)-(2.10) with f (0) = (2, 1, 1). ( 1) first and (2) second order Asymptotic Preserving and explicit Runge-Kutta schemes

Hyperbolic systems with relaxations

In this section, we propose and study the method for hyperbolic system with (stiff) relaxations. We propose the following temporal approximation to (1.12):

(3.1)          U n+1 -U n ∆t + f 1 (U n , V n ) x = 0 , V n+1 -V n ∆t + f 2 (U n , V n ) x = 1 ε [R(U n , V n ) + β(V n -g(U n ))] - β ε [V n+1 -g(U n+1 )] .
Assume all functions are smooth. Some simple mathematical manipulations on (3.1) give

V n+1 -g(U n+1 ) = -f 2 (U n , V n ) x + (g(U n+1 ) -g(U n ))/∆t ε∆t ε + β∆t + 1 + ∆t ε β + R(U n ,V n ) V n -g(U n ) 1 + β ∆t ε (V n -g(U n )) . (3.2) Note that R(U n , V n ) V n -g(U n ) = R(U n , V n ) -R(U n , g(U n )) V n -g(U n ) = ∂ v R(U n , W n ) < 0 for some W n , thus if β > 1 2 sup |∂ v R| ,
there exists a constant C, and 0 < r < 1 such that

|V n+1 -g(U n+1 )| ≤ C ε∆t ε + β∆t + r|V n -g(U n )| .
From here it is easy to see that

|V n -g(U n )| ≤ C 1 -r ε∆t ε + β∆t + r n |V 0 -g(U 0 )| This clearly gives (3.3) |V n -g(U n )| ≤ C (1 -r)β ε + r n |V 0 -g(U 0 )|
in which the first term on the right hand side is O(ε) independent of ∆t. For any ε ≪ 1, there exists an

N ε ≥ 1 such that r n |V 0 -g(U 0 )| ≤ ε , therefore (3.
3) implies the desired AP property (1.15).

Next we consider the linear penalty method (2.4):

(3.4)

         U n+1 -U n ∆t + f 1 (U n , V n ) x = 0 , V n+1 -V n ∆t + f 2 (U n , V n ) x = 1 ε [R(U n , V n ) + βV n ] - β ε V n+1 .
A simple mathematical manipulation on (3.4) gives

V n+1 -g(U n+1 ) = -f 2 (U n , V n ) x ε∆t ε + β∆t -[g(U n+1 ) -g(U n )] + 1 + ∆t ε µ + R(U n ,V n ) V n -g(U n ) 1 + µ ∆t ε (V n -g(U n )) . (3.5)
The first two terms on the right hand side of (3.5) can only be bounded by C(ε + ∆t), while the third term, under the condition

µ > 1 2 sup |∂ v R| ,
is similar to the second term on the right hand side of (3.2). In conclusion, corresponding to (3.3), here we can only obtain

(3.6) |V n -g(U n )| ≤ C(ε + ∆t) + r n |V 0 -g(U 0 )| which, if ∆t ≫ ε, gives only (1.16).
Another observation is the following. From (3.3), one sees that for prepared initial data

(3.7) V 0 = g(U 0 ) + O(ε) (3.3) implies that V n = g(U n ) + O(ε)
, for any n ≥ 1 Namely, if the data are within O(ε) of the local equilibrium, they remain so for all future times. However, for the linear penalty method, even if the initial data are well prepared as in (3.7), from (3.6) one sees that V 1 = g(U 1 ) + O(ε + ∆t) , so the deviation from the local equilibrium at later times is always of O(∆t) rather than O(ε). A similar analysis on method (2.3) gives a result as in (3.6). We omit the details here. Now to illustrate the effciency of our approach, we present numerical simulations on (1.12). We simply consider

(3.8)                  ∂u ∂t + ∂v ∂x = 0, ∂v ∂t + a ∂u ∂x = g(v) ε (f (u) -v) , u(t = 0, x) = 1 + 0.9 sin(π x), v(t = 0, x) = cos(π x), x ∈ (-1, 1), with g(v) = 1 + |v| 4 , f (u) = u 2 /2 and a = sup u |f ′ (u)| 2 .
In Figure 2, we represent the approximation of the solution at time t = 0.1 for different values of ε = 10 -1 ; 10 -2 and 10 -7 obtained with our AP scheme (3.1) and the linear penalty method (3.4). The number of points in space is n x = 800 and ∆t = 0.0006. Clearly, for the same time step, our scheme gives the correct behavior which corresponds to the well known solution to the Burgers equations when ε tends to zero, whereas the linear penalty method gives a stable approximation which is not accurate. Of course, when ∆t becomes smaller and ε is fixed, the linear penalty method is accurate. Here, the initial data is far from the equilibrium hence the linear penalty method is not appropriate since it does not have any mechanism of projection to the steady state when ε is small.

The Boltzmann equation

We now extend the method to the Boltzmann equation (1.1). To this aim, we rewrite the Boltzmann equation (1.1) in the following form (4.1)

     ∂f ∂t + v ∇ x f = Q(f ) -P (f ) ε + P (f ) ε , x ∈ Ω ⊂ R dx , v ∈ R dv , f (0, x, v) = f 0 (x, v), x ∈ Ω, v ∈ R dv ,
where the operator P is a "well balanced relaxation approximation" of Q(f ), which means that it satisfies the following (balance law)

R dv P (f ) φ(v) dv = 0, φ(v) = 1, v, |v| 2 ,
and preserves the steady state i.e. P (M ρ,u,T ) = 0 where M ρ,u,T is the Maxwellian distribution associated to ρ, u and T given by (1.5). Moreover, it is a relaxation operator in velocity (4.2)

P (f ) = β [M ρ,u,T (v) -f (v)] .
4.1. Choice of the free parameter β. For instance, P (f ) can be computed from an expansion of the Boltzmann operator with respect to M ρ,u,T :

Q(f ) ≃ Q(M ρ,u,T ) + ∇Q(M ρ,u,T ) [M ρ,u,T -f ] .
Thus, we choose β > 0 as an upper bound of the operator ∇Q(M ρ,u,T ). Other choices of β are also possible, for example

β = sup Q(f ) -Q(M) f -M = sup Q(f ) f -M ,
or, at time t n ,

β n = sup Q(f n ) -Q(f n-1 ) f n -f n-1 .
Then P (f ) given by (4.2) is just the BGK collisional operator [START_REF] Bhatnagar | A model for collision processes in gases[END_REF]. One can also choose β such that the operator P (f ) gives the same viscosity (of order to ε) as Q(f ) when applying a Chapman Enskog expansion. 4.2. Discretization to the Boltzmann equation. Since the convection term in (4.1) is not stiff, we will treat it explicitly. The source terms on the right hand side of (4.1) will be handled using the ODE solver in the previous section. For example, if the first order scheme (2.6) is used, then we have

(4.3)      f n+1 -f n ∆t + v • ∇ x f n = Q(f n ) -P (f n ) ε + P (f n+1 ) ε , f 0 (x, v) = f 0 (x, v) .
Using the relaxation structure of P (f ) given in (4.2), it can be written as

f n+1 = ε ε + β n+1 ∆t [f n -∆t v ∇ x f n ] + ∆t Q(f n ) -P (f n ) ε + β n+1 ∆t (4.4) + β n+1 ∆t ε + β n+1 ∆t M n+1 ,
where

β n+1 = β(ρ n+1 , T n+1
) and M n+1 is the Maxwellian distribution computed from f n+1 . Although (4.4) appears nonlinearly implicit, it can be computed explicitly. Specifically, upon multiplying (4.4) by φ(v) defined in (1.4), and use the conservation property of Q and P and the definition of M in (1.5), we define the macroscopic quantity U by U := (ρ, ρ u, T ) computed from f and get [START_REF] Coron | Numerical passage from kinetic to fluid equations[END_REF][START_REF] Puppo | Implicit-explicit schemes for BGK kinetic equations[END_REF]]

U n+1 = ε ε + β n+1 ∆t φ(v) (f n -∆t v • ∇ x f n ) dv + β n+1 ∆t ε + β n+1 ∆t U n+1 ,
or simply

U n+1 = φ(v) (f n -∆tv • ∇ x f n ) dv .
Thus U n+1 can be obtained explicitly, which defines M n+1 . Now f n+1 can be obtained from (4.4) explicitly. In summary, although (4.3) is nonlinearly implicit, it can be solved explicitly, thus satisfies the second condition of an AP scheme.

Clearly, scheme (4.3) satisfies the following properties Proposition 4.1. Consider the numerical solution given by (4.3). Then,

(i) If ε → 0 and f n = M n + O(ε), then f n+1 = M n+1 + O(ε)
. Thus, when ε → 0, the (moments of the) scheme becomes a consistent discretization of the Euler system (1.7). (ii) Assume ε ≪ 1 and f n = M n + ε g n . If there exists a constant C > 0 such that

(4.5) g n+1 -g n ∆t + ∇ x (vg n ) + g n ≤ C, and
(4.6) U n + U n+1 -U n ∆t ≤ C,
then the scheme (4.3) asymptotically becomes a first order in time approximation of the compressible Navier-Stokes (1.8).

Proof. We easily first check that for ε → 0 and

f n = M n + O(ε), we get f n+1 = M n+1 + O(ε).
Therefore, we multiply (4.3) by (1, v, |v| 2 /2) and integrate with respect to v, which yields that U n is given by a time explicit scheme of the Euler system (1.7). Now let us prove (ii). We apply the classical Chapman-Enskog expansion:

(4.7)

f n = M n + ε g n
and integrate (4.3) with respect to v ∈ R dv . By using the conservation properties of the Boltzmann operator (1.4) and of the well-balanced approximation P (f ),

(4.8) U n+1 -U n ∆t + ∇ v • R dv    1 v |v| 2 2    v (M n + εg n )dv = 0.
For εg = 0, this is the compressible Euler equations (1.7). Thus, a consistent approximation of the compressible Navier-Stokes is directly related to a consistent approximation of g n . Inserting decomposition (4.7) into the scheme (4.3) gives

M n+1 -M n ∆t + v ∇ x M n + ε g n+1 -g n ∆t + v ∇ x g n = Q(M n + εg n ) ε -β n g n -β n+1 g n+1 ,
Since Q is bilinear and Q(M) = 0, one has

Q(M + ε g) = Q(M) + ε L M (g) + ε 2 Q(g),
where L M is the linearized collision operator with respect to M. Thus, we get

M n+1 -M n ∆t -β n g n -β n+1 g n+1 + ε g n+1 -g n ∆t + v ∇ x g n -Q(g n ) = L M (g n ) -v ∇ x M n , (4.9)
It is well known that L M is a non-positive self-adjoint operator on L 2 M defined by the set

L 2 M := {ϕ : ϕ M -1/2 ∈ L 2 (R dv )} and that its kernel is N (L M ) = Span{M, v M, |v| 2 M}. Let Π M be the orthogonal projection in L 2 M onto N (L M ).
After easy computations in the orthogonal basis, one finds that [START_REF] Bouchut | Kinetic equations and asymptotic theory[END_REF]

Π M (ψ) = M ρ m 0 + v -u T m 1 + |v -u| 2 2T - d v 2 m 2
where

m 0 = R dv ψ dv, m 1 = R dv (v -u)ψ dv, m 2 = R dv |v -u| 2 2T - d v 2 ψ dv.
It is easy to verify that Π M n (M n ) = M n and

Π M n (g n ) = Π M n (g n+1 ) = Π M n (Q(g n )) = Π M n (L M n (g n )) = 0.
Then applying the orthogonal projection I -Π M n to (4.9), it yields

(I -Π M n ) M n+1 -M n ∆t -β n g n -β n+1 g n+1 + ε g n+1 -g n ∆t + (I -Π M n ) (v ∇ x g n ) -Q(g n ) = L M (g n ) -(I -Π M n ) (v ∇ x M n ).
Using the assumption (4.5) we get that the term

ε g n+1 -g n ∆t + (I -Π M n ) (v ∇ x g n ) -Q(g n )
is of order ε. Then, it remains to estimate the terms β n+1 g n+1β n g n and

(I -Π M n ) M n+1 -M n ∆t .
First, we have

β n+1 g n+1 -β n g n = β n+1 (g n+1 -g n ) + (β n+1 -β n ) g n .
Under the assumption (4.5) and (4.6), and since β n only depends on U n , we easily get (4.10)

β n+1 g n+1 -β n g n = O(∆t).
Next using a Taylor expansion we find that

M n+1 = M n 1 + ρ n+1 -ρ n ρ n + v -u n T n u n+1 -u n + |v -u n | 2 2T n - d 2 T n+1 -T n T n + O(∆t 2 )
and by definition of

Π M Π M n (M n+1 ) = M n 1 + ρ n+1 -ρ n ρ n + v -u n T n u n+1 -u n + |v -u n | 2 2T n - d 2 T n+1 -T n T n + M n |v -u n | 2 2T n - d 2 T n+1 -T n ρ n T n (ρ n+1 -ρ n ) + ρ n+1 dρ n T n (u n+1 -u n ) 2 + M n v -u n T n ρ n+1 -ρ n ρ n u n+1 -u n + O(∆t 2 ).
Thus, under assumption (4.6), we have

(4.11) (I -Π M n ) M n+1 -M n ∆t = O(∆t).
Gathering (4.10) and (4.11), the residual distribution function is given by

g n = L -1 M n (I -Π M n ) (v • ∇ x M n ) + O(ε) + O(∆t)
. Now, substituting this latter expression in (4.8), we get

U n+1 -U n ∆t + ∇ x • F (U ) = -ε ∇ x • R dv    v v ⊗ v v |v| 2 2    L -1 M n (Id -Π M n ) (v • ∇ x M n ) dv (4.12) + O ε∆t + ε 2 , (4.13)
where

F (U ) =   ρ u ρ u ⊗ u + p I (E + p) u   .
To complete the proof, it remains to compute the term in O(ε). After some computations, we first get

(I -Π M n ) (v • ∇ x M n ) = B ∇u + (∇u) T - d 2 ∇ • u I + A ∇T √ T M(v), with A = |v -u| 2 2T - d + 2 2 v -u √ T , B = 1 2 (v -u) ⊗ (v -u) 2T - |v -u| 2 dT I .
Therefore, it yields

L -1 M n (I -Π M n ) (v • ∇ x M n ) = L -1 M n (B M) ∇u + (∇u) T - d 2 ∇ • u I +L -1 M n (A M) ∇T √ T .
Substituting this expression in (4.8), we get a consistent time discretization scheme to the compressible Navier-Stokes system where the term of order of ε is given by

ε ∇ x •   0 µ ε σ(u ε ) µ ε σ(u ε ) u + κ ε ∇ x T ε   with σ(u) = ∇ x u + (∇ x u) T - 2 d v ∇ x • u I
while the viscosity µ ε = µ(T ε ) and the thermal conductivity κ ε = κ(T ε ) are defined according to the linearized Boltzmann operator with respect to the local Maxwellian [START_REF] Bardos | Fluid dynamic limits of kinetic equations. I. Formal derivations[END_REF].

At this stage, let us address several comments concerning Proposition 4.1.

• Note that the assumption (4.5) is very difficult to prove for our scheme. However, a similar assumption is done in [START_REF] Bennoune | Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible NavierStokes asymptotics[END_REF].

• We only prove theoretical results when the initial data is close enough to the local Maxwellian. However, it is expected that, as it is shown in Section 3, after the initial transient time the solution is only O(ε) distance from the local equilibrium. While this has not been proven for the Boltzmann equation (since it does not have a property similar to (1.13), our numerical results in Section 4 strongly suggest so. A rigorous proof remains an open question. • Under-resolved computations using AP schemes can only capture the solutions of the Euler equations. To capture the Navier-Stokes approximation that has O(ε) viscosity and heat conductivity, one needs the mesh size and c∆t to be o(ε) (c is a characteristic speed). Thus conclusion (ii) in the above proposition shows that the scheme is consistent to the Navier-Stokes equations provided that the viscous terms are resolved. In other words, one cannot expect to capture the Navier-Stokes solution with under-resolved ( ∆x c , ∆t ≫ ε) mesh sizes and time steps. On the other hand, if one has to resolve the viscous term using ∆x c , ∆t = o(ε) it will be more efficient to directly solve the Boltzmann equation directly. Thus we do not advocate an AP scheme for the compressible Navier-Stokes limit. Nevertheless, the result of Proposition 4.1 (ii) is still analytically interesting. If one directly compares the error of numerical solutions f with the solution of the Boltzmann equation by, say a first order method, one usually arrives at an error of O(∆t/ε) (see a related study in ( [START_REF] Golse | The Convergence of Numerical Transfer Schemes in Diffusive Regimes I: The Discrete-Ordinate Method[END_REF])), but if compared with the solutions of the Navier-Stokes equation, which are moments of f , (4.12) shows that the error is of order O(∆t + ε∆t). Here O(∆t) comes from the Euler time discretization of U t . 4.3. Second order IMEX scheme for the Boltzmann equation. In the following section, which is devoted to numerical simulations to the Boltzmann equation, we also have implemented a second order IMEX scheme :

(4.14)                                  f ⋆ = ε ε + β ⋆ ∆t [f n -∆t v ∇ x f n ] + ∆t Q(f n ) -P (f n ) ε + β ⋆ ∆t + β ⋆ ∆t ε + β ⋆ ∆t M ⋆ , f n+1 = ε ε + β n+1 ∆t/2 [f n -∆t v ∇ x f ⋆ ] + ∆t Q(f ⋆ ) -P (f ⋆ ) ε + β n+1 ∆t/2 + ∆t 2ε + β n+1 ∆t β n+1 M n+1 + β n (M n -f n ) ,
where β ⋆ = β(ρ ⋆ , T ⋆ ) and M ⋆ is the Maxwellian distribution computed from f ⋆ .

4.4. Space discretization. On the one hand, the approximation of the Boltzmann operator is performed by a fast spectral Fourier-Galerkin method already proposed in [START_REF] Filbet | Solving the Boltzmann equation in N log 2 N[END_REF]. On the other hand, the approximation in space is achived using a second order finite volume scheme. Let (x i+1/2 ) i∈I a set of points of the space domain and I a bounded set of integers, hence for ∆x =

x i+1/2 -x i-1/2 x i+1/2 x i-1/2 v x ∂f ∂x dx = F i+1/2 -F i-1/2 ∆x ,
where

F i+1/2 = v + x f l i+1/2 -v - x f r i+1/2 and v + x = max(v x , 0), v - x = max(-v x , 0), f l i+1/2 = f i + δf i+1/2 2 , f r i+1/2 = f i+1 - δf i+1/2 2 .
and δf represents a slope with a slope limiter (see for instance [START_REF] Jin | The Relaxation Schemes for Systems of Conservation Laws in Arbitrary Space Dimensions[END_REF]).

Numerical tests

In this section we perform several numerical simulations for the Boltzmann equation in different asymptotic regimes in order to check the performance (in stability and accuracy) of our methods. We have implemented the first order (2.6) and second order (2.8) scheme for the approximation of the Boltzmann equation. Here, the Boltzmann collision operator is discretized by a deterministic method [START_REF] Pareschi | Numerical solution of the Boltzmann equation I: Spectrally accurate approximation of the collision operator[END_REF][START_REF] Filbet | A numerical method for the accurate solution of the Fokker-Planck-Landau equation in the non homogeneous case[END_REF][START_REF] Filbet | High order numerical methods for the space non-homogeneous Boltzmann equation[END_REF][START_REF] Filbet | Accurate numerical methods for the collisional motion of (heated) granular flows[END_REF][START_REF] Filbet | Solving the Boltzmann equation in N log 2 N[END_REF][START_REF] Duclous | High order resolution of the Maxwell-Fokker-Planck-Landau model intended for ICF applications[END_REF], which gives a spectrally accurate approximation. A classical second order finite volume scheme with slope limiters is applied for the transport operator as sdescribed in section 4.4.

For all numerical simulations, we have considered Maxwellian molecules, that is α = 0 in (1.3). Hence, we take β = 2π ρ such that both operators P (f ) and the full Boltzmann operator Q(f ) have the same loss term corresponding to the dissipative part. 5.1. Approximation of smooth solutions. This test is used to evaluate the order of accuracy of our new methods. More precisely, we want to show that our methods (4.4) and (4.14) are uniformly accurate with respect to the parameter ε > 0. We consider the Boltzmann equation (1.1) in 1 d x ×2 d v . We take a smooth initial data

f 0 (x, v) = ρ 0 (x) 2π T 0 (x) exp - |v| 2 2 T 0 (x) , (x, v) ∈ [-L, L] × R 2 ,
with ρ 0 (x) = (11 -9 tanh(x))/10, T 0 (x) = (3tanh(x))/4, L = 1 and assume specular reflection boundary conditions in x. Numerical solutions are computed from different phase space meshes : the number of point in space is n x = 50, 100, 200,...,1600 and the number of points in velocity is n 2 v with n v = 8,...,64 (for which the spectral accuracy is achieved), the time step is computed such that the CFL condition for the transport is satisfied ∆t ≤ ∆x/v max , where ∆x is the space step and v max = 7 is the truncation of the velocity domain. Then different values of ε are considered starting from the fully kinetic regime ε = 1, up to the fluid limit ε = 10 -5 corresponding to the solution of the Euler system (1.7). The final time is T max = 1 such that the solution is smooth for the different regimes.

An estimation of the relative error in L p norm is given by

e 2 h = max t∈(0,T ) f h (t) -f 2 h (t) p f 0 p , 1 ≤ p ≤ +∞,
where f h represents the approximation computed from a grid of order h. The numerical scheme is said to be k-th order if e 2 h ≤ C h k , for all 0 < h ≪ 1. In Figure 3, the L 1 and L ∞ errors of the second order method (4.14) are presented. They show a uniformly second order convergence rate (an estimation of the slope is 1.9) in space and time (the velocity discretization is spectrally accuracy in v thus does not contribute much to the errors). The time step is not constrained by the value of ε, showing a uniform stability in time.

ε log(h) 1.9 x ε ε ε ε ε -5 -4 -3 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 = 1.E-5 = 1.E-2 = 1.E-1 = 2.E-1 = 5.E-1 = 1.E+0 ε log(h) 1.9 x ε ε ε ε ε -4 -3 -2 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 = 1.E-5 = 1.E-2 = 1.E-1 = 2.E-1 = 5.E-1 = 1.E+0 (1) (2)

5.2.

The Riemann problem. This test deals with the numerical solution to the 1d x ×2d v Boltzmann equation for Maxwellian molecules (γ = 0). We present numerical simulations for one dimensional Riemann problem and compute an approximation for different Knudsen numbers, from rarefied regime to the fluid regime.

Here, the initial data corresponding to the Boltzmann equations are given by the Maxwellian distributions computed from the following macroscopic quantities

   (ρ l , u l , T l ) = (1, 0, 1) , if 0 ≤ x ≤ 0.5 , (ρ r , u r , T r ) = (0.125, 0, 0.25) , if 0.5 < x ≤ 1 .
We perform several computations for ε = 1, 10 -1 , 10 -2 ,...,10 -4 . In Figure 4, we only show the results obtained in the kinetic regime (10 -2 ) using a spectral scheme for the discretization of the collision operator [START_REF] Filbet | Solving the Boltzmann equation in N log 2 N[END_REF] (with n v = 32 2 and a truncation of the velocity domain v max = 7) and second order explicit Runge-Kutta and second order method (4.14) for the time discretization with a time step ∆t = 0.005 satisfying the CFL condition for the transport part (with n x = 100). For such a value of ε, the problem is not stiff and this test is only performed to compare the accuracy of our second order scheme (4.14) with the classical (second order) Runge-Kutta method. We present several snapshots of the density, mean velocity, temperature and heat flux

Q(t, x) := 1 ε R dv (v -u ε ) |v -u ε | 2 f ε (t, x, v)dv
at different time t = 0.10 and 0.20. Both results agree well with only n x = 100 in the space domain and n v = 32 for the velocity space. Thus, in the kinetic regime our second order method (4.14) gives the same accuracy as a second order fully explicit scheme without any additional computational effort. (3) (4)

Figure 4. Riemann problem (ε = 10 -2 ), crosses (x) represent the numerical solution obtained with our second order method (4.14) and lines with the explicit Runge-Kutta method: evolution of (1) the density ρ, (2) mean velocity u, (3) temperature T and (4) heat flux Q at time t = 0.05, 0.1, 0.15 and 0.2. Now, we investigate the cases of small values of ε for which an explicit scheme requires the time step to be of order O(ε). In order to evaluate the accuracy of our method (4.14) in the Navier-Stokes regime (for small ε ≪ 1 but not negligible), we compared the numerical solution for ε = 10 -3 with one obtained with a small time step ∆t = O(ε) (for which the computation is still feasible). Note that a direct comparison with the numerical solution to the compressible Navier-Stokes system (1.8) is difficult since the viscosity µ ε = µ(T ε ) and the thermal conductivity κ ε = κ(T ε ) are not explicitly known. Therefore, in Figure 5, we report the numerical results for ε = 10 -3 and propose a comparison between the numerical solution obtained with the scheme (4.14) and the one obtained with a second order explicit Runge-Kutta method. In this case, the behavior of macroscopic quantities (density, mean velocity, temperature and heat flux) agree very well even if the time step is at least ten times larger with our method (4.4)

Figure 5. Riemann problem (ε = 10 -3 ), crosses (x) represent the numerical solution obtained with our second order method (4.14) and lines with the explicit Runge-Kutta method: evolution of (1) the density ρ, (2) mean velocity u, (3) temperature T and (4) heat flux Q at time t = 0.05, 0.1, 0.15 and 0.2.

Then in Figure 6, we compare the numerical solution of the Boltzmann equation (1.1) with the numerical solution to the compressible Navier-Stokes system derived from the BGK model since the viscosity and heat conductivity are in that case explicitly known [START_REF] Bennoune | Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible NavierStokes asymptotics[END_REF]. To approximate the compressible Navier-Stokes system, we apply a second order Lax-Friedrich scheme using a large number of points (n x = 1000) whereas we only used n x = 100, and 200 points in space and n 2 v = 32 2 points in velocity for the approximation of the kinetic equation (1.1). In this problem, the density, mean velocity and temperature are relatively close to the one obtained with the approximation of the Navier-Stokes system. Even the qualitative behavior of the heat flux agrees well with the heat flux corresponding to the compressible Navier-Stokes system κ ε ∇ x T ε , with κ ε = ρ ε T ε (see Figure 6), yet some differences can be observed, which means that the use of BGK models to derive macroscopic models has a strong influence on the heat flux.

Finally in Figure 7, we present a comparison to the numerical solution obtained with our AP scheme for a very small value of ε = 10 -8 with the numerical solution to the Euler system. The agreement on the density, mean velocity and temperature is very satisfying with only n x = 100 in the space domain for the solution to the kinetic model.

5.3.

A problem with mixing regimes. Now we consider the Boltzmann equation (1.1) with the Knudsen number ε > 0 depending on the space variable in a wide range of mixing scales. This kind of problem was already studied by several authors for the BGK model [START_REF] Degond | A smooth transition model between kinetic and hydrodynamic equations[END_REF] or the radiative transfer equation [START_REF] Jin | Uniformly accurate diffusive relaxation schemes for multiscale transport equations[END_REF]. In this problem, ε : R → R + is given by

ε(x) = ε 0 + 1 2 [tanh(1 -11 x) + tanh(1 + 11 x) ] ,
which varies smoothly from ε 0 to O(1). This numerical test is difficult because different scales are involved. It requires a good accuracy of the numerical scheme for all range of ε. In order to focus on the multi-scale nature we only consider periodic boundary conditions, even if the method has also been used with specular reflection in space. Furthermore, to increase the difficulty we consider an initial data which is far from the local equilibrium of the collision operator:

f 0 (x, v) = ρ 0 2 exp - |v -u 0 | 2 T + exp - |v + u 0 | 2 T 0 , x ∈ [-L, L], v ∈ R 2
with u 0 = (3/4, -3/4),

ρ 0 (x) = 2 + sin(k x) 2 , T 0 (x) = 5 + 2 cos(k x) 20
where k = π/L and L = 1/2.

Here we cannot compare the numerical solution with the one obtained by a macroscopic model. From the numerical simulations, we observe that the solution is smooth during a short time and some discontinuities are formed in the region where the Knudsen number ε is very small and then propagate into the physical domain.

On the one hand, we only take ε 0 = 10 -3 in order to propose a comparison of numerical solutions computed with a second order method using a time step ∆t = 0.001 (such that the CFL condition for the transport part is satisfied) and the one by the second order explicit Runge-Kutta method with a smaller time step ∆t = 0.0001) to get stability. The number of points in space is n x = 200 and in velocity is n 2 v = 32 2 . Clearly, in Figure 8, the results are in good agreement even if our new method does not solve accurately small time scales when the solution is far from the local equilibrium. Moreover in Figure 9, we present numerical results with only n x = 50 and n x = 200, and n 2 v = 32 2 to show the performance of the method with a small number of discretization points in space. With n x = 50 points the qualitative behavior of the macroscopic quantities (ρ, u, T ) is fairly good.

On the other hand, we have performed different numerical results when ε 0 = 10 -4 , then the variations of ε starts from 10 -4 to 1 in the space domain. In that case, the computational time of a fully explicit scheme would be more than one hundred times larger than the one required for the asymptotic preserving scheme (4.14). We observe that discontinuities appear on the density, mean velocity and temperature and then propagate accurately into the domain. The shock speed is roughly the same for the different numerical resolutions. Therefore, this method gives a very good compromise between accuracy and stability for the different regimes. Numerical results are not plotted since they are relatively close to the ones presented in Figures 8 and9.

Other applications: a nonlinear diffusion equation

In this section, we want to illustrate the efficiency of the asymptotic preserving scheme to treat high order differential operators. Such a scheme has already been applied to Willmore flow, a fourth order differential operator [START_REF] Smereka | Semi-implicit level set methods for curvature and surface diffusion motion. Special issue in honor of the sixtieth birthday of Stanley Osher[END_REF]. Here, we consider the flow of gas in a two dimensional porous medium with initial density g 0 (v) ≥ 0. The distribution function g(t, v) then satisfies the nonlinear degenerate

     ∂g ∂t = ∆ v g m , v ∈ R dv , g(t = 0, v) = g 0 (v), v ∈ R dv ,
where m > 1 is a physical constant. Assuming that R 2

(1 + |v| 2 ) g 0 (v)dv < +∞, J.A. Carrillo and G. Toscani [START_REF] Carrillo | Asymptotic L 1 -decay of solutions of the porous medium equation to self-similarity[END_REF] proved that g(t, v) behaves asymptotically in a self-similar way like the Barenblatt-Pattle solution, as t → +∞. More precisely, it is easy to see that if we consider the change of variables (6.2)

g(t, v) = 1 s(t) f log(s(t)), v s(t) ,
where s(t) := √ 1 + 2t, the new distribution function f is solution to

∂f ∂t = ∇ v • (v f + ∇ v f m ) ,
and converges to the Barenblatt-Pattle distribution where C is uniquely determined and depends on the initial mass g 0 but not on the "details" of the initial data. Instead of working on (6.1) directly, we will study the asymptotic decay towards its equilibrium. The key argument on the proof of J.A. Carrillo and G. Toscani is the control of the entropy functional

M(v) = C - m -1 2 m |v| 2 1/(m-1) + , 1
H(f ) = R 2 |v| 2 f (t, v) + m m -1 f m (t, v) dv, which satisfies dH(f ) dt = - R 2 f (t, v) v + m m -1 ∇f m-1 2 dv ≤ 0
or the control of the relative entropy H(f |M) = H(f ) -H(M) with respect to the steady state M. Numerical discretization of this problem leads to the following difficulty : explicit schemes are constrained by a CFL condition ∆t ≃ ∆v 2 whereas implicit schemes require the numerical resolution of a nonlinear problem at each time step (with a local constraint on the time step). We refer to [START_REF] Chainais-Hillairet | Asymptotic behavior of a finite volume scheme for the transient drift-diffusion model[END_REF][START_REF] Filbet | A finite volume scheme for the Patlak-Keller-Segel chemotaxis model[END_REF] for a fully implicit approximation preserving steady states for nonlinear Fokker-Planck type equations.

Here we do not focus on the velocity discretization, but only want to apply our splitting operator technique to remove this severe constraint on the time step. Here the parameter ε does not represent a physical time scale but is only related to the velocity space discretization ∆v. Therefore, we set

Q(f ) = ∇ v • (v f + ∇ v f m ) and P (f ) = ∇Q(M) (f -M), which leads to the following decomposition ∂f ∂t = ∇ v • v M + ∇ v (f m -m M m-1 (f -M)) non dissipative part + ∇ v • v (f -M) + m ∇ v M m-1 (f -M) stiff, dissipative linear part .
Then we apply a simple IMEX scheme which only requires the numerical resolution of a linear system at each time step. We choose m = 3 and a discontinuous initial datum far from the equilibrium f 0 (v) = l∈{1,2} k∈{0,...,n-1}

1 10 1 B(0,r0) (v -v k,l )
where n = 12, r 0 = 1/4 and v k,l = l e i θ k , with θ k = 2 k π/n, k = 0, . . . , n -1. We use a standard velocity discretization in the velocity space based on an upwind finite volume approximation for the transport term and a center difference for the diffusive part. We take n 2 v = 120 2 in velocity and a time step ∆t = 0.02 which is much larger than the time step satisfying a classical CFL condition for this problem ∆t ≃ O(∆v 2 ). The numerical scheme (2.6) is still stable and the numerical solution preserves nonegativity at each time step (see Figure 10)! For large time, the solution converges to an approximation of the steady state even if the present scheme is not exactly well-balanced (it does not preserve exactly the steady state). Moreover, to get a better idea on the behavior of the numerical solution, we plot the evolution of the entropy and its dissipation for different time steps. More surprisingly, the numerical entropy is decreasing and the dissipation converges towards zero when time goes to infinity.

Conclusion

We have proposed a new class of numerical schemes for physical problems with multiple time and spatial scales described by a stiff nonlinear source term. A prototype equation of this type is the Boltzmann equation for rarefied gas. When the Knudsen number is small, the stiff collision term of the Boltzmann equation drives the density distribution to the local Maxwellian, thus the macroscopic quantities such as mass, velocity and temperature evolve according to fluid dynamic equations such as the Euler or Navier-Stokes equations. Asymptotic-preserving (AP) schemes for kinetic equations have been successful since they capture the fluid dynamic behavior even without numerically resolving the small Knudsen number. However, the AP schemes need to treat the stiff collision terms implicitly, thus it yields a complicated numerical algebraic problem due to the nonlinearity and non-locality of the collision term. In this paper, we propose to augment the nonlinear Boltzmann collision operator by a much simpler BGK collision operator, and apply an implicit scheme only on the BGK operator which can be handled much more easily. For hyperbolic systems with relaxations We show that this method is AP in the Euler regime, after the initial transient time, and is also consistent to the Navier-Stokes approximations for suitably small time steps and mesh sizes. Numerical examples, including those with mixing scales and non-local-Maxwellian initial data, demonstrate the AP property as well as uniform convergence (in the Knudsen number) of this method.

This method can be extended to a wide class of PDEs (or ODEs) with stiff source terms that admit a stable and unique local equilibrium. One example is the hyperbolic system with relaxations which are studied in this paper. We also use the nonlinear Fokker-Planck equation as an example to illustrate this point, and will pursue more applications in the future.

It is worth to mention that the present method is essentially based on a decomposition of the nonlinear operator as the sum of a linear and dissipative part and a nonlinear part (2.5). Therefore, it does not need a specific velocity and space discretization and can be easily applied to different stochastic and deterministic schemes. Moreover, based on this decomposition, other schemes can be constructed [17]. In this paper we do not mention the numerical treatment of boundary conditions although the method naturally applies to periodic and specular reflection boundary conditions. However, for physical boundary conditions, as Maxwell diffusive conditions, boundary layers will be generated where the solution is very far away from Gaussian distributions [START_REF] Sone | Kinetic Theory and Fluid Dynamics[END_REF]. Therefore, adequate space and time discretizations deserve attention and will constitute a very interesting problem that we would like to deal with. We refer to [START_REF] Gamba | Shock and Boundary Structure formation by Spectral-Lagrangian methods for the Inhomogeneous Boltzmann Transport Equation[END_REF] for a numerical treatment of boundary conditions for the Boltzmann equation using deterministic method [START_REF] Gamba | Shock and Boundary Structure formation by Spectral-Lagrangian methods for the Inhomogeneous Boltzmann Transport Equation[END_REF], where the influence of boundary conditions is studied far away the boundary. This work was not aimed at the AP property. There have been very few studies on AP schemes for boundary value problems, except those on linear transport equation in the diffusive regimes with Dirichlet boundary conditions [START_REF] Jin | The Discrete-Ordinate Method in Diffusive Regime[END_REF][START_REF] Jin | Fully-discrete Numerical Transfer in Diffusive Regimes[END_REF][START_REF] Golse | The Convergence of Numerical Transfer Schemes in Diffusive Regimes I: The Discrete-Ordinate Method[END_REF]. This is an important subject for future research.

Figure 2 .

 2 Figure 2. Approximation of the solution to (3.8) obtained from our AP scheme (3.1) and the linear penalty method (3.4) for different values of the Knudsen number ε = 10 -1 , 10 -2 and 10 -7 .

Figure 3 .

 3 Figure 3. The L 1 and L ∞ errors of the second order method (4.14) for different values of the Knudsen number ε = 10 -5 , . . . , 1.

Figure 6 .

 6 Figure 6. Riemann problem (ε = 10 -4 ), comparison between the numerical solution to the Boltzmann equation with our second order method (4.14) represented with crosses (x) and the numerical solution to the compressible Navier-Stokes system (lines): evolution of (1) the density ρ, (2) mean velocity u, (3) temperature T and (4) heat flux Q at time t = 0.05, 0.1, 0.15 and 0.2.

Figure 7 .

 7 Figure 7. Riemann problem (ε = 10 -8 ), comparison with the solution to the Euler system: evolution of (1) the density ρ, (2) mean velocity u, and (3) temperature T at time t = 0.05, 0.1, 0.15 and 0.2.

Figure 8 .

 8 Figure 8. Mixing regime problem (ε 0 = 10 -3 ), comparison of the numerical solution to the Boltzmann equation with the second order method (4.14) represented with crosses (x) with the numerical solution obtained with the explicit Runge-Kutta method using a small time step (line): evolution of (1) the density ρ, (2) mean velocity u, (3) temperature T at time t = 0.25, 0.5 and 0.75.

Figure 9 .

 9 Figure 9. Mixing regime problem (ε 0 = 10 -3 ), comparison of the numerical solution to the Boltzmann equation obtained with the AP scheme (4.14) using n x = 50 (crosses x) and n x = 200 points (line): evolution of (1) the density ρ, (2) mean velocity u, (3) temperature T at time t = 0.25, 0.5 and 0.75.

Figure 10 .

 10 Figure 10. Nonlinear Fokker-Planck solution: convergence toward equilibrium (Barenblatt-Pattle distribution) obtained with the first order method (2.6) using n x = 100 at time t = 0.1, 0.4, 0.8, 1.0, 1.2 and 4 with a large time step.

Figure 11 .

 11 Figure 11. Nonlinear Fokker-Planck solution: convergence toward equilibrium (Barenblatt-Pattle distribution) obtained with the first order method (2.6) using n x = 100 with ∆t = 0.02 and 0.001.

  or (4.14).

	1 1 1 1 1 1 1 1				t=0.05				mean velocity u
	0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9				t=0.10		0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8
	0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8
	0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7		t=0.15				0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6		t=0.05
	0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6	t=0.20						t=0.10
	0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5						0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4		t=0.15
	0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4								t=0.20
	0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3						0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2
	0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2	density	ρ				0 0 0 0 0 0 0 0
	0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
	0 0 0 0 0 0 0 0		0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2	0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4	0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6	0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8	1 1 1 1 1 1 1 1	0 0 0 0 0 0 0 0	0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2	0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4	0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6	0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8	1 1 1 1 1 1 1 1

  1, 0.15 and 0.2.

	1 1 1 1 1 1 1 1				t=0.05				mean velocity u				1 1 1 1 1 1 1 1
	0.6 0.7 0.8 0.9 0.6 0.7 0.8 0.9 0.6 0.7 0.8 0.9 0.6 0.7 0.8 0.9 0.6 0.7 0.8 0.9 0.6 0.7 0.8 0.9 0.6 0.7 0.8 0.9 0.6 0.7 0.8 0.9	t=0.15 t=0.20		t=0.10		0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8		t=0.10 t=0.05			0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9			t=0.10	t=0.15	t=0.20
	0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5						0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4		t=0.15				0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6
	0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4								t=0.20				0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
	0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3						0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2						0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4			t=0.05
	0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2	density	ρ				0 0 0 0 0 0 0 0						0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3		temperature T
	0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1												0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2
	0 0 0 0 0 0 0 0		0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2	0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4	0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6	0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8	1 1 1 1 1 1 1 1	0 0 0 0 0 0 0 0	0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2	0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4	0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6	0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8	1 1 1 1 1 1 1 1	0 0 0 0 0 0 0 0	0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2	0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4	0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6	0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8	1 1 1 1 1 1 1 1
				(1)						(2)						(3)

Acknowledgments. F. Filbet thanks Ph. Laurençot, and L. Pareschi for interesting discussions on the topic.

F. Filbet is partially supported by the French ANR project "Jeunes Chercheurs" Méthodes Numériques pour les Équations Cinétiques (MNEC) and by the European Research Council ERC Starting Grant 2009, project 239983-NuSiKiMo. S. Jin was partially supported by NSF grant No. DMS-0608720, NSF FRG grant DMS-0757285, and a Van Vleck Distinguished Research Prize from University of Wisconsin-Madison.

1