Dynamics of ballistic annihilation - Laboratoire de Physique

14 Pages·2002·0.18 MB·English

PHYSICALREVIEWE66,066111 ~2002!
Dynamics of ballistic annihilation
Jarosław Piasecki,1 Emmanuel Trizac,2 and Michel Droz3
1InstituteofTheoreticalPhysics,UniversityofWarsaw,Hoza69,00-681Warsaw,Poland
2LaboratoiredePhysiqueThe´orique,*Baˆtiment210,Universite´ deParis–Sud,91405Orsay,France
3De´partementdePhysique,Universite´ deGene`ve,CH-1211Gene`ve4,Switzerland
~Received23May2002;published10December2002!
Theproblemofballisticallycontrolledannihilationisrevisitedforgeneralinitialvelocitydistributionsand
anarbitrarydimension.Ananalyticalderivationofthehierarchyequationsobeyedbythereduceddistributions
is given, and a scaling analysis of the corresponding spatially homogeneous system is performed. This ap-
proachpointstotherelevanceofthenonlinearBoltzmannequationfordimensionslargerthan1andprovides
expressionsfortheexponentsdescribingthedecayoftheparticledensityn(t)}t2jandtheroot-mean-square
velocity¯v}t2gintermsofaparameterrelatedtothedissipationofkineticenergy.TheBoltzmannequationis
thensolvedperturbativelywithinasystematicexpansioninSoninepolynomials.Analyticalexpressionsforthe
exponentsjandgareobtainedinarbitrarydimensionasafunctionoftheparametermcharacterizingthesmall
velocity behavior of the initial velocity distribution. Moreover, the leading non-Gaussian corrections to the
scaledvelocitydistributionarecomputed.Theseexpressionsforthescalingexponentsareingoodagreement
withthevaluesreportedintheliteratureforcontinuousvelocitydistributionsind51. Forthetwo-dimensional
case, we implement Monte Carlo and molecular dynamics simulations that turn out to be in excellent agree-
mentwiththeanalyticalpredictions.
DOI:10.1103/PhysRevE.66.066111 PACSnumber~s!: 05.70.Ln,05.20.Dd,73.23.Ad,82.20.Nk
I.INTRODUCTION ior near the origin of the initial velocity distribution. This
problem has been revisited by Rey etal. @6#. Based on the
Ballisticallycontrolledreactionsprovidesimpleexamples exact theoretical approach @2,3#, a dynamical scaling theory
ofnonequilibriumsystemswithcomplexkinetics.Theycon- was derived and its validity supported by numerical simula-
sist of an assembly of particles with a given velocity distri- tionsforseveralvelocitydistributions.Thisleadstothecon-
bution, moving freely between collisions in a d-dimensional jecturethatallcontinuousvelocitydistributions f(v) thatare
space.Inthesimplestversionofthesemodelsconsiderhere, symmetric,regular,andsuchthat f(0)(cid:222)0 areattractedinthe
when two particles meet, they instantaneously annihilate long-timeregimetowardsthesamedistribution,andthusbe-
each other and disappear from the system. Despite its appar- long to the same universality class. This conjecture was re-
ent simplicity, this problem is highly nontrivial and has at- inforced by numerical simulations in two dimensions @10#.
tracted substantial interest during the past years @1–10#.The The case of a continuous velocity distribution has also
one-dimensional case where the particles can only have two beenapproachedrecentlybyKrapivskyandSire@9#.Starting
velocities 6c has been studied in a pioneering work by El- from a Boltzmann equation, they investigated the decay of
skens and Frisch @1#. In particular, they proved that for a the particle density n(t);t2jand the root-mean-square ve-
symmetric initial velocity distribution, the particle density is locity¯v}t2g. They derived upper and lower bounds for the
decreasing, in the long-time limit, as n(t)}t2j}t21/2. The exponents as well as their leading expansion in 1/d, valid in
case of general distributions in dimension d51 was dis- high dimension. The main question with such an approach
cussed by Piasecki @2#, who reduced exactly the annihilation concerns the validity of a Boltzmann equation. This is not
dynamics to a single closed equation for the two-particle justified in the one-dimensional ~1D! case and remains an
conditional probability. Moreover, it was shown in the par- open problem in higher dimensions.
ticular bimodal situation of a discrete velocity distribution The purpose of this paper is to give a first-principles an-
(6c) that in one dimension, important dynamical correla- swer to this type of question. The paper is organized as fol-
tions are developing during the time evolution, invalidating lows.InSec.II,anoriginalanalyticalderivationoftheequa-
mean-field or Boltzmann-like approaches. This exact ap- tions governing the dynamics of ballistic annihilation is
proach was applied to the case of a three-velocity distribu- given.Thehierarchyequationsobeyedbythereduceddistri-
tion by Droz etal. @3#, with the result that the dynamical butions are obtained. It is shown that in the Grad limit, the
exponents were strongly depending upon the details of the hierarchy formally reduces to a Boltzmann-like form for d
initial velocity distribution. .1. If the initial reduced distributions factorize, the whole
Noanalyticalsolutionscouldbefoundforcontinuousve- hierarchy reduces to one nonlinear equation. In Sec. II, a
locity distributions. In one dimension, Ben-Naim etal. @4,5# scaling analysis of the exact spatially homogeneous hierar-
have shown that the exponent jcould depend on the behav- chy is performed. The exponents jand gare shown to de-
pend only on one parameter arelated to the dissipation of
energy. This scaling analysis turns out to be invalid for the
*Unite´ MixtedeRecherche8627duCNRS. case d51 with discrete velocity distributions, but correct in
1063-651X/2002/66~6!/066111~14!/$20.00 66066111-1 ©2002TheAmericanPhysicalSociety
PIASECKI,TRIZAC,ANDDROZ PHYSICALREVIEWE66,066111 ~2002!
the continuous case. Strong arguments are given in favor of where sdenotes the particle diameter, vˆ is a unit vector in
the validity of the Boltzmann approach for the case d.1 in the direction of the relative velocity, iajnd u denotes the
the long-time limit. The Boltzmann equation is then solved Heavisidedistribution.Indeed,movingbackwardintime,the
within a systematic approximation based on an expansion in particles collide during the time interval (0,t) if and only if
Soninepolynomials~Sec.III!.Thefirstnon-Gaussiancorrec- the following three conditions are simultaneously satisfied:
tions to the scaled velocity distribution are computed and ~i! r v .0 ~particles approach each other!, ~ii! s
predictions for the exponents jand gare explicitly worked A ij(cid:149) ij
. ur u22(r vˆ )2 ~the impact parameter is smaller than
out as functions of the dimension d and the parameter m ij ij(cid:149) ij A
s), ~iii!uv ut.r vˆ 2 s22ur u21(r vˆ )2 ~thetimetis
characterizing the small velocity behavior of the initial ve- ij ij(cid:149) ij ij ij(cid:149) ij
locity distribution: @f(v,t50)}uvumfor uvu!0]. These pre- long enough for the overlapping configuration to occur!.
dictions for j and gare asymptotically exact for large di- Hence, x(rij,vij;t)51 if and only if no overlapping takes
place during the time interval (0,t).
mensions,andreproducethe1/d correctiontothemean-field
Attimet,particles1,2,...k occupyinV theone-particle
values. In 1D, they are in very good agreement with the
exponents reported in the literature @9# at the Boltzmann states
level. In 2D, we implement extensive direct simulation
Monte Carlo methods ~DSMC!, where the nonlinear Boltz- ~r1,v1!,~r2,v2!, ...,~rk,vk!, ~4!
mann equation is solved, and molecular dynamics ~MD! withprobabilitydensity~1!.Usingthecharacteristicfunction
simulations, where the exact equations of motion are inte- ~3!, we can construct the probability density for finding the
grated~Sec.IV!.TheagreementbetweentheMDandDSMC
same particles in the phase space configuration
routes confirms the validity of the Boltzmann approach, and
the decay exponents measured are in exceptionally good ~r 1v dt,v !,~r 1v dt,v !, ...,~r 1v dt,v ! ~5!
1 1 1 2 2 2 k k k
agreement with the Sonine prediction. Conclusions are
drawn in Sec.V.Apreliminary account of part of the results at time (t1dt),dt.0. It reads
presented here has been published elsewhere @10#. F G
k
)
x~r 1v dt,v ;t1dt! mV~r ,v , ...,r ,v ;t!.
ij ij ij k 1 1 k k
II.EXACTRESULTS i,j
~6!
A.Derivationofthehierarchy
LetV bearegionoffinitemeasureinR2d. Wedenoteby In the limit dt!01, the above expression takes the
asymptotic form
mV~r ,v , ...,r ,v ;t! ~1! F S D
k 1 1 k k
k ] ]
(
mV~r ,v , ...,r ,v ;t! 11 1v
the probability density for finding ~at time t) exactly k par- k 1 1 k k i,j,i,j ]t ij(cid:149)]rij
ticles within V in the states (r ,v)PV, j51,2,...,k, G
j j
where r and v are the position and the velocity vectors,
respectivjely. Thje knowledge of the densities mV for all V 3x~rij,vij;t!dt . ~7!
k
PR2d andk50,1,2,... definesentirelythestateofthesys-
tem.ForagivenregionV, thenormalizationconditionreads Using definition ~3!, we find
S D
‘ E E ] ]
m0V~t!1k(51 Vdr1dv1(cid:149)(cid:149)(cid:149) Vdrkdvk ]t1vij(cid:149)]rij x~rij,vij;t!
3mV~r ,v , ...,r ,v ;t!51, ~2! 5~rˆij(cid:149)vij!d~uriju2s!@12u~rij(cid:149)vij!#, ~8!
k 1 1 k k
where rˆ 5r /ur u. We denote by Tv(i,j) the right-hand
where m0V(t) is the probability of finding the region-V void side of Eijq. ~8ij! anidj rewrite it in the form
of particles at time t.
E
Anecessary condition for the occurrence of a pair of par-
ticlesatthephasespacepoints(r ,v),(r ,v) attimet.0 is Tv~i,j!5sd21 dsˆ~sˆ v !u~2sˆ v !d~r 2ssˆ!. ~9!
j j i i (cid:149) ij (cid:149) ij ij
thatr 5r2r ,v 5v2v belongtotheregionofthephase
ij i j ij i j
space with the characteristic function
Here sˆ is the unit vector along the line passing through the
centersofthespheresatcontact.Theintegrationwithrespect
x~r ,v ;t!5u~ur u2s!$12u~r v !
ij ij ij ij(cid:149) ij to the measure dsˆ is thus the angular integration over the
3u~s2Aur u22~r vˆ !2! solidangleind-dimensionalspace.TheufunctioninEq.~9!
ij ij(cid:149) ij restricts this angular integral to the hemisphere correspond-
3u~uv ut2r vˆ ing to precollisional configurations.
ij ij(cid:149) ij OuraimistoconstructtheprobabilitydensitymV attime
A k
1 s22ur u21~r vˆ !2!%, ~3! (t1dt) for dt!01:
ij ij(cid:149) ij
066111-2
DYNAMICSOFBALLISTICANNIHILATION PHYSICALREVIEWE66,066111 ~2002!
S D
mkV~r11v1dt,v1, ...,rk1vkdtS,vk;t1dt! D ]1(k v ] 2 (k Tv~i,j! mV~1,...,k;t!
] (k ] ]t j51 j(cid:149)]rj i,j,i,j k
5mV~r ,v , ...,r ,v ;t!1 1 v E E
k 1 1 k k ]t j51 j(cid:149)]rj 52 d~k11! d~k12!Tv~k11,k12!
3mV~r ,v , ...,r ,v ;t!dt. ~10! V V
k 1 1 k k E
3mV ~1,2,...,k12;t!1 d~k11!
k12
Tothisend,wehavestilltoaddtotheterm~7!theprobabil- V
ity weights of two events. The first corresponds to the pres- ]
ence at time t of (k12) particles within V in the states 3vk11(cid:149)]rk11mkV11~1,...,k,k11;t!. ~14!
~r1,v1!,~r2,v2!, ...,~rk,vk!,~r11,vk11!,~rk12,vk12!~.11! Fmianlailzlayt,iothneceovnodliutitoionn~2e!q.uTahtiiosncofomrpmle0Vtefsotlhloewdserfirvoamtiotnheofntohre-
infinite hierarchy of equations satisfied by the probability
densities mV.
The state ~11! is then transformed into Eq. ~4! at time (t From Eqk. ~14!, one can derive in a straightforward way
1dt) as the result of an annihilating collision between the
the hierarchy satisfied by the reduced distributions
particles (k11) and (k12), during the time interval (t,t
f (1,2,...,k;t). They are relevant for the evaluation of
1dt). According to Eqs. ~7! and ~8!, the rate of the occur- k
physical parameters, as f (1,2,...,k;t)d1 dk represents
rence of binary collisions between pairs (i,j) is obtained by k (cid:149)(cid:149)(cid:149)
themeasureofthenumberofk-particlephasespaceconfigu-
applying the operator @2Tv(i,j)# defined in Eq. ~9! to the
rations, with k particles occupying the one-particle states
corresponding distribution. Hence, when dt!01, the (k
(r ,v ),(r ,v ), ...,(r ,v ) at time t. The distributions f
11,k12) annihilation process contributes to the density are1 re1lated2to2the probakbilikty densities mV by the equationk
~10! the term k
~see Ref. @11#!
E E
(‘ ~k1p!!
E E f ~1,2,...,k;t!5 d~k11! d~k1p!
2 d~k11! d~k12!Tv~k11,k12! k p50 p! V (cid:149)(cid:149)(cid:149) V
V V 3mV ~1,...,k,k11,...,k1p;t!.
k1p
3mkV12~1,2,...,k,k11,k12;t!dt, ~12! ~15!
Note that f (1,2,...,k;t) do not depend on V.
k
wheretheshorthandnotationdj[drdv for j51,2,... has In order to derive the evolution equation for f , one has
j j k
been used. thus to consider hierarchy equation ~14!, with k replaced by
Finally, we have to take into account the effects of the (k1p), and use relation ~15!. One finds
free flow of particles across the boundary ]V of the region S D
V. Indeed,thek-particlestatecanbecreatedordestroyedby ] (k ] (k (k
an additional particle (k11) leaving or entering the consid- ]t1j51 vj(cid:149)]rj2i,j i,j Tv~i,j! fk~1,...,k;t!
ered region. Denoting by nˆ the unit vector normal to ]V
E E
oriented outwards, we get the term (‘ ~k1p!!
5 d~k11! d~k1p!
p50 p! V (cid:149)(cid:149)(cid:149) V
HF
E E
k1p ] k k1p
dvk11 ]VdS~nˆ(cid:149)vk11!mkV11~1,...,k,k11;t!dt 3 2j5(k11 vj(cid:149)]rj1j(51 i5(k11 Tv~i,j!
E G
] k1p k1p
5 Vd~k11!vk11(cid:149)]rk11mkV11~1,...,k,k11;t!dt. 1k1(1<i (,j Tv~i,j! mk1p~1,...,k1p;t!
~13! E E
2 d~k1p11! d~k1p12!Tv~k1p11,
V V
Here dS is the measure of the surface area, and the equality
~13! follows from Gauss’theorem. 3k1p12!mk1p1S2~1,...,k1p11,Dk1p12;t!
E
The enumerated events combine together to create the ]
ceoqmuipvlaelteentraetveenotfschhaavnegteheofsatmheepprroobbaabbiilliittyymdeenassiutyre,mwkVe.cAans 1 Vd~k1p11! vk1p11(cid:149)]Jrk1p11
equate Eq. ~10! with the sum of contributions ~7!, ~12!, and
~13!, thus obtaining the hierarchy equations (k51,2,...) 3mV ~1,...k1p11;t! . ~16!
k1p11
066111-3
PIASECKI,TRIZAC,ANDDROZ PHYSICALREVIEWE66,066111 ~2002!
E E
It is then a question of inspection to see that on the right- k
1 (
hand side of Eq. ~16!, only the term l dvk11 dsˆ~sˆ(cid:149)vj(k11)!
j51
E E
(‘ ~k1p!! d~k11! d~k1p! 3u~2sˆ(cid:149)vj(k11)!fk11~1,...,j, ...,k,
p51 p! V (cid:149)(cid:149)(cid:149) V
k k1p rk115rj2ssˆ,vk11;t)/n~t!, ~22!
( (
3 j51 i5k11 Tv~i,j!mk1p~1,...,k1p;t! and its prefactor 1/l remains finite in the same limit. So,
E formallythehierarchyequation~18!atagiventimetreduces
k
( pointwise in the Grad limit to the Boltzmann-like hierarchy,
5 d~k11! Tv~j,k11!fk11~1,...,k,k11;t! S D
j51
] k ]
~17! 1( v fB~1,...,k;t!
]t j51 j(cid:149)]rj k
survives.All the remaining terms exactly cancel out. E
k
The hierarchy equations satisfied by the reduced distribu- 5 d~k11!( Tv~j,k11!
tions fk describing the annihilation dynamics thus read j51
S D
] (k ] (k (k 3fkB11~1,...,k,k11;t! ~k51,2,...!.
]t1j51 vj(cid:149)]rj2i,j i,j Tv~i,j! fk~1,...,k;t! ~23!
E
k The hierarchy ~23! propagates the factorization of the re-
(
5 d~k11! Tv~j,k11!fk11~1,...,k,k11;t!. duced distributions
j51
~18! )k
fB~1,...,k;t!5 fB~j;t!. ~24!
k k
Consider now Eqs. ~18!, supposing that the state of the sys- j51
temisspatiallyhomogeneous.Inthiscasethedistribution f
1 Hence, if the initial state is factorized, the whole hierarchy
doesnotdependontheparticleposition.Letusformallytake ~23! reduces to a nonlinear equation
the Grad limit
S D
E
] ]
s!0, n~t!!‘, n~t!sd215l215const, ~19! 1v fB~1;t!5 d2Tv~1,2!fB~1;t!fB~2;t!.
]t 1(cid:149)]r
1
~25!
where
E Equation~25!istheBoltzmannkineticequationcorrespond-
n~t!5 dvf ~v;t!. ingtotheannihilationdynamics.Inthefollowingsection,we
1 shall see that the formal Grad limit taken here, where n
!‘, is relevant for the description of the annihilation dy-
The fixed mean free path l introduces a relevant length namics at late times, even if the density n(t) decreases with
scale, so we pass to dimensionless positions, putting increase in time.
r5lx , j51,2,... . ~20!
j j B.Scalinganalysisofthehierarchy
Withthischangeofvariables,thecollisionoperator~9!takes The evolution of the annihilation kinetics shares a com-
mon feature with the Grad limit: The ratio of particle diam-
the form
etertomeanfreepathl51/(nsd21), whichisrelatedtothe
E
packingfraction,vanishesinbothcases.Tobemorespecific,
1
Tv~i,j!5@n~t!sd#d21 dsˆ~sˆ v !u~2sˆ v ! we perform a scaling analysis of the exact homogeneous hi-
l (cid:149) ij (cid:149) ij
erarchy equations and look for self-similar reduced distribu-
tions where the time dependence has been absorbed into the
3d@x 2n~t!sdsˆ#. ~21!
ij density n(t), with the velocities v and positions r renormal-
We conclude that the term on the left-hand side of Eq. ~18! ized by the typical ~root-mean-squared! velocity ¯v(t) and
involving the collision operators ~21! vanishes in the Grad meanfreepath,respectively.Wedefinethereducedvariables
limit for d.1, because the dimensionless parameter n(t)sd
v r
tendstozero.Notethatthisterminducesdynamicalcorrela- c5 and x5 ,
tions, hindering the propagation of the molecular chaos fac- ¯v l
torization. On the other hand, using the definition of Tv(j,k
11), we find that the term on the right-hand side equals with
066111-4
DYNAMICSOFBALLISTICANNIHILATION PHYSICALREVIEWE66,066111 ~2002!
E
1 scaling form
~¯v!25 v2f ~v;t!dv. ~26! S D
n~t! 1 k
n
f ~r ,v , ...,r ,v ;t!5 ˜f ~x ,c , ...,x ,c !.
k 1 1 k k ¯vd k 1 1 k k
For the one-body distribution, we therefore introduce the re- ~28!
duced function˜f , such that
1 The evolution equations of n(t) and kinetic energy den-
sity n¯v2(t) follow from integrating Eq. ~18! with weights
S D dv and v2dv , respectively, for k51. We obtain
1 1 1
n~t! n~t! v
f ~v;t!5 ˜f ~c!5 ˜f . ~27! dn
1 ¯v~t!d 1 ¯v~t!d 1 ¯v~t! 52v~t!n, ~29!
dt
d~n¯v2!
By definition, the moments of order 0 and 2 of˜f are con- 52av~t!n¯v2, ~30!
1 dt
strained to unity. Requiring that the k-body distribution
fk factorizes into )ik51f1(i) in the limit of infinite relative where the collision frequency vand kinetic energy dissipa-
separations between the particles, we consistently obtain the tion parameter aread
E
v~t!5n~t!¯v~t! dc dc dsˆ~2sˆ c !u~2sˆ c !˜f ~c ,c ,ssˆ!, ~31!
1 2 (cid:149) 12 (cid:149) 12 2 1 2
E
dc dc dsˆ~sˆ c !u~2sˆ c !c2˜f ~c ,c ,ssˆ!
1 2 (cid:149) 12 (cid:149) 12 1 2 1 2
a5FE GFE G. ~32!
c2˜f ~c!dc dc dc dsˆ~sˆ c !u~2sˆ c !˜f ~c ,c ,ssˆ!
1 1 2 (cid:149) 12 (cid:149) 12 2 1 2
S D
Equation ~32! is valid for general velocity rescalings. Defi- 2 11a
nition ~26! chosen here implies that the term *c2˜f1(c)dc in C511aln 11 2 v0t . ~35!
thedenominatorequalsunity.Thecoefficientamaybeseen
as the ratio of the kinetic energy dissipated in a typical col-
The corresponding time evolution is
lision normalized by the average kinetic energy, and is time
independent in the scaling regime. It is convenient to intro- F G
duce the internal ‘‘clock’’ C of the dynamics counting the n 11a 22/~11a!
numberofcollisions,suchthatdC5vdt. Withthisvariable, n 5 11 2 v0t , ~36!
0
Eqs. ~29! and ~30! integrate into
F G
n~t!5n exp@2C~t!# ¯v 11a (12a)/(11a)
0 5 11 vt . ~37!
¯v 2 0
0
and
Withoutknowingthedetailedformoftheone-particledistri-
¯v2~t!5¯v02exp@2~a21!C~t!#, ~33! bution function f , it is thus possible to conclude about the
1
time decay of n and¯v, which appear to obey algebraic laws
where the time origin with density n and kinetic energy
densityn ¯v2 hasbeenchosentocoincid0ewithC50. Knowl- in the long-time limit @n(t)}t2j and ¯v(t)}t2g, with j
0 0 52/(11a) and g5(a21)/(a11)]. The exponents jand
edge of the C dependence of n and ¯v allows one to relate gare consequently simply related to the unknown quantity
absolute time t to the number of accumulated collisions: a, forwhichaperturbativeexpansionwillbeputforwardin
From Eq. ~31!, we have
Sec. III before a numerical investigation in Sec. IV. More-
over, if the initial velocity distribution is of finite support
dC5v n ¯v 5v exp@2C~11a!/2#, ~34! ~i.e., vanishes outside a sphere of given velocity v*), ¯v ful-
dt 0n0 ¯v0 0 fills the bound ¯v<v* so that gis necessarily positive ~or
a>1). IntheframeworkofBoltzmann’sequation,itwillbe
where v[v(t50), so that showninAppendixAthatthequantityaisnecessarilylarger
0
066111-5
PIASECKI,TRIZAC,ANDDROZ PHYSICALREVIEWE66,066111 ~2002!
than 1. For the specific initial condition where all particles III.BOLTZMANNKINETICEQUATION
havethesamekineticenergyatt50, ¯v25^v2& istimeinde- This section is devoted to the analysis of the decay dy-
pendent, and Eq. ~37! implies that a51. From Eq. ~36!, we namics within the molecular chaos framework @13# of the
therefore obtain the time evolution for this situation, homogeneous nonlinear Boltzmann equation. No exact solu-
tion could be obtained, and our goal is to derive accurate
approximatepredictionsforthescalingexponentsjandgof
n 1 the density and root-mean-squared velocity.
5 , ~38!
n 11vt Before considering the kinetic equation obeyed by the re-
0 0
scaled distribution function, it is instructive to rewrite the
original homogeneous Boltzmann equation ~25! in the form
which is exact within the scaling theory. This relation a pri-
]f ~v;t!
ori holds in any dimension, except for d51 where the cor- 1 52n~v;t!f ~v;t!,
]t 1
responding initial velocity distribution is the symmetric dis-
crete bimodal function 6c, for which the scaling ansatz
with
underlying our approach fails ~see the discussion at the end
F E GE
of the present section!.
Inserting the scaling forms ~27! and ~28! into the first n~v ;t!5 sd21 dsˆ~sˆ vˆ !u~sˆ vˆ ! dv uv uf ~v ;t!,
1 (cid:149) 12 (cid:149) 12 2 12 1 2
equation of the hierarchy ~18! imposes the following con- ~39!
straint on the decay exponents: j1g51. This scaling rela-
tion may be simply obtained by elementary dimensional whereinthelastequation,theterminbracketsmaybecom-
analysis @4–6,9,10#, and may be considered as the compat- puted explicitly as a function of dimension ~it is understood
ibilityconditionofthehierarchywiththeself-similarscaling that the unit vector vˆ denotes an arbitrary direction!. For
12
solutions@12#.Itis,moreover,identicallyfulfilledbyexpres- our purpose, it is sufficient to notice that at all times, the
sions ~36! and ~37!. Under the constraint j1g51, the re- collision frequency n(v;t) of the population having velocity
maining equations of the hierarchy (k.1) turn out to be vremainsfiniteinthelimitv!0, providedthefirstmoment
compatiblewithEq.~28!withtheadditionalinformationthat of f exists. In this situation, Eq. ~39! implies that
1
the collision term on the left-hand side of Eq. ~18! decays f (v;t)/f (v;0) admits a finite limit for v!0, or equiva-
1 1
liket2g2dj, whereastheremainingtermsareassociatedwith lently, it may be stated that if the initial velocity distribution
m
a power 1/t. Given that g1dj511(d21)j>1, this colli- behaves like v near the velocity origin, this feature is pre-
sion term is asymptotically irrelevant except in one dimen- served by the Boltzmann dynamics. Previous works have
sionwhereitremainsofthesameorderasthedominantones shown accordingly that the scaling exponents j and gde-
(1/t). Wethereforerecovertheconclusionsobtainedbycon- pend on the exponent m@4,6,9,10#.
sidering the formal Grad limit, with distribution functions Making use of relations ~29! and ~30!, insertion of the
expected to obey a Boltzmann-like equation. This analysis scaling form ~27! into the Boltzmann equation leads to
points to the relevance of Boltzmann equation for d.1, a F S DS DG
12a d
point that is further corroborated by the numerical results 11 d1c ˜f ~c !
given in Sec. IV. 2 1dc 1 1
1
It is interesting to note that both Eq. ~38! and the relation E
j1g51 do not hold in 1D for discrete initial velocity dis- 5˜f ~c ! dc c12 ˜f ~c !, ~40!
tributions. In the symmetric situation of a bimodal distribu- 1 1 2^c12& 1 2
tion, the average kinetic energy per particle is conserved ~so
that g50), whereas the density decays as 1/At ~i.e., j where we have assumed an isotropic velocity distribu-
51/2 @1#!. The scaling assumption ~28! is consequently in- tion @˜f1(c)5˜f1(c)# and introduced the average ^((cid:149)(cid:149)(cid:149))&
correct in the specific situation of discrete distributions in 5*( )˜f (c )˜f (c )dc dc . ^c & is therefore the res-
(cid:149)(cid:149)(cid:149) 1 1 1 2 1 2 12
1D, but valid for continuous distributions @6#. To be more caled collision frequency. Once mhas been chosen, Eq. ~40!
specific, the scaling form ~28! implies that the collision fre- admits a solution for a unique value of a. We show in Ap-
quency scales with time, like v}n¯v. On the other hand, pendixAthat the inequality a.1 necessarily holds.
fromtheanalyticalsolutionofthebimodal6c situation@1#, Irrespectiveofa, thelargevelocitybehaviorof˜f maybe
1
we obtain v}n2¯v with¯v5c. This discrepancy is the signa- obtained following similar lines as in Refs. @4,5,9,14,15#: it
ture of dynamical correlations in the latter discrete case. is possible to integrate formally Eq. ~40! and cast˜f into
1
These correlations are responsible for the breakdown of Eq. S D
~28!, and in addition, violate molecular chaos. For continu- ˜f ~c! c 2[21d(12a)]/(12a)
ous velocity distributions, again in 1D, molecular chaos also 1 5
breaks down while the scaling ~28! is correct. As a conse- ˜f1~c8! c8 F G
quence,theexponentsobtainedattheBoltzmannleveldiffer E
from the exact ones ~see the discussion in the last paragraph 3exp 2 1 c ˜n~c9!dc9 . ~41!
of Sec. III!, whereas the relation j1g51 holds. 12a^c12& c8 c9
066111-6
DYNAMICSOFBALLISTICANNIHILATION PHYSICALREVIEWE66,066111 ~2002!
F G
In this equation,˜nis itself a functional of˜f1; ˜f ~c!5M~c! 11(‘ a S ~c2! . ~48!
E 1 n n
n51
˜n~c !5 c ˜f ~c !dc , ~42!
1 12 1 2 2
ThesepolynomialsareorthogonalwithrespecttotheGauss-
ian weight,
such that ˜n(c)/c goes to a finite limit for c!‘. We there- S D
fore obtain the large velocity tail d d/2
S D M~c!5 2p e2dc2/2, ~49!
2 c
˜f ~c!}c2[21d(12a)]/(12a)exp 2
1 a21 ^c & and the first few read
12
for c!‘. ~43! S ~x!51, ~50!
0
Inonedimension,werecovertheresultsofRefs.@4#and@9#. d
In Refs. @4,5#, an approximation was derived for a @or S ~x!5 ~2x11!, ~51!
1 2
equivalently,j52/(11a)] byassumingthatthelargeveloc-
ity behavior of ˜f1 could hold for the whole velocity spec- d2 d~d12! d~d12!
trum. In this picture, the power of c appearing on the right- S ~x!5 x22 x1 . ~52!
2 8 4 8
hand side of Eq. ~43! is equated to the exponent m
characteristicofthesmallvelocitybehavior~imposedbythe
The coefficients a follow from the orthogonality relation
initial distribution chosen, see above!, with the result *S (c2)S (c2)M(nc)dc}d . In particular,
n m nm
2 2d12m
a511m1d or j52d12m11. ~44! a152d^S1~c2!&5d2~12^c2&!50, ~53!
This prediction encodes the correct dependence on mand as obtained from the definition of rescaled velocities ~26!.
dimension(jincreaseswhenmordincrease!,andturnsout
The first non-Gaussian correction is thus embodied in a ,
2
to have an accuracy of order 10% when compared to the whichisproportionaltothefourthcumulant~kurtosis!ofthe
numerical results @10#. In the limit of large dimension, we
velocity distribution:
obtainfromEq.~44!j;12(2d)21, whereasKrapivskyand
Sire have shown that j;12d21(121/A2), also in the d2 d
frameworkoftheBoltzmannequation.Theremainderofthis a253 @^ci4&23^ci2&2#5d12^c4&21, ~54!
section is devoted to the derivation of a more precise value
for a, which furthermore coincides with the exact 1/d cor-
with c a Cartesian component of c. Upon truncating Eq.
rection for d!‘. ~48!atiafiniteordern,weobtainaregularvelocitydistribu-
Invoking the identity tion near c50. We consequently restrict our analysis to the
E S D casem50 ~thedependenceonmhasbeenconsideredinRef.
dccp d1c d ˜f ~c!52p^cp&, ~45! @10#!.
dc 1 It is also noteworthy that any truncation of Eq. ~48! at
arbitrary order n leads to a Gaussian high-energy behavior,
theenergydissipationparameteramaybegivenasthesetof incompatible with the result ~43! corresponding to an over-
equivalent expressions: populated tail with respect to the Maxwellian. However, it
S D
will be shown in Sec. IV that the difference between the
2 ^c cp& truncatedexpansion~48!andthenumericalvelocitydistribu-
a511 12 1 21 , ; p>0. ~46! tion becomes manifest far in the tail, where the distribution
p ^c &^cp&
12 1 has reached very low probabilities. Consequently, when the
moment involved in Eqs. ~46! and ~47! are evaluated from
Aparticularly useful relation between aand moments of˜f1 the truncation of Eq. ~48!, the accuracy of the result is ex-
followsfromconsideringthelimitc1!0 ofEq.~40!;weget pectedtobebetterforlowordersofpinEq.~46!.Hencethe
S D privilegedroleisplayedbyEq.~47!,whichisoflowerorder
a511 2 12 ^c1& . ~47! than any of the identities ~46!. In practice, upon truncating
m1d ^c & Eq. ~48! at order n, the n unknowns a,a , ...,a are com-
12 2 n
puted evaluating the corresponding moments appearing in
The ~infinite! family of relations ~47! and ~46! is equivalent Eq. ~47! and in n21 of the relation ~46!, among which it is
to the original integro-differential equation ~40!, and well convenient to retain the n21 even values of p (p50 ex-
suited to a perturbative analysis. To this end, a systematic cluded!. Truncation of Eq. ~48! at order n52 yields precise
approximation of the isotropic function˜f can be found by predictions for aand˜f , and already at Gaussian order, a
1 1
expanding it in a set of Sonine polynomials @16#: turns out to be very close to its numerical counterpart. Set-
066111-7
PIASECKI,TRIZAC,ANDDROZ PHYSICALREVIEWE66,066111 ~2002!
ting n50 ~or equivalently, n51 since a [0) in Eq. ~48!, tracted numerically from the exact dynamics (0.78560.005
1
we obtain immediately the zeroth-order approximation in Ref. @6#, and more recently 0.804 in @19#!.The difference
S D
between the exact exponents and those obtained assuming
A
2 2 molecular chaos is consistent with the conclusion of Sec. II:
a5a511 12 , ~55!
0 d 2 Inonedimension,thefactorizationofthetwo-bodydistribu-
tion˜f underlyingBoltzmannansatzisnotanexactproperty
2
which corresponds to of the distributions obeying hierarchy ~18!. On the other
hand, for d.1, the molecular chaos exponents are expected
2 2d
j5j0511a052~d11!2A2. ~56! tsoecbtieocno.meexact.Thispropertyisillustratedinthefollowing
For large dimensions, this estimation goes to unity, with the
exact 1/d correction computed in Ref. @9#, IV.SIMULATIONRESULTS
S D S D
Thenumericalresultspresentedinthissectioncorrespond
A
1 2 1 tothesituationm50, unlessstated.WerefertoRef.@10#for
j512 12 1O . ~57!
d 2 d2 the case of diverging (m,0) or depleted (m.0) velocity
distributions near v50.
This behavior turns out to be ‘‘universal,’’in the sense that
the mdependence does not appear at this order @10#.
A.Thenumericalmethods
Weshallalsobeinterestedinthenon-Gaussianfeaturesof
Wefollowtwocomplementarynumericalroutes.First,we
the velocity distribution, which we quantify by the fourth
cumulanta . The~cumbersome!calculationsatsecondorder solve the time-dependent homogeneous Boltzmann equation
2
by means of the direct simulation Monte Carlo method
in Sonine expansion are detailed inAppendix B. We obtain
~DSMC!, originally developed to study ordinary gases @20#.
d~2A223! This scheme, where a suitable Markov chain is constructed,
a2584d21d~62A2!, ~58! has been extended to deal with inelastic collisions @21,22#
and is easily modified to describe the situation under study
here, which does not conserve the total number of particles.
A
2 Restrictingtoaspatiallyhomogeneoussystem,thealgorithm
a5a1 a . ~59!
2 0 16d 2 is especially easy to implement, and may be summarized as
follows: among N initial particles having a given velocity
0
1Thae c)oarsrebsepfoonrde.inFgordden!si‘ty,eaxp;on2e(n3t2fo2lAlo2w)sd2fr1omirrejs2p5ec2t/i(v1e fdriosmtribthuetiosny,staempawirit(hi,aj)priosbcahboilsietynpartorpaonrdtioomna,latnoduvremu.oTvhede
2 2 ij
of m@10#, which reinforces the universal nature of large d. ~suitably renormalized! time variable is then incremented by
The correction to j carried by a behaves as 1/d2 in this the amount (N2uv u)21, where N is the number of particles
2 ij
limit, and does not affect the 1/d terms that still coincide remaining in the system before another pair is drawn, etc.
with the exact behavior ~57!. Both predictions ~55! and ~59! This scheme provides the numerical exact solution of Eq.
are such that a.1, which is required to obtain a normaliz- ~39!,andallowstotestthevalidityoftheapproximationsput
able distribution in Eqs. ~41! and ~43!. forward in Sec. III @essentially, truncation at second order in
The second-order expansion considered here may be im- expansion ~48! supplemented with calculations performed at
proved by consideration of higher-order Sonine terms and linear order in the fourth cumulant a ~seeAppendix B!#.A
2
inclusion of nonlinear terms in a . In the related context of preciseanalysisofthelatetimedynamics~andespeciallythe
2
inelastic hard spheres, the limitation of working at linear computation of velocity distributions! suffers from the con-
order in a with neglect of Sonine terms of order n>3 may comitant low number of particles left, and the statistical ac-
2
befoundinRef.@17#.Alternatively,keepingnonlinearterms curacy is improved by averaging over independent realiza-
ina andneglectingagainSoninetermsofordern>3 leads tions.
2
to multiple solutions.Astability analysis is then required to The second numerical method ~molecular dynamics @23#!
determine which one is stable, as discussed in Ref. @18#. consists of integrating the exact equations of motion for an
Here,thevalueobtainedfora isquitesmall~seebelow!. assembly of spheres confined in a ~hyper!cubic box with pe-
2
Ourapproximateexpressionsareaccuratewhencomparedto riodic boundary conditions. This route assesses the validity
the full numerical solution of Boltzmann’s equation, so that of the approach relying on the homogeneous Boltzmann
wedidnotcalculateanyhigher-ordercoefficients,nordidwe equation,butdoesnotofferthesameaccuracyasDSMC,nor
consider nonlinear terms in a . The existing literature re- the possibility to follow the evolution over comparable
2
ports, within the Boltzmann framework, numerical expo- times.Inparticular,oncethemeanfreepathl, whichgrows
nents in 1D, which are in excellent agreement with our pre- rapidly as tj, becomes of the order of the box size L, the
dictions, already at zeroth order. For the case m50, subsequent evolution suffers from finite size effects and
expressions ~56! and ~59! give j.0.773 and j.0.769, should be discarded. When l.L, the algorithm used is un-
0 2
whereas the numerical result obtained in Ref. @9# is j able to find collision events for those particles that make
.0.769. These exponents are close to their counterparts ex- more than one free flight round on the torus topologically
066111-8
DYNAMICSOFBALLISTICANNIHILATION PHYSICALREVIEWE66,066111 ~2002!
FIG.1. Inversedensityn /n versusn ¯v /(n¯v), ford52. The
0 0 0
dotted line has slope j50.87 @see Eq. ~56!#. Initial number of
2
particles:53106 ~withafurtheraverageover103 independentrep- FIG. 2. Plots of (n /n)21 @upper curve corresponding to
0
licas! for DSMC and 23105 for MD. In both cases, the initial DSMC, compared to its MD counterpart ~crosses!# and (¯v /¯v)2
0
velocitydistributionisGaussian(m50), andtheinitialconfigura- 21 ~lowerdashedcurveforDSMC,circlesforMD!asafunction
tion used for MD is that of an equilibrium hard disk fluid with oftime.Thedottedlineatshorttimeshasslope1.
packingfractionf50.1~chosenlowenoughtoavoidfinitepacking
effects!. those reported in the literature (j.0.89 in Ref. @4# and j
.0.87(5) using a multiparticle lattice gas method @24#!.
equivalent to the simulation box, which causes a spurious The time evolution of inverse density and inverse typical
slowing down of the dynamics. It is then tempting to reduce velocity square is shown in Fig. 2, where considering n0/n
l by increasing the particle diameter s, but then finite ~and 21 and¯v2/¯v221 instead of n /n and¯v2/¯v2 allows one to
0 0 0
necessarily transient! density effects—incompatible with the probe the short time behavior. Unless stated, the initial ve-
scaling assumption ~28!—may also arise if the packing frac- locity distribution is an isotropic Gaussian. From Eqs. ~29!
tion f, proportional to nsd, is not low enough. Simulating and~30!,nand¯v2 evolvelinearlywithtforvt!1 @seealso
0
explicitlythelimitofpointparticles,theDSMCschemecon- Eqs. ~36! and ~37!#; the same holds for inverse density and
sidered here is free of this defect.The initial number of par- inversetypicalvelocitysquared,whichisindeedobservedin
ticlesconsideredinMDneedstobelargetoallowthesystem Fig. 2. MD and DSMC results superimpose, except at late
to enter the scaling regime before finite size effects become times where MD suffers from the slowing down discussed
dominant, we considered systems with N5(53104) previously. For both numerical methods, the scaling relation
–(53105) spheres initially ~compared to N5106–108 in j1g51 iswellobeyed,inprincipleatlatetimesonly,inthe
DSMC!. scaling regime. Special combinations of n and ¯v can, how-
ever, be constructed with the requirement to match the short
B.Dynamicscalingbehavior timeevolutionwiththescalingbehavior.Oneofthesequan-
tities is displayed in Fig. 3, with a resulting scaling regime
The results of two-dimensional DSMC and MD simula-
extending over more than ten decades in time. In Fig. 4, we
tionsareshowninFig.1,wherethequantityonthexaxisis
notonlytestthevalidityofthetheoreticalscalingexponents,
expected to scale like real time t from the scaling relation
j1g51. This log-log plot is a direct probe of the exponent but also the full time dependence as predicted by Eqs. ~36!
j, from the slope measured. Both MD and DSMC methods and ~37!. In order to improve the agreement between theory
and simulation ~which holds over more than six orders of
give compatible results, with the possibility to follow the
magnitude in time!, the system has been left to time to enter
dynamics over a longer time interval in DSMC. The depar-
ture observed for n /n.200 corresponds to the slowing the scaling regime. The time origin t50 has been chosen
0
down of MD resulting from finite size effects ~see Fig. 2 when 80% of the particles originally present have disap-
below!. The theoretical predictions at zeroth and second or- peared.Thecorrespondingreferenceconfiguration~withsub-
der are very close (j50.872 and j50.870), and in excel- scripts 0! thus differs from the ones considered previously.
0 2
lent agreement with the simulation results over several de-
cades. On a similar graph as Fig. 1, the kinetic C.Velocitydistributioninthescalingregime
‘‘temperature’’¯v2 exhibitsapowerlawbehavior~notshown! Inordertounderstandthereasonsforthegoodagreement
withanexponent22g, ingoodagreementwiththetheoret- between our theoretical predictions and the simulations, we
ical prediction (g.0.13 in 2D!. Moreover, the exponents now consider the velocity distribution, restricting to Monte
obtained analytically and numerically are compatible to Carlo results ~leading to similar conclusions, MD is much
066111-9
PIASECKI,TRIZAC,ANDDROZ PHYSICALREVIEWE66,066111 ~2002!
FIG. 3. Plot of log10@n0/n21#1log10@n0¯v02/(n¯v2)21# on a FIG. 5. Probability distribution function˜f1(ci) of a given Car-
logarithmic time scale (d52). DSMC results are represented by tesian component c of the rescaled two-dimensional velocity c.
i
the continuous curve, while the crosses correspond to MD. The Thetime-independentdistributionobtainedinDSMCsimulationsat
dashed line has slope 2. At late times, the quantity displayed is latetimesiscomparedtotheGaussianMandtheSonineexpansion
expected to behave as 2log10(v0t) from the scaling relation j1g truncated at n52, with the fourth cumulant given by Eq. ~58! (a2
51. At short times, the same behavior is observed, for a different .0.109ford52 andm50).Alldistributionshavevariance1/2,so
reason@seeEqs.~36!and~37!#.TheultimateMDslowingdownis that^c2&51. Theresultshavebeenobtainedbyaveragingover50
againvisible. replicasofasystemwithinitiallyN5403106 particles.
theSoninepredictionpushedatsecondorder.Theagreement
more demanding on CPU time and does not allow to inves-
is remarkable, and it is also observed that the Gaussian ap-
tigatedetailedfeaturesof˜f withthesameaccuracy!.Aftera
1 proximation is already close to the asymptotic rescaled ve-
transient where the probability distribution function ˜f1 locity distribution, which is quite surprising in a kinetic pro-
evolves with time, a well-defined scaling regime is reached cess extremely far from equilibrium, with furthermore no
with a time independent˜f (c) shown in Fig. 5 together with conservation laws. Given that our perturbative analytical
1
work relies on the calculation on low order moments of˜f ,
1
this explains the accuracy of the zeroth-order predictions a
0
and j. From Eq. ~43!, we expect the differences between
0
the Sonine expansion and the numerical distribution to be-
come visible in the high-energy tail, which is confirmed in
Fig. 6.As predicted,˜f is overpopulated with respect to the
1
Gaussian, and displays a high velocity tail of the form ~43!
~see the inset of Fig. 6!.
D.Evolutiontowardstheasymptoticsolutions
Before the scaling regime is attained, ˜f is time depen-
1
dent, as shown in Fig. 7, where the distributions at different
timeshavebeenrenormalizedbyMtoemphasizethebuild-
ing up of non-Gaussianities.The evolution towards the scal-
ing solution 11a S can be observed. With respect to the
2 2
Gaussian,˜f is at all times overpopulated both at large and
1
smallvelocities~whichmayberelatedtothepositivesignof
a forthelattercase!;normalizationisensuredbyanunder-
2
population at intermediate velocities. Figures 5–7 show that
theindirectmeasureofa throughthenon-Gaussiancharac-
FIG. 4. Time dependence of n ~lower curve! and ¯v2 ~upper 2
curve! obtained in Monte Carlo, compared to the dashed curves ter of˜f1/M agrees with the theoretical prediction, but it is
correspondingtothetheoreticalpredictions~36!and~37!,wherethe also possible to compute directly a2 in the simulations
energy dissipation coefficient ais calculated at second order in throughitsdefinitionasafourthcumulant@Eq.~54!#.Itturns
[email protected].~59!#. out that both methods are numerically fully compatible.
2
066111-10