IOPscience

Journal of Physics A: Mathematical and General

Quick search Find article
Quick search
Find article

J. Phys. A: Math. Gen. 35 No 11 (23 March 2002) L147-L152
PII: S0305-4470(02)33044-0

LETTER TO THE EDITOR

Nontrivial velocity distributions in inelastic gases

P L Krapivsky1 and E Ben-Naim2

1 Center for Polymer Studies and Department of Physics, Boston University, Boston, MA 02215,USA
2 Theoretical Division and Center for Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, NM 87545, USA

Received 23 January 2002
Published 8 March 2002

Abstract. We study spatially homogeneous inelastic gases using the Boltzmann equation. We consider uniform collision rates and obtain analytical results valid for arbitrary spatial dimension d and arbitrary dissipation coefficient $\epsilon$. In the unforced case, we find that the velocity distribution decays algebraically, $P(v,t)\sim v^{-\sigma}$, for sufficiently large velocities. The exponent $\sigma(d,\epsilon)$ exhibits nontrivial dependence on the spatial dimension and the dissipation coefficient.

PACS Numbers: 05.20. Dd, 05.40.-a, 45.70. Mg, 47.70. Nd

Granular gases present novel challenges, previously not encountered in fluid dynamics [1]. Specifically, the strong underlying energy dissipation leads to clustering instabilities and strong velocity correlations [2,3,4,5,6,7]. A series of recent experimental and theoretical studies reveals a rich phenomenology. In particular, velocities are characterized by anomalous statistics, sensitive to the details of the driving conditions, the density and the degree of dissipation [8,9,10,11,12].

Kinetic theory provides a systematic framework for deriving macroscopic properties from the microscopic collision dynamics [13,14,15]. Yet, analysis of the corresponding Boltzmann equation often involves uncontrolled approximations or use of nearly Maxwellian distributions. Motivated by the latter issue, we examine both unforced and forced inelastic gases using a simplified Boltzmann equation. Specifically, we employ Maxwell's collision rate, which is proportional to the typical velocity rather than the relative velocity [16]. The resulting Boltzmann equation is analytically tractable, as reported in a few recent studies [17,18,19].

This kinetic theory leads to interesting behaviours in the freely evolving case. In one dimension, while moments of the velocity distribution exhibit multiscaling [17], the velocity distribution itself still approaches a scaling form with an algebraic large-velocity tail [18]. An algebraic tail was also found numerically in two dimensions [18]. Here, we show analytically that in an arbitrary spatial dimension the velocity distribution admits a scaling solution with an algebraic large-velocity tail. The corresponding exponent, a root of a transcendental equation, depends on the spatial dimension and the restitution coefficient. In the driven case, we find that the velocity distribution is non-Maxwellian.

Our starting point is a homogeneous system of identical inelastic spherical particles in arbitrary spatial dimension d. The mass and the cross-section are set to unity without loss of generality. When two particles collide, the normal component of the relative velocity is reduced by the restitution coefficient $r=1-2\epsilon$, while the tangential component remains the same. Denoting by n the impact direction, the unit vector connecting the centres of the two particles, the post-collision velocities v1,2 are given by a linear combination of the precollision velocities u1,2,

 \begin{equation} \bi{v}_{1,2}=\bi{u}_{1,2}\mp(1-\epsilon)~ (\bi{g}\bdot \bi{n})~\bi{n}, \end{equation} (1)

with the relative velocity g = u1 - u2. The energy dissipated in a collision equals $\Delta E=-\epsilon(1-\epsilon)(\bi{g}\bdot\bi{n})^2$. When $\epsilon=0$ collisions are elastic, while when $\epsilon=1/2$collisions are perfectly inelastic with maximal energy dissipation.

Let P(v,t) be the normalized density of particles with velocity v at time t. Assuming molecular chaos or perfect mixing, we arrive at the following Boltzmann equation for the velocity distribution:

 \begin{equation} \eqalign{ \fl {\partial P(\bi{v},t)\over \partial t} = \langle... ...psilon) (\bi{g}\bdot \bi{n})~\bi{n}] -\delta(\bi{v}-\bi{u}_1)\}. } \end{equation} (2)

The collision rate in the Maxwell approximation represents the typical velocity scale, $\langle g\rangle=\sqrt{T}$, with the `granular temperature' or the average kinetic energy per degree of freedom (for hard spheres, the collision rate equals the actual relative velocity $\bi{g}\bdot \bi{n}$). This evolution equation naturally describes a stochastic process where randomly chosen pairs of particles undergo inelastic collisions according to (1) with a randomly chosen impact direction n. In writing equation (2) we tacitly ignored the restriction $\bi{g}~\bdot~ \bi{n}\gt$ on the integration range, because the integrand obeys the reflection symmetry $\bi{n}\to -\bi{n}$. The integration measure should be normalized, $\int \mathrm{d}\bi{n}=1$.

In the absence of energy input the system `cools' indefinitely according to Haff's law [20]. Indeed, the time dependence of the temperature is found from the Boltzmann equation (2) to give $({\mathrm{d}/ \mathrm{d} t})T=-\lambda T^{3/2}$ with $\lambda=2\epsilon(1-\epsilon)\int \mathrm{d}\bi{n}~ n_1^2$ (the first axis was conveniently chosen to be parallel to g). Since $n_1^2+\cdots+n_d^2=1$ and $\int \mathrm{d}\bi{n}=1$ one has $\lambda=2\epsilon(1-\epsilon)/d$. Thus, the temperature decays according to T(t) = T0 (1+t/t0)-2, with the initial temperature T0 and the characteristic timescale $t_0=d/[\epsilon(1-\epsilon)\sqrt{T_0}]$. The temperature quantifies velocity fluctuations. In the following, we focus on the natural case of isotropic velocity distributions. We shall show that asymptotically, the temperature represents the only relevant velocity scale as the velocity distribution approaches the scaling form

 \begin{equation} P(\bi{v},t)\sim {1\over T^{d/2}}{\cal P}\biggl({v\over \sqrt{T}}\biggr). \end{equation} (3)

Given the convolution structure of the Boltzmann equation (2), we introduce F(k, t), the Fourier transform of the velocity distribution function, $F(\bi{k}, t)=\int \mathrm{d}\bi{v}~\mathrm{e}^{\mathrm{i}\bi{k}\bdot \bi{v}}~P(\bi{v},t)$. Applying the Fourier transform to equation (2) and integrating over the velocities gives

 \begin{equation} {1\over \sqrt{T}}{\partial \over \partial t} F(\bi{k},t)+F(\bi{k},t)= \int \mathrm{d}\bi{n}~F[\bi{k}-\bi{q},t] F[\bi{q},t], \end{equation} (4)

with $\bi{q}=(1-\epsilon)(\bi{k}\bdot \bi{n})~\bi{n}$. This rate equation reflects the momentum transferred between particles during collisions. We seek an isotropic scaling solution for the Fourier transform, the equivalent of (3),

 \begin{equation} F(\bi{k},t)=\Phi(k^2 T), \end{equation} (5)

with $k\equiv \vert\bi{k}\vert$. In the $k\to 0$ limit, the Fourier transform, $F(\bi{k},t)\cong 1-({1/ 2})~k^2~ T$, implies that the first two terms in the Taylor expansion of the corresponding scaling function $\Phi(x)$ are universal, $\Phi(x)\cong 1-({1/ 2})x$.

Let us first consider the simpler one-dimensional case. Equation (4) reduces to $({1/ \sqrt{T}})({\partial / \partial t})F(k,t)+F(k,t)= F[\epsilon k, t]~F[k-\epsilon k, t]$ and the scaling function (5) satisfies $-\lambda x\Phi'(x)+\Phi(x) =\Phi[\epsilon^2 x]~\Phi[(1-\epsilon)^2 x]~$[17]. This equation admits a very simple solution [18]

 \begin{equation} \Phi(x)=[1+\sqrt{x}]~\mathrm{e}^{-\sqrt{x}}. \end{equation} (6)

The inverse Fourier transform gives the scaling function of the velocity distribution as a squared Lorentzian

 \begin{equation} {\cal P}(w)={2\over \pi}{1\over (1+w^2)^2}. \end{equation} (7)

Therefore, the form of the scaling solution (7) is universal as it is independent of the dissipation degree. Another important feature is the algebraic tail of the velocity distribution, ${\cal P}(w)\sim w^{-4}$ as $w\to\infty$.

We now return to the general d-dimensional case. Substituting (5) into (4) and using the temperature cooling equation $({\mathrm{d}/ \mathrm{d} t})T=-\lambda T^{3/2}$ one finds that the scaling function $\Phi(x)$ satisfies $\Phi(x)-\lambda x\Phi'(x)=\int\! \mathrm{d}\bi{n}~\Phi(\xi x)\Phi(\eta x)$, where $\xi=1-(1-\epsilon^2)n_1^2$ and $\eta=(1-\epsilon)^2n_1^2$. To integrate over the impact direction n we use spherical coordinates and treat the first axis as the polar axis, $n_1=\cos \theta$. The $\theta$-dependent factor of the measure is $\mathrm{d}\bi{n}=N^{-1}(\sin \theta)^{d-2}~ \mathrm{d}\theta$ with the factor $N=\int_0^{\pi} (\sin \theta)^{d-2}~ \mathrm{d}\theta =B(({1/ 2}), {(d-1)/ 2})$ ensuring proper normalization (B(a,b) is the beta function). Using $\mu=\cos^2 \theta$ as the integration variable, the above governing equation for the scaling function reads

 \begin{equation} -\lambda ~x~\Phi'(x)+\Phi(x) =\int_0^1 {\cal D}\mu~~\Phi(\xi x)~~\Phi(\eta x), \end{equation} (8)

where $\xi\equiv \xi(\epsilon,\mu)=1-(1-\epsilon^2)\mu$ and $\eta\equiv\eta(\epsilon,\mu)=(1-\epsilon)^2\mu$. Additionally, the integration measure was re-written using ${\cal D}\mu$, defined via $B({1/ 2},{(d-1)/ 2}){\cal D}\mu= \mu^{-{(1/ 2)}}~(1-\mu)^{(d-3)/ 2}~\mathrm{d}\mu$ (it remains normalized, $\int_0^1{\cal D}\mu=1$). In the elastic case, the velocity distribution is Maxwellian. Indeed, $\xi+\eta=1$ and $\lambda=0$ when $\epsilon=0$, and thence $\Phi(x)=\mathrm{e}^{-x/2}$.

Our primary goal is to determine statistics of extremely fast particles, namely the tail of the velocity distribution. This can be accomplished by noting that the large-v behaviour of the velocity distribution is reflected by the small-k behaviour of its Fourier transform. For example, the small-x expansion of the one-dimensional solution (6) contains both regular and singular terms: $\Phi(x)=1-({1/ 2})~x+({1/ 3})~x^{3/2}+\cdots$, and the dominant singular x3/2 term reflects the w-4 tail of ${\cal P}(w)$. In general, an algebraic tail of the velocity distribution (3),

 \begin{equation} {\cal P}(w)\sim w^{-\sigma} \tqs \mbox{when } w\to\infty, \end{equation} (9)

indicates the existence of a singular component in the Fourier transform,

 \begin{equation} \Phi_{\mathrm{sing}}(x)\sim x^{(\sigma-d)/2} \tqs\mbox{when } x\to 0, \end{equation} (10)

and vice versa. This is seen by recasting the Fourier transform $\Phi(x)\propto \int_0^\infty \mathrm{d} w~ w^{d-1}{\cal P}(w)\mathrm{e}^{\mathrm{i} w\sqrt{x}}$, into the Laplace transform $I(s)\propto \int_0^\infty \mathrm{d} w~w^{d-1}{\cal P}(w)\mathrm{e}^{-ws}$ using x = -s2. The small-s expansion of the integral I(s) contains regular and singular components. For example, when $\sigma<d$, the integral I(s) diverges as $s\to 0$ and integration over large w yields the dominant contribution $I_{\mathrm{sing}}(s)\sim s^{\sigma-d}$. When $d<\sigma<d+1$, I(0) is finite, but the next term is the above singular term, so $I(s)=I(0)+I_{\mathrm{sing}}(s)+\cdots$. In general, the singular contribution is $I_{\mathrm{sing}}(s)\sim s^{\sigma-d}$, thereby leading to the singular term of equation (10).

The exponent $\sigma$ can now be obtained by inserting $\Phi(x)=\Phi_{\mathrm{reg}}(x)+\Phi_{\mathrm{sing}}(x)$ into equation (8) and equating the dominant singular terms. Combining $\Phi_{\mathrm{reg}}(0)=1$ with the anticipated leading singular term of equation (10), we find that the exponent $\sigma$ is a root of the following integral equation:

 \begin{equation} 1-\lambda~{\sigma-d\over 2} =\int_0^1 {\cal D}\mu~[\xi^{(\sigma-d)/2}+\eta^{(\sigma-d)/2}]. \end{equation} (11)

This equation has a trivial solution $\sigma=d~+~2$, following from the identity $\int{\cal D}\mu~(\xi+\eta)=1-\lambda$, where the singular and the regular components simply coincide, $x^{(\sigma-d)/2}=x$. Since we seek the leading singular term, the solution of equation (11) must therefore satisfy $\sigma\gt d+2$.

The integral equation (11) can be written explicitly in terms of special functions

 \begin{equation} \eqalign{ \fl 1-\epsilon(1-\epsilon)~{\sigma-d\over d}= {}_2F_1... ...)/2})\Gamma({d/ 2})\over \Gamma({\sigma/ 2})~\Gamma({1/ 2})},\cr } \end{equation} (12)

with 2F1(a,b;c;z) the hypergeometric function [21]. Interestingly, the exponent $\sigma\equiv\sigma(d,\epsilon)$ is a root of the transcendental equation (12) and thence it depends in a nontrivial fashion on spatial dimension d and the dissipation parameter $\epsilon$.

We first consider the dependence on the dissipation parameter by considering the quasi-elastic limit $\epsilon\to 0$. As discussed above, in the elastic case the velocity distribution is Maxwellian and the Fourier transform is simplyNote1 $\Phi(x)=\mathrm{e}^{-x/2}$. This implies a diverging exponent $\sigma\to\infty$ as $\epsilon\to 0$. Therefore, the right-hand side of equation (12) vanishes, and the leading behaviour is

 \begin{equation} \sigma\simeq {d\over \epsilon}\tqs\mbox{as } \epsilon\to 0. \end{equation} (13)

One can further expand $\sigma(d,\epsilon)$ in the $\epsilon\to 0$ limit to find $\sigma(d,\epsilon)=d~ \epsilon^{-1}+a_1(d)~\epsilon^{-1/2}+a_2(d)+\cdots$. We merely quote the leading correction in the physically relevant spatial dimensions $a_1(2)=-2(e^{-2}+1)/\sqrt{\pi}$ and $a_1(3)=-\sqrt{3\pi/2}$. Clearly, the quasi-elastic limit is singular. Dissipation, even if minute, seriously changes the nature of the system [6,22].

Next, we discuss the dependence on the dimension. First, one can verify that $\sigma=4$ when d = 1 using the identity 2F1(a,b;b;z) = (1-z)- a. The case d = 1 is unique in that the entire scaling function and in particular the exponent are not dependent on $\epsilon$. The case of $d\to\infty$ is similar to the $\epsilon=0$ case in that the inelastic nature of the collisions becomes irrelevant, the velocity distribution is Maxwellian and the exponent diverges, $\sigma\to\infty$ as $d\to\infty$. In this limit, the second integral in equation (11) is negligible as it vanishes exponentially with the dimension. The first integral can be evaluated by taking the limits $d\to\infty$ and $\mu\to 0$, with $z=\mu d/2$being finite. Then, the integration measure is transformed ${\cal D}\mu\to (\pi z)^{-1/2}\mathrm{e}^{-z}~\mathrm{d} z$, and equation (11) becomes $1-\epsilon(1-\epsilon)u=\int_0^{\infty} ~\mathrm{d} z~ (\pi z)^{-1/2}\mathrm{e}^{-[1+(1-\epsilon^2)u]z}$ where $u={\sigma/ d}-1$. Performing the integration yields $1-\epsilon(1-\epsilon) u=[1+(1-\epsilon^2)u]^{-1/2}$, from which we find u and

 \begin{equation} \sigma \simeq d~{1+({3/ 2})\epsilon-\epsilon^3-\epsilon^{1/2} (1+({5/ 4})\epsilon)^{1/2}\over \epsilon(1-\epsilon^2)}, \end{equation} (14)

as $d\to\infty$. In general, $\sigma\propto d$, and therefore, the algebraic decay becomes sharper as the dimension increases. The exponent $\sigma(d,\epsilon)$ increases monotonically with increasing dimension, and additionally, it increases monotonically with decreasing $\epsilon$ (see figure 1). Both features are intuitive as they mirror the monotonic dependence of the energy dissipation rate $\lambda=2\epsilon(1-\epsilon)/d$ on d and $\epsilon$. Hence, the completely inelastic case provides a lower bound for the exponent, $\sigma(d,\epsilon)\ge \sigma(d,\epsilon=1/2)$ with $\sigma(d,1/2)=6.287~53$, 8.329 37, for d = 2, 3, respectively. The former value should be compared with $\sigma\approx 5$, obtained from numerical simulations [18]. The algebraic tails are characterized by unusually large exponents which may be difficult to measure accurately in practice (for typical granular particles $\epsilon\approx 0.1$ yielding $\sigma\approx 30$). Figure 1 also shows that the quantity $\sigma/d$ weakly depends upon the dimension, and the large-d limit (14) provides a useful approximation.

\begin{figure} \leavevmode \epsfbox {3304401.eps} \end{figure}

Figure 1. The exact exponent $\sigma$, obtained from equation (11), versus the dissipation parameter $\epsilon$. The exponent was scaled by the dimension d. Also shown is the limiting large-dimension expression (14).

Thus far, we have discussed only freely cooling systems where the energy decreases indefinitely. In typical experimental situations, however, the system is supplied with energy to balance the energy dissipation [10,11]. Theoretically, it is natural to consider white-noise forcing [9], i.e. coupling to a thermal heat bath which leads to a nonequilibrium steady state. Specifically, we assume that in addition to changes due to collisions, velocities may also change due to an external forcing: $(\mathrm{d} v_j/ \mathrm{d} t)\vert _{\mathrm{heat}}=\xi_j$with j = 1,...,d. We use standard uncorrelated white noise $\langle\xi_i(t) \xi_j(t')\rangle =2D\delta_{ij}\delta(t-t')$with a zero average $\langle \xi_j\rangle=0$. The temperature rate equation is modified by the additional source term $({\mathrm{d}/ \mathrm{d} t})T+\lambda T^{3/2}=2D$, and the system approaches a steady state, $T_{\infty}=(2D/\lambda)^{2/3}$. The relaxation toward this state is exponential, $\vert T_{\infty}-T\vert\sim \mathrm{e}^{-t/\tau}$.

Uncorrelated white-noise forcing amounts to diffusion in velocity space, and equation (4) is modified as follows: $({1/ \sqrt{T}})({\partial/\partial t})\to ({1/ \sqrt{T}}) ({\partial/\partial t})+Dk^2$. In the steady state, the Fourier transform, $F(\bi{k},t=\infty)\equiv \Psi(y)$ with y = Dk2, obeys

 \begin{equation} (1+y)~\Psi(y)=\langle \Psi(\xi y)~ \Psi(\eta y)\rangle, \end{equation} (15)

where integration with respect to the measure ${\cal D}\mu$ is denoted by $\langle f\rangle =\int_0^1 {\cal D}\mu f(\mu)$.

Equation (15) is solved recursively by employing the cumulant expansion $\Psi(y)=\exp\big[{\sum}_{n\ge 1} (-y)^nF_n\big]$. Writing $1+y=\exp\big[{\sum}_{n\ge 1}(-y)^n/n\big]$, we recast equation (15) into

 \begin{equation} 1=\biggl\langle\exp\biggl[\sum_{n=1}^{\infty} (-y)^n (n^{-1}-G_n)\biggr]\biggr\rangle, \end{equation} (16)

where $G_n=F_n(1-\xi^n-\eta^n)$. The desired cumulants Fn are obtained by solving for $\langle G_n\rangle$ recursively and then using $F_n=\langle G_n\rangle/\langle 1-\xi^n-\eta^n\rangle$. In one dimension $\langle \mu^n\rangle=1$ and one immediately obtains $\langle G_n\rangle=n^{-1}$ [17]. In higher dimensions, the averages acquire non-trivial dependence on n though the qualitative nature of the solution remains the same since $\langle G_n\rangle$ and hence the cumulants Fn (as well as the moments of the distribution) are positive. This implies the following small-k behaviour of the Fourier transform at the steady state: $\ln F_{\infty}(\bi{k})\sim -(k/k_0)^2$. This behaviour agrees with the prediction of [23] derived via a small-$\epsilon$expansion, but differs from the stretched exponential behaviour exp (-v3/2) found for the driven inelastic hard sphere gas [9].

In summary, we have studied inelastic gases within the framework of the Boltzmann equation with a uniform collision kernel. In the freely evolving case, we have shown analytically that the density of high-energy particles is suppressed algebraically. The algebraic tails are characterized by remarkably large exponents, and may be hard to distinguish from (stretched) exponential tails. Our results, combined with previous kinetic theory studies which find exponential, stretched exponential, and Gaussian tails, indicate that the extremal characteristics can be very sensitive to parameters such as the restitution coefficient, and the dimension [8,9,13]. On the other hand, our results in the forced case support the near-Maxwellian assumptions typically used to obtain macroscopic transport coefficients from kinetic theory [14].

Acknowledgments

We thank A Baldassari and M H Ernst for fruitful correspondence, and G D Doolen and S Redner for useful comments. This research was supported by DOE (W-7405-ENG-36) and NSF(DMR9978902).

References

[1]
Pöschel T and Luding S (ed) 2000 Granular Gases (Berlin: Springer)
[2]
Bernu B and Mazighi R 1990 J. Phys. A: Math. Gen. 23 5745
IOPscience
[3]
McNamara S and Young W R 1992 Phys. Fluids A 4 496
CrossRef
[4]
Goldhirsch I and Zanetti G 1993 Phys. Rev. Lett. 73 1619
CrossRef
[5]
Soto R and Mareschal M 2001 Phys. Rev. E 63 041303
CrossRef
[6]
Ben-Naim E, Chen S Y, Doolen G D and Redner S 1999 Phys. Rev. Lett. 83 4069
CrossRef
[7]
Samadani A, Mahadevan L and Kudrolli A 2001 Preprint cond-mat/0110427
Preprint
[8]
Esipov S E and Pöschel T 1997 J. Stat. Phys. 86 1385
CrossRef
[9]
van Noije T P C and Ernst M H 1998 Granular Matter 1 57
CrossRef
[10]
Olafsen J S and Urbach J S 1998 Phys. Rev. Lett. 81 4369
CrossRef
[11]
Rouyer F and Menon N 2000 Phys. Rev. Lett. 85 3676
CrossRefPubMed
[12]
Barrat A, Biben T, Rácz Z, Trizac E and van Wijland F 2002 J. Phys. A: Math. Gen. 35 463
IOPscience
[13]
Jenkins J T and Richman M W 1985 Phys. Fluids 28 3485
CrossRef
[14]
Sela N and Goldhirsch I 1998 J. Fluid Mech. 361 41
CrossRef
[15]
Goldshtein A and Shapiro M 1995 J. Fluid Mech. 282 75
CrossRef
[16]
A review of elastic Maxwell molecules is given in:

Ernst M H 1981 Phys. Rep. 78 1
CrossRef
[17]
Ben-Naim E and Krapivsky P L 2000 Phys. Rev. E 61 R5
CrossRef
[18]
Baldassarri A, Marini Bettolo Marconi U and Puglisi A 2001 Preprint cond-mat/0111066
Preprint
[19]
Ernst M H and Brito R 2001 Preprint cond-mat/0111093
Preprint
Ernst M H and Brito R 2001 Preprint cond-mat/0112417
Preprint
[20]
Haff P K 1983 J. Fluid Mech. 134 401
CrossRef
[21]
Andrews G E, Askey R and Roy R 1999 Special Functions (New York: Cambridge University Press)
[22]
Puglisi A, Loreto V, Marini Bettolo Marconi U and Vulpiani A 1999 Phys. Rev. E 59 5582
CrossRef
[23]
Carrillo J A, Cercignani C and Gamba I M 2000 Phys. Rev. E 62 7700
CrossRef

Notes

Note1
Every initial velocity distribution evolves towards the Maxwellian distribution if d > 1 and $\epsilon=0$; for d = 1, there is no relaxation in the elastic case.


Please login to access our web services, or create an account if you don't yet have one.

You must have cookies enabled in your web browser to be able to login.

Username
Password

Forgotten your password? Get a new one here.