format: text

text:
# A model for the atomic-scale structure of a dense, nonequilibrium fluid: the homogeneous cooling state of granular fluids

## I Introduction

The atomic scale structure, e.g. the pair distribution function or, equivalently, the density-density equal time correlation function, of both equilibrium and nonequilibrium fluids is directly accessible to experiment by means of light scattering and has been used to study the behavior of complex systems such as sheared colloidal suspensions. However, the subject of nonequilibrium fluid structure remains obscure, particularly when compared to the structure of equilibrium fluids which is one of the most advanced areas of equilibrium statistical mechanics. While kinetic theory can provide some asymptotic results, see e.g. ,,, models of the small-scale structure remain phenomenological. This is not surprising: in equilibrium, one starts with the exact N-body distribution and the problem is to integrate out uninteresting degrees of freedom so as to get the pair distribution function (pdf). Away from equilibrium, even in steady states, the N-body distribution function is not known and the only starting point is either a complex dynamical equation, the Liouville equation, or equivalently a coupled set of equations, the BBGKY hierarchy relating the n-body to the (n+1)-body distribution functions. Another method used to derive equilibrium integral equations is based on the fact that from the BBGKY hierarchy it is easy, essentially by integrating over the velocities, to derive the YBG hierarchy which relates the pair distribution function to the triplet distribution function, the triplet to the quartic distribution function, etc. Some sort of closure hypothesis, like the neglect of three body correlations, then yields a closed equation for the pdf, see e.g.,. However, away from equilibrium, the distribution of velocities is not known and this procedure cannot be followed in detail. One approach to nonequilibrium fluid structure is that of Hess and coworkers,. The object there was to describe sheared colloidal suspensions as well as computer simulations of sheared simple fluids. The theory, based on a kind of relaxation approximation, is phenomenological and will be considered in more detail below. Fluctuating hydrodynamics offers another possibility,, but is restricted to length-scales much larger than atomic length scales. Another approach, developed about a decade later by Eu et al,, was to try to generalize the equilibrium integral equations, like the Kirkwood and the Percus-Yevick equations , to nonequilibrium fluids. For the reasons just mentioned, the only way to carry out this program was to start with an ansatz for the n-body distribution which allows the development to parallel that of the equilibrium equations. Whether or not this is a good approach is difficult to assess since the approximations made are uncontrolled, but it has been previously shown, and will be argued in more detail below, that velocity correlations are the crucial determinant of atomic-scale deviations from equilibrium structure and, it would follow, that approaches which do consider them may be missing an essential contribution to the structure. Indeed, the role played by velocity correlations in determining one aspect of the structure of sheared simple fluids, namely the value and angular dependence of the pair distribution function at contact has already been demonstrated and the purpose of this paper is to show that the knowledge of such velocity correlations can be used to extend equilibrium structural models into the nonequilibrium domain.

In order to avoid a number of complexities associated with the spatial anisotropy of sheared fluids, the system used here to illustrate this model is the dissipative hard sphere model used in the study of granular fluids. This fundamentally nonequilibrium system presents a unique opportunity for the application and development of nonequilibrium statistical mechanics in a system which is both of practical relevance and yet sufficiently simple to be amenable to theoretical analysis. The basis for such analysis has recently been formulated by Brey, Dufty and Santos who have constructed a Liouville description of the dynamics of the dissipative hard-sphere model. From this, Boltzmann- and Enskog-level kinetic equations for the one-body distribution function (the dense-fluid generalization of the Boltzmann equation) follow and their properties have been studied theoretically, in particular, the Chapman-Enskog solution to it has been developed, giving explicit expressions for the pressure and transport properties of the system, and it is possible to solve it numerically. Despite the simplicity of the dynamics, simulations have shown that this system exhibits a rich phenomenology which is only partially understood at present (see, e.g., ). The starting point for describing the phenomenology is the so-called homogeneous cooling state (HCS): since the collisions dissipate energy, the analog of the equilibrium state is one which is spatially homogeneous but in which the temperature decreases algebraically with time. However, this state is relatively unstable as it is subject to a number of forms of spontaneous symmetry breaking. The simplest example of this is the so-called clustering instability in which sufficiently large systems, subject to no external forces, exhibit a spatial clustering. This has been shown to be due to a small wave-vector hydrodynamic instability. Small systems for which the maximum wavevector (as determined by the size of the simulation cell) is too large, do not exhibit the transition. Other phenomena include the ”shearing state”, in which the system remains structurally homogeneous but long-ranged momentum correlations develop, and ”inelastic collapse” in which a few particles collide with one-another many times until all relative momentum is dissipated and the particles come to rest in contact with one another. One motivation for trying to understand the structure of the HCS is that it is needed as input to various kinetic theory calculations which could be useful in understanding some of these features.

The remainder of this paper is organized as follows. In Section II, the relevant elements of kinetic theory are reviewed and it is shown that the pdf at contact can in general be explicitly calculated with the same level of approximation as is used to get the Enskog equation. It is noted that both the pdf at contact, as well as the pressure, can be calculated in the HCS with no further approximations even though the exact solution for the one-body distribution is not known. In Section III, the dynamics of the dissipative hard sphere model is reformulated, through a change of variables, so that the HCS is mapped onto a steady state. Although mathematically, equivalent, the steady state formulation is convenient since it avoids numerical problems associated with the rapid decay of kinetic energy in the HCS and also allows statistical properties, like the pressure, to be computed by replacing ensemble averages by time averages (the assumption of ergodicity). Section IV consists of a presentation of a model for the structure of the HCS and Section V presents a comparison of this model to the results of simulations. The paper concludes with a discussion of the results and their bearing on the questions raised above.

## II Theory

The model granular fluid consists of a collection of $`N`$ classical hard-spheres of diameter $`\sigma `$ having positions $`𝐪`$ and momenta $`𝐩`$ which experience inelastic binary collisions. Between collisions, the atoms freely stream so the state evolves according to:

$`m_i{\displaystyle \frac{d}{dt}}𝐪_i`$ $`=`$ $`𝐩_i`$ (1)
$`{\displaystyle \frac{d}{dt}}𝐩_i`$ $`=`$ $`𝐅_i`$ (2)
where we allow for the possibility of an external force acting on the particles. Henceforth, we will consider only identical atoms and will use units in which the mass is equal to one. When two particles, say $`i`$ and $`j`$, collide, their total momentum is unaffected but their relative momentum is altered according to

$$𝐩_{ij}\widehat{b}_{ij}𝐩_{ij}𝐩_{ij}\left(1+\alpha \right)𝐩_{ij}\widehat{𝐪}_{ij}$$
(3)
where, in general, relative quantities are denoted as $`𝐩_{ij}=𝐩_i𝐩_j`$ , unit vectors as $`\widehat{𝐪}𝐪/\left|𝐪\right|`$ and where this equation serves to define the momentum transfer operator $`\widehat{b}_{ij}`$ . The parameter $`\alpha `$ is called the coefficient of restitution and takes on values between one (elastic hard spheres) and zero (completely dissipative inelastic collisions). Under this dynamics, the one-body distribution function, $`f_1(x,t)f_1(𝐪,𝐩,t)`$ and the two-body distribution function $`f_2(x_1,x_2,t)`$ are related by the first equation of the BBGKY hierarchy

$`{\displaystyle \frac{}{t}}f_1(x_1,t)+𝐩_1{\displaystyle \frac{}{𝐪_1}}f_1(x_1,t)`$ (4)
$`=`$ $`\sigma ^2{\displaystyle 𝑑x_2_\mathrm{\Omega }𝑑\widehat{\sigma }\delta \left(𝐪_{12}\sigma \widehat{\sigma }\right)\left(\widehat{\sigma }𝐩_{12}\right)\left[\alpha ^2\widehat{b}_{ij}^1+1\right]\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}\right)f_2(x_1,x_2,t)}`$ (5)
where the notation indicates that integration of $`\widehat{\sigma }`$ is over the unit sphere, $`\mathrm{\Theta }(x)=1`$ for $`x>0`$ and is otherwise zero and where $`\widehat{b}_{ij}^1`$ is the inverse of $`\widehat{b}_{ij}`$ and may be easily seen to be

$$\widehat{b}_{ij}^1𝐩_{ij}𝐩_{ij}\left(1+\alpha ^1\right)𝐩_{ij}\widehat{𝐪}_{ij}\widehat{𝐪}_{ij}$$
(6)
Finally, the two-body distribution must satisfy the identity, a kind of boundary condition,

$`\delta \left(𝐪_{12}\sigma \widehat{\sigma }\right)\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}\right)f_2(x_1,x_2,t)`$ $`=`$ $`\delta \left(𝐪_{12}\sigma \widehat{\sigma }\right)\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}\right){\displaystyle \frac{1}{\alpha ^2}}\widehat{𝐛}_{12}^1f_2(x_1,x_2,t)`$ (7)
$`=`$ $`\delta \left(𝐪_{12}\sigma \widehat{\sigma }\right){\displaystyle \frac{1}{\alpha ^2}}\widehat{b}_{12}^1\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}\right)f_2(x_1,x_2,t)`$ (8)
the origin of which - basically, the conservation of probability during a collision - is discussed in the appendix. Here, I only note that this identity is independent of, and provides information additional to, the first BBGKY hierarchy equation given above. In fact, as discussed in the appendix, it can actually be derived from the second equation of the BBGKY hierarchy. Using an abbreviated notation, $`W_{12}`$ $`\delta \left(𝐪_{12}\sigma \widehat{\sigma }\right)`$, it is easily seen that the full distribution at contact can be written as

$`W_{12}f_2(x_1,x_2,t)`$ $`=`$ $`W_{12}\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}\right)f_2(x_1,x_2,t)+W_{12}\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}\right)f_2(x_1,x_2,t)`$ (9)
$`=`$ $`W_{12}\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}\right)f_2(x_1,x_2,t)+W_{12}{\displaystyle \frac{1}{\alpha ^2}}\widehat{b}_{12}^1\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}\right)f_2(x_1,x_2,t)`$ (10)
The combination $`\delta \left(𝐪_{12}\sigma \widehat{\sigma }\right)\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}\right)f_2(𝐱_1,𝐱_2,t)`$ appearing here, as well as in the collisional term in eq.(5) refers to the probability for two atoms just prior to a collision and we refer to it as the pre-collisional part of the distribution. Equation (11) thus expresses the two-body distribution at contact solely in terms of the pre-collisional distribution. As discussed in ref. , the assumption of ”molecular chaos” used to obtain the Boltzmann equation and its generalization to dense hard-sphere fluids, the Enskog equation, is that this pre-collisional distribution can be approximated by neglecting velocity correlations between the particles so that one writes

$$W_{12}\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}\right)f_2(x_1,x_2,t)W_{12}\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}\right)g_0(𝐪_1,𝐪_2)f_1(x_1,t)f_1(x_2,t)$$
(11)
where $`g_0(𝐪_1,𝐪_2;t)`$ is the pre-collisional pair distribution function (pdf) at time $`t`$ which is normally taken to be the local-equilibrium form,. Using this in eq. (5), the generalized Enskog equation immediately results

$`{\displaystyle \frac{}{t}}f_1(x_1,t)+𝐩_1{\displaystyle \frac{}{𝐪_1}}f_1(x_1,t)`$ (12)
$`=`$ $`\sigma ^2{\displaystyle 𝑑x_2𝑑\widehat{\sigma }\delta \left(𝐪_{12}\sigma \widehat{\sigma }\right)\left(\widehat{\sigma }𝐩_{12}\right)\left[\alpha ^2\widehat{b}_{ij}^1+1\right]\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}\right)g_0(𝐪_1,𝐪_2)f_1(x_1,t)f_1(x_2,t)}`$ (13)
while from the boundary condition, we get

$`W_{12}f_2(x_1,x_2,t)`$ $`=`$ $`W_{12}\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}\right)g_0(𝐪_1,𝐪_2)f_1(x_1,t)f_1(x_2,t)`$ (15)
$`+W_{12}{\displaystyle \frac{1}{\alpha ^2}}\widehat{b}_{12}^1\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}\right)g_0(𝐪_1,𝐪_2)f_1(x_1,t)f_1(x_2,t)`$
$`=`$ $`W_{12}g_0(𝐪_1,𝐪_2)f_1(x_1,t)f_1(x_2,t)`$ (17)
$`+W_{12}\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}\right)\left({\displaystyle \frac{1}{\alpha ^2}}\widehat{b}_{12}^11\right)g_0(𝐪_1,𝐪_2)f_1(x_1,t)f_1(x_2,t)`$
thus expressing the distribution at contact as the sum of an uncorrelated term, the first on the right, and a term expressing the corrections due to momentum correlations. In equilibrium, the second term vanishes. Naturally, the momenta of the atoms are correlated after a collision and, in fact, (17) allows us to determine these correlations in a manner consistent with the degree of approximation of the generalized Enskog equation. Although such an evaluation might be used in a number of different ways, I will here focus one particular application of it to the problem of understanding the structure of the non-equilibrium state. Specifically, integrating over momenta gives

$`W_{12}n(𝐪_1;t)n(𝐪_2;t)g(𝐪_1,𝐪_2;t)`$ (18)
$`=`$ $`W_{12}n(𝐪_1)n(𝐪_2)g_0(𝐪_1,𝐪_2;t)`$ (20)
$`+W_{12}g_0(𝐪_1,𝐪_2){\displaystyle 𝑑𝐩_1𝑑𝐩_2\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}\right)\left(\frac{1}{\alpha ^2}\widehat{b}_{12}^11\right)f_1(x_1,t)f_1(x_2,t)}`$
where the nonequilibrium density is $`n(𝐪_1;t)=𝑑𝐩_1f_1(x_1,t)`$. Eq. (20) thus gives us an approximate evaluation of the contribution of momentum correlations to the structure of the fluid as expressed through the pdf. This relation was first used in ref.  to characterize velocity correlations in a sheared fluid and has recently been used in a study of velocity correlations in granular fluids near equilibrium.

The macroscopic balance equations for the density, momentum field $`m𝐔`$, and energy density $`\frac{3}{2}k_BT`$ follow from (5) by multiplying by $`1`$, $`𝐩_1`$, and $`\frac{1}{2m}p_1^2`$ respectively and integrating over $`𝐩_1`$ to get

$`{\displaystyle \frac{}{t}}n+n𝐔`$ $`=`$ $`0`$ (21)
$`{\displaystyle \frac{}{t}}U_i+𝐔U_i+(mn)^1_jP_{ij}`$ $`=`$ $`0`$ (22)
$`{\displaystyle \frac{}{t}}T+𝐔T+{\displaystyle \frac{2}{3nk_B}}(P_{ij}_jU_i+𝐪)`$ $`=`$ $`T^{3/2}\zeta `$ (23)
and explicit forms of the pressure tensor $`P_{ij}`$ and heat flux vector $`𝐪`$ are given in the literature. The source term in the temperature equation is due to the cooling caused by the inelastic collisions. It is easy to see that a spatially homogeneous solution to these equations is possible with

$`n(𝐫,t)`$ $`=`$ $`n_0`$ (24)
$`𝐔`$ $`=`$ $`\mathrm{𝟎}`$ (25)
$`{\displaystyle \frac{}{t}}T`$ $`=`$ $`T^{3/2}\zeta `$ (26)
where $`n_0`$ is a constant. Because there is no potential energy in hard-sphere systems, the time and energy scales are set by the temperature and so, on dimensional grounds, it is clear that $`\zeta `$ is independent of temperature: the temperature is therefore given by $`T(t)=T_0\left[1+\frac{1}{2}\zeta T_0^{1/2}t\right]^2`$ thus explicitly demonstrating the cooling. A linear stability analysis of the equations (23) expanded about the homogeneous cooling state shows that the state is unstable against small wavevector fluctuations.

The pdf at contact for HCS can also be evaluated without explicit knowledge of the one body distribution (see the Appendix) and is

$`\chi \left(\widehat{\sigma }\right)`$ $``$ $`{\displaystyle \frac{1}{V}}{\displaystyle 𝑑𝐪_1𝑑𝐪_2\delta \left(𝐪_{12}\sigma \widehat{\sigma }\right)g(𝐪_1,𝐪_2;t)}`$ (27)
$`=`$ $`{\displaystyle \frac{1+\alpha }{2\alpha }}\chi _0`$ (28)
where $`\chi _0`$ is defined by an analogous expression to this with $`g(𝐪_1,𝐪_2;t)`$ replaced by $`g_0(𝐪_1,𝐪_2;t)`$. Finally. since it will be of use below, we quote the expression for the collisional contribution to the pressure which is

$`{\displaystyle \frac{1}{nk_BT}}p^c`$ $`=`$ $`{\displaystyle \frac{1}{3}}{\displaystyle \underset{i=1}{\overset{3}{}}}{\displaystyle \frac{1}{nk_BT}}P_{ij}^c`$ (29)
$`=`$ $`{\displaystyle \frac{1+\alpha }{4}}\sigma ^3{\displaystyle _\mathrm{\Omega }}𝑑\widehat{\sigma }{\displaystyle 𝑑x_1𝑑x_2\delta \left(q_{12}\sigma \right)\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}\right)\left(\widehat{𝐪}_{12}𝐩_{12}\right)^2f_1(x_1,t)f_1(x_2,t)g_0(𝐪_1,𝐪_2)}`$ (30)
$`=`$ $`{\displaystyle \frac{1+\alpha }{2}}\chi _0`$ (31)
where the last line follows from the spatial isotropy and time invariance of the HCS and also does not require explicit knowledge of the one-body distribution. This is useful since an exact solution to the Enskog equation for the one-body distribution for HCS is not known.

## III Mapping HCS to a steady state

Since the time scale is determined by the temperature, we can simplify the description of this state by making a change of variables. Specifically, define a scaled time coordinate via $`ws=\mathrm{ln}(t/t_0)`$ where $`w`$ and $`t_0`$ are arbitrary constants, and the corresponding equations of motion are

$`{\displaystyle \frac{d}{ds}}𝐪`$ $`=`$ $`𝐜_i`$ (32)
$`{\displaystyle \frac{d}{ds}}𝐜_i`$ $`=`$ $`w𝐜_i`$ (33)
where $`\stackrel{}{c}_i=wt_0e^{ws}𝐩_i`$ is the momentum (velocity) in the new coordinates. To carry through the statistical description, we must define a new set of probability densities. If we denote the set of positions and velocities of atoms $`x_1\mathrm{}x_m`$ in the original system by $`\mathrm{\Gamma }_m`$ and of those in the scaled system by $`\mathrm{\Gamma }_m^{}`$ then the m-body distribution will transform according to

$$f_m(\mathrm{\Gamma }_m,t)=J\stackrel{~}{f}_m(\mathrm{\Gamma }_m^{},t)$$
(34)
where $`J=\left|\frac{\mathrm{\Gamma }_m^{}}{\mathrm{\Gamma }_m^{}}\right|=\left(wt_0e^{ws}\right)^{mD}`$ is the Jacobian of the transformation for a D-dimensional system. We then find, e.g., that eq.(5) becomes

$`{\displaystyle \frac{}{s}}\stackrel{~}{f}_1(x_1^{},s)+𝐜_1{\displaystyle \frac{}{𝐪_1}}\stackrel{~}{f}_1(x_1^{},s)+{\displaystyle \frac{}{𝐜_1}}w𝐜_1\stackrel{~}{f}_1(x_1^{},s)`$ (35)
$`=`$ $`\sigma ^2{\displaystyle 𝑑x_1^{}𝑑x_2^{}𝑑\widehat{\sigma }\delta \left(𝐪_{12}\sigma \widehat{\sigma }\right)\left(\widehat{\sigma }𝐜_{12}\right)\left[\alpha ^2\widehat{b}_{ij}^1+1\right]\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐜_{12}\right)\stackrel{~}{f}_2(x_1^{}x_2^{},t)}`$ (36)
so that the only change is that a new term appears in the streaming operator. This term has the same form as the artificial thermostats used in the study of shear flow: however, here it arises solely from a change of variables. The balance equations become

$`{\displaystyle \frac{}{t}}\stackrel{~}{n}+\stackrel{~}{n}𝐕`$ $`=`$ $`0`$ (37)
$`{\displaystyle \frac{}{t}}V_i+𝐕V_i+(m\stackrel{~}{n})^1_j\stackrel{~}{P}_{ij}wV_i`$ $`=`$ $`0`$ (38)
$`{\displaystyle \frac{}{t}}\stackrel{~}{T}+𝐕\stackrel{~}{T}+{\displaystyle \frac{2}{3nk_B}}(\stackrel{~}{P}_{ij}_jV_i+\stackrel{~}{𝐪})2w\stackrel{~}{T}`$ $`=`$ $`\stackrel{~}{T}^{3/2}\zeta `$ (39)
where $`𝐕(𝐪,s)=𝑑𝐜𝐜\stackrel{~}{f}_1(x_1^{},s)`$, $`\frac{3}{2}k_B\stackrel{~}{T}(𝐪,s)=𝑑𝐜\frac{1}{2}mc^2\stackrel{~}{f}_1(x_1^{},s)`$, etc. Now, the temperature of the homogeneous state is given by

$$\stackrel{~}{T}(t)=\left(\frac{2w}{\zeta }\right)^2\left(1+\left(\frac{2w}{\zeta \sqrt{\stackrel{~}{T}(0)}}1\right)e^{wt}\right)^2$$
(40)
so that any initial temperature will equilibrate to a final stable value of $`\left(\frac{2w}{\zeta }\right)^2`$. Thus, the original homogeneous cooling state is mapped by this change of variables onto a superficially steady state. It is worth noting that if one were to return to the original form of the dynamics and were to periodically rescale the momenta of the atoms so as to restore the temperature to its original value, then it is easy to see that as the time interval between rescalings goes to zero, the resulting equations of motion can be written in the form 33 with the constant replaced by a complicated function of the momenta. Periodic rescaling is commonly used in simulations of sheared fluids and has recently been employed in the simulation of HCS,.

This mapping also makes the hydrodynamic instability apparent. If we expand (39) about the steady state and retain terms to second order in the gradients (i.e., the Navier-Stokes equations) and to first order in the densities and transform to Fourier space (with wavevector $`𝐤`$), it is easy to see that the vorticity, $`\omega =\widehat{𝐤}\times 𝐕`$, satisfies

$$_t\omega +\left(\nu _0k^2w\right)\omega =0$$
(41)
where $`\nu _0`$ is the shear viscosity evaluated at density $`n_0`$ and temperature $`\stackrel{~}{T}_0=\left(\frac{2w}{\zeta }\right)^2`$. It is obviously unstable for sufficiently small wavevectors and for sufficiently large systems, we therefore expect a spontaneous shear to develop (note that $`\nu _0\sqrt{\stackrel{~}{T}_0}`$ so that the arbitrary constant, $`w`$, plays no role in the stability criterion). It should also be noted that there is a closely related instability in the total ($`𝐤=\mathrm{𝟎}`$) velocity in the mapped system which, in the linear stability analysis, obeys

$$\frac{}{t}V_i(\mathrm{𝟎})wV_i(\mathrm{𝟎})=0$$
(42)
and so is clearly unstable for all system sizes. Since we are at liberty to choose the initial conditions to be those for which $`V_i(\mathrm{𝟎})=0`$, we will find that this is represents a minor problem in the simulations and is of no theoretical significance.

## IV Modelling the structure of the HCS

In this section, we will assume spatial homogeneity so that the pdf depends only on the scalar separation between atoms. The simplest realistic model for the structure of the equilibrium hard sphere fluid is the Percus-Yevick approximation(see e.g., refs.,). This consists of the Ornstein-Zernike equation

$$h(𝐫)=c(𝐫)+\rho 𝑑𝐫^{}c(|𝐫𝐫^{}|)h(r^{})$$
(43)
where $`h(r)=g(r)1`$ and $`c(r)`$ is the direct correlation function, together with the boundary conditions

$`\mathrm{\Theta }(\sigma r)h(r)`$ $`=`$ $`1`$ (44)
$`\mathrm{\Theta }(r\sigma )c(r)`$ $`=`$ $`0`$ (45)
where the first condition is exact while the second defines the approximation. Comparison with computer simulation shows that the PY approximation for the pair distribution function is quite accurate for separations greater that about two hard-sphere diameters but less accurate near contact. The description of the small-separation structure can be significantly improved by considering the Yukawa closure for the Orstein-Zernicke equation which replaces the boundary condition on the direct correlation function by

$$\mathrm{\Theta }(r1)c(r)=\underset{i=1}{\overset{m}{}}K_i\frac{e^{v_ir}}{r}$$
(46)
and choosing the constants $`K_i`$ and $`v_i`$ to reproduce known properties: for example, taking $`m=1`$ and fitting the Carnahan-Starling equation of state as calculated by both the pressure equation and the compressibility equation. The original Mean Spherical Approximation, for arbitrary pair potentials $`\mathrm{\Phi }(r)`$, consists of requiring that $`\mathrm{\Theta }(r\sigma )c(r)=\mathrm{\Phi }(r)`$ where the effective hard sphere diameter is fit according to some criterion. The PY approximation is then seen to be the MSA for hard spheres. Equation (46) may therefore be viewed as either the MSA for a potential which is the sum of Yukawas or as a general expansion with coefficients to be fitted in which case it is termed the Generalized MSA or GMSA and can be viewed as being systematic since any function could be fitted as a sum of Yukawas.

Note that eq.(43) defines the direct correlation function and that the first of the boundary conditions in eq.(44) is an exact requirement. The only way in which this model uses the assumption that the system being modeled is in equilibrium is through the arguments that lead to the conclusion that the direct correlation function is short-ranged and hence the justification for the boundary conditions in eq.(44) and eq.(46). This connection to the equilibrium state is made even weaker in a reformulation of the model due to Yuste and Santos, ,. They begin by noting that the Laplace transform of the quantity $`rg(r)`$, in the PY approximation, is naturally written as

$`G(t)`$ $``$ $`{\displaystyle _0^{\mathrm{}}}𝑑re^{sr}rg(r)`$ (47)
$`=`$ $`{\displaystyle \frac{tF(t)e^t}{112\eta F(t)e^t}}`$ (48)
with

$$F(t)=\frac{1+A_1t}{S_0+S_1t+S_2t^2+S_3t^3}$$
(49)
They go on to point out that given the second equality of eq.(48) and making a Padé’ approximation for the function $`F(t)`$ one can deduce the correct order of the numerator and denominator of $`F(t)`$ as well as the PY expressions for the coefficients based solely on the asymptotic properties of the pdf. Specifically, they note that a)$`g(r)`$ at contact is given by $`g(\sigma )=`$ $`lim_t\mathrm{}t^2F(t)`$ thus fixing the relative number of terms in the numerator and denominator; b) $`lim_r\mathrm{}g(r)=1`$ implies that $`G(t)\underset{t0}{}t^2`$ and c) the fact that the static structure factor, given by $`S(q)=lim_{tiq}Re\left(tG(t)\right)`$, is finite at $`q=0`$ implies, given the previous condition, that for small $`t`$, $`G(t)=t^2+o(1)`$. (The last condition is equivalent to assuming that long-range correlations do not exist.) The minimal approximant satisfying these conditions is that given in eq.(49) with the PY value for the coefficients. They also note that the extension of the Padé’ approximant to include one more term in both the numerator and denominator is exactly equivalent to the one-Yukawa closure while the further extension of the approximation corresponds to a closure consisting of a sum of Yukawas. This method is also shown to give the PY solution for sticky hard spheres as well as the exact structure for both ordinary and sticky hard spheres in one-dimension. A straight forward extension of these ideas has also been used to model the square-well fluid. Thus, in this formulation, the PY form of $`G(t)`$ is taken as an ansatz characteristic of hard-core systems and the function $`F(t)`$ modeled as a Padé’ approximant subject to whatever knowledge exists about the structure.

With this justification in mind, we consider the application of this approach to the nonequilibrium HCS. In equilibrium, the next inclusion of an additional term in the numerator and denominator of $`F(z)`$ introduces two new parameters which are used to fit a known equation of state (normally the Carnahan-Starling equation of state) through both the pressure equation and the compressibility equation. The pressure equation for hard spheres in three dimensions reads

$$\frac{p}{nk_BT}=1+4\eta \chi _{eq}$$
(50)
and so allows to calculate the pdf at contact, recall $`\chi _{eq}=g(\sigma )`$, from the equation of state. In a nonequilibrium state, the collisional boundary condition can be used, together with assumption of molecular chaos, to give the same information. The equilibrium compressibility equation is

$`\left({\displaystyle \frac{}{\rho }}\beta P\right)^1`$ $`=`$ $`1+\rho {\displaystyle 𝑑\stackrel{}{r}(g(r)1)}`$ (51)
$`=`$ $`1+\underset{t0}{lim}\left(G(t){\displaystyle \frac{1}{t^2}}\right)`$ (52)
and for this, there is no obvious substitute for the nonequilibrium state. With nothing to use in its place, I will continue to apply this even in the nonequilibrium state, calculating the pressure from eq.(31), which might be viewed as a local-equilibrium approximation. In fact, the resulting model is relatively insensitive to the value used for the pressure since this only fixes the area of the structure function whereas the results are quite sensitive to the value of the pdf at contact. Using $`g_0(𝐪_1,𝐪_2;t)g_{eq}(q_{12})`$ as is normally done in Enskog theory, the model is given by

$$F(t)=\frac{1}{12\eta }\frac{1+A_1t+A_2t^2}{S_0+S_1t+S_2t^2+S_3t^3+S_4t^4}$$
(53)
with

$`S_0`$ $`=`$ $`1`$ (54)
$`S_1`$ $`=`$ $`A_11`$ (55)
$`S_2`$ $`=`$ $`A_2A_1+{\displaystyle \frac{1}{2}}`$ (56)
$`S_3`$ $`=`$ $`A_2+{\displaystyle \frac{1}{2}}A_1{\displaystyle \frac{1+2\eta }{12\eta }}`$ (57)
$`S_4`$ $`=`$ $`{\displaystyle \frac{1}{2}}A_2\left({\displaystyle \frac{1+2\eta }{12\eta }}\right)A_1+{\displaystyle \frac{2+\eta }{24\eta }}`$ (58)
$`A_1`$ $`=`$ $`{\displaystyle \frac{1}{2}}+\sqrt{{\displaystyle \frac{1}{12\eta }}{\displaystyle \frac{\left(\eta 1\right)^2Z\left(6\eta g(1;\alpha )+1\right)}{\left(2+\eta \right)2g(1;\alpha )\left(\eta 1\right)^2}}}`$ (59)
$`A_2`$ $`=`$ $`g(1;\alpha ){\displaystyle \frac{\left(1+2\eta \right)A_1\frac{1}{2}\left(2+\eta \right)}{1+6\eta g(1;\alpha )}}`$ (60)
and

$$Z\left(\frac{}{\rho }\beta P\right)^1=\left(1+\frac{1+\alpha }{2}\left(Z_{eq}1\right)\right)^1$$
(61)
with the Carnahan-Starling expression

$$\chi _{eq}=\frac{1\frac{1}{2}\eta }{\left(1\eta \right)^3}$$
(62)
together with eqs.(31) and (28) completes the model. This then reduces to the GMSA in equilibrium and can be seen as its natural generalization to a spatially isotropic nonequilibrium state.

## V Molecular Dynamics Simulations

To determine the structure of the model granular system, I have performed molecular dynamics simulations of a small systems of 108 and 500 particles in 3 dimensions governed by the steady-state dynamics described by eq.(33) and subject to periodic boundary conditions. The density was taken to be $`n^{}=0.5`$: high enough that finite density effects are important but low enough that the Enskog approximation is expected to be valid. The choice to simulate the steady-state dynamics, rather than to simulate the ”real” dynamics of the cooling system, was made on the basis that the systems cool very rapidly so that the time scales involved in the simulations become very large, the velocities and energies very small and numerical inaccuracies due to round-off error are a significant problem. This could be dealt with by periodically rescaling the velocities (i.e., redefining the time unit as in ,) but it is far more elegant and efficient to directly simulate the steady-state dynamics. Furthermore, the changes needed to implement this starting with a code for simulating equilibrium hard spheres are minimal. One point that does require attention is the instability with respect to the total momentum. Even if initial conditions are chosen so that the total momentum is zero at the start of the simulation, round-off errors lead to the spontaneous appearance of a non-zero total momentum which then quickly goes due to the instability. This effect is, however, benign and is easily suppressed by calculating the total momentum during each propagation step and subtracting $`(1/N)`$ of its value from the momentum of each particle.

The starting point for the simulations was an equilibrium configuration of velocities and positions. The value of the thermostat constant, $`w`$, was set arbitrarily. For each value of $`\alpha `$, the simulations were ”equilibrated” for a total of $`310^6`$ collisions and then statistics gathered over a second series of $`310^6`$ collisions . To obtain steady-state averages of one-body properties, quantities were time-averaged over periods of $`10^4`$ collisions throughout the simulations; these samples were then treated as statistically independent estimates and their average and standard error computed. The errors in all quantities reported below are found to be small, less than 1%. The determination of collisional effects, such as in the pressure and the pdf at contact, is somewhat different. For any collisional quantity of the form

$$A\underset{i<j}{}A(x_i,x_j)\delta (q_{ij}1)$$
(63)
the ergodic assumption gives

$`A`$ $`=`$ $`{\displaystyle \frac{1}{T}}{\displaystyle _0^T}𝑑tA(t)`$ (64)
$`=`$ $`{\displaystyle \underset{i<j}{}}{\displaystyle \frac{1}{T}}{\displaystyle _0^T}𝑑tA(x_i,x_j)\delta (q_{ij}1)`$ (65)
$`=`$ $`{\displaystyle \underset{i<j}{}}{\displaystyle \frac{1}{2T}}{\displaystyle _0^T}𝑑t\delta (\tau _{ij}1)\left({\displaystyle \frac{A(ij)}{\left|𝐪_{ij}𝐩_{ij}\right|}}+{\displaystyle \frac{A(x_i^{},x_j^{})}{\left|𝐪_{ij}𝐩_{ij}^{}\right|}}\right)`$ (66)
where the first two lines integrate the total time dependence of the function $`A`$. The third line follows from a change of variable in the delta-function $`\stackrel{}{p}_{ij}`$($`\stackrel{}{p}_{ij}^{}`$) is the relative momentum of the colliding pair immediately before(after) the collision and $`\tau _{ij}`$ is the time at which the pair (i,j) collide (which could be imaginary or outside the range of integration indicating in either case that they do not collide). The last expression obviously reduces to a sum over collisions:

$$A=\frac{1}{2T}\underset{collisions,\gamma }{}\left(\frac{A(\gamma )}{\left|𝐪_\gamma 𝐩_\gamma \right|}+\frac{A(\gamma ^{})}{\left|𝐪_\gamma 𝐩_\gamma ^{}\right|}\right)_{t_\gamma }$$
(67)
where $`\gamma `$ represents the colliding pair and which is the form used to evaluate the collisional part of the pdf at contact, $`\frac{N(N1)}{2V}4\pi g(1)=_{i<j}\delta (q_{ij}1)`$. Using $`\left|𝐪_\gamma 𝐩_\gamma ^{}\right|=\alpha \left|𝐪_\gamma 𝐩_\gamma \right|`$, this becomes

$$g(1)=\frac{1}{2\pi n}\left(\frac{\alpha +1}{2\alpha }\right)\frac{1}{T}\underset{collisions,\gamma }{}\frac{1}{\left|𝐪_\gamma \stackrel{}{p}_\gamma \right|}\mathrm{\Theta }\left(𝐪_\gamma 𝐩_\gamma \right)$$
(68)
while for the pressure one finds

$$p^c=\frac{\left(1+\alpha \right)}{2T}\underset{collisions,\gamma }{}\left|𝐪_\gamma 𝐩_\gamma \right|\mathrm{\Theta }\left(𝐪_\gamma 𝐩_\gamma \right)$$
(69)
and in both cases, the step function indicates that the expression is evaluated with the pre-collisional momenta. Finally, we present below determinations of the pdf for finite separations. These are determined in the obvious way by looping over all pairs of atoms and creating a histogram of the separations. The bin size used was 0.025 (hard sphere diameters) and these were compiled every $`10^4`$ collisions and the results averaged to obtain the final histogram.

Figure 1 shows the pdf at contact as determined from the simulations and from the collisional boundary condition. For $`\alpha >0.6`$, the agreement is seen to be good but it becomes worse values below this: furthermore, there appears to be a strong number dependence to the results with the larger system diverging more rapidly from the prediction. In both cases, simulations are only possible for $`\alpha `$ above some threshold: below this, the simulation code fails due to the time between collisions becoming smaller and smaller until the machine precision is reached. A detailed analysis of the sequence of collisions shows that this is due to a small number of particles with virtually no momentum relative to one another colliding over and over again - in other words, this is the phenomena of elastic collapse described by McNamara and Young. The threshold for the collapse is in the range $`0.3<\alpha <0.4`$ for the 500 particle system and $`0.2<\alpha <0.3`$ for the 108 particle system. In both cases, prior to the collapse, the value of the pdf at contact is several times the highest values shown in figure 1.

Similar behavior is seen in the pair distribution function. For reference, Fig. 2 shows the pdf at equilibrium as determined from the 500 atom simulations and from the model: the agreement is seen to be excellent as is also the case for the data coming from the 108 atom system. Figures 3 and 4 show the nonequilibrium part of the pdf (i.e., $`g(r)g_{eq}(r)`$) for the 108 atom system as determined by simulation and by from the nonequilibrium GMSA for $`\alpha =0.7`$ and $`0.5`$ respectively and in both cases, the extended MSA is seen to give a good quantitative description of all features of the nonequilibrium structure. (In fact, since the MD results are, by their nature, binned, the theoretical curve is obtained by integrating the model pdf over bins of the same size and position as used in the simulations.) Figures 5 and 6 show that, not surprisingly, the agreement is less good for the 500 atom system, particularly at the smaller value of $`\alpha `$.

The disparity between the results for the two systems is much greater than one finds in equilibrium and suggests that a qualitative difference between them. One obvious possible source of such a difference is the hydrodynamic instability discussed above. Using the values for the transport coefficients given by , one finds the critical size curve shown in Fig. 7 which indicates that the 108 atom system is always stable but that the 500 atom system becomes unstable around $`\alpha <0.6`$ however, knowledge of the transport coefficients at small $`\alpha `$ is only approximate so these numbers may only be indicative of the position of the instability. Nevertheless, the importance of the instability in the larger system is easily confirmed and Fig. 8, showing the velocities in one direction versus the positions along another as taken from a snapshot of the 500 atom system with $`\alpha =0.5`$, shows a spontaneously formed shearing profile. In a larger system, this would manifest itself in the form of vortices. Further evidence of the instability can be found in the kinetic contributions to the pressure tensor where, beginning at $`\alpha =0.7`$ in the 500 atom system, an oscillation develops whereby a large fraction, on the order of $`2/3`$, of the kinetic energy is concentrated in first one component of the pressure tensor and then another indicating that macroscopic flows are forming. The obvious interpretation is that around $`\alpha =0.7`$ the shear mode is soft or unstable. There is no evidence of such an unstable mode for any value of at $`\alpha `$ in the 108 atom system. This picture is thus in qualitative agreement with the predictions based on the calculations described above and has recently been observed in other studies.

We can suppress the instability in a crude way by periodically adjusting the velocities of the atoms. In these ”constrained” simulations we interpose correction whereby after every 100 collisions, we calculate the amplitude of the longest wavelength Fourier modes of the systems (i.e., $`𝐀_l=\frac{2}{N}_{i=1}^N𝐜_i\mathrm{cos}𝐤_l𝐪_i`$ for $`𝐤_1=\left(2\pi /L\right)\widehat{𝐱}`$, etc.) and we then subtract the mode from each atom’s velocity ( $`𝐜_i𝐜_i_{l=1}^3𝐀_l\mathrm{cos}𝐤_l𝐪_i`$). This is a crude procedure in that the amplitudes of the modes are only approximately set to zero and it also has the effect of removing kinetic energy from the system (which is, however, masked by the input of kinetic energy coming from the equations of motion). A more elegant procedure could be devised based on standard nonequilibrium molecular dynamics techniques such as Gauss’ principle of least constraint, but as the present purpose is only to control the unstable mode, the crude method was deemed sufficient. As shown in Fig. 9, the result is to give better agreement in the measured value of the pdf at contact between the two systems while having relatively little effect in the smaller system except at the highest values of alpha. Figure 10 shows that the pdf as determined from the constrained simulation of the larger system is in considerably better agreement with the model.

## VI Discussion

The main purpose of this paper has been to show that the GMSA can be extended to nonequilibrium systems by replacing the equilibrium input required by the GMSA with accessible nonequilibrium information coming from the collisional boundary condition and, incidently, to elucidate the atomic-scale structure of the HCS of granular fluids. In order to compare the predicted values of the pdf at contact and the pressure with the results of simulations, the equations of motion of the dissipative hard sphere system were mapped onto those that describe a steady state thus allowing us to use standard methods of steady state simulation such as the replacement of ensemble averages by time averages. The comparisons with molecular dynamics simulations also show that the pair distribution function at contact can be used as a signal of the onset of elastic collapse - its value steadily diverges from the predicted value as the elastic collapse threshold is approached and its value in the simulations that feature the collapse is very large.

The significant differences observed between the 108 and 500 atom systems were seen to be largely due to the soft hydrodynamic modes present in the larger system. Nevertheless, even when these are accounted for, there remains a significant deviation of the simulation results from the various predictions of the Enskog theory. It is tempting to conclude that this is due to a poor estimation of the quantity $`\chi _0`$, for which we use the equilibrium value, but examination of the pressure shows that the answer cannot be this simple: the pdf at contact would require a larger value for $`\chi _0`$ that increases with $`\alpha `$ in order to be in agreement with the simulations whereas Fig. 11 shows that pressure would require a smaller value that decreased with $`\alpha `$ (in both the unconstrained and constrained simulations). We thus conclude that the deviations are due to the Enskog approximation itself and could probably be described via mode coupling. Nevertheless, that $`\chi _0`$ should depend on $`\alpha `$ is intuitively clear: atoms moving slowly away from a collision will be more likely to be knocked, by a third atom, into a second collision with one another leading to such a dependence.

As mentioned in the Introduction, one phenomenological approach to the description of nonequilibrium structure is that of Hess and coworkers,. In its simplest form, this reduces to a relaxation model for the nonequilibrium contribution to the pdf:

$$\frac{}{t}g(𝐫,t)+𝐯(𝐫,t)g(𝐫,t)=\tau ^1\left(g(𝐫,t)g_0(𝐫)\right)$$
(70)
where $`\tau `$ is a relaxation time and $`g_0(𝐫)`$ is taken to be the equilibrium pdf. It is clear that for a homogeneous steady state with no flow this gives the trivial result that $`g(𝐫,t)=g_0(𝐫)`$. However, this model is intended only as a simplification of a more complex model given by

$$\frac{}{t}g(𝐫,t)+𝐯(𝐫,t)g(𝐫,t)+D\left(g_0(𝐫,t)\left(g(𝐫,t)/g_0(𝐫)\right)\right)=0$$
(71)
where $`D`$ is a kind of diffusion constant. Again, for a homogeneous steady state with no flow, only the last term survives so that this reduces to

$$g_0(r)\frac{1}{r^2}\frac{d}{dr}r^2\frac{d}{dr}\frac{g(r)}{g_0(r)}+\frac{dg_0(r)}{dr}\frac{d}{dr}\frac{g(r)}{g_0(r)}=0$$
(72)
which has the solution

$`g(r)`$ $`=`$ $`g_0(r)+\left({\displaystyle \frac{g(1)}{g_0(1)}}1\right){\displaystyle \frac{f(r)}{f(1)}}g_0(r)`$ (73)
$`f(r)`$ $`=`$ $`{\displaystyle _r^{\mathrm{}}}𝑑r{\displaystyle \frac{r^2}{g_0(r)}}`$ (74)
where we have used the boundary condition $`lim_r\mathrm{}g(r)=1`$. Although this has an interesting structure, the function $`f(r)`$ is positive-definite so that for HCS the difference $`g(r)g_0(r)`$ is also which precludes a description of the oscillatory nature of the differences found in the simulation.

Another approach is the theory of Eu et al,, also mentioned in the introduction. This theory is based on an ansatz for the N-body nonequilibrium distribution and takes the form of an integral equation:

$`\mathrm{ln}y(𝐪_1,𝐪_2;t)`$ $`=`$ $`n{\displaystyle 𝑑𝐪_3f_{NE}(𝐪_1,𝐪_3;t)y(𝐪_1,𝐪_3;t)}`$ (76)
$`\times \left[y(𝐪_2,𝐪_3;t)\left(1+f_{NE}(𝐪_2,𝐪_3;t)\right)1\right]`$
where

$`g(𝐪_1,𝐪_2;t)`$ $`=`$ $`\mathrm{exp}\left[V_{NE}(𝐪_1,𝐪_2;t)\right]y(𝐪_1,𝐪_2;t)`$ (77)
$`f_{NE}(𝐪_1,𝐪_3;t)`$ $`=`$ $`\mathrm{exp}\left[V_{NE}(𝐪_1,𝐪_2;t)\right]1`$ (78)
and where the nonequilibrium potential $`V_{NE}(𝐪_1,𝐪_2;t)`$ is a sum of the equilibrium potential and terms related to the moments of the velocity. This formulation has the desirable property that, if linearized in the density by replacing $`\mathrm{ln}y(𝐪_1,𝐪_2;t)y(𝐪_1,𝐪_2;t)1`$, it reduces to the PY approximation in equilibrium. It is difficult, however, to see how to incorporate the exact requirement of the information coming from the collisional boundary condition and it is therefore an open question whether this can give a model comparable to the nonequilibrium GMSA described above. Indeed, in the derivation of this model, the authors explicitly neglect velocity correlations of the kind used here to control the model of the structure. Nevertheless, the solution of this model to allow for such a comparison would be of some interest.

It is also appropriate to comment on the applicability of the Enskog-level description of the system. The density chosen for the simulations is one at which the Enskog description of elastic hard spheres is very good with most quantities being accurately predicted to within a few percent. Furthermore, neither the calculation of the pressure nor the pdf at contact requires explicit knowledge of the one-body distribution at contact which is good, because no exact solution to the Enskog equation for HCS exists. None the less, the results described above show that there are systematic deviations from the Enskog predictions at all values of the coefficient of restitution. Evidence has also been given that these are at least partly due to the presence of soft modes within the system. The conclusion is therefore mixed: while the Enskog description seems to be qualitatively accurate, all but the very smallest systems will contain soft or unstable modes that result in significant deviations from it. This, more than the separation of time scales discussed in ref. would seem to be the greatest obstacle to using a Boltzmann/Enskog description of the one body function or of using a hydrodynamic description of the macroscopic state. Indeed, the Enskog and hydrodynamic descriptions are successful in predicting fairly well the location of the hydrodynamic instability in the vorticity. A better test of these issues would be to study a related system, such as a sheared granular fluid, which may be more stable.

The nonequilibrium GMSA described here works surprisingly well. Encouraging results have also been found when this model was applied to simple sheared fluids and a systematic study of this system is in progress. It is also of interest to try to improve on the use of the compressibility equation in the nonequilibrium model: one substitute would be information coming from kinetic theory such as the asymptotic behavior of the pdf for which various approaches are possible.

## A Origin of the collisional boundary condition

### 1 From BBGKY hierarchy

It is instructive to derive the collisional boundary condition from the BBGKY hierarchy, the first equation of which appears above as eq.(5). The second BBGKY equation is

$`\left[{\displaystyle \frac{}{t}}+𝐩_1{\displaystyle \frac{}{𝐪_1}}+𝐩_2{\displaystyle \frac{}{𝐪_2}}\overline{T}_{}(x_1,x_2)\right]f^{(2)}(x_1,x_2;t)`$ (A1)
$`=`$ $`{\displaystyle 𝑑x_3\left[\overline{T}_{}(x_1,x_3)+\overline{T}_{}(x_2,x_3)\right]f^{(3)}(x_1,x_2,x_3;t)}`$ (A2)
where the collision operator is

$`\overline{T}_{}(x_1,x_2)`$ $`=`$ $`\sigma _0^2{\displaystyle _\mathrm{\Omega }}𝑑\widehat{\sigma }\mathrm{\Theta }\left(𝐩_{12}\widehat{\sigma }\right)\left(𝐩_{12}\widehat{\sigma }\right)\left[\delta \left(𝐪_{12}\sigma _0\widehat{\sigma }\right)\alpha ^2\widehat{b}_{12}^1\delta \left(𝐪_{12}+\sigma _0\widehat{\sigma }\right)\right]`$ (A3)
$`=`$ $`{\displaystyle 𝑑\sigma \delta \left(\sigma \sigma _0\right)\mathrm{\Theta }\left(𝐩_{12}\widehat{\sigma }\right)\left(𝐩_{12}\widehat{\sigma }\right)\left[\delta \left(𝐪_{12}\sigma \right)\alpha ^2\widehat{b}_{12}^1\delta \left(𝐪_{12}+\sigma \right)\right]}`$ (A4)
$`=`$ $`\delta \left(q_{12}\sigma _0\right)\left(𝐩_{12}𝐪_{12}\right)\left[\mathrm{\Theta }\left(𝐩_{12}𝐪_{12}\right)\alpha ^2\widehat{b}_{12}^1+\mathrm{\Theta }\left(𝐩_{12}𝐪_{12}\right)\right]`$ (A5)
$`=`$ $`\delta \left(q_{12}\sigma _0\right)\left(𝐩_{12}𝐪_{12}\right)\left[\alpha ^2\widehat{b}_{12}^1+1\right]\mathrm{\Theta }\left(𝐩_{12}𝐪_{12}\right)`$ (A6)
We now observe that no matter what the state, the atoms cannot overlap so we must be able to write the distribution as

$$f^{(2)}(x_1,x_2;t)=\mathrm{\Theta }\left(q_{12}\sigma _0\right)\stackrel{~}{f}^{(2)}(x_1,x_2;t)$$
(A7)
so that

$`\left[𝐩_1{\displaystyle \frac{}{𝐪_1}}+𝐩_2{\displaystyle \frac{}{𝐪_2}}\right]\mathrm{\Theta }\left(q_{12}\sigma _0\right)\stackrel{~}{f}^{(2)}(x_1,x_2;t)`$ (A8)
$`=`$ $`𝐩_{12}𝐪_{12}\delta \left(q_{12}\sigma _0\right)\stackrel{~}{f}^{(2)}(x_1,x_2;t)+\mathrm{\Theta }\left(q_{12}\sigma _0\right)\left[𝐩_1{\displaystyle \frac{}{𝐪_1}}+𝐩_2{\displaystyle \frac{}{𝐪_2}}\right]\stackrel{~}{f}^{(2)}(x_1,x_2;t)`$ (A9)
so that there are two sources of singularities in the second BBGKY equation: one coming from the collision operator and one from the streaming operator. If we integrate $`q_{12}`$ over a vanishingly small interval centered at $`\sigma _0`$, only the singular terms will contribute so that we conclude they must cancel independently of the regular terms. This gives

$$𝐩_{12}𝐪_{12}\delta \left(q_{12}\sigma _0\right)\stackrel{~}{f}^{(2)}(x_1,x_2;t)=\delta \left(q_{12}\sigma _0\right)\left(𝐩_{12}𝐪_{12}\right)\left[\alpha ^2\widehat{b}_{12}^1+1\right]\mathrm{\Theta }\left(𝐩_{12}𝐪_{12}\right)\stackrel{~}{f}^{(2)}(x_1,x_2;t)$$
(A10)
or

$$\delta \left(q_{12}\sigma _0\right)\mathrm{\Theta }\left(𝐩_{12}𝐪_{12}\right)\stackrel{~}{f}^{(2)}(x_1,x_2;t)=\delta \left(q_{12}\sigma _0\right)\alpha ^2\widehat{b}_{12}^1\mathrm{\Theta }\left(𝐩_{12}𝐪_{12}\right)\stackrel{~}{f}^{(2)}(x_1,x_2;t)$$
(A11)
which is the desired result relating the post-collisional distribution, on the left, to the pre-collisional distribution on the right. It is thus apparent that this identity carries part of the information of the second BBGKY hierarchy.

### 2 From Conservation of probability

It is clear that the probability to find two atoms moving towards each other with a given relative momentum, the left hand side, is the same as that to find two atoms in contact moving away from one another with the corresponding post-collision momentum:

$`\delta \left(𝐪_{12}\sigma \widehat{\sigma }\right)\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}\right)f_2(x_1,x_2,t)d^3x_1d^3x_2=\delta \left(𝐪_{12}\sigma \widehat{\sigma }\right)\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}^{}\right)f_2(x_1^{},x_2^{},t)d^3x_1^{}d^3x_2^{}`$

which simply states that the probability to find two atoms in contact moving towards each other with a given relative momentum, the left hand side, is the same as that to find two atoms in contact moving away from one another with the corresponding pre-collision momentum (the right hand side). Using

$`d^3x_1d^3x_2`$ $`=`$ $`{\displaystyle \frac{1}{\alpha ^2}}d^3x_1^{}d^3x_2^{}`$ (A12)
$`f_2(x_1,x_2,t)`$ $`=`$ $`\widehat{b}_{12}^1f_2(x_1^{},x_2^{},t)`$ (A13)
this can be written as

$$\delta \left(𝐪_{12}\sigma \widehat{\sigma }\right)\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}^{}\right)\frac{1}{\alpha ^2}\widehat{𝐛}_{12}^1f_2(x_1^{},x_2^{},t)=\delta \left(𝐪_{12}\sigma \widehat{\sigma }\right)\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}^{}\right)f_2(x_1^{},x_2^{},t)$$
(A14)
. The factors of $`\alpha `$ in the Jacobian arise because of the change from pre- to post-collisional momenta ($`𝐩_{12}^{}=\alpha 𝐩_{12}`$ gives one factor of $`\alpha `$) and because the relation between positions $`𝐪(t)`$ after the collision and positions before the collision also involves the momenta giving a second factor of $`\alpha `$.

### 3 Evaluating correlations at contact

From the identity and the assumption of molecular chaos, we have

$`\delta \left(𝐪_{12}\sigma \widehat{\sigma }\right)Vf_2(x_1,x_2,t)`$ $`=`$ $`\delta \left(𝐪_{12}\sigma \widehat{\sigma }\right)f_1(𝐱_1;t)f_1(𝐱_2;t)g_0(𝐪_1,𝐪_2;t)`$ (A16)
$`+\delta \left(𝐪_{12}\sigma \widehat{\sigma }\right)\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}\right)\left(\alpha ^2\widehat{b}_{12}^11\right)f_1(x_1;t)f_1(x_2;t)g_0(𝐪_1,𝐪_2;t)`$
so that if averages over $`f_1(𝐱_1;t)f_1(𝐱_2;t)g_0(𝐪_1,𝐪_2;t)`$ are denoted by $`\mathrm{}_0`$, then for any two-body function $`A=_{j<i}A(x_i,x_j)\delta \left(𝐪_{ij}\sigma \widehat{\sigma }\right)`$ one finds

$`A(\widehat{\sigma })`$ $`=`$ $`A`$ (A17)
$`=`$ $`A_0+{\displaystyle \frac{N(N1)}{2V}}{\displaystyle 𝑑x_1𝑑x_2\delta \left(𝐪_{12}\sigma \widehat{\sigma }\right)A(x_1,x_2)\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}\right)\left(\alpha ^2\widehat{b}_{12}^11\right)f_1(x_1;t)f_1(x_2;t)g_0(𝐪_1,𝐪_2;t)}`$ (A18)
$`=`$ $`A_0+{\displaystyle \frac{N(N1)}{2V}}{\displaystyle 𝑑x_1𝑑x_2\delta \left(𝐪_{12}\sigma \widehat{\sigma }\right)f_1(x_1;t)f_1(x_2;t)g_0(𝐪_1,𝐪_2;t)\left(\alpha ^1\widehat{b}_{12}1\right)\mathrm{\Theta }\left(\widehat{𝐪}_{12}𝐩_{12}\right)A(x_1,x_2)}`$ (A19)
$`=`$ $`A_0+{\displaystyle \underset{j<i}{}}\delta \left(𝐪_{ij}\sigma \widehat{\sigma }\right)\left(\mathrm{\Theta }\left(\widehat{𝐪}_{ij}𝐩_{ij}\right)\alpha ^1\widehat{b}_{ij}\mathrm{\Theta }\left(\widehat{𝐪}_{ij}𝐩_{ij}\right)\right)A(x_i,x_j)_0`$ (A20)
and, in particular, if $`A(x_i,x_j)=\frac{1}{4\pi n^2}`$ then

$`g(\widehat{\sigma })`$ $`=`$ $`g_0(\widehat{\sigma })+{\displaystyle \frac{1\alpha }{2\alpha }}g_0(\widehat{\sigma })`$ (A21)
$`=`$ $`{\displaystyle \frac{1+\alpha }{2\alpha }}g_0(\widehat{\sigma })`$ (A22)
as reported in the text.