Dynamics of annihilation. II. Fluctuations of global quantities

Dynamics of annihilation. II. Fluctuations of global quantities Pablo Maynar,1,2 María Isabel García de Soria,3Grégory Schehr,1Alain Barrat,1and Emmanuel Trizac3 1Laboratoire de Physique Théorique (CNRS UMR 8627), Bâtiment 210, Université Paris-Sud, 91405 Orsay Cedex, France 2Física Teórica, Universidad de Sevilla, Apartado de Correos 1065, E-41080, Sevilla, Spain 3LPTMS (CNRS UMR 8626), Université Paris-Sud, Orsay Cedex, F-91405, France 共Received 17 January 2008; published 28 May 2008兲 We develop a theory for fluctuations and correlations in a gas evolving under ballistic annihilation dynamics. Starting from the hierarchy of equations governing the evolution of microscopic densities in phase space, we subsequently restrict our attention to a regime of spatial homogeneity, and obtain explicit predictions for the fluctuations and time correlation of the total number of particles, total linear momentum, and total kinetic energy. Cross correlations between these quantities are worked out as well. These predictions are successfully tested against molecular dynamics and Monte Carlo simulations. This provides strong support for the theoretical approach developed, including the hydrodynamic treatment of the spectrum of the linearized Boltzmann operator. This paper makes use of the spectral analysis reported in the preceding paper 关Phys. Rev. E 77, 051127 共2008兲兴. DOI: 10.1103/PhysRevE.77.051128 PACS number共s兲: 02.50.⫺r, 51.10.⫹y, 05.20.Dd, 82.20.Nk I. INTRODUCTION Systems where particles may react, change chemical or physical nature, and ultimately disappear, model a rich variety of phenomena and provide prominent situations to develop and test the foundations of nonequilibrium statistical mechanics 共see, e.g., 关2–6兴and references therein兲. When reactions are controlled ballistically, the system can be modeled by an assembly of hard spheres or disks which annihilate with probability por collide elastically with probability 1−pevery time two particles meet each other 关7兴. Within the framework of this probabilistic ballistic annihilation 共PBA兲 model, most of the work carried out up to now has focused on the kinetic equations for the one-body distribution function and the information following from them 关7–11兴. In particular, the hydrodynamic equations, with explicit expressions for the transport coefficients, have been derived by using a generalization of the Chapman-Enskog expansion 关12兴. The preceding paper 关1兴, where we have established the hydrodynamic description in the context of the eigenfunctions and eigenvalues of the linearized Boltzmann collision operator, falls in this vein. Within this formalism, the conditions in which the hydrodynamic description is expected to apply are somewhat more transparent, and can be expressed in terms of some properties of the linearized Boltzmann collision operator. In the present paper, the goal is to go beyond the study of one-body quantities; the focus is on fluctuations and correlations. To this end, we use the tools and ideas developed in the context of the linearized Boltzmann equation. The dynamical behavior of the correlations in the dilute limit can indeed be obtained in terms of the linearized Boltzmann collision operator. The study of correlations in the PBA model allows one to go beyond the description at the level of average values, and to characterize how global magnitudes 共such as the total number of particles or the total energy兲fluctuate around their average. It has already been shown in other classes of dissipative systems, such as in granular systems, that the knowledge of fluctuations is relevant in order to understand the behavior of the system when vortices or clusters develop 关13,14兴, or even in simpler situations where the system is still homogeneous 关15,16兴. The main goal of this paper is to formulate a theory of fluctuations for the PBA model in the dilute limit and to apply it to one of the simplest possible states, namely, the homogeneous decay state, exploiting its scaling properties. This will allow us to obtain explicit expressions for the distributions characterizing the velocity correlations in the system, and to compute the statistics of the total number of particles, total momentum, and total energy, which decrease monotonically due to the annihilation process. The paper is organized as follows. In Sec. II, we present the general framework of the hierarchy method 关17兴which allows us to write the evolution equations of correlation functions. The specific case of the homogeneous decay state is considered in Sec. III, where the scaling properties of this state are used to simplify the equations. After briefly recalling in Sec. IV how fluctuations and correlations of global observables can be computed from the knowledge of the two-particle correlation functions, we first consider in Sec. V the correlation functions at equal time, which give access to the fluctuations of the total number of particles, total momentum, and total energy. We obtain theoretical predictions for the asymptotic scaling regime as well as for the shorttime behavior, and we test these predictions against numerical simulations 共both molecular dynamics and Monte Carlo兲. In Sec. VI, we generalize our results to the two-time correlation functions and compare also to numerical simulations. Finally, Sec. VII contains some discussions of the results and our conclusions. For the sake of readability, this paper contains some overlap with its companion 关1兴. Repetitions have been kept to a minimum though, and we therefore refer to 关1兴 for several technical details. II. GENERAL FRAMEWORK The system we consider consists of a dilute gas of identical smooth hard spheres or disks of mass mand diameter ␴ , PHYSICAL REVIEW E 77, 051128 共2008兲 1539-3755/2008/77共5兲/051128共14兲©2008 The American Physical Society051128-1 moving ballistically in dimension d. The particles are supposed to undergo only binary collisions. When two particles collide they annihilate with probability por collide elastically with probability 1−p. In this probabilistic ballistic annihilation, there is therefore no conserved quantity 共except for p=0兲and the number of particles decreases steadily. In this section, we will show how to obtain evolution equations for the correlation functions in this system. The general idea of the method is to derive a closed set of equations for the distribution functions describing the fluctuations by using the same kind of approximations as needed to derive the kinetic equation, in our case the Boltzmann equation. In this way, a unified formalism provides the usual kinetic equation as well as evolution equations for the one- and two-time correlations. Let Xj⬅兵Rj共t兲,Vj共t兲其 denote the position and velocity of particle jin the system at time t. Both, Rj共t兲and Vj共t兲are parametric functions of the initial positions and velocities of all particles. Microscopic one- and two-particle densities in the phase space are defined by F1共x1,t兲=兺 i=1 N ␦ 关x1−Xi共t兲兴,共1兲 F2共x1,x2,t兲=兺 i=1 N 兺 i⫽j N ␦ 关x1−Xi共t兲兴 ␦ 关x2−Xj共t兲兴,共2兲 and higher-order functions can similarly be defined. Here and in the following, the lower-case variables xi⬅兵ri,vi其are field variables referring to a particular point in phase space. The initial state of the system is characterized by a point in phase space, ⌫⬅兵X1,...,XN其, which is chosen at random according to a probability ␳ 共⌫,0兲. When the notation 具G典 ⬅兰d⌫G共⌫兲 ␳ 共⌫,0兲is introduced for the average over the initial conditions, the averages of the microscopic densities Fs共x1,...,xs,t兲over ␳ 共⌫,0兲are the usual one-time reduced distribution functions fs共x1, ...,xs,t兲=具Fs共x1, ...,xs,t兲典.共3兲 Similarly, two-time reduced distribution functions can also be defined in terms of the microscopic densities as fr,s共x1, ...,xr,t;x1 ⬘, ...,xs ⬘,t⬘兲 =具Fr共x1, ...,xr,t兲Fs共x1 ⬘, ...,xs ⬘,t⬘兲典.共4兲 For simplicity we will consider t⬎t⬘⬎0 in the following. We now introduce the two-particle correlation functions through the usual cluster expansion. The one-time correlation function g2and the two-time correlation function h1,1 are then defined by f2共x1,x2,t兲=f1共x1,t兲f1共x2,t兲+g2共x1,x2,t兲,共5兲 f1,1共x1,t;x2,t⬘兲=f1共x1,t兲f1共x2,t⬘兲+h1,1共x1,t;x2,t⬘兲.共6兲 It is easy to show from the definitions of f1,f2, and f1,1 that h1,1共x1,t;x2,t兲=g2共x1,x2,t兲+ ␦ 共x1−x2兲f1共x1,t兲.共7兲 The case of deterministic annihilation 共p=1兲was considered in Ref. 关10兴. The hierarchy of equations for the reduced distribution functions is then shown to be similar to the one describing elastic collisions, once the binary elastic collision operator is replaced by the operator describing annihilating collisions. In the PBA case, assuming molecular chaos, i.e., that no correlations exist between colliding particles, the equation for f1共x1,t兲is the Boltzmann equation 冉 ⳵ ⳵ t+L共0兲共x1兲 冊 f1共x1,t兲=J关f1,f1兴,共8兲 where L共0兲共x1兲=v1· ⳵ ⳵ r1 ,共9兲 J关f1,f1兴= 冕 dx2 ␦ 共r12兲T ¯ 0共v1,v2兲f1共x1,t兲f1共x2,t兲,共10兲 and T ¯ 0共v1,v2兲= ␴ d−1 冕 d ␴ ˆ⌰共v12 · ␴ ˆ兲共v12 · ␴ ˆ兲 ⫻关共1−p兲共b ␴ −1 −1兲−p兴,共11兲 is the PBA binary collision operator. In the above expressions v12=v1−v2is the relative velocity, r12=r1−r2the relative position, ⌰the Heaviside function, ␴ ˆa unit vector joining the centers of two particles at collision, and b ␴ −1 an operator that replaces all the velocities v1and v2appearing to its right by the precollisional values v1 ⴱand v2 ⴱ, b ␴ −1v1⬅v1 ⴱ=v1−共v12 · ␴ ˆ兲 ␴ ˆ,共12兲 b ␴ −1v2⬅v2 ⴱ=v2+共v12 · ␴ ˆ兲 ␴ ˆ.共13兲 The equation for the correlation function g2can be obtained under the same hypothesis required to derive the Boltzmann equation, following the same lines as in Ref. 关18兴in the case of inelastically colliding particles, as 冉 ⳵ ⳵ t+L共0兲共x1兲+L共0兲共x2兲−K共x1,t兲−K共x2,t兲 冊 g2共x1,x2,t兲 = ␦ 共r12兲T ¯ 0共v1,v2兲f1共x1,t兲f1共x2,t兲,共14兲 where we have introduced the linear operator K共xi,t兲, K共xi,t兲= 冕 dx3 ␦ 共ri3兲T ¯ 0共vi,v3兲共1+Pi3兲f1共x3,t兲,共15兲 and where the permutation operator Pab interchanges the labels of particles aand bin the quantities on which it acts. Finally, the evolution equation for h1,1 reads 冉 ⳵ ⳵ t+L共0兲共x1兲−K共x1,t兲 冊 h1,1共x1,t;x2,t⬘兲=0, t⬎t⬘, 共16兲 which has to be solved with the initial condition 共7兲, h1,1共x1,t⬘;x2,t⬘兲=g2共x1,x2,t⬘兲+ ␦ 共x1−x2兲f1共x1,t⬘兲. The equations for the correlation functions h1,1 and g2 contain a part corresponding to free streaming and another MAYNAR et al. PHYSICAL REVIEW E 77, 051128 共2008兲 051128-2 one that corresponds to collisions. In particular, the evolution of the one-time correlation function g2due to collisions has two parts: one due to collisions of particle 1 or 2 共corresponding to the indices of the correlation function兲with a third particle, which is governed by the Boltzmann collision operator; and a second one, due to collisions of particle 1 with particle 2, which can be written in terms of the oneparticle distribution function as a consequence of the molecular chaos hypothesis. In fact, as in the case of the inelastic granular gas 关18兴, the only difference between the evolution equations of the correlation functions for the PBA and for a system of elastic particles lies in the substitution of the elastic binary collision operator by the operator for the PBA model, T ¯ 0共v1,v2兲. However, this does not give any a priori guarantee on the range of validity of these equations. This prescription, i.e., how small the density of the system must be so that the above kinetic equations provide an accurate description, might depend on the parameter p, and also on the particular state being considered. We will see in Sec. IV how the knowledge of h1,1 and g2 allows us to obtain the correlation functions of any observable that is a function of the particle positions in phase space. III. HOMOGENEOUS DECAY STATE As recalled in the preceding paper 关1兴, the Boltzmann equation for the PBA model 共8兲admits a particular solution fH共v,t兲describing a spatially homogeneous decay state 共HDS兲, in which all the time dependence is contained in the evolution of the density nH共t兲and the temperature TH共t兲, which are defined as in standard kinetic theory as nH共t兲= 冕 dvfH共v,t兲,d 2nH共t兲TH共t兲= 冕 dvm 2v2fH共v,t兲. 共17兲 Although there exists no rigorous proof of its stability nor of the fact that such a state should be approached from arbitrary initial conditions, numerical results obtained by molecular dynamic simulations and by the direct Monte Carlo method strongly support the existence of such a homogeneous solution 关7,9兴. In this section, we review for completeness the evolution equation of the one-particle distribution function and obtain the equations for adequately rescaled correlation functions. All quantities concerning this homogeneous decay state will be labeled by an index H. In the HDS, the one-body distribution function does not depend on space and follows the scaling form 关10兴 fH共v,t兲=nH共t兲 vH d共t兲 ␹ H共c兲,共18兲 where nH共t兲is the uniform density, vH共t兲is the thermal 共rootmean-square兲velocity related to the granular temperature TH共t兲by vH共t兲= 冉 2TH共t兲 m 冊 1/2 ,共19兲 and ␹ 共c兲is an isotropic function depending only on the modulus c=兩c兩of the rescaled velocity c=v/vH共t兲. Moreover, the density and temperature fields evolve according to 关12兴 dnH共t兲 dt =−p ␯ H共t兲 ␨ nnH共t兲,共20兲 dTH共t兲 dt =−p ␯ H共t兲 ␨ TTH共t兲,共21兲 where ␯ H共t兲is the collision frequency of the corresponding hard-sphere fluid in equilibrium 共with same temperature and density兲 ␯ H共t兲=nH共t兲TH 1/2共t兲 ␴ d−1 m1/2 8 ␲ d−1/2 共d+2兲⌫共d/2兲,共22兲 and the dimensionless decay rates ␨ nand ␨ Tare functionals of the distribution function: p ␨ n=− ␥ 2 冕 dc1 冕 dc2T共c1,c2兲 ␹ H共c1兲 ␹ H共c2兲,共23兲 p ␨ T=− ␥ 2 冕 dc1 冕 dc2 冉 2c1 2 d−1 冊 T共c1,c2兲 ␹ H共c1兲 ␹ H共c2兲. 共24兲 In these expressions, ␥ , which does not depend on time, is given by ␥ =2nH共t兲vH共t兲 ␴ d−1/ ␯ H共t兲=共d+2兲冑2⌫共d/2兲/ 共4 ␲ 共d−1兲/2兲, and the binary collision operator T共c1,c2兲takes the form T共c1,c2兲= 冕 d ␴ ˆ⌰共c12 · ␴ ˆ兲共c12 · ␴ ˆ兲关共1−p兲b ␴ −1 −1兴. 共25兲 Finally, the scaled distribution function ␹ H共c兲obeys the equation p 冋 共d ␨ T−2 ␨ n兲+ ␨ Tc1· ⳵ ⳵ c1 册 ␹ H共c1兲 = ␥ 冕 dc2T共c1,c2兲 ␹ H共c1兲 ␹ H共c2兲.共26兲 The operator b ␴ −1 in the last equation is defined again by Eqs. 共12兲and 共13兲, but substituting 共v1,v2兲by 共c1,c2兲. The analytical form of ␹ His not known, but its behavior at large and small velocities has been determined 关9,10兴. As in the preceding paper 关1兴, we will use here the approximate form of the distribution function in the so-called first Sonine approximation, which is valid for velocities in the thermal region, and all the functionals of ␹ H共c兲, like the decay rates and the transport coefficients, will be evaluated in this approximation 关7,9兴. Considering the scaling form for the one-particle distribution function, it is convenient to introduce the rescaled correlation function g ˜ Hthrough DYNAMICS OF ANNIHILATION. II. FLUCTUATIONS OF…PHYSICAL REVIEW E 77, 051128 共2008兲 051128-3 g2,H共x1,x2,t兲=nH共t兲 vH 2d共t兲g ˜ H共 ␶ ,r12,c1,c2兲,共27兲 where we have taken into account that the system is invariant under space translation, so that g2,Hdepends on r12=r1−r2 and not on r1and r2separately. The dimensionless time scale ␶ , ␶ =1 2 冕 0 t dt⬘ ␯ H共t⬘兲,共28兲 is proportional to the number of collisions in the time interval 关0,t兴. The equation for the reduced function g ˜ Hin these units reads then 冉 − ⳵ ⳵ ␶ +⌳共c1兲+⌳共c2兲−2p ␨ n−lH共 ␶ 兲c12 · ⳵ ⳵ r12 冊 ⫻g ˜ H共 ␶ ,r12,c1,c2兲=− ␦ 共r12兲 ␥ T共c1,c2兲 ␹ H共c1兲 ␹ H共c2兲, 共29兲 where we have also introduced the length scale lH共t兲 =2vH共t兲/ ␯ H共t兲, which is proportional to the instantaneous mean free path, and the linearized Boltzmann operator 共see previous paper兲 ⌳共ci兲h共ci兲= ␥ 冕 dc3T共ci,c3兲共1+Pi3兲 ␹ H共c3兲h共ci兲 +p共2 ␨ n−d ␨ T兲h共ci兲−p ␨ Tci· ⳵ ⳵ ci h共ci兲.共30兲 Similarly, we define a rescaled two-time correlation function h ˜ Hthrough h1,1,H共x1,t;x2,t⬘兲=nH共t兲 vH d共t兲vH d共t⬘兲h ˜ H共r12;c1, ␶ ;c2, ␶ ⬘兲,共31兲 and obtain the following evolution equation: ⳵ ⳵ ␶ h ˜ H共r12;c1, ␶ ;c2, ␶ ⬘兲 = 冉 ⌳共c1兲−lH共 ␶ 兲c1· ⳵ ⳵ r1 冊 h ˜ H共r12;c1, ␶ ;c2, ␶ ⬘兲,共32兲 with h ˜ H共r12;c1, ␶ ;c2, ␶ 兲=g ˜ H共 ␶ ,r12,c1,c2兲+ ␦ 共c1−c2兲 ␦ 共r12兲 ␹ H共c1兲. It is interesting to note how, in this representation, all the time dependence due to the reference state is absorbed in the free streaming term through the function lH共 ␶ 兲, proportional to the mean free path. The evolution of the correlation functions, moreover, will be determined by the properties of the linearized Boltzmann operator ⌳, which we have already studied in the preceding paper 关1兴and which we will recall in Sec. V. IV. FROM PARTICLE CORRELATION FUNCTIONS TO CORRELATIONS AND FLUCTUATIONS OF GLOBAL MAGNITUDES In this section, we will show how the previously presented framework for correlation functions will allow us to study the fluctuations and correlations of global quantities for a PBA system in the homogeneous decay state. In particular, we will focus on the total number of particles N, the total momentum P, and the total energy E. Consider indeed two dynamical variables of the form A关⌫共t兲兴 =兺 i=1 N a共Vi兲= 冕 dx1a共v1兲F1共x1,t兲, B关⌫共t兲兴 =兺 i=1 N b共Vi兲= 冕 dx2b共v2兲F1共x2,t兲,共33兲 where aand bare functions of the particle velocities Vi, and F1共x1,t兲is the microscopic density in phase space 共1兲. Taking a=1, a=v, and a=mv2/2 yields for Athe total number of particles N, the total momentum P, and the total kinetic energy E, respectively. The deviations ␦ A共t兲=A共t兲−具A共t兲典H and ␦ B共t兲=B共t兲−具B共t兲典Hof Aor Bfrom their average values in the HDS 共denoted by 具¯典H兲, define the fluctuations of both magnitudes, which have average zero and correlations 具 ␦ A共t兲 ␦ B共t⬘兲典H=具A共t兲B共t⬘兲典H−具A共t兲典H具B共t⬘兲典H.共34兲 It is then straightforward to use the definition of the two-time correlation function h1,1 in Eq. 共6兲, to obtain 具 ␦ A共t兲 ␦ B共t⬘兲典H= 冕 dr1 冕 dv1 冕 dr2 冕 dv2a共v1兲b共v2兲 ⫻h1,1,H共r1,v1,t;r2,v2,t⬘兲.共35兲 In particular, for t=t⬘, this leads to 具 ␦ A共t兲 ␦ B共t兲典H=V 冕 dva共v兲b共v兲fH共v,t兲+V 冕 dv1 冕 dv2 ⫻a共v1兲b共v2兲 冕 dr12g2,H共r12,v1,v2,t兲,共36兲 where V=兰dr1is the total volume of the system. These formulas show how the correlations of two different global magnitudes are determined by the one-particle distribution function and by the correlation functions. The one-particle distribution function is known in the HDS in the first Sonine approximation 关7,9兴. V. FLUCTUATIONS IN THE HDS Let us focus in this section on the one-time correlation function g ˜ H. We will only need functions aand bwhich depend on the velocity degrees of freedom, so that it is convenient to integrate out the spatial dependence by introducing ␾ H共 ␶ ,c1,c2兲⬅ 冕 dr12g ˜ H共 ␶ ,r12,c1,c2兲,共37兲 whose evolution is obtained from 共29兲as MAYNAR et al. PHYSICAL REVIEW E 77, 051128 共2008兲 051128-4 冉 − ⳵ ⳵ ␶ +⌳共c1兲+⌳共c2兲−2p ␨ n 冊 ␾ H共 ␶ ,c1,c2兲 =− ␥ T共c1,c2兲 ␹ H共c1兲 ␹ H共c2兲.共38兲 Given an initial condition ␾ H共0,c1,c2兲, this Eq. 共38兲can be formally integrated as ␾ H共 ␶ ,c1,c2兲 =e关⌳共c1兲+⌳共c2兲−2p ␨ n兴 ␶ ␾ H共0,c1,c2兲+ 冕 0 ␶ d ␶ ⬘ ⫻e共 ␶ − ␶ ⬘兲关⌳共c1兲+⌳共c2兲−2p ␨ n兴 ␥ T共c1,c2兲 ␹ H共c1兲 ␹ H共c2兲 =e关⌳共c1兲+⌳共c2兲−2p ␨ n兴 ␶ 关 ␾ H共0,c1,c2兲− ␾ H s共c1,c2兲兴+ ␾ H s共c1,c2兲, 共39兲 where ␾ H s共c1,c2兲is the solution of 关⌳共c1兲+⌳共c2兲−2p ␨ n兴 ␾ H s共c1,c2兲=− ␥ T共c1,c2兲 ␹ H共c1兲 ␹ H共c2兲, 共40兲 where we implicitly assumed that ⌳is invertible. This happens to be the case; see below. The spectral properties of the linearized Boltzmann operator ⌳are thus crucial for the evaluation of ␾ H. We therefore start by recalling these properties. A. Spectral properties of ⌳ In our companion paper 关1兴, we have analyzed the eigenvalue problem associated with the linearized Boltzmann operator ⌳, ⌳共c兲 ␰ ␤ 共c兲=␭ ␤ ␰ ␤ 共c兲.共41兲 We have in fact restricted ourselves to the hydrodynamic part of ⌳, defined by those eigenvalues of the balance equations for the number density, momentum, and temperature following from the homogeneous linearized Boltzmann equation. Such eigenvalues are 关12兴 ␭1=0, ␭2=−p共 ␨ T+2 ␨ n兲,␭3=p ␨ T.共42兲 Although we were not able to prove that these eigenvalues are indeed the hydrodynamic ones, the self-consistency of the resulting description and the successful comparison with numerical simulations have validated this assumption 关1兴.In the previous paper we also obtained the corresponding eigenfunctions ␰ 1共c兲= ␹ H共c1兲+ ⳵ ⳵ c·关c ␹ H共c兲兴,共43兲 ␰ 2共c兲=z ␹ H共c兲− ⳵ ⳵ c·关c ␹ H共c兲兴,共44兲 ␰ 3共c兲=− ⳵␹ H共c兲 ⳵ c,共45兲 with z=2 ␨ n ␨ Ta function of the probability of annihilation, p. The eigenvalue ␭3is d-fold degenerate and we denote by ␰ 3i, i=1,...,d, the corresponding eigenvectors. The scalar product of two functions f共c兲and g共c兲is defined as 具f兩g典⬅ 冕 dc ␹ H −1共c兲fⴱ共c兲g共c兲,共46兲 fⴱbeing the complex conjugate of f. The eigenfunctions ␰ ␤ given in 共43兲–共45兲are not orthogonal, as a consequence of the operator ⌳being non-Hermitian in the associated Hilbert space. On the other hand, it is easily verified that the set of functions 兵 ␰ ¯ 1; ␰ ¯ 2; ␰ 3其=兵 ␹ H共c兲−z 1+z共1 2+c2 d兲 ␹ H共c兲; 1 1+z共1 2+c2 d兲 ␹ H共c兲;c ␹ H共c兲其 satisfy the biorthogonality condition 具 ␰ ¯ ␤ 兩 ␰ ␤ ⬘典= ␦ ␤ , ␤ ⬘, for ␤ , ␤ ⬘=1,2,3. The eigenfunctions of the operator 关⌳共c1兲+⌳共c2兲−2p ␨ n兴 that appear in Eq. 共40兲are then easily seen to be ␰ ␤ 1共c1兲 ␰ ␤ 2共c2兲, with 关⌳共c1兲+⌳共c2兲−2p ␨ n兴 ␰ ␤ 1共c1兲 ␰ ␤ 2共c2兲 =共␭ ␤ 1+␭ ␤ 2−2p ␨ n兲 ␰ ␤ 1共c1兲 ␰ ␤ 2共c2兲.共47兲 Since ␨ n⬎ ␨ T关12兴, and under the assumption that the norms of the “nonhydrodynamic” eigenvalues are always greater than the hydrodynamic ones, the eigenvalues of 关⌳共c1兲 +⌳共c2兲−2p ␨ n兴are therefore all negative. This has the important consequence that the exponential term in 共39兲decays to zero and that the large-time limit of ␾ H共 ␶ ,c1,c2兲is ␾ H s共c1,c2兲, the solution of Eq. 共40兲. B. Hydrodynamic part of the correlation functions Obtaining the full spectrum of ⌳is a formidable task. We will here assume, as in the preceding paper 关1兴, that the kinetic 共nonhydrodynamic兲modes have a fast decay, and work in the hydrodynamic subspace spanned by the functions ␰ ␤ defined in the previous section. With that purpose, we generalize the definition of the scalar product given in 共46兲to two-velocity functions by 具f共c1,c2兲兩g共c1,c2兲典 ⬅ 冕 dc1 冕 dc2 ␹ H −1共c1兲 ␹ H −1共c2兲fⴱ共c1,c2兲g共c1,c2兲. 共48兲 This allows us to define a projector operator P12 onto the space spanned by the functions ␰ ␤ 1共c1兲 ␰ ␤ 2共c2兲as P12f共c1,c2兲 ⬅兺 ␤ 1=1 3 兺 ␤ 2=1 3 具 ␰ ¯ ␤ 1共c1兲 ␰ ¯ ␤ 2共c2兲兩f共c1,c2兲典 ␰ ␤ 1共c1兲 ␰ ␤ 2共c2兲, 共49兲 and the hydrodynamic parts of ␾ H共 ␶ ,c1,c2兲and ␾ H s共c1,c2兲 are by definition ␾ H 共h兲共 ␶ ,c1,c2兲⬅P12 ␾ H共 ␶ ,c1,c2兲 =兺 ␤ 1, ␤ 2=1 3 a ␤ 1, ␤ 2共 ␶ 兲 ␰ ␤ 1共c1兲 ␰ ␤ 2共c2兲,共50兲 DYNAMICS OF ANNIHILATION. II. FLUCTUATIONS OF…PHYSICAL REVIEW E 77, 051128 共2008兲 051128-5 ␾ H s共h兲共c1,c2兲⬅P12 ␾ H s共c1,c2兲=兺 ␤ 1, ␤ 2=1 3 a ␤ 1, ␤ 2 s ␰ ␤ 1共c1兲 ␰ ␤ 2共c2兲. 共51兲 We can now obtain a closed equation for ␾ H 共h兲by applying the operator P12 on both sides of Eqs. 共39兲and 共40兲, under the additional assumption that P12⌳共ci兲=P12⌳共ci兲P12.共52兲 A theoretical estimation a priori of the accuracy of this approximation would require more knowledge than is available at present about the spectrum of ⌳共c兲and its adjoint. Therefore, it will be taken as a working hypothesis, to be validated later on by comparing the predictions it leads to with the results from numerical simulations of the system. It is worth emphasizing that, since ⌳leaves the hydrodynamic subspace invariant, Eq. 共52兲is equivalent to the commutation relation P12⌳=⌳P12. Proceeding further, we obtain from 共39兲 ␾ H 共h兲共 ␶ ,c1,c2兲 =e关⌳共c1兲+⌳共c2兲−2p ␨ n兴 ␶ P12关 ␾ H共0,c1,c2兲− ␾ H s共c1,c2兲兴 +P12 ␾ H s共c1,c2兲 =兺 ␤ 1, ␤ 2=1 3 共A ␤ 1, ␤ 2e共␭ ␤ 1+␭ ␤ 2−2p ␨ n兲 ␶ +a ␤ 1, ␤ 2 s兲 ␰ ␤ 1共c1兲 ␰ ␤ 2共c2兲, 共53兲 where we have introduced A ␤ 1, ␤ 2=a ␤ 1, ␤ 2共0兲−a ␤ 1, ␤ 2 s. We show in Appendix A how to obtain explicit formulas for the coefficients a ␤ 1, ␤ 2 sin terms of the cooling rates ␨ nand ␨ Tand other coefficients, which are also functionals of the one-time distribution function ␹ H. The values of a ␤ 1, ␤ 2共0兲depend on the initial condition ␾ H共0,c1,c2兲. For the specific case in which the variables N,P, and Edo not fluctuate at t=0, and taking into account that the system is in the HDS, the coefficients a ␤ 1, ␤ 2共0兲are calculated in Appendix B. C. Hydrodynamic approximation for global fluctuations In this section we compute the correlation functions of the global observables by replacing ␾ Himplicitly appearing in 共36兲by its hydrodynamic part ␾ H 共h兲. This can be done by invoking the relation 具 ␰ ¯ ␤ 1共c1兲 ␰ ¯ ␤ 2共c2兲兩f共c1,c2兲典 =具 ␰ ¯ ␤ 1共c1兲 ␰ ¯ ␤ 2共c2兲兩f共h兲共c1,c2兲兩典,共54兲 for ␤ i=1,2,3. However, it must be stressed that the theoretical prediction for ␾ H 共h兲in Eq. 共53兲has been calculated using the approximation 共52兲. If we substitute a共v兲=1 and b共v兲=1 in 共36兲, we obtain for the fluctuations of the number of particles 具 ␦ N2共 ␶ 兲典H=NH共 ␶ 兲 冋 冕 dc ␹ H共c兲+ 冕 dc1 冕 dc2 ␾ H 共h兲共 ␶ ,c1,c2兲 册 , 共55兲 where we have introduced the notation NH⬅具N典H. In order to calculate the fluctuations of the total momentum we substitute a共v兲=viand b共v兲=vjin Eq. 共36兲and obtain 具 ␦ Pi共 ␶ 兲 ␦ Pj共 ␶ 兲典H=NH共 ␶ 兲vH 2共 ␶ 兲 冋 冕 dccicj ␹ H共c兲 + 冕 dc1 冕 dc2c1ic2j ␾ H 共h兲共 ␶ ,c1,c2兲 册 . 共56兲 For the energy we substitute a共v兲=b共v兲=1 2mv2so that we have 具 ␦ E2共 ␶ 兲典H=m2 4NH共 ␶ 兲vH 4共 ␶ 兲 冋 冕 dcc4 ␹ H共c兲 + 冕 dc1 冕 dc2c1 2c2 2 ␾ H 共h兲共 ␶ ,c1,c2兲 册 .共57兲 Finally, we can calculate the correlation between ␦ Nand ␦ E by taking a共v兲=1 and b共v兲=1 2mv2 具 ␦ N共 ␶ 兲 ␦ E共 ␶ 兲典H=m 2NH共 ␶ 兲vH 2共 ␶ 兲 冋 冕 dcc2 ␹ H共c兲 + 冕 dc1 冕 dc2c2 2 ␾ H 共h兲共 ␶ ,c1,c2兲 册 .共58兲 After some algebra, we obtain 具 ␦ N2共 ␶ 兲典H=NH共 ␶ 兲共1+a1,1 s+2za1,2 s+z2a2,2 s+A1,1e−2p ␨ n ␶ +2zA1,2e−p共 ␨ T+4 ␨ n兲 ␶ +z2A2,2e−2p共 ␨ T+3 ␨ n兲 ␶ 兲,共59兲 具 ␦ Pi共 ␶ 兲 ␦ Pj共 ␶ 兲典H= ␦ ijNH共 ␶ 兲vH 2共 ␶ 兲 冉 1 2+a3i,3i s 冊 共1−e−2p共 ␨ n− ␨ T兲 ␶ 兲, 共60兲 具 ␦ E2共 ␶ 兲典H=m2 4NH共 ␶ 兲vH 4共 ␶ 兲 冋 d共d+2兲 4共1+a2兲+d2 4a1,1 s −d2 冉 1+ z 2 冊 a1,2 s+d2 冉 1+ z 2 冊 2 +d2 4A1,1e−2p ␨ n ␶ −d2 冉 1+ z 2 冊 A1,2e−p共 ␨ T+4 ␨ n兲 ␶ +d2 冉 1+ z 2 冊 2 A2,2e−2p共 ␨ T+3 ␨ n兲 ␶ 册 ,共61兲 and 具 ␦ N共 ␶ 兲 ␦ E共 ␶ 兲典H=m 2NH共 ␶ 兲vH 2共 ␶ 兲 冋 d 2−d 2a1,1 s+da1,2 s +dz 冉 1+ z 2 冊 a2,2 s−d 2A1,1e−2p ␨ n ␶ +dA1,2e−p共 ␨ T+4 ␨ n兲 ␶ +dz 冉 1+ z 2 冊 A2,2e−2p共 ␨ T+3 ␨ n兲 ␶ 册 ,共62兲 MAYNAR et al. PHYSICAL REVIEW E 77, 051128 共2008兲 051128-6 where a2is related to the fourth moment of ␹ H共c兲through 兰dcc4 ␹ H共c兲=d共d+2兲 4共1+a2兲, and has been evaluated in the first Sonine approximation in 关10兴. All the functions a ␣ , ␤ sand A ␣ , ␤ are evaluated in the Appendixes A and B. At this point, it is important to note that the Eqs. 共59兲–共62兲have been obtained under the assumption that the system is in the homogeneous decay state at all times, i.e., that the one-particle distribution function is ␹ H共c兲for the whole time evolution. For p⬍1, if we start with an arbitrary initial condition, numerical simulations show that, after a few collisions, the distribution function reaches the scaling regime given by Eq. 共18兲. Then, the evolution of ␾ H 共h兲is given by 共53兲and one expects that the same correlation functions 共59兲–共62兲will be obtained in the long-time limit, independently of the initial condition. This will be confirmed in the next section by numerical simulations. Equations 共59兲–共62兲lead to a certain number of theoretical predictions. In particular, they imply that the ratios 具 ␦ N2共 ␶ 兲典/NH共 ␶ 兲,具 ␦ E2共 ␶ 兲典/关NH共 ␶ 兲vH 4共 ␶ 兲兴,具 ␦ Pi 2共 ␶ 兲典/ 关NH共 ␶ 兲vH 2共 ␶ 兲兴, and 具 ␦ N共 ␶ 兲 ␦ E共 ␶ 兲典/关NH共 ␶ 兲vH 2共 ␶ 兲兴 reach stationary values at large times. The approach to these stationary values is exponential in ␶ , and is slower for the correlations of the total momentum, since the argument of the exponential is p共 ␨ n− ␨ T兲 ␶ , while the other quantities evolve on faster time scales. D. Numerical simulations We now compare our theoretical predictions with the results of molecular dynamics 共MD兲and direct simulation Monte Carlo 共DSMC兲calculations of a freely evolving system of Nhard disks of diameter ␴ which annihilate with probability por collide elastically with probability 1−pevery time two particles meet each other. In the MD case, the particles were localized in a square box of size Lwith periodic boundary conditions. The event-driven algorithm 关19兴 has been used and the initial density has been chosen low enough to be always in the dilute limit. The parameters for all the MD simulations were N共0兲=105,nH共0兲 ␴ 2=0.05, TH共0兲=1, and 0⬍pⱕ1. In the case of the DSMC simulations we have used Bird’s algorithm 关20兴with the same values of the parameters, except that the density that plays no role. The initial velocity distribution is a Maxwellian in both cases. We have measured the time evolution of the total number of particles and the total energy, averaging the data over various initial conditions 共the total momentum fluctuates around zero兲. We first checked that the Eqs. 共20兲and 共21兲, with the theoretical predictions derived in 关10兴for the cooling rates, correctly describe the decay of the average global quantities. In the same way, we have obtained the averaged values of N2共t兲,E2共t兲,Pi 2共t兲, and N共t兲E共t兲共the correlations between Piand Nor Eare zero兲. Figures 1–3show the time evolution of the various onetime correlation functions considered, for p=0.5 and 0.8. The DSMC results have been averaged over 4000 trajectories while the MD simulations have been averaged over 150 trajectories 共the DSMC method being computationally less expensive, it is then possible to average over a larger number of initial conditions than for the MD simulations兲. The dashed lines are the theoretical predictions, Eqs. 共59兲–共62兲. Note, however, that the system is not initially in the HDS; the initial distribution function is a Maxwellian and not ␹ H. Nevertheless, as the difference between these two distributions is very small 共at least for thermal velocities since a2 ⬃0.1兲and as the stationary values depend very weakly on p, Eqs. 共59兲–共62兲predict quite well the time evolution measured in the simulations. The ratios 具 ␦ N2共 ␶ 兲典/NH共 ␶ 兲, 具 ␦ E2共 ␶ 兲典/关NH共 ␶ 兲vH 4共 ␶ 兲兴, and 具 ␦ N共 ␶ 兲 ␦ E共 ␶ 兲典/关NH共 ␶ 兲vH 2共 ␶ 兲兴 reach stationary values as predicted. The fluctuations of the total momentum evolve more slowly, as also predicted, and the stationary value of 具 ␦ Pi 2共 ␶ 兲典/关NH共 ␶ 兲vH 2共 ␶ 兲兴 is barely reached. Note that ␶ =4 corresponds for p=0.5 to a total number of particles at the end of the simulation N⯝1700. We have performed simulations starting with other initial conditions further from the HDS. The initial velocity distribution function has been chosen as a constant function in a square centered in the origin in the velocity space such that the initial temperature is unity. As seen in Figs. 4and 5,we obtain a different short-time evolution but the scaled moments still converge toward the HDS values, which are represented by the dashed lines. The convergence is slower as we increase the value of p, and 具 ␦ E2共 ␶ 兲典/关NH共 ␶ 兲vH 4共 ␶ 兲兴, whose magnitude depends on the higher moments of the velocity distribution, is the most affected. Figures 6and 7show the comparison between the stationary values of the various ratios measured in the simulations and the theoretical predictions in Eq. 共59兲–共62兲at large ␶ . The agreement is very good for all values investigated. We have also computed the probability distribution for the number of particles, energy, and momentum. As we can see 00.5 11.5 22. 5 τ 0 0.5 1 1.5 2 2.5 3 <(δΝ)2>/<Ν> <(δΕ)2>/<Ν><Τ>2 <δΕδΝ>/<Ν><Τ> 00.5 11.5 22. 5 τ 0 0 .2 0 .4 0 .6 0.8 2<(δPx)2>/<N><T> 2<(δPy)2>/<N><T> ( b) ( a) FIG. 1. 共Color online兲Scaled second moment of the fluctuations as a function of ␶ for a system with p=0.8. These results are from DSMC simulations and have been averaged over 4000 trajectories. DYNAMICS OF ANNIHILATION. II. FLUCTUATIONS OF…PHYSICAL REVIEW E 77, 051128 共2008兲 051128-7 in Fig. 8, where we have considered a system with p=0.5, they are correctly described by a Gaussian distribution. The figure displays the distribution at four different times, showing that the shape of the probability distributions does not vary during the dynamical evolution. Similar results have been obtained for the probability distribution of the total momentum. VI. TWO-TIME CORRELATION FUNCTION IN THE HDS In this section, we study the two-time correlation function of the global quantities in the HDS. With this aim, we consider two dynamical variables A共t兲and B共t兲as in 共33兲, and compute the correlations 具 ␦ A共t兲 ␦ B共t⬘兲典Hfor t⬎t⬘, which are obtained from h1,1,Hthrough Eq. 共35兲. As in the previous section, we start by integrating out the spatial degrees of freedom and consider ␺ H共 ␶ , ␶ ⬘,c1,c2兲⬅ 冕 dr12h ˜ H共r12;c1, ␶ ;c2, ␶ ⬘兲,共63兲 whose evolution equation is obtained by integrating 共32兲 over space variables, ⳵ ⳵ ␶␺ H共 ␶ , ␶ ⬘,c1,c2兲=⌳共c1兲 ␺ H共 ␶ , ␶ ⬘,c1,c2兲.共64兲 This equation has to be solved with the initial condition ␺ H共 ␶ ⬘, ␶ ⬘,c1,c2兲= ␹ H共c1兲 ␦ 共c1−c2兲+ ␾ H共 ␶ ⬘,c1,c2兲,共65兲 where we have taken into account the scaling of fH共18兲and g2,H共27兲. Then, using the approximation 共52兲, we obtain P12 ␺ H共 ␶ , ␶ ⬘,c1,c2兲 =e⌳共c1兲共 ␶ − ␶ ⬘兲P12 ␺ H共 ␶ ⬘, ␶ ⬘,c1,c2兲 =具 ␰ ¯ 1共c1兲兩 ␺ H共 ␶ ⬘, ␶ ⬘,c1,c2兲典 ␰ 1共c1兲 +具 ␰ ¯ 2共c1兲兩 ␺ H共 ␶ ⬘, ␶ ⬘,c1,c2兲典 ␰ 2共c1兲e−p共 ␨ T+2 ␨ n兲共 ␶ − ␶ ⬘兲 +兺 i 具 ␰ ¯ 3i共c1兲兩 ␺ H共 ␶ ⬘, ␶ ⬘,c1,c2兲典 ␰ 3i共c1兲ep ␨ T共 ␶ − ␶ ⬘兲. 共66兲 Note that, assuming that the hydrodynamic eigenvalues are the biggest ones, we can deduce in the long-time limit that P12 ␺ H共 ␶ , ␶ ⬘兲is a linear combination of the exponentials written above. Hence, the approximation 共52兲affects only the coefficients that multiply the exponentials, i.e., the initial conditions. In the large-time limit ␶ , ␶ ⬘→⬁, ␶ − ␶ ⬘finite 共and positive兲, we can replace ␺ H共 ␶ ⬘, ␶ ⬘,c1,c2兲by ␹ H共c1兲 ␦ 共c1−c2兲 + ␾ H s共c1,c2兲, so that 具 ␦ N共 ␶ 兲 ␦ N共 ␶ ⬘兲典H=NH共 ␶ 兲共−A1,1 −zA1,2 −zA1,2 ⫻e−p共 ␨ T+2 ␨ n兲共 ␶ − ␶ ⬘兲−z2A2,2e−p共 ␨ T+2 ␨ n兲共 ␶ − ␶ ⬘兲兲, 共67兲 0 0.2 0.4 0.6 0.8 <(δΝ)2>/<Ν> 0123 4 τ MD DSMC 0 0.5 1 1.5 2 2.5 3 <(δΕ)2>/<Ν><Τ>2 0123 4 τ MD DSMC (b) (a) FIG. 2. 共Color online兲Second moment of the fluctuations of the number of particles 共left panel兲and of the total energy 共right panel兲 as a function of the number of collisions per particle ␶ , for a system with p=0.5 and initial number of particles N=105. The dashed lines are the theoretical predictions. 0 0.2 0.4 0.6 0 . 8 <δΕδΝ>/<Ν><Τ> 0123 4 τ MD DSMC 0 0.2 0.4 0.6 0 . 8 2<(δPy) 2 >/<N><T> 0123 4 τ MD DSMC (b) (a) FIG. 3. 共Color online兲Correlation between the fluctuations of the total energy and total number of particles 共left panel兲, and second moment of the fluctuations of the ycomponent of the total momentum 共right panel兲, as a function of the number of collisions per particle ␶ , for a system with p=0.5 and initial number of particles N=105. The dashed lines are the theoretical predictions. MAYNAR et al. PHYSICAL REVIEW E 77, 051128 共2008兲 051128-8 具 ␦ Pi共 ␶ 兲 ␦ Pj共 ␶ ⬘兲典H= ␦ ijNH共 ␶ 兲vH共 ␶ 兲vH共 ␶ ⬘兲 冉 a3i,3i s+1 2 冊 ep ␨ T共 ␶ − ␶ ⬘兲, 共68兲 具 ␦ E共 ␶ 兲 ␦ E共 ␶ ⬘兲典H= 冉 dm 4 冊 2 NH共 ␶ 兲vH 2共 ␶ 兲vH 2共 ␶ ⬘兲兵−A1,1 +共z+2兲A1,2 +关−共z+2兲2A2,2 +共z+2兲A1,2兴e−p共 ␨ T+2 ␨ n兲共 ␶ − ␶ ⬘兲其,共69兲 具 ␦ N共 ␶ 兲 ␦ E共 ␶ ⬘兲典H=md 4NH共 ␶ 兲vH 2共 ␶ ⬘兲兵A1,1 −共z+2兲A1,2 +关−z共z+2兲A2,2 +zA1,2兴e−p共 ␨ T+2 ␨ n兲共 ␶ − ␶ ⬘兲其. 共70兲 In the ␶ scale, it can be seen from Eqs. 共20兲and 共21兲that NH and vHdecay exponentially. For A,B=N,E,P, the normalized correlation functions CAB共 ␶ , ␶ ⬘兲=具 ␦ A共 ␶ 兲 ␦ B共 ␶ ⬘兲典/具 ␦ A共 ␶ ⬘兲 ␦ B共 ␶ ⬘兲典 共71兲 become therefore time-translationally invariant, i.e., functions of ␶ − ␶ ⬘, once the stationary regime for the ratios such as 具 ␦ N2共 ␶ 兲典/NH共 ␶ 兲has been reached 共see the previous section兲. We have checked numerically that this is indeed the case, and we compare in Figs. 9and 10 the evolution of CNN共 ␶ − ␶ ⬘兲,CEE共 ␶ − ␶ ⬘兲, and CNE共 ␶ − ␶ ⬘兲measured in DSMC simulations 共for p=1, averaged over 4000 trajectories兲with the theoretical predictions. The agreement is very good. Figure 11 also shows the theoretical prediction for the decay of the momentum correlation function CPP共 ␶ , ␶ ⬘兲 =具 ␦ Pi共 ␶ 兲 ␦ Pj共 ␶ ⬘兲典/具 ␦ Pi共 ␶ ⬘兲 ␦ Pj共 ␶ ⬘兲典 for p=0.5. The characteristic decay time of CPP is of the order of ␶ ⬃4. Because 00.5 11.5 22. 5 τ 0 0 .5 1 1.5 2 2.5 3 <(δΝ)2>/<Ν> <(δΕ)2>/<Ν><Τ>2 <δΕδΝ>/<Ν><Τ> FIG. 5. 共Color online兲Same as in Fig. 4but for a system with p=0.8. 0 0.2 0.4 0.6 0.8 1 p 0 0.5 1 1.5 2 <(δΕ) 2 >Η/<Τ>Η 2 <Ν>Η DSMC MD 0 0.2 0.4 0.6 0.8 1 p 0 0.2 0.4 0.6 0.8 <(δΝ)2>Η/<Ν>Η DSMC MD ( b ) (a) FIG. 6. 共Color online兲Average steady state of the scaled second moment of the fluctuations of the number of particles and of the energy fluctuations as a function of the annihilation probability p. The dashed lines show the large- ␶ predictions of Eqs. 共59兲and 共61兲. 0 0 .5 1 1.5 2 2 .5 3 02.5 5 7.5 1 0 τ <(δΝ)2>/<Ν> <(δΕ)2>/<Ν><Τ>2 <δΕδΝ>/<Ν><Τ> 0123 4 τ 0 0 .5 1 1.5 2 2.5 3 <(δΝ)2>/<Ν> <(δΕ)2>/<Ν><Τ>2 <δΕδΝ>/<Ν><Τ> (b) (a) FIG. 4. 共Color online兲DSMC results for the scaled second moment of the fluctuations as a function of ␶ for systems with p=0.2 共left panel兲and 0.5 共right panel兲. The dashed lines are the theoretical predictions for the stationary values. The initial velocity distribution at ␶ =0 is uniform in a square domain. DYNAMICS OF ANNIHILATION. II. FLUCTUATIONS OF…PHYSICAL REVIEW E 77, 051128 共2008兲 051128-9