Translational and rotational non-Gaussianities in homogeneous freely evolving granular gases

TL;DR

This paper investigates the non-Gaussian velocity distributions in homogeneous, freely evolving granular gases with rough particles, using a Sonine approximation and validating results through simulations.

Contribution

It introduces a theoretical framework for non-Gaussianities in rough granular gases considering translational and rotational degrees of freedom, improving upon previous Maxwellian assumptions.

Findings

Non-Gaussian velocity distributions characterized by cumulants.

High-velocity tails deviate from Maxwellian predictions.

Validation of theoretical results with simulations.

Abstract

The importance of roughness in the modeling of granular gases has been increasingly considered in recent years. In this paper, a freely evolving homogeneous granular gas of inelastic and rough hard disks or spheres is studied under the assumptions of the Boltzmann kinetic equation. The homogeneous base state reached by the system is studied from a theoretical point of view using a Sonine approximation, in contrast to a previous Maxwellian approach. A general theoretical description is done in terms of $d_t$ translational and $d_r$ rotational degrees of freedom, which accounts for the cases of spheres ($d_t=d_r=3$) and disks ($d_t=2$, $d_r=1$) within a unified framework. The non-Gaussianities of the velocity distribution function of this state are determined by means of the first nontrivial cumulants and by the derivation of non-Maxwellian high-velocity tails. The results are validated by computer simulations using direct simulation Monte Carlo and event-driven molecular dynamics algorithms.

Figures

Tables

Table 1. Table 1: Relevant collisional moments from the Sonine approximation in the hard-disk case.

$(p,q)$		${\mu }_{pq}/\sqrt{2\pi }$
$(2,0)$		$\left[2\overline{\alpha }\left(1-\overline{\alpha }\right)+\overline{\beta }\left(1-\overline{\beta }\right)\right]\left(1+\frac{3}{16}{a}_{20}\right)-\theta \frac{{\overline{\beta }}^2}{\kappa }\left(1-\frac{a_{20}}{16}+\frac{a_{11}}{4}\right)$
$(0,2)$		$\frac{\overline{\beta }}{\kappa }\left(1-\frac{\overline{\beta }}{\kappa }\right)\left(1-\frac{a_{20}}{16}+\frac{a_{11}}{4}\right)-\frac{{\overline{\beta }}^2}{\kappa \theta }\left(1+\frac{3}{16}{a}_{20}\right)$
$(4,0)$		$8{\overline{\alpha }}^3\left(2-\overline{\alpha }\right)\left(1+\frac{15}{16}{a}_{20}\right)+3{\overline{\beta }}^3\left(2-\overline{\beta }\right)\left(1+\frac{15}{16}{a}_{20}\right)+\left(\overline{\alpha }+\frac{\overline{\beta }}{2}\right)\left(9+\frac{223}{16}{a}_{20}\right)-{\overline{\alpha }}^2\left(17+\frac{327}{16}{a}_{20}\right)$
		$-{\overline{\beta }}^2(15+\frac{281}{16}{a}_{20})-4\overline{\alpha }\overline{\beta }(\overline{\alpha }\overline{\beta }-\overline{\alpha }-\overline{\beta })(1+\frac{15}{16}{a}_{20})-4\overline{\alpha }\overline{\beta }(1+\frac{23}{16}{a}_{20})-\frac{{\overline{\beta }}^2\theta }{2\kappa }\{9+\frac{35}{16}{a}_{20}+\frac{27}{4}{a}_{11}$
		$-4[2\overline{\alpha }(1-\overline{\alpha })+3\overline{\beta }(1-\overline{\beta })](1+\frac{3}{16}{a}_{20}+\frac{3}{4}{a}_{11})+6\frac{{\overline{\beta }}^2\theta }{\kappa }(1-\frac{1}{16}{a}_{20}+\frac{1}{2}{a}_{02}+\frac{1}{2}{a}_{11})\}$
$(0,4)$		$\frac{\overline{\beta }}{\kappa }[3(1-\frac{1}{16}{a}_{20}+a_{02}+\frac{1}{2}{a}_{11})-3\frac{\overline{\beta }}{\kappa }(2-2\frac{\overline{\beta }}{\kappa }+\frac{{\overline{\beta }}^2}{{\kappa }^2})(1-\frac{1}{16}{a}_{20}+\frac{1}{2}{a}_{02}+\frac{1}{2}{a}_{11})-\frac{3}{2}\frac{\overline{\beta }}{\kappa }{a}_{02}$
		$-3\frac{\overline{\beta }}{\theta }(1-2\frac{\overline{\beta }}{\kappa }+2\frac{{\overline{\beta }}^2}{{\kappa }^2})(1+\frac{3}{16}{a}_{20}+\frac{3}{4}{a}_{11})-3\frac{{\overline{\beta }}^3}{\kappa {\theta }^2}(1+\frac{15}{16}{a}_{20})]$
$(2,2)$		$\left[\overline{\alpha }\left(1-\overline{\alpha }\right)+\frac{\overline{\beta }}{2}\left(1-\overline{\beta }\right)\right]\left(1+\frac{3}{16}{a}_{20}+\frac{3}{4}{a}_{11}\right)+\left(\overline{\alpha }+\frac{\overline{\beta }}{2}\right)\frac{a_{11}}{2}-\frac{{\overline{\beta }}^2\theta }{\kappa }\left(1-\frac{1}{16}{a}_{20}+\frac{1}{2}{a}_{11}+\frac{3}{4}{a}_{02}\right)$
		$+\frac{\overline{\beta }}{\kappa }\left\{\frac{5}{4}+\frac{23}{64}{a}_{20}+\frac{27}{16}{a}_{11}-2\left[\overline{\alpha }\left(1-\overline{\alpha }\right)+\overline{\beta }\left(1-\frac{3}{2}\overline{\beta }\right)\right]\left(1+\frac{3}{16}{a}_{20}+\frac{3}{4}{a}_{11}\right)-\left(\overline{\alpha }+\overline{\beta }\right)a_{11}\right\}+3\frac{{\overline{\beta }}^3\theta }{{\kappa }^2}$
		$\times \left(1-\frac{a_{20}}{16}+\frac{a_{02}}{2}+\frac{a_{11}}{2}\right)-\frac{7}{4}\frac{{\overline{\beta }}^2}{\kappa \theta }\left(1+\frac{129}{112}{a}_{20}\right)-\frac{5}{4}\frac{{\overline{\beta }}^2}{{\kappa }^2}\left(1+\frac{23}{80}{a}_{20}+\frac{3}{4}{a}_{11}\right)+\frac{{\overline{\beta }}^2}{\kappa \theta }\left[2\overline{\alpha }\left(1-\overline{\alpha }\right)+3\overline{\beta }\left(1-\overline{\beta }\right)\right]$
		$\times \left(1+\frac{15}{16}{a}_{20}\right)+\frac{{\overline{\beta }}^2}{{\kappa }^2}\left[2\overline{\alpha }\left(1-\overline{\alpha }\right)+3\overline{\beta }\left(1-2\overline{\beta }\right)\right]\left(1+\frac{3}{16}{a}_{20}+\frac{3}{4}{a}_{11}\right)+3\frac{{\overline{\beta }}^4\theta }{{\kappa }^3}\left(1-\frac{a_{20}}{16}+\frac{a_{02}}{2}+\frac{a_{11}}{2}\right)$

Equations

$$(\hat{\mathbf{\sigma }}\cdot {\mathbf{g}}_{12}^{\prime })=-\alpha (\hat{\mathbf{\sigma }}\cdot {\mathbf{g}}_{12}),(\hat{\mathbf{\sigma }}\times {\mathbf{g}}_{12}^{\prime })=-\beta (\hat{\mathbf{\sigma }}\times {\mathbf{g}}_{12}),$$

$${\mathbf{c}}_{1,2}^{\prime }:={\mathcal{B}}_{12,\hat{\mathbf{\sigma }}}{\mathbf{c}}_{1,2}=$$

$${\mathbf{w}}_{1,2}^{\prime }:={\mathcal{B}}_{12,\hat{\mathbf{\sigma }}}{\mathbf{w}}_{1,2}=$$

$$\mathbf{c}=\frac{\mathbf{v}}{v_{\mathrm{th}}(t)},\mathbf{w}=\frac{\mathbf{\omega }}{{\omega }_{\mathrm{th}}(t)}.$$

$$\frac{d_t}{2}{T}_t(t)=\frac{m}{2}⟨v^2⟩,\frac{d_r}{2}{T}_r(t)=\frac{I}{2}⟨{\omega }^2⟩,$$

$$n=\int \mathrm{d}\mathbf{v}\int \mathrm{d}\mathbf{\omega }f(\mathbf{v},\mathbf{\omega };t),$$

$$T(t)=\frac{d_t{T}_t(t)+d_r{T}_r(t)}{d_t+d_r}.$$

$${\mathbf{\Delta }}_{12}=\overline{\alpha }({\mathbf{c}}_{12}\cdot \hat{\mathbf{\sigma }})\hat{\mathbf{\sigma }}+\overline{\beta }\left[{\mathbf{c}}_{12}-({\mathbf{c}}_{12}\cdot \hat{\mathbf{\sigma }})\hat{\mathbf{\sigma }}-2\sqrt{\frac{\theta }{\kappa }}\hat{\mathbf{\sigma }}\times {\mathbf{W}}_{12}\right]$$

$$\overline{\alpha }=\frac{1+\alpha }{2},\overline{\beta }=\frac{\kappa }{1+\kappa }\frac{1+\beta }{2}$$

$${\mathbf{c}}_{1,2}^{\prime \prime }:={\mathcal{B}}_{12,\hat{\mathbf{\sigma }}}^{-1}{\mathbf{c}}_{1,2}=$$

$${\mathbf{w}}_{1,2}^{\prime \prime }:={\mathcal{B}}_{12,\hat{\mathbf{\sigma }}}^{-1}{\mathbf{w}}_{1,2}=$$

$${\mathbf{\Delta }}_{12}^{-}=\overline{\alpha }\left(\frac{1}{\alpha }-\frac{1}{\beta }\right)({\mathbf{c}}_{12}\cdot \hat{\mathbf{\sigma }})\hat{\mathbf{\sigma }}+\frac{{\mathbf{\Delta }}_{12}}{\beta }.$$

$$\frac{\partial f(\mathbf{\Gamma };t)}{\partial t}={\sigma }^{d_t-1}{\mathcal{I}}_{\mathbf{\Gamma }}[f,f],$$

$${\mathcal{I}}_{{\mathbf{\Gamma }}_1}[f,f]=\int \mathrm{d}{\mathbf{\Gamma }}_2{\int }_{+}\mathrm{d}\hat{\mathbf{\sigma }}({\mathbf{v}}_{12}\cdot \hat{\mathbf{\sigma }})\left(\frac{f_1^{{}^{\prime \prime }}{f}_2^{{}^{\prime \prime }}}{\alpha J}-f_1{f}_2\right)$$

$$J=\left|\frac{\partial ({\mathbf{v}}_1^{\prime },{\mathbf{v}}_2^{\prime },{\mathbf{\omega }}_1^{\prime },{\mathbf{\omega }}_2^{\prime })}{\partial ({\mathbf{v}}_1,{\mathbf{v}}_2,{\mathbf{\omega }}_1,{\mathbf{\omega }}_2)}\right|=\left|\frac{\partial ({\mathbf{v}}_1,{\mathbf{v}}_2,{\mathbf{\omega }}_1,{\mathbf{\omega }}_2)}{\partial ({\mathbf{v}}_1^{\prime \prime },{\mathbf{v}}_2^{\prime \prime },{\mathbf{\omega }}_1^{\prime \prime },{\mathbf{\omega }}_2^{\prime \prime })}\right|=\alpha |\beta {|}^{2d_r/d_t}.$$

$$s(t)=\frac{1}{2}{\int }_0^t\mathrm{d}{t}^{\prime }\nu (t^{\prime }),$$

$$\nu (t)=Kn{\sigma }^{d_t-1}{v}_{\mathrm{th}}(t),K\equiv \frac{\sqrt{2}{\pi }^{\frac{d_t-1}{2}}}{\Gamma \left(\frac{d_t}{2}\right)}.$$

$$\varphi (\tilde{\mathbf{\Gamma }};s)=n^{-1}{v}_{\mathrm{th}}^{d_t}(t){\omega }_{\mathrm{th}}^{d_r}(t)f(\mathbf{\Gamma };t).$$

$$\frac{K}{2}{\partial }_s\varphi +\frac{{\mu }_{20}^{(0)}}{d_t}\frac{\partial }{\partial \mathbf{c}}\cdot \left(\mathbf{c}\varphi \right)+\frac{{\mu }_{02}^{(0)}}{d_r}\frac{\partial }{\partial \mathbf{w}}\cdot \left(\mathbf{w}\varphi \right)={\mathcal{I}}_{\tilde{\mathbf{\Gamma }}}[\varphi ,\varphi ],$$

$${\mu }_{pq}^{(r)}=$$

$$=$$

$$\times ({\mathcal{B}}_{12,\hat{\mathbf{\sigma }}}-1)\left[c_1^p{w}_1^q({\mathbf{c}}_1\cdot {\mathbf{w}}_1{)}^r+c_2^p{w}_2^q({\mathbf{c}}_2\cdot {\mathbf{w}}_2{)}^r\right]$$

$$\frac{K}{2}{\partial }_s{T}_t=-\frac{2}{d_t}{\mu }_{20}^{(0)}{T}_t,\frac{K}{2}{\partial }_s{T}_r=-\frac{2}{d_r}{\mu }_{02}^{(0)}{T}_r,$$

$$\frac{K}{2}{\partial }_s\ln \theta =2\left[\frac{{\mu }_{20}^{(0)}}{d_t}-\frac{{\mu }_{02}^{(0)}}{d_r}\right],$$

$$\frac{K}{2}{\partial }_{s}T=-{\zeta }^{\ast }T,$$

$$\frac{K}{2}\frac{\partial \ln M_{pq}^{(r)}}{\partial s}-\frac{(p+r){\mu }_{20}^{(0)}}{d_t}-\frac{(q+r){\mu }_{02}^{(0)}}{d_r}=-\frac{{\mu }_{pq}^{(r)}}{M_{pq}^{(r)}}.$$

$$\delta {\mathcal{B}}_{12,\hat{\mathbf{\sigma }}}\left(c_1^2+c_2^2\right)=$$

$$\times {\left(\hat{\mathbf{\sigma }}\times {\mathbf{c}}_{12}\right)}^2+8{\overline{\beta }}^2\frac{\theta }{\kappa }(\hat{\mathbf{\sigma }}\times {\mathbf{W}}_{12}{)}^2$$

$$-4\overline{\beta }(2\overline{\beta }-1)\sqrt{\frac{\theta }{\kappa }}{\mathbf{c}}_{12}\cdot (\hat{\mathbf{\sigma }}\times {\mathbf{W}}_{12}),$$

$$\delta {\mathcal{B}}_{12,\hat{\mathbf{\sigma }}}\left(w_1^2+w_2^2\right)=$$

$$\times {\left(\hat{\mathbf{\sigma }}\times {\mathbf{W}}_{12}\right)}^2+4\frac{\overline{\beta }}{\sqrt{\kappa \theta }}\left(2\frac{\overline{\beta }}{\kappa }-1\right)$$

$$\times {\mathbf{W}}_{12}\cdot (\hat{\mathbf{\sigma }}\times {\mathbf{c}}_{12}).$$

$$⟨⟨\psi ⟩⟩=\int \mathrm{d}{\tilde{\mathbf{\Gamma }}}_1\int \mathrm{d}{\tilde{\mathbf{\Gamma }}}_2\psi ({\tilde{\mathbf{\Gamma }}}_1,{\tilde{\mathbf{\Gamma }}}_2)\varphi ({\tilde{\mathbf{\Gamma }}}_1)\varphi ({\tilde{\mathbf{\Gamma }}}_2).$$

$${\mu }_{20}^{(0)}=$$

$$\displaystyle-2\overline{\beta}^{2}\frac{\theta}{\kappa}\left[3\langle\!\langle c_{12}W_{12}^{2}\rangle\!\rangle-\langle\!\langle c_{12}^{-1}(\mathbf{c}_{12}\cdot\mathbf{W}_{12})^{2}\rangle\!\rangle\right]\biggr{\}},$$

$${\mu }_{02}^{(0)}=$$

Full text

I Introduction

Granular systems are themselves worth studying from mechanical, physical, and mathematical points of view. They are very commonly observed in nature, where different geometries can take place. Grains, from a dynamical point of view, move in a three-dimensional space, but constraints make two-dimensional problems become real and of special interest [1, 2, 3, 4, 5, 6, 7, 8, 9, 10].

We will focus on the description of granular systems at low-density fluidized states, where the assumptions underlying the Boltzmann equation apply [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28]. The simplest collisional model for interactions in granular gaseous flows is the inelastic hard-sphere model, where the granular gas is assumed to be composed by inelastic and smooth identical hard disks, spheres, or hyperspheres in ${d_{t}}$ translational space dimensions [29, 30, 31, 32, 28]. However, this description might be limiting and can be improved by considering rotational degrees of freedom, which may play an important role in the dynamics of granular gases by means of surface roughness. Here, we will use the simplest collisional model that implements roughening, the inelastic and rough hard-sphere model. In the latter model, the inelasticity is parameterized by a constant coefficient of normal restitution, $\alpha$ (in common with the inelastic hard-sphere model), and roughness is introduced by means of a coefficient of tangential restitution, $\beta$. Although, in general, the effective coefficient of tangential restitution depends on the impact angle because of friction [33, 34], here we adopt the simplest model with constant $\beta$, as frequently done in the literature [12, 35, 32, 36, 37, 38, 39, 40, 41, 42, 28, 43, 44].

Whereas in the inelastic hard-sphere model, the ${d_{t}}$-dimensional description of an inelastic gas of smooth and spinless hard (hyper)spheres is straightforward, in the case of rough spheres, where angular velocities come into play, rotational degrees of freedom, ${d_{r}}$, need to be introduced. A description in terms of ${d_{t}}$ and ${d_{r}}$ becomes highly dependent on the geometry and constraints of the system as antecedently reported [43, 44, 45, 46]. As in Ref. [45], we will derive the general description to be valid just for the two relevant cases of hard disks and hard spheres. For hard disks, angular velocities are constrained to the direction orthogonal to the plane of motion, so ${d_{t}}=2$ and ${d_{r}}=1$. For hard spheres, on the other hand, angular and translational velocities are vectors of a Euclidean three-dimensional vector space, i.e., ${d_{t}}={d_{r}}=3$.

It is widely known and studied that the homogeneous Boltzmann equation—for both smooth and rough models—admits a scaling solution, in which the system cools down continuously and the whole evolution is driven by the granular temperature. This state is known as the homogeneous cooling state (HCS), and has been of interest for the granular gas community in the last three decades [32, 28, 47, 48, 49, 50, 51, 52, 53, 54, 55, 42, 43, 56, 39]. It is worth mentioning that, very recently, the HCS has been experimentally observed in microgravity experiments by Yu et al. [57]. In that paper, both Haff’s cooling law and the exponential high-velocity tail of the velocity distribution function (VDF) as predicted by kinetic theory together with the inelastic hard-sphere model are confirmed. In Ref. [57], while the results were compared with the constant and velocity-dependent models for the coefficient of normal restitution, it was concluded that the latter model had a negligible influence on the results, supporting the approximation of constant coefficients of restitution. In fact, the system in Ref. [57] is compatible with a constant $\alpha=0.66$, highlighting that the latter approximation is not subjected only to quasielastic systems. Moreover, the authors proposed a possible influence of surface roughness in the collisional rules due to an overestimate of the relaxation time from the inelastic hard-sphere model as compared with the experimental outcomes. Thus, they claimed that the rotational degrees of freedom could be an answer to these deviations.

Theoretically, some of the early attempts to study the Gaussian deviations of the HCS VDF of a granular gas of rough particles were done in Refs. [12, 35] using a Sonine expansion, that is, an isotropic expansion around a two-temperature—translational and rotational—Maxwellian VDF. However, although velocity correlations were not originally assumed for hard spheres, they were proved to be present [36]. More recently, the first nontrivial velocity cumulants were studied for freely evolving hard spheres [39]. Throughout the present paper, we will expose the results of those cumulants in a common frame for both disks and spheres from the collisional-moment point of view, and we will analyze results for hard disks.

In the case of freely cooling inelastic granular gases, deviations of the HCS VDF from a Maxwellian derived from the inelastic hard-sphere model are not only accounted for by the first nontrivial cumulants, but also an exponential high-velocity tail for this distribution was predicted by kinetic theory and computer simulations [58, 48, 49, 51, 59] and satisfactorily observed experimentally not only for freely evolving granular gases [57] as commented above, but also for uniformly heated systems [60]. Results in Ref. [39] put into manifest a highly populated tail for the marginal VDF of angular velocities in the inelastic and rough hard-sphere model, accompanied by high values of the fourth angular velocity cumulant for some values of the pair $(\alpha,\beta)$. However, this marginal distribution was interpreted as being consistent with an exponential form [39], similarly to what occurs in the inelastic hard-sphere model with the total (translational) VDF [58, 48, 49]. In this paper, we study the high-velocity tail for the marginal VDF of the translational and angular velocities, and for their product as well, where theory indicates algebraic tails for the two latter marginal distributions.

The study of the non-Gaussianities of the HCS VDF is also motivated by recent research in nonhomogeneous states [61, 45, 46] from Chapman–Enskog expansions around the HCS solution. Linear stability analyses of the homogeneous state hydrodynamics show that the Maxwellian approximation of the first-order VDF might not work for the hard-disk case [46], and for some very small region of the parameter space for hard spheres, yielding a wrong prediction of a completely unstable region in the parameter space. Therefore, it was conjectured that non-Gaussianities might be crucial in the cited region of parameters, this effect being more important for disks than for spheres [46]. This hypothesis was supported by high values of the first relevant cumulants in hard-sphere systems [39] and by results from the smooth case, where the homogeneous VDF is generally more disparate from a Maxwellian for disks than for spheres [62].

The present paper is structured as follows. In Sec. II, the inelastic and rough hard-sphere model and the binary collisional rules are introduced. Afterwards, the framework of the homogeneous Boltzmann equation is described in Sec. III and the hierarchy of evolution equations, the Sonine expansion of the VDF, and the definitions of cumulants are formally presented. Section IV is devoted to the Sonine approximation, where the infinite expansion is truncated beyond the first few nontrivial coefficients, the associated collisional moments are explicitly written for hard disks, and the HCS cumulants are obtained. Next, we study the forms of the marginal VDF from the Maxwellian and Sonine approximations, as well as their high-velocity tails in Sec. V in the context of the Boltzmann equation. In Sec. VI, we compare the theoretical predictions with direct simulation Monte Carlo (DSMC) and event-driven molecular dynamics (EDMD) computer simulations outcomes. Finally, concluding remarks and main results are summarized in Sec. VII.

II Inelastic and Rough Hard Particles

II.1 System

Let us consider a monodisperse dilute granular gas of hard disks or spheres, which are assumed to be inelastic and rough, their dynamics being described by their translational and angular velocities, $\mathbf{v}$ and $\bm{\omega}$, respectively. Whereas for spheres (${d_{t}}={d_{r}}=3$), both $\mathbf{v}$ and $\bm{\omega}$ are vectors in an Euclidean three-dimensional space, this is not the case for disks (${d_{t}}=2$, ${d_{r}}=1$), where $\bm{\omega}$ is a one-dimensional vector orthogonal to the two-dimensional vector space spanned by $\mathbf{v}$. In general, however, all vector relations will be written in a three-dimensional Euclidean embedding space.

The gas is considered to be formed by a fixed large number of identical hard ${d_{t}}$ spheres with mass $m$, diameter $\sigma$, reduced moment of inertia $\kappa=4I/m\sigma^{2}$ ($I$ being the moment of inertia), and whose inelasticity and roughness are characterized by a coefficient of normal restitution, $\alpha$, and a coefficient of tangential restitution, $\beta$, respectively, both assumed to be constant, and defined by

[TABLE]

where $\widehat{\bm{\sigma}}$ is the intercenter unit vector at contact, $\mathbf{g}_{12}=\mathbf{v}_{12}-\frac{\sigma}{2}\widehat{\bm{\sigma}}\times(\bm{\omega}_{1}+\bm{\omega}_{2})$ is the relative velocity of the contact points of particles 1 and 2 (with $\mathbf{v}_{12}=\mathbf{v}_{1}-\mathbf{v}_{2}$), and primed quantities refer to postcollisional values. Note that, while $\alpha$ is nonnegative, $\beta$ can be either positive or negative. A negative value means that the postcollisional tangential component of the relative velocity maintains the same sign as the precollisional one, implying that the effect of surface friction is not dramatic. On the other hand, if the particles are sufficiently rough, the sign of the tangential component is inverted upon collision. Figure 1 presents a sketch illustrating a collision between (a) two hard disks and (b) two hard spheres.

II.2 Direct collisional rules

The direct binary collisional rules are obtained from the assumption of conservation of linear and angular momenta at the point of contact in each collision. They can be expressed by [32, 39, 63, 42, 28, 43, 44, 45, 46]

[TABLE]

where $\mathcal{B}_{12,\widehat{\bm{\sigma}}}$ is the postcollisional operator acting on a dynamic quantity and giving the result after a collision, $\theta$ and $\bm{\Delta}_{12}$ are defined below, and $\{\mathbf{c},\mathbf{w}\}$ are the velocities reduced by their thermal value, that is,

[TABLE]

Here, $v_{\mathrm{th}}(t)=\sqrt{{2T_{t}(t)}/{m}}$, $\omega_{\mathrm{th}}(t)=\sqrt{{2T_{r}(t)}{/I}}$ are the thermal translational and angular velocities, $T_{t}$ and $T_{r}$ being the translational and rotational granular temperatures, respectively, which are defined by [32, 39, 63, 42, 28, 43, 44, 45, 46]

[TABLE]

where $\langle\cdots\rangle=n^{-1}\int\mathrm{d}\mathbf{v}\int\mathrm{d}\bm{\omega}\medspace(\cdots)f(\mathbf{v},\bm{\omega};t)$ represents a one-body average value with respect to the VDF $f(\mathbf{v},\bm{\omega};t)$ normalized as

[TABLE]

$n$ being the particle number density. In Eq. (2b), $\theta\equiv T_{r}/T_{t}$ is the rotational-to-translational granular temperature ratio. Moreover, the mean granular temperature is

[TABLE]

Finally, the quantity

[TABLE]

is the reduced impulse. In Eq. (7), $\mathbf{W}_{12}\equiv\frac{1}{2}(\mathbf{w}_{1}+\mathbf{w}_{2})$ and

[TABLE]

are effective coefficients of restitution.

Notice that the sets of vectors $\{\mathbf{c},\mathbf{w}\}$ and $\{\mathbf{v},\bm{\omega}\}$ span the same vector spaces, respectively. Therefore, the use of reduced quantities will be algebraically equivalent to the original velocity description.

II.3 Inverse collisional rules

The inverse collisional rules relating precollisional velocities $\{\mathbf{c}_{1}^{\prime\prime},\mathbf{w}_{1}^{\prime\prime},\mathbf{c}_{2}^{\prime\prime},\mathbf{w}_{2}^{\prime\prime}\}$ to postcollisional velocities $\{\mathbf{c}_{1},\mathbf{w}_{1}\,\mathbf{c}_{2},\mathbf{w}_{2}\}$ are [32, 39, 63, 42, 28, 43, 44, 45, 46]

[TABLE]

with

[TABLE]

From now on, throughout this paper, we will adopt the notation ${\mathbf{\Gamma}}=\{\mathbf{v},\bm{\omega}\}$ and ${\widetilde{\mathbf{\Gamma}}}=\{\mathbf{c},\mathbf{w}\}$.

III Boltzmann equation

III.1 Basics

We will carry out a description of the system under the assumption of molecular chaos or Stosszahlansatz [64], basing the analytical treatment on the homogeneous Boltzmann equation. As said before, we will generally derive the results keeping a dependence on the number of degrees of freedom, ${d_{t}}$ and ${d_{r}}$. The homogeneous Boltzmann equation reads

[TABLE]

where

[TABLE]

is the collision operator. Here, $f_{1,2}\equiv f({\mathbf{\Gamma}}_{1,2})$, $f_{1,2}^{\prime\prime}\equiv f({\mathbf{\Gamma}}_{1,2}^{\prime\prime})$, the subscript $+$ designates the constraint $\mathbf{v}_{12}\cdot\widehat{\bm{\sigma}}>0$, and $J$ is the Jacobian due to the collisional change of velocities [43], i.e.,

[TABLE]

Since the temporal change of the VDF is subjected only to collisions, it is convenient to change from laboratory time, $t$, to collisional time, $s$, as given by

[TABLE]

where $\nu(t)$ is the (nominal) collision frequency, defined by

[TABLE]

This variable $s(t)$ quantifies the accumulated average number of collisions per particle up to time $t$. Furthermore, the treatment based on the reduced velocities, ${\widetilde{\mathbf{\Gamma}}}$, allows us to define the reduced one-body VDF:

[TABLE]

The homogeneous Boltzmann equation for the reduced VDF then reads

[TABLE]

where

[TABLE]

are (reduced) collisional moments. Note that, in the particular case of disks on a plane, the index $r$ is meaningless due to the orthogonality between the vector spaces spanned by translational and angular velocities [see Fig. 1(a)]. However, from a general point of view, the three-dimensional vector forms will be maintained.

Upon derivation of Eq. (17), use has been made of the evolution equations for the translational and rotational temperatures,

[TABLE]

which imply

[TABLE]

where $\zeta^{*}\equiv 2(\mu_{20}^{(0)}+\mu_{02}^{(0)}\theta)/({d_{t}}+{d_{r}}\theta)$ is the (reduced) cooling rate, and thus Eq. (20b) represents Haff’s cooling law [65] for the inelastic and rough hard-sphere model.

From Eq. (17), one can directly derive the hierarchy equations for the evolution of the velocity moments $M_{pq}^{(r)}\equiv\langle c^{p}w^{q}(\mathbf{c}\cdot\mathbf{w})^{r}\rangle$:

[TABLE]

III.2 Collisional moments

The collisional change of a certain velocity function can be obtained by application of the operator $\delta\mathcal{B}_{12,\widehat{\bm{\sigma}}}\equiv\mathcal{B}_{12,\widehat{\bm{\sigma}}}-1$ on the function. For instance,

[TABLE]

The results for $\delta\mathcal{B}_{12,\widehat{\bm{\sigma}}}\left(c_{1}^{4}+c_{2}^{4}\right)$, $\delta\mathcal{B}_{12,\widehat{\bm{\sigma}}}\left(w_{1}^{4}+w_{2}^{4}\right)$, $\delta\mathcal{B}_{12,\widehat{\bm{\sigma}}}\left(c_{1}^{2}w_{1}^{2}+c_{2}^{2}w_{2}^{2}\right)$, and $\delta\mathcal{B}_{12,\widehat{\bm{\sigma}}}\left[(\mathbf{c}_{1}\cdot\mathbf{w}_{1})^{2}+(\mathbf{c}_{2}\cdot\mathbf{w}_{2})^{2}\right]$ can be found in the Supplemental Material [66].

Inserting the collisional changes into Eq. (III.1), the collisional moments $\mu_{pq}^{(r)}$ can be formally expressed in terms of two-body averages of the form

[TABLE]

In particular,

[TABLE]

where the factor $B_{3}=\pi^{\frac{{d_{t}}-1}{2}}/\Gamma\left(\frac{{d_{t}}+3}{2}\right)$ comes from an angular integral. The formally exact expressions of the collisional moments $\mu_{40}^{(0)}$, $\mu_{04}^{(0)}$, and $\mu_{22}^{(0)}$ in terms of two-body averages are given in the Supplemental Material [66], where also some related tests for the simulation data are included.

III.3 Sonine expansion

Assuming isotropy, $\phi({\widetilde{\mathbf{\Gamma}}};s)$ must depend on velocity only through three scalars: $c^{2}$, $w^{2}$, and $(\mathbf{c}\cdot\mathbf{w})^{2}$. This can be made explicit by the polynomial expansion

[TABLE]

where

[TABLE]

is the (two-temperature) Maxwellian distribution, $a_{jk}^{(\ell)}$ are Sonine coefficients, and the functions

[TABLE]

form a complete set of orthogonal polynomials [39]. Here, $L_{j}^{(\ell)}(x)$ are associated Laguerre polynomials, $u\equiv(\mathbf{c}\cdot\mathbf{w})/cw$ is the cosine of the angle formed by the vectors $\mathbf{c}$ and $\mathbf{w}$, and $P_{\ell}(u)$ are Legendre polynomials [67]. The orthogonality condition is

[TABLE]

where the inner product of two arbitrary real functions $\Phi_{1}({\widetilde{\mathbf{\Gamma}}})$ and $\Phi_{2}({\widetilde{\mathbf{\Gamma}}})$ is defined as

[TABLE]

Note that $\langle\Phi\rangle=\innerproduct{\Phi}{\phi/\phi_{\text{M}}}$. Using Eq. (28a) in Eq. (25), one can express the Sonine coefficients as

[TABLE]

In particular, $a_{00}^{(0)}=1$, $a_{10}^{(0)}=a_{01}^{(0)}=0$, while

[TABLE]

are fourth-order cumulants. Notice that ${a_{00}^{(1)}}$ is only meaningful in the hard-sphere case and thus it is not expressed in terms of the number of degrees of freedom.

The evolution equations for the cumulants defined by Eqs. (31) can be easily obtained from the moment hierarchy, Eq. (21), as

[TABLE]

where henceforth we simplify the notation as ${a_{jk}^{(0)}}\to a_{jk}$ and $\mu_{jk}^{(0)}\to\mu_{jk}$.

III.4 Homogeneous cooling state

The scaling method in the description of the kinetic equation suggests that a stationary solution, $\phi=\phi^{\mathrm{H}}$, of Eq. (17) applies for long times (hydrodynamic limit). This is the HCS, in which the temperature ratio, $\theta^{\mathrm{H}}$, is constant and the whole time dependence of the unscaled VDF $f^{\mathrm{H}}({\mathbf{\Gamma}};t)$ takes place through the mean temperature $T(t)$ only. On the other hand, this stationary solution $\phi=\phi^{\mathrm{H}}$ is not exactly known.

From Eqs. (20a) and (32), it follows that, in the HCS,

[TABLE]

Notice that, as expected, Eq. (33c) is only meaningful for spheres (${d_{t}}={d_{r}}=3$).

IV Approximate schemes

All the equations presented in Sec. III are formally exact within the framework of the Boltzmann equation. However, no explicit results can be obtained unless one makes use of approximations.

IV.1 Maxwellian approximation

The simplest approximation is the Maxwellian one, i.e., $\phi({\widetilde{\mathbf{\Gamma}}})\to\phi_{\text{M}}({\widetilde{\mathbf{\Gamma}}})$. In that case [45],

[TABLE]

In this Maxwellian approximation, Eqs. (19) and (20) can be solved to get the evolution of the partial and mean temperatures, as well as the HCS value of the temperature ratio $\theta^{\mathrm{H}}$. However, by construction, the Maxwellian approximation is unable to account for the non-Gaussianities of the VDF, either in the transient evolution to the HCS or in the HCS itself.

IV.2 Sonine approximation

The basic quantities measuring non-Gaussianities are the cumulants defined in Eqs. (31). Therefore, as the simplest scheme to capture those cumulants, we introduce the Grad–Sonine methodology [68, 39, 28] and truncate the infinite Sonine expansion, Eq. (25), after $j+k+2\ell\geq 3$, i.e.,

[TABLE]

where the term $a_{00}^{(1)}\Psi_{00}^{(1)}$ is not present in the hard-disk case. With the replacement given by Eq. (IV.2), the two-body averages appearing in the collisional moments [see, for instance, Eqs. (24)] can be explicitly calculated as linear and quadratic functions of the cumulants. Next, our Sonine approximation is constructed by neglecting quadratic terms, so only linear terms are retained.

By particularizing to the hard-sphere case (${d_{t}}={d_{r}}=3$), previous results are recovered [39]. Moreover, we obtain novel expressions for hard disks (${d_{t}}=2$, ${d_{r}}=1$), which are displayed in Table 1. Further details about some of the computations are available in the Supplemental Material [66].

For consistency with the truncation and linearization steps carried out in the Sonine approximation, the evolution equations in the hard-disk case are obtained by inserting the expressions in Table 1 into Eqs. (20a) and (32a)–(32c), and linearizing the bracketed quantities. This gives a closed set of four differential equations, which are linear in the cumulants and nonlinear in the temperature ratio. Likewise, the HCS values are obtained by linearizing Eqs. (33a) and (33b) with respect to the cumulants. The linear stability of the HCS versus uniform and isotropic perturbations is proved in Appendix A.

Figures 2 and 3 show the HCS quantities $\theta^{\mathrm{H}}$, $a_{20}^{\mathrm{H}}$, $a_{02}^{\mathrm{H}}$, and $a_{11}^{\mathrm{H}}$, obtained from the Sonine approximation for uniform disks ($\kappa=\frac{1}{2}$) and spheres ($\kappa=\frac{2}{5}$), respectively, as functions of the coefficients of restitution $\alpha$ and $\beta$. In the hard-sphere case, the cumulant $a_{00}^{(1)\mathrm{H}}$ is also included. It can be observed that, typically, hard-disk systems depart from the Maxwellian state more than hard-sphere systems. Interestingly, both hard-disk and hard-sphere systems present relatively large values of $a_{02}^{\mathrm{H}}$ and $a_{11}^{\mathrm{H}}$, thus signaling a possible quantitative breakdown of the Sonine approximation, which implicitly assumes small deviations from the Maxwellian VDF.

V Marginal distribution functions and High-velocity tails in the Homogeneous Cooling State

V.1 Marginal distribution functions

As said before, the reduced VDF $\phi({\widetilde{\mathbf{\Gamma}}})$ in isotropic states depend on the three scalars $c^{2}$, $w^{2}$, and $c^{2}w^{2}$ [plus $(\mathbf{c}\cdot\mathbf{w})^{2}$ only for spheres]. To disentangle those dependencies, it is convenient to define the following marginal distributions [39, 41]:

[TABLE]

where $x$ represents the product $c^{2}w^{2}$. Note that, by isotropy, $\phi_{\mathbf{c}}(\mathbf{c})$ and $\phi_{\mathbf{w}}(\mathbf{w})$ depend only on the moduli $c$ and $w$, respectively. Moreover, the marginal distributions in Eqs. (36) are directly related to the cumulants $a_{20}^{(0)}$, $a_{02}^{(0)}$, and $a_{11}^{(0)}$ defined by Eqs. (31), namely

[TABLE]

The Maxwellian expressions for these functions are

[TABLE]

where $\Omega_{d}={2\pi^{d/2}}/{\Gamma\left(\frac{d}{2}\right)}$ is the $d$-dimensional solid angle and $K_{a}(x)$ is the modified Bessel function of the second kind. In the Sonine approximation defined by Eq. (IV.2), one has

[TABLE]

While Eqs. (39) may reproduce the correct behavior of the HCS in the thermal domain, it is known from the smooth case [58, 48, 49] and from hard-sphere results [39] that they are unable to account for the high-velocity tail.

V.2 High-velocity tails

Let us now study the high-velocity tail for the marginal VDF in the HCS, in analogy to previous works for the smooth case [58, 48].

To carry out this asymptotic analysis, we start from the homogeneous Boltzmann equation, Eq. (17), and split the collisional operator into a loss and a gain term, that is [58, 48],

[TABLE]

where the loss term can be written as

[TABLE]

with $B_{1}=\pi^{\frac{{d_{t}}-1}{2}}/\Gamma\left(\frac{{d_{t}}+1}{2}\right)$. The gain term accounts for all the particles that after a collision have velocities ${\widetilde{\mathbf{\Gamma}}}_{1}$. In contrast, the loss term takes into account the amount of particles with ${\widetilde{\mathbf{\Gamma}}}_{1}$ that, after a collision, are not contributing any more to these velocities.

Intuitively, one would expect that escaping from the rapid regime is easier than entering the high-velocity limit, given the low likelihood of encountering rapid particles compared to thermal ones. Thus, the main assumption we will use is that, for high velocities of the HCS, the loss term prevails over the gain term. From Eq. (12), and following the case of smooth particles [48], the assumption above can be expressed as

[TABLE]

V.2.1 Tail of $\phi_{\mathbf{c}}^{\mathrm{H}}(\mathbf{c})$

Integrating over $\mathbf{w}$ on both sides of the stationary version of Eq. (17), neglecting the gain term, replacing $c_{12}\to c_{1}$ in Eq. (41), and taking the limit $c\gg 1$, we get the linear differential equation

[TABLE]

whose solution is

[TABLE]

where $\mathcal{A}_{c}$ is an integration constant. This is equivalent to the result in the smooth case [58, 48], except that now $\mu_{20}^{\mathrm{H}}$ takes into account the influence of surface roughness. The exponential decay of $\phi_{\mathbf{c}}^{\mathrm{H}}(\mathbf{c})$ implies that all the cumulants of the form $a_{j0}^{\mathrm{H}}$ are finite.

V.2.2 Tail of $\phi_{\mathbf{w}}^{\mathrm{H}}(\mathbf{w})$

Now we integrate over $\mathbf{c}$ on both sides of Eq. (17) and neglect again the gain term. This yields

[TABLE]

where

[TABLE]

Here, $\phi_{\mathbf{c}|\mathbf{w}}(\mathbf{c}|\mathbf{w})$ is a conditional probability distribution function defined as

[TABLE]

The quantity $\overline{c_{12}}(\mathbf{w})$ represents the average relative translational velocity of those particles with an angular velocity $\mathbf{w}$. It is a functional of the whole VDF $\phi({\widetilde{\mathbf{\Gamma}}})$, so Eq. (45) is not a closed equation for the marginal distribution $\phi_{\mathbf{w}}^{\mathrm{H}}(\mathbf{w})$.

The positive values observed in Figs. 2 and 3 for the cumulant $a_{11}^{\mathrm{H}}$ imply that high angular velocities are positively correlated to high translational velocities, so $\overline{c_{12}}^{\mathrm{H}}(\mathbf{w})$ is expected to increase with $w$. However, to estimate the tail of $\phi_{\mathbf{w}}^{\mathrm{H}}(\mathbf{w})$, we further assume that the dependence of $\overline{c_{12}}^{\mathrm{H}}(\mathbf{w})$ on $w$ is weak enough as to take $\overline{c_{12}}^{\mathrm{H}}(\mathbf{w})\approx\langle\!\langle c_{12}\rangle\!\rangle^{\mathrm{H}}$. With this adiabaticlike approximation, Eq. (45) becomes a closed linear equation whose solution is

[TABLE]

where $\mathcal{A}_{w}$ is the associated integration constant. In the expression of $\gamma_{w}$, we have made use of the HCS condition $\mu_{02}^{\mathrm{H}}/{d_{r}}=\mu_{20}^{\mathrm{H}}/{d_{t}}$ [see Eqs. (33a)].

While, in principle, Eqs. (48) are approximate because of the ansatz $\overline{c_{12}}^{\mathrm{H}}(\mathbf{w})\approx\langle\!\langle c_{12}\rangle\!\rangle^{\mathrm{H}}$, it accounts for an algebraic decay of $\phi_{\mathbf{w}}^{\mathrm{H}}(\mathbf{w})$ explaining the relatively high values attained by $a_{02}^{\mathrm{H}}$. In fact, Eqs. (48) imply that the coefficients of the form $a_{0k}^{\mathrm{H}}$ diverge if $2k\geq\gamma_{w}-1$.

V.2.3 Tail of $\phi_{cw}(x)$

Whereas the derivation of the high-velocity tail for $\phi_{\mathbf{c}}$ is clean, and the one for $\phi_{\mathbf{w}}$, although approximate, is reasonable, in the case of the distribution $\phi_{cw}$ the reasoning is somewhat more speculative. Let us start by introducing the marginal probability distribution function of the variable $w^{2}$, $\phi_{w^{2}}(w^{2})=(\Omega_{{d_{r}}}/2)w^{{d_{r}}-2}\phi_{\mathbf{w}}(\mathbf{w})$. According to Eqs. (48), the high-velocity tail of $\phi_{w^{2}}^{\mathrm{H}}(w^{2})$ is

[TABLE]

As can be inferred from Eqs. (44) and (48), the tail of angular velocities is much more populated than that of translational velocities. Therefore, it is reasonable to expect that the main contribution to $\phi_{cw}(c^{2}w^{2})$ comes essentially from particles with thermal translational velocities ($c\sim 1$) and high angular velocities ($w\gg 1$). Thus, in view of Eq. (49), we conjecture that

[TABLE]

This algebraic decay would be responsible for the relatively large values of $a_{11}^{\mathrm{H}}$ and implies the divergence of the coefficients of the form $a_{jj}^{\mathrm{H}}$ if $j\geq\gamma_{cw}-1$. An alternative justification of Eqs. (50) is provided in the Supplemental Material [66].

While, according to Eq. (44), the asymptotic decay of $\phi_{\mathbf{c}}(\mathbf{c})$ is governed by a velocity scale $c\sim\gamma_{c}^{-1}$, Eqs. (48) and (50) show that the decays of $\phi_{\mathbf{w}}(\mathbf{w})$ and $\phi_{cw}(x)$ are scale-free. It can be checked that the exponents $\gamma_{c}$, $\gamma_{w}$, and $\gamma_{cw}$ are generally smaller for disks than for spheres, meaning that the high-velocity tails are fatter in the former case than in the latter. Apart from that, they exhibit a similar qualitative dependence on the coefficients of restitution.

The consistency of Eq. (42) with the tails obtained here is discussed in Appendix B.

VI Simulation Results

To test the theoretical results, we have run two types of computer simulation algorithms for a dilute and homogeneous granular gas of inelastic and rough hard disks (${d_{t}}=2$, ${d_{r}}=1$) with different values of the coefficients of restitution $\alpha$ and $\beta$. In all cases, the disks are assumed to have a uniform mass distribution, so the reduced moment of inertia is $\kappa=\frac{1}{2}$.

First, we used DSMC, as proposed by Bird [69, 70] and conveniently adapted to the granular case [51, 39], to simulate a homogenenous and dilute granular gas of inelastic and rough hard disks, using $N=10^{4}$ representative particles. Additionally, we carried out EDMD computer simulations with $N=1600$ disks in a square box of side length $L/\sigma=565.7$, which correspond to a number density $n\sigma^{2}=0.005$, thus avoiding spatial instabilities [46]. Whereas the EDMD system has nonzero density, the solid fraction $\varphi=\frac{\pi}{4}n\sigma^{2}\simeq 3.9\times 10^{-3}$ is small enough to expect good agreement with the diluteness assumption. We ran $100$ and $50$ replicas for DSMC and EDMD, respectively, for each pair $(\alpha,\beta)$, not observing instabilities in the EDMD simulations. In addition to averaging over replicas, the stationary HCS values were measured by averaging over instantaneous values at $s=s_{\text{ini}},s_{\text{ini}}+\delta s,s_{\text{ini}}+2\delta s,\ldots,s_{\text{fin}}$ with $(s_{\text{ini}},s_{\text{fin}},\delta s)=(500,1500,5)$ and $(150,200,1)$ for DSMC and EDMD, except in the case $\alpha=0.9$, $\beta=-0.8$, where we took $(450,500,1)$ in the EDMD simulations. In the construction of histograms for the marginal distributions, we considered $2^{8}$ bins in the associated velocity variable.

In the Supplemental Material [66], we present a comparison between the Sonine-approximation results [see Eqs. (20a) and (32)] and simulation data for the temporal evolution toward the HCS of the temperature ratio and the cumulants, starting from an equipartioned Maxwellian state. A generally good agreement is observed, except for $a_{02}$ near the HCS if $a_{02}^{\mathrm{H}}$ reaches relatively high values. Now we present results for the relevant quantities in the HCS.

VI.1 Temperature ratio and cumulants

Figure 4 shows the HCS values of $\theta^{\mathrm{H}}$, $a_{20}^{\mathrm{H}}$, $a_{02}^{\mathrm{H}}$, and $a_{11}^{\mathrm{H}}$ versus $\beta$ for some representative values of $\alpha$. Figure 5 presents the same quantities versus $\alpha$ for some illustrative values of $\beta$. We observe that the Maxwellian approximation provides a good description of $\theta^{\mathrm{H}}$, although it tends to overestimate it if $\alpha\lesssim 0.7$ [see Figs. 5(a)–5(c)]. Those deviations are satisfactorily corrected by the Sonine approximation.

In the case of the cumulants, their qualitative shape as functions of both $\alpha$ and $\beta$ are well accounted for by the Sonine approximation. The quantitative agreement is good as long as the magnitude of the cumulants is small, thus validating the Sonine approximation in those cases. On the other hand, whenever the Sonine approximation predicts values $a_{ij}^{\mathrm{H}}=O(1)$, the approximation is itself signaling its breakdown. This situation, which is similar to that already reported in the case of HS [39]. is especially noteworthy in the cases of $a_{02}^{\mathrm{H}}$ and, to a lesser extent, $a_{11}^{\mathrm{H}}$, and is clearly indicative of the high-velocity tails discussed in Sec. V.2 and confirmed below.

VI.2 High-velocity tails

To further observe the non-Gaussianities of the HCS state, Fig. 6 displays the histograms from simulation data of $\phi^{\mathrm{H}}_{\mathbf{c}}$, $\phi^{\mathrm{H}}_{\mathbf{w}}$, and $\phi^{\mathrm{H}}_{cw}$, for nine combinations of coefficients of restitution ($\alpha=0.9$, $0.7$, $0.2$, and $\beta=0.5$, [math], $-0.5$). Except in the case of $\phi^{\mathrm{H}}_{\mathbf{c}}$ for $\alpha=0.9$ (where $a_{20}^{\mathrm{H}}$ is small), the deviations from the Maxwellian tail are quite apparent. In fact, the high-velocity tails observed in Fig. 6 are consistent with an exponential tail for $\phi^{\mathrm{H}}_{\mathbf{c}}$ and power-law tails for $\phi^{\mathrm{H}}_{\mathbf{w}}$ and $\phi^{\mathrm{H}}_{cw}$, in agreement with the analysis in Sec. V.2.

A more quantitative test is presented in Fig. 7, where the three cases with $\alpha=0.7$ have been selected and the straight lines representing the asymptotic tails are included. The theoretical predictions for the exponents derived in Sec. V.2 (with additional Maxwellian estimates for $\mu_{20}^{\mathrm{H}}$ and $\langle\!\langle c_{12}\rangle\!\rangle^{\mathrm{H}}$) agree reasonably well with the fitted values, except for $\beta=-0.5$, in which case the actual decays are slower than predicted.

Figure 8 shows the exponents $\gamma_{c}$, $\gamma_{w}$, and $\gamma_{cw}$ as functions of $\beta$ (for $\alpha=0.9$, $0.7$, $0.2$) and $\alpha$ (for $\beta=0.5$, [math], $-0.5$). There exists very good agreement between the theoretical estimates and the fitting simulation values in the case of the exponent $\gamma_{w}$, which seems to worsen as $\beta$ decreases. However, in the case of $\gamma_{c}$ and $\gamma_{cw}$, the agreement is mainly qualitative. This might be due to the fact that the tails of $\phi_{\mathbf{c}}^{\mathrm{H}}$ and $\phi_{cw}^{\mathrm{H}}$ are much less populated than that of $\phi_{\mathbf{w}}^{\mathrm{H}}$ (see Figs. 6 and 7) and, therefore, it is much more difficult to reach values of $c$ and $c^{2}w^{2}$ high enough to accurately measure the exponents $\gamma_{c}$ and $\gamma_{cw}$ in the simulations. If that were the case, then the values of $\gamma_{c}$ and $\gamma_{cw}$ empirically determined would characterize an intermediate velocity regime previous to the true asymptotic behavior. Of course, one cannot discard that our analysis becomes more limited as $\beta$ decreases.

Before closing this section, it is worth remarking the excellent mutual agreement between DSMC and EDMD results. There are, however, some small deviations for low values of $\alpha$, which might be a consequence of the smaller number of disks in the EDMD simulations and also a reflection of possible violations of the molecular chaos ansatz in those highly dissipative systems [71].

VII Conclusions

In this paper, we have studied low-density, monodisperse, and homogeneous granular gases of hard disks and hard spheres from a kinetic-theory point of view, using a general framework to express the results in terms of the number of translational and rotational degrees of freedom, ${d_{t}}$ and ${d_{r}}$, respectively. Special attention has been paid to the non-Gaussian features of the HCS, as measured by the fourth-order cumulants and the high-velocity tails of the marginal distributions. The theory has been complemented by DSMC and EDMD computer simulations.

The theoretical approach is based on the Boltzmann equation. First, we have expressed the collisional moments as formally exact functions of the parameters of the system ($\alpha$, $\beta$, and $\kappa$) and two-body averages. Next, we have employed a Grad–Sonine expansion of the complete one-body VDF, Eq. (25). Then, in analogy to Ref. [39], we have defined the Sonine approximation from the truncation of the Sonine expansion beyond the first nontrivial cumulants defined in Eq. (31). This contrasts with the Maxwellian approximation, which is based on approximating the VDF by a two-temperature Maxwellian distribution, i.e., $\phi\approx\phi_{\mathrm{M}}$.

Within the Sonine approximation, and neglecting quadratic terms, the relevant collisional moments have been evaluated, thus recovering previous results for hard spheres (${d_{t}}={d_{t}}=3$) [39], and obtaining results for hard disks (${d_{t}}=2$, ${d_{r}}=1$), as presented in Table 1. Cumulant-linearization in Eqs. (20a) and (32) allows us to deal with a closed set of differential equation for the evolution of the rotational-to-translational temperature ratio ($\theta$) and the fourth-order cumulants. Analogously, the stationary HCS values in the Sonine approximation have been obtained by linearization in Eqs. (33). As a consistency test, we have checked in Appendix A that the HCS is linearly stable with respect to homogeneous and isotropic perturbations. The HCS quantities have been shown in Figs. 2 and 3 for uniform disks ($\kappa=\frac{1}{2}$) and uniform spheres ($\kappa=\frac{2}{5}$), respectively. At a qualitative level, their dependence on $\alpha$ and $\beta$ is very similar for disks and spheres, but the values are generally more extreme in the former system than in the latter. In both cases, the kurtosis for the angular velocity, $a_{02}^{\mathrm{H}}$, reaches values of $O(1)$ in a lobular region of the parameter space with a vertex at $(\alpha,\beta)=(1,-1)$, thus announcing a breakdown of the Sonine approximation in that region.

Moreover, the non-Gaussianities of the HCS have been studied not only in the context of the first nontrivial cumulants, but also analyzing the tails of the marginal VDF $\phi_{\mathbf{c}}^{\mathrm{H}}(\mathbf{c})$, $\phi_{\mathbf{w}}^{\mathrm{H}}(\mathbf{w})$, and $\phi_{cw}^{\mathrm{H}}(c^{2}w^{2})$ defined by Eqs. (36). Using previous methods developed for the smooth case [58, 48, 49], which are based on the prevalence of the collisional loss term with respect to the gain term, we have obtained the expected exponential tail $\phi_{\mathbf{c}}^{\mathrm{H}}(\mathbf{c})\sim e^{-\gamma_{c}c}$, with formally the same expression for the exponent coefficient $\gamma_{c}$ as in the smooth case [see Eq. (44)]. On the other hand, we have found much slower scale-free decays $\phi_{\mathbf{w}}^{\mathrm{H}}(\mathbf{w})\sim w^{-\gamma_{w}}$ and $\phi_{cw}^{\mathrm{H}}(c^{2}w^{2})\sim(c^{2}w^{2})^{-\gamma_{cw}}$, with exponents given by Eqs. (48) and (50), respectively. These algebraic tails of $\phi_{\mathbf{w}}^{\mathrm{H}}(\mathbf{w})$ and $\phi_{cw}^{\mathrm{H}}(c^{2}w^{2})$ explain the relatively large values attained by the cumulants $a_{02}^{\mathrm{H}}$ and $a_{11}^{\mathrm{H}}$, especially in the hard-disk case, and predict divergences in higher-order cumulants which recall the ones already observed in the case of the three-dimensional inelastic and rough Maxwell model [72].

To test the theoretical predictions, we have run DSMC and EDMD computer simulations for hard disks (with $\kappa=\frac{1}{2}$), as described in Sec. VI. First, the quantities $\theta^{\mathrm{H}}$, $a_{20}^{\mathrm{H}}$, $a_{02}^{\mathrm{H}}$, and $a_{11}^{\mathrm{H}}$ have been studied for different values of $\alpha$ and $\beta$, as depicted in Figs. 4 and 5. The agreement between the Sonine approximation and simulation is rather good, except when the values of the cumulants are not small. Even in those cases, it is remarkable that the Sonine approximation reproduces qualitatively well the shape of the curves. Second, we have extended the comparison to the three marginal distributions in Figs. 6 and 7, finding that the predicted exponential tail of $\phi_{\mathbf{c}}^{\mathrm{H}}(\mathbf{c})$ and algebraic tails of $\phi_{\mathbf{w}}^{\mathrm{H}}(\mathbf{w})$ and $\phi_{cw}^{\mathrm{H}}(c^{2}w^{2})$ are supported by simulation data. The theoretical and fitting exponents have been compared in Fig. 8, where a good agreement for $\gamma_{w}$ has been observed, while the agreement is more qualitative for $\gamma_{c}$ and $\gamma_{cw}$. This might be due to a lack of statistically reliable simulation data in the high-velocity regime.

To sum up, the HCS VDF of a monodisperse granular gas of inelastic and rough hard particles is, in general, strongly non-Maxwellian. Moreover, the non-Gaussianities exposed in this paper might solve some inconsistencies reported in the stability analysis of Navier–Stokes hydrodynamics from a Maxwellian approximation in hard-disk systems [46] and improve the predictions of the inelastic hard-sphere model for real experimental systems, such as the one of Ref. [57]. As a follow-up of the study presented in this paper, we plan to extend it to driven hard-disk systems (in analogy to a previous study for hard spheres [41]), whose dynamics has very interesting implications [73, 74]. Finally, we hope this paper could stimulate further research in all these issues, not only from theoretical and simulation points of view, but also from experimental setups.

The data that support the findings of this study are openly available in Ref. [75].

Acknowledgements.

The authors acknowledge financial support from Grant No. PID2020-112936GB-I00 funded by MCIN/AEI/10.13039/501100011033, and from Grant No. IB20079 funded by Junta de Extremadura (Spain) and by ERDF, “A way of making Europe.” A.M. is grateful to the Spanish Ministerio de Ciencia, Innovación y Universidades for a predoctoral fellowship No. FPU2018-3503. The authors are grateful to the computing facilities of the Instituto de Computación Científica Avanzada of the University of Extremadura (ICCAEx), where the simulations were run.

Appendix A Linear stability analysis of the homogeneous cooling state

In this Appendix we show that, within the Sonine approximation, the HCS for hard disks is linearly stable under uniform and isotropic perturbations. Let us define the time-dependent set of perturbed quantities:

[TABLE]

Insertion into the Sonine approximation versions of Eqs. (20a) and (32a)–(32c), and linearization around the HCS values, yield

[TABLE]

where $\mathsf{L}$ is a constant matrix, its four eigenvalues, $\{\ell_{i};i=1,2,3,4\}$, determining the evolution of $\delta\mathbf{Y}(s)$ from an arbitrary initial perturbation $\delta\mathbf{Y}(0)$.

The dependence of the four eigenvalues on the coefficients of restitution is displayed in Fig. 9 for uniform disks ($\kappa=\frac{1}{2}$). As can be seen, the real parts are always positive, thus signaling the linear stability and attractor character of the HCS under uniform perturbations, as expected on physical grounds. In turn, since the cumulant-linearization scheme within the Sonine approximation is not univocally defined [51, 55], the fact that we get $\mathrm{Re}(\ell_{i})>0$ reinforces the reliability of the linearization criterion applied to the right-hand sides of Eqs. (32) and (33).

The imaginary parts plotted in Fig. 9(e) show the regions of the parameter space where the decay toward the HCS is oscillatory. In this respect, the plane $(\alpha,\beta)$ turns out to be split into three disjoint regions: a region where $(\ell_{1},\ell_{2})$ make a pair of complex conjugates but $(\ell_{3},\ell_{4})$ are real, a region where $(\ell_{3},\ell_{4})$ make a pair of complex conjugates but $(\ell_{1},\ell_{2})$ are real, and, finally, the blank region in Fig. 9(e), where the four eigenvalues are real.

Appendix B Consistency of the high-velocity tails

In Sec. V.2, the high-velocity tails of the HCS marginal distributions $\phi_{\mathbf{c}}^{\mathrm{H}}(\mathbf{c})$, $\phi_{\mathbf{w}}^{\mathrm{H}}(\mathbf{w})$, and $\phi_{cw}^{\mathrm{H}}(c^{2}w^{2})$ were obtained by assuming Eq. (42). Here, we test the self-consistency of that assumption.

B.1 $\phi_{\mathbf{c}}^{\mathrm{H}}(\mathbf{c})$

Let us insert Eqs. (44) into the ratio resulting from the replacement $\phi^{\mathrm{H}}({\widetilde{\mathbf{\Gamma}}})\to\phi_{\mathbf{c}}^{\mathrm{H}}(\mathbf{c})$ in Eq. (42):

[TABLE]

Assuming $c_{1}\gg\{1,c_{2},w_{1},w_{2}\}$ in the inverse binary collisional rules, Eq. (9a), and after some algebra, one gets

[TABLE]

where $\vartheta_{c}=\cos^{-1}\left|\widehat{\mathbf{c}}_{1}\cdot\widehat{\bm{\sigma}}\right|$. Therefore, the exponent in Eq. (53) is strictly negative, except for smooth particles ($\beta=-1$) and grazing collisions ($\cos\vartheta_{c}=0$). Thus, apart from those cases with zero Lebesgue measure, $\lim_{c_{1}\to\infty}{\phi_{\mathbf{c}}^{\mathrm{H}}(\mathbf{c}_{1}^{\prime\prime})\phi_{\mathbf{c}}^{\mathrm{H}}(\mathbf{c}_{2}^{\prime\prime})}/{\phi_{\mathbf{c}}^{\mathrm{H}}(\mathbf{c}_{1})\phi_{\mathbf{c}}^{\mathrm{H}}(\mathbf{c}_{2})}=0$.

B.2 $\phi_{\mathbf{w}}^{\mathrm{H}}(\mathbf{w})$

In the rotational case, from Eqs. (48) we have

[TABLE]

Let us take $w_{1}\gg\{1,w_{2},c_{1},c_{2}\}$. Then,

[TABLE]

where $\vartheta_{w}=\cos^{-1}\left|\widehat{\mathbf{w}}_{1}\cdot\widehat{\bm{\sigma}}\right|$. Therefore, $\lim_{w_{1}\to\infty}\phi^{\mathrm{H}}_{\mathbf{w}}(\mathbf{w}_{1}^{\prime\prime})\phi^{\mathrm{H}}_{\mathbf{w}}(\mathbf{w}_{2}^{\prime\prime})/\phi^{\mathrm{H}}_{\mathbf{w}}(\mathbf{w}_{1})\phi^{\mathrm{H}}_{\mathbf{w}}(\mathbf{w}_{2})=0$, except if $\sin\vartheta_{w}=0$, which has zero Lebesgue measure in its continuous domain.

B.3 $\phi_{cw}^{\mathrm{H}}(x)$

From Eqs. (50), one has

[TABLE]

If both $c_{1}$ and $w_{1}$ are much larger than $\{1,c_{2},w_{2}\}$, it is possible to obtain

[TABLE]

where $\vartheta_{cw}=\cos^{-1}\left[\widehat{\mathbf{c}}_{1}\cdot\left(\widehat{\bm{\sigma}}\times\widehat{\mathbf{w}}_{1}\right)\right]$. Thus, $\lim_{x_{1}\to\infty}{\phi_{cw}^{\mathrm{H}}(x_{1}^{\prime\prime})\phi_{cw}^{\mathrm{H}}(x_{2}^{\prime\prime})}/{\phi_{cw}^{\mathrm{H}}(x_{1})\phi_{cw}^{\mathrm{H}}(x_{2})}=0$, except if $\sin\vartheta_{c}=\sin\vartheta_{w}=\cos\vartheta_{cw}=0$, which again have zero Lebesgue measure.