format: text

text:
# Asymptotic solutions of the nonlinear Boltzmann equation for dissipative systems

## 1 Introduction

Classic kinetic theory Waldm58 ; Chapman ; resibois ; cercignani-book ; Uhl+Fo63  deals with elastic particles with energy conserving dynamics. The system is described by the single particle distribution function, whose time evolution is governed by the nonlinear Boltzmann equation. The asymptotic states of such systems follow the universal laws of thermodynamics, and the distribution function is the Maxwell Boltzmann distribution. This scenario does not apply to dissipative systems, where energy is lost in inelastic interactions.

In elastic systems the approach to asymptotic states is characterized by a kinetic stage of rapid relaxation in velocity space to a locally homogeneous equilibrium state, followed by a hydrodynamic stage of slow approach to a globally homogeneous equilibrium state. The time scale in the kinetic stage is the mean free time $`t_{mf}`$ between collisions. In the kinetic theory of inelastic systems campbell ; jenkins ; savage ; haff ; goldhirsch ; goldstein-shapiro ; goldhirsch-zanetti  the type of decay depends on the energy supply to the dissipative system. Without energy supply there is first a kinetic stage of rapid relaxation on the time scale $`t_{mf}`$ to a locally homogeneous adiabatic state, the homogeneous cooling state, described by scaling or similarity solutions with a slowly changing parameter, at least for weakly inelastic systems. With energy supply the evolution is more similar to the elastic case with, however, equilibrium states replaced by non-equilibrium steady states. The velocity distributions in these adiabatic or steady states are very different from a Maxwell Boltzmann distribution. The subsequent stage of evolution involves transport phenomena and complex hydrodynamic phenomena of clustering and pattern formation goldhirsch-zanetti ; Nagel .

The interest in granular matter in general has strongly stimulated new developments in the kinetic theory of granular fluids and gases, which show surprising new physics. A granular fluid is a collection of small or large macroscopic particles, with short range repulsive hard core interactions, in which energy is lost in inelastic collisions, and the system cools when not driven. When rapidly driven, gravity can be neglected. The dynamics is based on binary collisions and ballistic motion between collisions, which conserve total momentum. So these systems can be considered to be a granular fluid or gas. The prototypical model for these so-called rapid granular flows is a fluid or gas of perfectly smooth mono-disperse inelastic hard spheres, and its non-equilibrium behavior can be described by the nonlinear Boltzmann equation campbell ; jenkins ; savage ; haff ; goldhirsch ; goldstein-shapiro ; goldhirsch-zanetti . The inelastic collisions are modeled by a coefficient of restitution $`\alpha \left(0<\alpha <1\right)`$, where $`\left(1\alpha ^2\right)`$ measures the degree of inelasticity.

This review focuses on the first stage of evolution, and studies the velocity distribution $`F(v,t)`$ in spatially homogeneous states of inelastic systems. For that reason most of the citations, given in this article, only refer to kinetic theory studies of $`F(v,t)`$, and not to studies of transport properties. The revival WM96 ; brey-dufty-santos ; esipov ; TvN+ME-granmat ; BBRTvW01 ; piasecki ; brey ; MS00 ; brito+huthman ; ignacio-II ; BN-PK-00 ; Bobyl-00 ; CCG-00  in kinetic theory of inelastic systems has been strongly stimulated by the increasing sophistication of experimental techniques SciAm-01 ; GM-exp , which make direct measurements of velocity distributions feasible in non-equilibrium steady states. In this review we also include inelastic generalizations ETB  of the classical repulsive power law interactions Chapman ; resibois ; cercignani-book , which embed both the inelastic hard spheres $`\left(\nu =1\right)`$, as well as the recently much studied Rome1 ; Rome2 ; BN+PK-6-11 ; BN+PK-02 ; ME+RB-EPL ; ME+RB-rapid ; ME+RB-fest ; AB+CC-proof ; Nienhuis ; Droz ; AS+ME-condmat  inelastic Maxwell models $`\left(\nu =0\right)`$ in a single class of models, parametrized by an exponent $`\nu `$. This exponent characterizes the dependence of the collision frequencies on the energy of impact at collision.

In fact, the kinetic theory for such models is of interest in its own right, as the majority of inter-particle interactions in macroscopic systems involve some effects of inelasticity. Our goal is to expose the generic and universal features of the velocity distributions in dissipative fluids, and to compare them with conservative fluids to highlight the differences.

A classical problem in kinetic theory is the possibility of overpopulated high energy tails in velocity distributions bobylev-BKW ; krook+wu , as many physical and chemical processes only occur above a certain energy threshold. Consequently, this old problem has received a new stimulus through a large amount of recent theoretical and experimental studies on tail distributions in many particle systems with inelastic interactions. From the point of view of kinetic theory the intriguing question is, what is the generic feature causing overpopulated tails, possibly even power law tails, in systems of inelastic particles, how does the overpopulation depend on the scattering cross sections, and on the different forms of energy input. The generic feature is the mechanism for overpopulation, and not the specific shape of the tails.

Finally, from the point of view of nonequilibrium steady states, the structure of velocity distributions in elastic and dissipative systems, including the high energy tail, is a subject of continuing research, as the universality of the Gibbs’ state of thermal equilibrium is lacking outside thermal equilibrium, and a possible classification of generic structures would be of great interest in many fields of non-equilibrium statistical mechanics.

The plan of the paper is as follows: in Section 2 we discuss a simple inelastic BGK- or single-relaxation-time model BGK  to illustrate the phenomenon of power law tails. The exponent in the algebraic tail depends qualitatively in the same manner on the degree of inelasticity as in 2- and 3-dimensional Maxwell models. In Section 3 the nonlinear Boltzmann equation is constructed for inelastic generalizations of the classical repulsive power law potentials, to which we refer as Inelastic Repulsive Scatterers or IRS-models. This is done for freely cooling as well as for systems driven by thermostats or heat sources. Section 4 gives a systematic analysis for the energy balance equation, it derives the nonlinear integral equation for the scaling or similarity solution, denoted by $`f\left(c\right)`$, and presents an asymptotic analysis of the high energy tails in the form of stretched Gaussians, $`f\left(c\right)\mathrm{exp}\left[\beta c^b\right]`$ where $`b<2`$. The method used can only be applied to IRS-models with $`\nu >0`$, where the exponent $`b=b\left(\nu \right)`$ is found as a simple function of $`\nu `$. The case of freely cooling Maxwell models ($`\nu =0`$) forms a borderline case, discussed in Section 5. Here algebraic tails, $`f\left(c\right)1/c^{d+a}`$, are found, where the exponent $`a`$ is determined by a transcendental equation. It yields $`a=a\left(\alpha \right)`$ as a function of the degree of inelasticity. We end with some perspectives and conclusions.

## 2 Inelastic BGK Model

### 2.1 Kinetic equations

The goal of this section is to present in the nutshell of a simple Bhatnagar-Gross-Krook (BGK) model BGK  a preview of many of the qualitative features of velocity relaxation in homogeneous systems.

A crude scenario for the relaxation without energy input suggests that the system will cool down due to inelastic collisions, and the velocity distribution $`F(v,t)`$ will approach a Dirac delta function $`\delta \left(\stackrel{}{v}\right)`$ as $`t\mathrm{}`$, while the width or r.m.s. velocity $`v_0\left(t\right)`$ of this distribution, defined as $`v^2=\frac{1}{2}dv_0^2`$, is shrinking. With a constant supply of energy, the system can reach a non-equilibrium steady state (NESS).

To model this evolution we use a simple BGK-type kinetic equation, introduced in 1996 by Brey et al. BMD ,

$$_tF(v,t)D_\stackrel{}{v}^2F(v,t)=\omega _\nu \left(t\right)\left[F(v,t)F_0(v,t)\right].$$
(1)
We have added a heating term, $`D_\stackrel{}{v}^2F(v,t)`$ to the usual BGK equation which represents the heating by a white noise of strength $`D`$ 111For a systematic discussion of driven systems, see Section 3.3. Here the mean collision frequency, $`\omega _\nu \left(t\right)=1/t_{mf}`$, is a function of the r.m.s. velocity $`v_0\left(t\right)`$, chosen as $`\omega _\nu =v_0^\nu `$, in preparation of section 3.1. The kinetic equation describes the relaxation of $`F(v,t)`$ with a time-dependent rate $`\omega _\nu \left(t\right)`$ towards a Maxwellian with a width proportional to $`v_0\left(t\right)`$, defined by

$$F_0(v,t)=\left(\sqrt{\pi }\alpha v_0\right)^d\mathrm{exp}\left[\left(v/\alpha v_0\right)^2\right]\left(\alpha v_0\right)^d\varphi \left(c/\alpha \right),$$
(2)
where $`c=v/v_0`$. The constant $`\alpha `$ ($`0<\alpha <1`$) is related to the inelasticity, $`\gamma =\frac{1}{2}\left(1\alpha ^2\right)`$, of the model, and the totally inelastic limit $`(\alpha 0)`$ is ill-defined in this model, as the mean energy is divergent for $`\alpha =0`$.

### 2.2 Free Cooling ($`D=0`$)

The cooling law of the mean square velocity $`v^2=\frac{1}{2}dv_0^2`$, or the granular temperature $`T=v_0^2`$, is obtained by applying $`𝑑\stackrel{}{v}v^2\left(\mathrm{}\right)`$ to (1) with $`\omega _\nu =v_0^\nu `$. The result is $`_tv_0=\gamma v_0^{\nu +1}`$, yielding

$$v_0\left(t\right)=v_0\left(0\right)/\left[1+\nu \gamma tv_0^\nu (0)\right]^{1/\nu }.$$
(3)
The result is a homogeneous cooling law, $`Tt^{2/\nu }`$, which agrees with Haff’s law haff  for $`\nu =1`$, corresponding to inelastic hard spheres. Note that for negative $`\nu `$ the homogeneous cooling law takes the form,

$$v_0\left(t\right)=v_0\left(0\right)\left[1|\nu |\gamma t/v_0^{|\nu |}(0)\right]^{1/|\nu |},$$
(4)
i.e. for $`t>t_sv_0^{|\nu |}/|\nu |\gamma `$ the r.m.s. velocity and the mean kinetic energy become negative, which is unphysical, and so are the BGK-models with $`\nu <0`$.

As indicated in the introduction, the long time behavior of $`F(v,t)`$ in free cooling is determined by a scaling or similarity solution of the form $`F(v,t)=v_0^d\left(t\right)f\left(v/v_0\left(t\right)\right)`$. We insert this ansatz in (1), eliminate $`_tv_0=\gamma v_0^{\nu +1}`$, and obtain the scaling equation,

$$c\frac{d}{dc}f+\left(d+a\right)f=\frac{a}{\alpha ^d}\varphi \left(\frac{c}{\alpha }\right),$$
(5)
where $`a`$ is defined as,

$$a=1/\gamma =2/\left(1\alpha ^2\right).$$
(6)
We also note that the scaling equation is independent of $`\nu `$. The exact solution of this equation is:

$$f\left(c\right)=\frac{A}{c^{d+a}}+\frac{a}{\alpha ^d\pi ^{d/2}}\left(\frac{1}{c^{d+a}}\right)_0^c𝑑uu^{d+a1}\mathrm{exp}\left[u^2/\alpha ^2\right].$$
(7)
The integration constant $`A`$ is fixed by the normalizations,

$$𝑑\stackrel{}{c}\{1,c^2\}f\left(c\right)=\{1,\frac{1}{2}d\},$$
(8)
where $`d\stackrel{}{c}=\mathrm{\Omega }_dc^{d1}dc`$ with $`\mathrm{\Omega }_d=2\pi ^{d/2}/\mathrm{\Gamma }\left(d/2\right)`$ being the surface area of a $`d`$dimensional hyper-sphere. The normalization integral converges only near $`c0`$ if $`A=0`$. Then the solution (7) with $`A=0`$ is identical to the scaling form, obtained in BMD  for $`\nu =1`$.

For velocities far above thermal, i.e. $`c1`$, the solution has a power law tail,

$$f\left(c\right)\frac{a\alpha ^a}{\pi ^{d/2}}\left(\frac{1}{c^{d+a}}\right)_0^{\mathrm{}}𝑑uu^{d+a1}e^{u^2}=\frac{a\alpha ^a\mathrm{\Gamma }\left(\frac{d+a}{2}\right)}{2\pi ^{d/2}}\left(\frac{1}{c^{d+a}}\right),$$
(9)
with exponent $`a=1/\gamma =2/\left(1\alpha ^2\right)`$, such that $`c^2`$ is bounded for $`\alpha >0`$. The exact solution (7), including its high energy tail, is independent of the exponent $`\nu `$, that determines the energy dependence of the mean collision frequency $`\omega _\nu =v_0^\nu `$ in the BGK model.

As we shall see in Section 5, a similar heavily overpopulated tail, $`f\left(c\right)1/c^{d+a}`$ with $`d>1`$, will also be found in freely cooling Maxwell model with $`\omega _0=1\left(\nu =0\right)`$. There the exponent $`a\left(\alpha \right)`$ takes in the elastic limit ($`\alpha 1`$) the form $`a1/\gamma _0=4d/\left(1\alpha ^2\right)`$. However, in the general class of IRS-models with collision frequency $`\omega _\nu v_0^\nu `$ with $`\nu >0`$ the tails are not given by power laws, but by stretched Gaussians, $`f\left(c\right)\mathrm{exp}\left[\beta c^b\right]`$ with $`0<b=b\left(\nu \right)<2`$. These models will be introduced and discussed in Section 3.2.

### 2.3 NESS ($`D0`$)

Next we extend the result of BMD  to an inelastic BGK-model with $`\nu 0`$, driven by white noise. By applying $`𝑑\stackrel{}{v}v^2`$ to (1) we obtain the temperature balance equation at stationarity,

$$_tv_0^2=4D2\gamma _0v_0^{\nu +2}=0,$$
(10)
where the collisional dissipation, $`2\gamma _0v_0^{\nu +2}`$, is compensated by energy input from the external white noise. Also note that for $`\nu <2`$ the fixed point solution $`v_0\left(\mathrm{}\right)`$ in (10) still exists, but it is unstable ETB . If $`v_0\left(0\right)<v_0\left(\mathrm{}\right)`$, then $`v_0\left(t\right)`$ vanishes as $`t\mathrm{}`$, and if $`v_0\left(0\right)>v_0\left(\mathrm{}\right)`$ then $`v_0\left(t\right)`$ diverges 222Thanks are due to E. Trizac and A. Barrat for pointing out to us that the energy balance equations in this article have stability thresholds that are $`\nu `$dependent..

To obtain the solution $`F(v,\mathrm{})`$ of (1) in the NESS we rescale $`F(v,\mathrm{})=v_0^d\left(\mathrm{}\right)f\left(v/v_0\left(\mathrm{}\right)\right)`$ to the standard width $`c^2=\frac{1}{2}d`$, substitute the rescaled form in the kinetic equation, and eliminate $`D`$, using the stationarity condition (10) as well as (6). This yields the rescaled equation in universal form,

$$\frac{1}{c^{d1}}\frac{d}{dc}c^{d1}\frac{d}{dc}f\left(c\right)2af\left(c\right)=\frac{2a}{\alpha ^d}\varphi \left(\frac{c}{\alpha }\right),$$
(11)
where the normalizations (8) are imposed. The O.D.E. shows that $`f\left(c\right)`$ is independent of the noise strength, $`D`$, and does not contain any dependence on the initial data. This equation can be solved exactly, and more details will be published in ETB . However, for the purpose of this section, we only want to extract from the differential equation (11) the asymptotic form of $`f\left(c\right)`$. In that case we may neglect in (11) the inhomogeneity $`\varphi \left(c/\alpha \right)\mathrm{exp}\left[\left(c/\alpha \right)^2\right]`$ for $`c\alpha `$, and find the asymptotic solution for the BGK-model driven by white noise, i.e.

$$\begin{array}{cc}& f\left(c\right)\mathrm{exp}\left[\beta c^b\right]\hfill \\ & b=1;\beta =\sqrt{2a}=2/\sqrt{1\alpha ^2}.\hfill \end{array}$$
(12)
The constant $`\beta `$ is independent of the parameter $`\nu `$. The exponentially decaying high energy tail is also a ’stretched’ Gaussian, which is overpopulated when compared to a Maxwellian, but the overpopulation is much less heavy than is the freely cooling case (9) with an algebraic tail.

As we shall see in Section 4.3, a similar exponential high energy tail will be found in the white noise driven Maxwell model ($`\nu =0;b=1;\beta =\sqrt{8/\left(1\alpha ^2\right)}`$), but not in the general class of IRS-models, where the stretching exponent $`b`$ takes a value in the interval $`0<b=\overline{b}\left(\nu \right)2`$.

## 3 Basics of inelastic scattering models

### 3.1 Boltzmann equation as a stochastic process

The nonlinear Boltzmann equation for dissipative interactions in the homogeneous cooling state can be put in a broader perspective, that covers both elastic and inelastic collisions, as well as interactions where the scattering of particles is described either by conservative (deterministic) forces, or by stochastic ones. To do so it is convenient to interpret the Boltzmann equation as a stochastic process, similar to the presentations in the classical articles of WaldmannWaldm58 , and Uhlenbeck and Ford Uhl+Fo63 , or in Ulam’s stochastic model Ulam  showing the basics of the approach of a one-dimensional gas of elastic particles towards a Maxwellian distribution.

Consider a spatially homogeneous fluid of elastic or inelastic particles in $`d`$dimensions, specified by their velocities ($`\stackrel{}{v},\stackrel{}{w},\mathrm{}`$), and interacting through binary collisions, $`(\stackrel{}{v},\stackrel{}{w})(\stackrel{}{v}^{},\stackrel{}{w}^{})`$, that are described by transition probabilities. To describe fluids out of equilibrium the total momentum $`\stackrel{}{G}=\frac{1}{2}\left(\stackrel{}{v}+\stackrel{}{w}\right)`$ needs to be conserved in a binary collision. The outgoing or post-collision velocities $`(\stackrel{}{v}^{},\stackrel{}{w}^{})`$ can be parametrized in terms of the incoming velocities $`(\stackrel{}{v},\stackrel{}{w})`$, and an impact (unit) vector $`\stackrel{}{n}`$, that is chosen on the surface of a unit sphere with a certain probability, proportional to the collision frequency $`\text{a}(\stackrel{}{g},\stackrel{}{n})=\text{a}(\mathrm{g},\widehat{\stackrel{}{g}}\stackrel{}{n})`$, that in general depends on the relative speed $`g=\left|\stackrel{}{v}\stackrel{}{w}\right|`$ of the colliding particles, and the angle between $`\stackrel{}{g}`$ and $`\stackrel{}{n}`$, where $`\widehat{\stackrel{}{a}}`$ denotes a unit vector.

In this article we consider the simplest case of modeling the inelastic collisions through a velocity-independent coefficient of restitution $`\alpha \left(0<\alpha <1\right)`$, where the component $`\mathrm{g}_{}\stackrel{}{g}\stackrel{}{n}`$ is not only reflected, as in elastic collisions, but also reduced in size by a factor $`\alpha `$, i.e.

$$\mathrm{g}_{}^{}{}_{}{}^{}=\alpha \mathrm{g}_{}.$$
(13)
The components $`\stackrel{}{g}_{}=\stackrel{}{g}\mathrm{g}_{}\stackrel{}{n}`$, orthogonal to $`\stackrel{}{n}`$, remain unchanged. More explicitly, with the help of momentum conservation we find for the post-collision velocities resulting from the direct collisions $`(\stackrel{}{v},\stackrel{}{w})(\stackrel{}{v}^{},\stackrel{}{w}^{})=(\stackrel{}{v}^{},\stackrel{}{w}^{})`$,

$$\begin{array}{cc}\hfill \stackrel{}{v}^{}=& \stackrel{}{v}\frac{1}{2}\left(1+\alpha \right)\stackrel{}{g}\stackrel{}{n}\stackrel{}{n}\hfill \\ \hfill \stackrel{}{w}^{}=& \stackrel{}{w}+\frac{1}{2}\left(1+\alpha \right)\stackrel{}{g}\stackrel{}{n}\stackrel{}{n}.\hfill \end{array}$$
(14)
The corresponding energy loss in such a collision is,

$$\mathrm{\Delta }E=\frac{1}{2}\left(v^2+w^2v^2w^2\right)=\frac{1}{4}\left(1\alpha ^2\right)\left(\stackrel{}{g}\stackrel{}{n}\right)^2,$$
(15)
where $`\left(1\alpha ^2\right)`$ measures the degree of inelasticity. The value $`\alpha =1`$ describes elastic collisions. In the special case of hard spheres the impact vector $`\stackrel{}{n}`$ is the unit vector along the line of centers of the colliding spheres at contact, as illustrated in Figure 1.

The system above can be described by an isotropic velocity distribution, $`F(v,t)=F(\left|\stackrel{}{v}\right|,t)`$, as long as one only considers isotropic initial distributions. Its time evolution is given by the Boltzmann equation,

$$\begin{array}{c}_tF(v,t)=I(v|F)d\stackrel{}{w}d\stackrel{}{v}^{}d\stackrel{}{w}^{}d\stackrel{}{n}[W(\stackrel{}{v},\stackrel{}{w}|\stackrel{}{v}^{},\stackrel{}{w}^{};\stackrel{}{n})\hfill \\ \hfill \times F(v^{},t)F(w^{},t)W(\stackrel{}{v}^{},\stackrel{}{w}^{}|\stackrel{}{v},\stackrel{}{w};\stackrel{}{n})F(v,t)F(w,t)].\end{array}$$
(16)
Here $`W(\stackrel{}{v}^{},\stackrel{}{w}^{}|\stackrel{}{v},\stackrel{}{w};\stackrel{}{n})`$ is the transition probability per unit time that the incoming pair state $`(\stackrel{}{v},\stackrel{}{w})`$ at impact vector $`\stackrel{}{n}`$ is scattered into the outgoing pair state $`(\stackrel{}{v}^{},\stackrel{}{w}^{})`$. The loss term is the sum over all parameters of the direct collisions $`(\stackrel{}{v},\stackrel{}{w})(\stackrel{}{v}^{},\stackrel{}{w}^{})`$ at fixed $`\stackrel{}{v}`$, including a sum over impact vectors $`\stackrel{}{n}`$. Similarly the gain term is the sum over all parameters of the restituting collisions, $`(\stackrel{}{v}^{},\stackrel{}{w}^{})(\stackrel{}{v},\stackrel{}{w})`$ at fixed $`\stackrel{}{v}`$. The transition probability for the scattering event $`(\stackrel{}{v},\stackrel{}{w})(\stackrel{}{v}^{},\stackrel{}{w}^{})`$, obeying the inelastic reflection law (13) and momentum conservation, is in general proportional to the collision frequency $`\text{a}(\stackrel{}{g},\stackrel{}{n})`$, and contains delta functions, selecting the allowed collisions,

$$W(\stackrel{}{v}^{},\stackrel{}{w}^{}|\stackrel{}{v},\stackrel{}{w};\stackrel{}{n})=\text{a}(\stackrel{}{g},\stackrel{}{n})\delta ^{\left(d\right)}\left(\stackrel{}{G}^{}\stackrel{}{G}\right)\delta ^{\left(d1\right)}\left(\stackrel{}{g}_{}^{}\stackrel{}{g}_{}\right)\delta \left(\mathrm{g}_{}^{}{}_{}{}^{}+\alpha \mathrm{g}_{}\right).$$
(17)
The faster the collision frequency increases at large impact velocities, the more rapidly the high energy tail of $`F(v,t)`$ relaxes relative to the bulk values of $`F(v,t)`$ with $`v`$ in the thermal range.

For conservative interactions the total energy is conserved, as well as total momentum and total number of particles. Moreover (17) shows that the transition probability is symmetric for $`\alpha =1`$,

$$W(\stackrel{}{v}^{},\stackrel{}{w}^{}|\stackrel{}{v},\stackrel{}{w};\stackrel{}{n})=W(\stackrel{}{v},\stackrel{}{w}|\stackrel{}{v}^{},\stackrel{}{w}^{};\stackrel{}{n}).$$
(18)
This means that the transition probabilities for elastic binary collisions satisfy the condition of detailed balance. Energy conservation and the detailed balance relation in combination with the $`H`$theorem guarantee that the Maxwellian velocity distribution is approached at large times.

Once the detailed balance relation is obeyed, it is trivial to prove the $`H`$theorem. Defining the $`H`$function or entropy $`S\left(t\right)=H\left(t\right)`$ as,

$$H\left(t\right)=𝑑\stackrel{}{v}F(\stackrel{}{v},t)\mathrm{ln}F(\stackrel{}{v},t),$$
(19)
we obtain its decay rate by applying $`𝑑\stackrel{}{v}\mathrm{ln}F(\stackrel{}{v},t)`$ to (16), then symmetrizing over $`\stackrel{}{v}\stackrel{}{w}`$ and $`\stackrel{}{v}^{}\stackrel{}{w}^{}`$, and subsequently over $`(\stackrel{}{v},\stackrel{}{w})(\stackrel{}{v}^{},\stackrel{}{w}^{})`$. The result is,

$$_tH\left(t\right)=\frac{1}{4}𝑑\stackrel{}{v}𝑑\stackrel{}{w}𝑑\stackrel{}{v}^{}𝑑\stackrel{}{w}^{}𝑑\stackrel{}{n}W\left[YX\right]\mathrm{ln}\left(X/Y\right)0,$$
(20)
where $`X=F\left(v\right)F\left(w\right)`$ and $`Y=F\left(v^{}\right)F\left(w^{}\right)`$ and the inequality follows from $`\left[YX\right]\mathrm{ln}\left(X/Y\right)0`$. The equality sign holds if and only if $`X=Y`$. This implies that $`H`$ decreases monotonically and becomes stationary only when $`F`$ approaches the Maxwellian.

The reason for reviewing these ’obvious’ properties in elastic systems is that several of them, such as energy conservation, the detailed balance relation, the $`H`$theorem, and the approach to a stationary Maxwellian velocity distribution no longer hold in dissipative systems $`\left(\alpha <1\right)`$, as we shall see.

### 3.2 Boltzmann equation in standard form

The standard form of the collision term $`I\left(v|F\right)`$ in the Boltzmann equation for elastic interactions ($`\alpha =1`$) contains the differential scattering cross-section $`\sigma (g,\chi )`$. It can be calculated from the pair potential $`V\left(r\right)`$ cercignani-book , and depends on the relative speed $`g=\left|\stackrel{}{g}\right|`$ and the scattering angle $`\chi =\pi 2\phi `$, where $`\chi `$ and $`\phi `$ are defined in Figure 1. The collision frequency is then given by $`\text{a}(\stackrel{}{g},\stackrel{}{n})=g\sigma (g,\chi )`$. For repulsive power law potentials $`V\left(r\right)r^s`$ the collision frequency is $`\text{a}(\stackrel{}{g},\stackrel{}{n})=g^\nu \text{a}\left(\widehat{\stackrel{}{g}}\stackrel{}{n}\right)`$ with $`\nu =12\left(d1\right)/s`$ in $`d`$dimensions ME-PhysRep . By definition Maxwell molecules have a collision frequency, which is independent of $`g`$, corresponding to $`\nu =0`$ or $`s=2\left(d1\right)`$. For hard spheres, $`s\mathrm{}`$ or $`\nu =1`$ with $`\text{a}(\stackrel{}{g},\stackrel{}{n})=\left|\stackrel{}{g}\stackrel{}{n}\right|\sigma ^{d1}=\left|\mathrm{g}_{}\right|\sigma ^{d1}`$, and $`\sigma `$ is the hard sphere diameter. Using these relations the collision term in (16) can be reduced to the classical Boltzmann equation for elastic energy-conserving collisions. For inelastic hard spheres the collision rate is again given by $`\text{a}(\stackrel{}{g},\stackrel{}{n})=\left|\mathrm{g}_{}\right|\sigma ^{d1}`$. In general, a given positive function $`\text{a}(\stackrel{}{g},\stackrel{}{n})`$ defines a stochastic scattering model, and in particular the choice, $`\text{a}(\stackrel{}{g},\stackrel{}{n})=\text{a}\left(\widehat{\stackrel{}{g}}\stackrel{}{n}\right)`$, defines inelastic Maxwell models BN-PK-00 ; Bobyl-00 ; CCG-00 . As a simple realization of a dissipative scattering model we consider the IRS-models with a collision frequency $`\text{a}(\stackrel{}{g},\stackrel{}{n})=\text{a}_0\left|\stackrel{}{g}\stackrel{}{n}\right|^\nu `$ with $`\nu 0`$. This class includes inelastic hard spheres $`\left(\nu =1\right)`$ and inelastic Maxwell models $`\left(\nu =0\right)`$. In the remainder of this article we restrict ourselves to this class of models.

To reduce the Boltzmann equation to its standard form we use the transition probabilities for dissipative interactions (17) with $`\text{a}(\stackrel{}{g},\stackrel{}{n})=\text{a}_0\left|\mathrm{g}_{}\right|^\nu `$. Consider first the loss term in (16), insert the relation $`d\stackrel{}{v}^{}d\stackrel{}{w}^{}=d\stackrel{}{g}_{}^{}d\mathrm{g}_{}^{}{}_{}{}^{}d\stackrel{}{G}^{}`$, and carry out the integrations over the delta functions. The result is,

$$\left(_tF\left(v\right)\right)_{\text{loss}}=𝑑\stackrel{}{w}𝑑\stackrel{}{n}\text{a}_0\left|\mathrm{g}_{}\right|^\nu F\left(w\right)F\left(v\right).$$
(21)
The gain term involves the transition probability for the restituting collisions $`(\stackrel{}{v}^{},\stackrel{}{w}^{})(\stackrel{}{v},\stackrel{}{w})`$, obtained from (17) by interchanging primed and unprimed velocities,

$$\begin{array}{cc}\hfill \left(_tF\left(v\right)\right)_{\text{gain}}=& 𝑑\stackrel{}{w}𝑑\stackrel{}{g}_{}^{}𝑑\mathrm{g}_{}^{}{}_{}{}^{}𝑑\stackrel{}{G}^{}𝑑\stackrel{}{n}\text{a}_0\left|\mathrm{g}_{}^{}{}_{}{}^{}\right|^\nu \delta ^{\left(d\right)}\left(\stackrel{}{G}\stackrel{}{G}^{}\right)\hfill \\ & \times \delta ^{\left(d1\right)}\left(\stackrel{}{g}_{}\stackrel{}{g}_{}^{}\right)\frac{1}{\alpha }\delta \left(\mathrm{g}_{}^{}{}_{}{}^{}+\frac{1}{\alpha }\mathrm{g}_{}\right)F\left(w^{}\right)F\left(v^{}\right)\hfill \\ \hfill =& 𝑑\stackrel{}{w}𝑑\stackrel{}{n}\left(1/\alpha \right)\text{a}_0\left|\mathrm{g}_{}/\alpha \right|^\nu F\left(w^{}\right)F\left(v^{}\right).\hfill \end{array}$$
(22)
In the second integral we have carried out the integrations over the primed velocities, and used the following relations for the restituting velocities,

$$\begin{array}{cc}\hfill \stackrel{}{v}^{}=& \stackrel{}{G}^{}+\frac{1}{2}\stackrel{}{g}_{}^{}+\frac{1}{2}\mathrm{g}_{}^{}{}_{}{}^{}\stackrel{}{n}=\stackrel{}{G}+\frac{1}{2}\stackrel{}{g}_{}\frac{1}{2\alpha }\mathrm{g}_{}\stackrel{}{n}\hfill \\ \hfill =& \stackrel{}{v}\frac{1}{2}\left(1+\frac{1}{\alpha }\right)\stackrel{}{g}\stackrel{}{n}\stackrel{}{n}\stackrel{}{v}^{}\hfill \\ \hfill \stackrel{}{w}^{}=& \stackrel{}{w}+\frac{1}{2}\left(1+\frac{1}{\alpha }\right)\stackrel{}{g}\stackrel{}{n}\stackrel{}{n}\stackrel{}{w}^{}.\hfill \end{array}$$
(23)
In the first equality $`\stackrel{}{v}^{}`$ has been expressed in center of mass and relative velocities. In the second equality we have used the inelastic collision law (13) and conservation of total momentum, and the very last equality defines the restituting velocities, $`(\stackrel{}{v}^{},\stackrel{}{w}^{})`$. They are the incoming velocities that result in the scattered velocities $`(\stackrel{}{v},\stackrel{}{w})`$, described by the inverse of the transformation (14).

The space-homogeneous Boltzmann equation in its standard form is then obtained by combining (21) and (22) with $`\text{a}(\stackrel{}{g},\stackrel{}{n})=\left|\stackrel{}{g}\stackrel{}{n}\right|^\nu `$ to yield,

$$\begin{array}{cc}\hfill _tF\left(v\right)=& I\left(v|F\right)\hfill \\ \hfill I\left(v|F\right)=& _\stackrel{}{n}𝑑\stackrel{}{w}\left|\stackrel{}{g}\stackrel{}{n}\right|^\nu \left[\frac{1}{\alpha ^{\nu +1}}F\left(v^{}\right)F\left(w^{}\right)F\left(v\right)F\left(w\right)\right],\hfill \end{array}$$
(24)
where $`_\stackrel{}{n}\left(\mathrm{}\right)=\left(1/\mathrm{\Omega }_d\right)𝑑\stackrel{}{n}\left(\mathrm{}\right)`$ is an average over a $`d`$dimensional unit sphere, and we have absorbed constant factors in the time scale. Here $`\nu =1`$ corresponds to inelastic hard spheres and $`\nu =0`$ to inelastic Maxwell models. Velocities and time have been dimensionalized in terms of the width and the mean free time of the initial distribution. Moreover, the Boltzmann equation obeys conservation of particle number and total momentum, but the average kinetic energy or granular temperature, $`Tv^2`$, decreases in time on account of the dissipative collisions, i.e.

$$𝑑\stackrel{}{v}(1,\stackrel{}{v},v^2)F(v,t)=(1,\mathrm{\hspace{0.25em}0},\frac{1}{2}dv_0^2\left(t\right)),$$
(25)
where $`v_0\left(t\right)`$ is the r.m.s. velocity. The inelastic scattering models with collision frequency $`\text{a}g^\nu \left(\nu >0\right)`$, are the inelastic analogs of the deterministic models with repulsive power law potentials, $`V\left(r\right)r^s`$ with $`2\left(d1\right)<s<\mathrm{}`$, introduced at the start of this section. The inelastic case with $`\nu =2`$ corresponds to an exactly solvable stochastic scattering model with energy conservation, known as Very Hard Particle model ME-PhysRep .

The most basic and most frequently used model for dissipative systems with short range hard core repulsion is the Enskog-Boltzmann equation for inelastic hard spheres in $`d`$dimensions campbell , which simplifies in the spatially homogeneous case to (24) with $`\nu =1`$. Recently, inelastic Maxwell models have been studied extensively. Ben-Naim and Krapivky BN-PK-00  introduced the one-dimensional version of (24) with $`\left(d=1,\nu =0\right)`$, and Bobylev et al Bobyl-00  have introduced a three-dimensional Maxwell model with $`\text{a}(\stackrel{}{g},\stackrel{}{n})=\text{a}_0\left|\widehat{\stackrel{}{g}}\stackrel{}{n}\right|`$. The Maxwell models in (24) for general $`d`$ were first considered in BN+PK-6-11 ; ME+RB-EPL .

### 3.3 Cooling and driven systems

Homogeneous cooling and scaling:
An inelastic fluid without energy input will cool down due to the collisional dissipation in (15). In experimental studies of granular fluids energy has to be supplied at a constant rate to keep the system in a non-equilibrium steady state, while in analytic, numerical and simulation work freely cooling systems can be studied directly. Without energy input the velocity distribution $`F(v,t)`$ will approach a Dirac delta function $`\delta \left(\stackrel{}{v}\right)`$ as $`t\mathrm{}`$, and all moments approach zero, including the width $`v_0\left(t\right)`$.

However, an interesting structure is revealed when velocities, $`\stackrel{}{c}=\stackrel{}{v}/v_0\left(t\right)`$, are measured in units of the instantaneous width $`v_0\left(t\right)`$, and the long time limit is taken while keeping $`\stackrel{}{c}`$ constant, the so-called scaling limit. Monte Carlo simulations Rome1  of the Boltzmann equation suggest that in this limit the rescaled velocity distribution of the homogeneous cooling state can be collapsed on a scaling form or similarity solution $`f\left(c\right)`$. These observations seem to indicate that the long time behavior of $`F(v,t)`$ in freely cooling systems approaches a simple, and to some extent universal, form $`f\left(c\right)`$, which is the same for different initial distributions. Such scaling or similarity solutions have the structure,

$$F(v,t)=\left(v_0\left(t\right)\right)^df\left(v/v_0\left(t\right)\right),$$
(26)
where $`\stackrel{}{c}=\stackrel{}{v}/v_0\left(t\right)`$ is the scaling argument. Then $`f\left(c\right)`$ satisfies the normalizations,

$$𝑑\stackrel{}{c}f\left(c\right)=1;𝑑\stackrel{}{c}c^2f\left(c\right)=\frac{1}{2}d,$$
(27)
on account of (25). Substitution of the scaling ansatz (26) in the Boltzmann equation (24) then leads to a separation of variables, i.e.

$$\begin{array}{cc}\hfill I\left(c|f\right)=& \gamma \left(df+cdf/dc\right)=\gamma \stackrel{}{}_\stackrel{}{c}\stackrel{}{c}f\hfill \\ \hfill \dot{v}_0=& \gamma v_0^{\nu +1},\hfill \end{array}$$
(28)
where $`\gamma `$ is a separation constant.

The question is then: Can one determine physically acceptable solutions? Here a short history. The first kinetic model with dissipative interactions, that has been solved exactly for the case of free cooling, and exhibits a heavily overpopulated high energy tail, $`f\left(c\right)1/c^{d+a}`$, is the inelastic BGK model, discussed in Section 2. Asymptotic solutions $`\left(c1\right)`$ of the scaling equation for inelastic hard spheres esipov ; TvN+ME-granmat  have predicted the existence of exponential high energy tails, $`f\left(c\right)\mathrm{exp}\left[\beta c\right]`$. We talk about tails, over-populations or stretched Gaussians in $`f\left(c\right)`$, when the ratio of $`f\left(c\right)`$ and a Gaussian is an increasing function of $`c`$ at $`c1`$. The predictions about high energy tails were later confirmed in great detail by Monte Carlo simulations of the long time solutions of the nonlinear Boltzmann equation for inelastic hard spheres brey ; MS00 , as well as in further analytic work BBRTvW01 .

Similarity solutions for freely cooling inelastic Maxwell models BN-PK-00 ; Bobyl-00  where first studied in terms of the scaling variable, $`\left(1\alpha ^2\right)t`$, relating large times and small inelasticities Bobyl-00 , but the solutions obtained turned out to be unphysical, i.e. non-positive. The first exact positive similarity solution was found by Baldassarri et al. Rome1  for the one-dimensional Maxwell model. It shows a surprisingly strong high energy tail of algebraic type, $`1/c^4`$. Using Monte Carlo simulations of the Boltzmann equation, these authors also showed that rather general initial distributions approach towards this scaling form for long times. Algebraic tails, $`f\left(c\right)1/c^{d+a}`$, for Maxwell models in higher dimensions have been obtained analytically in BN+PK-6-11 ; BN+PK-02 ; ME+RB-EPL ; ME+RB-rapid ; ME+RB-fest , and the actual approach in time towards such scaling forms has been studied analytically in ME+RB-fest ; AB+CC-proof  for general classes of initial distributions.

Driving and steady states:
Heating may be described by applying an external stochastic force to the particles in the system, or by connecting the system to a thermostat, which may be modeled by a frictional force. For example, the friction force $`\gamma \stackrel{}{v}`$ – here with a negative friction coefficient $`\left(\gamma \right)`$ – is called a Gaussian thermostat. Complex fluids (e.g. granular) subject to such forces can be described by the microscopic equations of motion for the particles, $`\dot{\stackrel{}{r}}_i=\stackrel{}{v}_i`$, and $`\dot{\stackrel{}{v}}_i=\stackrel{}{a}_i+\stackrel{~}{\stackrel{}{\xi }}_i`$ $`\left(i=1,2,\mathrm{}\right)`$, where $`\stackrel{}{a}_i`$ and $`\stackrel{~}{\stackrel{}{\xi }}_i`$ are respectively the systematic and random forces per unit mass. Here $`\stackrel{}{a}_i`$ contains frictional forces, which may depend on velocity, and in the present case $`\stackrel{~}{\stackrel{}{\xi }}_i`$ represents external noise (modeling energy input), which is taken to be Gaussian white noise with zero mean, and variance,

$$\overline{\stackrel{~}{\xi }_{i,\alpha }\left(t\right)\stackrel{~}{\xi }_{j,\beta }\left(t^{}\right)}=2D\delta _{ij}\delta _{\alpha \beta }\delta \left(tt^{}\right),$$
(29)
where $`\alpha ,\beta `$ denote Cartesian components, and $`D`$ is the noise strength. The Boltzmann equation for system driven in this manner takes the form,

$$\begin{array}{cc}& _tF\left(\stackrel{}{v}\right)+F\left(\stackrel{}{v}\right)=I\left(v|F\right)\hfill \\ & =\stackrel{}{}_\stackrel{}{v}\stackrel{}{a}D_\stackrel{}{v}^2.\hfill \end{array}$$
(30)
A detailed derivation on how to include frictional and stochastic forces in kinetic equations can be found, for instance, in ME-TvN-badhonnef ; WM96 ; TvN+ME-granmat .

When energy is supplied at a constant rate, driven systems can reach a NESS. Again there is the question about universality of these asymptotic states as $`t\mathrm{}`$. Does the NESS depend on the inelasticity, on the type of thermostat, and on the initial distribution? The basic idea to show universality is always essentially the same. Rescale the velocity distribution $`F(v,\mathrm{})`$ by measuring velocities in terms of their typical size, i.e. the width $`v_0\left(\mathrm{}\right)`$,

$$F(v,\mathrm{})=\left(v_0\left(\mathrm{}\right)\right)^df\left(v/v_0\left(\mathrm{}\right)\right),$$
(31)
and analyze the scaling form $`f\left(c\right)`$. The scaling equation for the NESS function $`f\left(c\right)`$ will be analyzed in Sections 4.2 and 4.3. Here we make some comments about what is known. The scaling solutions for the NESS show again overpopulated tails in the form of stretched Gaussians, albeit with a larger stretching exponent $`b`$ than in free cooling. The scaling solution for inelastic hard spheres driven by white noise was first analyzed in TvN+ME-granmat , predicting $`f\left(c\right)\mathrm{exp}\left[\beta c^b\right]`$ with $`b=3/2`$. The experiments of Rouyer and Menon GM-exp ; SciAm-01 , and of Aranson and Olafsen GM-exp  seem to confirm the stretched Gaussian behavior with $`b=3/2`$. Maxwell models, driven by white noise, were first studied in BN-PK-00 ; CCG-00 , but did not give any predictions about overpopulated high energy tails. The high energy tails for Maxwell models, driven by white noise, were predicted in ME+RB-rapid  to have exponential tails, $`f\left(c\right)\mathrm{exp}\left[\beta c\right]`$. A number of more detailed analytical and numerical studies about this driven Maxwell model in one dimension Nienhuis ; Droz ; AS+ME-condmat ; BN-PK-Springer  have appeared as well. Moreover, as will be shown in Section 4.2, the integral equation for the scaling form $`f\left(c\right)`$ in free cooling is identical to the integral equation for the NESS distribution for a special thermostat, provided $`F(v,\mathrm{})`$ is also rescaled to the same width as in (27). How to determine the scaling form $`f\left(c\right)`$ in free cooling and NESS will be described in subsequent sections.

### 3.4 Qualitative analysis

In order to illustrate the rich behavior of the inelastic systems, we start with the Boltzmann equation for the one-dimensional inelastic Maxwell model (24) ($`d=1,\nu =0`$), where $`F(v,t)`$ satisfies,

$$_tF\left(v\right)D_v^2F\left(v\right)=I\left(v|F\right),$$
(32)
and the collision term has the form,

$$\begin{array}{cc}\hfill I\left(v|F\right)=& 𝑑w\left[\frac{1}{\alpha }F\left(v^{}\right)F\left(w^{}\right)F\left(v\right)F\left(w\right)\right]\hfill \\ \hfill =& F\left(v\right)+\frac{1}{p}𝑑uF\left(u\right)F\left(\frac{vqu}{p}\right).\hfill \end{array}$$
(33)
All velocity integrations extend over the interval $`(\mathrm{},+\mathrm{})`$. Here the outgoing velocities $`(v,w)`$, and the incoming ones $`(v^{},w^{})`$ are according to (23) related by,

$$v=qv^{}+pw^{};w=pv^{}+qw^{}$$
(34)
with $`p=1q=\frac{1}{2}\left(1+\alpha \right)`$. By changing integration variables $`wv^{}=u`$ with $`dw=\left(\alpha /p\right)du`$, and using the relation $`w^{}=\left(vqu\right)/p`$ one obtains the second equality in (33). The normalization of mass and mean square velocity $`v^2=\frac{1}{2}v_0^2`$ are given by (25) for $`d=1`$. The temperature balance equation is obtained from (32) as,

$$_tv_0^2=4D2pqv_0^2,$$
(35)
and describes the approach to the non-equilibrium steady state (NESS) with a width $`v_0^2\left(\mathrm{}\right)=2D/pq`$, where the heating through random forces, $`D`$, is balanced by the collisional losses.

To understand the physical processes involved we first discuss in a qualitative way the relevant limiting cases. Without the heating term $`\left(D=0\right)`$, equation (32) reduces to the freely cooling inelastic Maxwell model.

If one takes in addition the elastic limit $`(\alpha 1`$ or $`q0)`$, the collision laws reduce in the one-dimensional case to $`v^{}=w,w^{}=v`$, i.e. an exchange of particle labels; the collision term vanishes identically; every $`F(v,t)=F\left(v\right)`$ is a solution; there is no randomization or relaxation of the velocity distribution through collisions, and the model becomes trivial at the Boltzmann level of description, whereas the distribution function in the presence of infinitesimal dissipation $`\left(\alpha 1\right)`$ approaches a Maxwellian.

If we turn on the noise $`\left(D0\right)`$ at vanishing dissipation $`\left(q=0\right)`$, the collision term in (32) vanishes, and the granular temperature follows from (35) as $`v_0^2\left(t\right)=v_0^2\left(0\right)+4Dt`$, and increases linearly with time. With stochastic heating and dissipation (even in infinitesimal amounts) the system reaches a NESS.

### 3.5 Comments

Spatial dependence:
When spatial dependence of the distribution function $`F(\stackrel{}{r},\stackrel{}{v};t)`$ is relevant, the collision term in the Enskog-Boltzmann equation must be slightly modified, namely the angular integration $`_\stackrel{}{n}`$ over the full solid angle should be replaced by $`2_\stackrel{}{n}\theta \left(\stackrel{}{g}\stackrel{}{n}\right)`$, where the unit step function $`\theta \left(x\right)`$ restricts the $`\stackrel{}{n}`$integration to the pre-collision hemisphere with $`\stackrel{}{g}\stackrel{}{n}<0`$ (see Figure 1). In the spatially uniform case both representations are identical because the restituting velocities (23) are even functions of $`\stackrel{}{n}`$.

Origin of Maxwell models:
In several papers Bobyl-00 ; BN+PK-6-11 ; ME+RB-fest ; Nienhuis  inelastic Maxwell models have been introduced as more or less ad hoc mathematical simplifications of the nonlinear collision term for inelastic hard spheres. This has been done by replacing the relative velocity $`g`$ in the hard sphere collision frequency $`\text{a}(\stackrel{}{g},\stackrel{}{n})=\left|\stackrel{}{g}\stackrel{}{n}\right|`$ by its mean value, $`g\text{a}\left(\stackrel{}{n}\right)v_0\left(t\right)\text{a}\left(\stackrel{}{n}\right)`$, where $`v_0\left(t\right)`$ is the root mean square velocity.

This procedure guarantees that the homogeneous cooling law for inelastic Maxwell model constructed in this way, is identical to the one for inelastic hard spheres, and given by Haff’s law haff , $`T\left(t\right)t^2`$. The construction of inelastic Maxwell and IRS models, followed in the present article, is more in the spirit of BN-PK-00 , i.e. by defining the collision term through transition probabilities for the scattering process $`(\stackrel{}{v},\stackrel{}{w})(\stackrel{}{v}^{},\stackrel{}{w}^{})`$ with the proper constraints.

Particles with pseudo-power law repulsion:
A system of $`N`$ inelastic hard spheres is a model with microscopic particles, characterized by positions and velocities, and interacting via well-defined force laws, that can be studied by means of MD simulations. However microscopic particles with dissipative inter-particle forces, such as inelastic Maxwell molecules and IRS-models, are stochastic models, only defined in $`N`$particle velocity space. Molecular Dynamics simulations can not be performed for such models.

Violation of $`H`$theorem:
As the detailed balance symmetry and the $`H`$theorem are lacking in dissipative interaction models, there is no guarantee that the entropy $`S\left(t\right)=H\left(t\right)`$ is non-decreasing. In fact the solutions of the Boltzmann equation for such dissipative interaction models approach ME+RB-fest  for long times to a scaling form, defined in (26). By inserting such solutions into (19) and anticipating the cooling law (3) in the next section, one easily verifies that the entropy in the scaling state keeps decreasing as $`t`$ becomes large, i.e.

$$S\left(t\right)=H\left(t\right)\left(d/\nu \right)\mathrm{ln}t.$$
(36)
This result is typical for pattern forming mechanisms in configuration space, where spatial order or correlations are building up, as well as in chaos theory, where the rate of irreversible entropy production is negative on an attractor Dorfman . Moreover, there is no fundamental objection against decreasing entropies in an open system in contact with a reservoir, which is here the energy sink formed by the dissipative collisions. The dynamics in $`N`$particle velocity space corresponds to a contracting flow $`d\stackrel{}{v}^{}d\stackrel{}{w}^{}=\alpha d\stackrel{}{v}d\stackrel{}{w}`$ where $`\alpha <1`$, where the probability is contracting onto an attractor. This is a well known phenomenon in chaos theory.

Modeling dissipation:
There are many ways to model inelastic collisions that dissipate the relative kinetic energy of colliding particles; e.g. by using Hertz’ contact law Brill-TP-badhonnef ; chen-doolen-phys Lett , visco-elastic media Hutter , or coefficients of normal and tangential restitution campbell ; goldstein-shapiro . The coefficients of restitution may also depend on the relative speed of the colliding particles Brill-TP-badhonnef .

For the IRS-models, studied in this article, the scattering laws (scattering angle, collisional energy loss) are independent of $`\nu `$, i.e. they are the same for inelastic hard spheres, inelastic Maxwell models and for general IRS-models. Only the collision frequency $`\text{a}(𝐠,𝐧)`$ of the IRS-models has the same energy dependence as the collision frequency of elastic particles interacting through repulsive power law potentials.

## 4 Analysis of inelastic scattering models

### 4.1 Homogeneous cooling laws

We start with the freely evolving case and compute the cooling rate, defined through, $`_tv_0^2=\zeta _\nu \left(t\right)v_0^2`$. This can be done by applying $`\left(𝑑\stackrel{}{v}v^2\right)`$ to the Boltzmann equation (24), changing integration variables $`(\stackrel{}{v},\stackrel{}{w})(\stackrel{}{v}^{},\stackrel{}{w}^{})`$ in the gain term, and using the relations $`d\stackrel{}{v}d\stackrel{}{w}=\alpha d\stackrel{}{v}^{}d\stackrel{}{w}^{}`$ together with the energy loss $`\mathrm{\Delta }E`$ per collision in (15). In general the above cooling equation is a formal identity, and does not provide a closed equation for $`v_0(t)`$. However, if $`F(v,t)`$ rapidly relaxes to a scaling form, then the subsequent time evolution of $`v_0^2=(2/d)v^2`$ is described by a closed equation for $`v_0(t)`$, as we will see, i.e.

$$\begin{array}{c}_tv_0^2=\frac{2}{d}𝑑\stackrel{}{v}v^2I\left(v|F\right)=\frac{2}{d}_\stackrel{}{n}𝑑\stackrel{}{v}𝑑\stackrel{}{w}\left|\stackrel{}{g}\stackrel{}{n}\right|^\nu \mathrm{\Delta }EF(v,t)F(w,t)\hfill \\ \hfill =\frac{1}{2d}\left(1\alpha ^2\right)v_0^{\nu +2}_\stackrel{}{n}𝑑\stackrel{}{c}𝑑\stackrel{}{c}_1\left|\left(\stackrel{}{c}\stackrel{}{c}_1\right)\stackrel{}{n}\right|^{\nu +2}f\left(c\right)f\left(c_1\right),\end{array}$$
(37)
from which the cooling rate can be identified as,

$$\zeta _\nu \left(t\right)_tv_0^2/v_0^2=2\gamma _0\kappa _{\nu +2}v_0^\nu \left(t\right).$$
(38)
Here the coefficient $`\gamma _0=\frac{1}{4d}\left(1\alpha ^2\right)`$ measures the inelasticity. Moreover, the constant $`\kappa _\nu `$ is defined as,

$$\kappa _\nu =_\stackrel{}{n}𝑑\stackrel{}{c}𝑑\stackrel{}{c}_1\left|\left(\stackrel{}{c}\stackrel{}{c}_1\right)\stackrel{}{n}\right|^\nu f\left(c\right)f\left(c_1\right)=\beta _\nu 𝑑\stackrel{}{c}𝑑\stackrel{}{c}_1|\stackrel{}{c}\stackrel{}{c}_1|^\nu f\left(c\right)f\left(c_1\right),$$
(39)
and the $`\left(d1\right)`$dimensional angular integral $`\beta _\nu `$ is evaluated as,

$$\beta _\nu =_\stackrel{}{n}\left|\widehat{\stackrel{}{a}}\stackrel{}{n}\right|^\nu =\frac{_0^{\pi /2}𝑑\theta \left(\mathrm{sin}\theta \right)^{d2}\left(\mathrm{cos}\theta \right)^\nu }{_0^{\pi /2}𝑑\theta \left(\mathrm{sin}\theta \right)^{d2}}=\frac{\mathrm{\Gamma }\left(\frac{\nu +1}{2}\right)\mathrm{\Gamma }\left(\frac{d}{2}\right)}{\mathrm{\Gamma }\left(\frac{\nu +d}{2}\right)\mathrm{\Gamma }\left(\frac{1}{2}\right)}.$$
(40)
In one-dimensional systems $`\beta _\nu =1`$ for all $`\nu `$, and in higher dimensions for the $`\nu `$-values of interest here $`\beta _\nu `$ is convergent ($`\nu >0`$). As a consequence of the scaling ansatz the coefficient $`\kappa _\nu `$ is a time independent constant, that does depend on the unknown scaling form. For Maxwell models $`\left(\nu =0\right)`$ the cooling rate $`\zeta _0`$ can be calculated explicitly from (40), and the constant $`\kappa _2`$ is $`\kappa _2=\left(1/d\right)\left|\stackrel{}{c}\stackrel{}{c}_1\right|^2=1`$ on account of (27). The result is,

$$\zeta _0_tv_0^2/v_0^2=2\gamma _0.$$
(41)
The r.m.s. velocity decays then as $`v_0\left(t\right)=v_0\left(0\right)\mathrm{exp}\left[\gamma _0t\right]`$.

The most important time scale in kinetic theory is the mean free time $`t_{mf}`$, which equals the inverse of the mean collision rate $`\omega _\nu \left(t\right)`$, defined as the average of the collision frequency $`\text{a}(\stackrel{}{g},\stackrel{}{n})=\left|\left(\stackrel{}{c}\stackrel{}{c}_1\right)\stackrel{}{n}\right|^\nu `$ over the velocities of the colliding pair in the scaling state, i.e.

$$\omega _\nu \left(t\right)𝑑\stackrel{}{v}I_{\text{loss}}\left(v|F\right)=\kappa _\nu v_0^\nu \left(t\right).$$
(42)
Comparison of (38) and (42) shows that both frequencies are related as,

$$\zeta _\nu \left(t\right)=2\gamma _0\frac{\kappa _{\nu +2}}{\kappa _\nu }\omega _\nu \left(t\right)2\gamma _\nu \omega _\nu \left(t\right),$$
(43)
where $`\gamma _\nu \kappa _\nu =\gamma _0\kappa _{\nu +2}`$ are time independent constants. Using (42) and (43) the equation for the r.m.s. velocity becomes, $`_tv_0=\gamma _\nu \kappa _\nu v_0^{\nu +1}`$, yielding a solution, identical to (3) with $`\gamma `$ replaced by $`\gamma _\nu `$. Consequently the granular temperature at large times decays as $`T=v_0^2t^{2/\nu }`$. For $`\nu =1`$ one recovers the well known law of Haff haff , describing the long time decay of the granular temperature in the homogeneous cooling state of inelastic hard spheres. A homogeneous cooling law similar to (3) with $`\nu =6/5`$ has been derived in Brill-TP-badhonnef , not on the basis of the Boltzmann equation (24) with an energy dependent collision rate $`g^{5/6}`$, but by using the Boltzmann equation for inelastic hard spheres $`\left(\nu =1\right)`$ with an energy dependent coefficient of restitution $`\alpha \left(g\right)=1Ag^{2/5}`$. For Maxwell models $`\left(\nu 0\right)`$ the r.m.s. velocity in (3) reduces to $`v_0\left(t\right)=v_0\left(0\right)\mathrm{exp}\left[\gamma _0t\right]`$, in agreement with (41).

The homogeneous cooling law for the general class of IRS-models discussed here, can be cast into a universal form by changing to a new time variable, the collision counter or internal time $`\tau `$ of a particle, which represents the total number of collisions that a particle has suffered in the (external) time $`t`$. It is defined through the mean collision frequency,

$$d\tau =\omega _\nu \left(t\right)dt,$$
(44)
where $`\omega _\nu (t)v_0^\nu (t)`$. After inserting (3) in this equation the differential equation can be solved to yield,

$$\nu \gamma _\nu \tau =\mathrm{ln}\left[1+\nu \gamma _\nu \omega _\nu \left(0\right)t\right],$$
(45)
valid for all IRS-models. Combination of (45) and (3) shows that the r.m.s. velocity follows the universal homogeneous cooling law,

$$v_0\left(\tau \right)=v_0\left(0\right)\mathrm{exp}\left[\gamma _\nu \tau \right],$$
(46)
where $`\gamma _\nu `$ is according to (43) and (39) proportional to the inelasticity $`\gamma _0`$, and depends on $`\nu `$ through the collision integrals $`\kappa _\nu `$. We further observe that the relations (45) and (46) also holds for the inelastic BGK-model in Section 2 with $`\gamma _\nu `$ replaced by $`\gamma =\frac{1}{2}\left(1\alpha ^2\right)`$ in (6).

So far we have been dealing with freely cooling systems. Next we address the balance equation for the granular temperature in driven cases, where the external input of energy counterbalances the collisional cooling, and may lead to a NESS. We proceed in the same manner as for the free case, and apply $`\left(𝑑\stackrel{}{v}v^2\right)`$ to the Boltzmann equation in (30) with the result,

$$_tv^2=\zeta _\nu \left(t\right)v^2+2\stackrel{}{v}\stackrel{}{a}+2dD,$$
(47)
where the first term, $`𝑑\stackrel{}{v}v^2I\left(v|F\right)=\zeta _\nu \left(t\right)v^2`$, is obtained from (37), (38) and (25). The next two terms are obtained from the driving term in (30), i.e. $`𝑑\stackrel{}{v}v^2F\left(v\right)`$, by performing partial integrations.

The most common ways of driving dissipative fluids TvN+ME-granmat ; MS00 ; BN-PK-00 ; CCG-00  is by Gaussian white noise (WN) $`\left(\stackrel{}{a}=0;D0\right)`$, or by a Gaussian thermostat (GT) $`\left(\stackrel{}{a}=\gamma \stackrel{}{v};D=0\right)`$, yielding for the balance equations,

$$_tv_0^2\left(t\right)=\{\begin{array}{cc}\left(2\gamma \zeta _\nu \left(t\right)\right)v_0^2\left(t\right)\hfill & \text{(GT)}\hfill \\ 4D\zeta _\nu \left(t\right)v_0^2\left(t\right)\hfill & \text{(WN).}\hfill \end{array}$$
(48)
Here the collisional loss, $`\zeta _\nu v_0^2`$, is counterbalanced by the heat, $`2\gamma v_0^2`$, generated by the negative friction of the Gaussian thermostat, or by the heat, $`4D`$, generated by randomly kicking the particles.

As $`t\mathrm{}`$ the granular temperature will reach a NESS with $`T\left(\mathrm{}\right)=v_0^2\left(\mathrm{}\right)`$, where the r.m.s. velocity follows from (48) using (38),

$$v_0\left(\mathrm{}\right)=\{\begin{array}{cc}\left(\gamma /\gamma _0\kappa _{\nu +2}\right)^{1/\nu }\hfill & \text{(GT)}\hfill \\ \left(2D/\gamma _0\kappa _{\nu +2}\right)^{1/\left(\nu +2\right)}\hfill & \text{(WN)}.\hfill \end{array}$$
(49)
With reference to the discussion in Section 2.3 we note that the GT-fixed point solution $`v_0\left(\mathrm{}\right)`$ is attracting for $`\nu >0`$, leading to a stable NESS, and unstable for $`\nu <0`$ with $`v_0\left(t\right)`$ vanishing at $`t\mathrm{}`$ if $`v_0\left(0\right)<v_0\left(\mathrm{}\right)`$, and diverging if $`v_0\left(0\right)>v_0\left(\mathrm{}\right)`$. The case $`\nu =0`$ is marginally stable. For the WN-fixed point similar observations apply with the stability threshold $`\nu =0`$ replaced by $`\nu =2`$.

### 4.2 Scaling and non-equilibrium steady states

Free cases:
To investigate the existence of scaling solutions of the Boltzmann equation for freely cooling inelastic systems, we substitute the scaling ansatz (26) into (24), to obtain an integral equation for the scaling form $`f\left(c\right)`$. With the help of (44), (46) and (28) the left hand side of (24) becomes,

$$\begin{array}{c}\text{l.h.s.}=\dot{v}_0v_0^{d1}\left(\frac{d\tau }{dt}\right)\left\{df\left(c\right)+c\frac{d}{dc}f\left(c\right)\right\}\hfill \\ \hfill =\gamma _\nu \omega _\nu \left(\tau \right)v_0^d\left(\tau \right)_\stackrel{}{c}\stackrel{}{c}f\left(c\right)=\gamma _0\kappa _{\nu +2}_\stackrel{}{c}\stackrel{}{c}f\left(c\right)v_0^{\nu d}\left(\tau \right).\end{array}$$
(50)
In the last equality we have used the relation $`\gamma _\nu \kappa _\nu =\gamma _0\kappa _{\nu +2}`$, implied by (43). The resulting integral equation for $`f\left(c\right)`$ becomes,

$$I\left(c|f\right)=\gamma _0\kappa _{\nu +2}_\stackrel{}{c}\stackrel{}{c}f=\frac{1}{d}\varpi _2_\stackrel{}{c}\stackrel{}{c}f.$$
(51)
With the help of (37) and (38) the second moment of the collision term, $`\varpi _2`$, can be expressed as,

$$\varpi _2𝑑\stackrel{}{c}c^2I\left(c|f\right)=d\gamma _0\kappa _{\nu +2}.$$
(52)
Driven cases:
We consider equation (30) for the NESS distribution in its rescaled form (31) with $`v_0=v_0\left(\mathrm{}\right)`$, driven by a Gaussian thermostat (GT), $`\{\stackrel{}{a}=\gamma \stackrel{}{v};D=0\}`$ or by white noise (WN), $`\{\stackrel{}{a}=0;D0\}`$. The scaling equations for $`f\left(c\right)`$ take the form,

$$I\left(c|f\right)=\{\begin{array}{cc}\frac{\gamma }{v_0^\nu }_\stackrel{}{c}\stackrel{}{c}f=\frac{1}{d}\varpi _2_\stackrel{}{c}\stackrel{}{c}f\hfill & \text{(GT)}\hfill \\ \frac{D}{v_0^{\nu +2}}_\stackrel{}{c}^2f=\frac{1}{2d}\varpi _2_\stackrel{}{c}^2f\hfill & \text{(WN)}.\hfill \end{array}$$
(53)
The first equality in (GT) and (WN) suggests that $`f\left(c\right)`$ depends explicitly on $`\gamma `$ or $`D`$. This is however not the case. The stationarity relation (49) combined with (52) shows in fact that the following expressions,

$$\varpi _2=d\gamma v_0^\nu =2dDv_0^{\left(\nu +2\right)}=d\gamma _0\kappa _{\nu +2},$$
(54)
are independent of $`\gamma `$ or $`D`$. So, we have used (49) to eliminate $`\gamma `$ and $`D`$, and to put it in the universal form, containing $`\varpi _2`$.

### 4.3 Comments

Equivalence free cooling - Gaussian thermostat:
Comparison of integral equation (51) for free cooling and (53) for the Gaussian thermostat shows that both equations are identical, as first observed by Montanero and Santos MS00 . This implies that the scaling form $`f\left(c\right)`$ in free cooling ($`=0`$) is identical to the NESS distribution $`f\left(c\right)`$ of the same system, driven by a Gaussian thermostat, provided both forms are rescaled to the same constant width $`𝑑\stackrel{}{c}c^2f\left(c\right)=d/2`$. This also implies that the scaling form for the free case can be measured by performing Monte Carlo simulations in a steady state brey ; MS00 . The same idea of systematically rescaling the velocities has also been used in molecular dynamics simulations of a freely cooling system of $`N`$ inelastic hard spheres Soto .

Universality:
The scaling equations (53) for $`f\left(c\right)`$ in systems driven by a Gaussian thermostat and by white noise, have a universal form, because $`\varpi _2`$ is independent of the friction constant $`\gamma `$, the noise strength $`D`$ and the width $`v_0\left(\mathrm{}\right)`$, which may depend on the initial distribution. The scaling equations reduce for $`\nu =1`$ to those for hard spheres TvN+ME-granmat , and for $`\nu =0`$ to those for Maxwell models ME+RB-rapid .

Perturbative approach:
The moment of the collision term $`\varpi _2=d\gamma _0\kappa _{\nu +2}`$ in the scaling equations (53), and the cooling rate $`\zeta _\nu =\gamma _0\kappa _{\nu +2}v_0^\nu `$ in (48) contain for all models the quantity $`\kappa _\nu `$, as given by (38). It depends on the unknown function $`f\left(c\right)`$, except for $`\nu =0`$ where $`\kappa _2=1`$ and $`\varpi _2=d\gamma _0=\frac{1}{4}\left(1\alpha ^2\right)`$. For the case of inelastic hard spheres, a perturbative method has been developed in goldstein-shapiro ; TvN+ME-granmat  for small inelasticities to solve both integral equations in (53) by expanding $`f\left(c\right)`$ in a series of Sonine polynomials $`S_p\left(c^2\right)`$, i.e.

$$f\left(c\right)=\varphi \left(c\right)\{1+\underset{p=2}{\overset{\mathrm{}}{}}a_pS_p\left(c^2\right)\},$$
(55)
where $`\varphi \left(c\right)=\pi ^{d/2}\mathrm{exp}\left(c^2\right)`$ is the Maxwellian. Here $`a_2`$ is essentially the fourth cumulant which has been calculated explicitly in TvN+ME-granmat . For dimensions $`d>1`$ it turns out to be proportional to the inelasticity $`\left(1\alpha ^2\right)`$. However, the one-dimensional case is exceptional because $`a_2`$ approaches a finite value for $`\alpha 1`$ BBRTvW01 . The same method has been successfully applied by Cercignani et al. CCG-00  to an inelastic Maxwell model, and by several authors brey ; MS00 ; brito+huthman ; Brill-TP-badhonnef  to inelastic hard spheres and related problems. The method can be applied to the inelastic IRS-models as well.

The method focuses on the lower moments of $`f\left(c\right)`$, and the polynomial approximation, cut off after $`S_2\left(c^2\right)`$, gives a fair representation of $`f\left(c\right)`$ for velocities in the thermal range, $`c2`$. The expansion (55) can then be used to calculate $`\kappa _\nu `$ in (52) in the form $`\kappa _\nu =\kappa _\nu ^0+\left(1\alpha ^2\right)\kappa _\nu ^1+\mathrm{}`$. For more details we refer to TvN+ME-granmat . As an illustration we calculate the lowest approximation $`\kappa _\nu ^0`$ to $`\kappa _\nu `$ in (38) using $`f\left(c\right)\varphi \left(c\right)`$. The result is,

$`\kappa _\nu ^0`$ $`=`$ $`{\displaystyle _\stackrel{}{n}}{\displaystyle 𝑑\stackrel{}{c}𝑑\stackrel{}{c}_1\left|g_{}\right|^\nu \varphi \left(c\right)\varphi \left(c_1\right)}`$ (56)
$`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x\left|x\right|^\nu e^{x^2/2}={\displaystyle \frac{2^{\nu /2}}{\sqrt{\pi }}}\mathrm{\Gamma }\left({\displaystyle \frac{\nu +1}{2}}\right).`$
The lowest approximation to (54) is then $`\varpi _2^0=d\gamma _0\kappa _\nu ^0`$. The damping rate $`\gamma _\nu `$ in (43) of the r.m.s. velocity becomes then $`\gamma _\nu =\left(\nu +1\right)\gamma _0`$, in agreement with the results for $`\nu =1`$ in TvN+ME-granmat .

### 4.4  High energy tails

The polynomial expansion (55) describes $`f\left(c\right)`$ only in the thermal range, but contains no meaningful information about velocities $`c`$ in the asymptotic range $`\left(c1\right)`$. However the high energy tails in the IRS-models with $`\nu >0`$ can be determined by a procedure similar to the one used successfully for inelastic hard spheres systems TvN+ME-granmat , as well as in inelastic Maxwell models driven by white noise ME+RB-rapid . To do so we make the ansatz of stretched Gaussian behavior for the high energy tail, i.e. $`f\left(c\right)B\mathrm{exp}\left[\beta c^b\right]`$ with $`0<b<2`$ and $`\beta >0`$, and determine $`b`$ and $`\beta `$ by inserting this ansatz in the scaling equation (53) and requiring self-consistency. The border line case $`b2`$ corresponds to Gaussian tails, and $`b0`$ suggests power law tails with negative exponents.

An estimate of the rescaled collision term $`I\left(c|f\right)`$ in (24) is made in TvN+ME-granmat ; ME+RB-rapid . This suggests that the loss term is asymptotically dominant over the gain term as long as the exponent $`b`$ in $`\mathrm{exp}\left[\beta c^b\right]`$ is restricted to $`b>0`$. This estimate also applies to the IRS-models as long as the exponent $`\nu >0`$ and the inelasticity $`\gamma _0`$ is non-vanishing. So the gain term is neglected. Moreover the loss term can be simplified in the asymptotic velocity range. Its dominant contribution comes from collisions with particles having velocities $`c_1`$ that are typically in the thermal range $`\left(c_1=𝒪\left(1\right)\right)`$. Consequently the collision rate $`\left|\left(\stackrel{}{c}\stackrel{}{c}_1\right)\stackrel{}{n}\right|^\nu `$ for asymptotic dynamics may be replaced by $`c^\nu \left|\widehat{\stackrel{}{c}}\stackrel{}{n}\right|^\nu `$, and the total collision term simplifies to,

$$I\left(c|f\right)I_{\text{loss}}\left(c|f\right)c^\nu \beta _\nu f\left(c\right),$$
(57)
with $`\beta _\nu `$ given by (40). After these preparations we insert (57) and the stretched Gaussian form, $`f\left(c\right)\mathrm{exp}\left[\beta c^b\right]`$, into the Boltzmann equation (53), and follow the procedure sketched above. This gives the following universal results (see item 2 below (54)) for the asymptotic high energy tail in $`d`$dimensional inelastic $`\nu `$models, $`f\left(c\right)\mathrm{exp}\left[\beta c^b\right]`$ with

$$\begin{array}{ccc}b=0\hfill & \text{inconsistent}\hfill & \left(\text{GT};\nu =0\right)\hfill \\ b=\nu \hfill & \beta =d\beta _\nu /\nu \varpi _2\hfill & \left(\text{GT};\nu >0\right)\hfill \\ b=\frac{1}{2}\left(\nu +2\right)\hfill & \beta =\frac{2}{\nu +2}\sqrt{\frac{2d\beta _\nu }{\varpi _2}}\hfill & \left(\text{WN};\nu >0\right),\hfill \end{array}$$
(58)
where the scaling functions and high energy tails of GT-driven and freely cooling systems are equivalent.

We conclude this subsection with some comments.
driven systems:
For $`0<b<2`$ both GT- and WN-driving lead to consistent asymptotic solutions of the scaling equations for $`\nu >0`$ with overpopulated high energy tails of stretched Gaussian type. We also note, the larger the interactions (larger $`\nu `$-values) at large impact energies, the smaller the overpopulation of the high energy tails. For $`\nu =2`$, corresponding to the inelastic version of the Very Hard Particle model (see end of Section 3.2) we obtain Gaussian behavior $`(b=2)`$ for both types of thermostats, and there are no longer over-populated tails. In the case of stretched Gaussian tails all moments $`𝑑\stackrel{}{c}c^nf\left(c\right)<\mathrm{}`$. This would not be the case for power law tails.

In the case of white noise driving, the above results with $`b=1+\frac{1}{2}\nu `$, include inelastic hard spheres $`\left(\nu =1\right)`$, as well as Maxwell models ($`\nu =0`$), and the results coincide with the detailed predictions for $`\nu =1`$ TvN+ME-granmat  and $`\nu =0`$ ME+RB-rapid . For GT-driven or freely cooling IRS-models with $`\nu >0`$ the tail distributions have a stretching exponent $`b=\nu `$. The property $`\beta b1/\gamma _0`$ as $`\nu 0`$ would suggest power law behavior at the stability threshold, $`\nu =0`$ (Maxwell models), with an exponent $`a=1/\gamma _0`$ (see however next item). The result (58) also includes the exponential decay, $`f(c)\mathrm{exp}[\beta c]`$ for inelastic hard spheres, found in TvN+ME-granmat .

Stability thresholds:
In the previous comments we have speculated on power law tails, in agreement with the exact solution Rome1  for the one-dimensional freely cooling Maxwell models with $`f\left(c\right)c^4`$ for $`c1`$. As will be shown in section 5, the scaling form for $`d`$-dimensional Maxwell models in free cooling does indeed have an algebraic tail $`f\left(c\right)1/c^{d+a}`$, but with an exponent $`a(\alpha )1/\gamma _0`$. However, neither the existence of algebraic tails, nor the possible value of the power law exponents can be obtained legitimately from the results (58). The reason is that the result is not consistent with the a priori  assumption $`I_{\text{gain}}I_{\text{loss}}`$ at $`c1`$. In fact, gain and loss terms are of the same order of magnitude for $`\nu =0`$.

Maxwell models as approximations to hard spheres:
As discussed in item 2 of Section 3.5, inelastic Maxwell models have also been introduced as a sensibly looking mathematical simplification of the Boltzmann equation for inelastic hard spheres. The results of the analysis in the previous section show that the effect of this simplification on the shape of the tail may be very drastic. In the free cooling case the hard sphere tail is exponential, $`f\left(c\right)\mathrm{exp}\left[\beta c\right]`$, whereas in the Maxwell model it is a power law tail, $`f\left(c\right)1/c^{d+a}`$. A smaller difference exists in the WN-driven case, where the hard sphere tail is $`f\left(c\right)\mathrm{exp}\left[\beta c^{3/2}\right]`$, and the Maxwell tail is $`f\left(c\right)\mathrm{exp}\left[\beta c\right]`$.

MC-simulations at high and low inelasticity:
The results (58) predict more than the stretching exponent $`b`$. In fact the coefficient $`\beta `$ for $`\nu >0`$ in $`f\left(c\right)\mathrm{exp}\left[\beta c^b\right]`$ is given in terms of $`\varpi _2`$ in (54). For these $`\nu `$values we have described a perturbative calculation for $`\varpi _2`$, which converges rapidly for $`\alpha 1`$, but gives poor results for $`\alpha <0.6`$ TvN+ME-granmat ; MS00 . Nevertheless, it is possible to test the predictions of (58) for the GT- and the WN-thermostats, including the coefficient $`\beta `$ for all values of the restitution coefficient $`\alpha `$. This can be done by measuring the mean square velocity $`v_0^2\left(\mathrm{}\right)`$ in MC simulations on driven systems in a NESS for a given $`\gamma `$ or $`D`$. From these data $`\varpi _2`$ can be calculated using (54).

## 5 Inelastic Maxwell models

### 5.1 Fourier transform method

As discussed in Section 4.3, the behavior of the high energy tails of the scaling form in the IRS-models is only controlled by the loss term in the Boltzmann collision term, and generates the stretched Gaussian tails. However, in the borderline case ($`\nu =0)`$ the tail behavior is determined by an interplay between gain and loss term, which leads to algebraic tails. This makes the analysis more complicated. In this article we only discuss the borderline case in free cooling, formed by the Maxwell models, and we refer to ETB  for a more comprehensive discussion.

In this section we want to demonstrate that the Boltzmann equation (24) for inelastic Maxwell models ($`\nu =0`$),

$$_tF(v,t)=F(v,t)+\frac{1}{\alpha }_\stackrel{}{n}𝑑\stackrel{}{w}F(v^{},t)F(w^{},t),$$
(59)
has a scaling solution with a power law tail. To do so we first consider the Fourier transform of the distribution function, $`\mathrm{\Phi }(k,t)=\mathrm{exp}\left[i\stackrel{}{k}\stackrel{}{v}\right]`$, which is the characteristic function or generating function of the velocity moments. Because $`F(v,t)`$ is isotropic, $`\mathrm{\Phi }(k,t)`$ is isotropic as well.

As an auxiliary step we first Fourier transform the gain term in (59), i.e.

$$\begin{array}{c}𝑑\stackrel{}{v}\mathrm{exp}\left[i\stackrel{}{k}\stackrel{}{v}\right]I_{\text{gain}}\left(v|F\right)=\hfill \\ \hfill _\stackrel{}{n}𝑑\stackrel{}{v}𝑑\stackrel{}{w}\mathrm{exp}\left[i\stackrel{}{k}\stackrel{}{v}^{}\right]F(v,t)F(w,t)=_\stackrel{}{n}\mathrm{\Phi }(k\eta _+,t)\mathrm{\Phi }(k\eta _{},t).\end{array}$$
(60)
The transformation needed to obtain the first equality is the same as in (37). Then we use (14) to write the exponent as $`\stackrel{}{k}\stackrel{}{v}_1^{}=\stackrel{}{k}_{}\stackrel{}{v}_1+\stackrel{}{k}_+\stackrel{}{v}_2`$, where

$$\begin{array}{cc}\hfill \stackrel{}{k}_+& =p\stackrel{}{k}\stackrel{}{n}\stackrel{}{n}\left|\stackrel{}{k}_+\right|=kp\left|\left(\widehat{\stackrel{}{k}}\stackrel{}{n}\right)\right|=k\eta _+\left(\stackrel{}{n}\right)\hfill \\ \hfill \stackrel{}{k}_{}& =\stackrel{}{k}\stackrel{}{k}_+\left|\stackrel{}{k}_{}\right|=k\sqrt{1z\left(\widehat{\stackrel{}{k}}\stackrel{}{n}\right)^2}=k\eta _{}\left(\stackrel{}{n}\right),\hfill \end{array}$$
(61)
with $`p=1q=\frac{1}{2}\left(1+\alpha \right)`$ and $`z=1q^2`$. The Fourier transform of (59) then becomes,

$$_t\mathrm{\Phi }(k,t)=\mathrm{\Phi }(k,t)+_\stackrel{}{n}\mathrm{\Phi }(k\eta _+\left(\stackrel{}{n}\right),t)\mathrm{\Phi }(k\eta _{}\left(\stackrel{}{n}\right),t),$$
(62)
where $`\mathrm{\Phi }(0,t)=1`$ because of (25). In one-dimension $`\eta _+\left(\stackrel{}{n}\right)=p`$ and $`\eta _{}\left(\stackrel{}{n}\right)=q`$, and $`_\stackrel{}{n}`$ can be replaced by unity. Moreover this equation simplifies to BN-PK-00 ,

$$_t\mathrm{\Phi }(k,t)=\mathrm{\Phi }(pk,t)\mathrm{\Phi }(qk,t)\mathrm{\Phi }(k,t).$$
(63)
Because $`F(v,t)`$ is isotropic, only the even moments are non-vanishing, and the moment expansion of the characteristic function takes the form,

$$\mathrm{\Phi }(k,t)=\underset{n}{}^{}\frac{\left(ik\right)^n}{n!}\left(\widehat{\stackrel{}{k}}\stackrel{}{v}\right)^n=\underset{n}{}^{}\left(ik\right)^nm_n\left(t\right),$$
(64)
where the prime indicates that $`n=even`$, and the moment $`m_n\left(t\right)`$ is defined as,

$$m_n\left(t\right)=\beta _nv^n/n!,$$
(65)
where $`\beta _n=_\stackrel{}{n}\left(\widehat{\stackrel{}{k}}\widehat{\stackrel{}{v}}\right)^n`$ is given in (40). Moreover, the normalizations (25) give $`m_0\left(t\right)=1`$ and $`m_2\left(t\right)=\frac{1}{2}\beta _2v^2=\frac{1}{4}v_0^2`$.

Maxwell models have the unusual property that the system of moment equations for $`m_n\left(t\right)`$ is closed, and can be solved sequentially. The reason is that the rate $`\dot{m}_n`$ is only a function of lower moments $`m_s\left(t\right)`$ with $`sn`$.

The moment equations are readily obtained by inserting the expansion (64) in the Fourier transformed Boltzmann equation (62), and equating the coefficients of equal powers of $`k`$. The result is,

$$\dot{m}_n+\lambda _nm_n=\underset{l=2}{\overset{n2}{}}h(l,nl)m_lm_{nl}\left(n>2\right),$$
(66)
where all labels $`\{n,l,s\}`$ take even values only. Following BN+PK-02  or Appendix A of ME+RB-fest  the functions can be calculated for real positive values of $`\mathrm{}`$ and $`s`$ with the result,

$$\begin{array}{cc}\hfill h(l,s)=& _\stackrel{}{n}\eta _+^l\left(\stackrel{}{n}\right)\eta _{}^s\left(\stackrel{}{n}\right)=p^l\beta _l{}_{2}{}^{}F_{1}^{}(\frac{s}{2},\frac{l+1}{2};\frac{l+d}{2}|z)\hfill \\ \hfill \lambda _s=& 1h(s,0)h(0,s)=_\stackrel{}{n}\left[1\eta _+^s\left(\stackrel{}{n}\right)\eta _{}^s\left(\stackrel{}{n}\right)\right]\hfill \\ \hfill =& 1p^s\beta _s{}_{2}{}^{}F_{1}^{}(\frac{s}{2},\frac{1}{2};\frac{d}{2}|\mathrm{\hspace{0.25em}1}q^2).\hfill \end{array}$$
(67)
Here $`{}_{2}{}^{}F_{1}^{}(\alpha ,\beta ;\gamma |z)`$ is a hypergeometric function, and $`\lambda _2=2pq/d=2\gamma _0`$ is easily obtained from the above results.

Next we consider the Fourier transform of the scaling relation (26), yielding $`\mathrm{\Phi }(k,t)=\varphi \left(kv_0\left(t\right)\right)`$ with $`v_0\left(t\right)=v_0\left(0\right)\mathrm{exp}\left(\gamma _0t\right)`$ according to (41). Inserting $`\mathrm{\Phi }`$ in (62) gives the integral equation for the scaling form $`\varphi \left(k\right)`$, which reads

$$\gamma _0k\frac{d}{dk}\varphi \left(k\right)+\varphi \left(k\right)=_\stackrel{}{n}\varphi \left(k\eta _+\right)\varphi \left(k\eta _{}\right).$$
(68)
For $`d=1`$ it reduces to,

$$pqk\frac{d}{dk}\varphi \left(k\right)+\varphi \left(k\right)=\varphi \left(pk\right)\varphi \left(qk\right).$$
(69)
The scaling form $`\varphi \left(k\right)`$ is the generating function for the moments of $`f\left(c\right)`$, i.e.

$$\begin{array}{cc}\hfill \varphi \left(k\right)=& \underset{n}{}^{}\frac{\left(ik\right)^n}{n!}\beta _nc^n\underset{n}{}^{}\left(ik\right)^n\mu _n\hfill \\ \hfill & 1\frac{1}{4}k^2+k^4\mu _4k^6\mu _6+\mathrm{},\hfill \end{array}$$
(70)
where $`n=even`$, $`\mu _0=1`$ and $`\mu _2=\frac{1}{2}\beta _2c^2=1/4`$ on account of the normalizations (27) and $`\beta _2=1/d`$ (see (40)). By inserting (70) into (68) one obtains the recursion relation,

$$\mu _n=\frac{1}{\lambda _n\frac{1}{2}n\lambda _2}\underset{l=2}{\overset{n2}{}}h(l,nl)\mu _l\mu _{nl}\left(n>2\right)$$
(71)
with initialization $`\mu _2=1/4`$ and all labels $`n,l=even`$. How these moments behave as a function of $`\alpha `$ has been calculated in ME+RB-fest  by numerically solving the recursion relation.

### 5.2 Small-$`k`$ singularity of the characteristic function

The asymptotic analysis in Section 4.4 of the high energy tail $`f\left(c\right)\mathrm{exp}\left[\beta c^b\right]`$ for Maxwell models ($`\nu =0`$) yields $`b=0`$, which suggests a power law tail $`f\left(c\right)1/c^{d+a}`$ as the leading large$`c`$ behavior in $`d`$dimensional systems. If this is indeed the case, then the moments $`\mu _n`$ of the scaling form $`f\left(c\right)`$ are convergent if $`n<a`$ and are divergent if $`n>a`$. As we are interested in physical solutions which can be normalized, and have a finite energy, a possible value of the power law exponent must obey $`a>2`$.

The characteristic function is in fact a very suitable tool for investigating this problem. Suppose the moment $`\mu _n`$ with $`n>a`$ diverges, then the $`n`$-th derivative of the corresponding generating function also diverges at $`k=0`$, i.e. $`\varphi \left(k\right)`$ has a singularity at $`k=0`$. Then a simple rescaling argument of the inverse Fourier transform shows that $`\varphi \left(k\right)`$ has a dominant small-$`k`$ singularity of the form $`\varphi \left(k\right)\left|k\right|^a`$, where $`aeven`$. On the other hand, whenever the exponent $`b`$ in $`f\left(c\right)\mathrm{exp}\left[\beta c^b\right]`$ is positive, – as is the case for the inelastic Maxwell model driven by white noise, where $`b=1`$ and $`\beta =\sqrt{2/pq}`$ – , then all moments are finite, and the characteristic function $`\varphi \left(k\right)`$ is regular at the origin, i.e. can be expanded in powers of $`k^2`$.

We first illustrate our analysis for the one-dimensional case. As the requirement of finite total energy imposes the lower bound $`a>2`$ on the exponent, we make the ansatz, consistent with (70) that the dominant small-$`k`$ singularity has the form,

$$\varphi \left(k\right)=1\frac{1}{4}k^2+A\left|k\right|^a,$$
(72)
insert this in (69), and equate the coefficients of equal powers of $`k`$. This yields,

$$\frac{1}{2}a\lambda _2=\lambda _a\mathrm{or}apq=1p^aq^a.$$
(73)
The equation has two roots, $`a=2,3`$, of which $`a=3`$ is the one larger than 2. Here $`A`$ is left undetermined. Consequently the one-dimensional scaling solution has a power law tail, $`f\left(c\right)1/c^4`$, in agreement with the exact solution in Rome1 .

For general dimension we proceed in the same way as in the one-dimensional case, insert the ansatz (72) into (68), and equate the coefficients of equal powers of $`k`$. This yields for the coefficient of $`k^2`$ the identity $`2\gamma _0=\lambda _2`$, and for the coefficient of $`k^a`$ the transcendental equation,

$$\frac{1}{2}a\lambda _2=\lambda _a=_\stackrel{}{n}\left[1\eta _+^a\eta _{}^a\right].$$
(74)
The equation above obviously has the solution $`a=2`$. We are however interested in the solution with $`a>2`$. In the elastic limit ($`\alpha 1`$) the solution is simple. There $`\gamma _00`$ and $`a`$ diverges. The contributions of $`\eta _\pm ^a`$ on the right hand side vanish because $`\eta _\pm <1`$, and the result is,

$$a\frac{1}{\gamma _0}=\frac{4d}{1\alpha ^2}.$$
(75)
For general values of $`\alpha `$ one can conveniently use an integral representation of $`{}_{2}{}^{}F_{1}^{}`$ to evaluate $`\lambda _a`$ and solve the transcendental equation (74) numerically. We illustrate the solution method of (74) with the graphical construction in Figure 2, where we look for intersections of the line $`y=\frac{1}{2}\lambda _2=\gamma _0s`$ with the curve $`y=\lambda _s`$ for different values of $`\alpha `$.

The relevant properties of $`\lambda _s`$ are: (i) $`lim_{s0}\lambda _s=1`$; (ii) $`\lambda _s`$ is a concave function, monotonically increasing with $`s`$, and (iii) all eigenvalues for positive integers $`n`$ are positive (see Figure 2). As can be seen from the graphical construction, the transcendental equation (74) has two solutions, the trivial one $`\left(s_0=2\right)`$ and the solution $`s_1=a`$ with $`a>2`$. The numerical solutions for $`d=2,3`$ are shown in Figure 3 as a function of $`\alpha `$, and the $`\alpha `$-dependence of the root $`a\left(\alpha \right)`$ can be understood from the graphical construction. In the elastic limit as $`\alpha 1`$ the eigenvalue $`\lambda _2\left(\alpha \right)0`$ because of energy conservation. In that limit the transcendental equation (74) no longer has a solution with $`a>2`$, and $`a\left(\alpha \right)\mathrm{}`$ according to (75), as it should be. This is consistent with a Maxwellian tail distribution in the elastic case. Krapivsky and Ben-Naim have in fact solved the transcendental equation asymptotically for large $`d`$, which gives qualitatively the same results as shown in Figure 3 for two and three dimensions.

These results establish the existence of scaling solutions $`f\left(c\right)1/c^{d+a}`$ with algebraic tails, where the exponent $`a`$ is the solution of the transcendental equation (74) with $`a>2`$. The exponent $`a\left(\alpha \right)`$ behaves as a function of $`\alpha `$ qualitatively the same as in the simple inelastic BGK model of Section 2.

Using a somewhat different analysis Krapivsky and Ben-Naim BN+PK-6-11  independently obtained the same results for the algebraic tails in freely cooling Maxwell models.

### 5.3 Beyond asymptotic analysis

Free cooling:
As mentioned already in Section 3.3, Baldassarri et al. Rome1  have obtained an exact solution $`f\left(c\right)`$ of the scaling equation for the one-dimensional Maxwell model in free cooling, and they demonstrated the importance of this scaling solution by means of Monte Carlo simulations. In doing so they found that $`F(v,t)`$ for different classes of initial distributions can be collapsed for long  times $`t`$ on this exact scaling solution $`f\left(c\right)`$, when $`v_0\left(t\right)F(v,t)`$ is plotted versus $`c=v/v_0\left(t\right)`$.

We briefly illustrate here how this solution is obtained from the one-dimensional scaling equation (69). One first verifies that the following function,

$$\varphi \left(k\right)=\left(1+\vartheta \left|k\right|\right)\mathrm{exp}\left(\vartheta \left|k\right|\right),$$
(76)
with arbitrary positive $`\vartheta `$ is a solution, that can be Fourier inverted. It gives the scaling form,

$$f\left(c\right)=\frac{2}{\pi \vartheta }\frac{1}{\left(1+c^2/\vartheta ^2\right)^2}.$$
(77)
To determine the scaling form that satisfies the normalizations (27), we expand (76) in powers of $`\left|k\right|`$ to obtain,

$$\varphi \left(k\right)=1\frac{1}{2}\vartheta ^2k^2+\frac{1}{3}\vartheta ^3\left|k\right|^3+\mathrm{}.$$
(78)
Comparison of this result with (70) shows that $`\vartheta =1/\sqrt{2}`$. Moreover, it confirms the ansatz (72) used to find solutions with small-$`k`$ singularities. Comparison also shows the value of the coefficient $`A=1/\left[6\sqrt{2}\right]`$ in (72), and the high energy tail is $`f\left(c\right)𝒜/c`$ with $`𝒜=1/\left[\pi \sqrt{2}\right]`$. The coefficients $`A`$ and $`𝒜`$ can not be determined within the asymptotic method.

White noise driving:
For the one-dimensional Maxwell model driven by white noise the steady state solution has also been found exactly with the help of the Fourier transform method BN-PK-00 . The characteristic function satisfies in that case a nonlinear finite difference equation, which can be solved by iteration. However, the analytic structure is rather complex, and makes it difficult to extract analytic information from that solution.

The observation that $`f\left(c\right)\mathrm{exp}\left[\beta \left|c\right|\right]`$ has an exponential high energy tail, was probably first made in numerical work of van der Hart and Nienhuis Nienhuis . A more detailed analytic prediction about the asymptotic tail was given in ME+RB-rapid , where it was shown that $`\beta =\sqrt{8/\left(1\alpha ^2\right)}`$ for all Maxwell models, independent of the dimensionality $`\left(d=1,2,\mathrm{}\right)`$, at least with the normalization, $`c^2=\frac{1}{2}d`$, used in this article. This result follows directly from (58) and the value $`\varpi _2=d\gamma _0=\frac{1}{4}\left(1\alpha ^2\right)`$, given in item 3 below (54). Furthermore, additional numerical and analytical work was also published by Marconi and Puglisi Rome3 , and by Antal et al. Droz . Only recently more detailed analytic results have been extracted from the rather complex structure of the exact solution BN-PK-Springer ; AS+ME-condmat .

## 6 Conclusions and perspectives

We have studied asymptotic properties of scaling or similarity solutions, $`F(v,t)=\left(v_0\left(t\right)\right)^df\left(v/v_0\left(t\right)\right)`$, of the nonlinear Boltzmann equation in spatially homogeneous systems composed of particles with inelastic interactions for large times and large velocities. The large $`t`$ and $`v`$scales are relevant because on such scales the universal features of the solutions survive, while details of the initial distributions, of interaction strength, and degree of inelasticity, are mostly lost. The behavior of these scaling states, which describe nonequilibrium steady states (NESS), is less universal than the state of thermal equilibrium, because the form of the NESS distribution $`f(c)`$ depends on the way of driving the dissipative systems. Scaling solutions are very well suited to expose the universal features of the velocity distribution functions, because the velocities, $`c=v/v_0\left(t\right)`$, are measured in units of the r.m.s. velocity or instantaneous width $`v_0\left(t\right)`$ of the distribution.

The real importance of the scaling solutions is that the actual solutions $`F(v,t)`$ for large classes of initial distributions $`F(v,0)`$ (essentially initial data without over-populated tails -see ME+RB-fest ) rapidly approach these scaling solutions in the sense that after a short transient time the data $`v_0^d\left(t\right)F(v,t)`$ can be collapsed on a single scaling form $`f\left(c\right)`$, as first observed by Baldassarri et al. in MC simulations of the nonlinear Boltzmann equation for a one-dimensional inelastic Maxwell model. This rapid approach to a universal scaling function applies both to systems driven by a Gaussian or by a white noise thermostat, as well as for freely cooling systems, which show scaling behavior, identical to systems driven by a Gaussian thermostat.

Originally this scenario had the status of a conjecture for systems of inelastic hard spheres ($`\nu =1`$) and inelastic Maxwell models($`\nu =0`$), as formulated in ME+RB-fest  where also some analytical evidence for the approach to a scaling form has been presented. For Maxwell models the conjecture has been rigorously proven in the mean time by Bobylev et al. AB+CC-proof . The analysis in Section 4.1 also suggests what the basic criterion us for the approach of IRS-models to the scaling form. If the energy balance equation (48) has a stable/attractive fixed point solution, $`v_0(\mathrm{})`$, then the distribution function $`F(v,t)`$ approaches – in the sense detailed in the paragraphs above – a stable NESS described by the scaling solution.

The analysis in Section 4.1 of the energy balance equation also suggests that cases of marginal stability are candidates for power law tails, where the freely evolving inelastic Maxwell model is a well known example. A further candidate system for power law tails is an IRS-model with $`\nu =2`$, driven by white noise. These properties, and the relation between marginal stability and power law tails is further explored in ETB , where in addition the analytic predictions of the present paper are confirmed by MC simulations.

In this article we have focussed on the properties of the scaling solution, and in particular on its high energy tail. Here we have introduced the Boltzmann equation for new classes of inelastic interactions, named Inelastic Repulsive Scatterers or IRS-models, corresponding to pseudo-repulsive power law potentials, with collision probabilities proportional to $`\text{a}\left|\stackrel{}{g}\stackrel{}{n}\right|^\nu `$, covering hard scatterers like inelastic hard spheres ($`\nu =1`$) and soft scatterers like pseudo-Maxwell molecules ($`\nu =0`$). The energy loss in an inelastic interaction is proportional to the inelasticity, $`\gamma _0\left(1\alpha ^2\right)`$, where $`\alpha `$ is the coefficient of restitution. We have studied two typical cases: freely evolving systems in homogeneous cooling states, without energy supply, and systems with energy supply, driven by Gaussian thermostats with negative friction or driven by Gaussian white noise.

The homogeneous cooling laws in these systems are described by the granular temperature $`T=v_0^2`$, with a long time decay as $`t^{2/\nu }`$ or $`\mathrm{exp}\left[2\gamma _\nu \tau \right]`$, where $`t`$ is the external (laboratory) time, and $`\tau `$ is the internal time or collision counter, and $`\gamma _\nu `$ is proportional to the inelasticity $`\gamma _0`$.

An interesting feature of all IRS-models is that their scaling forms generically have a stretched Gaussian high energy tail, $`f\left(c\right)\mathrm{exp}[\beta c^b]`$. In freely cooling IRS-models with $`\nu >0`$ the stretching exponent $`b=\nu `$ satisfies $`0<b=\nu 2`$. Among the IRS-models at the stability threshold only the freely cooling Maxwell model has been analyzed in the present paper, and it does yield power law tails, $`f\left(c\right)1/c^{d+a}`$, and the exponent $`a`$ has been calculated from a transcendental equation. In the IRS-models with $`\nu >0`$, driven by white noise, there exist again stretched Gaussian tails with $`b=\frac{1}{2}\left(\nu +2\right)`$. All known results for higher dimensional inelastic hard spheres ($`\nu =0`$) and inelastic Maxwell models ($`\nu =1`$) are recovered in the present analysis.

In IRS-models with positive $`\nu `$ the tails are always stretched Gaussians, and they are determined only by the loss term in the Boltzmann equation. In the freely cooling Maxwell model, and more generally at a stability threshold – which are the cases leading to power law tails ETB  – the loss and gain term in the nonlinear Boltzmann equation are of comparable size, and partially balancing each other.

We also have analyzed in Section 2 an inelastic BGK model, introduced by Brey et al BMD , and generalized it to an energy dependent collision frequency, and we added energy source terms to this kinetic equation as well. In all cases the scaling solution can be calculated exactly. The freely cooling model shows an algebraic tail, $`f\left(c\right)1/c^{d+a}`$ with an exponent $`a1/\gamma _02/\left(1\alpha ^2\right)`$. For white noise driving one finds asymptotically $`f\left(c\right)\mathrm{exp}\left[\beta c\right]`$ with $`\beta \left(1\alpha ^2\right)^{1/2}`$.

It is interesting to observe that the overpopulated high energy tails in the inelastic BGK-models, both for free cooling as well as for white noise driving, are essentially the same as for the more complicated Maxwell models. However, the present extension to an energy dependent collision frequency does not capture the generic features of the less-singular stretched Gaussian tails for the IRS-models with $`\nu >0`$ , as predicted by the nonlinear Boltzmann equation for these models.

## Acknowledgements

The authors want to thank J.W. Dufty, A. Santos, E. Trizac, A. Barrat and K. Shundyak for helpful discussions and correspondence. This work is supported by DGES (Spain), Grant No BFM-2001-0291.