text:
# Dynamics of ballistic annihilation

## I Introduction

Ballistically controlled reactions provide simple examples of non-equilibrium systems with complex kinetics. They consist of an assembly of particles with a given velocity distribution, moving freely between collisions in a $`d`$-dimensional space. In the simplest version of these models we consider here, when two particles meet, they instantaneously annihilate each other and disappear from the system. Despite its apparent simplicity, this problem is highly non trivial and has attracted substantial interest during the past years Elskens ; Piasecki ; Droz ; BenNaim1 ; BenNaim2 ; Rey ; Kafri ; Blythe ; Krapivsky ; Trizac . The one-dimensional case where the particles can only have two velocities $`\pm c`$ has been studied in a pioneering work by Elskens and Frisch Elskens . In particular, they proved that for a symmetric initial velocity distribution, the particle density is decreasing, in the long time limit, as $`n(t)t^\xi t^{1/2}`$. The case of general distributions in dimension $`d=1`$ was discussed by Piasecki Piasecki , who reduced exactly the annihilation dynamics to a single closed equation for the two-particle conditional probability. Moreover, it was shown in the particular bimodal situation of a discrete velocity distribution ($`\pm c`$) that in one dimension, important dynamical correlations are developing during the time evolution, invalidating mean-field or Boltzmann-like approaches. This exact approach was applied to the case of a three velocity distribution by Droz et al Droz , with the result that the dynamical exponents were strongly depending upon the details of the initial velocity distribution.

No analytical solutions could be found for continuous velocity distributions. In one dimension, Ben-Naim et al.BenNaim1 ; BenNaim2  have shown that the exponent $`\xi `$ could depend on the behavior near the origin of the initial velocity distribution. This problem has been revisited by Rey et al Rey . Based on the exact theoretical approach Piasecki ; Droz , a dynamical scaling theory was derived and its validity supported by numerical simulations for several velocity distributions. This leads to the conjecture that all continuous velocity distributions $`f(v)`$ that are symmetric, regular, and such that $`f(0)0`$ are attracted, in the long-time regime, towards the same distribution and thus belong to the same universality class. This conjecture was reinforced by numerical simulations in two dimensions Trizac .

The case of a continuous velocity distribution has also been approached recently by Krapivsky and Sire Krapivsky . Starting from a Boltzmann equation, they investigated the decay of the particle density $`n(t)t^\xi `$ and the root mean-square velocity $`\overline{v}t^\gamma `$. They derived upper and lower bounds for the exponents as well as their leading expansion in $`1/d`$, valid in high dimension. The main question with such an approach concerns the validity of a Boltzmann equation. This is not justified in 1D and remains an open problem in higher dimensions.

The purpose of this paper is to give a first principle answer to this type of question. The article is organized as follows. In section II, an original analytical derivation of the equations governing the dynamics of ballistic annihilation is given. The hierarchy equations obeyed by the reduced distributions are obtained. It is shown that in the Grad limit, the hierarchy formally reduces to a Boltzmann-like form for $`d>1`$. If the initial reduced distributions factorize, the whole hierarchy reduces to one non-linear equation. In the second part of section II, a scaling analysis of the exact spatially homogeneous hierarchy is performed. The exponents $`\xi `$ and $`\gamma `$ are showed to depend only on one parameter $`\alpha `$ related to the dissipation of energy. This scaling analysis turns out to be invalid for the case $`d=1`$ with discrete velocity distributions, but correct in the continuous case. Strong arguments are given in favor of the validity of the Boltzmann approach for the case $`d>1`$ in the long time limit. The Boltzmann equation is then solved within a systematic approximation based on an expansion in Sonine polynomials (section III). The first non-Gaussian corrections to the scaled velocity distribution are computed and predictions for the exponents $`\xi `$ and $`\gamma `$ are explicitly worked out as a function of the dimension $`d`$ and the parameter $`\mu `$ characterizing the small velocity behaviour of the initial velocity distribution: \[$`f(𝐯,t=0)|𝐯|^\mu `$ for $`|𝐯|0`$\]. These predictions for $`\xi `$ and $`\gamma `$ are asymptotically exact for large dimensions, and reproduce the $`1/d`$ correction to the mean-field values. In 1D, they are in very good agreement with the exponents reported in the literature Krapivsky  at the Boltzmann level. In 2D, we implement extensive Direct Simulation Monte Carlo methods (DSMC) where the non-linear Boltzmann equation is solved, and Molecular Dynamics (MD) simulations where the exact equations of motion are integrated (section IV). The agreement between the MD and DSMC routes confirms the validity of the Boltzmann approach, and the decay exponents measured are in exceptionally good agreement with the Sonine prediction. Conclusions are drawn in section V. A preliminary account of part of the results presented here has been published elsewhere Trizac .

## II Exact results

### II.1 Derivation of the hierarchy

Let $`\mathrm{\Omega }`$ be a region of finite measure in $`R^{2d}`$. We denote by

$$\mu _k^\mathrm{\Omega }(𝐫_1,𝐯_1,\mathrm{},𝐫_k,𝐯_k;t)$$
(1)
the probability density for finding at time $`t`$ exactly $`k`$ particles within $`\mathrm{\Omega }`$ in the states $`(𝐫_j,𝐯_j)\mathrm{\Omega },j=1,2,\mathrm{},k`$, where $`𝐫_j`$ and $`𝐯_j`$ are the position and the velocity vectors, respectively. The knowledge of the densities $`\mu _k^\mathrm{\Omega }`$ for all $`\mathrm{\Omega }R^{2d}`$ and $`k=0,1,2,\mathrm{}`$, defines entirely the state of the system. For a given region $`\mathrm{\Omega }`$ the normalization condition reads

$$\mu _0^\mathrm{\Omega }(t)+\underset{k=1}{\overset{\mathrm{}}{}}_\mathrm{\Omega }𝑑𝐫_1𝑑𝐯_1\mathrm{}_\mathrm{\Omega }𝑑𝐫_k𝑑𝐯_k\mu _k^\mathrm{\Omega }(𝐫_1,𝐯_1,\mathrm{},𝐫_k,𝐯_k;t)=1$$
(2)
where $`\mu _0^\mathrm{\Omega }(t)`$ is the probability of finding the region $`\mathrm{\Omega }`$ void of particles at time $`t`$.

A necessary condition for the occurrence of a pair of particles at the phase space points $`(𝐫_j,𝐯_j),(𝐫_i,𝐯_i)`$ at time $`t>0`$ is that $`𝐫_{ij}=𝐫_i𝐫_j,𝐯_{ij}=𝐯_i𝐯_j`$ belong to the region of the phase space with the characteristic function

$$\chi (𝐫_{ij},𝐯_{ij};t)=\theta (|𝐫_{ij}|\sigma )\{1\theta (𝐫_{ij}𝐯_{ij})\theta (\sigma \sqrt{|𝐫_{ij}|^2(𝐫_{ij}\widehat{𝐯}_{ij})^2})$$
(3)
$$\times \theta (|𝐯_{ij}|t𝐫_{ij}\widehat{𝐯}_{ij}+\sqrt{\sigma ^2|𝐫_{ij}|^2+(𝐫_{ij}\widehat{𝐯}_{ij})^2})\}$$

where $`\sigma `$ denotes the particle diameter, $`\widehat{𝐯}_{ij}`$ is a unit vector in the direction of the relative velocity, and $`\theta `$ denotes the Heaviside distribution. Indeed, moving backward in time the particles collide during the time interval $`(0,t)`$ if and only if the following three conditions are simultaneously satisfied

$$\begin{array}{cccc}\hfill \text{(i) }& 𝐫_{ij}𝐯_{ij}>0\hfill & \text{ particles approach each other}\hfill & \\ \hfill \text{(ii) }& \sigma >\sqrt{|𝐫_{ij}|^2(𝐫_{ij}\widehat{𝐯}_{ij})^2}\hfill & \text{ the impact parameter is smaller than }\sigma \hfill & \\ \hfill \text{(iii) }& |𝐯_{ij}|t>𝐫_{ij}\widehat{𝐯}_{ij}\sqrt{\sigma ^2|𝐫_{ij}|^2+(𝐫_{ij}\widehat{𝐯}_{ij})^2}\hfill & \text{ the time }t\text{ is long enough}\hfill & \\ & & \text{ for the overlapping configuration to occur}\hfill & \end{array}$$

Hence, $`\chi (𝐫_{ij},𝐯_{ij};t)=1`$, if and only if no overlapping takes place during the time interval $`(0,t)`$.

At time $`t`$ particles $`1,2,\mathrm{}k`$ occupy in $`\mathrm{\Omega }`$ the one-particles states

$$(𝐫_1,𝐯_1),(𝐫_2,𝐯_2)\mathrm{}(𝐫_k,𝐯_k)$$
(4)
with probability density (1). Using the characteristic function (3) we can construct the probability density for finding the same particles in the phase space configuration

$$(𝐫_1+𝐯_1dt,𝐯_1),(𝐫_2+𝐯_2dt,𝐯_2)\mathrm{}(𝐫_k+𝐯_kdt,𝐯_k)$$
(5)
at time $`(t+dt),dt>0`$. It reads

$$\left[\underset{i<j}{\overset{k}{}}\chi (𝐫_{ij}+𝐯_{ij}dt,𝐯_{ij};t+dt)\right]\mu _k^\mathrm{\Omega }(𝐫_1,𝐯_1,\mathrm{},𝐫_k,𝐯_k;t)$$
(6)
In the limit $`dt0^+`$ the above expression takes the asymptotic form

$$\mu _k^\mathrm{\Omega }(𝐫_1,𝐯_1,\mathrm{},𝐫_k,𝐯_k;t)\left[1+\underset{i<}{\overset{k}{}}\underset{j}{\overset{k}{}}\left(\frac{}{t}+𝐯_{ij}\frac{}{𝐫_{ij}}\right)\chi (𝐫_{ij},𝐯_{ij};t)dt\right]$$
(7)
Using the definition (3) we find

$$\left(\frac{}{t}+𝐯_{ij}\frac{}{𝐫_{ij}}\right)\chi (𝐫_{ij},𝐯_{ij};t)=(\widehat{𝐫}_{ij}𝐯_{ij})\delta (|𝐫_{ij}|\sigma )[1\theta (𝐫_{ij}𝐯_{ij})]$$
(8)
where $`\widehat{𝐫}_{ij}=𝐫_{ij}/|𝐫_{ij}|`$. We denote by $`T^v(i,j)`$ the right hand side of (8) and rewrite it in the form

$$T^v(i,j)=\sigma ^{d1}𝑑\widehat{𝝈}(\widehat{𝝈}𝐯_{ij})\theta (\widehat{𝝈}𝐯_{ij})\delta (𝐫_{ij}\sigma \widehat{\sigma })$$
(9)
Here $`\widehat{𝝈}`$ is the unit vector along the line passing through the centers of the spheres at contact. The integration with respect to the measure $`d\widehat{𝝈}`$ is thus the angular integration over the solid angle in $`d`$-dimensional space. The $`\theta `$ function in (9) restricts this angular integral to the hemisphere corresponding to pre-collisional configurations.

Our aim is to construct the probability density $`\mu _k^\mathrm{\Omega }`$ at time $`(t+dt)`$ for $`dt0^+`$ :

$$\mu _k^\mathrm{\Omega }(𝐫_1+𝐯_1dt,𝐯_1,\mathrm{},𝐫_k+𝐯_kdt,𝐯_k;t+dt)=\mu _k^\mathrm{\Omega }(𝐫_1,𝐯_1,\mathrm{}𝐫_k,𝐯_k;t)$$
(10)
$$+\left(\frac{}{t}+\underset{j=1}{\overset{k}{}}𝐯_j\frac{}{𝐫_j}\right)\mu _k^\mathrm{\Omega }(𝐫_1,𝐯_1,\mathrm{}𝐫_k,𝐯_k;t)dt$$

To this end we have still to add to the term (7) the probability weights of two events. The first corresponds to the presence at time $`t`$ of $`(k+2)`$ particles within $`\mathrm{\Omega }`$ in the states

$$(𝐫_1,𝐯_1),(𝐫_2,𝐯_2),\mathrm{},(𝐫_k,𝐯_k),(𝐫_{+1},𝐯_{k+1}),(𝐫_{k+2},𝐯_{k+2})$$
(11)
The state (11) is then transformed into (4) at time $`(t+dt)`$ as the result of an annihilating collision between the particles $`(k+1)`$ and $`(k+2)`$, during the time interval $`(t,t+dt)`$. According to equations (7) and (8), the rate of the occurrence of binary collisions between pairs $`(i,j)`$ is obtained by applying the operator $`[T^v(i,j)]`$ defined in (9) to the corresponding distribution. Hence, when $`dt0^+`$, the $`(k+1,k+2)`$ annihilation process contributes to the density (10) the term

$$_\mathrm{\Omega }d(k+1)_\mathrm{\Omega }d(k+2)T^v(k+1,k+2)\mu _{k+2}^\mathrm{\Omega }(1,2,\mathrm{},k,k+1,k+2;t)𝑑t$$
(12)
where the shorthand notation $`djd𝐫_jd𝐯_j`$ for $`j=1,2,\mathrm{}`$ has been used.

Finally, we have to take into account the effects of the free flow of particles across the boundary $`\mathrm{\Omega }`$ of the region $`\mathrm{\Omega }`$. Indeed, the $`k`$-particle state can be created or destroyed by an additional particle $`(k+1)`$ leaving or entering the considered region. Denoting by $`\widehat{𝐧}`$ the unit vector normal to $`\mathrm{\Omega }`$ oriented outwards, we get the term

$$𝑑𝐯_{k+1}_\mathrm{\Omega }𝑑S(\widehat{𝐧}𝐯_{k+1})\mu _{k+1}^\mathrm{\Omega }(1,\mathrm{},k,k+1;t)𝑑t$$
(13)
$$=_\mathrm{\Omega }d(k+1)𝐯_{k+1}\frac{}{𝐫_{k+1}}\mu _{k+1}^\mathrm{\Omega }(1,\mathrm{},k,k+1;t)𝑑t$$

Here $`dS`$ is the measure of the surface area, and the equality (13) follows from Gauss’ theorem.

The enumerated events combine together to create the complete rate of change of the probability density $`\mu _k^\mathrm{\Omega }`$. As equivalent events have the same probability measure we can equate (10) with the sum of contributions (7),(12) and (13) obtaining thus the hierarchy equations $`(k=1,2,\mathrm{})`$

$$\left(\frac{}{t}+\underset{j=1}{\overset{k}{}}𝐯_j\frac{}{𝐫_j}\underset{i<}{\overset{k}{}}\underset{j}{\overset{k}{}}T^v(i,j)\right)\mu _k^\mathrm{\Omega }(1,\mathrm{},k;t)=$$
(14)
$$_\mathrm{\Omega }d(k+1)_\mathrm{\Omega }d(k+2)T^v(k+1,k+2)\mu _{k+2}^\mathrm{\Omega }(1,2,\mathrm{},k+2;t)+_\mathrm{\Omega }d(k+1)𝐯_{k+1}\frac{}{𝐫_{k+1}}\mu _{k+1}^\mathrm{\Omega }(1,\mathrm{},k,k+1;t).$$

Finally, the evolution equation for $`\mu _0^\mathrm{\Omega }`$ follows from the normalization condition (2). This completes the derivation of the infinite hierarchy of equations satisfied by the probability densities $`\mu _k^\mathrm{\Omega }`$.

From (14) one can derive in a straightforward way the hierarchy satisfied by the reduced distributions $`f_k(1,2,\mathrm{},k;t)`$. They are relevant for the evaluation of physical parameters, as $`f_k(1,2,\mathrm{},k;t)d1\mathrm{}dk`$ represents the measure of the number of $`k`$-particle phase space configurations, with $`k`$ particles occupying the one-particle states $`(𝐫_1,𝐯_1),(𝐫_2,𝐯_2)\mathrm{}(𝐫_k,𝐯_k)`$ at time $`t`$. The distributions $`f_k`$ are related to the probability densities $`\mu _k^\mathrm{\Omega }`$ by the equation Ruelle

$$f_k(1,2,\mathrm{},k;t)=\underset{p=0}{\overset{\mathrm{}}{}}\frac{(k+p)!}{p!}_\mathrm{\Omega }d(k+1)\mathrm{}_\mathrm{\Omega }d(k+p)\mu _{k+p}^\mathrm{\Omega }(1,\mathrm{},k,k+1,\mathrm{},k+p;t).$$
(15)
Note that the $`f_k(1,2,\mathrm{},k;t)`$ do not depend on $`\mathrm{\Omega }`$.

In order to derive the evolution equation for $`f_k`$ one has thus to consider the hierarchy equation (14) with $`k`$ replaced by $`(k+p)`$, and use the relation (15). One finds

$$\left(\frac{}{t}+\underset{j=1}{\overset{k}{}}𝐯_j\frac{}{𝐫_j}\underset{i<}{\overset{k}{}}\underset{j}{\overset{k}{}}T^v(i,j)\right)f_k(1,\mathrm{},k;t)=\underset{p=0}{\overset{\mathrm{}}{}}\frac{(k+p)!}{p!}_\mathrm{\Omega }d(k+1)\mathrm{}_\mathrm{\Omega }d(k+p)$$
(16)
$$\{[\underset{j=k+1}{\overset{k+p}{}}𝐯_j\frac{}{𝐫_j}+\underset{j=1}{\overset{k}{}}\underset{i=k+1}{\overset{k+p}{}}T^v(i,j)+\underset{k+1i}{\overset{k+p}{}}\underset{<j}{\overset{k+p}{}}T^v(i,j)]\mu _{k+p}(1,\mathrm{},k+p;t)$$

$$_\mathrm{\Omega }d(k+p+1)_\mathrm{\Omega }d(k+p+2)T^v(k+p+1,k+p+2)\mu _{k+p+2}(1,\mathrm{},k+p+1,k+p+2;t)$$

$$+_\mathrm{\Omega }d(k+p+1)(𝐯_{k+p+1}\frac{}{𝐫_{k+p+1}})\mu _{k+p+1}^\mathrm{\Omega }(1,\mathrm{}k+p+1;t)\}$$

It is then a question of inspection to see that on the right hand side of (16) only the term

$$\underset{p=1}{\overset{\mathrm{}}{}}\frac{(k+p)!}{p!}_\mathrm{\Omega }d(k+1)\mathrm{}_\mathrm{\Omega }d(k+p)\underset{j=1}{\overset{k}{}}\underset{i=k+1}{\overset{k+p}{}}T^v(i,j)\mu _{k+p}(1,\mathrm{},k+p;t)$$
(17)
$$=d(k+1)\underset{j=1}{\overset{k}{}}T^v(j,k+1)f_{k+1}(1,\mathrm{},k,k+1;t)$$

survives. All the remaining terms exactly cancel out.

The hierarchy equations satisfied by the reduced distributions $`f_k`$ describing the annihilation dynamics thus read

$$\left(\frac{}{t}+\underset{j=1}{\overset{k}{}}𝐯_j\frac{}{𝐫_j}\underset{i<}{\overset{k}{}}\underset{j}{\overset{k}{}}T^v(i,j)\right)f_k(1,\mathrm{},k;t)=d(k+1)\underset{j=1}{\overset{k}{}}T^v(j,k+1)f_{k+1}(1,\mathrm{},k,k+1;t).$$
(18)
Consider now equations (18) supposing that the state of the system is spatially homogeneous. In this case the distribution $`f_1`$ does not depend on the particle position. Let us formally take the Grad limit

$$\sigma 0,n(t)\mathrm{},n(t)\sigma ^{d1}=\lambda ^1=const$$
(19)
where

$$n(t)=𝑑𝐯f_1(𝐯;t)$$

The fixed mean free path $`\lambda `$ introduces a relevant length scale, so we pass to dimensionless positions putting

$$𝐫_j=\lambda 𝐱_j,j=1,2,\mathrm{}$$
(20)
With this change of variables the collision operator (9) takes the form

$$T^v(i,j)=[n(t)\sigma ^d]^{d1}\frac{1}{\lambda }𝑑\widehat{𝝈}(\widehat{𝝈}𝐯_{ij})\theta (\widehat{𝝈}𝐯_{ij})\delta [𝐱_{ij}n(t)\sigma ^d\widehat{\sigma }]$$
(21)
We conclude that the term on the left hand side of equation (18) involving the collision operators (21) vanishes in the Grad limit for $`d>1`$, because the dimensionless parameter $`n(t)\sigma ^d`$ tends to zero. Note that this term induces dynamical correlations, hindering the propagation of the molecular chaos factorization. On the other hand, using the definition of $`T^v(j,k+1)`$ we find that the term on the right hand side equals

$$\frac{1}{\lambda }\underset{j=1}{\overset{k}{}}d𝐯_{k+1}d\widehat{𝝈}(\widehat{𝝈}𝐯_{j(k+1)})\theta (\widehat{𝝈}𝐯_{j(k+1)})f_{k+1}(1,\mathrm{},j,\mathrm{},k,𝐫_{k+1}=𝐫_j\sigma \widehat{\sigma },𝐯_{k+1};t)/n(t)$$
(22)
and its prefactor $`1/\lambda `$ remains finite in the same limit. So, formally the hierarchy equations (18) at a given time $`t`$ reduce pointwise in the Grad limit to the Boltzmann-like hierarchy

$$\left(\frac{}{t}+\underset{j=1}{\overset{k}{}}𝐯_j\frac{}{𝐫_j}\right)f_k^B(1,\mathrm{},k;t)=d(k+1)\underset{j=1}{\overset{k}{}}T^v(j,k+1)f_{k+1}^B(1,\mathrm{},k,k+1;t),(k=1,2,\mathrm{})$$
(23)
The hierarchy (23) propagates the factorization of the reduced distributions

$$f_k^B(1,\mathrm{},k;t)=\underset{j=1}{\overset{k}{}}f_k^B(j;t)$$
(24)
Hence, if the initial state is factorized, the whole hierarchy (23) reduces to one non-linear equation

$$\left(\frac{}{t}+𝐯_1\frac{}{𝐫_1}\right)f^B(1;t)=d2T^v(1,2)f^B(1;t)f^B(2;t)$$
(25)
Equation (25) is the Boltzmann kinetic equation corresponding to the annihilation dynamics. In the following section, we shall see that the formal Grad limit taken here, where $`n\mathrm{}`$, is relevant for the description of the annihilation dynamics at late times, even if the density $`n(t)`$ decreases with time.

### II.2 Scaling analysis of the hierarchy

The evolution of the annihilation kinetics shares a common feature with the Grad limit: The ratio of particle diameter to mean-free-path $`\lambda =1/(n\sigma ^{d1})`$, which is related to the packing fraction, vanishes in both cases. To be more specific, we perform a scaling analysis of the exact homogeneous hierarchy equations and look for self-similar reduced distributions where the time dependence has been absorbed into the density $`n(t)`$, with the velocities $`𝐯`$ and positions $`𝐫`$ renormalized by the typical (root mean squared) velocity $`\overline{v}(t)`$ and mean-free-path respectively: We define the reduced variables

$$𝐜=\frac{𝐯}{\overline{v}}\text{and}𝐱=\frac{𝐫}{\lambda },\text{with}(\overline{v})^2=\frac{1}{n(t)}v^2f_1(𝐯;t)𝑑𝐯.$$
(26)
For the one-body distribution, we therefore introduce the reduced function $`\stackrel{~}{f}_1`$ such that

$$f_1(𝐯;t)=\frac{n(t)}{\overline{v}(t)^d}\stackrel{~}{f}_1(𝐜)=\frac{n(t)}{\overline{v}(t)^d}\stackrel{~}{f}_1\left(\frac{𝐯}{\overline{v}(t)}\right).$$
(27)
By definition, the moments of order 0 and 2 of $`\stackrel{~}{f}_1`$ are constrained to unity. Requiring that the $`k`$-body distribution $`f_k`$ factorizes into $`_{i=1}^kf_1(i)`$ in the limit of infinite relative separations between the particles, we consistently obtain the scaling form:

$$f_k(𝐫_1,𝐯_1,\mathrm{},𝐫_k,𝐯_k;t)=\left(\frac{n}{\overline{v}^d}\right)^k\stackrel{~}{f}_k(𝐱_1,𝐜_1,\mathrm{},𝐱_k,𝐜_k).$$
(28)
The evolution equations of $`n(t)`$ and kinetic energy density $`n\overline{v}^2(t)`$ follow from integrating (18) with weights $`d𝐯_1`$ and $`v_1^2d𝐯_1`$ respectively, for $`k=1`$. We obtain

$`{\displaystyle \frac{dn}{dt}}=\omega (t)n`$ (29)
$`{\displaystyle \frac{d(n\overline{v}^2)}{dt}}=\alpha \omega (t)n\overline{v}^2,`$ (30)
where the collision frequency $`\omega `$ and kinetic energy dissipation parameter $`\alpha `$ read

$`\omega (t)=n(t)\overline{v}(t){\displaystyle 𝑑𝐜_1𝑑𝐜_2𝑑\widehat{𝝈}(\widehat{𝝈}𝐜_{12})\theta (\widehat{𝝈}𝐜_{12})\stackrel{~}{f}_2(𝐜_1,𝐜_2,\sigma \widehat{𝝈})}`$ (31)
$`\alpha ={\displaystyle \frac{𝑑𝐜_1𝑑𝐜_2𝑑\widehat{𝝈}(\widehat{𝝈}𝐜_{12})\theta (\widehat{𝝈}𝐜_{12})c_1^2\stackrel{~}{f}_2(𝐜_1,𝐜_2,\sigma \widehat{𝝈})}{\left[c^2\stackrel{~}{f}_1(𝐜)𝑑𝐜\right]\left[𝑑𝐜_1𝑑𝐜_2𝑑\widehat{𝝈}(\widehat{𝝈}𝐜_{12})\theta (\widehat{𝝈}𝐜_{12})\stackrel{~}{f}_2(𝐜_1,𝐜_2,\sigma \widehat{𝝈})\right]}}.`$ (32)
Equation (32) is valid for general velocity rescalings; The definition (26) chosen here implies that the term $`c^2\stackrel{~}{f}_1(𝐜)𝑑𝐜`$ in the denominator equals unity. The coefficient $`\alpha `$ may be seen as the ratio of the kinetic energy dissipated in a typical collision normalized by the average kinetic energy, and is time-independent in the scaling regime. It is convenient to introduce the internal “clock” $`𝒞`$ of the dynamics counting the number of collisions, such that $`d𝒞=\omega dt`$. With this variable, Eqs. (29) and (30) integrate into

$$n(t)=n_0\mathrm{exp}[𝒞(t)]\text{and}\overline{v}^2(t)=\overline{v}_0^2\mathrm{exp}[(\alpha 1)𝒞(t)],$$
(33)
where the time origin with density $`n_0`$ and kinetic energy density $`n_0\overline{v}_0^2`$ has been chosen to coincide with $`𝒞=0`$. Knowledge of the $`𝒞`$-dependence of $`n`$ and $`\overline{v}`$ allows to relate absolute time $`t`$ to the number of accumulated collisions: From Eq. (31), we have

$$\frac{d𝒞}{dt}=\omega _0\frac{n}{n_0}\frac{\overline{v}}{\overline{v}_0}=\omega _0\mathrm{exp}[𝒞(1+\alpha )/2],$$
(34)
where $`\omega _0\omega (t=0)`$, so that

$$𝒞=\frac{2}{1+\alpha }\mathrm{ln}\left(1+\frac{1+\alpha }{2}\omega _0t\right).$$
(35)
The corresponding time evolution is

$`{\displaystyle \frac{n}{n_0}}=\left[1+{\displaystyle \frac{1+\alpha }{2}}\omega _0t\right]^{\frac{2}{1+\alpha }}`$ (36)
$`{\displaystyle \frac{\overline{v}}{\overline{v}_0}}=\left[1+{\displaystyle \frac{1+\alpha }{2}}\omega _0t\right]^{\frac{1\alpha }{1+\alpha }}.`$ (37)
Without knowing the detailed form of the one-particle distribution function $`f_1`$, it is thus possible to conclude about the time decay of $`n`$ and $`\overline{v}`$, which appear to obey algebraic laws in the long time limit \[$`n(t)t^\xi `$ and $`\overline{v}(t)t^\gamma `$, with $`\xi =2/(1+\alpha )`$ and $`\gamma =(\alpha 1)/(\alpha +1)`$\]. The exponents $`\xi `$ and $`\gamma `$ are consequently simply related to the unknown quantity $`\alpha `$, for which a perturbative expansion will be put forward in section III before a numerical investigation in section IV. Moreover, if the initial velocity distribution is of finite support (i.e. vanishes outside a sphere of given velocity $`v^{}`$), $`\overline{v}`$ fulfills the bound $`\overline{v}v^{}`$ so that $`\gamma `$ is necessarily positive (or $`\alpha 1`$). In the framework of Boltzmann’s equation, it will be shown in appendix A that the quantity $`\alpha `$ is necessarily larger than 1. For the specific initial condition where all particles have the same kinetic energy at $`t=0`$, $`\overline{v}^2=v^2`$ is time-independent, and (37) implies that $`\alpha =1`$. From Eq. (36), we therefore obtain the time evolution for this situation

$$\frac{n}{n_0}=\frac{1}{1+\omega _0t},$$
(38)
which is exact within the scaling theory. This relation a priori holds in any dimension, except for $`d=1`$ where the corresponding initial velocity distribution is the symmetric discrete bimodal function $`\pm c`$, for which the scaling ansatz underlying our approach fails (see the discussion at the end of the present section).

Inserting the scaling forms (27) and (28) into the first equation of the hierarchy (18) imposes the following constraint on the decay exponents: $`\xi +\gamma =1`$. This scaling relation may be simply obtained by elementary dimensional analysis BenNaim1 ; BenNaim2 ; Rey ; Krapivsky ; Trizac , and may be considered as the compatibility condition of the hierarchy with the self-similar scaling solutions Trizac\_coal . It is moreover identically fulfilled by the expressions (36) and (37). Under the constraint $`\xi +\gamma =1`$, the remaining equations of the hierarchy ($`k>1`$) turn out to be compatible with (28) with the additional information that the collision term on the lhs of (18) decays like $`t^{\gamma d\xi }`$ whereas the remaining terms are associated with a power $`1/t`$. Given that $`\gamma +d\xi =1+(d1)\xi 1`$, this collision term is asymptotically irrelevant except in one dimension where it remains of the same order as the dominant ones ($`1/t`$). We therefore recover the conclusions obtained by considering the formal Grad limit, with distribution functions expected to obey a Boltzmann-like equation. This analysis points to the relevance of Boltzmann equation for $`d>1`$, a point which is further corroborated by the numerical results given in section IV.

It is interesting to note that both Eq. (38) and the relation $`\xi +\gamma =1`$ do not hold in 1D for discrete initial velocity distributions. In the symmetric situation of a bimodal distribution, the average kinetic energy per particle is conserved (so that $`\gamma =0`$), whereas the density decays as $`1/\sqrt{t}`$ (i.e. $`\xi =1/2`$ Elskens ). The scaling assumption (28) is consequently incorrect in the specific situation of discrete distributions in 1D, but valid for continuous distributions Rey . To be more specific, the scaling form (28) implies that the collision frequency scales with time like $`\omega n\overline{v}`$. On the other hand, from the analytical solution of the bimodal $`\pm c`$ situation Elskens , we obtain $`\omega n^2\overline{v}`$ with $`\overline{v}=c`$. This discrepancy is the signature of dynamical correlations in the latter discrete case. These correlations are responsible for the breakdown of (28), and in addition violate molecular chaos. For continuous velocity distributions, again in 1D, molecular chaos also breaks down while the scaling (28) is correct. As a consequence, the exponents obtained at the Boltzmann level differ from the exact ones (see the discussion in the last paragraph of section III), whereas the relation $`\xi +\gamma =1`$ holds.

## III Boltzmann kinetic equation

This section is devoted to the analysis of the decay dynamics within the molecular chaos framework Resibois  of the homogeneous non-linear Boltzmann equation. No exact solution could be obtained, and our goal is to derive accurate approximate predictions for the scaling exponents $`\xi `$ and $`\gamma `$ of the density and root mean squared velocity.

Before considering the kinetic equation obeyed by the rescaled distribution function, it is instructive to rewrite the original homogeneous Boltzmann equation (25) in the form

$$\frac{f_1(𝐯;t)}{t}=\nu (𝐯;t)f_1(𝐯;t)\text{with}\nu (𝐯_\mathrm{𝟏};t)=\left[\sigma ^{d1}𝑑\widehat{𝝈}(\widehat{𝝈}\widehat{𝐯}_{12})\theta (\widehat{𝝈}\widehat{𝐯}_{12})\right]𝑑𝐯_2|𝐯_{12}|f_1(𝐯_2;t)$$
(39)
where in the last equation, the term in brackets may be computed explicitly as function of dimension (it is understood that the unit vector $`\widehat{𝐯}_{12}`$ denotes an arbitrary direction). For our purpose, it is sufficient to notice that at all times, the collision frequency $`\nu (𝐯;t)`$ of the population having velocity $`𝐯`$ remains finite in the limit $`v0`$, provided the first moment of $`f_1`$ exists. In this situation, Eq. (39) implies that $`f_1(𝐯;t)/f_1(𝐯;0)`$ admits a finite limit for $`v0`$, or equivalently, it may be stated that if the initial velocity distribution behaves like $`v^\mu `$ near the velocity origin, this feature is preserved by the Boltzmann dynamics. Previous works have shown accordingly that the scaling exponents $`\xi `$ and $`\gamma `$ depend on the exponent $`\mu `$ BenNaim1 ; Rey ; Krapivsky ; Trizac .

Making use of relations (29) and (30), insertion of the scaling form (27) into the Boltzmann equation leads to

$$\left[1+\left(\frac{1\alpha }{2}\right)\left(d+c_1\frac{d}{dc_1}\right)\right]\stackrel{~}{f}_1(c_1)=\stackrel{~}{f}_1(c_1)𝑑𝒄_2\frac{c_{12}}{c_{12}}\stackrel{~}{f}_1(c_2),$$
(40)
where we have assumed an isotropic velocity distribution \[$`\stackrel{~}{f}_1(𝐜)=\stackrel{~}{f}_1(c)`$\] and introduced the average $`(\mathrm{})=(\mathrm{})\stackrel{~}{f}_1(c_1)\stackrel{~}{f}_1(c_2)𝑑𝒄_1𝑑𝒄_2`$. $`c_{12}`$ is therefore the rescaled collision frequency. Once $`\mu `$ has been chosen, Eq. (40) admits a solution for a unique value of $`\alpha `$. We show in appendix A that the inequality $`\alpha >1`$ necessarily holds.

Irrespective of $`\alpha `$, the large velocity behaviour of $`\stackrel{~}{f}_1`$ may be obtained following similar lines as in BenNaim1 ; BenNaim2 ; Krapivsky ; vanNoije ; BBRTW : it is possible to integrate formally (40) and cast $`\stackrel{~}{f}_1`$ into

$$\frac{\stackrel{~}{f}_1(c)}{\stackrel{~}{f}_1(c^{})}=\left(\frac{c}{c^{}}\right)^{\frac{2+d(1\alpha )}{1\alpha }}\mathrm{exp}\left[\frac{2}{1\alpha }\frac{1}{c_{12}}_c^{}^c\frac{\stackrel{~}{\nu }(c^{\prime \prime })}{c^{\prime \prime }}𝑑c^{\prime \prime }\right].$$
(41)
In this equation, $`\stackrel{~}{\nu }`$ is itself a functional of $`\stackrel{~}{f}_1`$:

$$\stackrel{~}{\nu }(c_1)=c_{12}\stackrel{~}{f}_1(c_2)𝑑𝐜_2,$$
(42)
such that $`\stackrel{~}{\nu }(c)/c`$ goes to a finite limit for $`c\mathrm{}`$. We therefore obtain the large velocity tail

$$\stackrel{~}{f}_1(c)c^{\frac{2+d(1\alpha )}{1\alpha }}\mathrm{exp}\left(\frac{2}{\alpha 1}\frac{c}{c_{12}}\right)\text{for}c\mathrm{}.$$
(43)
In one dimension, we recover the results of references BenNaim1  and Krapivsky . In BenNaim1 ; BenNaim2 , an approximation was derived for $`\alpha `$ \[or equivalently $`\xi =2/(1+\alpha )`$\] by assuming that the large velocity behaviour of $`\stackrel{~}{f}_1`$ could hold for the whole velocity spectrum. In this picture, the power of $`c`$ appearing on the rhs of Eq. (43) is equated to the exponent $`\mu `$ characteristic of the small velocity behaviour (imposed by the initial distribution chosen, see above), with the result

$$\alpha =\mathrm{\hspace{0.17em}1}+\frac{2}{\mu +d}\text{or}\xi =\frac{2d+2\mu }{2d+2\mu +1}.$$
(44)
This prediction encodes the correct dependence on $`\mu `$ and dimension ($`\xi `$ increases when $`\mu `$ or $`d`$ increase), and turns out to have an accuracy of order 10% when compared to the numerical results Trizac . In the limit of large dimension, we obtain from (44) $`\xi 1(2d)^1`$, whereas Krapivsky and Sire have shown that $`\xi 1d^1(11/\sqrt{2})`$, also in the framework of the Boltzmann equation. The remainder of this section is devoted to the derivation of a more precise value for $`\alpha `$, which furthermore coincides with the exact $`1/d`$ correction for $`d\mathrm{}`$.

Invoking the identity

$$𝑑𝐜c^p\left(d+c\frac{d}{dc}\right)\stackrel{~}{f}_1(𝐜)=pc^p,$$
(45)
the energy dissipation parameter $`\alpha `$ may be given the set of equivalent expressions:

$$\alpha =\mathrm{\hspace{0.17em}1}+\frac{2}{p}\left(\frac{c_{12}c_1^p}{c_{12}c_1^p}1\right),p0.$$
(46)
A particularly useful relation between $`\alpha `$ and moments of $`\stackrel{~}{f}_1`$ follows from considering the limit $`c_10`$ of (40): we get

$$\alpha =1+\frac{2}{\mu +d}\left(1\frac{c_1}{c_{12}}\right).$$
(47)
The (infinite) family of relations (47) and (46) is equivalent to the original integro-differential equation (40), and well suited to a perturbative analysis. To this end, a systematic approximation of the isotropic function $`\stackrel{~}{f}_1`$ can be found by expanding it in a set of Sonine polynomials Landau :

$$\stackrel{~}{f}_1(c)=(c)\left[1+\underset{n=1}{\overset{\mathrm{}}{}}a_nS_n(c^2)\right].$$
(48)
These polynomials are orthogonal with respect to the Gaussian weight

$$(c)=\left(\frac{d}{2\pi }\right)^{d/2}e^{dc^2/2},$$
(49)
and the first few read

$`S_0(x)=\mathrm{\hspace{0.17em}1}`$ (50)
$`S_1(x)={\displaystyle \frac{d}{2}}\left(x+1\right)`$ (51)
$`S_2(x)={\displaystyle \frac{d^2}{8}}x^2{\displaystyle \frac{d(d+2)}{4}}x+{\displaystyle \frac{d(d+2)}{8}}.`$ (52)
The coefficients $`a_n`$ follow from the orthogonality relation $`S_n(c^2)S_m(c^2)(c)𝑑𝐜\delta _{nm}`$: In particular,

$$a_1=\frac{2}{d}S_1(c^2)=\frac{2}{d}\left(1c^2\right)=0$$
(53)
from the definition of rescaled velocities (26). The first non-Gaussian correction is thus embodied in $`a_2`$, that is proportional to the fourth cumulant (kurtosis) of the velocity distribution:

$$a_2=\frac{d^2}{3}\left[c_i^43c_i^2^2\right]=\frac{d}{d+2}c^4\mathrm{\hspace{0.17em}1},$$
(54)
with $`c_i`$ a Cartesian component of $`𝐜`$. Upon truncating Eq. (48) at a finite order $`n`$, we obtain a regular velocity distribution near $`c=0`$. We consequently restrict our analysis to the case $`\mu =0`$ (the dependence on $`\mu `$ has been considered in Trizac ).

It is also noteworthy that any truncation of (48) at arbitrary order $`n`$ leads to a Gaussian high energy behaviour, incompatible with the result (43) corresponding to an overpopulated tail with respect to the Maxwellian. However, it will be shown in section IV that the difference between the truncated expansion (48) and the numerical velocity distribution becomes manifest far in the tail, where the distribution has reached very low probabilities. Consequently, when the moment involved in (46) and (47) are evaluated from the truncation of (48), the accuracy of the result is expected to be better for low orders $`p`$ in (46). Hence the privileged role played by (47), which is of lower order than any of the identities (46). In practice, upon truncating (48) at order $`n`$, the $`n`$ unknowns $`\alpha `$, $`a_2`$, …, $`a_n`$ are computed evaluating the corresponding moments appearing in (47) and in $`n1`$ of the relations (46), among which it is convenient to retain the $`n1`$ even values of $`p`$ ($`p=0`$ excluded). Truncation of (48) at order $`n=2`$ yields precise predictions for $`\alpha `$ and $`\stackrel{~}{f}_1`$, and already at Gaussian order, $`\alpha `$ turns out to be very close to its numerical counterpart: Setting $`n=0`$ (or equivalently $`n=1`$ since $`a_10`$) in (48), we obtain immediately the zeroth order approximation

$$\alpha =\alpha _0=\mathrm{\hspace{0.17em}1}+\frac{2}{d}\left(1\frac{\sqrt{2}}{2}\right)$$
(55)
which corresponds to

$$\xi =\xi _0=\frac{2}{1+\alpha _0}=\frac{2d}{2(d+1)\sqrt{2}}.$$
(56)
For large dimensions, this estimation goes to unity, with the exact $`1/d`$ correction computed in Krapivsky

$$\xi =1\frac{1}{d}\left(1\frac{\sqrt{2}}{2}\right)+𝒪\left(\frac{1}{d^2}\right).$$
(57)
This behaviour turns out to be “universal”, in the sense that the $`\mu `$ dependence does not appear at this order Trizac .

We shall also be interested in the non-Gaussian features of the velocity distribution, that we quantify by the fourth cumulant $`a_2`$. The (cumbersome) calculations at second order in Sonine expansion are detailed Appendix B. We obtain

$`a_2=\mathrm{\hspace{0.17em}8}{\displaystyle \frac{d(2\sqrt{2}3)}{4d^2+d(6\sqrt{2})}}`$ (58)
$`\alpha _2=\alpha _0+{\displaystyle \frac{\sqrt{2}}{16d}}a_2.`$ (59)
The corresponding density exponent follows from $`\xi _2=1/(1+\alpha _2)`$ as before. For $`d\mathrm{}`$, $`a_22(32\sqrt{2})d^1`$, irrespective of $`\mu `$ Trizac  which reinforces the “universal” nature of large $`d`$. The correction to $`\xi `$ carried by $`a_2`$ behaves as $`1/d^2`$ in this limit, and does not affect the $`1/d`$ terms which still coincide with the exact behaviour (57). Both predictions (55) and (59) are such that $`\alpha >1`$, which is required to obtain a normalizable distribution in (41) and (43).

The second order expansion considered here may be improved by consideration of higher order Sonine terms and inclusion of non-linear terms in $`a_2`$. In the related context of inelastic hard spheres, the limitation of working at linear order in $`a_2`$ with neglect of Sonine terms of order $`n3`$ may be found in Santos . Alternatively, keeping non linear terms in $`a_2`$ and neglecting again Sonine terms of order $`n3`$ leads to multiple solutions. A stability analysis is then required to determine which one is stable, as discussed in Brilliantov .

Here, the value obtained for $`a_2`$ is quite small (see below). Our approximate expressions are accurate when compared to the full numerical solution of Boltzmann’s equation, so that we did not calculate any higher order coefficients, nor consider non linear terms in $`a_2`$. The existing body of literature reports, within Boltzmann framework, numerical exponents in 1D that are in excellent agreement with our predictions, already at zeroth order. For the case $`\mu =0`$, expressions (56) and (59) give $`\xi _00.773`$ and $`\xi _20.769`$ whereas the numerical result obtained in Krapivsky  is $`\xi 0.769`$. These exponents are close to their counterparts extracted numerically from the exact dynamics (0.785 $`\pm 0.005`$ in Rey , and more recently 0.804 Com ). The difference between the exact exponents and those obtained assuming molecular chaos is consistent with the conclusion of section II: In 1 dimension, the factorization of the two-body distribution $`\stackrel{~}{f}_2`$ underlying Boltzmann ansatz is not an exact property of the distributions obeying hierarchy (18). On the other hand, for $`d1`$, the molecular chaos exponents are expected to become exact. This property is illustrated in the next section.

## IV Simulation results

The numerical results presented in this section correspond to the situation $`\mu =0`$, unless stated. We refer to Trizac  for the case of diverging ($`\mu <0`$) or depleted ($`\mu >0`$) velocity distributions near $`v=0`$.

### IV.1 The numerical methods

We follow two complementary numerical routes. First, we solve the time-dependent homogeneous Boltzmann equation by means of the Direct Simulation Monte Carlo method (DSMC), originally developed to study ordinary gases Bird . This scheme, where a suitable Markov chain is constructed, has been extended to deal with inelastic collisions Montanero ; Frezzotti  and is easily modified to describe the situation under study here, which does not conserve the total number of particles. Restricting to a spatially homogeneous system, the algorithm is especially easy to implement, and may be summarized as follows: among $`N_0`$ initial particles having a given velocity distribution, a pair $`(i,j)`$ is chosen at random, and removed from the system with a probability proportional to $`|𝐯_{ij}|`$. The (suitably renormalized) time variable is then incremented by the amount $`(N^2|𝐯_{ij}|)^1`$, where $`N`$ is the number of particles remaining in the system, before another pair is drawn etc…This scheme provides the numerical exact solution of Eq. (39), and allows to test the validity of the approximations put forward in section III \[essentially: Truncation at second order in expansion (48) supplemented with calculations performed at linear order in the fourth cumulant $`a_2`$, see Appendix B\]. A precise analysis of the late time dynamics (and especially the computation of velocity distributions) suffers from the concomitant low number of particles left, and the statistical accuracy is improved by averaging over independent realizations.

The second numerical method (Molecular Dynamics Allen ) consists of integrating the exact equations of motion for an assembly of spheres confined in a (hyper)-cubic box with periodic boundary conditions. This route assesses the validity of the approach relying on the homogeneous Boltzmann equation, but does not offer the same accuracy as DSMC, nor the possibility to follow the evolution over comparable times. In particular, once the mean-free-path $`\lambda `$, which grows rapidly as $`t^\xi `$, becomes of the order of the box size $`L`$, the subsequent evolution suffers from finite size effects and should be discarded: When $`\lambda >L`$, the algorithm used is unable to find collision events for those particles which make more than one free flight round on the torus topologically equivalent to the simulation box, which causes a spurious slowing down of the dynamics. It is then tempting to reduce $`\lambda `$ by increasing particles diameter $`\sigma `$, but then finite (and necessarily transient) density effects –incompatible with the scaling assumption (28)– may also arise if the packing fraction $`\varphi `$, proportional to $`n\sigma ^d`$, is not low enough. Simulating explicitly the limit of point particles, the DSMC scheme considered here is free of this defect. The initial number of particles considered in MD needs to be large to allow the system to enter the scaling regime before finite-size effects become dominant: we considered systems with $`N=5.10^4`$ to $`N=5.10^5`$ spheres initially (compared to $`N=10^6`$ to $`10^8`$ in DSMC).

### IV.2 Dynamic scaling behaviour

The results of two-dimensional DSMC and MD simulations are shown in Figure 1, where the quantity on the $`x`$ axis is expected to scale like real time $`t`$ from the scaling relation $`\xi +\gamma =1`$. This log-log plot is a direct probe of the exponent $`\xi `$, from the slope measured. Both MD and DSMC methods give compatible results, with the possibility to follow the dynamics over a longer time interval in DSMC. The departure observed for $`n_0/n200`$ corresponds to the slowing down of MD dynamics resulting from finite-size effects (see Figure 2 below). The theoretical predictions at zeroth and second order are very close ($`\xi _0=0.872`$ and $`\xi _2=0.870`$), and in excellent agreement with the simulation results over several decades. On a similar graph as Fig. 1, the kinetic “temperature” $`\overline{v}^2`$ exhibits a lower law behaviour (not shown) with an exponent $`2\gamma `$ in good agreement with the theoretical prediction ($`\gamma 0.13`$ in 2D). Moreover, the exponents obtained analytically and numerically are compatible to those reported in the literature ($`\xi 0.89`$ in BenNaim1  and $`\xi 0.87(5)`$ using a multi-particle lattice gas method Chopard ).

The time evolution of inverse density and inverse typical velocity square is shown in Figure 2, where considering $`n_0/n1`$ and $`\overline{v}_0^2/\overline{v}^21`$ instead of $`n_0/n`$ and $`\overline{v}_0^2/\overline{v}^2`$ allows to probe the short time behaviour. Unless stated, the initial velocity distribution is an isotropic Gaussian. From Eqs. (29) and (30), $`n`$ and $`\overline{v}^2`$ evolve linearly with $`t`$ for $`\omega _0t1`$ \[see also (36) and (37)\]; the same holds for inverse density and inverse typical velocity squared, which is indeed observed in Fig. 2. MD and DSMC result super-impose, except at late times where MD dynamics suffers from the slowing down discussed previously. For both numerical methods, the scaling relation $`\xi +\gamma =1`$ is well obeyed, in principle at late times only, in the scaling regime. Special combinations of $`n`$ and $`\overline{v}`$ can however be constructed with the requirement to match the short time evolution with the scaling behaviour. One of these quantities is displayed in Fig. 3, with a resulting scaling regime extending over more than 10 decades in time. In Fig. 4, we not only test the validity of the theoretical scaling exponents, but also the full time dependence as predicted by Eqs. (36) and (37). In order to improve the agreement between theory and simulation (which holds over more than 6 orders of magnitude in time), the system has been left time to enter the scaling regime: The time origin $`t=0`$ has been chosen when 80% of the particles originally present have disappeared. The corresponding reference configuration (with subscripts 0) thus differs from the ones considered previously.

### IV.3 Velocity distribution in the scaling regime

In order to understand the reasons for the good agreement between our theoretical predictions and the simulations, we now consider the velocity distribution, restricting to Monte Carlo results (leading to similar conclusions, MD is much more demanding on CPU time and does not allow to investigate detailed features of $`\stackrel{~}{f}_1`$ with the same accuracy). After a transient where the probability distribution function $`\stackrel{~}{f}_1`$ evolves with time, a well defined scaling regime is reached with a time independent $`\stackrel{~}{f}_1(𝐜)`$ shown in Fig. 5 together with the Sonine prediction pushed at second order. The agreement is remarkable, and it is also observed that the Gaussian approximation is already close to the asymptotic rescaled velocity distribution, which is quite surprising in a kinetic process extremely far from equilibrium, with furthermore no conservation laws. Given that our perturbative analytical work relies on the calculation on low order moments of $`\stackrel{~}{f}_1`$, this explains the accuracy of the zeroth order predictions $`\alpha _0`$ and $`\xi _0`$. From Eq. (43), we expect the differences between the Sonine expansion and the numerical distribution to become visible in the high energy tail, which is confirmed by Fig. 6. As predicted, $`\stackrel{~}{f}_1`$ is overpopulated with respect to the Gaussian, and displays a high velocity tail of the form (43) (see the inset).

### IV.4 Evolution toward the asymptotic solutions

Before the scaling regime is attained, $`\stackrel{~}{f}_1`$ is time-dependent, as shown in Fig. 7, where the distributions at different times have been renormalized by $``$ to emphasize the building-up of non-Gaussianities. The evolution towards the scaling solution $`1+a_2S_2`$ can be observed. With respect to the Gaussian, $`\stackrel{~}{f}_1`$ is at all times overpopulated both at large and small velocities (which may be related to the positive sign of $`a_2`$ for the latter case); normalization is ensured by an under-population at intermediate velocities. Figures 5, 6 and 7 show that the indirect measure of $`a_2`$ through the non-Gaussian character of $`\stackrel{~}{f}_1/`$ agrees with the theoretical prediction, but it is also possible to compute directly $`a_2`$ in the simulations through its definition as a fourth cumulant \[Eq. (54)\]. It turns that both methods are numerically fully compatible. Moreover, the Sonine expansion (48) truncated at $`n=2`$ holds at any time, even in the transient regime, with the time dependent fourth cumulant $`a_2`$ measured from (54) (see Fig. 8). This result is not a priori expected and points to the relevance of the expansion (48). We did not try to solve analytically the homogeneous time-dependent Boltzmann equation within the same framework as in the scaling regime, so that we do not have any prediction for the (transient) time-dependence of $`a_2`$ and $`\alpha `$. However, as shown in the inset of Fig. 8, the relation (72) (which reads in 2D $`\alpha =5/4+7a_2/16`$) remarkably holds for all time. Here, the energy dissipation parameter $`\alpha `$ has been computed through the ratio $`c_{12}c_1^2/(c_{12}c_1^2)=c_1^2_{\text{coll}}/c_1^2`$, where $`c_1^2_{\text{coll}}`$, the mean energy dissipated in a collision is computed in the simulations and normalized by the time-dependent mean kinetic energy per particle $`c_1^2`$.

### IV.5 A final remark: Identification of “isobestic” points

For $`\mu =0`$, Figure 7 indicates that during the transient evolution toward scaling, the distributions of reduced velocities have fixed points (for a given initial velocity distribution, the curves corresponding to $`\stackrel{~}{f}_1`$ at different times all pass through common points, that we called isobestic points). This feature has been observed for all initial distributions investigated and appears to be a systematic property of the dynamics, which still holds for non vanishing values of $`\mu `$ (see Fig. 9). We have no analytical explanation for this observation.

## V Conclusion

An analytical derivation of the equations governing the dynamics of an infinite system of spherical particles in a $`d`$-dimensional space, moving freely between collisions and annihilating in pairs when meeting, has been obtained. The hierarchy equations obeyed by the reduced distributions $`f_k(1,2,\mathrm{},k;t)`$ have been derived. In the Grad limit, this hierarchy formally reduces to a Boltzmann-like hierarchy for dimensions $`d>1`$. If these reduced distributions $`f_k^B`$ factorize at the initial time, this factorization remains valid for all times and the whole hierarchy reduces to one non-linear equation.

In the long time limit, the ratio of the particle radius to the mean-free path vanishes. A scaling analysis of the exact homogeneous hierarchy equations has been performed. Self-similar reduced distributions in which the time dependence has been absorbed into the density $`n(t)`$ and the root mean-square velocity $`\overline{v}(t)`$ were introduced. As a result, the exponents $`\xi `$ and $`\gamma `$ describing the decay with time of $`n(t)`$ and $`\overline{v}(t)`$ depend only upon one single parameter $`\alpha `$, related to the dissipation of energy. Moreover, it turns out that in dimension higher than $`1`$, the terms responsible for the violation of molecular chaos are asymptotically irrelevant. Therefore, we recover the conclusions reached in the formal Grad limit and thus, the Boltzmann equation becomes exact in the long time limit in dimensions higher than $`1`$. The above arguments give a first principle justification for the use of the Boltzmann equation approach for $`d>1`$, as well as its limitations, in ballistic annihilation problems, an issue that has been overlooked so far. However, as discussed above, this scaling analysis is incorrect for a $`1`$ dimensional system with discrete initial velocity distribution, a situation for which the scaling assumption underlying our approach fails.

The Boltzmann equation has been solved within an expansion in Sonine polynomials $`S_n`$. Truncation to order $`n=2`$, provides the first non-Gaussian corrections to the scaled velocity distribution, and leads to analytical predictions for the exponents $`\xi `$ and $`\gamma `$ as a function of $`d`$. For large dimension $`d`$, these predictions coincide with the exact $`1/d`$ correction to the naive mean field values ($`\xi =1`$ and $`\gamma =0`$), calculated in Krapivsky . The above analytical predictions are in remarkable agreement with the results of extensive numerical simulations we have performed for two dimensional systems (implementing the complementary Monte Carlo and Molecular Dynamics techniques). In 1D for regular continuous velocity distributions, it is noteworthy that they are in excellent agreement with the numerical solution of the Boltzmann equation, and quite close to the exact values obtained with Molecular Dynamics (4% difference). This last point was unexpected since molecular chaos breaks in 1D. Finally, the time dependence of the reduced velocity distribution function shows an unexplained and remarkable feature, with the existence of fixed (“isobestic”) points, irrespective of initial conditions.

Acknowledgments: We would like to thank F. Coppex for a careful reading of the manuscript, T. van Noije, A. Barrat, F. van Wijland and A. Santos for useful discussions. E.T. acknowledges the hospitality of the Theoretical Physics Department (Genève) where part of this work was performed. J.P. and E.T. benefited from the financial support of the “Programme d’Action Integrée Franco-Polonais Polonium” (contract 03990WF). M.D. acknowledges the support of the Swiss National Science Foundation and of the CNRS through the attribution of a “poste de chercheur associé”.

## Appendix A

Within Boltzmann’s kinetic equation, we show in this appendix that $`\alpha >1`$. To this end, Eq. (40) is rewritten in the form

$$\stackrel{~}{f}_1(c_1)+\left(\frac{1\alpha }{2}\right)\text{div}_{𝒄_1}\left(𝒄_1\stackrel{~}{f}_1(c_1)\right)=\stackrel{~}{f}_1(c_1)𝑑𝒄_2\frac{c_{12}}{c_{12}}\stackrel{~}{f}_1(c_2),$$
(60)
and integrated with weight $`\psi (c_1)d𝐜_1`$, where we choose $`\psi (c)=c^{2/(1\alpha )}`$. Assuming that the moment $`\psi \stackrel{~}{f}_1`$ exists, we obtain after an integration by parts (neglecting surface terms)

$$0=𝑑𝒄_1𝑑𝒄_2\psi (c_1)\frac{c_{12}}{c_{12}}\stackrel{~}{f}_1(c_1)\stackrel{~}{f}_1(c_2).$$
(61)
The right hand side of Eq. (61) is a strictly positive quantity \[except when $`\stackrel{~}{f}_1(c)=\delta (c)`$\], which leads to a contradiction. We therefore conclude that the quantity $`\psi \stackrel{~}{f}_1`$ does not exist. Remembering that $`\psi (c)\stackrel{~}{f}_1(c)c^{\mu +2/(1\alpha )}`$ near the velocity origin, the divergence of $`\psi \stackrel{~}{f}_1`$ implies

$$\frac{2}{1\alpha }+d1+\mu <1.$$
(62)
Supplementing this condition with the normalization constraint $`\mu >d`$, we obtain $`\alpha >1`$. The physical meaning of this condition is that the typical kinetic energy dissipated per particle in a collision is larger than the average kinetic energy of the system at the same time. The temperature $`\overline{v}^2`$ is therefore a decreasing function of time.

## Appendix B Second order truncated Sonine expansion

In this Appendix, we calculate the dominant deviation of $`\stackrel{~}{f}_1`$ from the Gaussian shape, and the associated energy dissipation parameter $`\alpha `$ by setting $`\stackrel{~}{f}_1(c)=(c)\left\{1+a_2S_2(c^2)\right\}`$. The moments appearing in (47) and (46) with $`p=2`$ are then computed as a function of $`a_2`$. Upon writing

$$\alpha =\mathrm{\hspace{0.17em}1}+\frac{2}{\mu +d}\left(1\frac{c_1}{c_{12}}\right)=\frac{c_{12}c_1^2}{c_{12}c_1^2},$$
(63)
the last equality provides an equation for $`a_2`$ which is solved so that $`\alpha `$ is finally explicitly known as a function of input parameters $`\mu `$ and dimension $`d`$. As in the main text, the angular brackets denote averages with weight $`\stackrel{~}{f}_1(c_1)\stackrel{~}{f}_1(c_2)`$

$`\mathrm{}`$ $`=`$ $`{\displaystyle (\mathrm{})\stackrel{~}{f}_1(c_1)\stackrel{~}{f}_1(c_2)𝑑𝐜_1𝑑𝐜_2}`$ (64)
$`=`$ $`{\displaystyle (\mathrm{})(c_1)(c_2)\left\{1+a_2\left[S_2(c_1^2)+S_2(c_2^2)\right]\right\}𝑑𝐜_1𝑑𝐜_2}+𝒪\left(a_2^2\right).`$ (65)
In the following, non-linear terms of order $`a_2^2`$ will be neglected. In the spirit of Pre2 , it is convenient to introduce center-of-mass and relative velocities $`𝐜_1=𝐂+𝐜_{12}/2`$; $`𝐜_2=𝐂𝐜_{12}/2`$ and to define the generic moments

$`M_{np}`$ $`=`$ $`c_{12}^nC^p`$ (66)
$`=`$ $`{\displaystyle 𝑑𝐜_{12}𝑑𝐂c_{12}^nC^p\left(\frac{d}{2\pi }\right)^de^{dc_{12}^2/4dC^2}\left\{1+a_2\left[\frac{d^2}{8}\left(c_1^4+c_2^4\right)\frac{d(d+2)}{4}\left(c_1^2+c_2^2\right)+\frac{d(d+2)}{4}\right]\right\}}`$ (67)
From $`c_1^4+c_2^4=2C^4+2(𝐂𝐜_{12})^2+c_{12}^4/8+C^2c_{12}^2`$, the term $`(𝐂𝐜_{12})^2`$ appearing under the integral sign in (67) becomes $`C^2c_{12}^2/d`$ for symmetry reasons, and making use of $`c_1^2+c_2^2=2C^2+c_{12}^2/2`$, the variables $`𝐜_{12}`$ and $`𝐂`$ decouple in (67). The resulting integrals yield

$$\frac{M_{np}}{M_{np}^0}=\frac{c_{12}^nC^p_{}}{c_{12}^nC^p_0}=\mathrm{\hspace{0.17em}1}+\frac{a_2}{16d}\left\{d(n^2+p^2)2d(n+p)+2np(d+2)\right\}+𝒪\left(a_2^2\right).$$
(68)
In this equation, the subscript $`0`$ refers to averages with Gaussian measure (formally $`a_2=0`$):

$$c_{12}^nC^p_0=(\sqrt{d})^{np}\mathrm{\hspace{0.17em}2}^n\frac{\mathrm{\Gamma }\left(\frac{d+n}{2}\right)\mathrm{\Gamma }\left(\frac{d+p}{2}\right)}{\mathrm{\Gamma }^2(d/2)},$$
(69)
where $`\mathrm{\Gamma }`$ is the Euler function. The moments $`c_{12}`$ and $`c_{12}c_1^2=c_{12}(c_1^2+c_2^2)/2=c_{12}^3/4+c_{12}C^2`$ appearing in (63) are then known

$`{\displaystyle \frac{c_{12}_{}}{c_{12}_0}}=\mathrm{\hspace{0.17em}1}{\displaystyle \frac{1}{16}}a_2+𝒪\left(a_2^2\right)`$ (71)
$`{\displaystyle \frac{c_{12}c_1^2_{}}{c_{12}_0}}=\mathrm{\hspace{0.17em}1}+{\displaystyle \frac{1}{2d}}+{\displaystyle \frac{a_2}{32}}\left(6+{\displaystyle \frac{11}{d}}\right)+𝒪\left(a_2^2\right)`$
$``$ $`{\displaystyle \frac{c_{12}c_1^2}{c_{12}c_1^2}}=\mathrm{\hspace{0.17em}1}+{\displaystyle \frac{1}{2d}}+{\displaystyle \frac{a_2}{8}}\left(2+{\displaystyle \frac{3}{d}}\right)+𝒪\left(a_2^2\right).`$ (72)
For an elastic hard sphere fluid at equilibrium (with thus $`a_2=0`$), this last quantity equals $`1+1/(2d)`$ and represents the ratio of the mean kinetic energy of colliding particles (averaged over successive collision events) to the mean kinetic energy of the population. As expected, this ratio exceeds 1, since typical colliding partners are “hotter” than the mean background. This quantity is easily measured in Molecular Dynamics or Monte Carlo simulations (see e.g. Fig. 8). We also note that the ratio (71) has been computed in vanNoije  at the same level of approximation in the context of rapid granular flows, with the same result. van Noije and Ernst also reported non Gaussian corrections to the cooling rate $`\gamma `$ of an inelastic hard sphere fluid vanNoije

$$\frac{\gamma }{\gamma _0}=\mathrm{\hspace{0.17em}1}+\frac{3}{16}a_2+𝒪\left(a_2^2\right)d,$$
(73)
where $`\gamma _0`$ denotes the cooling rate evaluated assuming Maxwellian velocity distributions. In terms of the moments $`M_{np}`$ introduced in (66), it can be shown that $`\gamma /\gamma _0=M_{30}/M_{30}^0`$ and it is then easily checked that expression (68) reduces to (73) for $`n=3`$ and $`p=0`$, for arbitrary dimensionality.

The remaining unknown quantity is $`c_1`$ which may be calculated following similar lines as above:

$`{\displaystyle \frac{c_1^n}{c_1^n_0}}=\mathrm{\hspace{0.17em}1}+{\displaystyle \frac{a_2}{8}}n(n2)+𝒪\left(a_2^2\right)`$ (74)
$`c_1^n_0=\left({\displaystyle \frac{2}{d}}\right)^{n/2}{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{d+n}{2}\right)}{\mathrm{\Gamma }(\frac{d}{2})}},`$ (75)
from which we extract $`c_1=c`$. As expected, the $`a_2`$ correction in (75) vanishes for $`n=0`$ and $`n=2`$, which follows respectively from the normalization constraint and the definition (26) of $`𝐜`$ implying $`c^2=1`$. Gathering results, we obtain (58) and (59) from Eqs. (63).


format: text

text:
# Cold deconfined matter EOS through an HTL quasi-particle model

## 1 Introduction

The equation of state (EOS) of deconfined matter was determined perturbatively to order $`\alpha _s^{5/2}`$ in QCD; however, due to the large coupling and strong renormalization scheme dependencies, perturbative methods are plagued by bad convergence properties that make them fail in the region of physical interest; clearly, non-perturbative approaches are needed.

For zero quark chemical potential $`\mu =0`$ lattice QCD calculations are seemingly up to the task of determining the EOS for a quark-gluon plasma at finite temperature, and although recently there has been some progress also for non-vanishing $`\mu `$, the determination of the EOS for cold dense matter, which is of importance in astrophysical situations, is currently not possible by means of lattice calculations.

As a remedy for the situation, Peshier *et al.* suggested a method which can be used to map the available lattice data for $`\mu =0`$ to finite $`\mu `$ and small temperatures by describing the interacting plasma as a system of massive quasi-particles (QPs); section 2 gives a short review of this technique. While in their simple QP model the (thermal) masses of the quarks and gluons are approximated by the asymptotic limit of the hard-thermal-loop (HTL) self-energies, we make use of the full, momentum-dependent HTL self-energies in our model (which will be presented in section 3). This amounts to including more of the physically important plasmon-effect than in the simple QP model, so we expect our result to show the general trend of the correction that a full next-to-leading order (NLO) calculation of the self-energies (which is work in progress) would add to the result from Peshier *et al*. In section 4 we compare the results obtained from the two models and present our conclusions as well as an outlook in section 5.

## 2 Mapping lattice data to finite chemical potential using QP models

Assuming a plasma of gluons and $`N_f`$ light quarks in thermodynamic equilibrium can be described as an ideal gas of massive quasi-particles with residual interaction $`B`$, the pressure of the system is given by

$$P(T,\mu )=\underset{i=g,q}{}p_i(T,\mu _i(\mu ),m_i^2)B(m_g,m_q),$$
(1)
where the sum runs over gluons (g), quarks and anti-quarks (q) with respective chemical potential $`0,\pm \mu `$; $`p_i`$ are the model dependent QP pressures and $`m_i`$ are the QP masses, which are functions of the effective strong coupling $`G^2(T,\mu )`$.

Using the stationarity of the thermodynamic potential under variation of the self-energies and Maxwell’s relations, one obtains a partial differential equation for $`G^2`$,

$$a_T\frac{G^2}{T}+a_\mu \frac{G^2}{\mu }=b,$$
(2)
where $`a_T`$, $`a_\mu `$ and $`b`$ are coefficients that are given by integrals depending on $`T,\mu `$ and $`G^2`$. Given a valid boundary condition, a solution for $`G^2(T,\mu )`$ is found by solving the above flow equation by the method of characteristics; once $`G^2`$ is thus known in the $`T,\mu `$ plane, the QP pressure is fixed completely. The residual interaction $`B`$ is then given by the integral

$$B=\underset{i}{}\frac{p_i}{m_i^2}\left(\frac{m_i^2}{\mu }d\mu +\frac{m_i^2}{T}dT\right)+B_0,$$
(3)
where $`B_0`$ is an integration constant that has to be fixed by lattice data (usually by requiring $`P(T_c,\mu =0)=P_{lattice}(T_c)`$). Motivated by the fact that at $`\mu =0`$ and $`TT_c`$ the coupling should behave as predicted by the perturbative QCD beta-function, ref.  used the ansatz

$$G^2(T,0)=\frac{48\pi ^2}{(11N_c2N_f)\mathrm{ln}\frac{T+T_s}{T_c}\lambda }.$$
(4)
The parameters $`\lambda `$ and $`T_s`$ are determined by fitting the entropy of the model (which is independent of $`B`$) to available lattice data at $`\mu =0`$. Using (4) as boundary condition for (2), the above procedure allows one to map the lattice EOS from $`\mu =0`$ to the whole $`T,\mu `$-plane.

## 3 HTL QP model

Doing a perturbative expansion of the QP contribution to the pressure at $`\mu =0`$ for the simple QP model and comparing to the known result one finds that while the Stefan-Boltzmann and leading-order interaction terms are correctly reproduced, only $`1/(4\sqrt{2})`$ of the NLO term (the plasmon effect) is included.

Whereas in a simple QP model $`p_g`$ and $`p_q`$ in (1) are just the free pressure of massive (scalar) bosons and fermions, respectively, the HTL-resummed entropy leads to

$`p_g`$ $`=`$ $`d_g{\displaystyle }{\displaystyle \frac{d^3k}{(2\pi )^3}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{d\omega }{2\pi }}n(\omega )[2\mathrm{I}\mathrm{m}\mathrm{ln}(\omega ^2+\mathrm{k}^2+\widehat{\mathrm{\Pi }}_\mathrm{T})2\mathrm{I}\mathrm{m}\widehat{\mathrm{\Pi }}_\mathrm{T}\mathrm{Re}\widehat{\mathrm{D}}_\mathrm{T}`$
$`+\mathrm{Im}\mathrm{ln}(\mathrm{k}^2+\widehat{\mathrm{\Pi }}_\mathrm{L})+\mathrm{Im}\widehat{\mathrm{\Pi }}_\mathrm{L}\mathrm{Re}\widehat{\mathrm{D}}_\mathrm{L}]`$
$`p_q`$ $`=`$ $`d_q{\displaystyle }{\displaystyle \frac{d^3k}{(2\pi )^3}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{d\omega }{2\pi }}(f_+(\omega )+f_{}(\omega ))[\mathrm{Im}\mathrm{ln}(\mathrm{k}\omega +\widehat{\mathrm{\Sigma }}_+)`$ (5)
$`\mathrm{Im}\widehat{\mathrm{\Sigma }}_+\mathrm{Re}\widehat{\mathrm{\Delta }}_++\mathrm{Im}\mathrm{ln}(\mathrm{k}+\omega +\widehat{\mathrm{\Sigma }}_{})+\mathrm{Im}\widehat{\mathrm{\Sigma }}_{}\mathrm{Re}\widehat{\mathrm{\Delta }}_{}],`$
where $`d_g=2(N_c^21)`$, $`d_q=2N_cN_f`$ for gluons and quarks/anti-quarks, respectively; $`n(\omega )`$ and $`f_\pm (\omega )`$ are the bosonic and fermionic distribution functions and $`\widehat{D}_{T,L}`$, $`\widehat{\mathrm{\Delta }}_\pm `$ are the HTL propagators with $`\widehat{\mathrm{\Pi }}_{T,L}`$ and $`\widehat{\mathrm{\Sigma }}_\pm `$ the corresponding self-energies. This increases the included plasmon effect to $`1/4`$ of the known value, the remainder coming from NLO corrections to $`\widehat{\mathrm{\Pi }}`$ and $`\widehat{\mathrm{\Sigma }}`$.

## 4 Comparisons

To obtain the input parameters $`T_s`$ and $`\lambda `$ we fitted the entropy expressions from both models to lattice data for $`N_f=2`$, normalized and scaled. For the HTL-model one obtains a slightly higher value of $`\lambda ^{HTL}=19.4`$ than in the simple QP model, while $`T_s`$ turns out to be equal. Requiring $`P(T_c)=P_{lattice}(T_c)`$ then fixes $`B_0^{HTL}=0.82`$; comparisons for $`B(T,0)`$ and the effective coupling $`G^2`$ are shown in figures 2 and 2, respectively.

While the numerical values of the coefficients $`a_\mu `$ and $`a_T`$ in (2) are approximately equal in the two models, the coefficient $`b`$ turns out to differ noticeably. This difference is reflected in the shape of the characteristics (figure 4). Interestingly enough, for characteristics of the simple QP model starting very near $`T_c`$ intersections occur, while for the HTL model this is not the case (see figure 4).

Integrating the function $`B`$ along the characteristics (3) and adding the QP contributions one obtains the pressure in the whole $`T,\mu `$-plane. The result for $`T=0`$ is shown in figure 6, where the two curves are simple fits to the numerical results. Calculating also the energy density $``$ for $`T=0`$ one finds that the EOS $`(P)`$ for the HTL-model is nearly the same as for the simple model. Finally, a comparison between the QP model characteristics (which are nearly isobars for small $`\mu `$) starting at $`T=T_cT_c|_{\mu =0}`$ and recent lattice data for the phase transition line is shown in figure 6; the curvature of the slope for the QP model characteristics, $`T_cdT/d\mu ^20.07`$ coincides with the central value obtained from one lattice study.

## 5 Conclusions and Outlook

We have investigated the difference between a simple and an HTL QP model at finite chemical potential: although the latter incorporates only slightly more of the plasmon term, there are considerable differences in the characteristics of the flow equation that eliminate the ambiguities of the solutions for the simple QP model. The differences for the thermodynamical quantities are small (even for large $`\mu `$), which seems to reflect the fact that these are protected by the principle of stationarity. However, they should indicate the trend that an inclusion of the full plasmon term will give. It is encouraging that a first comparison of the QP models and recent lattice results for finite $`\mu T`$ shows some agreement (more extensive studies are in progress).