\documentstyle[aps,multicol,epsf]{revtex}
\begin{document}
\title{Decay Kinetics  in Ballistic Annihilation} 
\author{E.~Ben-Naim, S.~Redner}
\address{Center for Polymer Studies and Department of Physics, 
Boston University, Boston, MA 02215}
\author{F.~Leyvraz}
\address{Instituto de Fisica, Laboratorio der Cuernavaca, UNAM, MEXICO}
\maketitle
\begin{abstract}
We study the kinetics of ballistic
annihilation, $A+A\to 0$, with continuous initial particle velocity
distributions.  The concentration and the rms velocity are found to
decay as $c\sim t^{-\alpha}$ and $v_{\rm rms}\sim t^{-\beta}$ respectively, with the relation $\alpha+\beta=1$
holding in any spatial dimension.  A ``mean-field'' Boltzmann equation
for the evolution of the velocity distribution predicts that $\alpha$ and
$\beta$ depend strongly on the initial condition.  This non-universal
behavior is confirmed numerically in one and two dimensions.  
\end{abstract}
PACS. Numbers: 68.70+w, 03.20.+i, 05.20.Dd, 05.40.+j
\begin{multicols}{2}
For irreversible diffusion-controlled reactions, it is now widely
appreciated that the density decays more slowly than the predictions of
mean-field theory in sufficiently low spatial dimension.  For
two-species annihilation, this behavior is accompanied by the dynamic
formation of large-scale spatial heterogeneities in an initially
homogeneous system [1]. The contrasting situation where the
reactants move ballistically has received much less attention, however,
and relatively little is known.

A number of interesting results have been recently reported for the
kinetics of irreversible aggregation, $A_i+A_j\to A_{i+j}$, with
ballistic trajectories for the aggregates and with momentum conserving
collisions [2,3]. Here the subscript refers to the (conserved) mass
of the aggregates.  This model has been invoked as an idealization of
processes such as the coalescence of fluid vortices [4] and planet
formation by accretion [5]. For the ballistic aggregation model, a
scaling argument suggests that the concentration of aggregates decays as
$t^{-\alpha}$, with $\alpha=2d/(d+2)$, where $d$ is the spatial
dimension [2]. This dimension dependence for all $d$ is atypical of
the behavior pattern exhibited by diffusion-controlled reactions.
Furthermore, microscopic considerations show that the decay of the
density of fixed-mass aggregates disagrees with the predictions of the
scaling argument [3].

Motivated in part by these intriguing features, we investigate the decay
kinetics of the more elementary single-species annihilation process,
$A+A\to 0$, for arbitrary continuous initial velocity distributions.  We
find that the decay of the density depends non-universally on the
initial velocity distribution.  A Boltzmann equation for the evolution
of the velocity distribution accounts for the dependence of the decay
exponent $\alpha$ on the form of the velocity distribution and on the
spatial dimension.  Our predictions are verified in one and two
dimensions by numerical integration of the Boltzmann equation and by
Monte Carlo simulations.  It is worth noting that for one-dimensional
single-species annihilation with a discrete bimodal initial velocity
distribution, $P(v,t=0)\propto p\delta(v-1)+(1-p)\delta(v+1)$, the
density decays as $t^{-1/2}$ for $p=1/2$, while the minority velocity
species decays exponentially for $p\ne 1/2$ [6]. These results can
be inferred by mapping the kinetics onto a first-passage process for a
one-dimensional random walk.  When the velocities are continuously
distributed, this line of reasoning is inadequate to account for the wide
range of possible kinetic behaviors.

At time $t=0$, the system consists of identical particles which are
distributed in space with $P(v,t)i$ the initial concentration of particles
of velocity $v$.  Without loss of generality the average initial
velocity can be chosen to be zero.  The decay kinetics appears to be
independent of the initial spatial distribution of particles and for
simplicity we focus on a random initial distribution.
Particles move according to their initial velocity until a collision
occurs, which results in the removal of both colliding particles.  We
are interested in determining the time dependence of the macroscopic
concentration, $c(t)=\int dvP(v,t)$, and the moments of the velocity
distribution, $\langle v^n\rangle^{1/n}=\big(\int dv\, v^nP(v,t)/c(t)\big)^{1/n}$.

A simple power counting argument relates the density decay exponent
$\alpha$ with the exponent $\beta$ which characterizes the decay of the
typical velocity, $v_{\rm rms}\sim t^{-\beta}$.  Consider a system of
identical particles of fixed radius $r$ at concentration $c$ which move
with a velocity of the order of $v_{\rm rms}$.  From an elementary mean-free
path argument, the time between collisions is $t\sim 1/cv_{\rm rms} r^{d-1}$.
If one assumes the following power law forms for the concentration and
$v_{\rm rms}$, 
\begin{equation}
c\sim t^{-\alpha}\qquad v_{\rm rms}\sim t^{-\beta},
\end{equation} 
then the mean-free path argument indicates that the relation $\alpha+\beta=1$
should hold for all spatial dimension $d$.

Since the lifetime of particles with velocity $v$ is proportional to
$1/v$, faster particles tend to annihilate more quickly, and the typical
velocity should decay in time.  By the relation between $\alpha$ and
$\beta$, a value of $\alpha$ less than unity is therefore implied.  We
further argue that there is a strong dependence of the exponent
$\alpha$ on the form of the initial velocity distribution.  This behavior
differs from the mean field prediction of $c\sim t^{-1}$ which arises
from the naive rate equation $\dot c\propto -kc^2$.

A useful approach for determining the decay kinetics is to write a
Boltzmann equation for the time evolution of the velocity distribution.
For simplicity, consider the case of one dimension; generalization to
higher dimensions follows naturally.  Let $P(x,v,t)$ be the density of
particles with velocity $v$ at position $x$ and at time $t$.  At time
$t+\Delta t$, the velocity distribution changes both because of translation
of particles and because of reactions.  We treat the reaction term in a
mean-field approximation by assuming that a particle at $x'<x$ and
velocity $v'>v$ will necessarily react with the target particle at
$(x,v)$ when $x-x'<(v'-v)\Delta t$.  There is a complementary contribution
due to collisions between the target and a particle located at $x'>x$
with $v'<v$.

\end{multicols}
The sum of these two contributions leads to the Boltzmann equation for
ballistic annihilation
\begin{equation}
P(x+v\Delta t,v,t+\Delta t)-P(x,v,t)=-kP(x,v,t) 
\\\Bigl[\int\limits_v^{\infty}
      \int\limits_{x-(v'-v)\Delta t}^{x} dx'\, P(x',v',t)+
      \int\limits_{-\infty}^v dv'
      \int\limits_{x}^{x+(v-v')\Delta t} dx' \,P(x',v',t)\Bigr],
\end{equation} 
where $k$ is a dimensionless reaction constant.  Since a collision leads
to particle annihilation, there is no collision-induced gain term in the
equation.  This approximate equation overcounts collisions, since the
incident particle at $x'$ may react with a third particle rather than
with the target particle.  We anticipate that such three-body
interactions will have a relatively small effect on the kinetics.

\begin{multicols}{2}
To analyze the Boltzmann equation, we expand Eq.~(1) to first order in $\Delta t$
to arrive at
\begin{eqnarray}
{\partial P(x,v,t)\over \partial t}&=&-v{\partial P(x,v,t)\over \partial x}\\
&&-kP(x,v,t)\int\limits_{-\infty}^{\infty} dv'|v-v'|P(x,v',t)\,. \nonumber
\end{eqnarray} 
Since the initial velocity distribution is spatially homogeneous, and
because of the homogeneity of the reaction process, we assume that the
velocity distribution remains spatially homogeneous.  Thus we write
$P(v,t)$ to signify the time-dependent and spatially homogeneous
concentration of particles with velocity $v$.  This implies that the
convective term $\partial P/\partial x$ in the Boltzmann equation may be set to zero,
leading to
\begin{equation}
{\partial P(v,t)\over\partial t}=
-kP(v,t)\int\limits_{-\infty}^{\infty} dv'\,|v-v'|P(v',t)\,.
\end{equation} 
Thus the $|v-v'|$ dependence of the integral kernel controls the
reaction rate between two particles.  Eq.~(3) is also strongly
reminiscent of the Smoluchowski equation for ballistic
aggregation [3]. Despite the uncontrolled nature of the
approximations underlying Eq.~(3), this formulation gives a useful
quantitative description of the decay kinetics.


The first step in our analysis of the Boltzmann equation is to apply
dimensional analysis, together with the assumed asymptotic behaviors,
$c\sim t^{-\alpha}$ and $v_{\rm rms}\sim t^{-\beta}$, to reduce Eq.~(3) to a 
single variable equation.  From
these considerations, we expect that the velocity distribution will have
the following scaled form 
\begin{mathletters}
\begin{equation}
P(v,t)={c_0\over v_0}\left({t\over t_0}\right)^{\beta-\alpha}f(z),
\end{equation} 
where the scaling function $f(z)$ depends only on the dimensionless
velocity 
\begin{equation} 
z={v\over v_0}\left({t\over t_0}\right)^{\beta}.  
\end{equation} 
Here $t_0=1/(kc_0v_0)$ is the initial characteristic time between reactions,
with $k$ the dimensionless reaction constant, $c_0$ the
initial concentration, and $v_0$ the initial rms velocity.
\end{mathletters}
Upon substituting this scaling form into Eq.~(3), we immediately confirm the
exponent relation $\alpha+\beta=1$.  Additionally, we obtain a rescaled Boltzmann 
 equation for the scaling function 
\begin{equation}
(2\beta-1)f(z)+\beta z f'(z)=-f(z)\int\limits_{-\infty}^{\infty} dz'\,|z-z'|\,f(z').
\end{equation}
Notice that this equation is invariant under the transformation $f(z)\to
a^2f(az)$, so that a unit normalization of $f(z)$ can be achieved by a
scale change in $z$.  To find the large-$z$ tail of the scaled velocity
distribution, we approximate
$|z-z'|\sim|z|$ in the integral and use the fact that $f(z)$ vanishes
as $z\to\infty$.  These steps reduce Eq.~(5) To a simple differential 
equation whose solution is 
\begin{equation}
f(z)\sim|z|^{(1-2\beta)/\beta}\exp(-|z|/\beta)\qquad |z|\gg1. 
\end{equation}

Near the origin, Eq.~(5) admits different solutions for $f(z)$, and
correspondingly, the exponent $\beta$ depends on the form of the initial
velocity distribution.  To treat the small-$z$ limit, we first divide
Eq.~(5) by $f(z)$ to yield
\begin{equation} 
(2\beta-1)+\beta \lim_{z\to 0}z\, (\ln
f(z))'=-2\int_0^\infty dz'\,|z'|\,f(z'). 
\end{equation} 
Consider an initial velocity distribution with a power-law tail in the
small velocity limit, $P(v,t=0)\propto |v|^\mu$.  Our Monte Carlo
simulations (discussed below) reveal that the scaled distribution
retains the same power-law form, $f(z)\sim z^{\mu}$, in the small-$z$
limit.  Adopting this form in Eq.~(5), then the second term is simply
equal to $\beta\mu$.  The resulting equation then predicts that
$\beta$ is a monotonically decreasing function of $\mu$ (for $\mu>-1$)
whose precise form depends on the first moment of $f(z)$.  This
moment, in turn, depends on the full details of the velocity
distribution.  For example, if we take as a trial function
$f(z)\propto z^\mu e^{-z/\beta}$ in Eq.~(6), {\it i.e.}, the product
of the asymptotic behaviors, we obtain $\beta={1\over{3+2\mu}}$. Thus,
by tuning $\beta$, $\alpha=1-\beta$ can be set to any value between 0
and 1. This estimate for $\beta$ manifests the strong dependence of
the decay kinetics on the initial conditions. As the velocity
distribution contains slower particles the concentration decays slower
and conversely, the velocity decays faster.  This estimated form
qualitatively mirrors values for $\beta$ obtained by numerical
integration of Eq.~(3) and on Monte Carlo simulations for various
values of $\mu$.
% THE FOLLOWING TWO SENTENCES CAN BE OMITTED TO MY OPINION.
%More importantly, the relation between $\beta$ and $\mu$ depends only weakly on
%the full details of $f(z)$.  
%Thus any physically reasonable form for
%$f(z)$ provides a functional dependence of $\beta$ on $\mu$ which
%qualitatively agrees with our numerical results.

The generalization of the scaling approach for the
Boltzmann equation to higher spatial dimensions is straightforward.  The
scaled velocity distribution function now takes the form
\begin{equation}
P(v,t)={c_0\over v_0}\left({t\over t_0}\right)^{\beta d-\alpha}f(\vec z\,)\, ,
\quad{\rm with\ \ } \vec z={\vec v\over v_0}\left(t\over t_0\right)^{\beta},
\end{equation}
where the exponent combination $\beta d$ originates from an
integration over $d$-dimensional velocity space. 
The corresponding rescaled equation for $f(\vec z\,)$ becomes
\begin{equation}
\big( (d+1)\beta-1\big) f(\vec z\,)+\beta \vec z\cdot\nabla f(\vec z\,)=-f(\vec z\,)
\int\limits_{-\infty}^{\infty} d \vec z\,'|\vec z-\vec z\,'|f(\vec z\,').
\end{equation} 

To determine the relation between the exponents $\mu$ and $\beta$ we
again assume a small-$z$ power law form $f(\vec z\,)\propto|\vec
z\,|^{\mu}$, with $\mu >-d $, and a relatively sharp cutoff in $f(\vec
z)$ for large $|\vec z\,|$.  When employed in the Boltzmann equation,
these assumptions again lead to a qualitatively correct
$\mu$-dependence of $\beta$.  Namely, $\beta$ monotonically decreases
with $\mu$, and has the limits $\beta\to 1$ for $\mu\to -d$, and
$\beta\to 0$ for $\mu\to\infty$.  As in the case of one dimension, if
we adopt $f(\vec z\,)\propto |\vec z\,|^\mu e^{-|\vec z\,|/\beta}$,
then we find $\beta={1\over{1+2d+2\mu}}$.  Conversely, as the spatial
dimension increases, $\beta$ decreases systematically, and the limit
$\alpha=1$ is approached but never reached.  This behavior corresponds
to the naive rate equation $\dot c = -kc^2$.  Thus only in the
$d=\infty$ limit are particle trajectories sufficiently transparent
that the typical velocity does not decrease.  This is in contrast to
the situation of many diffusion-controlled reactions for which
``transparent'' behavior occurs for only $d>=d_c$ with $d_c$ finite
[1].

To test our approximate analysis, we have performed direct numerical
integration of the Boltzmann equations, Eq.~(3) and Eq.~(5) (Table 1).  In
one dimension, a typical integration was based on dividing the velocity
range $[-1,1]$ into 200 bins with a time step of $\Delta t=0.15$.  A
finer level of resolution gives essentially identical results.  The
integration was performed to 1000 time steps.  To estimate the value of
$\beta$, we computed the ``test'' scaling function $f(z;t)_{\rm
test}\propto t^{1-2\beta_{\rm test}}P(v,t)$ at different times, and
adjusted $\beta_{\rm test}$ to achieve the best data collapse between
different data sets.  We quantitatively made this determination by minimizing
the rms deviation between pairs of data sets for which the time differed
by a factor of 2.  By this method, we obtained estimates for the
exponent $\beta$ with an uncertainty of less than 0.005.

We also performed Monte Carlo simulations for ballistic annihilation in
one and two dimensions which are based on two independent approaches.
One method simulates the ballistic motion as a biased random.  In two
dimensions, an particle at $\vec r_i$ is assigned a velocity $(v_{ix},
v_{ix})$, with $|v_{ix}|, |v_{iy}|<1$, according to an initial velocity
distribution with support in $[-1,1]^2$.  A move attempt consists of
picking each occupied site $\vec r_i$ at random and the particle is
moved by an amount $({\rm sgn}(v_{ix}),0)$ with probability $|v_{ix}|$
and by an amount $(0,{\rm sgn}(v_{iy}))$ with probability $|v_{iy}|$.
If a particle lands on an occupied site, both particles are removed from
the system.  The time is then incremented by the inverse of the current
number of particles in the system.  The attraction of this stochastic
method is that it is easily implemented in any spatial dimension.
However, the stochastic nature of the individual particle moves
introduces diffusion in addition to the primary ballistic motion.  This
may lead to a crossover associated with the transition between diffusion
and drift in a biased random walk.

Our second method is an exact time evolution by synchronous
dynamics, an approach which is restricted to one spatial dimension.
Particle velocities and positions are initialized on a periodic
one-dimensional chain.  The collision time for each nearest-neighbor
pair is computed and the minimum pair collision time $\tau_{\rm min}$ is
retained.  The particles then move ballistically over a time $\tau_{\rm
min}$, so that the particle pair whose collision time equals $\tau_{\rm
min}$ is removed, and the time is incremented by $\tau_{\rm min}$.
The determination of $\tau_{\rm min}$ and subsequent update of particle
positions by this minimum time interval is then iterated.

\vspace{.15in}
\centerline{
\begin{tabular}{|c|c|c|c|c|}\hline
         &                   &    MC      &     MC    &   NI           \\ 
dimension&  $P(v,t=0)/c_0$   &   $\beta$  &   $\alpha$&  $\alpha$      \\ 
\hline
   1     &  uniform          &   0.20     &   0.77    &    0.77         \\
   1     &$|v|^{-1/2}/4 $    &   0.42     &   0.56    &    0.60         \\
   1     &$|v|^{-4/5}/10$    &   0.66     &   0.32    &    0.37         \\
   2     &  uniform          &   0.10     &   0.89    &    0.91         \\
\hline
\end{tabular}}
{\small Table 1 Numerical values for the decay exponents $\alpha$ and $\beta$ 
based on numerical integration of the Boltzmann equation, Eq.~(4), and
on Monte Carlo (MC) simulations.  Results are given for several
representative initial velocity distributions.}
\vspace{.1in}

Our two simulation methods give essentially identical results and we
quote exponent estimates based on the biased random walk algorithm,
since it can be applied in both one and two dimensions (Table 1).  The
exponents are determined by measuring the slopes of successive pairs of
data points when time-dependent quantities are plotted on a double
logarithmic scale.  Typically there is a non-negligible temporal range
for which the value of the slope is most stable, and we adopt the
average value of the slope in this range as the exponent estimate.  The
accuracy of the simulation can be inferred from the deviation of the
numerical estimate for $\alpha +\beta$ from its expected value of unity.  The
basic conclusion from our numerics is that the decay exponents $\alpha$ and
$\beta$ are indeed non-universal and depend on the nature of the initial
velocity distribution.  The numerical integration of the Boltzmann
equation also provides an excellent approximation for the simulation
results.

The non-universality displayed by ballistic single-species annihilation
with continuous velocity distributions suggests several interesting
avenues for further investigation.  One such situation is ballistic
annihilation with a trimodal initial velocity distribution,
$P(v,t=0)\propto p_+\delta(v-1)+p_0\delta(v)+p_-\delta(v+1)$, with
$p_++p_0+p_-=1$.  This system exhibits considerably richer kinetics than
that of ballistic annihilation with a bimodal velocity
distribution [6,7]. For the symmetric situation of $p_+=p_-\equiv
p=(1-p_0)/2$, numerical simulations reveal a decay which depends
non-universally on $p_0$.  For $p_0\to 0$, the density of stationary
particles decays as $c^{(0)}\sim t^{-\alpha_0}$, with $\alpha_0\cong 1$,
while the density of mobile particles decays as $c^{(\pm)}\sim
t^{-\alpha_\pm}$, with $\alpha_\pm\cong 1/2$, as might be expected.
However when the value of $p_0$ is increased, there is a systematic
decrease in $\alpha_0$ and a corresponding increase in $\alpha_\pm$.
When $p_0$ reaches 0.25, we find $\alpha_0\cong\alpha_\pm\cong 2/3$.
For larger values of $p_0$, $c^{(0)}$ saturates at a finite limiting value
while $c^\pm$ decays faster than a power law.

Another interesting variation is ballistic fusion, $A+A\to A$, in which
the collision product takes on the velocity of one of the incident
particles according to a specified rule (such as retaining the smaller,
or the larger of the two incident velocities).  Numerical simulations of
this process give results which are distinct from those of the
stoichiometrically identical momentum conserving aggregation process.

Finally, it may prove interesting to study single species annihilation
for which the diffusion coefficient of each reactant is drawn from a
continuous distribution.  Particles with a larger diffusion coefficient
will explore a larger area and thus may be expected to decay more
rapidly in time.  Hence, it is reasonable to assume that the average
diffusion coefficient of the surviving particles will decay according to
the power law $\langle D\rangle\sim t^{-\beta}$.  When this is used in an estimate of
the mean collision time between particles, one finds the relation
$2\alpha/d+\beta=1$ between the diffusion coefficient decay exponent and the
concentration decay exponent.  Based on the behavior observed in
ballistic reactions, we anticipate that non-universality in $\alpha$ and
$\beta$ may also occur for reactions in which the particles possess
continuously distributed diffusion coefficients.

In summary, the kinetics of ballistic annihilation with general
distributions of particle velocities exhibits a rich variety of decay
kinetics.  Both numerical and analytical approaches indicate that there
is non-universality in the exponents that describe the time dependence
of the concentration and the typical velocity.  An approximate theory,
based on a mean-field Boltzmann equation, successfully accounts for the
dependence of these exponents on the initial velocity distribution.  It
is intriguing that an initial velocity distribution with a large
component of slower particles gives a weak decay of the concentration
and relatively faster decay of the typical velocity.  As the spatial
dimension is increased, the ``transparent'' limit $\alpha=1$ is approached
but apparently never reached.

We gratefully acknowledge the Army Research Office and CONACYT for
partial support of this research.

\begin{thebibliography}{99}


\bibitem{} See {\it e. g.}, 
Ya. B. Zeldovich and A. A. Ovchinnikov, {\sl Chem. Phys.} {\bf 28},
215 (1978);D. Toussaint and F. Wilczek, {J. Chem. Phys.} {\bf 78} 2642
(1983); K. Kang and S. Redner, {\sl Phys. Rev. A} {\bf 32}, 435
(1985); G. Zumofen, A. Blumen, and J. Klafter, {\sl J. Chem. Phys.}
{\bf 82}, 3198 (1985); R. Kopelman, {\sl Science} {\bf 241} 1620
(1988); M. Bramson and J. L. Lebowitz, {J. Stat. Phys.} {\bf 65}, 941
(1991); F.~Leyvraz and S.~Redner, {\sl Phys. Rev. A} {\bf 46}, 3132
(1992).
\bibitem{} G.~F.~Carnevale, Y.~Pomeau and W.~R.~Young, {\sl Phys. Rev. Lett.} {\bf 64}, 2913 (1990).
\bibitem{} Y.~Jiang and F.~Leyvraz, (to be published).
\bibitem{} J. C. McWilliams, {\sl J. Fluid Mech.} {\bf 146}, 21 (1984).
\bibitem{} G. Wetherill, in {\sl The Formation and Evolution of Planetary
Systems} (Cambridge Univ.\ Press, 1988).
\bibitem{} Y.~Elskens and H.~L.~Frisch, {\sl Phys. Rev. A} {\bf 31}, 3812 
(1985).
\bibitem{} However a restricted trimodal aggregation-annihilation process gives
essentially the same results as bimodal annihilation (Ref.~2).  See,
W.~S.~Sheu, C.~Van den Broeck, and K.~Lindenberg, {\sl Phys. Rev. A} 
{\bf 43}, 4401 (1991).

\end{thebibliography} 


\section*{Figure Captions} 

Fig 1 Representative Monte Carlo simulation results for the concentration and
the rms velocity vs./ time.  (a) An average over 5000 realizations for
an initial velocity distribution $P(v,t=0)=c_0|v|^{-1/2}/4$ on a 1000-site
lattice at initial concentration $c_0=0.3$.  (b) An average over 500
realizations  for a uniform distribution of velocities on a
$100\times100$ lattice at initial concentration $c_0=0.3$.


\end{multicols}
\end{document}