From Boltzmann Equation for Granular Gases to a Modified Navier–Stokes–Fourier System

In this paper, we give an overview of the results established in Alonso (http://arxiv.org/org/abs/2008.05173, 2020) which provides the first rigorous derivation of hydrodynamic equations from the Boltzmann equation for inelastic hard spheres in 3D. In particular, we obtain a new system of hydrodynamic equations describing granular flows and prove existence of classical solutions to the aforementioned system. One of the main issue is to identify the correct relation between the restitution coefficient (which quantifies the rate of energy loss at the microscopic level) and the Knudsen number which allows us to obtain non trivial hydrodynamic behavior. In such a regime, we construct strong solutions to the inelastic Boltzmann equation, near thermal equilibrium whose role is played by the so-called homogeneous cooling state. We prove then the uniform exponential stability with respect to the Knudsen number of such solutions, using a spectral analysis of the linearized problem combined with technical a priori nonlinear estimates. Finally, we prove that such solutions converge, in a specific weak sense, towards some hydrodynamic limit that depends on time and space variables only through macroscopic quantities that satisfy a suitable modification of the incompressible Navier–Stokes–Fourier system.


INTRODUCTION
In this paper, we report on some recent results obtained in [3] about the problem of deriving rigorously some hydrodynamic limit from the Boltzmann equation for inelastic hard spheres with small inelasticity. Our aim here is to give an account of the main aspects of our work [3] in a shorterreader-friendly -version that includes the main results as well as the main ideas and arguments. We shall only sketch the proofs of our results, referring the reader to [3] for complete versions and details.

The problem.
The kinetic model. We consider here the (freely cooling) Boltzmann equation which provides a statistical description of identical smooth hard spheres suffering binary and inelastic collisions: is the density of granular gases having position x ∈ T d ℓ and velocity v ∈ R d at time t 0. We consider here for simplicity the case of flat torus for some typical length-scale ℓ > 0. The so-called restitution coefficient α belongs to (0, 1] and the collision operator Q α is defined in weak form aŝ (1.4) and the post-collisional velocities (v ′ , v ′ * ) are given by where q = v − v * ,q = q/|q|. (1.5) Here, dσ denotes the Lebesgue measure on S d−1 and the angular part b = b(σ ·q) of the collision kernel appearing in (1.4) is a non-measurable mapping integrable over S d−1 . There is no loss of generality assuminĝ Notice that one can also give a strong formulation of the collision operator Q α (see [3,Appendix A]). This strong formulation is simpler in the elastic case (α = 1), we here give it for later use: The true definition actually involves pre-collisional velocities and not post-collisional velocities v ′ and v ′ * but they match in the elastic case, which explains the formula (1.6).
The fundamental distinction between the classical elastic Boltzmann equation and the associated to granular gases lies in the role of the parameter α ∈ (0, 1), the coefficient of restitution that we suppose constant. This coefficient is given by the ratio between the magnitude of the normal component (along the line of separation between the centers of the two spheres at contact) of the relative velocity after and before the collision. The case α = 1 corresponds to perfectly elastic collisions where kinetic energy is conserved. However, when α < 1, part of the kinetic energy of the relative motion is lost since Notice that the microscopic description (1.5) preserves the momentum v ′ + v ′ * = v + v * and, taking ψ = 1 and then ψ = v in (1.3) yields the following conservation of macroscopic density and bulk velocity defined as Consequently, there is no loss of generality in assuming that R(t) = R(0) = 1, U (t) = U (0) = 0, ∀ t 0.
The main contrast between elastic and inelastic gases is that in the latter the granular temperature, where D α ( · , · ) denotes the normalised energy dissipation associated to Q α , see [16], given by where γ b is a positive constant depending only on the angular kernel b.

The problem of hydrodynamic limits.
To capture some hydrodynamic behaviour of the gas, we need to write the above equation in nondimensional form introducing the dimensionless Knudsen number which is proportional to the mean free path between collisions. We then introduce the classical Navier-Stokes rescaling of time and space (see [5]) to capture the hydrodynamic limit and introduce the particle density (1. 8) In this case, we choose for simplicity ℓ = ε in (1.2) which ensures now that F ε is defined on Under such a scaling, F ε satisfies the rescaled Boltzmann equation supplemented with the initial condition Conservation of mass and density is preserved under this scaling, if F ε solves (1.9a), then whereas the cooling of the granular gas is given by the equation x, v) dv dx and we recall that D α is defined in (1.7). The conservation properties of the equation imply that there is no loss of generality assuming that In order to understand the free-cooling inelastic Boltzmann equation (1.9a)-(1.9b), we perform a selfsimilar change of variables, which allows us to introduce an intermediate asymptotic and ensures that our equation has a non trivial steady state (see [15,16,17] for more details). After this change of variables, we are led to study the equation . Note that the drift term acts as an energy supply which prevents the total cooling down of the gas. It has been shown that there exists a spatially homogeneous steady state G α to (1.11). More specifically, there exists α 0 ∈ (0, 1) (where α 0 is an explicit threshold value) such that for α ∈ (α 0 , 1), there exists a unique distribution G α = G α (v) satisfying Moreover, there exists some constant C > 0 independent of α such that for some explicit temperature ϑ 1 > 0. The Maxwellian distribution M is a steady solution for α = 1 and its prescribed temperature ϑ 1 (which ensures (1.13) to hold) will play a role in the rest of the analysis. It is important to emphasize that, in all the sequel, all the threshold values on ε and the various constants involved are actually depending only on this initial choice.
In order to reach some incompressible Navier-Stokes type equation in the limit ε → 0, we introduce the following fluctuation h ε around the equilibrium G α : Our problem boils down to look at the following equation on h ε : where L α is the linearized collision operator (local in the x-variable) defined as We also denote by L 1 the linearized operator around G 1 = M, that is, From now on, we will always assume that The choice of prescribing as initial energy some constant E ε > 0 satisfying ε −1 (E ε − dϑ 1 ) → 0 as ε → 0 for our problem is natural because dϑ 1 is the energy of the Maxwellian M introduced in (1.14) and as we shall see later on, the restitution coefficient α is intended to tend to 1 as ε goes to 0 in our analysis (see (1.22)). It is also worth noticing that assumption (1.18) and (1.12) result in Moreover, equation (1.15) preserves mass and vanishing momentum since, if h ε solves (1.15), then one formally has Consequently, there is no loss of generality assuming that (1.20)

Relation between the restitution coefficient and the Knudsen number.
The central underlying assumption in our study is the following relation between the restitution coefficient and the Knudsen number.
Notice that under this assumption, the hypothesis made on the energy of the initial data in (1.18) implies thatˆT Indeed, using (1.18) and Assumption 1.1 combined with (1.13), we obtain Still under Assumption 1.1, we formally obtain that if ε → 0 in (1.15), then h ε → h with h ∈ Ker L 1 where L 1 is defined in (1.17). We recall that when seeing L 1 as an operator acting only on velocity on the and the projection π 0 onto Ker L 1 is given by We deduce formally that h takes the following form It is worth mentioning that a careful spectral analysis of the linearized collision operator L α defined in (1.16) shows that unless one assumes 1 − α at least of order ε 2 , the eigenfunction associated to the energy dissipation would explode and prevent some exponential stability for (1.15) to hold true (see Theorem 2.1). Actually, in our study, we will require λ 0 to be relatively small with respect to the spectral gap associated to the elastic linearized operator to ensure stability in the inelastic case. If one assumes λ 0 = 0 (for example, one could assume 1 − α of order ε q with q > 2), the effect of the inelasticity is too weak in the hydrodynamic scale and the expected model is the classical Navier-Stokes-Fourier system. In short, we are left with two cases: Case 1: If λ 0 = 0, the expected model is the classical Navier-Stokes-Fourier system. Case 2: If 0 < λ 0 < ∞ is small enough (compared to some explicit quantities), the cumulative effect of inelasticity is visible in the hydrodynamic scale and we expect a different model to the Navier-Stokes-Fourier system accounting for that. In this nearly elastic regime, the energy dissipation rate in the system happens in a controlled fashion since the inelasticity parameter is compensated accordingly to the number of collisions per time unit. Other regimes can be considered depending on the rate at which kinetic energy is dissipated; for example, an interesting regime is the mono-kinetic one which considers the extreme case of infinite energy dissipation rate. In this way, the limit is formally described by enforcing a Dirac mass solution in the kinetic equation yielding the pressureless Euler system (corresponding to sticky particles). Such a regime has been rigorously addressed in the one-dimensional framework in the interesting contribution [11]. It is an open question to extend such analysis to higher dimensions since the approach of [11] uses the so-called Bony functional which is a tool specifically tailored for 1D kinetic equations.
1.2. Notations and definitions. Let us introduce some useful notations for functional spaces. For any nonnegative weight function m : R d → R + , we define, for all p > 1 the space L p (m) through the norm We also define, for p 1 with the usual norm, i.e., for k ∈ N: For m ≡ 1, we simply denote the associated spaces by L p and W k,p . Notice that all the weights we consider here will depend only on velocity, i.e. m = m(v). We will also use the notation ξ := 1 + |ξ| 2 for ξ ∈ R d . On the complex plane, for any a ∈ R, we set and, for any r > 0, we set D(r) = {z ∈ C ; |z| r}.
We also introduce the following notion of hypo-dissipativity in a general Banach space (X, · ). A closed (unbounded) linear operator A : D(A) ⊂ X → X is said to be hypo-dissipative on X if there exists a norm, denoted by ||| · |||, equivalent to the · -norm such that A is dissipative on the space (X, ||| · |||), that is, Given two Banach spaces X and Y , we denote with · X→Y the operator norm on the space of B(X, Y ) linear and continuous operators from X to Y . Note also that in what follows, for two positive quantities A and B, we denote by A B if there exists a universal positive constant C (which is in particular independent of the parameters α and ε) such that A CB.  (1.15). The first one is the following Cauchy theorem regarding the existence and uniqueness of close-to-equilibrium solutions to (1.15). The functional spaces at stake are

Remark 1.3.
It is worth pointing out that the close-to-equilibrium solutions we construct are shown to decay with an exponential rate as close as we want to (which is the energy eigenvalue of the linearized operator, see Theorem 2.1 hereafter). The rate of convergence can thus be made uniform with respect to the Knudsen number ε (notice that if λ 0 = 0, we obtain a rate of decay as close as we want to η(ε), we thus obtain a uniform bound in time but not a uniform rate of decay).
The estimates on the solution h ε provided by Theorem 1.2 are enough to prove that the solution h ε (t) converges towards some hydrodynamic solution h which depends on (t, x) only through macroscopic quantities (̺(t, x), u(t, x), θ(t, x)) which are solutions to a suitable modification of the incompressible Navier-Stokes system. This is done under an additional assumption on the initial datum that is lightly restrictive. Before stating our main convergence result, we introduce the notation and we furthermore assume that in the definition of the functional spaces (1.27), the following conditions are satisfied: Theorem 1.4. We suppose that the assumptions of Theorem 1.2 are satisfied. We assume furthermore that there exists Then, for any T > 0, the family of solutions {h ε } ε constructed in Theorem 1.2 converges in some weak sense to where P is the Leray projection on divergence-free vector fields and (̺ 0 , u 0 , θ 0 ) have been introduced in (1.28).
The viscosity ν > 0 and heat conductivity γ > 0 are explicit and λ 0 > 0 is the parameter appearing in Assumption 1.1. The parameterc > 0 is depending on the collision kernel b( · ).
Remark 1.5. The data that we consider here are actually quite general. Indeed, the assumption that we make only tells that the macroscopic projection of h in ε converges towards some macroscopic distribution and we do not make any assumption on the macroscopic quantities of this distribution. Namely, we do not suppose that the divergence free and the Boussinesq relations are satisfied by (̺ 0 , u 0 , θ 0 ), the initial layer that could be created by such a lack of assumption is actually absorbed in our notion of weak convergence, the precise notion of which being very peculiar and strongly related to the a priori estimates used for the proof of Theorem 1.2 (see Theorem 4.2 for more details on the type of convergence).
To prove Theorem 1.4, our approach is reminiscent of the program established in [5,6,9,19] but simpler because our solutions are stronger than the renormalized ones that are used in [9]. It is based on computations and compactness arguments that were already used in the elastic case. Let us point out that in our case, additional terms appear due to the inelasticity and they can be handled in the framework of Assumption 1.1. In Section 4, we present the proof but only mention its main steps and arguments (details can be found in [3, Section 6]).

STUDY OF THE KINETIC LINEARIZED PROBLEM
2.1. Main result on the linearized operator. The first step in the proof of Theorem 1.2 is the spectral analysis of the linearized problem associated to (1.15). To that end, we introduce We are going to state our main result on G α,ε in the space E defined in (1.27) but our analysis actually allows to treat even larger spaces (namely, we can obtain the same result under the softer constraints m k 0 and q > 2) but we only state the linear result in this case because it is the only one that will be used in the rest of the paper. Let us also recall that, in any reasonable space (in particular in E and Y j for j = −1, 0, 1 defined in (2.6)-(2.8)), the elastic operator has a spectral gap: there exists µ ⋆ > 0 such that where 0 is an eigenvalue of algebraic multiplicity d + 2 of G 1,ε associated to the eigenfunctions (recall that C µ⋆ is defined in (1.26)). This can be proven by an enlargement argument due to [10] based on the fact that in the Hilbert space a result of hypocoercivity has been proven in [7] (the constraint m 1 would actually be enough but we will only make use of this result for m > d in the sequel). More precisely, introducing the other Hilbert space there exists µ ⋆ > 0 and a norm equivalent to the usual one uniformly in ε (we still denote it by · H and · , · H its associated scalar product to lighten the notations) such that As we shall see in the following result, the scaling (1.21) in Assumption 1.1 is precisely the one which allows to preserve exactly d + 2 eigenvalues in the neighborhood of zero for G α,ε . Let us now state our main spectral result (see Figure 1 for an illustration where we have denoted λ ε := −λ d+2 (ε)): Theorem 2.1. Assume that Assumption 1.1 is met. For µ close enough to µ ⋆ defined in (2.1) (in an explicit way), there are some explicit ε > 0 and λ > 0 depending only on χ := µ ⋆ − µ such that, for all ε ∈ (0, ε) and λ 0 ∈ [0, λ), the linearized operator G α(ε),ε has the following spectral property in E:

Remark 2.2.
It is worth noticing that the eigenvalue λ 1 (ε) = 0 corresponds with the property of mass conservation of the operator G α(ε),ε . Concerning the intermediate eigenvalues λ j (ε) for j = 2, . . . , d+ 1, as one can see on their definition, they may be positive, this is due to the fact that the collision operator Q α preserves momentum while the drift operator ∇ v (v · ) does not. However, using (1.19), one can prove that the vanishing momentum is preserved by the whole operator G α(ε),ε , consequently, those eigenvalues won't affect the long-time analysis of our problem. Finally, the eigenvalue λ d+2 (ε) is directly linked to the non-preservation of energy property of G α(ε),ε .
We are going to prove Theorem 2.1 in two stages. First, we perform a perturbative argument (reminiscent of [21]) in a L 2 v,x -based Sobolev space, namely in with σ 0 and σ 1 defined in (2.9) and κ > d/2. The key point of our approach is to see G α,ε as a perturbation of the elastic linearized operator G 1,ε . We then use an enlargement argument (from [10]) to extend the result from Y to the space E defined in (1.27). Several remarks are in order: (i) First, let us remind that the global equilibrium of our equation G α defined in (1.12) has some exponential fat tail and in particular, decays more slowly than a standard Maxwellian distribution (see [17]). As a consequence, we can not rely on classical works on the elastic linearized operator which are developed in spaces of type L 2 v,x (M −1/2 ) with M defined in (1.14). To overcome this difficulty, we exploit results coming from [10] in which an enlargement theory has been carried out. The results proven in [10] include a spectral analysis of the elastic Boltzmann operator G 1,1 in larger spaces (in particular of type L 2 v,x ) with "soft weights" that can be polynomial or stretched exponential. In the same line of ideas, these results have been extended to the rescaled elastic operator G 1,ε in [8].
(ii) Let us also point out that the perturbation at stake does not fall into the realm of the classical perturbation theory of unbounded operators as described in [12] because the perturbation is not relatively bounded. Indeed, the domain of G 1,ε in Y is given by W s+1,2 v W ℓ+1,2 x ( v r+1 ) while if one wants to be sharp in terms of rate, the best estimate in terms of functional spaces that we are able to get is where the spaces Y j are defined through with κ > d/2. These estimates (whose proofs can be found in [3,Lemma 3.3]) are a generalization and optimisation of estimates obtained in [17] and are sharp in terms of rate, this sharpness being needed in our analysis since it allows us to deal with the case λ 0 > 0 in Assumption 1.1. (iii) As a consequence, we have to use refined perturbation arguments whose key insights come from [21].
Note however that we drastically simplify the analysis performed in [21] by remarking that the difference operator G α,ε − G 1,ε does not involve any spatial derivative and that we "only" need to develop a spectral analysis of G α,ε without being able to obtain decay properties on the associated semigroup. As a consequence, we don't need to use a spectral mapping theorem, nor do we need to use an iterated version of Duhamel formula and this is crucial in order to reach the optimal scaling (1.21) for our restitution coefficient. (iv) Let us finally mention that we perform our perturbative argument in Y which is a L 2 v,x -based Sobolev space instead of performing it in E (which is L 1 v L 2 x -based) directly. This intermediate step seems necessary because even if L α − L 1 satisfies nice estimates in L 1 v , the use of Fubini theorem is actually crucial to get the rate (1 − α)/ε 2 in estimates of type (2.7).

2.2.
Elements of proof of Theorem 2.1. As mentioned above, the basis of the proof of this theorem is to see L α as a perturbation of L 1 .
We start by giving a splitting of it into two parts: one which has some good regularizing properties (in the velocity variable) and another one which is dissipative. For any δ > 0, one can write L 1 = A (δ) + B where we have used the shorthand notations g = g(v), g * = g(v * ), g ′ = g(v ′ ), g ′ * = g(v ′ * ). We define is an appropriate truncature function. The dissipativity property of B 1 comes from the fact that the truncature function Θ δ is defined so that the first term is small as δ goes to 0 and the fact that there exist σ 0 > 0 and σ 1 > 0 such that As a consequence, B 1 is going to be dissipative for δ small enough. This leads to the following decomposition of L α : and then the following decomposition of G α,ε : Our analysis of this splitting and then of the spectrum of G α,ε relies on several elements: the nice properties of the above-mentioned splitting of L 1 = A (δ) + B (δ) 1 coming from [10], some refined bilinear estimates on the collision operator coming from [1], new estimates on G α −M that are reminiscent of estimates proven in [4] (see [3,Lemma 2.3]) and also new estimates on Q α − Q 1 (see [3, Lemmas 2.1 and 2.2]). Concerning the latter point, we exploit ideas developed in [17] but our situation is more involved because we work in polynomially weighted spaces whereas in [17], the authors were working with stretched exponential weights.
Let us now explain how we develop our perturbative argument to prove Theorem 2.1. The following proposition (which is an adaptation of [21, Lemma 2.16]) is the first step in the development of the perturbative argument and its proof relies on Corollary 2.4 and Lemma 2.5.
To estimate now R(λ, G ε ) B(Y) , one simply notices that from which, as soon as λ ∈ C µ \ D(µ ⋆ − µ), One checks, using the previous computations, that for λ ∈ C µ \ D(µ ⋆ − µ), and deduces (2.15). This achieves the proof. A first obvious consequence of Proposition 2.6 is that, for any µ ∈ (0, µ ⋆ ), there is ε 3 > 0 depending only We denote by P ε (resp. P 0 ) the spectral projection associated to the set One can deduce then the following lemma whose proof is similar to [21,Lemma 2.17].
With Lemma 2.7, we can now end the proof of Theorem 2.1. Sketch of the proof of Theorem 2.1. The structure of S(G ε ) ∩ C µ in the space Y comes directly from Lemma 2.7 together with Proposition 2.6. To describe more precisely the spectrum, one first remarks that This comes from the fact that for each eigenvalue of L α(ε) , the eigenfunction depends only on v and thus remains an eigenfunction for the operator G ε . Since, for ε small enough, the same perturbative argument that we developed above implies that the spectral projection Π L α(ε) associated to S(L α ) ∩ C µ satisfies dim(Range(Π L α(ε) )) = dim(Range(Π L1 )) = d + 2 = dim(Range(P ε )), that is, the eigenvalues λ j (ε) are actually eigenvalues of L α(ε) . The development of the energy eigenvalue λ d+2 (ε) comes from [17]. The conservation of mass gives us that 0 is an eigenvalue for our problem. The intermediate eigenvalues λ j (ε) for j = 2, . . . , d + 1 are obtained thanks to the fact that Notice that all this allows us to find eigenfunctions (that depend only on v) in L 2 v,x ( v r ). Using once more the splitting L α = A (δ) + B δ α defined in (2.10) and the regularizing properties of A (δ) , one can actually prove that our eigenfunctions lie in Y, which yields the conclusion of Theorem 2.1 in the space Y. To extend the result to the space E, we use an enlargement argument coming from [10], we omit the details here and just mention that this argument is based on the splitting G ε = A ε + B ε introduced in (2.11).

STUDY OF THE KINETIC NONLINEAR PROBLEM
Let us recall that the spaces E and E 1 are defined in (1.27). In this section, we assume that Assumption 1.1 is met and consider ε ∈ (0, ε), λ 0 ∈ 0, λ where ε and λ are defined in Theorem 2.1. As in Section 2, to lighten the notations, we write G ε = G α(ε),ε as well as B ε = B α(ε),ε .
3.1. Splitting of the nonlinear inelastic Boltzmann equation. Now that the spectral analysis of the linearized operator G ε in the space E has been performed, in order to prove Theorem 1.2, we are going to prove several a priori estimates for the solutions to (1.15). The crucial point in the analysis lies in the splitting of (1.15) into a system of two equations mimicking a spectral enlargement method from a PDE perspective (see [18,Section 2.3] and [8] for pioneering ideas on such a method). More precisely, using (2.11), the splitting amounts to look for a solution of (1.15) of the form and In order to lighten the notations, in this section, we will write h in , h, h 0 and h 1 instead of h in ε , h ε , h 0 ε and h 1 ε . The goal is to obtain nice nested a priori estimates on h 0 and h 1 . Notice first that our splitting is more complicated than the one of [8] because it relies on perturbative considerations around the elastic case that come out in the equation satisfied by h 0 . As a consequence, our a priori estimates are more intricate and require the use of non standard Gronwall lemma. Notice also that since the initial datum of h 1 is vanishing, we can study the equation on h 1 in any functional space. In particular, we can study it in the Hilbert space H = W m,2 x,v M −1/2 in which we have a good understanding of the elastic linearized operator G 1,ε . Indeed, in this type of spaces, the symmetries of the collision operator Q 1 allow to get some nice hypocoercive estimates (see (2.4)).
Remark 3.1. In [10], the authors treat the elastic case (α = 1) of the non-rescaled equation (ε = 1) and they do not resort to such a splitting method to study the nonlinear equation, their approach is based on the use of a norm which is equivalent to the usual one and is such that G 1,1 is dissipative in this norm in large spaces. Such an approach is no longer usable when one wants to deal with rescaled equations and obtain uniform in ε estimates. Indeed, the definition of the equivalent norm in [10] does not take into account the different behaviors of microscopic and macroscopic parts of the solution with respect to ε: typically, the microscopic part of the solution vanishes as ε → 0 whereas the macroscopic one does not. Conversely, in the splitting method, the equation that defines h 1 is treated thanks to hypocoercivity tricks that allow to distinguish microscopic and macroscopic behaviors.

3.2.
Estimating h 0 . Concerning h 0 , let us first mention that the dissipativity properties of B ε stated in Lemma 2.3 can actually be improved a bit. More precisely, one can show that there exist norms on the spaces E and E 1 that are equivalent to the standard ones (with multiplicative constants independent of ε) that we still denote · E and · E1 and that satisfy: where we have denoted by (S Bε (t)) t 0 the semigroup generated by B ε and ν 0 is defined in Lemma 2.3. Let us also introduce the Banach space E 2 which satisfies the following continuous embeddings: H ֒→ E 2 ֒→ E 1 (recall that E 1 is defined in (1.27)). Let us point out that the spaces E 1 and E 2 allow us to get the following estimates (see [3,Remark 3.5] and [1,2]): where the multiplicative constants are uniform in α. One can then obtain the following proposition: For ν ∈ (0, ν 0 ) (where ν 0 is defined in Lemma 2.3), there exists an explicit ε 5 ∈ (0, ε) (where ε is defined in Theorem 2.1) such that: Sketch of the proof. Using (3.3) as well as (3.4) and recalling that h 0 solves (3.1), we can compute the evolution of h 0 (t) E and estimate it: Using that the embedding E 2 ֒→ H is continuous, recalling that h 0 (0) = h in and choosing ε 5 small enough so that C ε 5 ∆ 0 ν 0 − ν, we obtain We conclude to (3.5) by assuming furthermore that ε 5 ∆ 0 1.

3.3.
Estimating h 1 . We now comment and study the equation satisfied by h 1 . Let us point out that getting estimates on h 1 is trickier than in [8], indeed, in the latter paper, the idea is to estimate separately P 0 h 1 and (Id − P 0 )h 1 where P 0 is the projector onto Ker(G 1,ε ) defined by where the functions Ψ i have been defined in (1.24), and thanks to the properties of preservation of mass, momentum and energy of the whole equation, one could write that P 0 h = 0 so that P 0 h 1 = −P 0 h 0 and directly get an estimate on P 0 h 1 from the one on h 0 . In our case, the energy is no longer preserved which induces additional difficulties. However, we keep the same strategy and start by estimating P 0 h 1 (see Remark 3.4 for a comment on this choice of strategy).
For the sequel, we also introduce Notice that thanks to Cauchy-Schwarz inequality in velocity, one can easily prove that we have P 0 ∈ B(E, H).
One can then obtain the following proposition: For ε ∈ (0, ε 5 ) (ε 5 is defined in Proposition 3.2), Sketch of the proof. Due to the properties of preservation of mass and vanishing momentum of our equation, we have P 0 h = 0 which implies that P 0 h 1 = −P 0 h 0 . Consequently, we easily get an estimate on P 0 h 1 using that P 0 ∈ B(E, H): (3.10) It now remains to estimate Π 0 h 1 . To this end, we first notice that where we used that P 0 h = Π 0 h due to the preservation of mass and vanishing momentum so, using (3.5) and (3.10), we only need to estimate Π 0 h to get an estimate on P 0 h 1 . To this end, we start by computing the evolution of Π 0 h: By direct inspection, using the definition of Π 0 given in (3.8) and the dissipation of energy (1.7) (see [3,Lemmas 4.2 and 4.5]), we obtain: as ε → 0, Similarly, we have by direct computation that so that, using Minkowski's inequality to estimate D α(ε) (h, h), we obtain Gathering previous estimates, we are able to deduce that With this, inequality (3.9) holds by using ε 5 ∆ 0 1 from the proof of Proposition 3.2.
nothing guarantees that Π ε ε −1 Q α(ε) (h, h) remains of order 1 with respect to ε whereas we have seen in (3.11) that due to the dissipation of kinetic energy, Π 0 ε −1 Q α(ε) (h, h) is actually of order ε. This explains our choice of strategy.
Let us now focus on the estimate of (Id − P 0 )h 1 . We can proceed similarly as in [8], using in particular that P 0 Q 1 = 0. Another crucial point is that the source term A ε h 0 can be bounded in H using the fact that A ε ∈ B(E, H) (see Lemma 2.3). Moreover, it is important to mention that the fact that the bound on A ε induces a rate of ε −2 will be counterbalanced by the fact that the semigroup associated with B ε has an exponential decay rate of type e −νt/ε 2 (see (3.3)). We recall that the Hilbert space H 1 is defined in (2.3) and is such that For µ ∈ (0, µ ⋆ ) (where µ ⋆ is defined in (2.1)) and for ∆ 0 small enough, we have that: Sketch of the proof. From (3.2), the fact that P 0 Q 1 (g, g) = 0 and the fact that P 0 commutes with G 1,ε , we can compute the evolution of Φ(t) := (Id − P 0 )h 1 : We now use the hypocoercive norm on H for G 1,ε introduced in (2.4) and also denote by Φ ⊥ the microscopic part of Φ, namely Φ ⊥ := (Id − π 0 )Φ where we recall that π 0 is the projection onto the kernel of L 1 that has been introduced in (1.23). We compute the evolution of Φ(t) 2 H : Notice that we have been able to replace Φ by Φ ⊥ in the second term due to the conservation laws satisfied by Q 1 and the fact that π 0 is orthogonal in H. Then, from the properties of the hypocoercive norm (see (2.4)), using (3.12) and the facts that P 0 ∈ B(H), A ∈ B(E, H) (from Lemma 2.3) as well as Cauchy-Schwarz inequality, we obtain that Making an appropriate use of Young inequality to treat the third term, we obtain that for µ ∈ (0, µ ⋆ ), In the second term, we decompose h 1 into two parts: h 1 = P 0 h 1 + Φ and use that P 0 = P 2 0 together with the fact that P 0 ∈ B(E, H) to obtain We can thus conclude the proof by taking ∆ 0 small enough and integrating the above differential inequality.
Notice that the inequality stated in the Proposition holds for the equivalent "hypocoercive norm" introduced above and thus also holds for the usual norm on H because of the equivalence (uniformly in ε) between those two norms. Combining estimates from Propositions 3.2 and 3.5, one can obtain that Corollary 3.6. Assume that h 0 ∈ E, h 1 ∈ H are such that For µ ∈ (0, µ ⋆ ) (where µ ⋆ is defined in (2.4)), for ∆ 0 small enough and for any δ > 0, we have that: is problematic when applying Gronwall lemma if one hopes to recover some decay in time but the extra small constant that appears in the estimate of (Id − P 0 )h 1 in (3.14) allows us to circumvent this difficulty.

Estimates on the kinetic problem.
Combining the previous corollary with Proposition 3.2, we are able to get our final a priori estimates on h in the space E: Proposition 3.9. Let r ∈ (0, 1). Assume that h 0 ∈ E, h 1 ∈ H are such that where ∆ 0 is small enough so that the conclusion of Corollary 3.6 holds. There exists ε † ∈ (0, ε 6 ), λ † ∈ (0, λ 6 ) (where ε 6 and λ 6 are defined in Proposition 3.8) such that for any ε ∈ (0, ε † ) and any λ 0 ∈ [0, λ † ) (where λ 0 is defined in Assumption 1.1), (1 − α(ε))/ε 2 has been defined in Theorem 2.1 and the constant C depends on r, ∆ 0 , µ ⋆ (defined in (2.1)) and ν 0 (defined in Lemma 2.3). Thanks to the above a priori estimates, we can prove Theorem 1.2 by introducing a suitable iterative scheme that is stable and convergent. We refer to [3,Section 5] for the details of the proof. We can actually prove the following more precise estimates (which will be useful in what follows) on h 0 ε and h 1 ε that are respectively solutions to (3.1) and (3. where we recall that the spaces H and H 1 are respectively defined in (2.2) and (2.3). Notice that in the previous inequalities, the multiplicative constants only involve quantities related to the initial data of the problem and are independent of ε.

DERIVATION OF THE FLUID LIMIT SYSTEM
The Cauchy theory developed in the previous results give all the a priori estimates that will allow to prove Theorem 1.4. To this end, we make additional assumptions in the definition of the spaces E and E 1 , namely, in this section, those spaces are defined through: We assume that Assumption 1.1 is met, consider ε, λ 0 and η 0 sufficiently small so that the conclusion of Theorem 1.2 holds in those spaces and consider {h ε } ε a family of solutions to (1.15) constructed in this theorem that splits as h ε = h 0 ε + h 1 ε with h 0 ε and h 1 ε defined in Section 3. We also fix T > 0 for the rest of the section.
4.1. Weak convergence. We start by the following lemma which in particular tells that the microscopic part of h ε vanishes in the limit ε → 0: Lemma 4.1. For any 0 t 1 t 2 T , there holds:

2)
where we recall that π 0 is the projection onto the kernel of L 1 defined in (1.23).

Incompressibility condition and Boussinesq relation.
Using (4.8) in the equations (4.7a)-(4.7b), using also that the restitution coefficient satisfies Assumption 1.1, we can easily obtain the incompressibility condition as well as the Boussinesq relation: div x u = 0 and ∇ x (̺ + ϑ 1 θ) = 0 (4.10) where we recall that ̺, u and θ are defined in (4.5). Using furthermore that the global mass of h ε vanishes (see (1.20)), we have that and thus that´T d ̺(t, x) dx = 0. It implies that we have the following strengthened Boussinesq relation: for almost every (t, x) ∈ (0, T ) × T d , (4.11) Remark 4.5. Notice here that the derivation of the strong Boussinesq relation ̺ + ϑ 1 θ = 0 is not as straightforward as in the elastic case. In the elastic case, the classical Boussinesq relation ∇ x (̺+ϑ 1 θ) = 0 straightforwardly implies the strong form of Boussinesq because the two functions ̺ and θ have zero spatial averages. This cannot be deduced directly in the granular context due to the dissipation of energy and we will see later on how to obtain it (see Proposition 4.8).

Equations of motion and temperature.
In order to identify the equations satisfied by u and θ, as in the elastic case, we start by studying the convergence of quantities that are related to (4.12) More precisely, we inverstigate the convergence of u ε := exp −t 1 − α(ε) ε 2 Pu ε and θ ε := 1 2 (|v| 2 − (d + 2)ϑ 1 )h ε where P is the Leray projection on divergence-free vector fields. Notice that if we compare our approach to the elastic case, we have added the exponential term in the definition of u ε in order to absorbe the term in the RHS in (4.7b). We compute the evolution of u ε and θ ε (by applying the Leray projector P to (4.7b) and by making an appropriate linear combination of (4.7a) and (4.7c)) and obtain: where A is defined in (4.6) and 14) The study of the limit ε → 0 in those equations is more favorable because compared to (4.7a)-(4.7b)-(4.7c), the gradient term in (4.7b) has been eliminated thanks to the Leray projector and also because A and b belong to the range of Id − π 0 so that thanks to Lemma 4.1, we know that the quantities ε −1 Div x A h ε and ε −1 div