Singular Energy Distributions in Granular Media
E. Ben-Naim1, ∗ and A. Zippelius2, †
Theoretical Division and Center for Nonlinear Studies,
Los Alamos National Laboratory, Los Alamos, New Mexico 87545 USA
Institut für Theoretische Physik, Georg-August-Universität, 37077 Göttingen, Germany
We study the kinetic theory of driven granular gases, taking into account both translational and
rotational degrees of freedom. We obtain the high-energy tail of the stationary bivariate energy
distribution, depending on the total energy E and the ratio x =
Ew/E of rotational energy Ew to
total energy. Extremely energetic particles have a unique and well-defined distribution f(x) which
has several remarkable features: x is not uniformly distributed as in molecular gases; f(x) is not
smooth but has multiple singularities. The latter behavior is sensitive to material properties such
as the collision parameters, the moment of inertia and the collision rate. Interestingly, there are
preferred ratios of rotational-to-total energy. In general, f(x) is strongly correlated with energy
and the deviations from a uniform distribution grow with energy. We also solve for the energy
distribution of freely cooling Maxwell Molecules and find qualitatively similar behavior.
PACS numbers: 45.70.Mg, 47.70.Nd, 05.40.-a, 81.05.Rm
I. INTRODUCTION
Energy dissipation has profound consequences in granular materials, especially in dilute gases, where the dynamics
are controlled by collisions [1, 2, 3]. Dissipation is responsible for many interesting collective phenomena including
clustering [4, 5, 6, 7, 8], formation of shocks [9, 10, 11, 12, 13], and hydrodynamic instabilities [14, 15]. Another
consequence is the anomalous statistical physics that includes the non-Maxwellian velocity distributions [16, 17, 18,
19, 20, 21] and the breakdown of energy equipartition in mixtures [22, 23].
For an elastic gas in equilibrium, the temperature, defined as the average kinetic energy, characterizes the entire
distribution function including all of the moments, the bulk of the distribution, as well as the tail of the distribution.
Outside of equilibrium, the temperature is not sufficient to characterize the energy distribution. Granular gases are
inherently out of equilibrium and a complete characterization must therefore include the behavior of typical particles,
the behavior of energetic particles, as well as the moments of the distribution. For example, the energy distribution
may have power-law tails with divergent high-order moments [24, 25, 26] and consequently, the moments exhibit
multiscaling [27]. Generally, nonequilibrium effects are pronounced in the absence of energy input to balance the
dissipation but can be suppressed by injection of energy where the deviation from a Maxwellian distribution affects
only extremely energetic particles [17, 28, 29, 30].
While there is substantial understanding of the energy distribution of frictionless granular gases, much less is known
theoretically [31, 32, 33, 34, 35, 36, 37, 38] and experimentally [39, 40, 41] when the rotational degrees of freedom are
taken into account. It is difficult to measure the rotational motion experimentally, and the few available measurements
are restricted to two-dimensions. Surface roughness and friction have important consequences and the hydrodynamic
theory [42, 43, 44, 45] must be modified, if the particles have spin [46]. Equipartition does not hold for the average
rotational and translational temperature – neither in the free cooling case [33, 34, 35, 36] nor for a driven system [37].
In general, rotational and linear degrees of freedom are correlated in direction [47].
In this paper, we investigate the nature of the full energy distribution, that is, the bivariate distribution of rotational
and translational energy. Motivated by the fact that on average the total energy is not partitioned equally between
rotational and translational degrees of freedom, we focus on the bivariate distribution P (E, x) of total energy E and
the modified ratio x =
Ew/E of rotational to total energy. We thereby generalize the understanding of frictionless
granular matter in terms of the energy distribution to rough grains.
Our starting point is the nonlinear Boltzmann equation with a collision rule that accounts for the coupling of
translational and rotational motion due to tangential restitution. We study stationary solutions of the inelastic
Boltzmann equation that describe steady states achieved through a balance between energy injections that are powerful
but rare and energy dissipation through inelastic collisions. For high-energy particles we derive a linear equation for
∗Electronic address: ebn@lanl.gov
†Electronic address: annette@theorie.physik.uni-goettingen.de
http://arxiv.org/abs/0704.1520v1
mailto:ebn@lanl.gov
mailto:annette@theorie.physik.uni-goettingen.de
the bivariate energy distribution. The latter can be shown to factorize – P (E, x) = p(E)f(x) – into a product of the
distribution of the total energy, p(E), and the distribution of the fraction of energy stored in the rotational degrees
of freedom, f(x). The former distribution decays algebraically with energy: p(E) ∼ E−ν . The fraction of energy
stored in rotational motion is universal for energetic particles in the sense that f(x) approaches a limiting distribution
independent of energy. Furthermore, this quantity has a number of interesting features. First, the distribution is not
uniform, as it would be, if equipartition were to hold. Second, the distribution is not analytic but has singularities
at special energy ratios. Third, the distribution and in particular its singularities depend sensitively on the moment
of inertia and the collision parameters. Only for energetic particles is this distribution well defined. In general, the
partition of energy into rotational and translational motion depends on the magnitude of the energy. This paper
specifically addresses two-dimensions, although the theoretical approach and the reported qualitative behavior are
generic.
We also develop a general framework for describing high-energy collisions and we use this framework to study freely
cooling Maxwell Molecules where the moments of the energy distribution can be found in a closed form. For example,
the two granular temperatures corresponding to the rotational and translational motions are coupled and generally,
they are not equal. The high-energy behavior found for driven steady-states extends to freely cooling gases.
The rest of this paper is organized as follows. We review the collision rules and introduce the nonlinear kinetic
theory in section II. We then derive the linear kinetic theory for high-energy particles in section III. Next, in section
IV, we study driven steady states and solve for the stationary energy distribution. Freely cooling Maxwell molecules
are discussed in section V and we conclude in section VI. The Appendices detail technical derivations.
II. THE NONLINEAR KINETIC THEORY
Our system consists of an infinite number of identical particles with mass m = 1, radius R, and moment of inertia
I = qR2 where 0 ≤ q ≤ 1 is a dimensionless quantity. Each particle has a linear velocity v and an angular velocity
w. Its total energy is shared by the linear and the rotational motion, E = Ev + Ew, or explicitly,
v2 + qR2w2
where v ≡ |v| and w ≡ |w|.
In a collision between two particles, their velocities (vi,wi) with the labels i = a, b, change according to
(va,wa) + (vb,wb) → (v′a,w′a) + (v′b,w′b) (2)
where the postcollision velocities are denoted by primes. In a binary collision, rotational and translational energy
are exchanged, while the total energy decreases. In this study, we consider tangential restitution in addition to the
standard normal restitution. Let ri be the position of particle i, then the directed unit vector connecting the centers
of the colliding particles is n̂ = (ra − rb)/|ra − rb|. We term this vector the impact direction. The collision rules are
most transparent in terms of ui the particle velocity at the contact point
ua = va +R n̂×wa (3a)
ub = vb −R n̂×wb. (3b)
The inelastic collision laws state that the normal component of the relative velocity U = ua − ub is reversed and
reduced by the normal restitution coefficient 0 ≤ rn ≤ 1. The tangential component is either reversed (rough particles)
or not (smooth particles) and in any case reduced by the tangential restitution coefficient −1 ≤ rt ≤ 1, according to
the following collision rules:
′ · n̂ = −rnU · n̂, (4a)
′ × n̂ = −rt U× n̂. (4b)
Inelastic collisions conserve linear and angular momentum. Conservation of linear momentum implies that the
total linear velocity does not change, and conservation of angular momentum enforces that the angular momentum
of each particle with respect to the point of contact remains the same, because there is no torque acting at the point
of contact. The collision laws (4) combined with these conservation laws specify the postcollision velocities as linear
combinations of the precollision velocities [33]
a = va − ηnV · n̂ n̂− ηt (V −V · n̂ n̂)− ηtR n̂×W w′a = wa +
n̂×V + ηt
n̂× n̂×W (5a)
b = vb + ηnV · n̂ n̂+ ηt (V −V · n̂ n̂) + ηtR n̂×W w′b = wb +
n̂×V + ηt
n̂× n̂×W (5b)
where the shorthand notations V = va − vb and W = wa +wb were introduced. These collision rules involve the
normal and tangential collision parameters, defined as
1 + rn
, and ηt =
1 + q
1 + rt
. (6)
Their range of values is bounded by 1/2 ≤ ηn ≤ 1 and 0 ≤ ηt ≤ q/(1 + q). Details of the derivation of the collision
rules are given in Appendix A, as they are relevant for our discussion. The energy dissipation, ∆E = Ea+Eb−E′a−E′b,
is given by
1− r2n
(V · n̂)2 + q
1 + q
1− r2t
(V −V · n̂ n̂+R n̂×W)2. (7)
The energy dissipation is always positive, except when the collisions are elastic, rn = 1 and rt = −1 (perfectly smooth
spheres) or rt = 1 (perfectly rough spheres).
The collision rate K(va,vb) is the rate by which the two particles approach each other. For hard spheres, this rate
is simply the normal component of the relative velocity, but we study the general case
K(va,vb) = |(va − vb) · n̂|γ (8)
with 0 ≤ γ ≤ 1. Of course, the collision rate vanishes, K = 0, when the particles are moving away from each other,
(va − vb) · n̂ > 0. When particles interact via the central potential r−κ then γ = 1 − 2 d−1κ [48]. The two limiting
cases are hard spheres (γ = 1) and Maxwell molecules (γ = 0) where the collision rate is independent of the velocity
[49, 50, 51, 52].
The central quantity in kinetic theory is the probability P (v,w, t) that a particle has the velocities (v,w) at time
t. We study spatially homogeneous situations where this velocity distribution function is independent of position.
Under the strong assumption that the velocities of the two colliding particles are completely uncorrelated, the velocity
distribution obeys the Boltzmann equation
∂P (v,w)
dvadwadvbdwb |(va − vb) · n̂|γ P (va,wa)P (vb,wb) (9)
δ(v − v′a)δ(w −w′a) + δ(v − v′b)δ(w −w′b)− δ(v − va)δ(w −wa)− δ(v − vb)δ(w −wb)
We integrate over all impact directions with
dn̂ = 1 [53] and over the precollision velocities weighted by the respective
probability distributions. There are two gain terms and two loss terms, because the velocities of interest (v,w) can
be identified with any one of the four velocities in the collision rule (2) and the kernel is simply the collision rate (8).
III. THE LINEAR KINETIC THEORY
The focus of this study is the energy distribution that generally depends only on two variables: Ev and Ew. It
is our aim to compute the distribution P (Ev, Ew) for asymptotically large energies. This will be done for a system
which is driven at very high energies as well as for an undriven system.
As a first step to this goal, we simplify the Boltzmann equation in the limit of large energies. Extremely energetic
particles are rare and as a result it is unlikely that such particles will encounter each other. Hence, energetic particles
typically collide with much slower particles. Since the collision rules are linear, the velocity of the slower particle barely
affects the outcome of the collision. We can therefore neglect the slower velocity. Substituting (va,wa) = (v0,w0)
and (vb,wb) = (0,0) or (va,wa) = (0,0) and (vb,wb) = (v0,w0) into (5) gives the cascade process [54, 55]
(v0,w0) → (v1,w1) + (v2,w2) (10)
where (v0,w0) is the precollision velocity of the energetic particle and (vi,wi) with i = 1, 2 are the consequent
postcollision velocities. With these definitions, the collision rules for extremely energetic particles are
v1 = (1− ηn)v0 · n̂ n̂+ (1− ηt)(v0 − v0 · n̂ n̂)− ηtn̂×w0 w1 =
1− ηt
n̂× v0 (11a)
v2 = ηnv0 · n̂ n̂+ ηt(v0 − v0 · n̂ n̂) + ηtn̂×w0 w2 = −
n̂× v0, (11b)
where we have set R = 1, so that the moment of inertia, I = q, is dimensionless. A collision between a high-energy
particle and a typical-energy particle produces two energetic particles with an energy total that is smaller than the
initial energy. This cascade process transfers energy from large scales to small scales.
Since the cascade process (10) involves only one particle, the tail of the probability distribution P (v,w) obeys the
linear equation
∂P (v,w)
dn̂dv0dw0|v0 · n̂|γP (v0,w0)
δ(v−v1)δ(w−w1)+δ(v−v2)δ(w−w2)−δ(v−v0)δ(w−w0)
. (12)
There are two gain terms and one loss term according to the cascade process (10). Formally, this linear rate equation
can be obtained from the full nonlinear equation (9) by treating either one of the precollision velocities as negligible
and then integrating over this small velocity. This procedure leads to four gain terms and two loss terms and thus,
the factor 1/2 in (9) drops out. We stress that the linear equation (12) is valid only in the high-energy limit.
We also comment that the linear equation (12) for the high-energy tail of the velocity distribution may be valid
in cases where the full nonlinear equation is not. Whereas the nonlinear equation requires that all possible velocities
are uncorrelated, the linear equation merely requires that energetic particles are uncorrelated with typical particles.
This is a much weaker condition.
In this paper, we restrict ourselves to two space dimensions, i.e. rotating disks. In that case the rotational velocities
are always perpendicular to the linear velocities. Thus, we conveniently denote the unit vector in the tangential
direction by t̂ and the unit vector coming out of the plane by ẑ, such that n̂ · t̂ = 0 and n̂× t̂ = ẑ. The precollision
velocities of the energetic particle v0 = vn n̂ + vt t̂ and w0 = w ẑ are compactly written as [vn, vt, w]. With this
notation, the postcollision velocities specified in (12) are
(1− ηn)vn, (1− ηt)vt + ηtw, (ηt/q)vt + (1− ηt/q)w
, and
ηnvn, ηt(vt − w), (ηt/q)(vt − w)
, (13)
respectively. We now treat the three velocity components, namely the normal component of the velocity vn, the
tangential component of the velocity vt, and the scaled angular velocity
qw as a three dimensional vector with
magnitude V0, polar angle θ0, and azimuthal angle φ0:
(vn, vt,
qw) = (V0 sin θ0 cosφ0, V0 sin θ0 sinφ0, V0 cos θ0). (14)
The magnitude V0 gives the energy E0 =
V 20 =
(v2n + v
t + qw
2) while the polar angle characterizes the fraction
of energy stored in the rotational degree of freedom, 1
qw2/E0 = cos
2 θ. In this representation, the postcollision
velocities are three-dimensional vectors with magnitude Vi, polar angle 0 ≤ θi ≤ π, and azimuthal angle 0 ≤ φi ≤ 2π.
The collision rules (11) allow us to express these quantities in terms of V0, θ0, φ0:
(Vi sin θi cosφi, Vi sin θi sinφi, Vi cos θi) = (V0Ai, V0Bi, V0Ci) (15)
where i = 1, 2. The magnitudes of the postcollision velocities are proportional to the magnitude of the precollision
velocity. The three velocity components are scaled by three dimensionless constants Ai, Bi and Ci, that depend on
the angles θ0 and φ0 of the energetic particle, the collision parameters ηn and ηt, and the moment of inertia q,
A1 = (1− ηn) sin θ0 cosφ0 (16a)
B1 = (1− ηt) sin θ0 sinφ0 + (ηt/
q) cos θ0 (16b)
C1 = (ηt/
q) sin θ0 sinφ0 + (1− ηt/q) cos θ0 (16c)
A2 = ηn sin θ0 cosφ0 (16d)
B2 = ηt sin θ0 sinφ0 − (ηt/
q) cos θ0 (16e)
C2 = (ηt/
q) sin θ0 sinφ0 − (ηt/q) cos θ0. (16f)
The new energies are proportional to the precollision energies
Ei = αiE0, with αi = A
i + C
i . (17)
We term the parameters 0 < αi < 1 the contraction parameters. Since the collisions are dissipative, these parameters
satisfy the inequality α1 + α2 ≤ 1. The equality α1 + α2 = 1 holds only for elastic collisions (rn = |rt| = 1). The
energy dissipation is ∆E = E0 − E1 − E2 = ΛE with Λ = 1− α1 − α2 or explicitly,
1− r2n
sin2 θ0 cos
2 φ0 +
1 + q
1− r2t
sin2 θ0 sin
2 φ0 +
cos2 θ0
. (18)
The polar and azimuthal angles are given by
cos θi =
A2i +B
i + C
and tanφi =
, (19)
respectively.
Let us represent solid angles by Ω ≡ cos θ, φ. With this definition, the cascade process (11) is
(E0,Ω0) → (E1,Ω1) + (E2,Ω2) (20)
with Ei and Ωi given by (17) and (19). Energetic particles have an important property: the solid angle is not coupled
to the energy! Indeed, the postcollision angles depend only on the precollision angle. The cascade process has the
following geometric interpretation: a three dimensional vector is duplicated into two vectors. Subsequently, these
two vectors are scaled down by the contraction parameters (17), and rotated according to the angular transformation
(19).
We can now write the linear Boltzmann equation for P (E,Ω), the distribution of energy and solid angle, in a closed
∂P (E,Ω)
dE0dΩ0
E0sin θ0cosφ0
P (E0,Ω0)
δ(E−E1)δ(Ω−Ω1)+δ(E−E2)δ(Ω−Ω2)−δ(E−E0)δ(Ω−Ω0)
Time was rescaled, t → 2γ/2t, to absorb the constant which arises from replacing velocity by energy in the collision
rate (8). Henceforth, we implicitly assume that the distribution P (E,Ω) is independent of φ because the distribution
of linear velocities must be isotropic. The integration over the energy is performed using the collision rule (17), leading
to the linear rate equation for the tail of the energy distribution
∂P (E,Ω)
=Eγ/2
∣ sin θ0 cosφ0
δ(Ω− Ω1)
1+γ/2
δ(Ω− Ω2)
1+γ/2
−P (E,Ω0)δ(Ω− Ω0)
. (21)
This is a non-local equation as the density of particles with energy E is coupled to the density of particles with the
higher energies E/α1 and E/α2. We stress that this equation is a straightforward consequence of the cascade process
(20) and that it can also be derived from the full nonlinear Boltzmann equation. Yet, there may be situations where
the linear equation (21) is valid, while the nonlinear equation (9) is not valid. The bivariate energy distributions
P (E,Ω) and P (Ev, Ew) are completely equivalent but we analyze the former because the cascade process (20) is
transparent in terms of the total energy and the solid angle.
IV. DRIVEN STEADY-STATES
The inelastic Boltzmann equation admits stationary solutions for frictionless particles. These stationary solutions
describe driven steady-states with rare but powerful injection of energy. The injected energy cascades from high-
energies down to small energies, thereby balancing the energy lost in collisions. At energies below the injection scale,
Eqs. (9), (12) and (21) are not altered by the energy source and consequently, the stationary solution of the inelastic
Boltzmann equation holds up to this large energy scale [54, 55]. Here, we seek a corresponding stationary solution
for particles with rotational degrees of freedom in the high energy limit.
The stationary solution has to fulfill Eq. (21) with the left hand side set to zero
∣ sin θ0 cosφ0
1+γ/2
δ(Ω− Ω1)+
1+γ/2
δ(Ω− Ω2)−P (E,Ω0)δ(Ω− Ω0)
. (22)
At high-energies, the solid angle is not coupled to the energy, as follows from Eq. (19). This fact has a ma-
jor consequence: the bivariate energy distribution P (E,Ω) takes the form of a product of the energy distribution
p(E) =
dΩP (E,Ω) and the distribution of solid angle, g(Ω),
P (E,Ω) → p(E) g(Ω) (23)
as E → ∞. The angle distribution is normalized,
dΩ g(Ω) = 1. It does not depend on the azimuthal angle, because
on average the two components of the linear velocity are equivalent. Due to the equi-dimensional (in E) structure of
the steady-state equation (22), the product ansatz (23) is a solution when the distribution p(E) decays algebraically
p(E) ∼ E−ν , (24)
-1 -0.5  0  0.5  1  0
FIG. 1: The exponent ν for hard spheres (γ = 1) as a function of the coefficients of normal, rn, and tangential, rt, restitution
coefficients. The numerical procedure for solving (25) is detailed below.
as E → ∞ [54, 55]. We obtain a closed equation for the distribution g(Ω) by substituting the product ansatz (23)
with the power-law form (24) into the steady-state equation (22)
dΩ0 g(Ω0)
∣ sin θ0 cosφ0
ν−1−γ/2
1 δ(Ω− Ω1) + α
ν−1−γ/2
2 δ(Ω− Ω2)− δ(Ω− Ω0)
. (25)
This equation is linear in g(Ω). However, it is nonlinear in ν and moreover, the solid angles Ωi ≡ Ωi(Ω0) in (19) and
the contraction parameters αi ≡ αi(Ω0) in (17) are complicated functions of the solid angle Ω0.
Equation (25) involves two unknowns quantities, the exponent ν and the distribution function g(Ω). A solution
does not exist for arbitrary values of ν. In fact, there is one and only one value of ν for which there is a solution for
g(Ω). This is the value selected by the cascade dynamics! In other words, (25) is an eigenvalue equation: ν is the
eigenvalue and g(Ω) is the eigenfunction. This eigenvalue equation circumvents the full nonlinear equation (9) and
thus, represents a significant simplification.
The physical interpretation of (25) involves a cascade process in which the solid angle undergoes a creation-
annihilation process
∅ with rate β0,
Ω1 with rate β1,
Ω2 with rate β2.
Here, βi = | sin θ0 cosφ0
αi for i = 0, 1, 2 and α0 = 1. There is one annihilation process and two creation processes.
These processes have relative weights that reflect the powerlaw decay of p(E). At the steady-state, the creation and
the annihilation terms balance (see Appendix B), as reflected in the integrated form of (25)
dΩ0 g(Ω0)
∣ sin θ0 cosφ0
ν−1−γ/2
1 + α
ν−1−γ/2
2 − 1
. (27)
To achieve a steady-state, βi < β0 for i = 1, 2 and therefore α
ν−1−γ/2
i < 1. Since αi < 1, we have the lower bound
ν > 1 + γ/2.
We can immediately check that for elastic collisions, ν = 2 + γ/2 [56, 57] because α1 + α2 = 1, and therefore, we
conclude the bounds 1 + γ/2 ≤ ν ≤ 2 + γ/2. The exponent ν varies continuously with the restitution coefficients
rn and rt and the normalized moment of inertia q. This quantity must coincide with the value found for frictionless
particles where tangential restitution is irrelevant (rt = −1) [54, 55, 58], but otherwise the exponent is distinct, as
shown in Fig. 1. Also, the exponent ν increases monotonically with rn and |rt|. We conclude that the rotational
degrees of motion do affect the power-law behavior (24).
The azimuthal angle θ characterizes the fraction of energy stored in the rotational mode, cos2 θ = Ew/E with
qw2. The angle distribution g(Ω) = (2π)−1f̃(cos θ) therefore captures the partition of energy into rotational
and translational energies. We introduce the natural variable 0 ≤ x ≤ 1 defined by x = | cos θ| so that
and present results for the angle distribution f(x) = 2f̃(cos θ). In equilibrium, energy is partitioned equally into all
degrees of freedom and therefore geq(Ω) = (4π)
−1 or equivalently,
feq(x) = 1 (29)
for 0 ≤ x ≤ 1. In particular, 〈x2〉 = 1/3.
A. Simulation Methods
We numerically studied the angle distribution f(x) by solving the linear eigenvalue equation (25) for the “angular”
process (26) and by solving the full nonlinear Boltzmann equation (9) for the collision process (2). Both of these
equations are solved using Monte Carlo simulations.
The eigenvalue equation is solved by mimicking the angular process. Throughout the simulation, the value ν is
fixed. There are N particles, each with a given polar angle. A particle with polar angle θ0 is picked at random and
then, a random azimuthal angle φ0 is drawn. The polar angles θ1 and θ2 are then calculated according to (19). The
original particle is annihilated with probability β0 and simultaneously, a new particle with angle θ1 is created with
probability β1 and similarly, a second particle with angle θ2 is created with probability β2. Therefore, the number of
particles may increase by one, remain unchanged, or decrease by one. The exponent ν is the value that keeps the total
number of particles constant in the long time limit. The eigenvalue ν is calculated by repeating this simulation for
various values of ν and then using the bisection method [59]. We present Monte Carlo simulations of 100 independent
realizations with N = 107 particles.
Driven steady-states are obtained by simulating the two competing processes of inelastic collisions and energy
injection. In an inelastic collision, two particles are picked at random and also, the impact direction is chosen at
random. The particle velocities are updated according to the collision law (5). Collisions are executed with probability
proportional to the collision rate. Throughout this process, we keep track of the total energy loss. With a small rate,
we augment the energy of a randomly selected particle by an amount equal to the loss total and subsequently, reset
the total energy loss to zero. A fraction of the injected energy is rotational and the complementary fraction is
translational. We draw this fraction according to the equilibrium distribution (29). We experimented with different
angle distributions and found that the resulting stationary state did not change.
Obtaining the distribution f(x) is generally challenging as it requires excellent statistics. The simulations are
most efficient for Maxwell molecules because all possible collisions are equally likely. Therefore, for the full nonlinear
Boltzmann equation (9), we present the angle distribution of the energetic particles only for the case γ = 0.
For Maxwell molecules, the injection rate is 10−4 and the system size is N = 107. The corresponding values for
hard spheres are 10−2 and N = 105. In all cases, the simulation results represent an average over 102 independent
realizations. Unless noted otherwise, the simulation results are for maximally dissipative (rn = rt = 0) disks (q = 1/2).
B. The Distribution of Total Energy
The numerical simulations confirm several of our theoretical predictions. First, the energy distribution approaches a
steady-state with a power-law high-energy tail. Second, the distribution of the total energy p(E) decays algebraically
as in (24). Third, the exponent ν is in excellent agreement with the predictions of the eigenvalue equation. For
Maxwell molecules, Monte Carlo simulation of the full nonlinear equation yields ν = 1.570± 0.005 whereas numerical
solution of the eigenvalue equation (25) gives ν = 1.569 ± 0.005 (Fig. 2). For hard-spheres, where the simulation
results are slightly less accurate, the corresponding values are ν = 2.065 ± 0.005 and ν = 2.060 ± 0.005 (Fig. 3).
The behavior of the distribution of total energy is therefore qualitatively similar to the behavior in the no-rotation
case [54, 55]. However, the quantitative behavior is different because the exponent ν does depend on the tangential
restitution coefficient and the moment of inertia (Fig. 1).
C. The Angle distribution
The numerical simulations also confirm several of our theoretical predictions concerning the angle distribution.
Extremely energetic particles have a universal distribution f(x). This distribution is independent of the energy,
theory
simulation
FIG. 2: The tail of the energy distribution for driven
Maxwell molecules. Shown are simulation results (solid line)
and a line with the slope predicted by the theory (dashed
line). The energy is normalized by the typical energy 10−4.
theory
simulation
FIG. 3: The tail of the energy distribution for driven hard
spheres. Shown are simulation results (solid line) and a line
with the slope predicted by the theory (dashed line).
0 0.2 0.4 0.6 0.8 1
Angular Process
Collision Process
FIG. 4: The angle distribution f(x) obtained by Monte
Carlo simulation of the angular process (26) (solid line) and
the collision process (2) (dashed line) for Maxwell molecules.
The special values x1, x2, and x3 discusses in the text are
indicated by arrows.
0 0.2 0.4 0.6 0.8 1
=0, r
=-0.9
=0.9, r
=0.9, r
FIG. 5: The angle distribution f(x) for various collision
parameters (rn and rt) for Maxwell molecules.
provided that the energy is sufficiently large. We had to probe only the most energetic particle out of roughly 103
particles to measure this distribution. For this reason, the linear analysis and the resulting eigenvalue equation are
valuable because they allow for an accurate and efficient determination of the angle distribution of the energetic
particles. We also verified that the distribution f(x) obeys the eigenvalue equation (25), as demonstrated in Fig. 4,
where the simulations are compared to the solution of the angular process.
The distribution f(x) has several noteworthy features. First, it is not uniform, implying the breakdown of en-
ergy equipartition in a granular gas. Furthermore, this distribution is nonanalytic. It contains singularities and
discontinuous derivatives. There are notable peaks in the distribution so that special values x and special ratios
Ew/E are strongly preferred. The reason for these peaks is the fact that the polar angle is limited. For example,
cos2 θ2 < 1/(1 + q) as seen by substituting cos θ0 = ±1 into (16) and (19). Consequently, there is a special ratio
1 + q
with the corresponding special energy ratio Ew/E = x
1. This is the most pronounced peak in Fig. 4,
0 0.2 0.4 0.6 0.8 1
FIG. 6: The angle distribution f(x) for hard spheres.
0 0.2 0.4 0.6 0.8 1
FIG. 7: The angle distribution fall(x) of all particles for
Maxwell molecules (solid line). Also shown for reference is
the uniform equilibrium distribution (broken line).
2/3 = 0.81649. Numerically, we observe that the peak becomes more pronounced as the distribution is mea-
sured at a finer scale, indicating that the distribution function diverges at this point.
Similarly, there is another special ratio that corresponds to θ1 when cos θ0 = ±1, and unlike (30), this location
depends on the tangential restitution,
1− ηt/q
η2t /q + (1− ηt/q)2
. (31)
Indeed, there is a barely noticeable cusp at x2 =
8/9 = 0.942809. Singularities may induce less pronounced
“echo”-singularities. For example, using cos θ0 = x1 and φ0 = π/2 yields the special ratio
1 + ηt(1− 1/q)
q[1− ηt(1− 1/q)]2 + [1 + ηt(1− 1/q)]2
. (32)
There is a noticeable peak at the corresponding value x3 =
50/99 = 0.710669 in Fig. 4. We anticipate that as
the transformation (19) is iterated, the strength of the singularities weakens and as a result there are discontinuous
derivatives of increasing order, a subtle behavior that is difficult to measure.
The location of the singularities varies with the collision parameters rn and rt and the moment of inertia q. In
fact, the angle distribution is extremely sensitive to material properties as its shape changes dramatically with these
parameters, see Fig. 5. The angle distribution also depends on the collision rate and it is much smoother for hard
spheres, see Fig. 6. Since the collision rate vanishes for grazing collisions, φ = π/2, the associated singularities
including in particular (32) are suppressed. Nevertheless, there is a pronounced jump at the special ratio given by
(30) and there are also noticeable cusps.
The angle distribution of all particles fall(x) ∝
dE P (E,Ω) is shown in Fig. 7. It is substantially different from
f(x). Therefore, the energy distribution P (E,Ω) does not factorize in general and there are correlations between the
solid angle and the total energy. Only for energetic particles does (23) hold. Moreover, fall(x) is much smoother in
comparison with f(x) although there is a jump in the first derivative at the special ratio (30) showing that the angle
distribution of all particles is also non-analytic, see Fig. 7. Generally, the angle distribution depends on energy and
the deviation from a uniform distribution grows with energy.
We also comment that lone measurement of the moment 〈x2〉 can be misleading. The angle distribution may very
well have a value close to the equipartition value 〈x2〉eq = 1/3 but still, be very far from the equilibrium distribution.
Indeed, in Fig. 4, 〈x2〉 ∼= 0.318, a value that barely differs from the equilibrium value, even though the corresponding
distribution is far from uniform. The second moment may also differ substantially from the equipartition value and
for example, 〈x2〉 = 0.202 when rn = 0.9 and rt = 0 (Fig. 5).
We argue that the qualitative features of the angle distribution should be generic in granular materials. Colli-
sions involving energetic particles must follow the linear cascade rules (20) with the angular transformations (19).
The singularities are a direct consequence of these transformations and therefore should be generic. Measuring the
parameter-sensitive distribution f(x) experimentally is challenging because a huge number of particles must be probed
and the measurement has to be accurate. The distribution fall(x) provides a detailed probe of the partition of energy
into rotational and translation motion.
V. FREE COOLING
We now consider freely cooling granular gases that evolve via purely collisional dynamics. Without energy input,
all energy is eventually dissipated and the particles come to rest. This system has been studied extensively [1] for
hard spheres with [33, 47] and without rotation [60].
We consider Maxwell molecules where in the absence of rotation an exact treatment is possible [24, 25, 27, 61, 62].
When γ = 0 the Boltzmann equation (9) simplifies
∂P (v,w)
dvadwadvbdwb P (va,wa)P (vb,wb) (33)
× [δ(v − v′a)δ(w −w′a) + δ(v − v′b)δ(w −w′b)− δ(v − va)δ(w −wa)− δ(v − vb)δ(w −wb)].
Consequently, the equations for the moments 〈vnwm〉 =
dvdwP (v,w)vnwm close.
A. The Temperatures
Here, we consider only the translational temperature defined as the average translational energy, Tv = 〈Ev〉, and
the rotational temperature, defined as the average rotational energy Tw = 〈Ew〉. These two temperatures are coupled
through the linear equation
λvv λvw
λwv λww
. (34)
Appendix C details the derivation of the matrix of coefficients
λvv = ηn(1− ηn) + ηt(1 − ηt) (35a)
λvw = −2η2t /q (35b)
λwv = −η2t /q (35c)
λww = 2(ηt/q)(1− ηt/q). (35d)
The two temperatures are coupled as long as ηt 6= 0 or alternatively, rt 6= −1.
The solution of (34) is a linear combination of the two eigenvectors
e−λ−t + C+
e−λ+t (36)
with the constants C− and C+ set by the initial conditions, and c± = (λ± − λvv)/λvw. The eigenvalues are
λvv + λww
λvv − λww
+ λvwλwv . (37)
The larger eigenvalue is irrelevant in the long time limit and therefore,
e−λt (38)
such that both temperatures decay with the same rate λ ≡ λ−. Of course, the total temperature also follows the
same exponential decay, T = Tv + Tw ∼ e−λt. In this regime, the fraction of rotational energy is on average
1 + c−
λ− λvv
λ+ λvw − λvv
. (39)
The approach toward this value is exponentially fast and the relaxation time is inversely proportional to the difference
in eigenvalues τ = 1/(λ+ − λ−).
In equilibrium, Tw/T = 1/3 but for nonequilibrium granular gases the ratio varies. In Fig. 8 we plot the ratio of
the average rotational energy to the total energy as a function of the coefficients of restitution. In accordance with
our findings for driven steady-states, energy is not partitioned equally between all the degrees of freedom.
 1  0
FIG. 8: The ratio of average rotational energy to total energy
as a function of the coefficients of normal, rn, and tangential,
rt, restitution.
theory
simulation
FIG. 9: The scaling function underlying the energy distri-
bution (solid line). The distribution was obtained using a
Monte Carlo simulation with N = 107 particles. A dashed
line with the slope predicted by the theory is also shown for
reference.
B. The Energy Distribution
To study the full energy distribution, it is again convenient to make a transformation of variables from the velocity
pair (v,w) to the total energy and the solid angle (E,Ω). The energy distribution is now time dependent and assuming
that the temperature – T ∼ e−λt – is the characteristic energy scale we postulate the self-similar form
P (E,Ω, t) → eλtΦ(Eeλt,Ω) (40)
with the prefactor ensuring proper normalization,
dz dΩΦ(z,Ω) = 1. We focus on the high-energy behavior where
the linear equation (21) holds. By substituting the scaling form (40) into this linear equation and setting γ = 0, we
find the integro-differential equation governing the scaling function
λΦ(z,Ω) + λz
Φ(z,Ω)=
δ(Ω− Ω1) +
δ(Ω− Ω2)− Φ(z,Ω0)δ(Ω− Ω0)
.(41)
We again write the multivariate energy distribution as a product Φ(z,Ω) → ψ(z)g(Ω) of the distribution of the total
energy ψ(z) =
dΩΦ(z,Ω) and the distribution of the solid angle g(Ω). This form is a solution of the equi-dimensional
equation (41) when the distribution of the total energy decays as a power-law
ψ(z) ∼ z−ν (42)
at large energies, z → ∞. The angle distribution satisfies the eigenvalue equation
dΩ0 g(Ω0)
αν−11 δ(Ω− Ω1) + α
2 δ(Ω− Ω2)− [1− λ(ν − 1)]δ(Ω− Ω0)
. (43)
Of course, setting λ = 0, one recovers the steady-state equation (25) reflecting that the similarity solution is stationary.
The factor 1 is replaced by the smaller factor 1 − λ(ν − 1) that accounts for the constant decrease in the number
of particles at any given energy because of dissipation. Again, we have a nonlinear eigenvalue equation with the
eigenvalue ν and the eigenfunction g(Ω).
We solve this eigenvalue equation by performing a Monte Carlo simulation of the same angular process as described
by (26) but with a different annihilation rate β0 = 1−λ(ν − 1). We compare the angle distribution predicted by (43)
with the behavior of the energetic particles in the freely cooling gas.
The numerical simulations of the inelastic collision process confirm the theoretical predictions. First, the energy
distribution is self-similar as in (40) and the characteristic scale is proportional to the temperature. Second, the
distribution of the total energy has a power-law tail, as displayed in Fig. 9 and the exponent ν is very close to the
theoretical prediction (numerical simulations of the collision process gives ν = 2.98±0.05 while the eigenvalue equation
yields ν = 2.92± 0.05).
0 0.2 0.4 0.6 0.8 1
Angular Process
Collision Process
FIG. 10: The angle distribution of the energetic particles.
Shown are results for the collision process (solid line) and
for the angular process (dashed line).
0 0.2 0.4 0.6 0.8 1
FIG. 11: The angle distribution of all particles for a freely
cooling gas (solid line). Also shown for reference is the uni-
form equilibrium distribution.
The angle distribution deviates even more strongly from the uniform distribution with a very pronounced peak
(see Fig. 10) because the dynamics are purely collisional. The singularities are weaker although the one at x1 given
by (30) is clear. The agreement between the solution of the angular process and the Monte Carlo simulations is
slightly worse than for driven systems because the statistics become prohibitive: now it is necessary to probe the most
energetic out of roughly 106 particles to obtain the asymptotic angle distribution! The sharper power-law decay is
responsible for this three order of magnitude increase: the cumulative distribution of total energy decays according
dE′p(E′) ∼ E−µ with µ = ν − 1 about three times larger than before. Finally, the angle distribution of all
particles deviates only slightly from a uniform distribution (see Fig. 11). We conclude that the behavior of the freely
cooling gas is qualitatively similar to that found in driven steady-states.
VI. CONCLUSIONS AND OUTLOOK
The complete description of granular media with translational and rotational degrees of freedom requires the full
bivariate distribution of energies. It is not sufficient to consider only the average kinetic energy of translations and
rotations. Instead the full bivariate distribution is highly nontrivial. We have shown that in the limit of large particle
energy, this distribution obeys a linear equation. Its solution can be written as a product of two distributions, one
for the total energy, E = Ev + Ew, and one for the variable x =
Ew/E, which captures the partition of the
total energy between rotational and translational motion. The distribution of the total energy decays algebraically
and the characteristic exponent depends on the collision parameters and the moment of inertia. The variable x is
not uniformly distributed as in equilibrium. Instead the distribution f(x) is not analytic and displays a series of
singularities of varying strengths. Remarkably, there are special preferred ratios of rotational-to-total energy. This
violation of energy equipartition among different degrees of freedom is a direct consequence of the energy dissipation.
The total energy and the variable x are correlated in general with the deviations from equilibrium increasing with
energy. These two variable become uncorrelated only at extremely high-energies.
We have studied both, the system which is driven at extremely high energies and displays a stationary energy
cascade on energy scales below the driving one, and a freely cooling gas. In the latter gas the bivariate energy
distribution is time dependent, reflecting the overall decrease of energy. Nevertheless, scaling the total energy with
temperature, one finds a self-similar form for the distribution, which again factorizes in the high-energy limit. As in
the driven system, the distribution of the total energy decays as a power law with, however, different exponents for
the driven and the free cooling system. The angular distribution deviates even more from the uniform (equipartition)
one in the cooling system.
It should be straightforward to extend these results to three dimensions where the angular process takes place in
three dimensions. In the limit of high energies one would again expect a limiting distribution for the partition angle
Ew/E. Another possible extension refers to a more realistic law of friction, including Coulomb friction [39, 40].
Finally, it would be of interest to extend the analysis to other systems, where equipartition is violated. An example
is a binary mixture, where the energy is shared unequally between the two components.
Acknowledgments
We thank the Kavli Institute for Theoretical Physics in University of California, Santa Barbara where this work
was initiated. We acknowledge financial support from DOE grant DE-AC52-06NA25396.
[1] Kinetic theory of granular gases, N. Brilliantov and T. Pöschel, (Oxford, Oxford, 2003).
[2] Granular Gases, T. Pöschel and S. Luding (editors), (Springer, Berlin, 2000).
[3] Granular Gas Dynamics, T. Pöschel and N. Brilliantov (editors), (Springer, Berlin, 2003).
[4] S. McNamara and W. R. Young, Phys Fluids A 4, 496 (1992).
[5] J. S. Olafsen and J. S. Urbach, Phys. Rev. Lett. 81, 4369 (1998).
[6] S. Luding and H. J. Herrmann, Chaos 9, 673 (1999).
[7] X. Nie, E. Ben-Naim, and S. Y. Chen, Phys. Rev. Lett. 89, 204301 (2002).
[8] D. van der Meer, K. van der Weele, and D. Lohse, Phys. Rev. Lett 88, 174302 (2002).
[9] E. Ben-Naim, S. Y. Chen, G. D. Doolen, and S. Redner, Phys. Rev. Lett. 83, 4069 (1999).
[10] E. Efrati, E. Livne, and B. Meerson Phys. Rev. Lett. 94, 088001 (2005)
[11] V. Yu. Zaburdaev, M. Brinkmann, and S. Herminghaus Phys. Rev. Lett. 97, 018001 (2006).
[12] E. C. Rericha, C. Bizon, M. D. Shattuck, and H. L. Swinney, Phys. Rev. Lett. 88, 014302 (2002).
[13] A. Samadani, L. Mahadevan, and A. Kudrolli, J. Fluid Mech. 452, 293 (2002).
[14] I. Goldhirsch, and G. Zanetti, Phys. Rev. Lett. 70, 1619 (1993).
[15] E. Khain and B. Meerson, Europhys. Lett. 65, 193 (2004).
[16] W. Losert, D. G. W. Cooper, J. Delour, A. Kudrolli, and J. P. Gollub, Chaos 9, 682 (1999).
[17] F. Rouyer and N. Menon, Phys. Rev. Lett. 85, 3676 (2000)
[18] I. S. Aranson and J. S. Olafsen Phys. Rev. E 66, 061302 (2002).
[19] Y. Du, H. Li, and L. P. Kadanoff, Phys. Rev. Lett. 74, 1268 (1995).
[20] A. Kudrolli, M. Wolpert and J. P. Gollub, Phys. Rev. Lett. 78, 1383 (1997).
[21] E. L. Grossman, T. Zhou, and E. Ben-Naim, Phys. Rev. E 55, 4200 (1997).
[22] R. D. Wildman and D. J. Parker, Phys. Rev. Lett. 88, 064301 (2002).
[23] K. Feitosa and N. Menon, Phys. Rev. Lett. 88, 198301 (2002).
[24] P. L. Krapivsky and E. Ben-Naim, J. Phys. A 58, 182 (2002).
[25] M. H. Ernst and R. Brito, Europhys. Lett. 58, 182 (2002).
[26] A. Baldassarri, U. M. B. Marconi, and A. Puglisi, Europhys. Lett. 58, 14 (2002).
[27] E. Ben-Naim and P. L. Krapivsky, Phys. Rev. E 61, R5 (2000).
[28] T. P. C. van Noije and M. H. Ernst, Gran. Matt. 1, 57 (1998).
[29] E. Ben-Naim and P. L. Krapivsky, Phys. Rev. E 66, 011309 (2002).
[30] K. Kohlstedt, A. Snezhko, M. V. Sapoznikov, I. S. Aranson, J. S. Olafsen, and E. Ben-Naim, Phys. Rev. Lett. 95, 068001
(2005).
[31] N. Schorghofer and T. Zhou, Phys. Rev. E 54, 5511 (1996).
[32] I. Goldhirsch and M. L. Tan, Phys. Fluids 7, 1752 (1996).
[33] M. Huthmann and A. Zippelius, Phys. Rev. E 56, R6275 (1997).
[34] S. Luding, M. Huthmann, S. McNamara and A. Zippelius, Phys. Rev. E 58, 3416 (1998).
[35] O. Herbst, M. Huthmann, and A. Zippelius, Granular Matter 2, 211 (2000).
[36] J. T. Jenkins and C. Zhang, Physics of Fluids 14, 1228 (2002).
[37] O. Herbst, R. Cafiero, A. Zippelius, H. J. Herrmann and S. Luding, Physics of Fluids 17, 107102 (2005).
[38] M. Huthmann, J. Orza, and R. Brito, Granular matter 2, 189 (2000).
[39] O. R. Walton, In ”Particle Two-Phase Flow”, ed. M. C. Rocco, Butterworth, London 1993, p884.
[40] S. F. Foerster, M. Y. Louge, H. Chang, and K. Allia, Physics of Fluids 6, 1108 (1994).
[41] J. C. Tsai, F. Ye, J. P. Gollub, and T. C. Lubensky, Phys. Rev. Lett. 94, 214301 (2005).
[42] J. T. Jenkins and M. W. Richman, Phys. Fluids 28, 3485 (1985).
[43] J. J. Brey, J. W. Dufty, C. S. Kim, and A. Santos, Phys. Rev. E 58, 4638 (1998).
[44] J. Lutsko, J. J. Brey, and J. W. Dufty, Phys. Rev. E 65, 051304 (2002).
[45] I. Goldhirsch, Ann. Rev. Fluid. Mech. 35, 267 (2003).
[46] I. Goldhirsch, S. H. Noskowicz, and O. Bar-Lev, Phys. Rev. Lett. 95, 068002 (2005).
[47] N. V. Brilliantov, T. Pöschel, W. T. Kranz, and A. Zippelius, Phys. Rev. Lett. 98, 128001 (2007).
[48] P. Résibois and M. de Leener, Classical Kinetic Theory of Fluids (John Wiley, New York, 1977).
[49] J. C. Maxwell, Phil. Trans. R. Soc. 157, 49 (1867).
[50] C. Truesdell and R. G. Muncaster, Fundamentals of Maxwell’s Kinetic Theory of a Simple Monoatomic Gas (Academic
Press, New York, 1980).
[51] R. S. Krupp, A nonequilibrium solution of the Fourier transformed Boltzmann equation, M.S. Thesis, MIT (1967); Inves-
tigation of solutions to the Fourier transformed Boltzmann equation, Ph.D. Thesis, MIT (1970).
[52] M. H. Ernst, Phys. Reports 78, 1 (1981).
[53] We tacitly ignore the condition (va−vb)·n̂ < 0 because the collision rules are invariant under exchange of the two particles.
[54] E. Ben-Naim and J. Machta, Phys. Rev. Lett. 94, 138001 (2005).
[55] E. Ben-Naim, B. Machta, and J. Machta, Phys. Rev. E 72, 021302 (2005).
[56] A. V. Bobylev, Sov. Sci. Rev. C. Math. Phys. 7, 111 (1988).
[57] L. Acedo, A. Santos, and A. V. Bobylev, J. Stat. Phys. 109, 1027 (2002).
[58] The high-energy behavior (24) is equivalent to the large-velocity tail P (v) ∼ v−σ with σ = 2(ν − 1) + d.
[59] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes, (Cambridge University Press,
Cambridge, 1992).
[60] S. E. Esipov and T. Pöschel, J. Stat. Phys. 86, 1385 (1997).
[61] A. V. Bobylev, J. A. Carrillo, and I. M. Gamba, J. Stat. Phys. 98, 743 (2000).
[62] A. V. Bobylev and C. Cercignani, J. Stat. Phys. 106, 547 (2002).
APPENDIX A: THE COLLISION RULES
The total linear momentum v′a + v
b = va + vb is conserved in the collision. The angular momenta of the two
particles with respect to the point of contact, ω′i, are given by
Iωa = Iwa +mR n̂× va (A1a)
Iωb = Iwb −mR n̂× vb. (A1b)
These are conserved, ω′i = ωi with i = a, b, because there is no torque at the point of contact. In inelastic collisions,
the normal and tangential components of the relative velocity at the point of contact obey the collision law (4) where
U = V +R n̂×W.
It is convenient to introduce the momentum transfer δ, defined as follows: v′a = va−δ and v′b = vb+δ. Conservation
of the angular velocity with respect to the point of contact and Eq. (A1) gives w′i = wi +
n̂ × δ. In terms of δ,
the difference in velocity at the point of contact is U′ = U− 2δ + 2
n̂× n̂× δ. Substituting this expression into the
collision laws (4), the normal and the tangential components of δ are simply
δ · n̂ = ηn U · n̂ (A2a)
δ × n̂ = ηt U× n̂. (A2b)
Consequently, the momentum transfer is δ = ηnU · n̂ n̂+ ηt(U−U · n̂ n̂) or explicitly,
δ = ηn V · n̂ n̂+ ηt (V −V · n̂ n̂) + ηt R n̂×W (A3)
We now have the explicit collision rules (5).
APPENDIX B: PARTICLE NUMBER CONSERVATION
In this appendix, we verify that the stationary solution is consistent with particle number conservation. Maxwell
Molecules are considered for simplicity. It is straightforward to generalize this calculation to all γ and to free cooling.
Our starting point is Eq. (21), specialized to Maxwell molecules, i.e. γ = 0,
∂P (E,Ω)
δ(Ω− Ω1) +
δ(Ω− Ω2)− P (E,Ω0)δ(Ω− Ω0)
. (B1)
As a first step we integrate this equation over the solid angle
∂p(E)
− P (E,Ω0)
. (B2)
The power-law behavior (24) typically holds in a restricted energy range, El ≤ E ≤ Eu, where El and Eu are upper
and lower cutoffs. In the driven case, the upper cutoff is set by the energy injection scale. Let N =
dE p(E) be
the total number of particles in this range. With the powerlaw decay (24), then
N ∼ 1
ν − 1
E1−νl − E
. (B3)
To evaluate this time evolution of N , we substitute the product form (23) into (B2) and integrate over the energies
in the aforementioned power-law range,
dΩ0 g(Ω0)
αν−11 + α
2 − 1
. (B4)
Using Eq. (27), we confirm that the total number of particles is conserved, ∂N/∂t = 0.
APPENDIX C: THE MATRIX COEFFICIENTS
In an inelastic collision, the translational energy loss is ∆Ev = Ev − E′v with Ev = 12 (v
a + v
b ) and similarly, the
rotational energy loss is ∆Ew = Ew −E′w with Ew = 12q(w
b ). We can conveniently calculate these quantities by
using v′a = va − δ, v′b = vb + δ, and w′i = wi + (1/qR)n̂× δ, and by expressing the momentum transfer δ using the
natural coordinate system, δ = ηn Vnn̂+ ηt(Vt −W )t̂,
∆Ev = ηn(1− ηn)V 2n + ηt(1− ηt)V 2t − η2tW 2 + ηt(2ηt − 1)VtW (C1a)
∆Ew = −(η2t /q)V 2t + ηt(1− ηt/q)W 2 − ηt(1− 2ηt/q)VtW. (C1b)
The rate of change of the respective temperatures equals 1/2 the average of this quantities, d
〈∆Ev〉 and
〈∆Ew〉. This is seen by multiplying (33) by 12v
2 and by integrating over the velocity. The averaging is with
respect to the probability distribution functions of the two colliding particles. The cross-term vanishes, 〈VtW 〉 = 0,
by symmetry. Using 〈V 2n 〉 = 2〈v2n〉 = 〈v2〉 = 2Tv and 〈W 2〉 = 2〈w2〉 = 4Tw/q we obtain the matrix elements (35).
	Introduction
	The Nonlinear Kinetic Theory
	The Linear Kinetic Theory
	Driven Steady-states
	Simulation Methods
	The Distribution of Total Energy
	The Angle distribution
	Free Cooling
	The Temperatures
	The Energy Distribution
	Conclusions and Outlook
	Acknowledgments
	References
	The Collision rules
	Particle Number Conservation
	The matrix coefficients

  We study the kinetic theory of driven granular gases, taking into account
both translational and rotational degrees of freedom. We obtain the high-energy
tail of the stationary bivariate energy distribution, depending on the total
energy E and the ratio x=sqrt{E_w/E} of rotational energy E_w to total energy.
Extremely energetic particles have a unique and well-defined distribution f(x)
which has several remarkable features: x is not uniformly distributed as in
molecular gases; f(x) is not smooth but has multiple singularities. The latter
behavior is sensitive to material properties such as the collision parameters,
the moment of inertia and the collision rate. Interestingly, there are
preferred ratios of rotational-to-total energy. In general, f(x) is strongly
correlated with energy and the deviations from a uniform distribution grow with
energy. We also solve for the energy distribution of freely cooling Maxwell
Molecules and find qualitatively similar behavior.

Introduction
	The Nonlinear Kinetic Theory
	The Linear Kinetic Theory
	Driven Steady-states
	Simulation Methods
	The Distribution of Total Energy
	The Angle distribution
	Free Cooling
	The Temperatures
	The Energy Distribution
	Conclusions and Outlook
	Acknowledgments
	References
	The Collision rules
	Particle Number Conservation
	The matrix coefficients
