Recent development in kinetic theory of granular materials: analysis and numerical methods

Jose Antonio Carrillo Jose Antonio Carrillo Department of Mathematics, Imperial College London SW7 2AZ London, UK Email: [email protected] , Jingwei Hu Jingwei Hu Department of Mathematics, Purdue University West Lafayette, IN 47907, USA Email: [email protected] , Zheng Ma Zheng Ma Department of Mathematics, Purdue University West Lafayette, IN 47907, USA Email: [email protected] and Thomas Rey Thomas Rey Univ. Lille, CNRS, UMR 8524, Inria – Laboratoire Paul Painlevé F-59000 Lille, France Email: [email protected]

Abstract.

Over the past decades, kinetic description of granular materials has received a lot of attention in mathematical community and applied fields such as physics and engineering. This article aims to review recent mathematical results in kinetic granular materials, especially for those which arose since the last review [91] by Villani on the same subject. We will discuss both theoretical and numerical developments. We will finally showcase some important open problems and conjectures by means of numerical experiments based on spectral methods.
Keywords: Granular gases equation, inelastic Boltzmann equation, fast spectral method, asymptotic behavior, review, hydrodynamic equations, granular flows, GPU.
2010 Mathematics Subject Classification: 82C40, 76P05, 76T25, 65N35.

JAC acknowledges support by the Engineering and Physical Sciences Research Council (EPSRC) under grant no. EP/P031587/1 and by the National Science Foundation (NSF) under grant no. RNMS11-07444 (KI-Net). JH was partially funded by NSF grant DMS-1620250 and NSF CAREER grant DMS-1654152. TR was partially funded by Labex CEMPI (ANR-11-LABX-0007-01) and ANR Project MoHyCon (ANR-17-CE40-0027-01).

1. The Boltzmann equation for granular gases

Granular gases have been initially introduced to describe the nonequilibrium behavior of materials composed of a large number of non necessarily microscopic particles, such as grains or sand. These particles form a gas, interacting via energy dissipating inelastic collisions. Although a statistical mechanics description of particle systems through inelastic collisions faces basic derivation problems as the inelastic collapse [72], i.e. infinite many collisions in finite time, the kinetic description of rapid granular flows [65, 66, 57] has been able to compute transport coefficients for hydrodynamic descriptions used successfully in situations that are a long way from their supposed limits of validity, to describe, for instance, shock waves in granular gases [82, 28], clustering [61, 32, 37], and the Faraday instability for vibrating thin granular layers [46, 73, 89, 88, 90, 28, 27, 36]. A large amount of practical systems can be described as a granular gas, such as for example spaceship reentry in a dusty atmosphere (Mars for instance), planetary rings [70, 7] and sorting behavior in vibrating layers of mixtures. A lot of other examples can be found in the thesis manuscript [40], and in the seminal book of Brilliantov and Pöschel [31].

Usually, a granular gas is composed of $10^{6}$ to $10^{16}$ particles. The study of such a system will then be impossible with a direct approach, and we shall adopt a kinetic point of view, studying the behavior of a one-particle distribution function $f$ , depending on time $t\geq 0$ , space $x\in\Omega\subset\mathbb{R}^{d_{x}}$ and velocity $v\in\mathbb{R}^{d}$ , for $d_{x}\leq d\in\{1,2,3\}$ . The statistical mechanics description of the system has been then admitted in the physical community as the tool to connect the microscopic description to macroscopic system of balance laws in rapid granular flows [65, 66, 57, 59, 30] as in the classical rarefied gases [39]. In this first section, we shall review some basics on the inelastic Boltzmann equation, and present the mathematical state of the art since the previous review paper on the subject [91].

Microscopic dynamics.

The microscopic dynamics can be summarized with the following hypotheses:

(1)

The particles interact via binary collisions. More precisely, the gas is rarefied enough so that collisions between 3 or more particles can be neglected.
(2)

These binary collisions are localized in space and time. In particular, all the particles are considered as point particles, even if they describe macroscopic objects.

(3)

Collisions preserve mass and momentum, but dissipate a fraction $1-e$ of the kinetic energy in the impact direction, where the inelasticity parameter $e\in[0,1]$ is called restitution coefficient:

(1)

\left\{\begin{aligned} &v^{\prime}\,+\,v_{*}^{\prime}\,=\,v\,+\,v_{*},\\ &|v^{\prime}|^{2}\,+\,|v_{*}^{\prime}|^{2}\,-\,|v|^{2}\,-\,|v_{*}|^{2}\,=\,-\frac{1-e^{2}}{2}\,|(v-v_{*})\cdot\omega|^{2}\,\leq\,0,\end{aligned}\right.

with $\omega\in\mathbb{S}^{d-1}$ being the impact direction. Using these conservation, one has the following two possible parametrizations (see also Fig. 2) of the post-collisional velocities, as a function of the pre-collisional ones:

•

The $\omega$ -representation or reflection map, given for $\omega\in\mathbb{S}^{d-1}$ by

	$\displaystyle v^{\prime}$	$\displaystyle=v-\frac{1+e}{2}\left((v-v_{*})\cdot\omega\right)\omega,$
(2)		$\displaystyle v_{*}^{\prime}$	$\displaystyle=v_{}+\frac{1+e}{2}\left((v-v_{})\cdot\omega\right)\omega.$

•

The $\sigma$ -representation or swapping map, given for $\sigma\in\mathbb{S}^{d-1}$ by

	$\displaystyle v^{\prime}$	$\displaystyle=\frac{v+v_{}}{2}+\frac{1-e}{4}(v-v_{})+\frac{1+e}{4}\|v-v_{*}\|\sigma,$
(3)		$\displaystyle v_{*}^{\prime}$	$\displaystyle=\frac{v+v_{}}{2}-\frac{1-e}{4}(v-v_{})-\frac{1+e}{4}\|v-v_{*}\|\sigma.$

Remark 1.

Taking $e=1$ in both (• ‣ 3) and (• ‣ 3) yields the classical energy-conservative elastic collision dynamics, as illustrated in Fig. 1.

Figure 1. Geometry of the inelastic collision in the physical space (green is elastic, red is inelastic).

The geometry of collisions is more complex than the classical elastic one. Indeed, fixing $v,v_{*}\in\mathbb{R}^{d}$ , denote by

\Omega_{\pm}:=\frac{v+v_{*}}{2}\pm\frac{1-e}{2}(v_{*}-v),\qquad O:=\frac{v+v_{*}}{2}=\frac{v^{\prime}+v_{*}^{\prime}}{2}.

Then if $u:=v-v_{*}$ is the relative velocity, one has

|\Omega_{-}-v^{\prime}|=|\Omega_{+}-v_{*}^{\prime}|=\frac{1+e}{2}|u|,

namely $v^{\prime}\in\mathcal{S}\left(\Omega_{-},|u|(1+e)/2\right)$ and $v_{*}^{\prime}\in\mathcal{S}\left(\Omega_{+},|u|(1+e)/2\right)$ , where $\mathcal{S}(x,r)$ is the sphere centered in $x$ and of radius $r$ (see also Fig. 2).

Refer to caption — Figure 2. Geometry of the inelastic collision in the phase space (dashed lines represent the elastic case).

Restitution coefficient.

The physics literature is quite divided on the question of whether the restitution coefficient $e$ should be a constant or not [31]. Although most of the early mathematical results on the topic consider a constant $e$ [91], it seems that this case is only realistic in dimension 1 of velocity (the so-called “collisional cannon” described in the chapter 4 of [31] is a famous counter-example). The true realistic case considers that $e$ depends on the relative velocity $|v-v_{*}|$ of the colliding particles. Even more precisely, it must be close to the elastic case $1$ for small relative velocities (namely no dissipation, elastic case), and decay towards $0$ when this relative velocity is large. The first mathematical result on this direction can be found in [85], where

(4)

e(|v-v_{*}|)=\dfrac{1}{1+c\,|v-v_{*}|^{\gamma}},

for a nonnegative constant $c$ characterizing the inelasticity strength ( $c=0$ being elastic), and $\gamma\in\mathbb{R}$ .

Another important case is the so-called viscoelastic hard spheres one, thoroughly studied mathematically in a series of papers [6, 4, 11, 5], where $e$ is given by the implicit relation

(5)

e(|v-v_{*}|)+a|v-v_{*}|^{1/5}e(|v-v_{*}|)^{3/5}=1,

for $a>0$ . More details on the derivation of this expression can be found in [31, Chapter 3].

Another quite rigorous study has been made in [81], with a threshold-dependent restitution coefficient:

e(r)=\left\{\begin{aligned} 1&&\text{ if }r<r_{*},\\ \bar{e}&&\text{ if }r\geq r_{*},\end{aligned}\right.

$\bar{e}<1$ and $r_{*}>0$ being fixed.

Finally, the case $e=0$ describes sticky collisions: the normal component of the kinetic energy being completely dissipated during impact, the particles stick and travel together in the tangent direction after impact. A derivation of the model from the microscopic dynamics on the line can be found in [29, 42].

Remark 2.

This model is meaningful even in dimension $1$ , which is not the case for elastic collisions. Indeed, such monodimensional collisions are only

\{v^{\prime},v_{*}^{\prime}\}=\{v,v_{*}\},

meaning that the particle velocities are either swapped or preserved. The particles being indistinguishable, nothing happens¹¹1Because of that, the elastic collision operator is simply equal to $0$ for a one-dimensional velocity space, the Boltzmann equation reducing only to the free transport equation.. In the $1$ d inelastic case, the collisions are given using (• ‣ 3) by

\{v^{\prime},v_{*}^{\prime}\}=\{v,v_{*}\}\quad\text{ or }\quad\left\{\frac{v+v_{*}}{2}\pm\frac{e}{2}(v-v_{*})\right\}

depending on the value of $\sigma\in\{\pm 1\}$ .

The granular gases operator: Weak form.

Using the microscopic hypotheses (1–2–• ‣ 3), one can derive the granular gases collision operator $\mathcal{Q}_{\mathcal{I}}$ , by following the usual elastic procedure (see [91] for more details). Its weak form in the $\sigma$ -representation is given by

(6)

\int_{\mathbb{R}^{d}}\mathcal{Q}_{\mathcal{I}}(f,f)(v)\,\psi(v)\,dv=\frac{1}{2}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{S}^{d-1}}f_{*}\,f\,\left(\psi^{\prime}+\psi_{*}^{\prime}-\psi-\psi_{*}\right)B(|v-v_{*}|,\cos\theta,E(f))\,d\sigma\,dv\,dv_{*},

where the collision kernel is typically of the form $B(|u|,\cos\theta,E(f))=\Phi(|u|)b(\cos\theta,E(f))$ , and $E(f)$ is the kinetic energy of $f$ , namely its second moment in velocity, the postcollisional velocities are computed by (• ‣ 3), and $\theta$ is the angle between $\sigma$ and $u$ . We shall assume in all the following of this section that the collision kernel is of generalized hard sphere type, namely

(7)

B(|u|,\cos\theta,E))=\Phi(|u|)\widehat{b}(\cos\theta,E)=|u|^{\lambda}\,b(\cos\theta)\,E^{\gamma},

where $\lambda\in[0,1]$ ( $\lambda=0$ being the simplified Maxwellian pseudo-molecules case and $\lambda=1$ the more relevant hard sphere case), $\gamma\in\mathbb{R}$ and the angular cross section $b$ verifies

(8)

0<\beta_{1}\leq b(x)\leq\beta_{2}<\infty,\quad\forall x\in[-1,1].

Remark 3.

Note that we assumed that the collision kernel $B$ in (6) depends on the relative velocity, the angle of collision, and on $E(f)$ . These former dependencies are quite classical, but the latter is not. Nevertheless, it makes a lot of sense physically speaking, as one can see in [83].

The weak form in the $\omega$ -representation can be written analogously as

(9)

\int_{\mathbb{R}^{d}}\mathcal{Q}_{\mathcal{I}}(f,f)(v)\,\psi(v)\,dv=\frac{1}{2}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{S}^{d-1}}f_{*}\,f\,\left(\psi^{\prime}+\psi_{*}^{\prime}-\psi-\psi_{*}\right)\widetilde{B}(|u|,\cos\theta,E(f))\,d\omega\,dv\,dv_{*},

where the postcollisional velocities are computed by (• ‣ 3), $\theta$ is the angle between $\omega$ and $u$ , and

\widetilde{B}(|u|,\cos\theta,E)=|u|^{\lambda}\widetilde{b}(\cos\theta)E^{\gamma}

with $\widetilde{b}(t)=3|t|b(1-2t^{2})$ for $-1\leq t\leq 1$ by the change of variables between the $\sigma$ - and the $\omega$ -representation, see [34] for details.

The granular gases operator: Strong form.

Deriving a strong form of $\mathcal{Q}_{\mathcal{I}}$ with the reflection map in the $\omega$ -representation is a matter of a change of variables. However, deriving a strong form of $\mathcal{Q}_{\mathcal{I}}$ is not as easy as in the elastic case in the $\sigma$ -representation since the collisional transform $(v,v_{*},\sigma)\to(v^{\prime},v_{*}^{\prime},\sigma)$ is not an involution and we have to go through the $\omega$ -representation, see [34] for details.

More precisely, given the restitution coefficient $e=e(|u|)$ depending on the relative velocity of the particles $u=v-v_{*}$ , we assume the collisional transform’s Jacobian for (• ‣ 3) is $J(|u|,\cos\theta)\neq 0$ for all $z$ . Notice $J=e$ in the constant restitution case. It is in general a complicated expression of the relative speed $r=|u|$ and $s=\cos\theta$ involving $e$ and its derivative. Then, the precollisional velocities read as

	$\displaystyle^{\prime}v$	$\displaystyle=v-\frac{1+e}{2e}\left((v-v_{*})\cdot\omega\right)\omega,$
(10)		${}^{\prime}v_{*}$	$\displaystyle=v_{}+\frac{1+e}{2e}\left((v-v_{})\cdot\omega\right)\omega.$

The final strong from of the operator is $\mathcal{Q}_{\mathcal{I}}(f,f)(v)=\mathcal{Q}_{\mathcal{I}}^{+}(f,f)(v)-f(v)\,L(f)(v)$ with the loss part of the operator given by

L(f)(v)=\int_{\mathbb{R}^{d}\times\mathbb{S}^{d-1}}\widetilde{B}(|v-v_{*}|,\cos\theta,E(f))f_{*}\,d\omega\,dv_{*}

and the gain part of the operator in strong form written as

(11)

\displaystyle\mathcal{Q}_{\mathcal{I}}^{+}(f,f)(v)

\displaystyle=\int_{\mathbb{R}^{d}\times\mathbb{S}^{d-1}}\widetilde{\Phi}_{e}^{+}(|u|,\cos\theta)\widetilde{b}_{e}^{+}(\cos\theta)\frac{E^{\gamma}}{J(|u|,\cos\theta)}\,\,^{\prime}f\,^{\prime}f_{*}\,d\omega\,dv_{*},

with $\widetilde{\Phi}_{e}^{+}(r,s)$ and $\widetilde{b}_{e}^{+}(s)$ given by

(12)

\widetilde{b}_{e}^{+}(s)=\widetilde{b}\left(\frac{s}{\sqrt{e^{2}+(1-e^{2})s^{2}}}\right)\ ,

and

(13)

\widetilde{\Phi}_{e}^{+}(r,s)=\Phi\left(\frac{r}{e}\sqrt{e^{2}+(1-e^{2})s^{2}}\right)=\left(\frac{r}{e}\sqrt{e^{2}+(1-e^{2})s^{2}}\right)^{\lambda}.

We can derive now the following strong form of the collision operator also in the $\sigma$ -representation by just changing variable in the operator from $\omega$ to $\sigma$ , see [34], to find the expressions of the loss and the gain terms in the $\sigma$ -representation:

L(f)(v)=\int_{\mathbb{R}^{d}\times\mathbb{S}^{d-1}}B(|v-v_{*}|,\cos\theta,E(f))f_{*}\,d\sigma\,dv_{*}

and

(14)

\displaystyle\mathcal{Q}_{\mathcal{I}}^{+}(f,f)(v)

\displaystyle=\int_{\mathbb{R}^{d}\times\mathbb{S}^{d-1}}\Phi_{e}^{+}(|u|,\cos\theta)b_{e}^{+}(\cos\theta)\frac{E^{\gamma}}{J(|u|,\cos\theta)}\,\,^{\prime}f\,^{\prime}f_{*}\,d\sigma\,dv_{*},

with $\Phi_{e}^{+}(r,s)$ and $b_{e}^{+}(s)$ given by

(15)

b_{e}^{+}(s)=b\left(\frac{(1+e^{2})s-(1-e^{2})}{(1+e^{2})-(1-e^{2})s}\right)\frac{\sqrt{2}}{\sqrt{{(1+e^{2})-(1-e^{2})s}}}\ ,

and

(16)

\Phi_{e}^{+}(r,s)=\Phi\left(\frac{r}{\sqrt{2}e}\sqrt{(1+e^{2})-(1-e^{2})s}\right)=\left(\frac{r}{\sqrt{2}e}\sqrt{(1+e^{2})-(1-e^{2})s}\right)^{\lambda}.

In these expressions, the precollisional velocities are given in the $\sigma$ -representation by

	$\displaystyle^{\prime}v$	$\displaystyle=\frac{v+v_{}}{2}+\frac{1-e}{4e}(v-v_{})+\frac{1+e}{4e}\|v-v_{*}\|\sigma,$
(17)		${}^{\prime}v_{*}$	$\displaystyle=\frac{v+v_{}}{2}-\frac{1-e}{4e}(v-v_{})-\frac{1+e}{4e}\|v-v_{*}\|\sigma.$

The granular gases collision operator has then the same structure of the elastic Boltzmann operator under Grad’s cutoff assumption, namely the difference between an inelastic gain term $\mathcal{Q}_{\mathcal{I}}^{+}(f,f)$ and a loss operator $f\,L(f)$ , which depends only on the chosen collision kernel, but not on the inelasticity.

We shall call the following granular gases equation, or the inelastic Boltzmann equation:

(18)

\frac{\partial f}{\partial t}+v\cdot\nabla_{x}f=\mathcal{Q}_{\mathcal{I}}(f,f).

We shall denote its first fluid moments (resp. density, mean momentum, and kinetic energy) by

\left(\rho(f),\,\rho(f)\bm{u}(f),\,E(f)\right):=\int_{\mathbb{R}^{d}}\left(1,v,|v|^{2}/2\right)f(v)\,dv.

Remark 4.

There is another popular approach to describe granular gases, which uses an Enskog-type collision operator. It is more relevant physically because it allows to keep the particles’ radii $\delta$ positive, hence delocalizing the collision²²2Note that using a BBGKY approach [50] to derive (11) is not expected to succeed, because among other problems the macroscopic size of the particles composing a granular gas is incompatible with the Boltzmann-Grad scaling assumption.. The microscopic hypothesis 2 is in particular not valid anymore. The strong form of the collision operator in the constant restitution coefficient case is given by

(19)

\mathcal{Q}_{\mathcal{E}}(f,f)(x,v)=\delta^{d-1}\,\int_{\mathbb{R}^{d}\times\mathbb{S}^{d-1}}\left(\widetilde{\Phi}_{e}^{+}(|u|,\cos\theta)\widetilde{b}_{e}^{+}(\cos\theta)\frac{G(\rho_{+})}{e}\,^{\prime}f_{+}\,^{\prime}f_{*}-G(\rho_{-})f_{-}\,f_{*}\right)\,d\omega\,dv_{*},

where $\rho$ is the local density of $f$ , $\pm$ denotes for a given function $g=g(x)$ the shorthand notation

g_{\pm}(x):=g(x\pm\delta\,\omega),

and $G$ is the local collision rate (also known as the correlation rate, see [91]). The global existence of renormalized solutions for the full granular gases equation (18) with the collision operator (19) has been established for both elastic and inelastic collisions in [45]. Existence and $L^{1}(dx\,dv)$ stability of such solutions has been proved in [95], for close to vacuum initial datum.

1.1. Cauchy theory of the granular gases equation

The space homogeneous setting.

Most of the rigorous mathematical results concerning the granular gases equation are obtained in the space homogeneous setting, where $f=f(t,v)$ is the solution to

(20)

\left\{\begin{aligned} &\partial_{t}f=\frac{\mathcal{Q}_{\mathcal{I}}(f,f)}{\varepsilon},\\ &f(0,v)=f_{in}(v),\end{aligned}\right.

for a given scaling parameter $\varepsilon>0$ .

The first existence results for solution to (20) can be found in [15, 16, 12, 13, 21, 38]. These works deal with the generalized Maxwellian pseudo-molecule kernel (7) $\lambda=0$ , $b\equiv 1$ and $\gamma=1/2$ , with a velocity dependent restitution coefficient $e=e(|v-v_{*}|)$ . Such a model allows to use some Fourier techniques to deal with the collision operator, altogether with the correct large time behavior for the kinetic energy, the so-called Haff’s cooling Law (26), and the correct hydrodynamic limit (31). The main result of these works is the global well-posedness of the solutions to (20) in $L_{2}^{1}(\mathbb{R}^{3})$ , the convergence towards equilibrium, and the contraction in different metrics for the equation. The proof relies in the careful study of the self-similar solutions to (20). Some extensions of these results, using the same collision kernel, can be found in the works [19, 18, 17], where in particular the uniqueness is obtained.

The physically relevant case of the hard sphere kernel $\lambda=1$ , $\gamma=0$ was first considered in [85] in the unidimensional case. This work establishes the global existence of measure solutions with finite kinetic energy for this problem. The proof relies on a priori estimates of the solution to (20) in the Monge-Kantorovich-Wasserstein metrics $W_{2}$ . This work also investigates the quasi-elastic $1-e^{2}\sim\varepsilon\to 0$ limit of the model, a nonlinear McNamara-Young-like friction equation. This friction model was later investigated in [71], where its global wellposedness was shown in the space of measure of finite energy.

The tail behavior of the equilibrium solution to the granular gases equation with a thermal bath $\Delta_{v}f$ was investigated in many papers, the main ones being [44, 20, 51]. They all proved the existence of non-gaussian, overpopulated tails for diffusively excited granular gases, namely:

Theorem 1 (From [20] and [44]).

Let $F(v)\geq 0$ for $v\in\mathbb{R}^{d}$ be a solution to the stationary equation

\mathcal{Q}_{\mathcal{I}}(F,F)+\Delta_{v}F=0

with all polynomial moments in velocity. Then,

F(v)\sim_{|v|\to\infty}\exp\left(-|v|^{\alpha}\right),

with $\alpha=1$ in the Maxwellian molecules case and $\alpha=3/2$ in the hard spheres case.

Indeed, the thermal bath gives an input of kinetic energy, preventing the appearance of trivial Dirac delta equilibria. The propagation of the Sobolev norms of the space dependent version of this equation was then established in [51]. It uses a careful estimate of the inelastic entropy production (23), and a fixed point argument for the existence and uniqueness of solutions.

Finally, the work [77] establishes the global well posedness of the granular gases equation without a thermal bath, for a general case of collision kernel (which contains (7)) and velocity dependent restitution coefficients:

Theorem 2 (Theorem 1.4 of [77]).

Let $0\leq f_{in}\in L^{1}_{3}\cap\in BV_{4}$ . Then for any $T\in(0,T_{c})$ , where $T_{c}:=\sup\,\{T>0:\mathcal{E}(f)(t)>0,\,\forall\,t<T\}$ is the so-called blowup time, there exists an unique nonnegative solution $f\in\mathcal{C}(0,T;L_{2}^{1})\cap L^{\infty}(0,T;L_{3}^{1})$ of (20). It preserves mass and momentum, and converges in the weak-* topology of measures towards a Dirac delta.

Their proof relies on careful estimates of the collision operator $\mathcal{Q}_{\mathcal{I}}$ in Orlicz space (specially the $L\log L$ space of finite entropy measures).

Remark 5.

The related (but still mostly open) problem of the propagation of chaos was considered in [78] for a very simplified inelastic collision operator with a thermal bath.

Cauchy problem in the space dependent setting.

The case of the space inhomogeneous setting³³3Physically more realistic, in part because of the spontaneous loss of space homogeneity that has been observed in [58]. has been much less investigated.

The first result can be found in [9] for the model introduced in [85] with a restitution coefficient given by (4), in one dimension of space and velocity. This work establishes the existence and uniqueness of mild (perturbative) solutions, first for small $L^{1}(dx\,dv)$ initial data, and then for compactly supported initial data. Their main argument is reminiscent from a work due to Bony in [23] concerning discrete velocity approximation of the Boltzmann equation in dimension $1$ .

The global existence of mild solutions in the general $\mathbb{R}_{x}^{3}\times\mathbb{R}_{v}^{3}$ setting, for a large class of velocity-dependent restitution coefficient, but for initial data close to vacuum, was obtained in [6]. The proof is based on a Kaniel-Shinbrot iteration on a very small functional space. The stability in $L^{1}(\mathbb{R}^{3}_{x}\times\mathbb{R}^{3}_{v})$ under the same assumptions was established in [94]. Finally the existence and convergence to equilibrium in $\mathbb{T}_{x}^{3}\times\mathbb{R}_{v}^{3}$ for a weakly inhomogeneous granular gas⁴⁴4Namely, the initial condition is chosen with a lot of exponential moments in velocity, and close to a space homogeneous profile. with a thermal bath was proved in [87], using a perturbative approach.

1.2. Large time behavior

Macroscopic properties of the granular gases operator.

Modeling-wise, the main microscopic difference between a granular gas and a perfect molecular gas is the dissipation of the kinetic energy. Using the weak form (6) among with the microscopic relations (1) of the inelastic collision operator, this yields

\int_{\mathbb{R}^{d}}\mathcal{Q}_{\mathcal{I}}(f,f)(v)\begin{pmatrix}1\\ v\\ |v|^{2}\end{pmatrix}dv\,=\,\begin{pmatrix}0\\ 0\\ -D(f)\end{pmatrix},

where $D(f)\geq 0$ is the energy dissipation functional, which depends only on the collision kernel:

(21)

D(f)\,:=\,\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}f\,f_{*}\,\Delta\left(|v-v_{*}|,E(f)\right)dv\,dv_{*}.

The quantity $\Delta\left(|u|,E\right)$ is the so-called energy dissipation rate, given using (1) by

(22)

\Delta\left(|u|,E\right)\,:=\,\frac{1-e^{2}}{4}\int_{\mathbb{S}^{d-1}}|u\cdot\omega|^{2}\,B(|u|,\cos\theta,E)\,d\omega\geq 0,\quad\forall e\in[0,1].

This dissipation of kinetic energy has a major consequence on the behavior of the solutions to the granular gases equation. Indeed, combined with the conservation of mass and momentum, it implies (at least formally) an explosive behavior, namely convergence in the weak-* topology of solutions to (18) towards Dirac deltas, centered in the mean momentum $\bm{u}$ :

f(t,\cdot)\rightharpoonup\delta_{v=\bm{u}},\quad t\to\infty.

As for the entropy, it is not possible to obtain any entropy dissipation for this equation, in order to precise this large time behavior. Indeed, as noticed in [51], taking $\psi(v)=\log f(v)$ in (6) yields

	$\displaystyle\int_{\mathbb{R}^{d}}\mathcal{Q}_{\mathcal{I}}(f,f)(v)\,\log f(v)\,dv=$	$\displaystyle\ \frac{1}{2}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{S}^{d-1}}f_{}\,f\log\left(\frac{f^{\prime}\,f_{}^{\prime}}{f\,f_{}}\right)B\,d\sigma\,dv\,dv_{}$
	$\displaystyle=$	$\displaystyle\ \frac{1}{2}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{S}^{d-1}}f_{}\,f\left[\log\left(\frac{f^{\prime}\,f_{}^{\prime}}{f\,f_{}}\right)-\frac{f^{\prime}\,f_{}^{\prime}}{f\,f_{}}+1\right]B\,d\sigma\,dv\,dv_{}$
(23)			$\displaystyle+\frac{1}{2}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{S}^{d-1}}\left(f_{}^{\prime}\,f^{\prime}-f_{}\,f\right)B\,d\sigma\,dv\,dv_{*}.$

The first term in (23), the elastic contribution, is nonpositive because $\log\lambda-\lambda+1\leq 0$ (this is Boltzmann’s celebrated H Theorem). Nevertheless, the second term has no sign a priori: it is 0 only in the elastic case (because of the involutive collisional transformation $(v,v_{*},\sigma)\to(v^{\prime},v_{*}^{\prime},\sigma)$ ). Boltzmann’s entropy

\mathcal{H}(f):=\int_{\mathbb{R}^{d_{v}}}f(v)\log f(v)\,dv

is then not dissipated by the solution of the granular gases equation if $e<1$ . Some work has been done on that direction in the numerical side. Indeed, adding a drift term or a thermal bath in velocity can yield numerical entropy dissipation, as noticed in [53].

Kinetic energy dissipation and the Haff’s cooling Law.

Let us assume in this subsection that the granular gas considered is space homogeneous, namely $f$ is solution to (20). Having no known entropy, one has then to use other macroscopic quantities to study the large time behavior of solutions to (18). Because of its explicit dissipation functional, kinetic energy is a good candidate for this. Moreover, being related to the variance, it allows to measure the concentration in velocity of the solution.

In order to have an explicit bound for the energy dissipation, let us assume that the collision kernel is of the general type (7). Using polar coordinates, it is straightforward to compute the dissipation rate (22):

(24)

\Delta\left(|u|,E\right)=b_{1}\,\frac{1-e^{2}}{4}\,|u|^{\lambda+2}E^{\gamma},

where thanks to (8)

b_{1}=\left|\mathbb{S}^{d-2}\right|\int_{0}^{\pi}\cos^{2}(\theta)\sin^{d-3}(\theta)\,b(\cos(\theta))\,d\theta<\infty.

Using the conservation of mass and momentum, one can always assume that the initial condition is of unit mass and zero momentum. Plugging (24) into (21) then yields using Hölder and Jensen inequalities

	$\displaystyle\frac{d}{dt}E(f)(t)$	$\displaystyle=-b_{1}\,\frac{1-e^{2}}{4}\,E(f)^{\gamma}(t)\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}f\,f_{}\,\|v-v_{}\|^{\lambda+2}\,dv\,dv_{*}$
		$\displaystyle\leq-b_{1}\,\frac{1-e^{2}}{4}\,E(f)^{\gamma}(t)\int_{\mathbb{R}^{d}}f(v)\,\|v\|^{\lambda+2}\,dv$
(25)			$\displaystyle\leq-b_{1}\,\frac{1-e^{2}}{4}\,E(f)^{1+\gamma+\lambda/2}(t).$

In particular, one will have the following large time behaviors: Setting $C_{e}=b_{1}\,\rho\,(1-e^{2})/4$ and $\alpha:=\gamma+1/2$ ,

•

Maxwellian pseudo-molecules ( $\lambda=\gamma=0$ ) decays exponentially fast towards the Dirac delta :

$E(f)(t)=E\left(f_{in}\right)e^{-C_{e}\,t};$

Notice that the inequality in (25) is an identity for this case.
•

Hard spheres ( $\lambda=1$ , $\gamma=0$ ) exhibits the seminal quadratic Haff’s cooling Law [60]:

(26) $E(f)(t)\leq{\left(E\left(f_{in}\right)^{-1/2}+C_{e}\,t/2\right)^{-2}}.$

•

Anomalous gases ( $\gamma\neq 0$ ) exhibits more general behaviors:

E(f)(t)\leq\left\{\begin{aligned} &{\left(E\left(f_{in}\right)^{\alpha}+C_{e}\alpha\,t\right)^{-\frac{1}{\alpha}}}&&\text{if }\gamma>{-1/2}\ (\alpha>0,\ \text{ finite time extinction)};\\ &E\left(f_{in}\right)e^{-C_{e}\,t}&&\text{if }\gamma={-1/2};\\ &\left(E\left(f_{in}\right)^{\alpha}-C_{e}\alpha\,t\right)^{-\frac{1}{\alpha}}&&\text{if }\gamma<{-1/2}\ (\alpha<0).\end{aligned}\right.

All of these formal results have been proven to be rigorous and sharp, with explicit lower bounds, in [15, 13] for the Maxwellian and hard sphere cases [74], and in [83] for the anomalous cases. Extension to the viscoelastic case can be found in [6, 5], where the energy is shown to behave as

E(f)(t)\sim_{t\to\infty}C\,(1+t)^{-5/3}.

All these papers share a common approach of proof, using the fact that the space homogeneous granular gases equation admits a self-similar behavior. Hence, introducing some well chosen time-dependent scaling function $\omega$ and $\tau$ , the distribution $f$ is written as

f(t,v)=\omega(t)^{d}g\left(\tau(t),\omega(t)\,v\right),

to take into account the concentration in the velocity variables⁵⁵5One can see the velocity scaling function $\omega$ as the inverse of the variance of the distribution $f$ . This scaling is then a continuous “zoom” on the blowup, and can be used to develop numerical methods for solving the full granular gases equation, see [49].. The rescaled function $g$ is then solution to the granular gases equation, with an anti-drift term in velocity:

\partial_{t}g+\nabla_{v}\cdot(v\,g)=\mathcal{Q}_{\mathcal{I}}(g,g).

Using some regularity estimates of the gain term of $\mathcal{Q}_{\mathcal{I}}$ “à la” Lions/Bouchut-Desvillettes [24] and some new Povzner-like estimates [3], [74] then obtains a lower bound for the energy of $g$ , yielding the generalized Haff’s law by coming back to $f$ .

Remark 6.

In the viscoelastic case, note that the rescaling in velocity induces a time dependency on the restitution coefficient, complicating the proof of the Haff’s cooling law [4]. It is also the case in the anomalous setting, where the rescaling function depends nonlinearly on the solution $f$ .

The question of the uniqueness, stability and exponential return to an universal equilibrium profile (hypocoercivity, see [92]) of the self-similar solutions has then been fully addressed in the series of work [75, 76], for a constant restitution coefficient, with and without a thermal bath. Extension of these results to the viscoelastic case has been done in [5].

Remark 7.

These results are reminiscent of the works [8, 35, 12, 13, 21, 22], where the exponential convergence to equilibrium of the solution to the nonlocal granular media equation or the Maxwellian case has been shown in the $W_{2}$ or Fourier metrics.

1.3. Compressible hydrodynamic limits

Let us consider in this subsection the following hyperbolic scaling of the granular gases equation:

(27)

\frac{\partial f_{\varepsilon}}{\partial t}+v\cdot\nabla_{x}f_{\varepsilon}=\frac{1}{\varepsilon}\mathcal{Q}_{\mathcal{I}}(f_{\varepsilon},f_{\varepsilon}).

Determining the precise hyperbolic limit $\varepsilon\to 0$ of equation (27) is a fundamental, yet very difficult question.

Indeed, for the elastic case $e=1$ , one simply has to use the fact that the equilibria of the collision operator are at the thermodynamical equilibrium (gaussian distributions) and the conservation of mass, momentum and kinetic energy to obtain the classical compressible Euler-Fourier system [39]. Because of the trivial Dirac equilibria, this question is more intricate for the true inelastic case.

Pressureless Euler dynamics.

Adopting the same approach as in the elastic case, one can formally plug the “equilibrium” Dirac deltas in the pressure to obtain the following pressureless Euler system:

(28)

\left\{\begin{aligned} &\frac{\partial\rho}{\partial t}\,+\,\nabla\cdot(\rho\,\bm{u})\,=\,0,\\ &\frac{\partial(\rho\bm{u})}{\partial t}\,+\,\nabla\cdot\left(\rho\,\bm{u}\otimes\bm{u}\right)\,=\,\bm{0}.\end{aligned}\right.

This system can describe various interesting physical situations, such as galactic clusters, but is notoriously difficult to study mathematically. Its solution are in general ill-posed, as classical solutions cannot exists for large times and weak solutions are not unique.

In the unidimensional case, it is however possible to recover a well posed theory by imposing a semi-Lipschitz condition on $u$ . This theory was introduced in [25], and later extended in [26] and [63]. We cite below the main result of [63], where $M^{1}(\mathbb{R})$ denotes the space of Radon measures on $\mathbb{R}$ and $L^{2}(\rho)$ for $\rho\geq 0$ in $M^{1}(\mathbb{R})$ denotes the space of functions which are square integrable against the density $\rho$ .

Theorem 3 (From [63]).

For any $\rho^{0}\geq 0$ in $M^{1}(\mathbb{R})$ and any $u^{0}\in L^{2}(\rho^{0})$ , there exists $\rho\in L^{\infty}(\mathbb{R}_{+},M^{1}(\mathbb{R}))$ and $u\in L^{\infty}(\mathbb{R}_{+},L^{2}(\rho))$ solution to (28) in the sense of distribution and satisfying the semi-Lipschitz Oleinik-type bound

(29)

u(t,x)-u(t,y)\leq\frac{x-y}{t},\quad\mbox{for}\ a.e.\;x>y.

Moreover the solution is unique if $u^{0}$ is semi-Lipschitz or if the kinetic energy is continuous at $t=0$

\int_{\mathbb{R}}\rho(t,dx)\,|u(t,x)|^{2}\longrightarrow\int_{\mathbb{R}}\rho^{0}(dx)\,|u^{0}(x)|^{2},\qquad\mbox{as}\ t\rightarrow 0.

The proof of Th. 3 is quite delicate, relying on duality solutions. For this reason, we only explain the rational behind the bound (29), which can be seen very simply from the discrete sticky particles dynamics. We refer in particular to [29] for the limit of this sticky particles dynamics as $N\rightarrow\infty$ .

Consider $N$ particles on the real line. We describe the $i^{\text{th}}$ particle at time $t>0$ by its position $x_{i}(t)$ and its velocity $v_{i}(t)$ . Since we are dealing with a one dimensional dynamics, we can always assume the particles to be initially ordered

x_{1}^{in}<x_{2}^{in}<\ldots<x_{N}^{in}.

The dynamics is characterized by the following properties

(1)

The particle $i$ moves with velocity $v_{i}(t)$ : $\frac{d}{dt}x_{i}(t)=v_{i}(t)$ .
(2)

The velocity of the $i^{\text{th}}$ particle is constant, as long as it does not collide with another particle: $v_{i}(t)$ is constant as long as $x_{i}(t)\neq x_{j}(t)$ for all $i\neq j$ .
(3)

The velocity jumps when a collision occurs: if at time $t_{0}$ there exists $j\in\{1,\ldots,N\}$ such that $x_{j}(t_{0})=x_{i}(t_{0})$ and $x_{j}(t)\neq x_{i}(t)$ for any $t<t_{0}$ , then all the particles with the same position take as new velocity the average of all the velocities

$v_{i}(t_{0}+)=\frac{1}{|{j|x_{j}(t_{0})=x_{i}(t_{0})}|}\sum_{j|x_{j}(t_{0})=x_{i}(t_{0})}v_{j}(t_{0}-).$

Note in particular that particles having the same position at a given time will then move together at the same velocity. Hence, only a finite number of collisions can occur, as the particles aggregates.

This property also leads to the Oleinik regularity. Consider any two particles $i$ and $j$ with $x_{i}(t)>x_{j}(t)$ . Because they occupy different positions, they have never collided and hence $x_{i}(s)>x_{j}(s)$ for any $s\leq t$ . If neither had undergone any collision then $x_{i}(0)=x_{i}(t)-v_{i}(t)\,t>x_{j}(0)=x_{j}(t)-v_{j}(t)\,t$ or

(30)

\frac{\left(v_{i}-v_{j}\right)_{+}}{\left(x_{i}-x_{j}\right)_{+}}<\frac{1}{t},

where $x_{+}:=\max(x,0)$ . It is straightforward to check that (30) still holds if particles $i$ and $j$ had some collisions between time $0$ and $t$ .

As one can see this bound is a purely dispersive estimate based on free transport and the exact equivalent of the traditional Oleinik regularization for Scalar Conservation Laws, see [80]. It obviously leads to the semi-Lipschitz bound (29) as $N\rightarrow\infty$ .

This result was extended to the one dimensional (in space and velocity) granular gases equation (27) in [64]. The basic idea of the proof of this work is to compare the granular gases dynamics to the pressureless gas system (28). The main difficulty is to show that $f_{\varepsilon}$ becomes monokinetic at the limit (see also the very recent paper [68]). This is intimately connected to the Oleinik property (29), just as this property is critical to pass to the limit from the discrete sticky particles dynamics.

Unfortunately it is not possible to obtain (29) directly. Contrary to the sticky particles dynamics, this bound cannot hold for any finite $\varepsilon$ (or for any distribution that is not monokinetic). This is the reason why it is very delicate to obtain the pressureless gas system from kinetic equations (no matter how natural it may seem). Indeed we are only aware of one other such example in [69].

One of the main contributions of [64] is a complete reworking of the Oleinik estimate, still based on dispersive properties of the free transport operator $v\,\partial_{x}$ but compatible with kinetic distributions that are not monokinetic, through the introduction of a new, global nonlinear energy. The main result in this paper is the following:

Theorem 4 (from [64]).

Consider a sequence of weak solutions $f_{\varepsilon}\in L^{\infty}([0,\ T],\ L^{p}(\mathbb{R}^{2}))$ for some $p>2$ and with total mass $1$ to the granular gases Eq. (27) such that all initial $v$ -moments are uniformly bounded in $\varepsilon$ , some moment in $x$ is uniformly bounded, and $f^{0}_{\varepsilon}$ is, uniformly in $\varepsilon$ , in some $L^{p}$ for $p>1$ . Then any weak-* limit $f$ of $f_{\varepsilon}$ is monokinetic, $f=\rho(t,x)\,\delta(v-u(t,x))$ for $a.e.\ t$ , where $\rho,\;u$ are a solution in the sense of distributions to the pressureless system (28) while $u$ has the Oleinik property (29).

Quasi-elastic limit.

The physical community usually considers another approach, that is assuming that the granular gas is in a quasi-elastic $1-e^{2}\sim\varepsilon\rightarrow 0$ setting. This was first proposed in [65, 66], using an approach very similar to the seminal Grad’s 13 moments closure for rarefied gas dynamics. The difficulty of a hydrodynamic description of granular materials has been addressed in well reasoned terms in [55, 56], and as already discussed in the introduction, the hydrodynamic equations obtained with the kinetic theory of granular gases have been shown to be insightful well beyond their supposed limit of validity, i.e., away from the quasi-elastic limit assumption with external sources of energy. In fact, assuming that solutions of the kinetic problem do not deviate from being Gaussians, one can then obtain in the hard sphere case the following quasi-elastic compressible Euler system

(31)

\left\{\begin{aligned} &\frac{\partial\rho}{\partial t}\,+\,\nabla\cdot(\rho\,\bm{u})\,=\,0,\\ &\frac{\partial(\rho\bm{u})}{\partial t}\,+\,\nabla\cdot\left(\rho\,(\bm{u}\otimes\bm{u})\,+\,\rho\,T\,{\rm\bf I}\right)\,=\,\bm{0},\\ &\frac{\partial W}{\partial t}\,+\,\nabla\cdot\left(\bm{u}\,(W\,+\,\rho\,T)\right)\,=\,-K\,\rho\,T^{3/2},\end{aligned}\right.

representing the evolution of number density of particles $\rho(t,x)$ , velocity field $\bm{u}(t,x)$ and the total energy $W$ , which is given by $W=\frac{3}{2}\rho T+\frac{1}{2}\rho|\bm{u}|^{2}=\rho\epsilon+\frac{1}{2}\rho|\bm{u}|^{2}$ , with $T$ being the granular temperature (3D). Here $K=K(d,B)$ is an explicit nonnegative constant, see below. This is a compressible Euler-type system, which dissipates the kinetic energy thanks to its nonzero right hand side. The particular expression of this RHS allows, after integration in space, to recover the correct Haff’s cooling Law.

The assumption that the solutions are not far from Gaussians obviously degenerates in a free cooling granular gas at some point leading to the so-called clustering instability studied by means of (31), see for instance [37] and the references therein. In fact, this assumption can be shown to be valid in the quasi-elastic limit, see [76] for a rigorous justification of this property. Physicists argue that this assumption is also generically true in practical experiments with external sources of energy such as the shock waves in granular flows under gravity [82, 28], clustering [61, 32, 37], the Faraday instability for vibrating thin granular layers [46, 73, 89, 88, 90, 28, 27, 36], and many other applications, see [36, 67] and the references therein.

Passing from the granular gases equation (18) to (31) has not been established properly. It can still be done formally under the weak inelasticity hypothesis $1-e^{2}\sim\varepsilon$ , see [86]. This particular scaling insures that the granular gases operator converges towards the elastic Boltzmann operator, as was shown rigorously in [75] in the space homogeneous setting. Moreover, it allows to characterize the equilibrium distribution of the limit operator, which is Gaussian.

A first step towards a rigorous compressible hydrodynamic limit is available in [84], where the study of the spectrum of the heated granular gases operator $\mathcal{Q}_{\mathcal{I}}+\tau\Delta_{v}$ , linearized with respect to the equilibrium described in [76], is done. For small inelasticity $1-e^{2}\sim 0$ , this work provides a spectral decomposition, and more importantly the existence and computation of eigenvalue branches. This work generalizes the seminal paper [43] on the spectrum of the linearized Boltzmann operator in $L^{2}$ to the $L^{1}$ and inelastic setting.

Other types of fluid limits (such as viscous limits) of the granular gases equation has been described in the review paper [41] and in the recent survey [54] for many different physical regimes, but none has been rigorously established. To illustrate the kind of equations obtained through these procedure, we write the generalized Navier-Stokes compressible equations for granular media in conservative form, see [36], as

(32)

\begin{array}[]{l l}&{\displaystyle\frac{\partial\rho}{\partial t}}+\nabla\cdot(\rho\bm{u})=0\vspace{6pt}\\ &{\displaystyle\frac{\partial(\rho\bm{u})}{\partial t}}+\nabla\cdot\left[\rho\,(\bm{u}\otimes\bm{u})\right]=\nabla\cdot\bm{P}+\rho\,\bm{F}\vspace{6pt}\\ &{\displaystyle\frac{\partial W}{\partial t}}+\nabla\cdot\left[\bm{u}W\right]=-\nabla\cdot q+\bm{P}:\bm{E}+\bm{u}\cdot(\nabla\cdot\bm{P})-\gamma+\rho\bm{u}\cdot\bm{F},\end{array}

representing the evolution of number density of particles $\rho(t,x)$ , velocity field $\bm{u}(t,x)$ and the total energy $W$ , which is given by $W=\frac{3}{2}\rho T+\frac{1}{2}\rho|\bm{u}|^{2}=\rho\epsilon+\frac{1}{2}\rho|\bm{u}|^{2}$ , with $T$ being the granular temperature (3D). The symbol $F$ stands for external forces applied to the system. The constitutive relations for the momentum and heat fluxes write, as usual,

P_{ij}=\left[-p+\left(\lambda-\frac{2}{3}\mu\right)\sum_{i}E_{ii}\right]\delta_{ij}+2\mu E_{ij}

for the stress tensor, with $E_{ij}=\frac{1}{2}\left(\frac{\partial U_{i}}{\partial x_{j}}+\frac{\partial U_{j}}{\partial x_{i}}\right)$ . The thermal conductivity relates linearly the heat flux $q$ to the temperature gradient, $q=-\chi\nabla T$ .

The equation of state is relevant here: we use the expression $G(\nu)=[1-(\nu/\nu_{max})^{\frac{4}{3}\nu_{max}}]^{-1}$ for the contact value of the pair correlation function $G(\nu)$ , which accounts for high densities, and the equation of state $p=\rho T(1+2(1+e)\nu G(\nu))$ , where $\nu$ is the packing fraction and $e$ the constant restitution coefficient. Random close-packing is achieved in 3D at $\nu_{max}=0.65$ ; we do not allow any packing fraction higher than 99.99% of this value. The transport coefficients for constant restitution given in [65, 14] write,

\gamma=\frac{12}{\sigma\sqrt{\pi}}(1-e^{2})\rho T^{3/2}G(\nu),

for the cooling coefficient $\gamma$ , which models energy dissipation through collisions. Other kinetic coefficients are the shear viscosity,

\mu=\frac{\sqrt{\pi T}}{6}\rho\sigma\left[\frac{5}{16G(\nu)}+1+\frac{4}{5}\left(1+\frac{12}{\pi}\right)G(\nu)\right],

the bulk viscosity,

\lambda=\frac{8}{3\sqrt{\pi}}\rho\sigma\sqrt{T}G(\nu)

and the thermal conductivity

\chi=\frac{15}{16}\sqrt{\pi T}\rho\sigma\left[\frac{5}{24G(\nu)}+1+\frac{6G(\nu)}{5}\left(1+\frac{32}{9\pi}\right)\right],

implemented directly with no further assumptions.

A first attempt to derive equations of compressible Navier-Stokes type was done in the paper [33] using singular perturbations of the collision operator $\mathcal{Q}_{\mathcal{I}}$ and a central manifold approach inspired from [47]. The fact is that transport coefficients for compressible Navier-Stokes like equations can be derived by moment closures under different assumptions and these equations are able to recover realistic phenomena in granular gases, see [54] for a very recent review. For instance, the Faraday instability for vertically oscillated granular layer was analysed in [36, 2] comparing molecular dynamics simulations, and different closures proposed in the literature. The conclusion, as the reader can check in the numerical experiments in [1], is that the qualitative behavior of the numerical experiments using equations (32) fully recovers the physical expected outcome. Thus, this is a physical validation of these kind of approximations that are totally opened for mathematicians to be derived in full rigor.

2. Fourier spectral methods for the inelastic Boltzmann collision operator

Due to the complexity of the inelastic Boltzmann collision operator $\mathcal{Q}_{\mathcal{I}}(f,f)$ (a five-fold nonlinear integral operator), numerical simulation of granular gases is challenging and mostly done at the particle level (direct simulation Monte Carlo method [10] or its variants). Over the past decade, a class of deterministic numerical methods – the Fourier-Galerkin spectral method – has received a lot of popularity for its high accuracy and relatively low computational cost. The first attempt was made in [79] for the one-dimensional model. Later in [48, 52, 49], both two and three dimensional cases were considered. Although the implementation details may differ, the essential ideas in these works are the same, that is, utilizing the translational invariance of the collision operator and convolution property of the Fourier basis to write the collision operator as a weighted convolution in the Fourier space. In this way, the $O(N^{3d})$ cost per evaluation of the collision operator in the Galerkin framework (since $\mathcal{Q}_{\mathcal{I}}(f,f)$ is quadratic) is readily reduced to $O(N^{2d})$ , where $N$ is the number of basis used in each velocity dimension. Even though this reduction is dramatic compared to other spectral basis, numerical implementation of the “direct” Fourier spectral method is still computationally demanding; what makes it worse is that the method also requires $O(N^{2d})$ memory to store the precomputed weights, which quickly becomes a bottleneck as $N$ increases. Recently, a fast Fourier spectral method was introduced in [62], wherein the key idea is to shift some offline precomputed items to online computation so that the weighted convolution in the original method can be rendered into a few pure convolutions to be evaluated efficiently by the fast Fourier transform (FFT). As a result, both the computational complexity and the memory requirement in the direct Fourier method are reduced.

In this section, we briefly review the original direct Fourier spectral method proposed in [48] and then its fast version introduced in [62]. To this end, let us introduce the following weak form of the inelastic Boltzmann collision operator $\mathcal{Q}_{\mathcal{I}}(f,f)$ using the $\sigma$ -representation (• ‣ 3):

(33)

\int_{\mathbb{R}^{d}}\mathcal{Q}_{\mathcal{I}}(f,f)(v)\psi(v)\mathop{}\!\mathrm{d}{v}=\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}\int_{\mathbb{S}^{d-1}}B(|v-v_{*}|,\sigma\cdot\widehat{(v-v_{*})})ff_{*}\left(\psi^{\prime}-\psi\right)\mathop{}\!\mathrm{d}{\sigma}\mathop{}\!\mathrm{d}{v}\mathop{}\!\mathrm{d}{v_{*}}\,,

where $\sigma\cdot\widehat{(v-v_{*})}=\cos\theta$ ( $\hat{u}$ denotes the unit vector along $u$ ), and for simplicity we assume $B$ does not depend on the kinetic energy $E$ (compare with (7)). From the numerical point of view, this does not impose any limitation since the $E$ part can always be factored out from the integral sign.

2.1. The direct Fourier spectral method

We first perform a change of variable $v_{*}\rightarrow g=v-v_{*}$ in (33) to obtain

\int_{\mathbb{R}^{d}}\mathcal{Q}_{\mathcal{I}}(f,f)(v)\psi(v)\mathop{}\!\mathrm{d}{v}=\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}\int_{\mathbb{S}^{d-1}}B(|g|,\sigma\cdot\hat{g})f(v)f(v-g)\left(\psi(v^{\prime})-\psi(v)\right)\mathop{}\!\mathrm{d}{\sigma}\mathop{}\!\mathrm{d}{g}\mathop{}\!\mathrm{d}{v},

where

v^{\prime}=v-\frac{1+e}{4}(g-|g|\sigma).

We then assume that $f$ has a compact support: $\text{Supp}_{v}(f)\approx B_{S}$ , where $B_{S}$ is a ball centered at the origin with radius $S$ . Hence it suffices to truncate the infinite integral in $g$ to a larger ball $B_{R}$ with radius $R=2S$ :

(34)

\int_{\mathbb{R}^{d}}\mathcal{Q}_{\mathcal{I}}(f,f)(v)\psi(v)\mathop{}\!\mathrm{d}{v}=\int_{\mathbb{R}^{d}}\int_{B_{R}}\int_{\mathbb{S}^{d-1}}B(|g|,\sigma\cdot\hat{g})f(v)f(v-g)\left(\psi(v^{\prime})-\psi(v)\right)\mathop{}\!\mathrm{d}{\sigma}\mathop{}\!\mathrm{d}{g}\mathop{}\!\mathrm{d}{v}.

Next we restrict $v$ to a cubic computational domain $D_{L}=[-L,L]^{d}$ , and approximate $f$ by a truncated Fourier series:

f(v)\approx f_{N}(v)=\sum_{k=-\frac{N}{2}}^{\frac{N}{2}-1}\hat{f}_{k}e^{i\frac{\pi}{L}k\cdot v},\quad\hat{f}_{k}=\frac{1}{(2L)^{d}}\int_{D_{L}}f(v)e^{-i\frac{\pi}{L}k\cdot v}\mathop{}\!\mathrm{d}{v}.

Here $k=(k_{1},\dots,k_{d})$ is a multidimensional index, and $\sum_{k={-N/2}}^{N/2-1}:=\sum_{k_{1},\dots,k_{d}=-N/2}^{N/2-1}$ . The choice of $L$ should be chosen at least as $L\geq(3+\sqrt{2})S/2$ to avoid aliasing, see [48] for more details. Now substituting $f_{N}$ into (34) and choosing $\psi(v)=e^{-i\frac{\pi}{L}k\cdot v}$ , we can obtain the $k$ -th mode of the collision operator as

(35)

\hat{Q}_{k}=\sum_{\begin{subarray}{c}l,m=-\frac{N}{2}\\ l+m=k\end{subarray}}^{\frac{N}{2}-1}G(l,m)\hat{f}_{l}\hat{f}_{m},

where the weight $G(l,m)$ is given by

G(l,m)=\int_{B_{R}}e^{-i\frac{\pi}{L}m\cdot g}\left[\int_{\mathbb{S}^{d-1}}B(|g|,\sigma\cdot\hat{g})\left(e^{i\frac{\pi}{L}\frac{1+e}{4}(l+m)\cdot(g-|g|\sigma)}-1\right)\mathop{}\!\mathrm{d}{\sigma}\right]\mathop{}\!\mathrm{d}{g}.

In the original spectral method [48], the weight $G(l,m)$ is precomputed and stored since it is independent of the solution $f$ which leads to a memory requirement of $O(N^{2d})$ . During the online computation, the weighted sum (35) is directly evaluated whose complexity is $O(N^{2d})$ .

2.2. The fast Fourier spectral method

To reduce the complexity of the direct spectral method as well as to alleviate its memory requirement, the key idea introduced in [62] is to render the weighted convolution (35) into a pure convolution so that it can be computed efficiently by the FFT. One way to achieve this is through a low-rank approximation of $G(l,m)$ , namely,

(36)

G(l,m)\approx\sum_{p=1}^{N_{p}}\alpha_{p}(l+m)\beta_{p}(m),

where $\alpha_{p}$ and $\beta_{p}$ are some functions to be determined and the number of terms $N_{p}$ in the expansion is small. Then (35) becomes

(37)

\hat{Q}_{k}\approx\sum_{p=1}^{N_{p}}\alpha_{p}(k)\sum_{\begin{subarray}{c}l,m=-\frac{N}{2}\\ l+m=k\end{subarray}}^{\frac{N}{2}-1}\hat{f}_{l}\left(\beta_{p}(m)\hat{f}_{m}\right),

where the inner summation is a pure convolution of two functions $\hat{f}_{l}$ and $\beta_{p}(m)\hat{f}_{m}$ . Hence the total complexity to evaluate $\hat{Q}_{k}$ (for all $k$ ) is brought down to $O(N_{p}N^{d}\log N)$ , i.e., a few number of FFTs.

Specifically, we first split $G(l,m)$ into a gain term and a loss term:

G(l,m)=G_{\text{gain}}(l,m)-G_{\text{loss}}(m),

where

G_{\text{gain}}(l,m):=\int_{B_{R}}e^{-i\frac{\pi}{L}m\cdot g}\left[\int_{\mathbb{S}^{d-1}}B(|g|,\sigma\cdot\hat{g})e^{i\frac{\pi}{L}\frac{1+e}{4}(l+m)\cdot(g-|g|\sigma)}\mathop{}\!\mathrm{d}{\sigma}\right]\mathop{}\!\mathrm{d}{g},

and

G_{\text{loss}}(m):=\int_{B_{R}}e^{-i\frac{\pi}{L}m\cdot g}\left[\int_{\mathbb{S}^{d-1}}B(|g|,\sigma\cdot\hat{g})\mathop{}\!\mathrm{d}{\sigma}\right]\mathop{}\!\mathrm{d}{g}.

Note that the loss term is readily a function of $m$ , hence no approximation/decomposition is actually needed. This suggests to evaluate the loss term of the collision operator by a precomputation of $G_{\text{loss}}(m)$ and then compute

\hat{Q}^{-}_{k}=\sum_{\begin{subarray}{c}l,m=-\frac{N}{2}\\ l+m=k\end{subarray}}^{\frac{N}{2}-1}\hat{f}_{l}\left(G(m)\hat{f}_{m}\right)

by FFT. For the gain term, to get a decomposition of form (36), we introduce a quadrature rule to discretize $g$ , then $G_{\text{gain}}(l,m)$ can be approximated as

(38)

G_{\text{gain}}(l,m)\approx\sum_{\rho,\hat{g}}w_{\rho}w_{\hat{g}}\,\rho^{d-1}e^{-i\frac{\pi}{L}\rho\,m\cdot\hat{g}}F(l+m,\rho,\hat{g}),

where $\rho:=|g|\in[0,R]$ is the radial part of $g$ and $\hat{g}\in\mathbb{S}^{d-1}$ is the angular part, and $w_{\rho}$ and $w_{\hat{g}}$ are the corresponding quadrature weights. The function $F$ is given by

(39)

F(l+m,\rho,\hat{g}):=\int_{\mathbb{S}^{d-1}}B(\rho,\sigma\cdot\hat{g})e^{i\frac{\pi}{L}\rho\frac{1+e}{4}(l+m)\cdot(\hat{g}-\sigma)}\mathop{}\!\mathrm{d}{\sigma}.

With this approximation, the gain term of the collision operator can be evaluated as

\hat{Q}_{k}^{+}\approx\sum_{\rho,\hat{g}}w_{\rho}w_{\hat{g}}\,\rho^{d-1}\,F(k,\rho,\hat{g})\sum_{\begin{subarray}{c}l,m=-\frac{N}{2}\\ l+m=k\end{subarray}}^{\frac{N}{2}-1}\hat{f}_{l}\left[e^{-i\frac{\pi}{L}\rho\,m\cdot\hat{g}}\hat{f}_{m}\right],

which is in the same form as explained in (37).

As for the quadratures, the radial direction $\rho$ can be approximated by the Gauss-Legendre quadrature. Since the integrand in (38) is oscillatory on the scale of $O(N)$ , the number of quadrature points needed for $\rho$ should be $O(N)$ . The angular direction in $2$ D can be discretized using simple rectangular rule which is expected to yield spectral accuracy due to the periodicity. While in 3D, we choose to use the spherical design [93] which is the near optimal quadrature on the sphere.

To summarize, the total complexity to evaluate $\hat{Q}_{k}$ is $O(MN^{d+1}\log N)$ , where $M$ is the number of points used on $\mathbb{S}^{d-1}$ and $M\ll N^{d-1}$ . Furthermore, the only quantity that needs to be precomputed and stored is (39), which in the worst scenario only requires $O(MN^{d+1})$ memory.

3. Numerical experiments and results

The accuracy and efficiency of the fast spectral method has been validated in [62]. In this section, we perform some additional tests to demonstrate the potential of the method in predicting some mathematical theories. We also introduce a GPU implementation of the method that significantly improves the CPU version in [62]. This is critical especially for long time simulation.

We consider the following spatially homogeneous equation with a heat bath:

(40)

\partial_{t}f=\mathcal{Q}_{\mathcal{I}}(f,f)+\tau\Delta_{v}f,

where $\tau$ is the parameter describing the strength of the heat bath. Notice that it is not necessarily related to the inelasticity parameter $e$ , contrarily to e.g. [76]. The heat bath $\Delta_{v}f$ will also be discretized using the Fourier spectral method and Runge-Kutta method is used for time marching.

For the collision operator, we consider the simplified variable hard sphere kernel

(41)

B(|g|,\sigma\cdot\hat{g},E)=C_{\lambda}|g|^{\lambda},\quad 0\leq\lambda\leq 1,

where $C_{\lambda}>0$ is some constant (namely (7) with $b=C_{\lambda}$ and $\gamma=0$ .

For Maxwell molecules, given the initial condition $f_{0}(v)$ whose macroscopic quantities are $\rho_{0}$ , $u_{0}$ and $T_{0}$ , the density and velocity are conserved so $\rho(t)=\rho_{0}$ , $u(t)=u_{0}$ and the temperature will evolve as

(42)

T(t)=\left(T_{0}-\frac{8\tau}{1-e^{2}}\right)\exp{\left(-\frac{\rho_{0}(1-e^{2})}{4}t\right)}+\frac{8\tau}{1-e^{2}},

We could see

\lim\limits_{t\rightarrow\infty}T(t)=\frac{8\tau}{1-e^{2}}.

As in [62], this analytical formula of temperature works as the reference solution to ensure the correctness of the numerical solution.

3.1. GPU parallelized implementation

From the implementation perspective, we dramatically improve the efficiency of the fast spectral method by using GPU via Nivida’s CUDA. GPU-parallelized implementations are run on $2$ Intel^® Xeon^® Silver 4110 2.10 GHz CPUs with $4$ NVIDIA Geforce GTX 2080 Ti (Turing) GPUs accompanying CUDA driver 10.0 and CUDA runtime 10.0. The operating system used is 64-bit Unbuntu 18.04. The CPU has 8 cores, 16 threads with max turbo frequency 3.00 Ghz equiped with 128GB DDR4 REG ECC memory. The GPU has 4352 CUDA cores, 11GB device memory. The algorithm has been written in python with packages Numpy and Scipy. The CPU implementation is based on PyFFTW which is a python wrapper of the C library FFTW. The GPU implementation is based on CuPy which is an implementation of NumPy-compatible multi-dimensional array on CUDA. As shown in Table 3.1, GPU version is up to $15$ times faster than CPU version depending on different $N$ s.