# Deviation from Maxwell Distribution in Granular Gases with Constant Restitution Coefficient

## Abstract

We analyze the velocity distribution function of force-free granular gases in the regime of homogeneous cooling when deviations from the Maxwellian distribution may be accounted only by leading term in the Sonine polynomial expansion. These are quantified by the magnitude of the coefficient $`a_2`$ of the second term of the expansion. In our study we go beyond the linear approximation for $`a_2`$ and observe that there are three different values (three roots) for this coefficient which correspond to a scaling solution to the Boltzmann equation. The stability analysis performed showed, however, that among these three roots only one corresponds to a stable scaling solution. This is very close to $`a_2`$, obtained in previous studies in a linear with respect to $`a_2`$ approximation.

PACS numbers: 81.05.Rm, 36.40.Sx, 51.20.+d, 66.30.Hs

The granular gases, i.e., rarefied systems composed of inelastically colliding particles have been of particular interest during the last decade (e.g. ). Compared to gases of elastically colliding particles, the dissipation of energy at inelastic collisions leads to some novel phenomena in these systems. One can mention clustering (e.g. ), formation of vortex patterns (e.g.), etc. Before clustering starts, the granular gas being initially homogeneous, keeps for some time its homogeneity, although its temperature permanently decreases. This regime is called the homogeneous cooling regime (HC).

In the present study we address the properties of the velocity distribution of granular particles in the regime of HC, such as the deviation from the Maxwellian distribution and the stability of the distribution function. We assume that the restitution coefficient $`ϵ`$ does not depend on the impact velocity, i.e. that $`ϵ=\mathrm{const}`$. The properties of the velocity distribution for the system with the impact-velocity dependent coefficient of restitution (e.g. ) will be addressed elsewhere .

It is well known that granular gases in the HC regime do not reveal Maxwellian distribution (e.g.). The high-velocity tail is overpopulated , while the main part of the distribution is described by the sum of the Maxwellian one and the correction to it, written in terms of the Sonine polynomial expansion (e.g. ). Usually only the leading, second term, in this expansion is taken into account , moreover in previous studies  only linear analysis with respect to the coefficient $`a_2`$, which refers to this second term has been performed. Finding that $`a_2`$, obtained within the linear approximation, is small, the authors of Ref. conclude a posteriori that the the linear approximation is valid.

In our approach we also assume that one can restrict oneself to the leading term in the Sonine polynomial expansion and ignore the other. However we go beyond the linear approximation with respect to the coefficient $`a_2`$ and perform complete analysis within this level of the system description. We observed that there are three different values of $`a_2`$ exist which correspond to the scaling solution of the Boltzmann equation. The stability analysis for the velocity distribution function shows, however, that only one value of $`a_2`$ corresponds to a physically acceptable stable scaling solution. The stable solution is close to the result previously obtained within the linear analysis .

To introduce notations and specify the problem we briefly sketch the derivation of the coefficient $`a_2`$ , in accordance with the approach developed in Ref..

So far we introduce the (time-dependent) temperature $`T(t)`$, and thermal velocity $`v_0(t)`$, which are related to the velocity distribution function $`f(𝐯,t)`$ as

$$\frac{3}{2}nT(t)=𝑑𝐯\frac{v^2}{2}f(𝐯,t)=\frac{3}{2}nv_0^2(t)$$
(1)   
here $`n`$ is the number density of the granular gas, the particles are assumed to be of a unit mass ($`m=1`$), and (1) is written for 3D-systems. The inelasticity of collisions is characterized by the coefficient of the normal restitution $`ϵ`$ (here we consider smooth particles), which relates after-collisional velocities $`𝐯_1^{}`$, $`𝐯_2^{}`$ to the pre-collisional ones, $`𝐯_1`$, $`𝐯_2`$ as:

$$𝐯_{1/2}^{}=𝐯_{1/2}\frac{1}{2}(1+ϵ)(𝐯_{12}𝐞)𝐞$$
(2)   
where $`𝐯_{12}=𝐯_1𝐯_2`$ is the relative velocity, the unit vector $`𝐞=𝐫_{12}/|𝐫_{12}|`$ gives the direction of the vector $`𝐫_{12}=𝐫_1𝐫_2`$ at the instant of the collision. The time-evolution of the velocity distribution function is subjected to the Enskog-Boltzmann equation, which for the force-free case reads :

$`{\displaystyle \frac{}{t}}f(𝐯,t)=g_2(\sigma )\sigma ^2{\displaystyle 𝑑𝐯_2𝑑𝐞\mathrm{\Theta }(𝐯_{12}𝐞)|𝐯_{12}𝐞|}`$ (3)   
$`\times \left\{{\displaystyle \frac{1}{ϵ^2}}f(𝐯_1^{},t)f(𝐯_2^{},t)f(𝐯_1,t)f(𝐯_2,t)\right\}`$ (4)   
where $`\sigma `$ is the diameter of particles, $`g_2(\sigma )=(2\eta )/2(1\eta )^3`$ ($`\eta =\frac{1}{6}\pi n\sigma ^3`$ is packing fraction) denotes the contact value of the two-particle correlation function , which accounts for the increasing collision frequency due to the excluded volume effects; $`\mathrm{\Theta }(x)`$ is the Heaviside function. The velocities $`𝐯_1^{}`$ and $`𝐯_2^{}`$ refer to the precollisional velocities of the so-called inverse collision, which results with $`𝐯_1`$ and $`𝐯_2`$ as the after-collisional velocities. The factor $`1/ϵ^2`$ in the gain term appears respectively from the Jacobian of the transformation $`d𝐯_1^{}d𝐯_2^{}d𝐯_1d𝐯_2`$ and from the relation between the lengths of the collisional cylinders $`ϵ|𝐯_{12}^{}𝐞|dt=|𝐯_{12}𝐞|dt`$ .

Assuming that the velocity distribution function is of a scaling form:

$$f(𝐯,t)=\frac{n}{v_0^3(t)}\stackrel{~}{f}(𝐜)$$
(5)   
one can show, that the scaling function satisfies the time-independent equation :

$$\frac{\mu _2}{3}\left(3+c_1\frac{}{c_1}\right)\stackrel{~}{f}(𝐜)=\stackrel{~}{I}(\stackrel{~}{f},\stackrel{~}{f})$$
(6)   
with the dimensionless collision integral:

$`\stackrel{~}{I}(\stackrel{~}{f},\stackrel{~}{f})={\displaystyle 𝑑𝐜_2𝑑𝐞\mathrm{\Theta }(𝐜_{12}𝐞)|𝐜_{12}𝐞|}`$ (7)   
$`\times \left\{ϵ^2\stackrel{~}{f}(𝐜_1^{})\stackrel{~}{f}(𝐜_2^{})\stackrel{~}{f}(𝐜_1)\stackrel{~}{f}(𝐜_2)\right\}`$ (8)   
and with the moments of the dimensionless collision integral :

$$\mu _p𝑑𝐜_1c_1^p\stackrel{~}{I}(\stackrel{~}{f},\stackrel{~}{f}),$$
(9)   
while the time-evolution of temperature reads:

$$dT/dt=(2/3)BT\mu _2$$
(10)   
where $`B=B(t)v_0(t)g_2(\sigma )\sigma ^2n`$.

To proceed we use the Sonine polynomial expansion for the velocity distribution function 

$$\stackrel{~}{f}(𝐜)=\varphi (c)\left\{1+\underset{p=1}{\overset{\mathrm{}}{}}a_pS_p(c^2)\right\}$$
(11)   
where $`\varphi (c)\pi ^{d/2}\mathrm{exp}(c^2)`$ is the Maxwellian distribution and the first few Sonine polynomials read: $`S_0(x)=1`$, $`S_1(x)=x^2+\frac{3}{2}`$, $`S_2(x)=\frac{x^2}{2}\frac{5x}{2}+\frac{15}{8}`$, etc. Multiplying both sides of Eq. (6) with $`c_1^p`$ and integrating over $`d𝐜_1`$, we obtain :

$$\frac{\mu _2}{3}pc^p=\mu _p$$
(12)   
where integration by parts has been performed and where we define

$$c^pc^p\stackrel{~}{f}(𝐜,t)𝑑𝐜.$$
(13)   
The odd moments $`c^{2n+1}`$ are zero, while the even ones, $`c^{2n}`$, may be expressed in terms of $`a_k`$ with $`0kn`$. Calculations show that $`c^2=\frac{3}{2}`$, implying $`a_1=0`$, according to the definition of the temperature (1) (e.g. ), and that $`c^4=\frac{15}{4}\left(1+a_2\right)`$.

Now we assume, that the dissipation is not large, so that the deviation from the Maxwellian distribution may be accurately described only by the second term in the expansion (11) with all high-order terms with $`p>2`$ discarded. Then (12) is an equation for the coefficient $`a_2`$. Using the above results for $`c^2`$ and $`c^4`$ it is easy to show that Eq. (12) converts for $`p=2`$ into identity, while for $`p=4`$ it reads:

$$5\mu _2\left(1+a_2\right)\mu _4=0$$
(14)   
The coefficients $`\mu _p`$ may be expressed in terms of $`a_2`$ due to the definition (9) and the assumption $`\stackrel{~}{f}=\varphi (c)[1+a_2S_2(c^2)]`$. Using the properties of the collision integral one obtains for $`\mu _p`$ :

$`\mu _p={\displaystyle \frac{1}{2}}{\displaystyle 𝑑𝐜_1𝑑𝐜_2𝑑𝐞\mathrm{\Theta }(𝐜_{12}𝐞)|𝐜_{12}𝐞|\varphi (c_1)\varphi (c_2)}`$ (15)   
$`\left\{1+a_2\left[S_2(c_1^2)+S_2(c_2^2)\right]+a_2^2S_2(c_1^2)S_2(c_2^2)\right\}\mathrm{\Delta }(c_1^p+c_2^p)`$ (16)   
where $`\mathrm{\Delta }\psi (𝐜_i)\left[\psi (𝐜_i^{})\psi (𝐜_i)\right]`$ denotes change of some function $`\psi (𝐜_i)`$ in a direct collision. Calculations, similar to that, described in , yield the following result (some detail are given in ):

$$\mu _2=\sqrt{2\pi }(1ϵ^2)\left(1+\frac{3}{16}a_2+\frac{9}{1024}a_2^2\right)$$
(17)   
and

$$\mu _4=4\sqrt{2\pi }\left\{T_1+a_2T_2+a_2^2T_3\right\}$$
(18)   
with

$`T_1={\displaystyle \frac{1}{4}}(1ϵ^2)\left({\displaystyle \frac{9}{2}}+ϵ^2\right)`$ (19)   
$`T_2={\displaystyle \frac{3}{128}}(1ϵ^2)(69+10ϵ^2)+{\displaystyle \frac{1}{2}}(1+ϵ)`$ (20)   
$`T_3={\displaystyle \frac{1}{64}}(1+ϵ)+{\displaystyle \frac{1}{8192}}(1ϵ^2)(930ϵ^2)`$ (21)   
The coefficients $`\mu _2`$ and $`\mu _4`$ were provided in Ref.  up to terms of the order of $`𝒪(a_2)`$. One obtains the coefficient $`a_2`$ in the Sonine polynomial expansion in this approximation by substituting (17,18) into (14) and discarding in Eqs. (17,18) all terms of the order of $`𝒪(a_2^2)`$:

$$a_2^{\mathrm{NE}}=\frac{16(1ϵ)(12ϵ^2)}{8117ϵ+30ϵ^2(1ϵ)}$$
(22)   
Calculations including the next order terms $`𝒪(a_2^2)`$ in the coefficients $`\mu _2`$ and $`\mu _4`$ show that Eq. (14) is a cubic equation, which for physical values of $`ϵ`$, $`0ϵ1`$, has three different real roots, as it shown on Fig. 1.    FIG. 1.: The left hand side of Eq. 14 over $`a_2`$ for $`ϵ=0.8`$. Obviously Eq. 14 has three real solutions. 

Although the cubic equation may be generally solved, the resultant expressions for the roots are too cumbersome to be written explicitly. However, one of the roots (the middle one) is rather small and close to that given by Eq. (22), obtained within the linear approximation. This suggests the perturbative solution of the cubic equation near this root:

$$a_2=a_2^{\mathrm{NE}}\left[1\frac{1005(1ϵ^2)4096T_3}{6080(1ϵ^2)4096T_2}a_2^{\mathrm{NE}}+\mathrm{}\right]$$
(23)   
where we do not write explicitly terms of the order $`𝒪\left([a_2^{\mathrm{NE}}]^3\right)`$ and high-order terms. In Fig. 2 the dependence of $`a_2^{\mathrm{NE}}`$ and of the corresponding improved value $`a_2`$ are shown as a function of the restitution coefficient $`ϵ`$. As one can see from Fig. 2 the maximal deviation between these is less than $`10\%`$ at small $`ϵ`$ and decreases as $`ϵ`$ tends to $`1`$.

The other two roots, shown on Fig. 3 are of the order of 1 or 10, i.e. are not small. Physically, this means that one can not cut the Sonine polynomial expansion in this case at the second term and next order terms are not negligible to be discarded. Taking into account the next order terms, i.e., releasing the assumption that $`a_p0`$ for $`p>2`$, breaks down the above analysis, since the coefficients $`\mu _2`$, $`\mu _4`$ occur to be dependent not only on $`a_2`$, but on $`a_3`$, $`a_4,\mathrm{}`$ as well. Thus the occurrence of several roots for the $`a_2`$, found within the above approach, which satisfy the conditions required by the scaling ansatz (5) does not imply the existence of several different scaling solutions. Nevertheless such possibility may not be completely excluded. If one assumes that few scaling distributions of the velocity may realize, depending on the initial conditions at which the HC state has been prepared, a natural question arises: Whether the particular scaling solution is stable with respect to small perturbations, and what is the domain of attraction of this particular scaling solution in some parametric space.    FIG. 3.: The other two solutions for second Sonine coefficient $`a_2`$ of Eq. 14 over the coefficient of restitution $`ϵ`$.  

Certainly, the stability problem is very complicated to be solved in general. Therefore, we restrict ourselves to the stability analysis of the scaling distribution (5) where the scaling function $`\stackrel{~}{f}(𝐜)`$ has nonzero value of the coefficient $`a_2`$, while the other coefficients $`a_p`$ with $`p>2`$ are negligibly small. (For this scaling solution our above results for the coefficients $`\mu _2`$, $`\mu _4`$ are valid). Moreover, we assume, that small perturbations of the (vanishingly small) coefficients $`a_p`$ with $`p>2`$ do not influence the stability of the distribution, and analyze the stability only with respect to variation of the coefficient $`a_2`$.

To analyze the stability of the velocity distribution we write it in a more general form:

$$f(𝐯,t)=\frac{n}{v_0^3(t)}\stackrel{~}{f}(𝐜,t)$$
(24)   
which leads, as it easy to show, to the following generalization of Eq. (6:

$$\frac{\mu _2}{3}\left(3+c_1\frac{}{c_1}\right)\stackrel{~}{f}(𝐜,t)+B^1\frac{}{t}\stackrel{~}{f}(𝐜,t)=\stackrel{~}{I}(\stackrel{~}{f},\stackrel{~}{f})$$
(25)   
with the collisional integral and coefficients $`\mu _p`$ being now time-dependent. The quantities $`c^p`$ also depend now on time. The temperature, however, evolves still according to (10).

Using $`\stackrel{~}{f}=\varphi (c)[1+a_2(t)S_2(c^2)]`$ and performing essentially the same manipulations which led before to Eq. (14), we arrive at the following equation for the coefficient $`a_2(t)`$:

$$\dot{a}_2(4/3)B\mu _2\left(1+a_2\right)+(4/15)B\mu _4=0$$
(26)   
with $`\mu _2`$, $`\mu _4`$ still given by (17,18), but with the time-dependent coefficient $`a_2(t)`$. Writing the above value $`B(t)`$ as

$`B(t)`$ $`=`$ $`(8\pi )^{1/2}\tau _c(0)^1u(t)^{1/2}`$ (27)   
$`\tau _c(0)^1`$ $``$ $`4\pi ^{1/2}g_2(\sigma )\sigma ^2nT_0^{1/2},`$ (28)   
where $`\tau _c(0)`$ is related to the initial mean-collision time at the initial temperature $`T_0`$, and $`u(t)T(t)/T_0`$ is the reduced temperature, we recast Eq. (26) into the form:

$$\frac{da_2}{d\widehat{t}}=\frac{\sqrt{2/\pi }}{15}u^{1/2}F(a_2)$$
(29)   
where $`\widehat{t}`$ is the reduced time, measured in units of $`\tau _c(0)`$, and where we define a function:

$$F(a_2)5\mu _2(1+a_2)\mu _4.$$
(30)   
The form of the function $`F(a_2)`$ for some particular value of $`ϵ`$ is shown on Fig. 1. This form of $`F(a_2)`$ persists for all physical values of the restitution coefficient, $`0ϵ1`$. This has three different roots, $`F(a_2^{(i)})=0`$, $`i=1,2,3`$, which makes $`da_2/dt`$ vanish yielding the scaling form for the solution of the Enskog-Boltzmann equation. The stability of the scaling solution, corresponding to $`a_2^{(i)}`$ requires for the derivative $`dF/da_2`$, taken at $`a_2^{(i)}`$ to be negative, since only in this case a small deviation $`a_2a_2^{(i)}`$ from $`a_2^{(i)}`$, corresponding to a scaling solution will decay with time. As one can see from Fig. 1 only the middle root, which corresponds to small values of $`a_2`$, and is close to $`a_2^{\mathrm{NE}}`$, predicted by linear theory , has negative $`dF/da_2`$, and thus is stable. We also observed that for any $`0ϵ1`$ the point $`a_2=0`$ belongs to the attractive interval of this stable root. Naturally, this means that initial Maxwellian distribution will relax to the non-Maxwellian with $`a_2a_2^{\mathrm{NE}}`$.

Note that relaxation of any (small) perturbation to this value of $`a_2`$ occurs, as it follows from Eq. (29), on the collision time-scale, i.e., practically “immediately” on the time-scale which describes the evolution of the temperature. Therefore we conclude, that the scaling solution of the Enskog-Boltzmann equation with $`a_2`$ corresponding to the middle root of the function $`F(a_2)`$, given with a high accuracy by Eqs. (23,22), and with negligibly small other coefficients $`a_3`$, $`a_4,\mathrm{}`$ of the Sonine polynomial expansion is a stable one with respect to (relatively) small perturbations.

In conclusion, we analyzed the velocity distribution function of the granular gas with a constant restitution coefficient at the regime of the homogeneous cooling. We assume that the deviations from the Maxwellian distribution may be described using only the leading term in the Sonine polynomial expansion, with all other high-order terms discarded. In this approach the deviations from the Maxwellian distribution are completely characterized by the magnitude of the coefficient $`a_2`$ of the leading term. We go beyond previous linear theories and perform a complete analysis (on the level of the description chosen), without discarding any nonlinear with respect to $`a_2`$ terms.

Performing the stability analysis of the scaling solution of the Enskog-Boltzmann equation we observe, that only one value of $`a_2`$, obtained within our nonlinear analysis corresponds to a stable scaling solution. We also report a corrections for this value of $`a_2`$ with respect to the previous result of the linear theory. This corrections are small (less than $`10\%`$) for all values of the restitution coefficient $`ϵ`$ and vanishes as $`ϵ`$ tends to unity in the elastic limit.