Front formation in a ballistic annihilation model

PHYSICA ELSEVIER

Physica A 229 (1996) 515-529

Front formation in a ballistic annihilation model Jaros~aw Piasecki a, l, Pierre-Antoine Rey b, 2, Michel Droz b. 2, * a htstitute o f Theoretical Physics, Warsaw UniversiO', Ho2a 69, P-O0 681 Warsau'. Poland b Departenwnt de Physique ThOorique, Universitk de GenOve, CH-1211 Geneva 4, Switzerhmd Received 18 January 1996.

Abstract

We study a simple one-dimensional model of ballistically-controlled annihilation in which the two annihilating species are initially spatially separated. The time dependent properties of the annihilation front are exactly derived. It is shown that the front wanders in a brownian fashion around its average value.

PACS: 05.20.Db, 05.40.÷j Keywords: Nonequilibrium statistical mechanics; Ballistic annihilation; Front

1. Introduction

The study of the kinetics of diffusion-controlled or ballisticaly-controlled annihilation provides a nice playground on which the tools of nonequilibrium statistical mechanics can be tested [1,2]. Particularly, in low dimensions, the kinetics of these systems is governed by the fluctuations and any mean-field like (or Boltzmann-like) approximation will not give a satisfactory description of the properties of the system. A lot of efforts have been expended in studying the diffusion-controlled case, both for the situation in which the reactants are uniformly distributed in space [3] and for the one in which they are initially spatially separated, leading to reaction-diffusion fronts [4, 5]. In both cases, the long time behavior of physical observables is characterized by power law. The associated exponents have some universal properties [2]. More recently, the case of ballistically-controlled annihilation has been investigated for spatially homogeneous conditions [6]. In particular, exact results have been obtained for one-dimensional, single species models by Piasecki [7] and Droz et al. [8, 9]. * Corresponding author. Fax: 4122 702 6870; e-mail: [email protected]. t Work partially supported by the Committee for Scientific Research (Poland), grant 2 P302 074 07. 2 Work partially supported by the Swiss National Science Foundation 0378-4371/96/$15.00 @ 1996 Elsevier Science B.V. All rights reserved PI[ S 0 3 7 8 - 4 3 7 I ( 9 6 ) 0 0 0 3 6 - 2

516

J. Piasecki et al./Physica A 229 (1996) 515-529

In these models, one considers point particles which move freely, with a given velocity. When two particles collide they instantaneously annihilate each other and disappear from the system. The system with only two possible velocities + c or - c has been studied in a pioneering work by Elskens and Frisch [10]. The case of an arbitrary discrete velocity distributions has been investigated by Droz et al. [9]. An exact equation for the survival probability until the time t of a particle moving with velocity v was derived, allowing to compute the density and the time dependent velocity distribution in the asymptotic regime t ~ ~ . For a symmetric three velocities distribution, different kinetic regimes were found as a function of the initial fraction of particles at rest. Such processes can model several physical situations as a recombination reaction in the gas phase or the fluorescence of laser excited gas atoms with quenching on contact (the one-dimensional aspect can be obtained by working in a suitable porous media [11]) or, the annihilation of kink-antikink pairs in solid state physics [12]. A related but different class of problems is the case in which the two species (called A and B), are initially separated in space and annihilate on contact. Such a process can model the situation in which chemical species incorporated in a gel move ballistically under the action of a drift [13]. As the two species cannot penetrate one into the other (as they annihilate on contact), a well defined reaction front is formed. We aim at computing exactly the time dependent properties of this front. The paper is organized as follows. In Section 2, the model is defined. In Section 3, we compute the average number of collisions that have occured before time t. An exact analytic expression for the probability density to find the front at a given point and a given time is derived. In Section 4, it is shown that the front makes a random walk around its (time dependent) average value. Possible extensions of this work are discussed in the conclusion.

2. The model We consider a one-dimensional system formed of two species of particles. Initially, particles A are spatially randomly distributed in the region (-cx~, 0) and particles B are spatially randomly distributed in the region (0, oo). The positions of the particles obey a Poissonian distribution. The initial linear average densities for the A and B particles are respectively pA and pS. The velocities of each particle is an independent random variable taking the value + c with even probability. Particles of the same kind suffer elastic collisions. When two particles A and B meet, they annihilate. Thus, practically the A particles with velocity - c and the B particles with velocity + c move freely. Accordingly, the relevant part of the dynamics concerns the A particles with velocity + c and the B ones with velocity - c . Let (Yl, Y2,..., Yk.... ) be the initial positions of the A particles with velocity +c and ( x l , x 2 . . . . . xk . . . . ) the initial positions of the B's with velocity - c (see Fig. 1). The relative velocity between the A and B particles being 2c, the pair of particles initially at ( x k , y k ) will collide at time tk = ~(xk - Y k ) .

J. P i a s e c k i et a l . / P h y s i c a A 2 2 9 ( 1 9 9 6 ) 515 529

B particles

A particles yk

517

,

x2

Y2

Xk

=

Yk-:

"" •

yx 0

Y3

Xl

X3

.."

Xk_

1

Fig. 1. T y p i c a l initial c o n d i t i o n . T h e particles are l a b e l e d such that Yk < " ' " < Yl < 0 < xl < . . .

< xt.

The position at which this collision will take place defines the position of the annihilation front. Thus, the position of the front at time t is ½(xl - Yl ), if 2ct <~xl + v (i.e. if no particle have yet collided), ½(x2 - Y2), if xl + Yl <~2ct <<,x2 + Y2 (i.e. only the first particles on the right and on the left have collided), and so on. The dynamics in itself is purely deterministic, the only stochastic aspect comes from the initial conditions. The properties of the front are completely defined by the probability density p(X; t) to find the front at the point X, at time t. Moreover, another quantity of interest is N(t), the number of particles that have been annihilated until a given time t.

3. Analytical solution We first investigate the behavior of N(t). For this purpose, we introduce the probability P(k; t) that exactly k collisions have occurred within the time interval [0, t]. For a given initial configuration, at least k collisions have occurred before the time t, if 2ct >~x~. y~., and at most k collisions have occured before t, if 2ct <~Xk+l Yk+l, with xa. < xk+t and Yk > Yk+l. Hence averaging over the initial conditions, we find:

P(k; t) = (0 (2ct - [xk - Yk]) 0 ([xk, 1 - Yk+l] - 2ct)} (O(2ct - [xk

Yk])) - (O(2ct - [x~+. - y k + l ] ) ) ,

(1)

where 0 is the usual Heaviside function and the brackets denote the average over the initial conditions (see Appendix A). Thus: (2}

N(t) = k k P ( k ; t ) , k=l

which, using (1), reduces to :x~

X(t) = Z(O(2ctk

(3)

[xt. - y , ] ) ) .

1

Thus, using the results of Appendix A, N(t) reads: 0

N(I) : ~ k= I

j

dyk pA

,?,c

dxk . / ~ e p'l.,',,I

(]Ykl(k -"A)k-'l )'

0 1

x e v~'~ ''(xkpB)k- O ( 2 c t - [ X k - - Y k ] ) (k

1)~

-

J. Piasecki et al./Physica A 229 (1996) 515-529

518

Noticing that the sum over k can be exchanged with the two integrals, we find I (where a = ~(xk - yk) and b = ½(xk + Yk))

N(t) = 2 pApB

da

0

db e-"(P%P%e6(P~-P%Io 2V/pAp ~ (a 2 - b2) ,

--a

I0 is the modified Bessel function. The integral over b is calculated in Appendix B and the integration on a is straightforward. We obtain

N(t) - pA + p--------~ ct

2(pA + pB)

"

In terms of the dimensionless time

z =--ct(p A + pB) and the parameter

4 pApB 2-- (pA + pS)Z ' we have: N(~:): 5 z-~ After a transient regime, N(T) grows linearly with time. The front position is defined as the point where the next collision will take place. Thus the probability density # (X; t) to find the front at the point X, at time t is given by

X, ÷ yl)O(Xl --Y, _ 2 c t ) l

+

=

6

_xk+j÷Yk+l 2

O(2ct--[xk--yk])

×0([xk+j + Yk+l] -- 2ct) >

xk+t +2 Yk+~)

2

(6)

519

J. Piasecki et al./Physica A 229 (1996) 515 529 The first term in the right-hand side reads simply

16(X

x l +2y l ) ) = p A p B

dxl/ 0

dv'e-lyx'e-v~v6 ~

X

-v, +2 y l

--::,c

__ 2pap" pA + t78

[O(X)e

z"~x +o(-X)e2V'x].

(7)

To calculate the second term in the right-hand side, we remark that one has to average over the position xk, xk+l, )'a and Yk+l. We shall then first average over particles k + 1, to obtain an expression containing only the positions xk and ya.. To do this, we need to calculate the following integral

S =pB

pA

dxk+ I va

dy/,-+l

a X

xa.+l+)'k+l 2

exp[-- pB

xa ÷I + P4 Ya -1 ],

oo

I I putting a = ~(xa-+l -Yk+J ) and b = g(xa-.i + Yk+l ), we have:

S = 2 pApB

/

da

~c

db O ( a + b - x k ) O ( a - b + y k ) d ( x - b ) -- oc

p")a + (pA _ p")b]

× exp [ (pA +

The integral over b is straightforward, and the integral over a gives after some algebra

2pAp" [o(Xk+Yk S - t7.4 + p. 2 2

X ) exp[--(pA+p")xk+2pAX] exp

[(pA +

)y~ --

2 p" X]

.

Hence, the second term in the right-hand side of Eq. (6) becomes

k=lL/{ DA2pApB -}- p"

[0 ( X k +2y k

xk +

X ) eO''+ve)v'92p'v .

x~ + yk ) /

) xO(2ct - [x, - Y~-]) ) • Now using the relation (see Appendix A)

da "dbexp[_(pA + pB)a + (pA

Z ( f ( x k + y~.,xk - Yk)) = 2 pap. k= I

0

pB)b ]

a

x f(2b, 2a)lo (2x/pAp" (a2 - b2)) ,

(8)

J. Piaseckiet al./Physica A 229 (1996) 515-529

520 one finally finds: /~(X;t) -

2 /+: pB pA

[O e-JX + O(-xSX]

+2 pA pB

/J da

0

x

{

db exp [_(pA + pe)a + (pA _ pS)b]

--a

-6(X-b)+

pA + p~----~O ( b - X ) e x p [ - 2 p A ( b - X ) ]

+O(X-b)exp[-2pe(X-b)]]}

(9)

In terms of the dimensionless variables z and 2, and the dimensionless coordinate

=_(pA + pS)X ' we can write pA + pe

.

U (~; r) -- _ _ e2~ ¢

2e-I¢1

with

p.4 _ pB pA + pB " Note that the function between the curly brackets is an even function of ¢; we shall call it/J(¢;z). The rest of the paper will be devoted to the study of this function.

4. Results and interpretation Starting from (10) and using the integrals quoted in Appendix B, it is possible to verify that y ( x ; t ) is correctly normalized to 1. Moreover, we can also calculate its different moments. However, for this purpose, it is simpler to consider the Fourier transform of y. Let us define the function fi(p; z) as: fi (p; z) = :

= ~

dXe-(P%P~)PXe-(p:-pB)XI2(X; z) d@-P~/~ (~; z).

(11)

521

Z Piasecki et al. IPhysica .4 229 (1996) 515-529

Here again, this integral can be calculated using the formulas quoted in Appendix B. One finds:

l [cosh(

sin

It is now easy to calculate the different moment of #, using the fact that

(~n)=

.J

l

dX[(pA +pB)x]n~(x;t)=(_l)n [?.~/~(p; F ?n ~ r) p : - ~ .

PC

(13)

Noticing that K 2 ~- ,~ 1, it is straightforward to check the normalization. The first moment ~(z) = (41 is given by =

2K [1 + N ( r ) ] , ~(r) = K [~ + z - ~ 1 (1 _ e _ a r ) l = -~-

(14)

and the variance is

--

4 2 3 5 2 [ 1 _2r(l_2)]e_2~ )2 2 + ~ -- T + 2 T + 2-- ~ 1 - 2e_4r 4

(15)

If pa¢pA, and in the long time limit, the front moves with a velocity proportional to the density difference. The case 2 -- 1 is of particular interest, because it corresponds to the symmetric case pA=pB. As expected: ~(z) ~ 0

(16)

{~s(r) = 2 [1 + N(z)] .

(17)

but

In the long time limit, we see that ~x r: the front moves essentially as an unbiased random walker. Thus, we can expect /~ to be essentially a gaussian. In fact, we shall prove that for any 2, the long time limit of # is a gaussian. Starting from the definition of an arbitrary moment o f / t (Eq. (13)), we may write for long time: ( ~ ) = ( - 1 ) ' 2 e -~

+ O(r "-1 ), k

~0

using

8p" we obtain (¢n) = (~cr)n + O(z.-I ).

Jp . . . .

(18)

J. Piasecki et al./Physica A 229 (1996) 515-529

522

If 2 < 1, K # 0 and we know the leading term for every moment of/~. If 2 = 1, the situation is different and will be considered later. It is known that for a gaussian distribution with a nonvanishing mean ~(~) and a variance ~2s(O, one has: ,/2 1 (¢2's(r)'~' (in)gauss = n!~n(z) ~ (n - 2k)!k! \ 2~2(~) J '

(20)

k=0

In the long time limit: 2~2(z) -- O

,

(21)

and (~n)gauss = (K'[) n ~ - O ( T n - - I ) . Thus, when 2 < 1 and for sufficiently long times (i.e. when the corrections to the leading term become negligible), the probabilty density # takes the gaussian form: p(~; r) o( exp [ ({ _~_~)2] •

22v

j •

(22)

The front wanders around its mean value with an amplitude of order of magnitude x/~, like a biased random walker with a diffusion coefficient of i , The 2 = 1 case is different because the amplitude of the zn term vanishes. However, we can still prove that for large time, the probability density # is a gaussian. More precisely, we shall give an upper and a lower bound, both of gaussian form. For this purpose, we shall come back to the Fourier transform fi(p; ~) defined by Eq. (11). It can be shown (see Appendix C), that the inverse Fourier transform of fi(p, z) is given by [ ( ) sinh (lx/TT~r) f i ( ~ ; v ) = 2e -I~1 ~ eosh x / T - ~ + v/l+)0

×

-2

sinh (x/1 + 2 z - I~]) +

/

dr cos(r~) [ 7~

1 +r 2

- 2,0(~ - I ~ l ) e - ~

c o s h ( v ' l + ,t, - I¢1)

( sinh ~ z

x/1 + ) ~

)

+

cosh ( 2x/Tf T ' -r2-c) 2X/~r

2

.

(23)

0

It is now easy to find an upper and a lower bound for ft. Indeed, for ]~] <~n/2x/~, one has:

cos(r~) cos (r~) 1 + ~ ~< ~1 + r ~< cos(r~)

J. Piasecki et al./ Physica A 229 (1996) 515-529

523

as r E [0, V~]. The remaining integral can be calculated exactly:

2 f dr cos(r~) rt ,

sinh (2V/-2~-r2"c) ÷

0

1' (V/)" (z2 - ~2))

( X/2('c2

V~7- r 2

32))

(24)

yielding, then two bounds for #. Two conclusions can then be drawn. First for ), the asymptotic time behavior of (24) is e 1~

1,

2~T e s: 2 / 2r

1 and r large, we have: and then, for { ~<2rt e ~: : 2~

2

e _~2,2r

(25)

In the long time limit, the probability density/~ becomes the probability density for a random walk, with a diffusion coefficient of 1_ 2" The second consequence concerns the case ~ = 0, for any value of ).: we have, in the long time limit:

fi(0; r) o(

{ e(l-,/~,)r xF - -

1

if 2 < 1,

(26)

if 2

(27)

z

1.

These results are not really surprising. Indeed, it can be verified that ~ = 0 is a maximum of I~, for 2 = 1 but not for 2 < 1. It is useful to get an image of what is going on by plotting the function/~ for various values of 2. However, for this purpose, both Eqs. (10) and (23) are inadequate for technical reasons. Eq. (10) is not well suited because it contains a double integral which requires a lot of CPU time and Eq. (23) because the integral diverges for r ,~, leading to numerical difficulties. To avoid these problems, we have obtained a different integral representation, starting from (10), and using the following relation:

[(~2

(')2 I0 (V/)~ (,[2 __ ~2)) = j[ 10 (~). (,C2:2)) -~f121

.

This relation follows directly from the integral representation for 10 1o

_ %

=

+

c

+

(28)

524

J. Piasecki et al./ Physica A 229 (1996) 515-529

o 8

mgiiTRN

l

-10

°o"

10 Fig. 2. Plot of #(~;z), for pA= 1 and ps= 1 (i.e. 2 = 1). where C is, for example, the circle of radius 1 centered in z = 0. Inserting (28) in (10), integrating by parts twice and using that ~2

0/~2

e-~-II~-~l = 2e-~6(/~ -

~)'

we get

lj I

/~ (~; "c) = 2 { e -I~1 + ~ o de e -~

e I~+~1)

I,

× ~

7--~--~ -I

-- 20(~ -- {)e ~-2~ -- 20(~ q- ~)e 2c~-~] / " ,.i

)

Performing the integral over the two theta functions, rearranging the terms in the remaining integral and using the integrals of Appendix B, leads to: /2(~;z)= Re-I:1-'

cosh ~ + - - 2 r

+

x/l+2

-20(~ - ,~,)e-~{ sinh(~ - ~) - /

d~ sinh(~ - I~,)

J. Piasecki et al./Physica A 229 (1996) 515 529

525

ZO

0.4

15

Fig. 3. Plot of/~(~;z), for p A : 1 + ½x/3 and p S : 1 - ½x/3 (i.e.),

0.5).

We have not been able to simplify further this expression and in particular to perform analytically the integration. However, this form can be easily integrated numerically. For example, the function/t(~; 3) is plotted on Figs. 2 and 3 for ). = 1 and for 2 = 0.5, respectively, pA > pB. One sees that for ), = 1, the probability density /~ spreads out and approaches a gaussian for long times. For the assymetric case )~ = 0.5, one sees clearly the drift of the front which becomes linear in the long time regime.

5. Conclusions Using a combinatorial analysis approach, we have been able to derive exact results for the number of time dependent collisions, the probability density # (X; t) to find the front at point X at time t and its first two moments. Although simple, this nonequilibrium model leads to nontrivial behavior as far as the dynamics of the annihilation front is concerned. Several extensions of the present work are possible. One is to consider a similar model in which the velocities are not restricted to be bimodal, but obey a more general (discrete or continuous) distribution. Another interesting aspect is given by the generalization of this model to higher dimensions. In this case, some A particles can invade the right part of the system without being annihilated (and vice versa for the B's). A new definition of the front in term of the annihilation rate should be introduced and different time dependent behavior may occur. These more complicated problems are under investigation.

6. Appendix A The average over the initial conditions is taken as follows: the B particles are distributed to the right of the origin, with a Poissonian distribution of density pS. Thus

J. Piasecki et al./Physica A 229 (1996) 515 529

526

the probability to find the first particle at position xl is

KI(Xl,pB) ~ e-Pexl pB . The probability to find k - ! particles between 0 and xk and a particle in xk is

rtk(xk,pB) = e_p%k (xk pB) k-I p8 . (k - 1)! The average o f a function depending only on xk is thus

(f(xk)) = f dxkf(xk)e -p'x~ (xk pS) k-I p8 (k - 1)!

d

0

(A.1)

"

Now, the average o f a function depending both on xk and xk+l is

dx k

(f(xk,Xk+l)) = 0

dxk+lj(Xk,Xk+l)e-p x~+, ( k -

1)t

(pC)2.

(A.2)

x~

Indeed, the probability to find exactly k - 1 particles between 0 and xk, one in xk and another in xk+l, but none between is

e-PS(X~+x-X~) pB e--penxl, (xk PB) k-I pB, (k -

1)!

as e -p%~+'-xk) p8 is the probability that the interval [xk,xk+l] is empty. Similar considerations apply for the A particles. Now, we propose to prove the formula (8). o~

~_~(f(Xk + yk,Xk

Yk))

k--I ~c

0

Xk pA (_yk)] k-I [ ( k - 1)!] 2

k--1

0

--oc

×f(xk + Yk, Xk -- Yk )

(A.3)

=2pApB j d a j d b O ( a + b ) --¢X3

O(a-b)e[-(p A + pB)a + (pA _ pS)b ]

- - OG

(A.4) oo

H

= 2 pap, f d a / d b e x p [ - ( p A + pe)a+(pA - p')b]f(2b, 2a) 0

--a

0 .

(A.5)

J. Piasecki et al./Physica A 229 (1996) 515 529

527

where in (A.4), we put 2b = xk + Yk and 2a = xk - Yk and use the series representation of I0.

Appendix B In this appendix, we quote two integrals that we often use in our derivation. The first is used, for example, in deriving the formula for N(t) (4), or in calculating the different moments of/~ without introducing the Fourier transform. I,,

),

(,',,

2' J d.,,,",'io .1

0

: sinh

~y~/fl~+fl~)

(B.2)

(see for example [14, Eq.(8)]). In calculating (28), we used the previous integral and also the following i •

I, (V/).('t"2 _ ~ 2 ) )

d~cosh(~K)2r

V/2(.c2 _ ~2)

=-

cosh(z) - cosh(Kz),

(B.3)

0

which can be found in [14, Eq.(6)].

Appendix C In this appendix, we shall prove formula (23), starting from Eq. (12). Remembering the definition (11) of fi(p; z), we may write ioo

K

-ioc

K

sinh( ~ r ) The integrant has two poles ( p = 4-1) and a cut on the imaginary axis between -ix/~. and i v/2~. For ~ > r, we can close the contour to the left, without changing the value of the integral. We get then a contribution from the pole p = - l and eventually from the cut (depending on the sign of •). Nevertheless, as cosh (x) + sinh (x)/x is an even function, the contribution of the cut vanishes and we are left with fi(~;r)=2e

-~-'

~[

cosh

(

~ z

) sinh(~/ +

,/1+2

'

~>~

~c,~

528

J. Piasecki et al./Physica A 229 (1996) 515-529

For ~ ~
,

~ < -z.

(C.2)

For ~ C ( - r , z), the situation is more complicated. First we write

cos

V/p 2 + 2

1

~

1+

+

1

For the first exponential, we can close the contour to the left and to the right for the second. Now, for both integrals, the integrant is not an even function of v / p 2 + 2 and one shall get a contribution of the cut. Integrating, one then eventually arrives to the formula /J(~;r) = 2e -~

e - l~T~;~ (1

--~) cosb(~)

drc°s(r~)[sinh( l+r 2

+2/ o

|

,

131 < ~.

(C.3)

Combining (C. 1)-(C.3) together gives (23).

Acknowledgements J.P. wishes to thank Dr. M. Droz for hospitality at the Department of Theoretical Physics of the University of Geneva where most of this research has been done. References [1] S. Redner, in: Proc 2nd Internat Colloquium on Quantum Field Theory and Stochastic Processes, (World Scientific, Singapore 1996), to appear; P.L. Krapivsky, S. Redner and F. Leyvraz, Phys. Rev. E 51 (1995) 3977. [2] M. Droz, Some aspects of pattern formation in reaction~tiffusion systems, in: Diffusion Processes: Experiments, Theory, Simulations, Lecture Notes in Physics, Vol. 438, Springer, Berlin, (1995) pp. 105-115. [3] D. Toussaint and F. Wilczek, J. Chem. Phys. 78 (1983) 2642. [4] L. G~ilfi and Z. Rficz, Phys.Rev A 38 (1988) 3151.

J. Piasecki et al./ Physica A 229 (1996) 515 ,529

529

[5] S. Cornell, M. Droz and B. Chopard, Phys. Rev. A 44 (1991) 4826. [6] E. Ben-Naim, S. Redner and F. Leyvraz, Phys. Rev. Lett. 70 (1993) 1890. [7] J. Piasecki, Phys. Rev. E 51 (1995) 5535. [8] M. Droz, P.-A. Rey, L. Frachebourg and J. Piasecki, Phys. Rev. E 51 (1995) 5541. [9] M. Droz, P.-A. Rey, L. Frachebourg and J. Piasecki, Phys. Rev. Lett. 75 (1995) 160. [10] Y. Elskens and H.L. Frisch, Phys. Rev. A 31 (1985) 3812. [ll] R. Kopelman, Sciences 241 (1988) 1620. [12] M. Bfittiker and T. Christen, Phys. Rev. Lett. 75 (1995) 1895. [13] E. Karp~ti-Smidr6czki, A. B/iri and M. Zrinyi, Colloid. Polm. Sci. 273 (1995) 857. [14] A.P. Prudnikov, Yu.A. Brychkov and O.1. Marichev, Inteyrals and Series, Vol. 2 (Gordon and Breach, London, 1988) pp. 311.