text:
# Kinetics and scaling in ballistic annihilation

## Abstract

We study the simplest irreversible ballistically-controlled reaction, whereby particles having an initial continuous velocity distribution annihilate upon colliding. In the framework of the Boltzmann equation, expressions for the exponents characterizing the density and typical velocity decay are explicitly worked out in arbitrary dimension. These predictions are in excellent agreement with the complementary results of extensive Monte Carlo and Molecular Dynamics simulations. We finally discuss the definition of universality classes indexed by a continuous parameter for this far from equilibrium dynamics with no conservation laws.

Systems with reacting particles model a rich variety of phenomena and provide prominent situations to develop and test the foundations of non-equilibrium statistical mechanics. In this context, the diffusion controlled first order annihilation process ($`A+A\mathrm{}`$) has been extensively studied and the corresponding decay kinetics is well understood. On the other hand much less is known in the contrasting case where the reactants move ballistically between the collision events, despite the relevance of such motion for growth and coarsening processes Krug ; Gerwinski . A few theoretical results are available in $`d=1`$ dimension for such irreversible kinetic processes with discrete initial velocity distributions. In a pioneering work, Elskens and Frisch show from combinatorial considerations that the particle density $`n(t)`$ decays like $`1/\sqrt{t}`$ for the simplest binary velocity distribution Elskens . Powerful generalizations of this result were obtained still in 1D, either for a larger class of stochastic ballistic annihilation and coalescence models Blythe ; Kafri  or from kinetic theory for discrete multi-velocity distributions Piasecki ; Droz . No exact results could be obtained for the generic case of continuous distributions, where the decay exponents have been computed numerically BenNaim ; Rey ; Krapivsky . Recently however, Krapivsky and Sire considered the latter situation in the framework of the Boltzmann equation (relying on the so-called “molecular chaos” factorization Resibois ) and derived bounds for the exponents as well as their leading large $`d`$ behavior. The existing body of literature has essentially focussed numerically on the one dimensional case, and no accurate predictions seem to be available for the decay exponents.

In this Letter, we obtain predictions for the decay exponents and velocity distribution (assumed initially continuous), revisiting Boltzmann kinetic theory in arbitrary dimension, with the explicit inclusion of non Gaussian corrections to velocity distributions. These predictions are compared both with the existing numerical results in 1D and the expressions derived in BenNaim ; Krapivsky , and further tested against extensive numerical simulations in dimension 2 and 3, following two complementary routes: we first solve the mean-field non-linear Boltzmann equation describing the annihilation process by means of a Monte Carlo scheme, which validates the analytical expressions obtained within the molecular chaos framework; second, we go beyond mean-field and investigate the exact decay kinetics by implementing Molecular Dynamics simulations. The two numerical approaches yield the same exponents, in excellent agreement with the analytical prediction. Finally, we address the question of universality in this process Rey  by partitioning the possible continuous velocity distributions into groups associated with the same asymptotic dynamic scaling behaviour, akin to equilibrium universality classes.

We consider an assembly of identical spherical particles with radius $`\sigma `$ in dimension $`d`$, with initial velocity distribution $`f(𝒗,t=0)`$ and random initial positions. Particles follow free flight motion until a collision occurs which results in the removal of both partners. We are interested in the time evolution of density $`n(t)=f(𝒗,t)𝑑𝒗`$ and typical velocity $`\overline{v}(t)`$, related to the kinetic temperature $`T(t)`$ defined as the variance of the velocity distribution

$$T(t)=\frac{1}{n(t)}v^2f(𝒗,t)𝑑𝒗=(\overline{v})^2.$$
(1)
Insight into the decay kinetics may be gained by writing the rate equations for $`n`$ and $`T`$

$`{\displaystyle \frac{dn}{dt}}=\omega (t)n`$ (2)
$`{\displaystyle \frac{d(nT)}{dt}}=\omega (t)nT_{\text{coll}}=\alpha \omega (t)nT,`$ (3)
where the first line stands for a definition of the instantaneous mean collision frequency $`\omega `$, while $`T_{\text{coll}}`$ is the time dependent total kinetic energy of a colliding pair, which is thus dissipated in a binary encounter, as stated by the rhs equality in Eq. (3). On dimensional grounds, the collision frequency is expected to scale like the inverse time, which together with Eqs. (2) and (3) implies an algebraic time decay for $`n`$ and $`\overline{v}`$, as well as a time-independent energy dissipation parameter $`\alpha `$ \[defined in Eq. (3) as $`\alpha =T_{\text{coll}}/T`$\]. We therefore introduce two exponents $`\xi `$ and $`\gamma `$ such that $`n(t)t^\xi `$ and $`\overline{v}t^\gamma `$ (and $`Tt^{2\gamma }`$). With a ballistic dynamics controlled by the mean-free-path $`\mathrm{}1/(n\sigma ^{d1}`$), the collision frequency may be written as the ratio $`\overline{v}/\mathrm{}`$. From $`\omega 1/t`$, we obtain the scaling relation $`\xi +\gamma =1`$ BenNaim ; Rey ; Krapivsky ; Dufty , which may be combined with the ratio of Eqs. (2) and (3) to give

$$\xi =\frac{2}{1+\alpha }\text{ and }\gamma =\frac{\alpha 1}{\alpha +1}.$$
(4)
Since particles with a higher velocity are likely to disappear with a higher rate than the average particle with temperature $`T`$, we expect $`\alpha =T_{\text{coll}}/T`$ to be larger than 1, so that the typical velocity should decrease with time \[$`\gamma >0`$ from Eq. (4)\]. This moreover explains the failure of the naive mean-field picture where the density decay rate is written $`\dot{n}n^2`$, so that $`n(t)1/t`$. This transparency limit would hold in the absence of collisional correlations ($`\alpha =1`$) which becomes only asymptotically exact in the limit of infinite dimension $`d`$.

We now turn to the computation of $`\alpha `$ within the molecular chaos framework, which is a priori an uncontrolled approximation. It will however be shown to capture the essential collisional correlations missed by the naive mean-field argument, and to provide decay exponents in excellent agreement with their numerical counterparts. The corresponding Boltzmann equation reads

$$\frac{f(𝒗,t)}{t}=f(𝒗,t)𝑑𝒘|𝒗𝒘|f(𝒘,t),$$
(5)
which implies that if the initial distribution behaves like a power law $`|v|^\mu `$ near the velocity origin, this property is preserved at subsequent times by the dynamics, which in turn should affect the exponents $`\xi `$ and $`\gamma `$, expected to depend explicitly on $`\mu `$ (as appears on the analytical predictions of Ben-Naim et al BenNaim  $`\xi =(2d+2\mu )/(2d+2\mu +1)`$, or on the bounds derived by Krapivsky and Sire Krapivsky ). Looking for a scaling solution of the kinetic equation (5), we introduce a rescaled velocity $`𝒄=𝒗/\overline{v}`$ and rescaled single particle distribution function $`\phi `$ through

$$f(𝒗,t)=\frac{n(t)}{\overline{v}^d}\phi (𝒄,t),$$
(6)
so that $`\phi (𝒄,t)`$ is the probability distribution function of the velocity $`𝒄`$ at time $`t`$, satisfying the constraints $`\phi 𝑑𝒄=1`$ and $`c^2\phi 𝑑𝒄=1`$ at any time. If $`f(𝒗,t)`$ evolves into a self-similar decay state, the only relevant time dependence occurs via $`n(t)`$ and $`\overline{v}(t)`$, so that $`\phi (𝒄,t)`$ no longer depends on time and the evolution equation for $`f`$ (assumed isotropic) translates into

$`\left[1+\left({\displaystyle \frac{1\alpha }{2}}\right)\left(d+c_1{\displaystyle \frac{d}{dc_1}}\right)\right]\phi (c_1)=`$
$`\phi (c_1){\displaystyle 𝑑𝒄_2\frac{c_{12}}{c_{12}}\phi (c_2)},`$ (7)
where $`(\mathrm{})=(\mathrm{})\phi (c_1)\phi (c_2)𝑑𝒄_1𝑑𝒄_2`$ so that $`c_{12}|𝒄_1𝒄_2|`$ denotes the rescaled collision frequency.

Equation (7) may be considered as an eigenvalue problem for $`\alpha `$, which has been computed numerically in 1D Krapivsky . However, it is useful to reformulate Eq. (7) into an infinite hierarchy of consistency relations obtained by computing the corresponding moment of order $`p`$:

$$\alpha =\mathrm{\hspace{0.17em}1}+\frac{2}{p}\left(\frac{c_{12}c_1^p}{c_{12}c_1^p}1\right).$$
(8)
Note that the special case $`p=2`$ coincides with the definition of $`\alpha `$ through the kinetic energy dissipation as expressed by Eq. (3): $`\alpha =T_{\text{coll}}/T=c_{12}c_1^2/(c_{12}c_1^2)`$. We look for explicit solutions by expanding $`\phi `$ in a basis of Sonine functions Landau

$$\phi (𝒄)=(𝒄)\left[1+\underset{n=1}{\overset{\mathrm{}}{}}a_nS_n(c^2)\right]$$
(9)
where the polynomials $`S_n`$ are orthogonal with respect to the Gaussian weight $`(𝒄)`$. Computing the averages involved in (8) from the functional expression (9) provides a system of equations for the coefficients $`a_n`$.

In practice, only a few terms are required in the expansion (9) in order to get a precise estimation for $`\alpha `$, provided relations of lowest order $`p`$ as possible are retained among the hierarchy (8). In this respect, taking the limit of vanishing velocity of (7) yields the “optimal” relation involving $`\alpha `$ and moments of $`\phi `$ of order 1:

$$\alpha =1+\frac{2}{\mu +d}\left(1\frac{c_1}{c_{12}}\right),$$
(10)
that we consider as the first equation of (8) corresponding to the limit $`p0^+`$. At Gaussian order for $`\phi `$ \[i.e. truncating (9) at order $`n=0`$\], it is straightforward to get

$$\alpha =\alpha _0=\mathrm{\hspace{0.17em}1}+\frac{2}{d+\mu }\left(1\frac{\sqrt{2}}{2}\right)$$
(11)
which, together with Eq. (4) yields the zeroth order estimation for $`\xi `$:

$$\xi _0=\frac{2d+2\mu }{2(d+\mu +1)\sqrt{2}}.$$
(12)
It is noteworthy that in the limit of large dimension, we obtain $`\xi _01d^1(11/\sqrt{2})+𝒪(1/d^2)`$ irrespective of $`\mu `$, which has been shown to be the exact $`1/d`$ behaviour within Boltzmann molecular chaos framework Krapivsky . The first non Gaussian correction is carried by $`a_2`$ ($`a_1`$ identically vanishes from the definition of temperature vanNoije ) and this coefficient is related to the kurtosis of the velocity distribution: $`a_2`$ is proportional to the fourth cumulant $`c_i^43c_i^2^2`$, where $`c_i`$ is a given Cartesian coordinate of $`𝒄`$. After a lengthy calculation performed at linear order in $`a_2`$, we obtain

$`a_2=\mathrm{\hspace{0.17em}8}{\displaystyle \frac{\mu +d(32\sqrt{2})}{4d^2+6\mu +d(6+4\mu \sqrt{2})}}`$ (13)
$`\alpha _2=\alpha _0+{\displaystyle \frac{\sqrt{2}}{16d}}a_2.`$ (14)
The above predictions rely on a perturbative expansion starting from the Maxwellian $``$ (regular at $`𝒗=\mathrm{𝟎}`$) and are therefore expected to be particularly relevant for $`\mu `$ close to 0. The agreement with the existing numerical data is excellent; an accurate estimation has been reported in 1D within molecular chaos for the much studied $`\mu =0`$ case Krapivsky : $`\xi =0.769(5)`$ whereas we obtain at zeroth order $`\xi _0=0.773`$ from (12) and at second order $`\xi _2=2/(1+\alpha _2)=0.769(3)`$ from Eq. (14). This exponent is compatible with its counterpart extracted from the exact dynamics (0.78 in Rey ). Moreover, we have investigated numerically the annihilation dynamics in higher dimensions by means of a) the Direct Simulation Monte Carlo procedure Bird  (DSMC) solving the non-linear homogeneous Boltzmann equation (5) and b) Molecular Dynamics simulations (MD) implementing the exact dynamics with periodic boundary conditions Allen . The DSMC technique provides precise data for the velocity distributions and decay exponents, and allows to test the validity of the analytical truncated expansion of the scaling form $`\phi `$, leading to (11) or (14). Alternatively, MD results assess the reliability of the molecular chaos ansatz, but are more demanding on computer resources: on the one hand, the system needs to reach very low densities in order to develop the self-similar decay stage where $`f(𝒗,t)`$ takes the scaling form (6), but on the other hand, the mean free path $`\mathrm{}`$ which increases with time like $`t^\xi `$ must remain smaller than the simulation box size $`L`$, which provides a lower bound for $`n(t)`$ or equivalently an upper bound for accessible times before finite size effects hinder the precise determination of $`\xi `$ and $`\gamma `$. In practice, we considered systems with $`N=10^5\mathrm{5\hspace{0.17em}10}^5`$ particles in MD and $`N=10^610^7`$ in DSMC where it is further possible to average over $`10^3`$ to $`10^4`$ replicas to increase the statistics of the velocity distributions, which is crucial for computations at large times with a concomitant low number of particles left.

The results of two dimensional simulations are shown in Fig. 1 where it appears that the MD data are fully compatible with DSMC, although less precise. For $`\omega _0t10^3`$, the MD density and temperature tend to saturate, which corresponds to the upper time limit where $`\mathrm{}L`$, and the subsequent evolution is discarded. The predictions $`\xi _0=0.872`$ and $`\xi _2=0.870`$ for $`\mu =0`$ (indistinguishable in Fig 1) are in good agreement with the simulations, irrespective of the initial $`f(𝒗)`$ chosen (we considered several distributions with the constraint $`\mu =0`$, see below the discussion concerning universality). The above exponent is compatible with that reported in the context of a multi-particle lattice gas method (0.87 Chopard ). Moreover, the initial spatial configuration is irrelevant (the long time dynamics and rescaled velocity distributions are the same starting from a fluid-like structure or from various crystalline arrays), and the scaling relation $`\xi +\gamma =1`$ is seen to be well obeyed in the asymptotic regime (inset of Fig. 1). The same scenario holds in dimensions 3 and 4, where the predictions at zeroth and second order are very close, and indistinguishable from the numerics ($`\xi _00.91`$ in 3D and 0.93 in 4D for $`\mu =0`$). However, the agreement is expected to become worse as $`\mu `$ deviates from 0 (with $`\mu >d`$ to ensure proper normalization). This is confirmed in Tables 1 and 2 which summarize the results obtained for various $`\mu `$, with comparison to the theoretical prediction of Ben-Naim et al. BenNaim  (coinciding with the lower bound for $`\xi `$ obtained in Krapivsky , the upper bound being 1). For $`\mu =0`$, the non Gaussian parameter $`a_2`$ is small \[with an even smaller correction to $`\alpha `$ due to the prefactor $`\sqrt{2}/(16d)`$ in (14)\]. This fourth cumulant however rapidly increases with $`\mu `$, so that inclusion of higher order terms \[$`n=3`$…in (9)\] would be required to obtain the same level of accuracy as for regular distributions near the velocity origin.

In the remainder, we consider the possibility to define universality classes for ballistic annihilation kinetics, in the following sense: does $`\mu `$ completely specify the asymptotic velocity distribution and decay exponents, irrespective of further details concerning the initial conditions Rey ? To answer this question we have run several simulations (MD and Monte Carlo) corresponding to different initial conditions sharing the same $`\mu `$, for several values of this parameter. The corresponding decay exponents $`\xi `$ and $`\gamma `$ are monitored, which provides a first test, however quite insensitive to possible non Gaussianities (see above the numerical proximity between $`\xi _0`$ and the non Gaussian corrected $`\xi _2`$). A more sensitive and severe probe is provided by the kurtosis $`a_2`$, which may be computed in two different ways: first from its definition involving the fourth cumulant $`c_i^43c_i^2^2`$, or alternatively from the direct computation of $`\phi (𝒄)/(𝒄)`$, which may further be compared to the analytical expansion $`1+a_2S_2(c^2)`$ with $`a_2`$ given by Eq. (13) (recall that $`a_10`$). The latter method is illustrated in Fig. 2 where the four initial distributions shown in the inset evolve after a transient towards the same attractor, that is furthermore in quantitative agreement with the Sonine prediction. Moreover the same values of $`\xi `$ and $`\gamma `$ are measured within statistical inaccuracy for the 4 distributions. We have observed the same phenomenology for $`\mu 0`$, which points to the relevance of defining universality classes of initial conditions as distributions having the same regularity exponent $`\mu `$, as conjectured in 1D for $`\mu =0`$ Rey .

In conclusion, we have shown that the non trivial dynamic scaling behaviour of ballistic annihilation may be investigated within Boltzmann kinetic theory, and accurate decay exponents have been explicitly worked out. Their evaluation (12) at zeroth order turns out to be straightforward, but follows from a kinetic equation and is therefore specific to the precise model considered here. A more versatile approach that would apply to any ballistically controlled reaction (including coalescence with arbitrary conservation laws, with or without stochasticity in the reactions) consists in reconsidering the rate equations (2) and (3), and identify the proper energy dissipation parameter $`\alpha `$ before approximating it assuming a Gaussian velocity distribution. This “model-independent” approach gives $`\alpha =1+1/(2d)`$ in the particular case of pure annihilation, which corresponds to $`\xi =4d/(4d+1)`$ (i.e. 0.8, 0.89 and 0.92 in dimensions 1, 2 and 3) in reasonable agreement with the exponents mentioned above (0.77, 0.87 and 0.91 respectively). We conjecture that the exponent $`\xi =4d/(4d+1)`$ becomes exact when the particles annihilate with probability $`p`$ (and collide elastically otherwise), in the limiting case $`p0^+`$ (whereas $`p=1`$ for “pure” annihilation). This hopefully provides an illustration of the central role played by the energy dissipation parameter $`\alpha `$ in ballistically controlled reactions, and calls for further investigations with more involved reactions.