We treat the case of an undriven gas of inelastic hard-spheres with short-ranged attractive potentials via an extension of the pseudo-Liouville operator formalism. New evolution equations for the granular temperature and coordination number are obtained. The granular temperature exhibits deviation from both Haff’s law and the case of long-ranged potentials. We verify this departure using soft-sphere discrete element method simulations. Excellent agreement is found for the duration of the simulation even beyond where exclusively binary collisions are expected. Simulations show the emergence of strong spatial-velocity correlations on the length scale of the last peak in the pair-correlation function but do not show strong correlations beyond this length scale. We argue that molecular chaos may remain an adequate approximation if the system is modelled as a Smoluchowski type equation with aggregation and break-up processes.

## I. INTRODUCTION

The use of micron-sized particles is ubiquitous in industry, and the properties of particles at this scale are of interest in a host of applications: fluidization of Geldart C powders,^{1,2} filtration,^{3} materials processing,^{4} powder mixing,^{5} rheology,^{6–9} etc. These particles also present themselves in extra-terrestrial problems, such as aggregation of planetesimals in the early solar system^{10} and contamination by Martian and Lunar regolith.^{11} These particles are interesting not only due to their ubiquity but also due to their unusual behavior. Energy is not conserved in collisions due to many internal degrees of freedom, and yet reversible short-range interactions remain important. The combination of these two interactions leads to particle trapping/sticking in binary collisions. That is, for collisions below a critical impact velocity, $ v crit $ (absolute value of the normal component of the relative velocity), a pair of particles does not retain adequate kinetic energy after a collision to re-emerge from a short-ranged potential well. They appear to stick after a binary collision, an impermissible result in both classical conservative and granular systems. The presence of this sticking mechanism qualitatively changes the behavior of many-particle systems.

Few kinetic theory approaches have been proposed and adequately developed for handling systems with both dissipation and attractive wells, where sticking behavior may present itself. One approach has been to focus on systems where sticking does not occur^{12} and treat weaker long-range interactions. An alternative approach is to treat the sticking systems with aggregation and break-up processes.^{2,13} This approach, though attractive, has two short-comings that arise from approximating the collision integral, which we stress is unknown *a priori*: $ 1 $ the dissipation from multi-particle collisions has been found to behave quite differently from binary collisions^{14} and $ 2 $ the evolution of the scattering cross-section is not known, i.e., aggregates develop unknown intraparticle structure. Often this lack of information is handled by treating systems in the limit where *all* particles stick to one another, producing larger spheres, the so-called ballistic coalescence.^{15–17} Recently, such theories have been applied to assessing the long-time cooling behavior of granular gases.^{18}

We will treat systems where sticking behavior is present using an extension of the pseudo-Liouville approach^{19} and subsequently explore the validity of such a theory. In an effort to understand these systems from a first principles approach, we study the simple problem of a gas of cohesive particles cooling via inelastic collisions. The problem of a homogeneous gas of inelastic hard-spheres was first considered by Haff.^{20} Through a heuristic treatment, Haff predicted the evolution of energy in a gas of such particles, known today as Haff’s law. More thorough treatments followed, employing the pseudo-Liouville operator theory^{21} and Sonine polynomial expansions.^{22} Other studies identified a number of instabilities and interesting behavior^{23,24} such as the growth of vortical structures and two-point density correlations, i.e., clustering. We expect that the nature of such instabilities is qualitatively changed in the presence of short-range attractive potentials, and we show evidence of interesting consequences of sticking by examining two-point correlations.

We solve the mean-field temperature evolution problem with the addition of short-range attractive potentials, building on the long-range results of Müller and Luding^{12} for dilute systems. The main difference between the approach developed in this work and that of Müller and Luding^{12} is the form of the final ordinary differential equation (ODE) governing the decay of kinetic energy. The discrepancy in the temperature evolution equations arises from differences in treatment of the collision dynamics: a single collision event^{12} vs. an infinite series of recollision events. Müller and Luding^{12} show that long-range attractive potentials accelerate cooling by effectively increasing the range of interaction, and hence, the number of particles with which a given particle will collide. The amount of dissipation in a collision is also affected. Corrections to the cooling law to account for higher density are also proposed by the authors.^{12}

For particles with short-range potentials, i.e., range of interaction much smaller than the mean free path and with collision durations much shorter than the free flight time, the treatment can be simplified. Collisions are again considered as discrete events, but the effective coefficient of restitution must be modified. A full population balance theory (presented in Appendix D) includes unclosed monomer–cluster and cluster–cluster terms in the collision integral. Our approach is thus able to treat the decay due to particle re-collision and aggregation. The theory may be simplified in the short-time limit, where the number of sticking collision events is small. A cooling law is also established in this limit by solving the final differential equation numerically, which we compare with Haff’s law and the results of Müller and Luding.^{12} Finally, we compare our results with soft-sphere discrete element method (DEM) simulations. For short times, where the number of sticking collision events is small, we quantitatively predict departure from both aforementioned cooling laws.

## II. ANALYTIC METHODS

### A. The binary system with short-ranged attraction

We begin our analytic treatment with a focus on the mechanics of a binary inelastic hard-sphere system. Interactions within the system are completely defined by so-called collision rules. Later, these rules will be used to construct collision operators appropriate for cases when cohesion is present. The simplest and most frequently used collision model is found in Eq. (1): a collision between two inelastic hard-spheres with constant coefficient of restitution. The most general form of the model relating pre-collisional velocities and post-collisional velocities is often given as

where subscripts denote vector components of velocities *v* and positions *r*. Einstein notation is also used, where repeated indices indicate an inner product. A single superscript indicates a particle index and denotes the quantity associated with that particle in the laboratory frame. Two superscripts indicate a relative vector constructed from the inertial reference frame of the particle with the first superscript, that is, $ v l j k = v l k \u2212 v l j $. Finally, a unit vector is indicated by a caret. The sole parameter determining the outcome of a collision in Eq. (1) is the coefficient of restitution, given by *ϵ*. The physical limits that *ϵ* may take are $\u03f5\u2208 0 , 1 $.

In a similar vein, collision rules may also be constructed for cases with short-ranged attraction so that we may apply the pseudo-Liouville operator theory. Whereas the instantaneous nature of hard-sphere collisions guarantees that we may straightforwardly construct collision rules for use in the pseudo-Liouville theory, particles with attractive wells are afforded no such luxury. We impose a few restrictions to set a context such that the collision events still remain essentially discrete and binary. The first restriction is that the length-scale of the potential well, *l _{well}*, must be much smaller than the particle diameter

*D*. This is true for micron sized particles, where interactions such as the van der Waals force remain important, and $D/ l w e l l =O 100 $. When considering the granular gas, we will also restrict ourselves to cases where the mean-free path λ is sufficiently large $ \lambda > l w e l l $. Finally, we assume that the time-scale for the duration of a collision, i.e., the time it takes for a particle to enter the well, collide, and exit the well, is much smaller than the free-flight time. The two latter assumptions restrict our cases to systems of moderate dilution, depending on solids volume fraction, particle properties, and thermal speed of the granules. Together, the three restrictions allow us to treat systems with the additional short-ranged attraction as approximately hard-sphere systems.

The attractive potential we consider is the van der Waals potential for macroscopic spheres,

where *A* and *R _{s}* denote the effective Hamaker constant for the materials that comprise the particles and the effective radius of curvature of the particle surfaces, respectively. The effective radius of curvature takes the value

*R*=

_{s}*D*/4 for smooth spheres. The effective separation distance between two surfaces of identical particles

*x*is defined as $x= r \u2225 j k \u22122a+ d 0 $, with the potential at contact

*x*=

*d*

_{0}being defined as Φ

_{c}. Details about using this potential to construct collision laws may be found in Appendix B.

Three stages of the collision process are identified: approach, collision, and separation (see Appendix B for details). If all three stages of a collision between the cohesive particles are considered, one obtains an equation relating the normal relative velocity of a pair of particles before entering the well, $ v \u2225 , 1 j k $ to the normal relative velocity of a pair exiting the well $ v \u2225 , 4 j k $. We refer to the relationship as an effective coefficient of restitution, given by

where $ Ha p =4 \Phi c / m v \u2225 j k 2 $ and the critical velocity is $ v crit =\u2212 \u2212 4 \Phi c 1 \u2212 \u03f5 2 / m \u03f5 2 $. The dimensionless quantity *Ha*_{p} is the ratio of potential energy at contact to the kinetic energy present in the normal component of the relative velocity at the edge of the potential well. Note that *Ha*_{p} and *v _{crit}* are negative quantities. This asymptotic expression was previously derived by Dahneke

^{25}and experimentally verified for 1.27

*μ*m latex spheres impacting a fused silica surface. In the dilute granular gas problem, the probability of particles beginning within the active region of the well is small, due to the short-ranged nature of the interaction, so we do not consider collisions between initially separating pairs.

The effective coefficient of restitution has interesting behavior, as shown in Fig. 1 for different values of *ϵ*. There exist two qualitatively different solutions, depending on the normal relative velocity at the edge of the well, $ v \u2225 , 1 j k = v i , 1 j k r \u02c6 i j k $. First, if the magnitude of the initial velocity is greater than the critical velocity, particles only experience an attenuation of the effective coefficient of restitution *ϵ*_{eff}. Particle pairs with more kinetic energy in the normal relative degree of freedom or alternatively weaker potentials, i.e., low |*Ha*_{p}|, experience almost no difference between *ϵ*_{eff} and *ϵ*. As |*Ha*_{p}| is increased, the discrepancy between *ϵ*_{eff} and *ϵ* also increases and collisions become more dissipative. The angle of separation for the collision partners also decreases (see Figs. 7(a) and 7(b) in Appendix A). Eventually, the initial kinetic energy is low enough that qualitatively different behaviors emerge, where particles no longer have sufficient kinetic energy to exit the well. In this case, particles cannot reach the edge of well and instead fall back into contact, and the pair experiences a rapid succession of recollisions. This quickly dissipates the energy in the $ r \u02c6 \u2225 j k $ degree of freedom. The conceptual picture for this collision outcome differs from the completely inelastic case (see Figs. 7(c) and 7(d) in Appendix A). For all practical purposes, these particles have stuck to one another forming a dimer. Additionally, the tangential components of relative velocity are transformed into rotating degrees of freedom for the newly formed dimer (see Fig. 7(d)). This result is fully characterized by $ Ha crit,p =\u2212 \u03f5 2 / 1 \u2212 \u03f5 2 $, which is only a function of *ϵ*.

We submit that much more complicated models for binary particle systems have been considered and thoroughly treated in the past. As an example, we point to Brilliantov *et al.*,^{26} who considered the addition of the van der Waals force to collisions between smooth visco-elastic particles. They also accounted for the increase in cohesive force at contact due to a substantial increase in contact area, which is consistent with the Johnson, Kendall, and Roberts (JKR) theory.^{27} Although, the model we consider is far less sophisticated, it maintains many of the key features obtained by the sophisticated model.^{26} For example, the relationship between *ϵ*_{eff} and impact velocity at low impact velocities remains qualitatively similar. Both theories also display dissipation and particle size-dependent critical impact velocities, under which particles adhere to one another, as observed in experiments.^{28,29} The advantage of our model is that we can explicitly express *ϵ*_{eff} as a function of the single parameter *Ha*_{p}, which allows for a cleaner analytic treatment in Sec. II B and Appendix D.

### B. Solving the cohesive problem

A thorough treatment using the pseudo-Liouville operator theory on a system of inelastic spheres is contained in Appendix C. In order to apply the pseudo-Liouville theory, as mentioned in Sec. II A, we must be able to approximate dilute systems of dissipative hard-spheres with short-ranged attractive potentials as hard-sphere systems with an altered hard-core distance, $| r i , c j k |=2a+ l w e l l $. In the treatment that follows, we neglect the additional length scale *l _{well}*, since

*l*≪

_{well}*a*. In Appendix D, we have developed an extension to the pseudo-Liouville operator that is applicable to systems with short-ranged attraction. The complete theory forms a set of population balance equations, which accounts for the enhanced dissipation from restituting and sticking collisions, as well as the aggregation and break-up of aggregates. However, because we do not know the source and sink terms in the aggregate–monomer and aggregate–aggregate difference operators

*a priori*, we must restrict ourselves to a theory for cooling that only accounts for the loss of energy due to monomer–monomer collisions, Eq. (D8). For this theory to be valid, we must restrict ourselves to a time interval, when these monomer–monomer interactions dominate. Note that this restriction leaves open the possibility for the development of aggregation and break-up kernels from simulation data, to extend the applicability of the theory to later times.

The pseudo-Liouville operators that we will use for the cohesive case are given by Eqs. (D2) and (D3) for restituting and sticking collisions, respectively. For restituting collisions, the difference operator is identical to that of Eq. (C2) with *ϵ*_{eff} replacing *ϵ*. The sticking collisions use the difference operator Eq. (D8). Following the procedure described in Appendix C up to Eq. (C5), one obtains a one-dimensional integral equation for the evolution of the ensemble-averaged kinetic energy per particle in terms of collision frequency $\omega =16 \pi 1 / 2 a 2 n g c T g 1 / 2 t $ and granular temperature, *T _{g}*,

The limits for the first term represent the sink for all pre-collisional energy, which is the same for both the restituting and sticking pseudo-Liouville operators (see Appendix D). The second integral arises only from the pseudo-Liouville operator for restituting collisions when $ v \u2225 \u2264 v crit / 2 $. Evaluating the integrals gives the evolution of the ensemble-averaged kinetic energy per particle,

Further simplifying by substituting for *v _{crit}*, we obtain the final evolution equation for the ensemble-averaged kinetic energy per particle,

where *Ha* = 4Φ_{c}/*mT _{g}* is the ratio of cohesive energy to average kinetic energy contained in the normal relative velocity, which is defined in terms of the granular temperature. We remind the reader that $ Ha crit,p = Ha p v \u2225 j k = v crit =\u2212 \u03f5 2 / 1 \u2212 \u03f5 2 $. This evolution equation appears strikingly similar to that obtained by Haff for inelastic spheres (see Eq. (C7)), with the exception of the exponential factor multiplying

*ϵ*

^{2}, and the parameter

*v*. This exponential factor determines the apparent inelasticity of the system. In the limit of very high granular temperature $ T g \u226b v crit 2 $ in Eq. (6), for example, we obtain cooling behavior that is nearly unchanged from that of a non-cohesive system. As Fig. 2(a) shows, as cooling progresses the apparent inelasticity increases, until in the limit of $ T g \u226a v crit 2 $, completely inelastic cooling behavior is expected.

_{crit}Müller and Luding^{12} obtained a cooling law for a granular gas undergoing weak long-range interactions using the pseudo-Liouville operator formalism that also looks similar to Haff’s law,

and it is useful to compare and contrast the result in Eq. (6) with theirs. Before discussing the quantitative differences between the two theories, we will comment on the differences in the approaches. *The main difference between the short-range and long-range theories is that the former considers the possibility of multiple dissipation events, while the latter considers only a single dissipation event.* More specifically, the theory developed herein says that attractive potentials increase the apparent inelasticity of particles. At sufficiently high values of *Ha*_{p} for a given *ϵ*, particles have the ability to aggregate, thus dissipating all energy in the normal relative velocity. Müller says that the effect of long-range interactions, in addition to increasing the effective inelasticity, is to permit particles that are separating to be pulled into contact. The effect of the long-range interactions is such that the higher the cohesion the greater the access to more potential collision partners at higher separating velocities. Their theory does not try to explicitly account for repeated collisions but does indirectly account for them by capturing particles in regions of phase space with insufficient separating velocities. Indeed, full treatment of recollision events requires higher-order ring-type equations,^{30} as noted in Müller and Luding.^{12} The potentials we consider, with a short-range of interaction allowing for temporal and spatial scale separation, allow us to account for these repeated collision/aggregation events in the context of the mean-field theory.

With high attractive potentials, we find that the rate of cooling is doubled in the Müller theory, Eq. (7). This is consistent with the increased access to collision partners at contact. Eventually, access is given to all pairs populating the velocity phase space at contact and is not just restricted to approaching collision partners. Figure 2(b) shows a comparison of cooling rates between the two theories. In the non-cohesive limit, both theories are identical to Haff’s law, and there is no significant difference between the theories for low values of $ Ha \u226a 1 $. The largest differences we find are in the highly cohesive regime for both low and high values of *ϵ*. In all but one limiting case, Müller’s theory predicts higher cooling rates than our current theory. The greatest difference in the predicted cooling rate between the theories is found for highly attractive and nearly elastic particles, where our current theory predicts the gas will cool much faster. The reason for this is that the majority of particles will still aggregate according to our theory, while the Müller theory has an upper limit of cooling only twice as fast as Haff’s law. Finally, it is interesting to note that there is never a large discrepancy for moderately inelastic cases, i.e., *ϵ* ∼ 0.6 − 0.8, since *Ha*_{crit,p} ∼ − 1.

## III. SOFT SPHERE DISCRETE ELEMENT METHOD

Previous treatment of granular gases often made use of event-driven codes, marching in collisions rather than time. However, in the attractive case, the microscopic sticking mechanism necessarily leads to an ever increasing number of repeated collision events within the well, and hence, soft-sphere DEM must be used. Use of DEM does not come without its own disadvantages; introducing such disparate length scales makes these simulations exceedingly expensive. This study concerns several soft-sphere systems with a linear spring dash-pot contact model.^{31} The parameter *k* is the spring stiffness in this model. First, dilute nearly elastic cases with *ϵ* = 0.97 are considered for several different initial values of *Ha*. Denser and more inelastic systems are also considered: one system without cohesion and two with cohesion. The two cohesive systems both have *ϵ* = 0.7 but start from different initial *Ha* values. The non-cohesive systems are simulated with *ϵ* = 0 and 0.7. Using the molecular dynamics package LAMMPS,^{32} particles are initialized according to a Mat$ e \u0301 rn$ hard-core process^{33} and with Maxwellian velocities with zero net linear momentum. Particles are then allowed to make approximately 20 collisions before dissipation and/or cohesion are activated. Cohesive forces are calculated using the gradient of Eq. (2), with a constant cohesive force $ F v d w =A R s /6 d 0 2 $ while particles are in contact. Table I contains all system parameters necessary to perform the simulations.

Parameters . | Dilute . | Moderately dilute . |
---|---|---|

ϵ | 0.97 | 0.7; 0 |

$ k \u2217 =k/ \rho D T g $ | 5.0 × 10^{7}, 4.0 × 10^{8} | 7.7 × 10^{12} |

Ha_{0} | −0.049, −0.398 | −0.24, −2.4 |

D/d_{0} | 10^{3} | 10^{4} |

R/_{s}D | 0.25 | 0.25 |

$\nu =N\pi D 3 / 6 L 3 $ | 0.01 | 0.084 |

L/D | 30 | 50 |

N | 515 | 20 146 |

Parameters . | Dilute . | Moderately dilute . |
---|---|---|

ϵ | 0.97 | 0.7; 0 |

$ k \u2217 =k/ \rho D T g $ | 5.0 × 10^{7}, 4.0 × 10^{8} | 7.7 × 10^{12} |

Ha_{0} | −0.049, −0.398 | −0.24, −2.4 |

D/d_{0} | 10^{3} | 10^{4} |

R/_{s}D | 0.25 | 0.25 |

$\nu =N\pi D 3 / 6 L 3 $ | 0.01 | 0.084 |

L/D | 30 | 50 |

N | 515 | 20 146 |

## IV. RESULTS AND DISCUSSION

This section is composed of three subsections. The first subsection explores cohesive systems in the very dilute and nearly elastic limit. The second subsection considers denser systems with higher inelasticity. Finally, in the third subsection, the appropriateness of mean-field closures used in our theory is explored.

### A. Dilute nearly elastic systems

In the spirit of the analysis done by Müller and Luding^{12} for systems with long-range attractive interactions, we look at the cooling predictions made for systems that are nearly elastic *ϵ* = 0.97 and also dilute *ν* = 0.01. The DEM results for two systems with different levels of initial non-dimensionalized cohesion *Ha*_{0} are plotted alongside the predictions made for the present short-ranged theory, the long-ranged theory,^{12} and Haff’s law^{20} (see Figs. 3(a) and 3(c)). It is seen that the results obtained from DEM most closely match the evolution of the present theory. However, after some time, the DEM abruptly diverges from all theories presented here. We find that this behavior occurs as the system rapidly aggregates into small clusters, eventually forming a single cluster as the simulation proceeds.

The accuracy of the predictions is more easily ascertained by looking at the quality factor obtained by normalizing the temperature by Haff’s law, *Q* = *T _{g}*/

*T*

_{g,Haff}. It is clear from Figs. 3(b) and 3(d) that the theory presented herein offers the best agreement with results from DEM. The long-range theory over predicts the rate of cooling (note that this was predicted in Fig. 2(b)) up until the rapid clustering event occurs. The cooling law that we present remains valid for these systems for a finite time interval, the extent of which is useful to ascertain.

We discussed in Sec. II B that the cooling law is expected to fail when the dynamics of the cooling from monomer–monomer interactions no longer dominate the system. We can estimate this time by applying the pseudo-Liouville theory to obtain the rate for number of sticking collision events per partner. To obtain this, we simply need to use the sticking piece of the extended pseudo-Liouville operator (see Eq. (D3)) and the difference operator in Eq. (D4). This operator gives the loss of monomers per sticking collision, and applying the procedure described in Appendix C, we can obtain the evolution equation for the ensemble-averaged number of sticking collisions per particle, i.e., the average coordination number $ Z $,

The average coordination number differs form the collision frequency by only the additional exponential term, which determines the number of collisions that may stick. Solving this equation, which is coupled to Eq. (6), one can predict when the temperature predictions should fail. We find that the neighborhood where $ Z \u21921$, in Figs. 3(b) and 3(d), correlates well with the time at which the DEM rapidly deviates from the predictions of our theory. This seems to confirm that the reason for the departure is that the dynamics are no longer dominated by monomer–monomer collisions.

### B. Inelastic systems at higher solid volume fraction

In order to assess the model’s capability to predict cooling outside of the range of validity, we also consider the accuracy of predictions for higher densities, *ν* = 0.087, and inelasticity, *ϵ* = 0.97. As an illustration of what aggregation in these systems looks like, we show in Fig. 4(a) a visualization of aggregates at the end of a simulation for a case where *Ha*_{0} = − 0.24. There are several large aggregates present in the system, with as many as 95 particles in the largest aggregate. A large number of smaller particle aggregates, e.g., dimers, trimers, and tetramers, also exist throughout the domain.

Figure 4(c) shows a comparison of the DEM behavior with all three cooling laws. For reference, we have also included the non-cohesive cases, where *ϵ* = 0 and 0.7 (see Fig. 4(b)). Our theory again *quantitatively* predicts the departure in cooling behavior from Haff’s law at short-times. The error in our prediction at the end of the strongly cohesive case (*Ha*_{0} = − 2.4) remains under 15% after ∼7 collisions per particle. The moderately cohesive case (*Ha*_{0} = − 0.24) is under 2% after ∼12 collisions per particle, which is approximately the same error in the non-cohesive case with *ϵ* = 0.7. The error in the completely inelastic case *ϵ* = 0 is also ∼15%. Similar findings regarding the accuracy of temperature predictions in denser systems were reported for the long-range attraction study for comparable solid volume fractions.^{12}

The discrepancy between Müller and Luding^{12} theory and ours is most visible in the *Ha*_{0} = − 0.24 case, where they over predict the cooling behavior. At short times, our theory also agrees better with the case of stronger cohesion. The system, however, is not well predicted by either theory at the final time, which we attribute to the large amount of aggregates in the system. We expect that the theories would show larger disagreement between one another if highly inelastic or slightly inelastic particles were considered, as mentioned in Sec. II B. Ultimately, the better agreement between DEM and our theory in this regime means that accounting for sticking, as discussed in Sec. II A, is essential for the accurate prediction of temperature evolution in systems with dissipation and short-ranged attractive potentials.

We also look at the number of contacts in the system to determine how far outside the range of expected validity our predictions have been made (see Fig. 4(d)). Comparing the average coordination number per particle $ Z $ from DEM to the numerical solution to Eq. (8), we find good agreement for both moderate (*Ha*_{0} = − 0.24) and strong cohesion (*Ha*_{0} = − 2.4) at early times (see Fig. 4(d)). Both cases show disagreement shortly after the times where $ Z =1$: *t*/*τ _{Haff}* ≈ 0.2 for

*Ha*

_{0}= − 2.4 and

*t*/

*τ*≈ 1 for

_{Haff}*Ha*

_{0}= − 0.24. In the strongly cohesive case,

*Ha*

_{0}= − 2.4, the theory appears to over-predict the number of contacts formed after this time, while the theory under-predicts the formation of contacts for the moderately cohesive case. We do not yet have a firm explanation for this difference. However, we note that the current theory does not fully account for aggregation events. Inevitably, any mechanistic explanation relies on information about the packing of monomers as more and more aggregation events occur. Examining packing information is beyond the scope of this study.

In comparing the prediction of temperature and coordination number, we observe quite unexpected behavior. Looking at Fig. 4(d) and comparing with more dilute systems of Sec. IV A, one would expect that the temperature predictions would be quite poor. Yet, if we look at the moderately cohesive case, where by the end of the simulation $ Z \u22433$, we find that the temperature prediction is no worse than the predictions made where Haff’s law holds (see Fig. 4(b)). It does appear as though sensitivity of the temperature on developing aggregate structure is not very important at these early times for systems that are sufficiently large. We believe that for too small a box size, an initial cluster can grow to dominate the system dynamics more readily.

### C. Accuracy of mean-field closures

The low error in the temperature evolution is promising. However, it is still of interest to explore in some detail the evolution of pair statistics which provide closure for the mean-field kinetic theory. We first examine the pdf of normal relative velocities, $f v \u2225 j k $, plotted in Fig. 5. The pdfs displayed are instantaneous pair statistics calculated from *N*(*N* − 1) samples in a single realization. The initial pdf is approximately that of an i.i.d. Maxwellian with variance 2*T _{g}* consistent with the theory. We have also plotted the result obtained for i.i.d. Laplace distributed particles, i.e., double exponential distribution. Note that the distribution is given by Laplace $ v i j ; \u2009 T g \u223c1/2 T g exp \u2212 v i j / T g $. The distribution of normal relative velocities $ v \u2225 j k $ given i.i.d Laplace distributed particle velocities $ v i j $ is found to be $f z =1/2 1 + 2 z exp \u2212 2 z $, where $z= v \u2225 j k / 2 T g $. This is the well-known asymptotic similarity solution

^{34}of the Boltzmann equation governing tail statistics in a granular gas obtained via the Lorentz gas model.

The comparison between statistics in the fully inelastic (*ϵ* = 0) case and the strongly cohesive case (*Ha*_{0} = − 2.4) is interesting. We find that for the inelastic case, an exponential-like tail begins to form even at these relatively short times, see inset of Fig. 5. At low velocities, the pdf remains essentially Maxwellian, in agreement with previous studies. On the other hand, the cohesive case has different and surprising behaviors. The entire pdf seems to agree well with the i.i.d. Laplace distribution, even for very low velocities. The mechanism leading to this nice result seems rooted in the selective inelastic behavior and microscopic sticking instability. We note that at initial time $ v crit 2 / 2 T g =\u22121.12$ and at the final time $ v crit 2 / 2 T g =\u22123.6$. Every collision with a velocity between the critical value and zero sticks, resulting in a growing population of particle pairs where $ v \u2225 j k =0$. Particles colliding in the tail portions experience relatively little change due to cohesion. Finally, particles that have stuck do not separate to collide again as in the completely inelastic case, but instead remain with the transverse portion of their relative velocities intact in the rotational degrees of freedom.

The more interesting behavior to observe is the emergence of spatial pair correlations, which serves as a measure to determine the adequacy of the molecular chaos assumption and the appropriateness of using the Enskog correction factor *g _{c}*. For example, Fig. 6(a) displays the pair correlation function

*g*(

*r*) at the initial time. The relative velocity correlations conditional on separation between a pair of particles, $ v \alpha 2 j k | r $, in the longitudinal $ \alpha = \u2225 $ and transverse directions $ \alpha = \u22a5 $ are also plotted. Note that these are all from a single snapshot from a single realization. We have plotted the correlations at initial time for reference in Fig. 6(a). At initialization, the pair correlation function displays typical behavior, with a small peak due to excluded volume effects. The velocity correlations also appear negligible; the velocity near contact appears well-predicted by the molecular chaos result $ v i j k v l j k =2 T g $.

The completely inelastic case gains a large peak in *g*(*r*) near contact, indicating that the Enskog correction factor may not provide an accurate closure. Velocity correlations close to contact likewise appear to be developing with longitudinal correlations extending further than transverse correlations, i.e., developing bias towards grazing collisions. At large distances, correlations do not appear to be present. This behavior has been observed and explored before and is thus not surprising. However, the correlations in the inelastic limit provide a nice comparison with the cohesive cases. Note that at early times, appreciable correlations do not develop for the non-cohesive case where *ϵ* = 0.7, and so, we have not plotted them.

We first observe that in the cohesive case, fine structure appears in *g*(*r*) not just at contact but further out. These fine structures indicate the formation of coordination shells, which gives rough detail about the average size of aggregates. The stronger cohesion case, Fig. 6(e), shows coordination structure out to ∼4*D* indicating that aggregate structures of at least 8*D* exist in the simulation box. The coordination peaks are followed by a slight trough extending to ∼6*D* indicating the area surrounding aggregates is underpopulated.

The correlations in normal relative velocities also display fine structure in the cohesive cases (see Figs. 6(d) and 6(f)). One finds that troughs develop in the longitudinal component at the same location as the coordination shells. This indicates a high population of particles moving neither away nor towards one another at these locations. The peaks in the longitudinal component occur in between these shells. The transverse component has much smoother structure. We attribute this structure primarily to the rotation of aggregates.

The most important feature of the velocity correlations found in Figs. 6(c)–6(f) is that the velocity correlations do not seem to extend beyond the aggregate length scale inferred from *g*(*r*). This implies that molecular chaos may still be an adequate model if a mean-field type theory of aggregation and break-up as presented in Appendix D is constructed. This result is counter-intuitive; the sticking behavior of inelastic systems with cohesion can serve to preserve molecular chaos. The sticking result does not allow pairs with more highly correlated velocities to separate and undergo further collisions with other particles, i.e., bias towards grazing collisions should be subdued. Local correlation/order is trapped in new larger aggregate structures, which do not diffuse as is the case for the traditional inelastic hard-sphere gas^{35} or undergo viscous heating.^{24} As a final remark, the addition of short-ranged attraction will have more of an effect on the granular gas than to simply hasten the formation of instabilities, which is also discussed independently in a paper by Gonzalez, Thornton, and Luding.^{36} For example, it is well-known that at late times, the traditional granular gas shows a strong correlation between regions of high density and low temperature.^{24,37} Even relatively weak cohesion would change the behavior of these regions, collapsing and arresting them in large aggregate structures with tensile strength. Indeed, there are many questions yet to be addressed and satisfactorily understood where short-ranged interactions are present with far-reaching implications, particularly for powder processing and fluidized reactor technologies.

## V. CONCLUSION

In this paper, we have considered an extension to the freely cooling granular gas in the presence of short-ranged attractive potentials. The theory developed herein makes use of spatial and temporal scale separation to treat collisions as discrete events via an extension to the pseudo-Liouville operator formalism. The predicted cooling behavior is qualitatively different from that predicted for long-range potentials in other work^{12} and agrees well with DEM over a relatively short simulation time, where the dynamics are dominated by monomer–monomer interactions. In addition, new criteria determining the validity of the mean-field theory are given based on the number of contacts formed, though the cooling predictions appear to surpass expectations for the larger boxes studied. We also note some differences in the formation and evolution of clusters and velocity correlations that arise due to sticking.

Finally, the behavior of pair statistics determining the applicability of the mean-field theory is explored. We find that spatial-velocity correlations extend about as far as density correlations at a given time. Molecular chaos appears to remain an adequate approximation for population balance models including aggregation and break-up processes. We form the foundation of this framework through the extended pseudo-Liouville operator theory. The theory developed here has implications for building transport models of highly sheared cohesive granular systems. Questions about the ultimate trajectory of these systems, aggregate structure, and aggregate collision outcomes remain and merit further study.

## Acknowledgments

We gratefully acknowledge the support provided for this work from NSF Grant No. CMMI 0927660 and partially from NSF Grant No. CBET 1134500.

### APPENDIX A: THE GEOMETRIC INTERPRETATION OF COLLISION LAWS

It is useful to consider the limiting behavior of the collision rule and how the post-collisional trajectory of the two-particle system with a short-ranged attractive potential differs from the completely inelastic case. For example, taking *ϵ* = 1, one recovers elastic hard-sphere behavior. If we employ the geometric interpretation of these collision laws plotted in the inertial frame of particle *j* (see Figs. 7(a)–7(c)), we find that *ϵ* serves to alter the angle of separation, the angle between the collision plane $ r \u02c6 \u22a5 j k $ and post-collisional relative velocity vector $ v i j k $. For *ϵ* = 1, this angle is equal to the angle of the approach. Lower values of *ϵ* attenuate this angle. For the limiting case of completely inelastic collisions with *ϵ* = 0, this angle is parallel to the collision plane $ r \u02c6 \u22a5 j k $. All energies in the relative velocity corresponding to the normal $ v \u02c6 \u2225 j k $ degree of freedom are dissipated. It is important to note that a completely inelastic collision does not imply that particles remain in post-collisional contact. This can only occur for the inelastic model if the tangential components of relative velocity, $ v \u02c6 \u22a5 j k $, are identically zero.

The effective coefficient of restitution discussed in Sec. II A has similar behavior as discussed above. Where particles do not stick, the attenuation of the angle of separation is the same as in the non-cohesive case. However, for sub-critical impact velocities, $ v \u2225 j k \u2208 v crit , 0 $ particles stick to one another, creating dimers (see Fig. 7(d)). The effect of adding short-ranged attraction causes particles to not only lose more energy but also to aggregate.

### APPENDIX B: CONSTRUCTING THE EFFECTIVE COEFFICIENT OF RESTITUTION

This section considers, in detail, how collision laws may be constructed for inelastic hard-sphere systems with an additional short-ranged attractive potential. The details that follow focus on the interaction contained in Eq. (2) but in principle may be straightforwardly extended to arbitrary short-ranged attractive potentials. Before forming the collision laws, we must establish some facts that are helpful in our construction. First, the van der Waals potential is not singular at contact but instead reaches a minimum at the inter-atomic separation *d*_{0}, i.e., the approximate diameter of the electron cloud of an atomic nucleus, typically ∼0.2 nm for most materials.^{38} Note that the interatomic length scale *d*_{0} is not physically meaningful for macroscopic particles, where such information has essentially been sacrificed due to coarse graining. Particles at this separation are in mechanical contact, and *d*_{0} simply serves as a parameter to set the attractive force between two particles at contact. This is the reason for setting the effective separation to *x* = *d*_{0} at contact, i.e., $ r \u2225 j k =2a$.

This model is applicable for particle diameters much larger than the length scale *d*_{0}. As a side note, because our particles are hard, i.e., rigid, we do not consider any increase in the attractive force at contact due to contact deformation, consistent with the theory of Derjaguin, Muller, and Toporov (DMT).^{39} We stress here that the analysis that follows is not restricted to the van der Waals interaction and is applicable to any generic short-ranged attractive potential.

The short-ranged attractive nature of these wells ensures that particles effectively do not feel one another at even relatively close separations. For example, two particles separated by 100*d*_{0} have 99% of the potential energy of two particles at infinite separation, $ \Phi c \u2212 \Phi 100 d 0 / \Phi c \u2212 \Phi \u221e =0.99$. Due to this fact, we are justified in using the approximation $\Phi l w e l l \u224a\Phi \u221e $. There is essentially no difference between a collision with two spheres initially separated by very large distances or initially separated by $O 10 nm $.

The outcome of a collision between two particles with short-ranged attractive potentials can now be determined. Collisions between cohesive spheres consists of three stages: $ 1 $ the approach where particles accelerate towards contact, $ 2 $ reflection and loss of energy due to the coefficient of restitution, and $ 3 $ slowing as particles separate. Equations (B1)–(B3) describe the collision process between two identical particles. The relationship between initial and final velocities is completely determined by a conserved Hamiltonian in stages one and three and by the coefficient of restitution in the second stage,

### APPENDIX C: THE PSEUDO-LIOUVILLE OPERATOR APPROACH WITHOUT SHORT-RANGED ATTRACTION

The backbone of the analytic approach for hard-sphere gases in this paper is the pseudo-Liouville operator.^{19,30} For inelastic hard-spheres, the pseudo-Liouville operator takes the form

The first term on the RHS is the streaming piece, which has no effect on the change in energy for statistically homogeneous problems. The second term, known as the collision operator, consists of two sifting terms: a Heaviside function Θ which works to select only pairs approaching one another and a Dirac delta function *δ* that selects pairs in contact to operate on, respectively. The last piece, $ b + j k \u2212 1 $ is a difference operator that replaces pre-collisional quantities with post-collisional quantities. One uses collision rules to derive the form of these difference operators. For example, using the collision rules in Eq. (1) for monodisperse particles, the difference operator for the change in kinetic energy from a binary collision between particles *j* and *k* takes the form

Here, the $ \u2032 $ superscript distinguishes a post-collisional value from a pre-collisional value. The main use of the pseudo-Liouville operator is to evolve functions of phase space variables, such a $ E kin \Gamma $, in time. By taking the ensemble average of the interaction portion of the pseudo-Liouville operator, we can obtain the cooling behavior,

The initial joint probability distribution $\rho \Gamma ; 0 $ for all particles is chosen to be i.i.d. Maxwellian in velocity space and uniform in position space. Note that the set of phase space variables is abbreviated as $\Gamma = { r i 1 , r i 2 , \u2026 , r i N ; v i 1 , v i 2 , \u2026 , v i N } $ and the differential *d*Γ is given as $d\Gamma = \u220f j N d v i j d r i j \u220f k < j \Theta r i j k \u2212 2 a $. We point the readers to previous works^{12,21} for the well-documented intermediate treatment, where the 6N dimensional integral is reduced to a 12 dimensional integral containing only the two-particle distribution. A few key assumptions must be made in order to perform this reduction, namely,

molecular chaos or the mean-field approximation, i.e., pre-collisional velocities, is uncorrelated;

the joint velocity distribution remains close to that of an equilibrium hard-sphere gas, which is implied by combining Eq. (C3), our initial distribution, and our previous assumption;

the collision velocities remain isotropic, implied by the same reasons as above;

excluded volume effects are the only spatial correlations present in the system and appear in the radial distribution function at contact, $ g c =g r 12 = 2 a $. We approximate the pair contact density using the model of Carnahan and Starling

^{40}$ g c = 1 \u2212 \nu / 2 / 1 \u2212 \nu 3 $ with volume fraction*ν*;the number of particles is sufficiently large that (

*N*− 1)/*V*≈*n*, where*N*,*n*, and*V*are the number of particles, number density, and volume of the system, respectively.

The resultant integral equation for the two-particle system is then found to be

where *T _{g}* is the granular temperature that is related to the average kinetic energy of a particle and the system via $ E kin , p = E kin /N=3m T g /2$. A change of variables is performed ($ r i 12 = r i 1 \u2212 r i 2 $, $ R i 1 = r i 1 $, $ v i = v i 1 \u2212 v i 2 / 2 $, and $ V i = v i 1 + v i 2 / 2 $) in order to simplify the final integro-differential equation. By integrating Eq. (C4) over all components except for the normal component of $ v \u2225 = v i r \u02c6 i 12 $, one obtains the following one-dimensional integral in terms of instantaneous collision frequency, $\omega =16 \pi 1 / 2 a 2 n g c T g 1 / 2 t $:

The problem has been reduced from a many-body problem to an ensemble of isolated particles interacting with a mean-field distribution. Hence, we expect that if correlations in two-point statistics remain small, this theory should provide a good model for the trajectory of the macroscopic system. Note that this is equivalent to saying that the molecular chaos approximation is adequate.

The ODE obtained by evaluating (C5) is

with solution

where the time scale *τ*_{Haff} is given by

The temperature evolution is commonly referred to as Haff’s law.

### APPENDIX D: THE EXTENDED PSEUDO-LIOUVILLE OPERATOR FOR SYSTEMS WITH SHORT-RANGE ATTRACTION

In Sec. II A, we observed that the nature of the hard-sphere inelastic system is qualitatively changed by the introduction of short-ranged attractive wells. Particles become sticky forming aggregates from low energy collisions. This causes problems when trying to straightforwardly apply the pseudo-Liouville operator,^{19,30} which was originally developed for colliding and separating particles. Instead, the operator must be altered to account for the fact that particles can stick and form aggregates. Here, we develop a Smoluchowski-type aggregation theory or birth-death process in the pseudo-Liouville operator framework, which will have the ability to account for annihilation of monomers and birth of dimers from collisions. Such a theory could be extended to larger size aggregates. However, the only aggregation kernel that we can compute explicitly using the collision laws is the source to the dimer phase from the monomer phase.

The first task in extending the pseudo-Liouville operator to the attractive problem is to split the operator into three pieces,

where the three terms on the right hand side of Eq. (D1) are (i) the free-streaming operator, (ii) a piece handling restituting collisions $i L + , r $, and (iii) one handling sticking collisions $i L + , s $, respectively.

The restituting operator is straightforwardly constructed as

The only part in Eq. (D2) that is altered is the argument for the Heaviside function, which controls the limits of integration. We alter notation here to distinguish between particle phases. For example, a superscript $ 1 p $ indicates the monomer phase and $ 2 p $ the dimer phase. The restituting operator maintains many of the same features of the pseudo-Liouville operator for hard-spheres. Both particle number and linear momentum are conserved. The source and sink term operating on the energy also take a familiar form given in Eq. (C2) with *ϵ* being replaced by *ϵ*_{eff}.

The second operator for sticking collisions takes a slightly different form,

We first notice that the limits of integration have again been changed so that only sticking collisions are operated upon. Next, we notice that the source and sink operators have changed here. The sink is denoted by $ b s j k , 1 p $, which accounts for the annihilation of quantities from the monomer phase. The source, given by $ b + , s j , 2 p $, serves as a source of some quantity to the dimer phase. The operator is non-zero only when operating on the denoted phase, e.g., the source term with a 2*p* superscript operating on the number of particles in the monomer phase is zero, $ b + , s j , 2 p N j k , 1 p =0$.

All that remains is to define these source and sink operators for the quantities of interest in this study, namely, phasic number and energy. The number equations for the monomer and dimer phases are

The number equation tells us that in a sticking collision, 2 monomers are annihilated to form a single dimer. Note that this is valid for monodisperse monomers in a statistically isotropic environment. Extension to account for both polydispersity and anisotropy is straight-forward but not as clean.

The energy equations for the monomer and dimer phases include the energies in the translational and rotational degrees of freedom for the dimer phase,

where the vectors given by $ w i j k , 1 p $ and $ R i j k , 1 p $ refer to the center of mass velocity and position of the collision partners. Next, all energies in the 6 degrees of freedom of the binary system are annihilated (Eq. (D5)). The energy in the center of mass degrees of freedom is straight-forwardly transformed into the energy in the translational degrees of freedom for the dimer (Eq. (D6)). The degrees of freedom in the tangential components of the relative velocity are transformed into rotational energy for the dimer (Eq. (D7)). Finally, the energy in the normal relative velocity degree of freedom is completely dissipated due to the enhanced dissipation from the attractive potential.

Despite having this nice extension, we cannot make adequate use of it without knowledge of how larger aggregates form and break-up. Indeed, this could only be meaningfully obtained through either experimental or, most likely, numerical means. So, we will restrict ourselves to the case where aggregation is ignored, i.e., there is no coupling to an *N*-particle phase. We account for energy dissipation only. In that case, the source and sink operators take the form