Dynamics of Annihilation I : Linearized Boltzmann Equation and Hydrodynamics
Physical Review E - Statistical, Nonlinear and Soft Matter Physics (2008)
- DOI: 10.1103/PhysRevE.77.051127
- PubMed: 18643046
- arXiv: 0801.2299
Available from arxiv.org
or
Abstract
We study the nonequilibrium statistical mechanics of a system of freely moving particles, in which binary encounters lead either to an elastic collision or to the disappearance of the pair. Such a system of ballistic annihilation therefore constantly loses particles. The dynamics of perturbations around the free decay regime is investigated using the spectral properties of the linearized Boltzmann operator, which characterize linear excitations on all time scales. The linearized Boltzmann equation is solved in the hydrodynamic limit by a projection technique, which yields the evolution equations for the relevant coarse-grained fields and expressions for the transport coefficients. We finally present the results of molecular dynamics simulations that validate the theoretical predictions.
Available from arxiv.org
Page 1
Dynamics of Annihilation I : Linearized Boltzmann Equation and Hydrodynamics
ar
X
iv
:0
80
1.
22
99
v1
[
co
nd
-m
at.
sta
t-m
ec
h]
1
5 J
an
20
08
Dynamics of Annihilation I : Linearized Boltzmann Equation and Hydrodynamics
Mar´ıa Isabel Garc´ıa de Soria,1 Pablo Maynar,2, 3 Gre´gory Schehr,2 Alain Barrat,2 and Emmanuel Trizac1
1Universite´ Paris-Sud, LPTMS, UMR 8626, Orsay Cedex, F-91405 and CNRS, Orsay, F-91405
2Laboratoire de Physique The´orique (CNRS UMR 8627),
Baˆtiment 210, Universite´ Paris-Sud, 91405 Orsay cedex, France
3F´ısica Teo´rica, Universidad de Sevilla, Apartado de Correos 1065, E-41080, Sevilla, Spain
(Dated: February 2, 2008)
We study the non-equilibrium statistical mechanics of a system of freely moving particles, in
which binary encounters lead either to an elastic collision or to the disappearance of the pair. Such
a system of ballistic annihilation therefore constantly looses particles. The dynamics of perturbations
around the free decay regime is investigated from the spectral properties of the linearized Boltzmann
operator, that characterize linear excitations on all time scales. The linearized Boltzmann equation
is solved in the hydrodynamic limit by a projection technique, which yields the evolution equations
for the relevant coarse-grained fields and expressions for the transport coefficients. We finally present
the results of Molecular Dynamics simulations that validate the theoretical predictions.
PACS numbers: 51.10.+y,05.20.Dd,82.20.Nk
I. INTRODUCTION
Understanding the differences and similarities between a flow of macroscopic grains and that of an ordinary liquid
is an active field of research [1, 2]. From a fundamental perspective, it is tempting to draw a correspondence between
the grains of the former and the atoms of the latter in order to make use of the powerful tools of statistical mechanics
to derive a large scale description for the various fields of interest, such as the local density of grains. A key difference
between a granular system and an ordinary liquid is that collisions between macroscopic grains dissipate energy, due to
the redistribution of translational kinetic energy into internal modes. This simple fact has far reaching consequences
[2, 3], but also poses an a priori serious problem concerning the validity of the procedure leading to the hydrodynamic
description. Indeed, the standard approach retains in the coarse-grained description only those fields associated with
quantities that are conserved in collisions (such as density and momentum). There is however good evidence –both
numerical and theoretical– that in the granular case, a relevant description should include the kinetic temperature
field, defined as the kinetic energy density ([2, 4] and references therein), which is therefore not associated to a
conserved quantity.
Our purpose here is to test a hydrodynamic description with suitable coarse-grained fields, for a model system where
not only the kinetic energy is not conserved during binary encounters, but also the number of particles and the linear
momentum. The ballistic annihilation model [5–10] provides a valuable candidate: in this model, each particle moves
freely (ballistically) until it meets another particle; such binary encounters lead to the annihilation of the colliding
pair of particles. In addition, we introduce a parameter 0 ≤ p ≤ 1 that may be thought of as a measure of the distance
to equilibrium, so that an ensemble of spherical particles in dimension d undergoing ballistic motion either annihilate
upon contact (with probability p) or scatter elastically (with probability 1− p). For the corresponding probabilistic
ballistic annihilation model, the Chapman-Enskog [11] scheme was applied recently [12]. The hydrodynamic equations
were derived and explicit formulas for the transport coefficients obtained. Our goal here is two-fold. First, we would
like to shed light on the context and limitations of the derivation, by obtaining the hydrodynamic description directly
from the linearized Boltzmann equation. Second, we aim at putting to the test the theoretical framework thereby
obtained by careful comparison with numerical simulations of the annihilation process. For granular gas dynamics,
the same program is quite complete, although challenges remain [1, 2]. The objective here is to initiate a similar
formulation for the ballistic annihilation model in view of a more stringent test of the hydrodynamic machinery.
The paper is organized as follows. We start in section II with a reminder of results derived in Refs [8, 9]. The
kinetic description adopted is that of the Boltzmann equation, since it has been shown that for p = 1 (all collision
events leading to annihilation), the underlying molecular chaos closure provides an exact description at long times,
provided space dimension d is strictly larger than 1 [9]. We may assume that the same holds for an arbitrary but
non vanishing value of p, since the density is then still a decreasing function of time. The focus is here on an
unforced system, which is characterized by an algebraic decay with time of the total density and kinetic energy
density (homogeneous decay state) [8, 9]. More precisely, we are interested in small perturbations around this state,
so that the Boltzmann equation will be subsequently linearized. After having identified the operator that generates
the dynamics of fluctuations, attention will be paid in section III to its spectral properties. This will provide the
basis for finding in section IV the evolution equations for the hydrodynamic fields (i.e. those chosen for the coarse-
X
iv
:0
80
1.
22
99
v1
[
co
nd
-m
at.
sta
t-m
ec
h]
1
5 J
an
20
08
Dynamics of Annihilation I : Linearized Boltzmann Equation and Hydrodynamics
Mar´ıa Isabel Garc´ıa de Soria,1 Pablo Maynar,2, 3 Gre´gory Schehr,2 Alain Barrat,2 and Emmanuel Trizac1
1Universite´ Paris-Sud, LPTMS, UMR 8626, Orsay Cedex, F-91405 and CNRS, Orsay, F-91405
2Laboratoire de Physique The´orique (CNRS UMR 8627),
Baˆtiment 210, Universite´ Paris-Sud, 91405 Orsay cedex, France
3F´ısica Teo´rica, Universidad de Sevilla, Apartado de Correos 1065, E-41080, Sevilla, Spain
(Dated: February 2, 2008)
We study the non-equilibrium statistical mechanics of a system of freely moving particles, in
which binary encounters lead either to an elastic collision or to the disappearance of the pair. Such
a system of ballistic annihilation therefore constantly looses particles. The dynamics of perturbations
around the free decay regime is investigated from the spectral properties of the linearized Boltzmann
operator, that characterize linear excitations on all time scales. The linearized Boltzmann equation
is solved in the hydrodynamic limit by a projection technique, which yields the evolution equations
for the relevant coarse-grained fields and expressions for the transport coefficients. We finally present
the results of Molecular Dynamics simulations that validate the theoretical predictions.
PACS numbers: 51.10.+y,05.20.Dd,82.20.Nk
I. INTRODUCTION
Understanding the differences and similarities between a flow of macroscopic grains and that of an ordinary liquid
is an active field of research [1, 2]. From a fundamental perspective, it is tempting to draw a correspondence between
the grains of the former and the atoms of the latter in order to make use of the powerful tools of statistical mechanics
to derive a large scale description for the various fields of interest, such as the local density of grains. A key difference
between a granular system and an ordinary liquid is that collisions between macroscopic grains dissipate energy, due to
the redistribution of translational kinetic energy into internal modes. This simple fact has far reaching consequences
[2, 3], but also poses an a priori serious problem concerning the validity of the procedure leading to the hydrodynamic
description. Indeed, the standard approach retains in the coarse-grained description only those fields associated with
quantities that are conserved in collisions (such as density and momentum). There is however good evidence –both
numerical and theoretical– that in the granular case, a relevant description should include the kinetic temperature
field, defined as the kinetic energy density ([2, 4] and references therein), which is therefore not associated to a
conserved quantity.
Our purpose here is to test a hydrodynamic description with suitable coarse-grained fields, for a model system where
not only the kinetic energy is not conserved during binary encounters, but also the number of particles and the linear
momentum. The ballistic annihilation model [5–10] provides a valuable candidate: in this model, each particle moves
freely (ballistically) until it meets another particle; such binary encounters lead to the annihilation of the colliding
pair of particles. In addition, we introduce a parameter 0 ≤ p ≤ 1 that may be thought of as a measure of the distance
to equilibrium, so that an ensemble of spherical particles in dimension d undergoing ballistic motion either annihilate
upon contact (with probability p) or scatter elastically (with probability 1− p). For the corresponding probabilistic
ballistic annihilation model, the Chapman-Enskog [11] scheme was applied recently [12]. The hydrodynamic equations
were derived and explicit formulas for the transport coefficients obtained. Our goal here is two-fold. First, we would
like to shed light on the context and limitations of the derivation, by obtaining the hydrodynamic description directly
from the linearized Boltzmann equation. Second, we aim at putting to the test the theoretical framework thereby
obtained by careful comparison with numerical simulations of the annihilation process. For granular gas dynamics,
the same program is quite complete, although challenges remain [1, 2]. The objective here is to initiate a similar
formulation for the ballistic annihilation model in view of a more stringent test of the hydrodynamic machinery.
The paper is organized as follows. We start in section II with a reminder of results derived in Refs [8, 9]. The
kinetic description adopted is that of the Boltzmann equation, since it has been shown that for p = 1 (all collision
events leading to annihilation), the underlying molecular chaos closure provides an exact description at long times,
provided space dimension d is strictly larger than 1 [9]. We may assume that the same holds for an arbitrary but
non vanishing value of p, since the density is then still a decreasing function of time. The focus is here on an
unforced system, which is characterized by an algebraic decay with time of the total density and kinetic energy
density (homogeneous decay state) [8, 9]. More precisely, we are interested in small perturbations around this state,
so that the Boltzmann equation will be subsequently linearized. After having identified the operator that generates
the dynamics of fluctuations, attention will be paid in section III to its spectral properties. This will provide the
basis for finding in section IV the evolution equations for the hydrodynamic fields (i.e. those chosen for the coarse-
Page 2
2grained description) and for obtaining explicit formulas for the transport coefficients. Finally, our predictions will be
confronted in section V against extensive Molecular Dynamics simulations. Such a comparison is an essential step in
testing the foundations of the hydrodynamic treatment.
II. THE BOLTZMANN EQUATION APPROACH TO THE HOMOGENEOUS DECAY STATE
A. Non-linear description
The Boltzmann equation describes the time evolution of the one particle distribution function f(r,v1, t). For a
system of smooth hard disks or spheres of mass m and diameter σ, which annihilate with probability p or collide
elastically with probability 1− p, it has the form
(
∂
∂t
+ v1 · ∇
)
f(r,v1, t) = pJa[f |f ] + (1− p)Jc[f |f ], (1)
where the annihilation operator Ja is defined by [9]
Ja[f |g] = −σd−1
∫
dv2
∫
dσˆΘ(v12 · σˆ)(v12 · σˆ)f(r,v1, t)g(r,v2, t). (2)
The elastic collision operator Jc reads [13, 14]
Jc[f |g] = σd−1
∫
dv2
∫
dσˆΘ(v12 · σˆ)(v12 · σˆ)(b−1σ − 1)f(r,v1, t)g(r,v2, t), (3)
with v12 = v1 − v2, Θ the Heaviside step function, σˆ a unit vector joining the centers of the two particles at contact
and b−1σ an operator replacing all the velocities v1 and v2 appearing in its argument by their precollisional values v
∗
1
and v∗2 , given by
b−1σ v1 = v
∗
1 = v1 − (v12 · σˆ)σˆ, (4)
b−1σ v2 = v
∗
2 = v2 + (v12 · σˆ)σˆ. (5)
We assume that the system can be characterized macroscopically by coarse grained (hydrodynamic-like) fields, that
we define as in standard Kinetic Theory in terms of the local velocity distribution function f(r,v, t)
n(r, t) =
∫
dvf(r,v, t), (6)
n(r, t)u(r, t) =
∫
dvvf(r,v, t), (7)
d
2
n(r, t)T (r, t) =
∫
dv
m
2
V 2f(r,v, t), (8)
where n(r, t), u(r, t), and T (r, t) are the local number density, velocity, and temperature, respectively. We have
introduced here V = v − u, the velocity of the particle relative to the local velocity flow. We stress that the
temperature defined has a kinetic meaning only, but lacks a thermodynamic interpretation. It seems natural to
consider these fields, as they are the usual hydrodynamical fields of the equilibrium system (with p = 0). It is however
not obvious at this point that restricting our coarse-grained description to the above three fields provides a relevant
and consistent framework. A major goal of this paper is to provide strong hints that this is indeed the case. We will
show in particular that closed equations can be obtained for these fields in the appropriate time and length scales,
under reasonable assumptions.
The Boltzmann equation (1) admits a homogeneous scaling solution fH in which all the time dependence is em-
bedded in the hydrodynamic fields, with the further simplification that those fields are position independent. The
existence of this regime could not be shown rigorously, but, numerically, such a scaling solution quickly emerges from
an arbitrary initial condition [8, 9]. It has the form [9]
fH(v, t) =
nH(t)
vH(t)d
χH(c), with c =
v
vH(t)
, (9)
confronted in section V against extensive Molecular Dynamics simulations. Such a comparison is an essential step in
testing the foundations of the hydrodynamic treatment.
II. THE BOLTZMANN EQUATION APPROACH TO THE HOMOGENEOUS DECAY STATE
A. Non-linear description
The Boltzmann equation describes the time evolution of the one particle distribution function f(r,v1, t). For a
system of smooth hard disks or spheres of mass m and diameter σ, which annihilate with probability p or collide
elastically with probability 1− p, it has the form
(
∂
∂t
+ v1 · ∇
)
f(r,v1, t) = pJa[f |f ] + (1− p)Jc[f |f ], (1)
where the annihilation operator Ja is defined by [9]
Ja[f |g] = −σd−1
∫
dv2
∫
dσˆΘ(v12 · σˆ)(v12 · σˆ)f(r,v1, t)g(r,v2, t). (2)
The elastic collision operator Jc reads [13, 14]
Jc[f |g] = σd−1
∫
dv2
∫
dσˆΘ(v12 · σˆ)(v12 · σˆ)(b−1σ − 1)f(r,v1, t)g(r,v2, t), (3)
with v12 = v1 − v2, Θ the Heaviside step function, σˆ a unit vector joining the centers of the two particles at contact
and b−1σ an operator replacing all the velocities v1 and v2 appearing in its argument by their precollisional values v
∗
1
and v∗2 , given by
b−1σ v1 = v
∗
1 = v1 − (v12 · σˆ)σˆ, (4)
b−1σ v2 = v
∗
2 = v2 + (v12 · σˆ)σˆ. (5)
We assume that the system can be characterized macroscopically by coarse grained (hydrodynamic-like) fields, that
we define as in standard Kinetic Theory in terms of the local velocity distribution function f(r,v, t)
n(r, t) =
∫
dvf(r,v, t), (6)
n(r, t)u(r, t) =
∫
dvvf(r,v, t), (7)
d
2
n(r, t)T (r, t) =
∫
dv
m
2
V 2f(r,v, t), (8)
where n(r, t), u(r, t), and T (r, t) are the local number density, velocity, and temperature, respectively. We have
introduced here V = v − u, the velocity of the particle relative to the local velocity flow. We stress that the
temperature defined has a kinetic meaning only, but lacks a thermodynamic interpretation. It seems natural to
consider these fields, as they are the usual hydrodynamical fields of the equilibrium system (with p = 0). It is however
not obvious at this point that restricting our coarse-grained description to the above three fields provides a relevant
and consistent framework. A major goal of this paper is to provide strong hints that this is indeed the case. We will
show in particular that closed equations can be obtained for these fields in the appropriate time and length scales,
under reasonable assumptions.
The Boltzmann equation (1) admits a homogeneous scaling solution fH in which all the time dependence is em-
bedded in the hydrodynamic fields, with the further simplification that those fields are position independent. The
existence of this regime could not be shown rigorously, but, numerically, such a scaling solution quickly emerges from
an arbitrary initial condition [8, 9]. It has the form [9]
fH(v, t) =
nH(t)
vH(t)d
χH(c), with c =
v
vH(t)
, (9)
Page 3
3where
vH(t) =
[
2TH(t)
m
]1/2
(10)
is the “thermal” (root-mean-square) velocity and χH(c) is an isotropic function depending only on the modulus c = |c|
of the rescaled velocity. By taking moments in the Boltzmann equation and using the scaling (9), it can be seen that
the homogeneous density and temperature obey the equations [12]
∂nH(t)
∂t
= −pνH(t)ζnnH(t), (11)
∂TH(t)
∂t
= −pνH(t)ζTTH(t), (12)
where we have introduced the collision frequency of the corresponding hard sphere fluid in equilibrium (with same
temperature and density)
νH(t) =
nH(t)T
1/2
H (t)σ
d−1
m1/2
8pi
d−1
2
(d + 2)Γ(d/2)
(13)
and the dimensionless decay rates ζn and ζT , that are functionals of the distribution function
pζn = −
γ
2
∫
dc1
∫
dc2T (c1, c2)χH(c1)χH(c2), (14)
pζT = −
γ
2
∫
dc1
∫
dc2
(
2c21
d
− 1
)
T (c1, c2)χH(c1)χH(c2). (15)
In these expressions, γ is a quantity that does not depend on time, which reads
γ =
2nH(t)vH(t)σd−1
νH(t)
=
(d + 2)
√
2Γ(d/2)
4pi(d−1)/2
, (16)
and the binary collision operator T (c1, c2), that should not be confused with the temperature, takes the form
T (c1, c2) =
∫
dσˆΘ(c12 · σˆ)(c12 · σˆ)[(1 − p)b−1σ − 1]. (17)
Finally, we can write an equation for the scaled distribution function χH(c) in terms of the coefficients and operators
defined above
p
[
(dζT − 2ζn) + ζT c1 ·
∂
∂c1
]
χH(c1) = γ
∫
dc2T (c1, c2)χH(c1)χH(c2). (18)
The operator b−1σ in the last equation is defined again by equation (4), but substituting (v1,v2) by (c1, c2).
Although an exact and explicit solution of equation (18) is not known, its behavior at large and small velocities has
been determined [8, 9]. In this work we will use the approximate form of the distribution function in the so-called
first Sonine approximation (an expansion around a Gaussian functional form, see Appendix A), which is valid for
velocities in the thermal region, and all the functionals of χH(c), that is the decay rates and the transport coefficients,
will be evaluated in this approximation [8, 15].
B. Linearized Boltzmann Equation
In the remainder, we consider a situation where the system is very close to the homogeneous decay state, so that
we can write
f(r,v1, t) = fH(v1, t) + δf(r,v1, t), |δf(r,v1, t)| ≪ fH(v1, t). (19)
Substitution of equation (19) into the Boltzmann equation (1), keeping only linear terms in δf , yields
(
∂
∂t
+ v1 · ∇
)
δf(r,v1, t)
= p {Ja[δf |fH ] + Ja[fH |δf ]}+ (1− p) {Jc[δf |fH ] + Jc[fH |δf ]} , (20)
vH(t) =
[
2TH(t)
m
]1/2
(10)
is the “thermal” (root-mean-square) velocity and χH(c) is an isotropic function depending only on the modulus c = |c|
of the rescaled velocity. By taking moments in the Boltzmann equation and using the scaling (9), it can be seen that
the homogeneous density and temperature obey the equations [12]
∂nH(t)
∂t
= −pνH(t)ζnnH(t), (11)
∂TH(t)
∂t
= −pνH(t)ζTTH(t), (12)
where we have introduced the collision frequency of the corresponding hard sphere fluid in equilibrium (with same
temperature and density)
νH(t) =
nH(t)T
1/2
H (t)σ
d−1
m1/2
8pi
d−1
2
(d + 2)Γ(d/2)
(13)
and the dimensionless decay rates ζn and ζT , that are functionals of the distribution function
pζn = −
γ
2
∫
dc1
∫
dc2T (c1, c2)χH(c1)χH(c2), (14)
pζT = −
γ
2
∫
dc1
∫
dc2
(
2c21
d
− 1
)
T (c1, c2)χH(c1)χH(c2). (15)
In these expressions, γ is a quantity that does not depend on time, which reads
γ =
2nH(t)vH(t)σd−1
νH(t)
=
(d + 2)
√
2Γ(d/2)
4pi(d−1)/2
, (16)
and the binary collision operator T (c1, c2), that should not be confused with the temperature, takes the form
T (c1, c2) =
∫
dσˆΘ(c12 · σˆ)(c12 · σˆ)[(1 − p)b−1σ − 1]. (17)
Finally, we can write an equation for the scaled distribution function χH(c) in terms of the coefficients and operators
defined above
p
[
(dζT − 2ζn) + ζT c1 ·
∂
∂c1
]
χH(c1) = γ
∫
dc2T (c1, c2)χH(c1)χH(c2). (18)
The operator b−1σ in the last equation is defined again by equation (4), but substituting (v1,v2) by (c1, c2).
Although an exact and explicit solution of equation (18) is not known, its behavior at large and small velocities has
been determined [8, 9]. In this work we will use the approximate form of the distribution function in the so-called
first Sonine approximation (an expansion around a Gaussian functional form, see Appendix A), which is valid for
velocities in the thermal region, and all the functionals of χH(c), that is the decay rates and the transport coefficients,
will be evaluated in this approximation [8, 15].
B. Linearized Boltzmann Equation
In the remainder, we consider a situation where the system is very close to the homogeneous decay state, so that
we can write
f(r,v1, t) = fH(v1, t) + δf(r,v1, t), |δf(r,v1, t)| ≪ fH(v1, t). (19)
Substitution of equation (19) into the Boltzmann equation (1), keeping only linear terms in δf , yields
(
∂
∂t
+ v1 · ∇
)
δf(r,v1, t)
= p {Ja[δf |fH ] + Ja[fH |δf ]}+ (1− p) {Jc[δf |fH ] + Jc[fH |δf ]} , (20)
Page 4
4Given the scaling form of f (Eq. (9)), it is convenient to introduce as well the scaled deviation of the distribution
function, δχ, as follows
δf(r,v1, t) =
nH(t)
vH(t)d
δχ(r, c1, τ) . (21)
Moreover, Eqs. (11) and (12) suggest to use the dimensionless time scale τ defined by
τ =
1
2
∫ t
0
dt′νH(t′), (22)
which counts the number of collision per particle in the time interval [0, t]. Combining Eq. (13) together with Eq.
(11, 12) yields immediately νH(t) = (1/νH(0) + p(ζn + ζT /2)t)−1 and thus
τ =
1
p(2ζn + ζT )
log[1 + νH(0)p(ζn + ζT /2)t] . (23)
In this time scale τ (22), these equations (11, 12) are easily integrated, yielding
nH(τ) = nH(0) exp(−2pζnτ), TH(τ) = TH(0) exp(−2pζT τ), (24)
and power law behaviors in time t, nH(t) ∝ t−2ζn/(2ζn+ζT ) and TH(t) ∝ t−2ζT /(2ζn+ζT ) at large time t ≫ 1 . It proves
also convenient to introduce Fourier components [with the notation hk =
∫
dr exp−ik·r h(r)] so that the evolution
equation for a general k component of δχ is, in the τ timescale,
∂
∂τ
δχk(c1, τ) = [Λ(c1)− ilH(τ)k · c1] δχk(c1, τ). (25)
In this equation, the time dependent length scale lH = 2vH(τ)/νH(τ) is proportional to the instantaneous mean free
path (lH(τ) ∝ n−1H (τ), see Eq. (13)) and the homogeneous scaled Boltzmann linear operator reads
Λ(c1)h(c1) = γ
∫
dc2T (c1, c2)(1 + P12)χH(c1)h(c2)
+p(2ζn − dζT )h(c1)− pζT c1 ·
∂
∂c1
h(c1). (26)
In this expression, the permutation operator P12 interchanges the labels of particles 1 and 2 and subsequently allows
for more compact notations. In the present representation, all the time dependence due to the reference state is
absorbed in the mean free path, obtained from lH(τ) ∝ n−1H (τ) as
lH(τ) = lH(0) exp(2pζnτ), (27)
which, as expected, is an increasing function of time.
C. Linearized Hydrodynamic Equations around the homogeneous decay state
Let us define the relative deviations of the hydrodynamic fields from their homogeneous values by
ρ(r, τ) ≡ δn(r, τ)
nH(τ)
=
∫
dcδχ(r, c, τ), (28)
w(r, τ) ≡ δu(r, τ)
vH(τ)
=
∫
dccδχ(r, c, τ), (29)
θ(r, τ) ≡ δT (r, τ)
TH(τ)
=
∫
dc
(
2c2
d
− 1
)
δχ(r, c, τ), (30)
where δy(r, τ) ≡ y(r, τ) − yH(t) denotes the deviation of a local macroscopic variable, y(r, τ), from its homogeneous
decay state value, yH(t). Taking velocity moments in the Boltzmann equation (25), we obtain the linearized balance
function, δχ, as follows
δf(r,v1, t) =
nH(t)
vH(t)d
δχ(r, c1, τ) . (21)
Moreover, Eqs. (11) and (12) suggest to use the dimensionless time scale τ defined by
τ =
1
2
∫ t
0
dt′νH(t′), (22)
which counts the number of collision per particle in the time interval [0, t]. Combining Eq. (13) together with Eq.
(11, 12) yields immediately νH(t) = (1/νH(0) + p(ζn + ζT /2)t)−1 and thus
τ =
1
p(2ζn + ζT )
log[1 + νH(0)p(ζn + ζT /2)t] . (23)
In this time scale τ (22), these equations (11, 12) are easily integrated, yielding
nH(τ) = nH(0) exp(−2pζnτ), TH(τ) = TH(0) exp(−2pζT τ), (24)
and power law behaviors in time t, nH(t) ∝ t−2ζn/(2ζn+ζT ) and TH(t) ∝ t−2ζT /(2ζn+ζT ) at large time t ≫ 1 . It proves
also convenient to introduce Fourier components [with the notation hk =
∫
dr exp−ik·r h(r)] so that the evolution
equation for a general k component of δχ is, in the τ timescale,
∂
∂τ
δχk(c1, τ) = [Λ(c1)− ilH(τ)k · c1] δχk(c1, τ). (25)
In this equation, the time dependent length scale lH = 2vH(τ)/νH(τ) is proportional to the instantaneous mean free
path (lH(τ) ∝ n−1H (τ), see Eq. (13)) and the homogeneous scaled Boltzmann linear operator reads
Λ(c1)h(c1) = γ
∫
dc2T (c1, c2)(1 + P12)χH(c1)h(c2)
+p(2ζn − dζT )h(c1)− pζT c1 ·
∂
∂c1
h(c1). (26)
In this expression, the permutation operator P12 interchanges the labels of particles 1 and 2 and subsequently allows
for more compact notations. In the present representation, all the time dependence due to the reference state is
absorbed in the mean free path, obtained from lH(τ) ∝ n−1H (τ) as
lH(τ) = lH(0) exp(2pζnτ), (27)
which, as expected, is an increasing function of time.
C. Linearized Hydrodynamic Equations around the homogeneous decay state
Let us define the relative deviations of the hydrodynamic fields from their homogeneous values by
ρ(r, τ) ≡ δn(r, τ)
nH(τ)
=
∫
dcδχ(r, c, τ), (28)
w(r, τ) ≡ δu(r, τ)
vH(τ)
=
∫
dccδχ(r, c, τ), (29)
θ(r, τ) ≡ δT (r, τ)
TH(τ)
=
∫
dc
(
2c2
d
− 1
)
δχ(r, c, τ), (30)
where δy(r, τ) ≡ y(r, τ) − yH(t) denotes the deviation of a local macroscopic variable, y(r, τ), from its homogeneous
decay state value, yH(t). Taking velocity moments in the Boltzmann equation (25), we obtain the linearized balance
Page 5
5equation for the k components of the hydrodynamic fields
(
∂
∂τ
− 2pζn
)
ρk + ilH(τ)k ·wk − p δζn[δχk] = 0, (31)
[
∂
∂τ
− p(2ζn + ζT )
]
wk
+
i
2
lH(τ)k(ρk + θk) + ilH(τ)k ·Π[δχk]− p δζu[δχk] = 0, (32)
[
∂
∂τ
− 2p(ζn + ζT )
]
θk
−2pζTρk + i
2
d
lH(τ)k · (wk + φ[δχk])− p δζT [δχk] = 0. (33)
Here, we have introduced the traceless pressure tensor and the heat flux as
Π[δχk] =
∫
dc∆(c)δχk(c, τ), (34)
φ[δχk] =
∫
dcΣ(c)δχk(c, τ), (35)
where ∆ and Σ are defined as
∆ij(c) = cicj −
c2
d
δij , (36)
Σ(c) =
(
c2 − d + 2
2
)
c, (37)
and the functionals
p δζn[δχ] = γ
∫
dc1
∫
dc2T (c1, c2)(1 + P12)χH(c1)δχk(c2, τ), (38)
p δζu[δχ] = γ
∫
dc1
∫
dc2c1T (c1, c2)(1 + P12)χH(c1)δχk(c2, τ), (39)
p δζT [δχ] = γ
∫
dc1
∫
dc2
(
2c21
d
− 1
)
T (c1, c2)(1 + P12)χH(c1)δχk(c2, τ). (40)
The previous analysis therefore amounts to obtaining a set of complicated equations expressing the evolution of the
hydrodynamic fields as a function of the rescaled homogeneous distribution function χH and the perturbation δχ. In
order to obtain a closed set of equations for the hydrodynamic fields (31), (32), (33), we need therefore to express the
functionals Π, φ, δζn, δζu and δζT , in terms of the hydrodynamic fields themselves. We will see in the next section
that, as long as we can treat lH(τ)k as a small parameter and if the linear Boltzmann operator has some specific
properties, it is possible to carry out this program and to close the linear hydrodynamic equations. However, since
the mean free path lH(τ) increases with time (Eq. (27)), the requirement of a small lH(τ)k is necessarily limited to a
time window depending on both k and the probability of annihilation p. An upper bound for this window is provided
by the time when the mean free path becomes of the order of the system size.
III. SOLUTION OF THE LINEARIZED BOLTZMANN EQUATION
In this Section we explore the solutions to the linearized Boltzmann equation (25) and establish some properties
of the homogeneous linear Boltzmann operator that will be essential for the coarse-grained description. From the
expression of the linearized Boltzmann equation, we can identify the operator Λ− ik · clH(τ) as the “generator of the
dynamic” of δχk. As we are interested in the solutions of this equation in the hydrodynamic regime (large enough
scales), it is convenient to study first the eigenvalue problem of the homogeneous linear Boltzmann operator. The
inhomogeneous term will be treated perturbatively later on.
(
∂
∂τ
− 2pζn
)
ρk + ilH(τ)k ·wk − p δζn[δχk] = 0, (31)
[
∂
∂τ
− p(2ζn + ζT )
]
wk
+
i
2
lH(τ)k(ρk + θk) + ilH(τ)k ·Π[δχk]− p δζu[δχk] = 0, (32)
[
∂
∂τ
− 2p(ζn + ζT )
]
θk
−2pζTρk + i
2
d
lH(τ)k · (wk + φ[δχk])− p δζT [δχk] = 0. (33)
Here, we have introduced the traceless pressure tensor and the heat flux as
Π[δχk] =
∫
dc∆(c)δχk(c, τ), (34)
φ[δχk] =
∫
dcΣ(c)δχk(c, τ), (35)
where ∆ and Σ are defined as
∆ij(c) = cicj −
c2
d
δij , (36)
Σ(c) =
(
c2 − d + 2
2
)
c, (37)
and the functionals
p δζn[δχ] = γ
∫
dc1
∫
dc2T (c1, c2)(1 + P12)χH(c1)δχk(c2, τ), (38)
p δζu[δχ] = γ
∫
dc1
∫
dc2c1T (c1, c2)(1 + P12)χH(c1)δχk(c2, τ), (39)
p δζT [δχ] = γ
∫
dc1
∫
dc2
(
2c21
d
− 1
)
T (c1, c2)(1 + P12)χH(c1)δχk(c2, τ). (40)
The previous analysis therefore amounts to obtaining a set of complicated equations expressing the evolution of the
hydrodynamic fields as a function of the rescaled homogeneous distribution function χH and the perturbation δχ. In
order to obtain a closed set of equations for the hydrodynamic fields (31), (32), (33), we need therefore to express the
functionals Π, φ, δζn, δζu and δζT , in terms of the hydrodynamic fields themselves. We will see in the next section
that, as long as we can treat lH(τ)k as a small parameter and if the linear Boltzmann operator has some specific
properties, it is possible to carry out this program and to close the linear hydrodynamic equations. However, since
the mean free path lH(τ) increases with time (Eq. (27)), the requirement of a small lH(τ)k is necessarily limited to a
time window depending on both k and the probability of annihilation p. An upper bound for this window is provided
by the time when the mean free path becomes of the order of the system size.
III. SOLUTION OF THE LINEARIZED BOLTZMANN EQUATION
In this Section we explore the solutions to the linearized Boltzmann equation (25) and establish some properties
of the homogeneous linear Boltzmann operator that will be essential for the coarse-grained description. From the
expression of the linearized Boltzmann equation, we can identify the operator Λ− ik · clH(τ) as the “generator of the
dynamic” of δχk. As we are interested in the solutions of this equation in the hydrodynamic regime (large enough
scales), it is convenient to study first the eigenvalue problem of the homogeneous linear Boltzmann operator. The
inhomogeneous term will be treated perturbatively later on.
Page 6
6A. Hydrodynamic Eigenfunctions of Λ
Let us consider the eigenvalue problem of the homogeneous linear Boltzmann operator Λ
Λ(c)ξβ(c) = λβξβ(c). (41)
Finding all the solutions of this equation is an insurmountable task. Nevertheless, it is possible to obtain some
particular solutions, which will turn out to be the relevant ones in the hydrodynamic regime. The problem will be
posed in a Hilbert space of functions of c with scalar product given by
〈g|h〉 =
∫
dcχ−1H (c)g
∗(c)h(c), (42)
where g∗ denotes the complex conjugate of g.
Of particular interest here are the eigenfunctions and eigenvalues associated with linear hydrodynamics. Following
[16, 17], we use the fact that the homogeneous decay state is parameterized by the hydrodynamic fields nH , TH and
uH . Writing the Boltzmann equation satisfied by χH and differentiating it with respect to these fields allows then to
obtain three exact relations from which one can extract eigenfunctions of the linearized Boltzmann collision operator.
In Appendix B, we show that the functions
ξ1(c) = χH(c) +
∂
∂c
· [cχH(c)] , (43)
ξ2(c) = zχH(c)−
∂
∂c
· [cχH(c)] , (44)
ξ3(c) = −
∂
∂c
χH(c), (45)
with z = 2ζn/ζT , are solutions of Eq. (41), with eigenvalues
λ1 = 0, λ2 = −p(ζT + 2ζn), λ3 = pζT , (46)
respectively, λ3 being d-fold degenerate. Although we cannot prove in general that these eigenvalues are indeed the
hydrodynamic ones (i.e the upper part of the spectrum), we will assume that this is the case ; the self-consistency of
the approach and comparison with numerical simulations will validate this assumption. Interestingly, in the particular
case of Maxwell molecules where the full spectrum of Λ may be computed exactly (see Appendix C), it appears that
the above “hydrodynamic” modes dominate at long times, provided that p < 1/4. For larger values of p, the “kinetic”
mode with largest eigenvalue decays slower than one of the three “hydrodynamic” modes.
As a consequence of the non-hermitian character of the operator Λ, the functions {ξβ}β=1,...,3 are not orthogonal
with respect to the scalar product defined in (42). They are nevertheless independent and, in order to define the
projection onto the subspace spanned by these functions, it is necessary to introduce a set of functions {ξ¯β}β=1,...,3
verifying the biorthonormality condition
〈ξ¯β |ξβ′〉 = δβ,β′. (47)
Although the set {ξ¯β}β=1,...,3 is not unique, a convenient choice is given by
ξ¯1(c) =
[
2 + z
2(1 + z)
− z
1 + z
c2
d
]
χH(c), (48)
ξ¯2(c) =
[
1
2(1 + z)
+
1
1 + z
c2
d
]
χH(c), (49)
ξ¯3(c) = cχH(c). (50)
Indeed, the functions {ξ¯β}β=1,...,3 have to be linear combinations of χH(c), cχH(c) and c2χH(c) to ensure that
projection of δχk onto the {ξ¯β} yields the coarse-grained fields ρk, θk and wk, or combinations thereof. The functions
{ξ¯β}β=1,...,3 span a dual subspace of that spanned by the eigenfunctions and for any linear combination of the
hydrodynamic modes
g(c) =
3
∑
β=1
aβξβ(c), (51)
Let us consider the eigenvalue problem of the homogeneous linear Boltzmann operator Λ
Λ(c)ξβ(c) = λβξβ(c). (41)
Finding all the solutions of this equation is an insurmountable task. Nevertheless, it is possible to obtain some
particular solutions, which will turn out to be the relevant ones in the hydrodynamic regime. The problem will be
posed in a Hilbert space of functions of c with scalar product given by
〈g|h〉 =
∫
dcχ−1H (c)g
∗(c)h(c), (42)
where g∗ denotes the complex conjugate of g.
Of particular interest here are the eigenfunctions and eigenvalues associated with linear hydrodynamics. Following
[16, 17], we use the fact that the homogeneous decay state is parameterized by the hydrodynamic fields nH , TH and
uH . Writing the Boltzmann equation satisfied by χH and differentiating it with respect to these fields allows then to
obtain three exact relations from which one can extract eigenfunctions of the linearized Boltzmann collision operator.
In Appendix B, we show that the functions
ξ1(c) = χH(c) +
∂
∂c
· [cχH(c)] , (43)
ξ2(c) = zχH(c)−
∂
∂c
· [cχH(c)] , (44)
ξ3(c) = −
∂
∂c
χH(c), (45)
with z = 2ζn/ζT , are solutions of Eq. (41), with eigenvalues
λ1 = 0, λ2 = −p(ζT + 2ζn), λ3 = pζT , (46)
respectively, λ3 being d-fold degenerate. Although we cannot prove in general that these eigenvalues are indeed the
hydrodynamic ones (i.e the upper part of the spectrum), we will assume that this is the case ; the self-consistency of
the approach and comparison with numerical simulations will validate this assumption. Interestingly, in the particular
case of Maxwell molecules where the full spectrum of Λ may be computed exactly (see Appendix C), it appears that
the above “hydrodynamic” modes dominate at long times, provided that p < 1/4. For larger values of p, the “kinetic”
mode with largest eigenvalue decays slower than one of the three “hydrodynamic” modes.
As a consequence of the non-hermitian character of the operator Λ, the functions {ξβ}β=1,...,3 are not orthogonal
with respect to the scalar product defined in (42). They are nevertheless independent and, in order to define the
projection onto the subspace spanned by these functions, it is necessary to introduce a set of functions {ξ¯β}β=1,...,3
verifying the biorthonormality condition
〈ξ¯β |ξβ′〉 = δβ,β′. (47)
Although the set {ξ¯β}β=1,...,3 is not unique, a convenient choice is given by
ξ¯1(c) =
[
2 + z
2(1 + z)
− z
1 + z
c2
d
]
χH(c), (48)
ξ¯2(c) =
[
1
2(1 + z)
+
1
1 + z
c2
d
]
χH(c), (49)
ξ¯3(c) = cχH(c). (50)
Indeed, the functions {ξ¯β}β=1,...,3 have to be linear combinations of χH(c), cχH(c) and c2χH(c) to ensure that
projection of δχk onto the {ξ¯β} yields the coarse-grained fields ρk, θk and wk, or combinations thereof. The functions
{ξ¯β}β=1,...,3 span a dual subspace of that spanned by the eigenfunctions and for any linear combination of the
hydrodynamic modes
g(c) =
3
∑
β=1
aβξβ(c), (51)
Page 7
7the coefficients aβ are given by
aβ = 〈ξ¯β |g〉 =
∫
dcχ−1H (c)ξ¯β(c)g(c). (52)
In particular, the projection of the distribution function δχk on the subspace spanned by the functions ξβ is given by
the coefficients
{〈ξ¯β |δχk〉} =
{
1
1 + z
ρk −
z
2(1 + z)
θk,
1
1 + z
ρk +
1
2(1 + z)
θk,wk
}
. (53)
Notably, these coefficients are simply linear combinations of the hydrodynamic fields linearized around the homoge-
neous decay state.
B. Projection of the Linearized Boltzmann Equation on the Hydrodynamic Subspace
In this section, we study the linearized Boltzmann equation on the hydrodynamic subspace. Let us define the
projectors
Ph(c) =
3
∑
β=1
〈ξ¯β |h〉ξβ(c), (54)
and
P⊥ = 1− P. (55)
so that any function can be decomposed as
h(c) = Ph(c) + P⊥h(c). (56)
In the definition (54) we are considering the functions (43)-(45) and (48)-(50) defined above.
Let us now consider the function δχk. If we apply the projectors P and P⊥ to equation (25), we obtain the following
relations
[
∂
∂τ
− P (Λ− ilHk · c)P
]
Pδχk = −PilHk · cP⊥δχk + PΛP⊥δχk, (57)
[
∂
∂τ
− P⊥(Λ− ilHk · c)P⊥
]
P⊥δχk = −P⊥ilHk · cPδχk, (58)
where we have used that
P⊥ΛP = 0 , (59)
which is obtained straightforwardly since ξβ are right-eigenfunctions of Λ. We note however that the {ξ¯β}β=1,...,3 are
not left-eigenfunctions of Λ, so that PΛP⊥ 6= 0. This also means that P and Λ do not commute.
Equations (57) and (58) for the functions Pδχk and P⊥δχk are coupled. Nevertheless, we shall see that, under
certain conditions, we can decouple the equation for Pδχk in the long time limit. If we solve formally equation (58),
we obtain
P⊥δχk(c, τ) = G0(τ)P⊥δχk(c, 0)−
∫ τ
0
dτ ′Gτ ′(τ − τ ′)P⊥ilH(τ ′)k · cPδχk(c, τ ′), (60)
where we have introduced the operator Gτ ′(τ − τ ′) defined from
d
dτ
Gτ ′(τ − τ ′) = P⊥ [Λ(c)− ilH(τ)k · c]P⊥Gτ ′(τ − τ ′), Gτ ′(0) = 1. (61)
If the hydrodynamic eigenvalues of the operator Λ are separated enough from the rest of the spectrum, the first term
on the right hand side of (60) decays with the “non hydrodynamic” modes, faster than the second one. We can then
write
P⊥δχk(c, τ) ≈ −
∫ τ
0
dτ ′Gτ ′(τ − τ ′)P⊥ilH(τ − τ ′)k · cPδχk(c, τ − τ ′). (62)
aβ = 〈ξ¯β |g〉 =
∫
dcχ−1H (c)ξ¯β(c)g(c). (52)
In particular, the projection of the distribution function δχk on the subspace spanned by the functions ξβ is given by
the coefficients
{〈ξ¯β |δχk〉} =
{
1
1 + z
ρk −
z
2(1 + z)
θk,
1
1 + z
ρk +
1
2(1 + z)
θk,wk
}
. (53)
Notably, these coefficients are simply linear combinations of the hydrodynamic fields linearized around the homoge-
neous decay state.
B. Projection of the Linearized Boltzmann Equation on the Hydrodynamic Subspace
In this section, we study the linearized Boltzmann equation on the hydrodynamic subspace. Let us define the
projectors
Ph(c) =
3
∑
β=1
〈ξ¯β |h〉ξβ(c), (54)
and
P⊥ = 1− P. (55)
so that any function can be decomposed as
h(c) = Ph(c) + P⊥h(c). (56)
In the definition (54) we are considering the functions (43)-(45) and (48)-(50) defined above.
Let us now consider the function δχk. If we apply the projectors P and P⊥ to equation (25), we obtain the following
relations
[
∂
∂τ
− P (Λ− ilHk · c)P
]
Pδχk = −PilHk · cP⊥δχk + PΛP⊥δχk, (57)
[
∂
∂τ
− P⊥(Λ− ilHk · c)P⊥
]
P⊥δχk = −P⊥ilHk · cPδχk, (58)
where we have used that
P⊥ΛP = 0 , (59)
which is obtained straightforwardly since ξβ are right-eigenfunctions of Λ. We note however that the {ξ¯β}β=1,...,3 are
not left-eigenfunctions of Λ, so that PΛP⊥ 6= 0. This also means that P and Λ do not commute.
Equations (57) and (58) for the functions Pδχk and P⊥δχk are coupled. Nevertheless, we shall see that, under
certain conditions, we can decouple the equation for Pδχk in the long time limit. If we solve formally equation (58),
we obtain
P⊥δχk(c, τ) = G0(τ)P⊥δχk(c, 0)−
∫ τ
0
dτ ′Gτ ′(τ − τ ′)P⊥ilH(τ ′)k · cPδχk(c, τ ′), (60)
where we have introduced the operator Gτ ′(τ − τ ′) defined from
d
dτ
Gτ ′(τ − τ ′) = P⊥ [Λ(c)− ilH(τ)k · c]P⊥Gτ ′(τ − τ ′), Gτ ′(0) = 1. (61)
If the hydrodynamic eigenvalues of the operator Λ are separated enough from the rest of the spectrum, the first term
on the right hand side of (60) decays with the “non hydrodynamic” modes, faster than the second one. We can then
write
P⊥δχk(c, τ) ≈ −
∫ τ
0
dτ ′Gτ ′(τ − τ ′)P⊥ilH(τ − τ ′)k · cPδχk(c, τ − τ ′). (62)
Page 8
8and we see, by substituting equation (62) in (57), that we obtain an involved but closed equation for Pδχk. It is
worth pointing out that we have not proved scale separation, but assumed it in order to derive (62). For an explicit
discussion of the scale separation assumption in a similar but somewhat simplified context, we refer to Appendix C,
already alluded to above.
The set of hydrodynamic equations (31)-(33) have been obtained through the projection of the Boltzmann equation
onto the hydrodynamic subspace. It now appears that, in the hydrodynamic time scale, the use of Eq. (62) will allow
us to close these equations by substituting the distribution function by its decomposition in terms of the projectors,
δχk = Pδχk + P⊥δχk. This is the aim of the next section.
IV. LINEAR HYDRODYNAMIC EQUATIONS IN NAVIER-STOKES ORDER
In this section we shall use the decomposition of δχk into its hydrodynamic part, Pδχk, and non-hydrodynamic
part, P⊥δχk, to close the linear hydrodynamic equations (31)-(33). We shall do so in Navier-Stokes order, that is, in
the long time limit and in second order in the gradients (order k2).
Let us first introduce Pδχk in the linear pressure tensor and in the heat flux vector. Here the calculation is
straightforward and we obtain
Π[Pδχk] = 0, φ[Pδχk] = 0, (63)
because the functions χH(c)∆(c) and χH(c)Σ(c) are orthogonal to the subspace spanned by the hydrodynamic
eigenfunctions {ξβ(c)}β=1,...3. Turning our attention to the other functionals δζn, δζu and δζT , the calculations
become somewhat lengthy, and we show the details in the Appendix D. We obtain
δζn[Pδχk] = −4ζnρk − ζnθk, (64)
δζu[Pδχk] = −2ζnwk, (65)
δζT [Pδχk] = −4ζTρk − (3ζT + 2ζn)θk. (66)
The negative signs occurring on the right hand side of these relations account for the fact that a fluctuation with a
local enhanced density will induce an increased collision rate, hence a faster density decay. The same remark holds
for temperature or local velocity flow fluctuations.
We now have to calculate the contribution of P⊥δχk to the same functionals, to second order in k. This requires the
knowledge of P⊥δχk to first order in k since the heat flux and pressure tensor enter the balance equations (31)-(33)
through their gradients and are already weighted by a factor k. However, it should be noted that for consistency, the
decay rates should be computed to second order in the gradients (see Eqs. (31)-(33)). We shall nevertheless restrict
to first order, henceforth neglecting the various terms of order two that symmetry allows (such as ∇2n and ∇2T for
δζn and δζT , or as ∇2u for δζu). We will further comment this approximation below. To leading order, we have that
Gτ−τ ′(τ ′) ≈ eP⊥ΛP⊥τ
′
, so that we obtain from equation (62)
P⊥δχk(c, τ) ≈
≈ −
∫ τ
0
dτ ′eP⊥ΛP⊥τ
′
P⊥ilH(τ − τ ′)k · cPδχk(c, τ − τ ′)
≈ −lH(τ)
∫ τ
0
dτ ′eP⊥ΛP⊥τ
′
e−2pζnτ
′
P⊥ik · cPδχk(c, τ − τ ′), (67)
where we have used that lH(τ) ∝ e2pζnτ (Eq. (27)). We now have to relate Pδχk(τ − τ ′) to Pδχk(τ), and to be
consistent with the approximation made above, we also have to restrict to leading order in k. In doing so, Markovian
equations for the fields will be derived. From Eq. (57), we get
Pδχk(c, τ − τ ′) ≈ e−PΛPτ
′
Pδχk(c, τ) =
3
∑
β=1
e−λβτ
′〈ξ¯β |δχk(τ)〉ξβ(c). (68)
Substituting (68) in (67), we obtain an equation for P⊥δχk to first order in k
P⊥δχ
(1)
k (c, τ) = −lH(τ)
3
∑
β=1
〈ξ¯β |δχk(τ)〉
∫ τ
0
dτ ′eP⊥(Λ−2pζn−λβ)P⊥τ
′
P⊥ik · cξβ(c), (69)
worth pointing out that we have not proved scale separation, but assumed it in order to derive (62). For an explicit
discussion of the scale separation assumption in a similar but somewhat simplified context, we refer to Appendix C,
already alluded to above.
The set of hydrodynamic equations (31)-(33) have been obtained through the projection of the Boltzmann equation
onto the hydrodynamic subspace. It now appears that, in the hydrodynamic time scale, the use of Eq. (62) will allow
us to close these equations by substituting the distribution function by its decomposition in terms of the projectors,
δχk = Pδχk + P⊥δχk. This is the aim of the next section.
IV. LINEAR HYDRODYNAMIC EQUATIONS IN NAVIER-STOKES ORDER
In this section we shall use the decomposition of δχk into its hydrodynamic part, Pδχk, and non-hydrodynamic
part, P⊥δχk, to close the linear hydrodynamic equations (31)-(33). We shall do so in Navier-Stokes order, that is, in
the long time limit and in second order in the gradients (order k2).
Let us first introduce Pδχk in the linear pressure tensor and in the heat flux vector. Here the calculation is
straightforward and we obtain
Π[Pδχk] = 0, φ[Pδχk] = 0, (63)
because the functions χH(c)∆(c) and χH(c)Σ(c) are orthogonal to the subspace spanned by the hydrodynamic
eigenfunctions {ξβ(c)}β=1,...3. Turning our attention to the other functionals δζn, δζu and δζT , the calculations
become somewhat lengthy, and we show the details in the Appendix D. We obtain
δζn[Pδχk] = −4ζnρk − ζnθk, (64)
δζu[Pδχk] = −2ζnwk, (65)
δζT [Pδχk] = −4ζTρk − (3ζT + 2ζn)θk. (66)
The negative signs occurring on the right hand side of these relations account for the fact that a fluctuation with a
local enhanced density will induce an increased collision rate, hence a faster density decay. The same remark holds
for temperature or local velocity flow fluctuations.
We now have to calculate the contribution of P⊥δχk to the same functionals, to second order in k. This requires the
knowledge of P⊥δχk to first order in k since the heat flux and pressure tensor enter the balance equations (31)-(33)
through their gradients and are already weighted by a factor k. However, it should be noted that for consistency, the
decay rates should be computed to second order in the gradients (see Eqs. (31)-(33)). We shall nevertheless restrict
to first order, henceforth neglecting the various terms of order two that symmetry allows (such as ∇2n and ∇2T for
δζn and δζT , or as ∇2u for δζu). We will further comment this approximation below. To leading order, we have that
Gτ−τ ′(τ ′) ≈ eP⊥ΛP⊥τ
′
, so that we obtain from equation (62)
P⊥δχk(c, τ) ≈
≈ −
∫ τ
0
dτ ′eP⊥ΛP⊥τ
′
P⊥ilH(τ − τ ′)k · cPδχk(c, τ − τ ′)
≈ −lH(τ)
∫ τ
0
dτ ′eP⊥ΛP⊥τ
′
e−2pζnτ
′
P⊥ik · cPδχk(c, τ − τ ′), (67)
where we have used that lH(τ) ∝ e2pζnτ (Eq. (27)). We now have to relate Pδχk(τ − τ ′) to Pδχk(τ), and to be
consistent with the approximation made above, we also have to restrict to leading order in k. In doing so, Markovian
equations for the fields will be derived. From Eq. (57), we get
Pδχk(c, τ − τ ′) ≈ e−PΛPτ
′
Pδχk(c, τ) =
3
∑
β=1
e−λβτ
′〈ξ¯β |δχk(τ)〉ξβ(c). (68)
Substituting (68) in (67), we obtain an equation for P⊥δχk to first order in k
P⊥δχ
(1)
k (c, τ) = −lH(τ)
3
∑
β=1
〈ξ¯β |δχk(τ)〉
∫ τ
0
dτ ′eP⊥(Λ−2pζn−λβ)P⊥τ
′
P⊥ik · cξβ(c), (69)
Page 9
9where P⊥δχk = P⊥δχ
(1)
k + O(k2). The pressure tensor and the heat flux up to first order in the gradients of the
fields are now obtained by substituting equation (69) into equations (34) and (35). Taking into account the symmetry
properties of the system, the resulting expressions can be written in the form
Πij [P⊥δχ
(1)
k ] = −ilH(τ)η˜(τ)
[
kjwi,k + kiwj,k +
2
d
k ·wkδij
]
, (70)
φ[P⊥δχ(1)k ] = −ilH(τ)k [κ˜(τ)θk + µ˜(τ)ρk] . (71)
Equation (70) is the expected Navier-Stokes expression for the pressure tensor, involving the shear viscosity coefficient
η˜, but equation (71) contains, besides the usual Fourier law characterized by the heat conductivity κ˜, an additional
contribution proportional to the density gradient and with an associated transport coefficient µ˜. This latter term is
analogous to the one appearing in granular gases [3, 18].
The expression of the (time dependent) transport coefficients are
η˜(τ) =
∫
dc∆xy(c)F3,xy(c, τ) =
1
d2 + d− 2
d
∑
i,j
∫
dc∆ij(c)F3,ij(c, τ), (72)
µ˜(τ) =
1
d(1 + z)
∫
dcΣ(c) [F1(c, τ) + F2(c, τ)] , (73)
κ˜(τ) =
1
2d(1 + z)
∫
dcΣ(c) [−zF1(c, τ) + F2(c, τ)] , (74)
where we have introduced the functions
F3,ij(c, τ) =
∫ τ
0
dτ ′eP⊥(Λ−2pζn−pζT )τ
′
P⊥ciξ3,j(c), (75)
F1(c, τ) =
∫ τ
0
dτ ′eP⊥(Λ−2pζn)τ
′
P⊥cξ1(c), (76)
F2(c, τ) =
∫ τ
0
dτ ′eP⊥(Λ+pζT )τ
′
P⊥cξ2(c), (77)
and in the second equality of equation (72), we have summed over all the i, j, taking into account the symmetry of
the linearized Boltzmann operator.
Similarly, we calculate the deviations of the decay rates to first order in the gradients of the fields by substituting
equation (69) into equations (38), (39) and (40). Taking into account the symmetry properties, we arrive at
δζn[P⊥δχ
(1)
k ] = 0, (78)
δζu[P⊥δχ(1)k ] = ilH(τ)k [ζu,ρ(τ)ρk + ζu,θ(τ)θk] , (79)
δζT [P⊥δχ
(1)
k ] = 0. (80)
The expression for the coefficients are
ζu,ρ(τ) =
γβ
d(1 + z)
∫
dc1
∫
dc2χH(c1)c12(c1 + c2) · [F1(c2, τ) + F2(c2, τ)] , (81)
ζu,θ(τ) =
γβ
2d(1 + z)
∫
dc1
∫
dc2χH(c1)c12(c1 + c2) · [−zF1(c2, τ) + F2(c2, τ)] , (82)
with β = pi(d−1)/2/Γ[(d + 1)/2] the d-dimensional solid angle.
At this point, it is important to note that the transport coefficients defined in equations (72)-(74) and (81)-(82)
are time-dependent, and this dependence is governed by the Fi functions. The (exponential) integrands appearing
in the definitions of the Fi decay with the non-hydrodynamic (kinetic) modes, as a consequence of the action of
the projector P⊥. From our assumption of hydrodynamic versus kinetic scale separation, all kinetic eigenvalues are
smaller than the smallest hydrodynamic eigenvalue of Λ, which is λ2 = −pζT − 2pζn. This ensures the convergence of
the integrals (75)-(77) for τ →∞. In order for the transport coefficients to reach their τ →∞ limit faster than any of
(1)
k + O(k2). The pressure tensor and the heat flux up to first order in the gradients of the
fields are now obtained by substituting equation (69) into equations (34) and (35). Taking into account the symmetry
properties of the system, the resulting expressions can be written in the form
Πij [P⊥δχ
(1)
k ] = −ilH(τ)η˜(τ)
[
kjwi,k + kiwj,k +
2
d
k ·wkδij
]
, (70)
φ[P⊥δχ(1)k ] = −ilH(τ)k [κ˜(τ)θk + µ˜(τ)ρk] . (71)
Equation (70) is the expected Navier-Stokes expression for the pressure tensor, involving the shear viscosity coefficient
η˜, but equation (71) contains, besides the usual Fourier law characterized by the heat conductivity κ˜, an additional
contribution proportional to the density gradient and with an associated transport coefficient µ˜. This latter term is
analogous to the one appearing in granular gases [3, 18].
The expression of the (time dependent) transport coefficients are
η˜(τ) =
∫
dc∆xy(c)F3,xy(c, τ) =
1
d2 + d− 2
d
∑
i,j
∫
dc∆ij(c)F3,ij(c, τ), (72)
µ˜(τ) =
1
d(1 + z)
∫
dcΣ(c) [F1(c, τ) + F2(c, τ)] , (73)
κ˜(τ) =
1
2d(1 + z)
∫
dcΣ(c) [−zF1(c, τ) + F2(c, τ)] , (74)
where we have introduced the functions
F3,ij(c, τ) =
∫ τ
0
dτ ′eP⊥(Λ−2pζn−pζT )τ
′
P⊥ciξ3,j(c), (75)
F1(c, τ) =
∫ τ
0
dτ ′eP⊥(Λ−2pζn)τ
′
P⊥cξ1(c), (76)
F2(c, τ) =
∫ τ
0
dτ ′eP⊥(Λ+pζT )τ
′
P⊥cξ2(c), (77)
and in the second equality of equation (72), we have summed over all the i, j, taking into account the symmetry of
the linearized Boltzmann operator.
Similarly, we calculate the deviations of the decay rates to first order in the gradients of the fields by substituting
equation (69) into equations (38), (39) and (40). Taking into account the symmetry properties, we arrive at
δζn[P⊥δχ
(1)
k ] = 0, (78)
δζu[P⊥δχ(1)k ] = ilH(τ)k [ζu,ρ(τ)ρk + ζu,θ(τ)θk] , (79)
δζT [P⊥δχ
(1)
k ] = 0. (80)
The expression for the coefficients are
ζu,ρ(τ) =
γβ
d(1 + z)
∫
dc1
∫
dc2χH(c1)c12(c1 + c2) · [F1(c2, τ) + F2(c2, τ)] , (81)
ζu,θ(τ) =
γβ
2d(1 + z)
∫
dc1
∫
dc2χH(c1)c12(c1 + c2) · [−zF1(c2, τ) + F2(c2, τ)] , (82)
with β = pi(d−1)/2/Γ[(d + 1)/2] the d-dimensional solid angle.
At this point, it is important to note that the transport coefficients defined in equations (72)-(74) and (81)-(82)
are time-dependent, and this dependence is governed by the Fi functions. The (exponential) integrands appearing
in the definitions of the Fi decay with the non-hydrodynamic (kinetic) modes, as a consequence of the action of
the projector P⊥. From our assumption of hydrodynamic versus kinetic scale separation, all kinetic eigenvalues are
smaller than the smallest hydrodynamic eigenvalue of Λ, which is λ2 = −pζT − 2pζn. This ensures the convergence of
the integrals (75)-(77) for τ →∞. In order for the transport coefficients to reach their τ →∞ limit faster than any of
Page 10
10
the hydrodynamic time scales, we need moreover the more stringent condition that the fastest kinetic mode is at least
separated by a pζT gap from λ2 = −pζT − 2pζn: under this condition, the time dependence of the exponential term
in the integral giving F2 is fast enough so that the transport coefficients, that depend on the Fi functions through
(72)-(74), can be considered as constants on the hydrodynamic time scale. With this proviso in mind, it is possible
to set τ → ∞ in the integrals (75)-(77) and the time-independent transport coefficients obtained in this section are
then equivalent to those calculated in reference [12] by the Chapman-Enskog method. We recall in Appendix A their
expressions in the first order Sonine approximation.
Finally, if we substitute the expressions derived above for the fluxes, equations (70)-(71), and the decay rates,
equations (78)-(80), and we take into account that in the hydrodynamic time scale we can substitute all the coefficients
by their τ →∞ limit, we obtain the following closed equations for the linear deviation of the hydrodynamic fields
(
∂
∂τ
+ 2pζn
)
ρk + ilH(τ)kwk|| + pζnθk = 0, (83)
[
∂
∂τ
− pζT + l2H(τ)η˜k2
]
wk⊥ = 0, (84)
[
∂
∂τ
− pζT +
2(d− 1)
d
l2H(τ)η˜k
2
]
wk||
+
i
2
lH(τ)k [(1 − 2pζu,ρ)ρk + (1 − 2pζu,θ)θk] = 0, (85)
[
∂
∂τ
+ pζT +
2
d
l2H(τ)κ˜k
2
]
θk
+
[
2pζT +
2
d
l2H(τ)µ˜k
2
]
ρk + i
2
d
lH(τ)kwk|| = 0, (86)
where wk|| and wk⊥ are the longitudinal and transversal parts of the velocity vector defined by
wk|| = wk · kˆ, wk⊥ = wk − wk||kˆ, (87)
and kˆ is the unit vector along the direction given by k.
Equation (84) for the shear mode is decoupled from the other equations and can be readily integrated. If we
introduce a non-dimensional wave number k˜ = lH(0)k, scaled by the mean free path at the time origin, we obtain the
explicit solution
wk⊥(τ) = exp
[
pζT τ −
η˜k˜2
4pζn
(
e4pζnτ − 1
)
]
wk⊥(0). (88)
Interestingly, depending on k˜, the perturbation may initially increase if pζT − η˜k˜2 > 0. For long times however, the
exponential e4pζnτ always dominates the linear term pζT τ and the perturbation decays.
The other three fields, namely, the density ρk, temperature θk, and the longitudinal velocity wk||, obey the system
of coupled linear equations
∂
∂τ
ρk
wk||
θk
= M(τ) ·
ρk
wk||
θk
, (89)
where the time-dependent matrix is
M(τ) =
−2pζn −ilH(τ)k −pζn
− i2 lH(τ)k(1 − 2pζu,ρ) pζT −
2(d−1)
d l
2
H(τ)η˜k
2 − i2 lH(τ)k(1 − 2pζu,θ)
−2pζT − 2d l2H(τ)µ˜k2 −i 2d lH(τ)k −pζT − 2d l2H(τ)κ˜k2
.
(90)
We note here that this matrix differs from Eq. (59) of Ref. [12], where the analysis amounts to overlooking the time
dependence of the mean free path, so that all entries of the hydrodynamic matrix exhibit the same time dependence.
The different time dependences present in Eq. (90) render the stability analysis more difficult. For long times however,
the eigenvalues of the matrix M(τ) are always negative and the perturbations a priori decay. A caveat is nevertheless
the hydrodynamic time scales, we need moreover the more stringent condition that the fastest kinetic mode is at least
separated by a pζT gap from λ2 = −pζT − 2pζn: under this condition, the time dependence of the exponential term
in the integral giving F2 is fast enough so that the transport coefficients, that depend on the Fi functions through
(72)-(74), can be considered as constants on the hydrodynamic time scale. With this proviso in mind, it is possible
to set τ → ∞ in the integrals (75)-(77) and the time-independent transport coefficients obtained in this section are
then equivalent to those calculated in reference [12] by the Chapman-Enskog method. We recall in Appendix A their
expressions in the first order Sonine approximation.
Finally, if we substitute the expressions derived above for the fluxes, equations (70)-(71), and the decay rates,
equations (78)-(80), and we take into account that in the hydrodynamic time scale we can substitute all the coefficients
by their τ →∞ limit, we obtain the following closed equations for the linear deviation of the hydrodynamic fields
(
∂
∂τ
+ 2pζn
)
ρk + ilH(τ)kwk|| + pζnθk = 0, (83)
[
∂
∂τ
− pζT + l2H(τ)η˜k2
]
wk⊥ = 0, (84)
[
∂
∂τ
− pζT +
2(d− 1)
d
l2H(τ)η˜k
2
]
wk||
+
i
2
lH(τ)k [(1 − 2pζu,ρ)ρk + (1 − 2pζu,θ)θk] = 0, (85)
[
∂
∂τ
+ pζT +
2
d
l2H(τ)κ˜k
2
]
θk
+
[
2pζT +
2
d
l2H(τ)µ˜k
2
]
ρk + i
2
d
lH(τ)kwk|| = 0, (86)
where wk|| and wk⊥ are the longitudinal and transversal parts of the velocity vector defined by
wk|| = wk · kˆ, wk⊥ = wk − wk||kˆ, (87)
and kˆ is the unit vector along the direction given by k.
Equation (84) for the shear mode is decoupled from the other equations and can be readily integrated. If we
introduce a non-dimensional wave number k˜ = lH(0)k, scaled by the mean free path at the time origin, we obtain the
explicit solution
wk⊥(τ) = exp
[
pζT τ −
η˜k˜2
4pζn
(
e4pζnτ − 1
)
]
wk⊥(0). (88)
Interestingly, depending on k˜, the perturbation may initially increase if pζT − η˜k˜2 > 0. For long times however, the
exponential e4pζnτ always dominates the linear term pζT τ and the perturbation decays.
The other three fields, namely, the density ρk, temperature θk, and the longitudinal velocity wk||, obey the system
of coupled linear equations
∂
∂τ
ρk
wk||
θk
= M(τ) ·
ρk
wk||
θk
, (89)
where the time-dependent matrix is
M(τ) =
−2pζn −ilH(τ)k −pζn
− i2 lH(τ)k(1 − 2pζu,ρ) pζT −
2(d−1)
d l
2
H(τ)η˜k
2 − i2 lH(τ)k(1 − 2pζu,θ)
−2pζT − 2d l2H(τ)µ˜k2 −i 2d lH(τ)k −pζT − 2d l2H(τ)κ˜k2
.
(90)
We note here that this matrix differs from Eq. (59) of Ref. [12], where the analysis amounts to overlooking the time
dependence of the mean free path, so that all entries of the hydrodynamic matrix exhibit the same time dependence.
The different time dependences present in Eq. (90) render the stability analysis more difficult. For long times however,
the eigenvalues of the matrix M(τ) are always negative and the perturbations a priori decay. A caveat is nevertheless
Page 11
11
0
0,01
0,02
0,03
0,04
0,05
wk
m
⊥
0 5 10 15 20
τ
A)
0
0,025
0,05
nσ
2
0 5 10 15 20
τ
-1,5
-1
-0,5
0
ln
(δu
k m
⊥(τ
)/δ
u
k m
⊥(0
))
0 100 200 300 400
(n(0)/n(τ))2
B)
Figure 1: A) Time evolution of a linear perturbation of the transversal velocity for a system with annihilation probability
p = 0.1. The solid lines are the Molecular Dynamics results and the dashed lines are the theoretical predictions (no adjustable
parameter). Note how the theoretical predictions correctly account for the increase of the perturbation at short times. The
inset shows the evolution of density with rescaled time τ , where σ is the discs’ radius. For τ = 20, nσ2 ≃ 9.10−4. B)
Evolution of δuk⊥(τ )/δuk⊥(0) as a function of n2(0)/n2(τ ). The dashed line is the exponential decay predicted by Eq. (92).
n2(0)/n2(τ ) = 400 corresponds to τ ≃ 15.
in order. It is worth pointing out that equations (88) and (89) break down at long times. Once the function lH(τ)k
exceeds unity indeed, the expansion in the gradients we have performed is no longer valid, and it would be necessary
to include terms of higher order in k. In addition, if the perturbation initially increases sufficiently to leave the linear
regime, our description breaks down and it becomes necessary to consider non linear terms.
Although quite involved, the evolution equations (89) can be numerically integrated, using for instance the transport
coefficients computed in the Sonine approximation in Appendix A. This is what we do in the next section, in order
to compare the theoretical predictions with Molecular Dynamics simulations.
V. MOLECULAR DYNAMICS SIMULATIONS
To put our theoretical predictions to the test, we have performed Molecular Dynamics (MD) simulations of a system
of N smooth hard disks (d = 2) which undergo ballistic flights punctuated by collisional events : at each collision, the
discs annihilate with probability p ; otherwise, they collide elastically. The particles are localized in a square box of
size L with periodic boundary conditions. An 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 are N(0) = 105,
reduced density n(0)σ2 = 0.05 where σ is the discs’ radius and 0 < p ≤ 1. The initial conditions we have considered
correspond to small amplitude perturbations around the homogeneous decay state, to enforce the validity of the
linearized hydrodynamic equations (83)-(86).
A. Perturbation of the transversal velocity
Since Eq. (84) for the shear mode is decoupled from the other equations, one of the simplest macroscopic pertur-
bation one can think of consists in an initial harmonic perturbation of the transversal component of the velocity field,
whose evolution is given by Eq. (88). We shall consider a small perturbation in real space of the form
ux(r, 0) = A sin(kmy), (91)
with A = 10−1vH(0) and km = 2pi/L, where L is the linear size of the system. The reason for choosing the smallest
possible value of k compatible with the boundary conditions is twofold. First, the corresponding mode is the most
unstable at short times (see Eq. (88)). For the parameters of the simulations, we indeed probe the region where
pζT > η˜k˜2, so that the hydrodynamic equation (88) predicts an initial increase of the perturbation. Second, the low
k regime is that where our large scale predictions are most likely to be relevant.
Figure 1 displays the evolution of wkm⊥ as a function of the number of collisions per particle τ for a system with
p = 0.1, averaging data over 50 different trajectories. The solid line is the simulation result and the dashed line is
the theoretical prediction, Eq. (88), where the shear viscosity and the decay rates have been computed using the
0
0,01
0,02
0,03
0,04
0,05
wk
m
⊥
0 5 10 15 20
τ
A)
0
0,025
0,05
nσ
2
0 5 10 15 20
τ
-1,5
-1
-0,5
0
ln
(δu
k m
⊥(τ
)/δ
u
k m
⊥(0
))
0 100 200 300 400
(n(0)/n(τ))2
B)
Figure 1: A) Time evolution of a linear perturbation of the transversal velocity for a system with annihilation probability
p = 0.1. The solid lines are the Molecular Dynamics results and the dashed lines are the theoretical predictions (no adjustable
parameter). Note how the theoretical predictions correctly account for the increase of the perturbation at short times. The
inset shows the evolution of density with rescaled time τ , where σ is the discs’ radius. For τ = 20, nσ2 ≃ 9.10−4. B)
Evolution of δuk⊥(τ )/δuk⊥(0) as a function of n2(0)/n2(τ ). The dashed line is the exponential decay predicted by Eq. (92).
n2(0)/n2(τ ) = 400 corresponds to τ ≃ 15.
in order. It is worth pointing out that equations (88) and (89) break down at long times. Once the function lH(τ)k
exceeds unity indeed, the expansion in the gradients we have performed is no longer valid, and it would be necessary
to include terms of higher order in k. In addition, if the perturbation initially increases sufficiently to leave the linear
regime, our description breaks down and it becomes necessary to consider non linear terms.
Although quite involved, the evolution equations (89) can be numerically integrated, using for instance the transport
coefficients computed in the Sonine approximation in Appendix A. This is what we do in the next section, in order
to compare the theoretical predictions with Molecular Dynamics simulations.
V. MOLECULAR DYNAMICS SIMULATIONS
To put our theoretical predictions to the test, we have performed Molecular Dynamics (MD) simulations of a system
of N smooth hard disks (d = 2) which undergo ballistic flights punctuated by collisional events : at each collision, the
discs annihilate with probability p ; otherwise, they collide elastically. The particles are localized in a square box of
size L with periodic boundary conditions. An 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 are N(0) = 105,
reduced density n(0)σ2 = 0.05 where σ is the discs’ radius and 0 < p ≤ 1. The initial conditions we have considered
correspond to small amplitude perturbations around the homogeneous decay state, to enforce the validity of the
linearized hydrodynamic equations (83)-(86).
A. Perturbation of the transversal velocity
Since Eq. (84) for the shear mode is decoupled from the other equations, one of the simplest macroscopic pertur-
bation one can think of consists in an initial harmonic perturbation of the transversal component of the velocity field,
whose evolution is given by Eq. (88). We shall consider a small perturbation in real space of the form
ux(r, 0) = A sin(kmy), (91)
with A = 10−1vH(0) and km = 2pi/L, where L is the linear size of the system. The reason for choosing the smallest
possible value of k compatible with the boundary conditions is twofold. First, the corresponding mode is the most
unstable at short times (see Eq. (88)). For the parameters of the simulations, we indeed probe the region where
pζT > η˜k˜2, so that the hydrodynamic equation (88) predicts an initial increase of the perturbation. Second, the low
k regime is that where our large scale predictions are most likely to be relevant.
Figure 1 displays the evolution of wkm⊥ as a function of the number of collisions per particle τ for a system with
p = 0.1, averaging data over 50 different trajectories. The solid line is the simulation result and the dashed line is
the theoretical prediction, Eq. (88), where the shear viscosity and the decay rates have been computed using the
Page 12
12
-0,3
-0,25
-0,2
-0,15
-0,1
-0,05
0
ln
(δu
k m
⊥(τ
)/δ
u
k m
⊥(0
))
0 100 200 300 400
(n(0)/n(τ))2
B)
Figure 2: A) Same as in Fig. 1 but for a system with p = 0.5, which –loosely speaking– may therefore be considered as being
more “distant” to equilibrium. The inset shows the evolution of density with rescaled time τ . For τ = 4, the rescaled density
is nσ2 ≃ 9.10−4 (while n(0)σ2 = 0.05). B) Evolution of δuk⊥(τ )/δuk⊥(0) as a function of n2(0)/n2(τ ). n2(0)/n2(τ ) = 400
corresponds to τ ≃ 3.
0,6
0,7
0,8
0,9
1
∼ η/
(ζ n
~ η e
)
0 0,2 0,4 0,6 0,8 1
p
Figure 3: Symbols: ratio η˜/ζn (normalized by its p→ 0 value) extracted from the exponential fit of δuk⊥(τ ) to Eq. (92). The
solid line is the theoretical prediction in the first Sonine approximation.
standard tools of kinetic theory (here, the first Sonine approximation [12], see Appendix A). An excellent agreement
is obtained without any adjustable parameter, including the predicted increase of the perturbation at short times, and
as also observed for a larger annihilation probability p = 0.5 (see Fig. 2, in which the MD data are obtained by an
average over 150 runs).
Recalling that nH(τ)/nH(0) = exp(−2pζnτ) and that vH(τ) = vH(0) exp(−pζT τ), it proves also convenient to
consider the actual velocity field δu(τ) = vH(τ)w(τ) instead of its dimensionless counterpart w, since the prediction
(88) then takes the form :
δuk⊥(τ) = exp
{
− η˜k˜
2
4pζn
[
n2H(0)
n2H(τ)
− 1
]
}
δuk⊥(0). (92)
The plot of δuk⊥(τ)/δuk⊥(0) as a function of (nH(0)/nH(τ))2 allows then by simple exponential fitting to extract
η˜/ζn (recall that k˜ = lH(0)k is known). Figure 3 compares the ratios η˜/ζn extracted from such fits for various values
of p with the theoretical prediction. We have plotted the value of the viscosity normalized by its elastic value, η˜e. For
p = 0.1 the agreement is quite good but for p = 1 we obtain discrepancies of the order of 15%. Such deviations could
be due to the limitation for high dissipation of the first Sonine approximation that has been used to compute the
numerical values of the various quantities involved in the description (in particular, the deviations of the homogeneous
decay state velocity distribution from its Gaussian form might be relevant). Additionally, the shear viscosity could
suffer from finite size effects. For a related discussion in the realm of granular gases, where both effects alluded to
are at work, see [20–22]. Finally, neglecting the k2 contribution to δζu might not be innocuous. From symmetry
-0,3
-0,25
-0,2
-0,15
-0,1
-0,05
0
ln
(δu
k m
⊥(τ
)/δ
u
k m
⊥(0
))
0 100 200 300 400
(n(0)/n(τ))2
B)
Figure 2: A) Same as in Fig. 1 but for a system with p = 0.5, which –loosely speaking– may therefore be considered as being
more “distant” to equilibrium. The inset shows the evolution of density with rescaled time τ . For τ = 4, the rescaled density
is nσ2 ≃ 9.10−4 (while n(0)σ2 = 0.05). B) Evolution of δuk⊥(τ )/δuk⊥(0) as a function of n2(0)/n2(τ ). n2(0)/n2(τ ) = 400
corresponds to τ ≃ 3.
0,6
0,7
0,8
0,9
1
∼ η/
(ζ n
~ η e
)
0 0,2 0,4 0,6 0,8 1
p
Figure 3: Symbols: ratio η˜/ζn (normalized by its p→ 0 value) extracted from the exponential fit of δuk⊥(τ ) to Eq. (92). The
solid line is the theoretical prediction in the first Sonine approximation.
standard tools of kinetic theory (here, the first Sonine approximation [12], see Appendix A). An excellent agreement
is obtained without any adjustable parameter, including the predicted increase of the perturbation at short times, and
as also observed for a larger annihilation probability p = 0.5 (see Fig. 2, in which the MD data are obtained by an
average over 150 runs).
Recalling that nH(τ)/nH(0) = exp(−2pζnτ) and that vH(τ) = vH(0) exp(−pζT τ), it proves also convenient to
consider the actual velocity field δu(τ) = vH(τ)w(τ) instead of its dimensionless counterpart w, since the prediction
(88) then takes the form :
δuk⊥(τ) = exp
{
− η˜k˜
2
4pζn
[
n2H(0)
n2H(τ)
− 1
]
}
δuk⊥(0). (92)
The plot of δuk⊥(τ)/δuk⊥(0) as a function of (nH(0)/nH(τ))2 allows then by simple exponential fitting to extract
η˜/ζn (recall that k˜ = lH(0)k is known). Figure 3 compares the ratios η˜/ζn extracted from such fits for various values
of p with the theoretical prediction. We have plotted the value of the viscosity normalized by its elastic value, η˜e. For
p = 0.1 the agreement is quite good but for p = 1 we obtain discrepancies of the order of 15%. Such deviations could
be due to the limitation for high dissipation of the first Sonine approximation that has been used to compute the
numerical values of the various quantities involved in the description (in particular, the deviations of the homogeneous
decay state velocity distribution from its Gaussian form might be relevant). Additionally, the shear viscosity could
suffer from finite size effects. For a related discussion in the realm of granular gases, where both effects alluded to
are at work, see [20–22]. Finally, neglecting the k2 contribution to δζu might not be innocuous. From symmetry
Page 13
13
-0,02
0
0,02
0,04
0,06
wk
m
//
0 5 10 15 20
τ
0 5 10 15 20
τ
0
1
2
3
4
lHk
Figure 4: Time evolution of the longitudinal velocity for a system with p = 0.1 as a function of τ . The solid line shows the
simulation results and the dashed line is the numerical solution of (89). For τ = 15 the rescaled density is nσ2 ≃ 2.5.10−3 . The
inset shows the increase of mean-free path with time and that for τ > 10, lHk is no longer a small quantity.
-0,06
-0,03
0
0,03
ek
m
0 5 10 15 20
0 5 10 15 20
τ
-0,03
-0,02
-0,01
0
0,01
ρk
m
Figure 5: Time evolution of ekm and ρkm as a function of τ for a system with p = 0.1. The solid lines show Molecular Dynamics
results while the dashed lines are for the numerical solution of (89). For τ = 15 the rescaled density is nσ2 ≃ 2.5.10−3 .
considerations, such a term must be of the form k2lH(τ)2wk, so that the equation (84) for the transversal velocity
has the same form as the one we considered, but with a “shifted” shear viscosity. It is worth pointing out here that in
the corresponding equation (88), putative order k2 corrections to δζn and δζT play no role (see Eq. (84)) : the decay
rates ζn and ζT appearing in (88) are fingerprints of the τ dependence of lH and of the rescaling procedure leading
to wk from the actual velocity flow. Those two decay rates are therefore properties of the homogeneous solution and
do not suffer any finite k correction.
B. Perturbation of the longitudinal velocity
In order to investigate further the validity of the hydrodynamic equations, we consider a perturbation of the
longitudinal velocity
ux(r, 0) = A sin(kmx), (93)
where A = 10−1vH(0) and km = 2pi/L. Since the hydrodynamic matrix M depends on time and dMdt does not
commute with M , we could not solve analytically the set of equations (89), and we turned to a numerical integration,
using the transport coefficients computed in Appendix A.
In Fig. 4, we have plotted the time evolution of the rescaled longitudinal velocity field, wkm‖, as a function of the
internal time clock τ (number of collisions per particle) for a system with p = 0.1. The results have been averaged over
16 trajectories. It can be seen that the theoretical framework is able to account for the non trivial time dependence
of the perturbation dynamics. Moreover, as the equations for θkm and ρkm are coupled with the equation for wkm‖, a
-0,02
0
0,02
0,04
0,06
wk
m
//
0 5 10 15 20
τ
0 5 10 15 20
τ
0
1
2
3
4
lHk
Figure 4: Time evolution of the longitudinal velocity for a system with p = 0.1 as a function of τ . The solid line shows the
simulation results and the dashed line is the numerical solution of (89). For τ = 15 the rescaled density is nσ2 ≃ 2.5.10−3 . The
inset shows the increase of mean-free path with time and that for τ > 10, lHk is no longer a small quantity.
-0,06
-0,03
0
0,03
ek
m
0 5 10 15 20
0 5 10 15 20
τ
-0,03
-0,02
-0,01
0
0,01
ρk
m
Figure 5: Time evolution of ekm and ρkm as a function of τ for a system with p = 0.1. The solid lines show Molecular Dynamics
results while the dashed lines are for the numerical solution of (89). For τ = 15 the rescaled density is nσ2 ≃ 2.5.10−3 .
considerations, such a term must be of the form k2lH(τ)2wk, so that the equation (84) for the transversal velocity
has the same form as the one we considered, but with a “shifted” shear viscosity. It is worth pointing out here that in
the corresponding equation (88), putative order k2 corrections to δζn and δζT play no role (see Eq. (84)) : the decay
rates ζn and ζT appearing in (88) are fingerprints of the τ dependence of lH and of the rescaling procedure leading
to wk from the actual velocity flow. Those two decay rates are therefore properties of the homogeneous solution and
do not suffer any finite k correction.
B. Perturbation of the longitudinal velocity
In order to investigate further the validity of the hydrodynamic equations, we consider a perturbation of the
longitudinal velocity
ux(r, 0) = A sin(kmx), (93)
where A = 10−1vH(0) and km = 2pi/L. Since the hydrodynamic matrix M depends on time and dMdt does not
commute with M , we could not solve analytically the set of equations (89), and we turned to a numerical integration,
using the transport coefficients computed in Appendix A.
In Fig. 4, we have plotted the time evolution of the rescaled longitudinal velocity field, wkm‖, as a function of the
internal time clock τ (number of collisions per particle) for a system with p = 0.1. The results have been averaged over
16 trajectories. It can be seen that the theoretical framework is able to account for the non trivial time dependence
of the perturbation dynamics. Moreover, as the equations for θkm and ρkm are coupled with the equation for wkm‖, a
Page 14
14
0
0,02
0,04
0,06
wk
m
//
0 1 2 3 4
τ
0 1 2 3 4
τ
0
1
2
3
4
lHk
Figure 6: Time evolution of the longitudinal velocity for a system with p = 0.5 as a function of τ . The solid line shows the
simulation results and the dashed line is for the numerical solution. We have also plotted with a dotted line the numerical
solution considering the elastic values of κ˜ and µ˜ (i.e. their limit when p → 0+). The inset shows the increase of mean-free
path with time.
-0,06
-0,04
-0,02
0
ek
m
0 1 2 3 4
0 1 2 3 4
τ
-0,04
-0,02
0
ρk
m
Figure 7: Time evolution of ekm and ρkm as a function of τ for a system with p = 0.5. The different lines have the same
meaning as in Fig. 6.
perturbation such as (93) induces a response of the above two other fields (at variance with the transversal velocity,
whose dynamics is decoupled from the other three modes, at least at the linear level of description adopted here). In
Fig. 5, we have plotted the energy density ekm = θkm + ρkm and ρkm as a function of τ . The agreement with theory
is very good both qualitatively and quantitatively, with an evolution that is well predicted till τ ≃ 15 ; for long times
the simulation data become somewhat noisy (less statistics can be achieved due to the smaller number of particles
left in the system. The number of particles at τ = 15 is N ≃ 5000).
In Fig.s 6 and 7, results for a system with p = 0.5 are shown. Data are averaged over 64 runs. In this case,
the agreement is still good qualitatively, with similar shapes for the theoretical and numerical curves, but some
discrepancies are observed. The most significant deviation from the theoretical prediction occurs for the density,
ρkm , where the value of the minimum and also its position are not predicted accurately. However, the hydrodynamic
framework still captures correctly the trends of the complex dynamics of the perturbations. The above discrepancies
could be ascribable to the failure of the first Sonine approximation for high dissipation (i.e. “high” p) or to finite size
effects. We mention here that within the first Sonine approximation, the dimensionless coefficients µ˜ and κ˜ exhibit a
divergent behavior in the vicinity of p = 0.8 [12] which is presumably unphysical, and is an indication of the limitation
of the method. For completeness and to assess the robustness of the predictions with respect to a modification of
the numerical values of the key parameters, we have also reported in Fig.s 6 and 7 the predictions obtained when
the transport coefficients take their elastic hard disc value (e.g. µ˜(p = 0.5) ≃ 0.505 while µ˜(p = 0) = 0, as required
by Fourier’s law ). The important point is that the time evolutions are not significantly affected, and that the main
features remain the same. Additionally, as mentioned after Eq. (66) and (92), a possible source of inaccuracy lies
in the truncation of decay rates to their first order (k1) in the gradients. While the corresponding terms have been
shown to be small for inelastic hard spheres [23] (with the notable simplification there that the velocity decay rate
0
0,02
0,04
0,06
wk
m
//
0 1 2 3 4
τ
0 1 2 3 4
τ
0
1
2
3
4
lHk
Figure 6: Time evolution of the longitudinal velocity for a system with p = 0.5 as a function of τ . The solid line shows the
simulation results and the dashed line is for the numerical solution. We have also plotted with a dotted line the numerical
solution considering the elastic values of κ˜ and µ˜ (i.e. their limit when p → 0+). The inset shows the increase of mean-free
path with time.
-0,06
-0,04
-0,02
0
ek
m
0 1 2 3 4
0 1 2 3 4
τ
-0,04
-0,02
0
ρk
m
Figure 7: Time evolution of ekm and ρkm as a function of τ for a system with p = 0.5. The different lines have the same
meaning as in Fig. 6.
perturbation such as (93) induces a response of the above two other fields (at variance with the transversal velocity,
whose dynamics is decoupled from the other three modes, at least at the linear level of description adopted here). In
Fig. 5, we have plotted the energy density ekm = θkm + ρkm and ρkm as a function of τ . The agreement with theory
is very good both qualitatively and quantitatively, with an evolution that is well predicted till τ ≃ 15 ; for long times
the simulation data become somewhat noisy (less statistics can be achieved due to the smaller number of particles
left in the system. The number of particles at τ = 15 is N ≃ 5000).
In Fig.s 6 and 7, results for a system with p = 0.5 are shown. Data are averaged over 64 runs. In this case,
the agreement is still good qualitatively, with similar shapes for the theoretical and numerical curves, but some
discrepancies are observed. The most significant deviation from the theoretical prediction occurs for the density,
ρkm , where the value of the minimum and also its position are not predicted accurately. However, the hydrodynamic
framework still captures correctly the trends of the complex dynamics of the perturbations. The above discrepancies
could be ascribable to the failure of the first Sonine approximation for high dissipation (i.e. “high” p) or to finite size
effects. We mention here that within the first Sonine approximation, the dimensionless coefficients µ˜ and κ˜ exhibit a
divergent behavior in the vicinity of p = 0.8 [12] which is presumably unphysical, and is an indication of the limitation
of the method. For completeness and to assess the robustness of the predictions with respect to a modification of
the numerical values of the key parameters, we have also reported in Fig.s 6 and 7 the predictions obtained when
the transport coefficients take their elastic hard disc value (e.g. µ˜(p = 0.5) ≃ 0.505 while µ˜(p = 0) = 0, as required
by Fourier’s law ). The important point is that the time evolutions are not significantly affected, and that the main
features remain the same. Additionally, as mentioned after Eq. (66) and (92), a possible source of inaccuracy lies
in the truncation of decay rates to their first order (k1) in the gradients. While the corresponding terms have been
shown to be small for inelastic hard spheres [23] (with the notable simplification there that the velocity decay rate
Page 15
15
vanishes identically due to momentum conservation), their relevance in the present case has not been assessed, apart
indirectly for the velocity decay rate δζu, by noting that second order corrections do not spoil the accuracy of the
prediction (88), see Figures 1 and 2.
VI. CONCLUSIONS
The objective here has been to explore the validity of a hydrodynamic description based on density, momentum,
and kinetic temperature fields for a gas composed of particles which annihilate with probability p or scatter elastically
otherwise, by a direct analysis of the spectrum of the linear Boltzmann equation. The motivation mainly was to study
the applicability of hydrodynamics to systems which a priori lack scale separation and in which there are no collisional
invariants. The analysis performed here has been shown to lead to results equivalent to those obtained previously
from the more formal Chapman-Enskog expansion [12], with the difference that the present approach considers linear
excitations only. However, the current spectral method is arguably more straightforward and explicitly shows that
the hydrodynamic description arises in the appropriate time scale, when the “kinetic modes” can be considered as
negligible against the hydrodynamic ones.
The eigenvalue problem of the Boltzmann operator linearized around the homogeneous decay state has been ad-
dressed and we have identified the hydrodynamic eigenfunctions. These eigenfunctions are not simply linear com-
binations of 1, v and v2 as it happens in the elastic case, but they are replaced by derivatives of the homogeneous
decay state velocity distribution function χH , that is not known analytically. As a consequence of the non-hermitian
character of the linearized Boltzmann operator, the eigenfunctions are not orthogonal. It is nevertheless possible to
construct a set of biorthonormal functions, {ξ¯β}β=1,...3, which are linear combinations of 1, v and v2, a crucial point in
order to obtain the hydrodynamic equations. The analysis is complicated by the fact that none of the {ξ¯β} functions
are left eigenfunctions of the linearized Boltzmann operator, since no quantity is conserved during binary encounters.
We have used these hydrodynamic eigenfunctions to derive to Navier-Stokes order the heat and momentum fluxes,
together with the various decay rates. To this end, we have decomposed the distribution function, δχk, into its
hydrodynamic and non-hydrodynamic parts. This decomposition enables us to close the hydrodynamic equations in
the long time limit and to order k2, and provides Green-Kubo formulas for the transport coefficients. We then arrived
at the linearized equations (around the homogeneous decay state) for the hydrodynamic fields, in the usual form of
partial differential equations with coefficients that are independent of the space variable but depend on time, since
the reference state considered is itself time dependent. If we analyze the stability of these equations, we may conclude
that small perturbations should decay in the long time limit. Nevertheless, it must be stressed that the perturbation
may increase at short times, thereby possibly leaving the linear domain where our analysis holds. The long time
dynamics in such a case remains an open question.
In section V, we have reported Molecular Dynamics simulations for the evolution of a perturbation of the transversal
and longitudinal velocity fields, which show a rich dynamics. The agreement between theory and simulations is
very good for moderate values of the annihilation probability p (say p < 0.5), which gives strong support to the
theoretical analysis developed here. The theoretical curves still agree qualitatively at larger p, with however some
quantitative discrepancies which might be a manifestation of the approximations underlying the computation of the
transport coefficients (namely κ and µ, evaluated to first order in a Sonine expansion [12]). Indeed, those coefficients
are predicted to exhibit a divergent behavior for p ∼ 0.8, but the simulations we have performed do not show a
qualitatively different behavior for these values of p. A further complication –that is also an a priori limitation for
the efficiency of a hydrodynamic approach– is that for values of p close to unity, the separation of time scales between
the kinetic and hydrodynamic modes is not clear cut (the density decay rate is on the order of the collision frequency,
comparable to the inverse typical time of the kinetic modes). To address this concern, one should study in detail the
spectrum of non hydrodynamic modes to find the slowest, and compare it to the fastest decay rate in our problem (i.e.
ζn). Such a program, left for future work, has been achieved in Appendix C within the Maxwell model framework,
with the conclusion that scale separation does not hold for p > p∗ = 1/4. While the threshold p∗ is a priori model
dependent, a similar phenomenon is to be expected in the original “hard-sphere” dynamics considered in this paper.
The coarse-grained description for p > p∗ is an open question.
Acknowledgments
We would like to thank the Agence Nationale de la Recherche for financial support. M. I. G. S. acknowledges
financial support from Becas de la Fundacio´n La Caixa y el Gobierno France´s and from the HPC-EUROPA project
(RII3-CT-2003-506079), with the support of the European Community Research Infrastructure Action.
vanishes identically due to momentum conservation), their relevance in the present case has not been assessed, apart
indirectly for the velocity decay rate δζu, by noting that second order corrections do not spoil the accuracy of the
prediction (88), see Figures 1 and 2.
VI. CONCLUSIONS
The objective here has been to explore the validity of a hydrodynamic description based on density, momentum,
and kinetic temperature fields for a gas composed of particles which annihilate with probability p or scatter elastically
otherwise, by a direct analysis of the spectrum of the linear Boltzmann equation. The motivation mainly was to study
the applicability of hydrodynamics to systems which a priori lack scale separation and in which there are no collisional
invariants. The analysis performed here has been shown to lead to results equivalent to those obtained previously
from the more formal Chapman-Enskog expansion [12], with the difference that the present approach considers linear
excitations only. However, the current spectral method is arguably more straightforward and explicitly shows that
the hydrodynamic description arises in the appropriate time scale, when the “kinetic modes” can be considered as
negligible against the hydrodynamic ones.
The eigenvalue problem of the Boltzmann operator linearized around the homogeneous decay state has been ad-
dressed and we have identified the hydrodynamic eigenfunctions. These eigenfunctions are not simply linear com-
binations of 1, v and v2 as it happens in the elastic case, but they are replaced by derivatives of the homogeneous
decay state velocity distribution function χH , that is not known analytically. As a consequence of the non-hermitian
character of the linearized Boltzmann operator, the eigenfunctions are not orthogonal. It is nevertheless possible to
construct a set of biorthonormal functions, {ξ¯β}β=1,...3, which are linear combinations of 1, v and v2, a crucial point in
order to obtain the hydrodynamic equations. The analysis is complicated by the fact that none of the {ξ¯β} functions
are left eigenfunctions of the linearized Boltzmann operator, since no quantity is conserved during binary encounters.
We have used these hydrodynamic eigenfunctions to derive to Navier-Stokes order the heat and momentum fluxes,
together with the various decay rates. To this end, we have decomposed the distribution function, δχk, into its
hydrodynamic and non-hydrodynamic parts. This decomposition enables us to close the hydrodynamic equations in
the long time limit and to order k2, and provides Green-Kubo formulas for the transport coefficients. We then arrived
at the linearized equations (around the homogeneous decay state) for the hydrodynamic fields, in the usual form of
partial differential equations with coefficients that are independent of the space variable but depend on time, since
the reference state considered is itself time dependent. If we analyze the stability of these equations, we may conclude
that small perturbations should decay in the long time limit. Nevertheless, it must be stressed that the perturbation
may increase at short times, thereby possibly leaving the linear domain where our analysis holds. The long time
dynamics in such a case remains an open question.
In section V, we have reported Molecular Dynamics simulations for the evolution of a perturbation of the transversal
and longitudinal velocity fields, which show a rich dynamics. The agreement between theory and simulations is
very good for moderate values of the annihilation probability p (say p < 0.5), which gives strong support to the
theoretical analysis developed here. The theoretical curves still agree qualitatively at larger p, with however some
quantitative discrepancies which might be a manifestation of the approximations underlying the computation of the
transport coefficients (namely κ and µ, evaluated to first order in a Sonine expansion [12]). Indeed, those coefficients
are predicted to exhibit a divergent behavior for p ∼ 0.8, but the simulations we have performed do not show a
qualitatively different behavior for these values of p. A further complication –that is also an a priori limitation for
the efficiency of a hydrodynamic approach– is that for values of p close to unity, the separation of time scales between
the kinetic and hydrodynamic modes is not clear cut (the density decay rate is on the order of the collision frequency,
comparable to the inverse typical time of the kinetic modes). To address this concern, one should study in detail the
spectrum of non hydrodynamic modes to find the slowest, and compare it to the fastest decay rate in our problem (i.e.
ζn). Such a program, left for future work, has been achieved in Appendix C within the Maxwell model framework,
with the conclusion that scale separation does not hold for p > p∗ = 1/4. While the threshold p∗ is a priori model
dependent, a similar phenomenon is to be expected in the original “hard-sphere” dynamics considered in this paper.
The coarse-grained description for p > p∗ is an open question.
Acknowledgments
We would like to thank the Agence Nationale de la Recherche for financial support. M. I. G. S. acknowledges
financial support from Becas de la Fundacio´n La Caixa y el Gobierno France´s and from the HPC-EUROPA project
(RII3-CT-2003-506079), with the support of the European Community Research Infrastructure Action.
Page 16
16
Appendix A: APPROXIMATE EXPRESSION FOR THE TRANSPORT COEFFICIENTS AND DECAY
RATES
For completeness, we recall here the explicit expressions used for the transport coefficients and decay rates, as
obtained within the first order Sonine scheme in Ref. [12]. The distribution function in this approximation reads
χH(c) =
1
pid/2
e−c
2
{
1 + a2
[
d(d + 2)
8
− d+ 2
2
c2 +
c4
2
]}
, (A1)
with the kurtosis
a2 =
8(3− 2
√
2)p
(4d + 6−
√
2)p+ 8
√
2(d− 1)(1− p)
. (A2)
This expression allows us to calculate the decay rates
ζn =
d + 2
4
(
1− a2
1
16
)
, (A3)
ζT =
d + 2
8d
(
1 + a2
8d + 11
16
)
, (A4)
and also the transport coefficients
η˜ =
1
4νη − 2pζT
, (A5)
κ˜ =
1
νκ − 2pζT
[
1
2
pζnµ˜+
d + 2
8
(2a2 + 1)
]
, (A6)
µ˜ =
1
2νκ − 3pζT − 2pζn
[
pζT κ˜+
d + 2
8
(2a2 + 1)
]
, (A7)
where the values of the coefficients νη and νκ are
νη =
p
8d
[
3 + 6d+ 2d2 − a2
278 + 375d+ 96d2 + 2d3
32(d+ 2)
]
+ (1− p)
(
1− a2
1
32
)
, (A8)
νκ =
p
32d
[
16 + 27d+ 8d2 + a2
2880 + 1544d− 2658d2 − 1539d3 − 200d4
32d(d+ 2)
]
+ (1− p)d− 1
d
(
1 + a2
1
32
)
. (A9)
Finally, the expressions for ζu,ρ and ζu,θ are
ζu,ρ =
8(d− 1)
d(d + 2)
µ˜ζu, (A10)
ζu,θ =
8(d− 1)
d(d + 2)
κ˜ζu, (A11)
with
ζu =
(d + 2)2
32(d− 1)
[
1 + a2
−86− 101d+ 32d2 + 88d3 + 28d4
32(d+ 2)
]
. (A12)
Appendix A: APPROXIMATE EXPRESSION FOR THE TRANSPORT COEFFICIENTS AND DECAY
RATES
For completeness, we recall here the explicit expressions used for the transport coefficients and decay rates, as
obtained within the first order Sonine scheme in Ref. [12]. The distribution function in this approximation reads
χH(c) =
1
pid/2
e−c
2
{
1 + a2
[
d(d + 2)
8
− d+ 2
2
c2 +
c4
2
]}
, (A1)
with the kurtosis
a2 =
8(3− 2
√
2)p
(4d + 6−
√
2)p+ 8
√
2(d− 1)(1− p)
. (A2)
This expression allows us to calculate the decay rates
ζn =
d + 2
4
(
1− a2
1
16
)
, (A3)
ζT =
d + 2
8d
(
1 + a2
8d + 11
16
)
, (A4)
and also the transport coefficients
η˜ =
1
4νη − 2pζT
, (A5)
κ˜ =
1
νκ − 2pζT
[
1
2
pζnµ˜+
d + 2
8
(2a2 + 1)
]
, (A6)
µ˜ =
1
2νκ − 3pζT − 2pζn
[
pζT κ˜+
d + 2
8
(2a2 + 1)
]
, (A7)
where the values of the coefficients νη and νκ are
νη =
p
8d
[
3 + 6d+ 2d2 − a2
278 + 375d+ 96d2 + 2d3
32(d+ 2)
]
+ (1− p)
(
1− a2
1
32
)
, (A8)
νκ =
p
32d
[
16 + 27d+ 8d2 + a2
2880 + 1544d− 2658d2 − 1539d3 − 200d4
32d(d+ 2)
]
+ (1− p)d− 1
d
(
1 + a2
1
32
)
. (A9)
Finally, the expressions for ζu,ρ and ζu,θ are
ζu,ρ =
8(d− 1)
d(d + 2)
µ˜ζu, (A10)
ζu,θ =
8(d− 1)
d(d + 2)
κ˜ζu, (A11)
with
ζu =
(d + 2)2
32(d− 1)
[
1 + a2
−86− 101d+ 32d2 + 88d3 + 28d4
32(d+ 2)
]
. (A12)
Page 17
17
Appendix B: EIGENFUNCTIONS OF Λ
In this Appendix some of the details leading to the solution of the eigenvalue problem (41) are given. Consider first
the function
Ψ1(c) = χH(c). (B1)
and let the linearized operator Λ act onto χH
Λ(c1)χH(c1) = γ
∫
dc2T (c1, c2)(1 + P12)χH(c1)χH(c2)
+ p(2ζn − dζT )χH(c1)− pζT c1 ·
∂
∂c1
χH(c1). (B2)
If we take into account the equation for χH(c1)
(dζT − 2ζn)pχH(c1) + pζT c1 ·
∂
∂c1
χH(c1) = γ
∫
dc2T (c1, c2)χH(c1)χH(c2), (B3)
we can rewrite equation (B2) as
Λ(c1)χH(c1) = (dζT − 2ζn)pχH(c1) + pζT c1 ·
∂
∂c1
χH(c1). (B4)
Consider now the function
Ψ2(c) = c ·
∂
∂c
χH(c). (B5)
In order to proceed, we perform the change of variables c1 = λc′1 in equation (B3)
(dζT − 2ζn)pχH(λc′1) + pζT c′1 ·
∂
∂c′1
χH(λc′1)
= γλd+1
∫
dc′2
∫
dσˆθ(c′12 · σˆ)c′12 · σˆ[(1− p)b−1σ − 1]χH(λc′1)χH(λc′2).
(B6)
Deriving with respect to λ we obtain
(dζT − 2ζn)p
∂
∂λ
χH(λc1) + pζT c1 ·
∂
∂c1
(
∂
∂λ
χH(λc1)
)
= (d + 1)λdγ
∫
dc2T (c1, c2)χH(λc1)χH(λc2)
+ γλd+1
∫
dc2T (c1, c2)
∂χH(λc1)
∂λ
χH(λc2)
+ γλd+1
∫
dc2T (c1, c2)χH(λc1)
∂χH(λc2)
∂λ
. (B7)
and taking λ = 1 we arrive at the equation for Ψ2(c)
(dζT − 2ζn)pΨ2(c1) + pζT c1 ·
∂
∂c1
Ψ2(c1)
= (d + 1)γ
∫
dc2T (c1, c2)χH(λc1)χH(λc2)
+ γ
∫
dc2T (c1, c2)Ψ2(c1)χH(λc2)
+ γ
∫
dc2T (c1, c2)χH(λc1)Ψ2(c2), (B8)
Appendix B: EIGENFUNCTIONS OF Λ
In this Appendix some of the details leading to the solution of the eigenvalue problem (41) are given. Consider first
the function
Ψ1(c) = χH(c). (B1)
and let the linearized operator Λ act onto χH
Λ(c1)χH(c1) = γ
∫
dc2T (c1, c2)(1 + P12)χH(c1)χH(c2)
+ p(2ζn − dζT )χH(c1)− pζT c1 ·
∂
∂c1
χH(c1). (B2)
If we take into account the equation for χH(c1)
(dζT − 2ζn)pχH(c1) + pζT c1 ·
∂
∂c1
χH(c1) = γ
∫
dc2T (c1, c2)χH(c1)χH(c2), (B3)
we can rewrite equation (B2) as
Λ(c1)χH(c1) = (dζT − 2ζn)pχH(c1) + pζT c1 ·
∂
∂c1
χH(c1). (B4)
Consider now the function
Ψ2(c) = c ·
∂
∂c
χH(c). (B5)
In order to proceed, we perform the change of variables c1 = λc′1 in equation (B3)
(dζT − 2ζn)pχH(λc′1) + pζT c′1 ·
∂
∂c′1
χH(λc′1)
= γλd+1
∫
dc′2
∫
dσˆθ(c′12 · σˆ)c′12 · σˆ[(1− p)b−1σ − 1]χH(λc′1)χH(λc′2).
(B6)
Deriving with respect to λ we obtain
(dζT − 2ζn)p
∂
∂λ
χH(λc1) + pζT c1 ·
∂
∂c1
(
∂
∂λ
χH(λc1)
)
= (d + 1)λdγ
∫
dc2T (c1, c2)χH(λc1)χH(λc2)
+ γλd+1
∫
dc2T (c1, c2)
∂χH(λc1)
∂λ
χH(λc2)
+ γλd+1
∫
dc2T (c1, c2)χH(λc1)
∂χH(λc2)
∂λ
. (B7)
and taking λ = 1 we arrive at the equation for Ψ2(c)
(dζT − 2ζn)pΨ2(c1) + pζT c1 ·
∂
∂c1
Ψ2(c1)
= (d + 1)γ
∫
dc2T (c1, c2)χH(λc1)χH(λc2)
+ γ
∫
dc2T (c1, c2)Ψ2(c1)χH(λc2)
+ γ
∫
dc2T (c1, c2)χH(λc1)Ψ2(c2), (B8)
Page 18
18
or equivalently
Λ(c1)Ψ2(c1) = −(d + 1)
[
(dζT − 2ζn)pχH(c1) + pζT c1 ·
∂
∂c1
χH(c1)
]
. (B9)
Consequently, equations (B4) and (B9) can be written as
Λ(c1)Ψ1(c1) = (dζT − 2ζn)pΨ1(c1) + pζTΨ2(c1), (B10)
Λ(c1)Ψ2(c2) = −(d+ 1)[(dζT − 2ζn)pΨ1(c1) + pζTΨ2(c1)]. (B11)
With equations (B10) and (B11) we can easily see that
Λ(c1)[(d + 1)Ψ1(c1) + Ψ2(c1)] = 0, (B12)
so that
ξ1(c1) ≡ (d + 1)Ψ1(c1) + Ψ2(c1), (B13)
is an eigenfunction of the operator Λ(c1) with eigenvalue λ1 = 0. It is also straightforward to see that
ξ2(c1) = (dζT − 2ζn)pΨ1(c1) + pζTΨ2(c1), (B14)
is an eigenfunction of Λ. We obtain that
Λ(c1)ξ2(c1) = (dζT − 2ζn)pΛ(c1)Ψ1(c1) + pζTΛ(c1)Ψ2(c1)
= −(ζT + 2ζn)pξ2(c1), (B15)
where we have used equations (B10) and (B11). Therefore, we have that ξ2(c1) is an eigenfunction of Λ(c1) with
eigenvalue λ2 = −(ζT + 2ζn)p.
Finally, let us consider the last function
Ψ3(c) = −
∂
∂c
χH(c). (B16)
Deriving the equation obeyed by χH(c−w), with respect to w and subsequently evaluating the result for w = 0, we
obtain
(dζT − 2ζn)pΨ3(c1) + pζTΨ3(c1) + pζT c1 ·
∂
∂c1
Ψ3(c1)
= γ
∫
dc2T (c1, c2)(1 + P12)Ψ3(c1)χH(c2), (B17)
or equivalently
Λ(c1)Ψ3(c1) = pζTΨ3(c1). (B18)
In other words, ξ3(c1) ≡ Ψ3(c1) is an eigenfunction of Λ(c1) with eigenvalue λ3 = pζT .
Appendix C: LINEARIZED BOLTZMANN OPERATOR FOR MAXWELL MOLECULES
The objective in this Appendix is to study the spectrum of the linearized Boltzmann operator for Maxwell molecules
with annihilation. It will be shown that for 0 ≤ p ≤ p∗, where p∗ depends on the specific Maxwell model under
consideration, the norm of the hydrodynamic eigenvalues are smaller than the rest of the spectrum.
The main characteristic of Maxwell models is that the differential cross section multiplied by the relative velocity is
independent of the relative velocity. We are going to assume that it is also independent of the angle between the two
colliding particles. Then, the Boltzmann equation for a system of Maxwell molecules which annihilate in a collision
with probability p and collide elastically otherwise (with probability 1− p) reads
(
∂
∂t
+ v1 · ∇
)
f(r,v1, t) = −pβΩ
∫
dv2f(r,v1, t)f(r,v2, t)
+ (1− p)β
∫
dv2
∫
dσˆ[b−1σ − 1]f(r,v1, t)f(r,v2, t),
(C1)
or equivalently
Λ(c1)Ψ2(c1) = −(d + 1)
[
(dζT − 2ζn)pχH(c1) + pζT c1 ·
∂
∂c1
χH(c1)
]
. (B9)
Consequently, equations (B4) and (B9) can be written as
Λ(c1)Ψ1(c1) = (dζT − 2ζn)pΨ1(c1) + pζTΨ2(c1), (B10)
Λ(c1)Ψ2(c2) = −(d+ 1)[(dζT − 2ζn)pΨ1(c1) + pζTΨ2(c1)]. (B11)
With equations (B10) and (B11) we can easily see that
Λ(c1)[(d + 1)Ψ1(c1) + Ψ2(c1)] = 0, (B12)
so that
ξ1(c1) ≡ (d + 1)Ψ1(c1) + Ψ2(c1), (B13)
is an eigenfunction of the operator Λ(c1) with eigenvalue λ1 = 0. It is also straightforward to see that
ξ2(c1) = (dζT − 2ζn)pΨ1(c1) + pζTΨ2(c1), (B14)
is an eigenfunction of Λ. We obtain that
Λ(c1)ξ2(c1) = (dζT − 2ζn)pΛ(c1)Ψ1(c1) + pζTΛ(c1)Ψ2(c1)
= −(ζT + 2ζn)pξ2(c1), (B15)
where we have used equations (B10) and (B11). Therefore, we have that ξ2(c1) is an eigenfunction of Λ(c1) with
eigenvalue λ2 = −(ζT + 2ζn)p.
Finally, let us consider the last function
Ψ3(c) = −
∂
∂c
χH(c). (B16)
Deriving the equation obeyed by χH(c−w), with respect to w and subsequently evaluating the result for w = 0, we
obtain
(dζT − 2ζn)pΨ3(c1) + pζTΨ3(c1) + pζT c1 ·
∂
∂c1
Ψ3(c1)
= γ
∫
dc2T (c1, c2)(1 + P12)Ψ3(c1)χH(c2), (B17)
or equivalently
Λ(c1)Ψ3(c1) = pζTΨ3(c1). (B18)
In other words, ξ3(c1) ≡ Ψ3(c1) is an eigenfunction of Λ(c1) with eigenvalue λ3 = pζT .
Appendix C: LINEARIZED BOLTZMANN OPERATOR FOR MAXWELL MOLECULES
The objective in this Appendix is to study the spectrum of the linearized Boltzmann operator for Maxwell molecules
with annihilation. It will be shown that for 0 ≤ p ≤ p∗, where p∗ depends on the specific Maxwell model under
consideration, the norm of the hydrodynamic eigenvalues are smaller than the rest of the spectrum.
The main characteristic of Maxwell models is that the differential cross section multiplied by the relative velocity is
independent of the relative velocity. We are going to assume that it is also independent of the angle between the two
colliding particles. Then, the Boltzmann equation for a system of Maxwell molecules which annihilate in a collision
with probability p and collide elastically otherwise (with probability 1− p) reads
(
∂
∂t
+ v1 · ∇
)
f(r,v1, t) = −pβΩ
∫
dv2f(r,v1, t)f(r,v2, t)
+ (1− p)β
∫
dv2
∫
dσˆ[b−1σ − 1]f(r,v1, t)f(r,v2, t),
(C1)
Page 19
19
where β is a constant representing the microscopic scattering collision frequency, Ω = 2pid/2/Γ(d/2) is the d-
dimensional solid angle and the operator b−1σ is defined in the main text, Eq. (4).
If we consider the homogeneous case, it is straightforward to see that the temperature is a constant in time,
TH(t) = TH(0), and that the density evolves as
nH(t) =
nH(0)
1 + pζnt
, (C2)
with ζn = βΩnH(0). Moreover, it was shown in [24] that there is an exact mapping between the homogeneous equation
for Maxwell molecules with annihilation (arbitrary p) and the usual elastic Maxwell molecules (p = 0). If we introduce
the time scale
s(t) =
∫ t
0
dt′
nH(t′)
nH(0)
, (C3)
it can be seen that the distribution function for arbitrary p is
f(v, t) =
nH(t)
nH(0)
fE [v, (1 − p)s(t)] , (C4)
where nH(t) is given by formula (C2), s(t) by (C3), and the function fE is the distribution function for elastic Maxwell
molecules. Note that this relation is also valid for p = 1 where fE is frozen in the initial condition.
In the elastic case, every homogeneous distribution tends to relax to a Maxwellian after a transient time. Then,
the same is going to happen for arbitrary p < 1 because of the mapping. In this sense, we can consider that the state
analogous to the homogeneous decay state introduced in the main text for hard particles and that will constitute the
appropriate reference will be characterized by the distribution function
fH(v, t) =
nH(t)
vdH
χM (v/vH), (C5)
where vH =
( 2TH
m
)1/2
and χM (v) = 1/pid/2e−v
2
is the Maxwellian distribution. Note that in the homogeneous decay
state, both the density and temperature decay, whether in the present case only the density decays. Then, as vH
plays no role, we will consider units with vH = 1 for simplicity.
Let us study now the linear response to an inhomogeneous small perturbation around the reference state as it was
done in section II. If we introduce the scaled distribution function
δχ(r,v, τ) =
nH(0)
nH(t)
[f(r,v, t) − fH(v, t)], (C6)
the equation for δχ in the s scale defined in (C3) is
[
∂
∂s
+ hH(s)v1 · ∇
]
δχ(r,v, s) = Λ(v1)δχ(r,v, s), (C7)
where we have introduced the function hH(s) = nH(0)/nH(s), and the linearized Boltzmann operator
Λ(v1)g(v1) = (1− p)
ζn
Ω
∫
dv2
∫
dσˆ[b−1σ − 1](1 + P12)χM (v1)g(v2)
− pζnχM (v1)
∫
dv2g(v2). (C8)
Let us stress that, although there is an exact mapping for the full non-linear homogeneous equation between p = 0
and arbitrary p, no such mapping exists for the linear inhomogeneous Boltzmann equation, Eq. (C7). Then, as in
the main text, the possibility of an hydrodynamic description depends on the properties of the linearized Boltzmann
operator. Here we will see that it is possible to calculate all the eigenfunctions and eigenvalues of this operator and
that there is a region of the parameter p in which we have an appropriate scale separation.
Let us write the linearized Boltzmann operator as
Λ(v1)g(v1) = (1− p)ΛE(v1)g(v1)− pζnχM (v1)
∫
dv2g(v2), (C9)
where β is a constant representing the microscopic scattering collision frequency, Ω = 2pid/2/Γ(d/2) is the d-
dimensional solid angle and the operator b−1σ is defined in the main text, Eq. (4).
If we consider the homogeneous case, it is straightforward to see that the temperature is a constant in time,
TH(t) = TH(0), and that the density evolves as
nH(t) =
nH(0)
1 + pζnt
, (C2)
with ζn = βΩnH(0). Moreover, it was shown in [24] that there is an exact mapping between the homogeneous equation
for Maxwell molecules with annihilation (arbitrary p) and the usual elastic Maxwell molecules (p = 0). If we introduce
the time scale
s(t) =
∫ t
0
dt′
nH(t′)
nH(0)
, (C3)
it can be seen that the distribution function for arbitrary p is
f(v, t) =
nH(t)
nH(0)
fE [v, (1 − p)s(t)] , (C4)
where nH(t) is given by formula (C2), s(t) by (C3), and the function fE is the distribution function for elastic Maxwell
molecules. Note that this relation is also valid for p = 1 where fE is frozen in the initial condition.
In the elastic case, every homogeneous distribution tends to relax to a Maxwellian after a transient time. Then,
the same is going to happen for arbitrary p < 1 because of the mapping. In this sense, we can consider that the state
analogous to the homogeneous decay state introduced in the main text for hard particles and that will constitute the
appropriate reference will be characterized by the distribution function
fH(v, t) =
nH(t)
vdH
χM (v/vH), (C5)
where vH =
( 2TH
m
)1/2
and χM (v) = 1/pid/2e−v
2
is the Maxwellian distribution. Note that in the homogeneous decay
state, both the density and temperature decay, whether in the present case only the density decays. Then, as vH
plays no role, we will consider units with vH = 1 for simplicity.
Let us study now the linear response to an inhomogeneous small perturbation around the reference state as it was
done in section II. If we introduce the scaled distribution function
δχ(r,v, τ) =
nH(0)
nH(t)
[f(r,v, t) − fH(v, t)], (C6)
the equation for δχ in the s scale defined in (C3) is
[
∂
∂s
+ hH(s)v1 · ∇
]
δχ(r,v, s) = Λ(v1)δχ(r,v, s), (C7)
where we have introduced the function hH(s) = nH(0)/nH(s), and the linearized Boltzmann operator
Λ(v1)g(v1) = (1− p)
ζn
Ω
∫
dv2
∫
dσˆ[b−1σ − 1](1 + P12)χM (v1)g(v2)
− pζnχM (v1)
∫
dv2g(v2). (C8)
Let us stress that, although there is an exact mapping for the full non-linear homogeneous equation between p = 0
and arbitrary p, no such mapping exists for the linear inhomogeneous Boltzmann equation, Eq. (C7). Then, as in
the main text, the possibility of an hydrodynamic description depends on the properties of the linearized Boltzmann
operator. Here we will see that it is possible to calculate all the eigenfunctions and eigenvalues of this operator and
that there is a region of the parameter p in which we have an appropriate scale separation.
Let us write the linearized Boltzmann operator as
Λ(v1)g(v1) = (1− p)ΛE(v1)g(v1)− pζnχM (v1)
∫
dv2g(v2), (C9)
Page 20
20
-1
-0,8
-0,6
-0,4
-0,2
0
0 0,2 0,4 0,6 0,8 1
p
λ10/ζn λ01/ζn
λ00/ζn
λ
κ
/ζ
n
Figure 8: Hydrodynamic eigenvalues, λ00, λ10, and λ01, and the slowest kinetic eigenvalue, λk, as a function of the dissipation
parameter p. The eigenvalues are normalized by ζn = βΩnH(0).
where we have introduced the linearized Boltzmann operator for elastic Maxwell particles
ΛE(v1) =
ζn
Ω
∫
dv2
∫
dσˆ[b−1σ − 1](1 + P12)χM (v1)g(v2), (C10)
whose spectral properties are well known [25, 26]. In particular, for d = 3 its eigenfunctions are
φrlm(v) = ArlχM (v)Srl+1/2(v
2)vlYlm(θ, ϕ), r = 0, 1, . . . (C11)
where Ylm(θ, ϕ) are the spherical harmonics, functions of the polar angles (θ, ϕ) of v with respect to an arbitrary
direction, Srl+1/2(v
2) are the Sonine polynomials which satisfy
∫ ∞
0
dxe−xSnq (x)S
n′
q (x) =
Γ(n+ q + 1)
n!
δnn′ , (C12)
and Arl are some constants that are introduced in order to normalize the eigenfunctions and that play no role in
the following analysis. The eigenvalues of ΛE , λErl, are also known. It can be seen that λ
E
00, λ
E
10, and λ
E
01 (which is
3 times degenerate) vanish, corresponding to the five hydrodynamic eigenvalues and that the slowest kinetic mode
corresponds to an eigenvalue λEk = −ζn/3.
The important point here is that the functions φrlm are also eigenfunctions of the second operator in (C9). Taking
into account the orthogonality properties of the spherical harmonics and of the Sonine polynomials, Eq. (C12), one
has
∫
dvφrlm(v) = δr0δl0. (C13)
Then, as φ000(v) = χM (v), the functions φrlm(v) are in fact the eigenfunctions of the total linearized Boltzmann
operator for arbitrary p. With the aid of (C13) is straightforward to see that the eigenvalues are
λrl = (1− p)λErl − pζnδr0δl0. (C14)
In Fig. 8 we plot the hydrodynamic eigenvalues as well as the slowest kinetic eigenvalue, λk, as a function of the
dissipation parameter p. It can be seen that for 0 ≤ p < 1/4, there is scale separation in the sense that the three
modes (density, linear momentum and kinetic energy) retained in the coarse-grained description decay slower that any
of the other “kinetic” modes. On the other hand, for 1/4 < p ≤ 1, the largest kinetic eigenvalue is slower than λ00.
We therefore conclude here that a conservative requirement for the validity of our approach in the case of Maxwell
molecules would be p < 1/4.
-1
-0,8
-0,6
-0,4
-0,2
0
0 0,2 0,4 0,6 0,8 1
p
λ10/ζn λ01/ζn
λ00/ζn
λ
κ
/ζ
n
Figure 8: Hydrodynamic eigenvalues, λ00, λ10, and λ01, and the slowest kinetic eigenvalue, λk, as a function of the dissipation
parameter p. The eigenvalues are normalized by ζn = βΩnH(0).
where we have introduced the linearized Boltzmann operator for elastic Maxwell particles
ΛE(v1) =
ζn
Ω
∫
dv2
∫
dσˆ[b−1σ − 1](1 + P12)χM (v1)g(v2), (C10)
whose spectral properties are well known [25, 26]. In particular, for d = 3 its eigenfunctions are
φrlm(v) = ArlχM (v)Srl+1/2(v
2)vlYlm(θ, ϕ), r = 0, 1, . . . (C11)
where Ylm(θ, ϕ) are the spherical harmonics, functions of the polar angles (θ, ϕ) of v with respect to an arbitrary
direction, Srl+1/2(v
2) are the Sonine polynomials which satisfy
∫ ∞
0
dxe−xSnq (x)S
n′
q (x) =
Γ(n+ q + 1)
n!
δnn′ , (C12)
and Arl are some constants that are introduced in order to normalize the eigenfunctions and that play no role in
the following analysis. The eigenvalues of ΛE , λErl, are also known. It can be seen that λ
E
00, λ
E
10, and λ
E
01 (which is
3 times degenerate) vanish, corresponding to the five hydrodynamic eigenvalues and that the slowest kinetic mode
corresponds to an eigenvalue λEk = −ζn/3.
The important point here is that the functions φrlm are also eigenfunctions of the second operator in (C9). Taking
into account the orthogonality properties of the spherical harmonics and of the Sonine polynomials, Eq. (C12), one
has
∫
dvφrlm(v) = δr0δl0. (C13)
Then, as φ000(v) = χM (v), the functions φrlm(v) are in fact the eigenfunctions of the total linearized Boltzmann
operator for arbitrary p. With the aid of (C13) is straightforward to see that the eigenvalues are
λrl = (1− p)λErl − pζnδr0δl0. (C14)
In Fig. 8 we plot the hydrodynamic eigenvalues as well as the slowest kinetic eigenvalue, λk, as a function of the
dissipation parameter p. It can be seen that for 0 ≤ p < 1/4, there is scale separation in the sense that the three
modes (density, linear momentum and kinetic energy) retained in the coarse-grained description decay slower that any
of the other “kinetic” modes. On the other hand, for 1/4 < p ≤ 1, the largest kinetic eigenvalue is slower than λ00.
We therefore conclude here that a conservative requirement for the validity of our approach in the case of Maxwell
molecules would be p < 1/4.
Page 21
21
Appendix D: EVALUATION OF δζn, δζu AND δζT
In this Appendix we calculate the contribution of the hydrodynamic part of δχk to the functionals δζn, δζu and
δζT defined in (38)-(40). To this end, we write explicitly Pδχk as
Pδχk(c1) = ρkξ1(c1) +
(
1
2
θk + ρk
)
ξ2(c1) +wk · ξ3(c1)
= ρkχH(c1)−
1
2
θk
∂
∂c1
· [c1χH(c1)]−wk ·
∂
∂c1
χH(c1). (D1)
Let us first evaluate δζn[Pδχk]. After some algebra it can be seen that
δζn[χH(c1)] = −4ζn, (D2)
δζn
[
∂
∂c1
· [c1χH(c1)]
]
= 2ζn, (D3)
δζn
[
∂
∂c1
χH(c1)
]
= 0, (D4)
where we have used equations (B5), (B8), the definition of the density decay rate, equation (14), and symmetry
considerations. Then, if we consider equations (D2)-(D4) we finally obtain
δζn[Pδχk] = −4ζnρk − ζnθk. (D5)
Now let us calculate δζu[Pδχk]. Using the definition of the density decay rate, equation (14), and symmetry
considerations, it appears that
δζu[χH(c1)] = 0, (D6)
δζu
[
∂
∂c1
· (c1χH(c1))
]
= 0, (D7)
δζui
[
∂
∂c1j
χH(c1)
]
= δijδζui
[
∂
∂c1i
χH(c1)
]
= 2ζn. (D8)
Then, we have
δζu[Pδχk] = −2ζnwk. (D9)
Finally, we turn to δζT [Pδχk]. Using the definitions of the decay rates, equations (14) and (15), we obtain
δζT [χH(c1)] = −4ζT , (D10)
δζT
[
∂
∂c1
· (c1χH(c1))
]
= 6ζT + 4ζn, (D11)
δζT
[
∂
∂c1
χH(c1)
]
= 0 (D12)
from which it follows that
δζT [Pδχk] = −4ζTρk − (3ζT + 2ζn)θk. (D13)
[1] I. Goldhirsch, Ann. Rev. Fluid Mech. 35, 267 (2005).
[2] J. Dufty, arXiv:0709.0479 (2007).
[3] A. Barrat, E. Trizac, and M. H. Ernst, J. Phys.: Condens. Matter 17, S2429 (2005).
[4] J. Brey and D. Cubero, in Granular gases, edited by S. Luding and T. Po¨schel (Springer, Berlin, 2001).
[5] E. Ben-Naim, P. Krapivsky, F. Leyvraz, and S. Redner, J. Chem. Phys. 98, 7284 (1994).
[6] R. Blythe, M. R. Evans, and Y. Kafri, Phys. Rev. Lett. 85, 3759 (2000).
[7] P. Krapivsky and C. Sire, Phys. Rev. Lett. 86, 2494 (2001).
Appendix D: EVALUATION OF δζn, δζu AND δζT
In this Appendix we calculate the contribution of the hydrodynamic part of δχk to the functionals δζn, δζu and
δζT defined in (38)-(40). To this end, we write explicitly Pδχk as
Pδχk(c1) = ρkξ1(c1) +
(
1
2
θk + ρk
)
ξ2(c1) +wk · ξ3(c1)
= ρkχH(c1)−
1
2
θk
∂
∂c1
· [c1χH(c1)]−wk ·
∂
∂c1
χH(c1). (D1)
Let us first evaluate δζn[Pδχk]. After some algebra it can be seen that
δζn[χH(c1)] = −4ζn, (D2)
δζn
[
∂
∂c1
· [c1χH(c1)]
]
= 2ζn, (D3)
δζn
[
∂
∂c1
χH(c1)
]
= 0, (D4)
where we have used equations (B5), (B8), the definition of the density decay rate, equation (14), and symmetry
considerations. Then, if we consider equations (D2)-(D4) we finally obtain
δζn[Pδχk] = −4ζnρk − ζnθk. (D5)
Now let us calculate δζu[Pδχk]. Using the definition of the density decay rate, equation (14), and symmetry
considerations, it appears that
δζu[χH(c1)] = 0, (D6)
δζu
[
∂
∂c1
· (c1χH(c1))
]
= 0, (D7)
δζui
[
∂
∂c1j
χH(c1)
]
= δijδζui
[
∂
∂c1i
χH(c1)
]
= 2ζn. (D8)
Then, we have
δζu[Pδχk] = −2ζnwk. (D9)
Finally, we turn to δζT [Pδχk]. Using the definitions of the decay rates, equations (14) and (15), we obtain
δζT [χH(c1)] = −4ζT , (D10)
δζT
[
∂
∂c1
· (c1χH(c1))
]
= 6ζT + 4ζn, (D11)
δζT
[
∂
∂c1
χH(c1)
]
= 0 (D12)
from which it follows that
δζT [Pδχk] = −4ζTρk − (3ζT + 2ζn)θk. (D13)
[1] I. Goldhirsch, Ann. Rev. Fluid Mech. 35, 267 (2005).
[2] J. Dufty, arXiv:0709.0479 (2007).
[3] A. Barrat, E. Trizac, and M. H. Ernst, J. Phys.: Condens. Matter 17, S2429 (2005).
[4] J. Brey and D. Cubero, in Granular gases, edited by S. Luding and T. Po¨schel (Springer, Berlin, 2001).
[5] E. Ben-Naim, P. Krapivsky, F. Leyvraz, and S. Redner, J. Chem. Phys. 98, 7284 (1994).
[6] R. Blythe, M. R. Evans, and Y. Kafri, Phys. Rev. Lett. 85, 3759 (2000).
[7] P. Krapivsky and C. Sire, Phys. Rev. Lett. 86, 2494 (2001).
Page 22
22
[8] E. Trizac, Phys. Rev. Lett. 88, 160601 (2002).
[9] J. Piasecki, E. Trizac, and M. Droz, Phys. Rev. E 66, 066111 (2002).
[10] A. Lipowski, D. Lipowska, and A. Feirrera, Phys. Rev. E 73, 032102 (2006).
[11] S. Chapman and T. G. Cowling, The mathematical theory of nonuniform gases (Cambridge University Press, London,
1960).
[12] F. Coppex, M. Droz, and E. Trizac, Phys. Rev. E 70, 061102 (2004).
[13] J. J. Brey and M. J. Ruiz-Montero, Phys. Rev. E 69, 011305 (2004).
[14] O. E. Lanford, Physica A 106, 70 (1981).
[15] F. Coppex, M. Droz, and E. Trizac, Phys. Rev. E 69, 011303 (2004).
[16] J. J. Brey, J. W. Dufty, and M. J. Ruiz-Montero, in Granular Gas Dynamics, edited by T. Po¨schel and N. Brilliantov
(Springer, Berlin, 2003).
[17] J. J. Brey and J. W. Dufty, Phys. Rev. E 72, 011303 (2005).
[18] J. Dufty, arXiv:0707.111, submitted to Journal of Physical Chemistry B (2007).
[19] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Oxford Science Publications, Bristol, 1987).
[20] J. M. Montanero, A. Santos, and V. Garzo´, Physica A 376, 75 (2007).
[21] V. Garzo´, A. Santos, and J. M. Montanero, Physica A 376, 94 (2007).
[22] J. J. Brey, M. I. Garc´ıa de Soria, and P. Maynar, (to be published).
[23] J. J. Brey, J. W. Dufty, C. S. Kim, and A. Santos, Phys. Rev. E 58, 4638 (1998).
[24] A. Santos and J. J. Brey, Phys. Fluids 29, 1750 (1986).
[25] P. Re´sibois and M. de Leener, Classical Kinetic Theory of Fluids (John Wiley, New York, 1977).
[26] J. A. McLennan, Introduction to Nonequilibrium Statistical Mechanics (Prentice-Hall, Englewood Cliffs, NJ, 1989).
[8] E. Trizac, Phys. Rev. Lett. 88, 160601 (2002).
[9] J. Piasecki, E. Trizac, and M. Droz, Phys. Rev. E 66, 066111 (2002).
[10] A. Lipowski, D. Lipowska, and A. Feirrera, Phys. Rev. E 73, 032102 (2006).
[11] S. Chapman and T. G. Cowling, The mathematical theory of nonuniform gases (Cambridge University Press, London,
1960).
[12] F. Coppex, M. Droz, and E. Trizac, Phys. Rev. E 70, 061102 (2004).
[13] J. J. Brey and M. J. Ruiz-Montero, Phys. Rev. E 69, 011305 (2004).
[14] O. E. Lanford, Physica A 106, 70 (1981).
[15] F. Coppex, M. Droz, and E. Trizac, Phys. Rev. E 69, 011303 (2004).
[16] J. J. Brey, J. W. Dufty, and M. J. Ruiz-Montero, in Granular Gas Dynamics, edited by T. Po¨schel and N. Brilliantov
(Springer, Berlin, 2003).
[17] J. J. Brey and J. W. Dufty, Phys. Rev. E 72, 011303 (2005).
[18] J. Dufty, arXiv:0707.111, submitted to Journal of Physical Chemistry B (2007).
[19] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Oxford Science Publications, Bristol, 1987).
[20] J. M. Montanero, A. Santos, and V. Garzo´, Physica A 376, 75 (2007).
[21] V. Garzo´, A. Santos, and J. M. Montanero, Physica A 376, 94 (2007).
[22] J. J. Brey, M. I. Garc´ıa de Soria, and P. Maynar, (to be published).
[23] J. J. Brey, J. W. Dufty, C. S. Kim, and A. Santos, Phys. Rev. E 58, 4638 (1998).
[24] A. Santos and J. J. Brey, Phys. Fluids 29, 1750 (1986).
[25] P. Re´sibois and M. de Leener, Classical Kinetic Theory of Fluids (John Wiley, New York, 1977).
[26] J. A. McLennan, Introduction to Nonequilibrium Statistical Mechanics (Prentice-Hall, Englewood Cliffs, NJ, 1989).
Sign up today - FREE
Mendeley saves you time finding and organizing research. Learn more
- All your research in one place
- Add and import papers easily
- Access it anywhere, anytime
Start using Mendeley in seconds!
Readership Statistics
2 Readers on Mendeley
by Discipline
100% Physics
by Academic Status
50% Ph.D. Student
50% Associate Professor
by Country
50% Brazil
50% France