Pair correlation function characteristics of nearly jammed disordered and ordered hard-sphere packings

We study the approach to jamming in hard-sphere packings and, in particular, the pair correlation function $g_2\left(r\right)$ around contact, both theoretically and computationally. Our computational data unambiguously separate the narrowing $\delta$-function contribution to $g_2$ due to emerging interparticle contacts from the background contribution due to near contacts. The data also show with unprecedented accuracy that disordered hard-sphere packings are strictly isostatic: i.e., the number of exact contacts in the jamming limit is exactly equal to the number of degrees of freedom, once rattlers are removed. For such isostatic packings, we derive a theoretical connection between the probability distribution of interparticle forces $P_f\left(f\right)$, which we measure computationally, and the contact contribution to $g_2$. We verify this relation for computationally generated isostatic packings that are representative of the maximally random jammed state. We clearly observe a maximum in $P_f$ and a nonzero probability of zero force, shedding light on long-standing questions in the granular-media literature. We computationally observe an unusual power-law divergence in the near-contact contribution to $g_2$, persistent even in the jamming limit, with exponent $-0.4$ clearly distinguishable from previously proposed inverse-square-root divergence. Additionally, we present high-quality numerical data on the two discontinuities in the split-second peak of $g_2$ and use a shared-neighbor analysis of the graph representing the contact network to study the local particle clusters responsible for the peculiar features. Finally, we present the computational data on the contact contribution to $g_2$ for vacancy-diluted fcc crystal packings and also investigate partially crystallized packings along the transition from maximally disordered to fully ordered packings. We find that the contact network remains isostatic even when ordering is present. Unlike previous studies, we find that ordering has a significant impact on the shape of $P_f$ for small forces.

Figure 1

(Color online) The polytope of allowed displacements, ${\mathcal{P}}_{\mathbf{\Delta }\mathbf{R}}$, for a locally jammed disk (light shade) trapped among three (left) or six (right, as in the triangular lattice) fixed disks (dark shade). The exclusion disks (dashed lines) of diameter twice the disk diameter are drawn around each of the fixed disks, along with their tangents (solid lines) and the polytope ${\mathcal{P}}_{\mathbf{\Delta }\mathbf{R}}$ they bound (dark). For the isostatic case on the left this polytope is a triangle (a simplex in two dimensions) and a hexagon for the hyperstatic case on the right.

Figure 2

A highly schematic plot of the subspace in the density-disorder $\left(\varphi \text{−}\psi \right)$ plane where strictly jammed three-dimensional packings exist. Point $A$ corresponds to the lowest-density jammed packing, and it is intuitive to expect that a certain ordering will be needed to produce low-density jammed packings. Point $B$ corresponds to the most dense jammed packing. Point MRJ represents the maximally random jammed state. This is the most disordered jammed packing in the given jammed category (locally, collectively or strictly jammed). We conjecture that the Lubachevsky-Stillinger packing algorithm [21, 22] typically produces packings along the right (maximally dense) branch, and we do not know of an algorithm that produces packings along the left (minimally dense) branch.

Figure 3

(Color online) The inverse of the “instantaneous” (averaged over several hundred collisions per particle) pressure of a nearly jammed (isostatic) packing of 1000 particles, as it is slowly diluted (using a negative expansion rate for the particles in the molecular dynamics algorithm $\gamma =-{10}^{-5}$) from ${\varphi }_J\approx 0.627$ until an unjamming particle rearrangement occurs. Up to this occurrence, the free-volume theoretical relation $p={\delta }^{-1}$ is satisfied to very high accuracy. There is a short transient region during the initial equilibration of the packing. Rattlers have been removed from the packing.

Figure 4

(Color online) The coefficient $C$ during a typical slow densification (expansion rate is ${10}^{-4}$) of a 10 000–particle system, starting from an equilibrated liquid at $\varphi =0.5$ up to jamming. The final packing has 259 rattlers, so the expected coefficient is $3\times 0.9741\approx 2.92$, a value which is shown with a red line. It is clear that close to the jamming point Eq. (4) is very accurate, but a marked lowering from a coefficient of 3 is seen for pressures lower than about ${10}^6$, likely explaining the coefficient 2.67 reported in works of Speedy [24, 25]. The inset shows the estimated “collisional” coordination, defined as the average number of different particles that a particle has collided with during a time interval of about 100 collisions per particle, during the same densification. The expected number $6\times 0.9741\approx 5.85$ is shown (this number is not asymptotically reached exactly since some of the $M$ contacts do not participate in collisions frequently enough to be registered during the time interval used), and we see that as many as eight contacts per particle are active at sufficiently large jamming gaps.

Figure 5

(Color online) The short-term (“instantaneous”) pressure versus number of events (mostly binary collisions) processed by the molecular dynamics algorithm [12], corresponding to a total run of about half a million collisions per particle. For the 1000-particle packing the pressure is stable, but for the larger packings a systematic pressure leak is observed.

Figure 6

(Color online) The cumulative coordination $Z\left(l\right)$ [i.e., the integral of $g_2\left(l\right)$] as a function of the gap tolerance $l$, for a sequence of snapshots of a 1000-particle packing during the final compression stages of the LS algorithm. Each snapshot is shown with a separate solid curve and only the last one is labeled in the figure legend. For a sufficiently slow expansion (expansion rate is ${10}^{-5}$ times the average thermal velocity), the packing is clearly seen to jam in an isostatic configuration. A subisostatic configuration is found for fast expansion (expansion rate is comparable to the thermal velocity). The inset shows the properly normalized derivative of $Z\left(l\right)$, right around contact, along with a comparison to our semitheoretical prediction for $g_2^{\left(\delta \right)}\left(l\right)$, for a packing with $\delta =2.5\times {10}^{-12}$.

Figure 7

(Color online) Computational data on the interparticle force distribution along with the best fit we could achieve. Packings of both 1000 and 10 000 particles, using either molecular dynamics to average the collisional forces or inversion of the rigidity matrix, were used, consistently producing the same probability distribution. Comparison to other data in the granular-media literature is beyond the scope of this work.

Figure 8

(Color online) The near-contact $Z^{\left(b\right)}\left(l\right)$ for a nearly jammed 10 000–particle packing, along with a power-law fit for small gaps, shown in both a linear-linear scale and a log-log scale (inset). In this inset we also show a line with slope 0.5 (i.e., a square-root dependence), which is clearly inconsistent with the numerical data.

Figure 9

(Color online) Computational data on the split-second peak of $g_2\left(r\right)$ averaged over five packings of 10 000 particles. The values $r=\sqrt{3}D$ and $r=2D$ are highlighted and match the two observed discontinuities. Also visible is the divergence near contact. The inset shows the probability distribution $P_{\theta }\left(\theta \right)$ of bond-pair angles in the contact network of the packings, also revealing two divergences at $\theta =\pi ∕3$ and at $\theta =2\pi ∕3$. No peaks are observed at $r=\sqrt{2}D$ or $r=\sqrt{5}D$, which are typical of crystal packings, indicating that there is no detectable crystal ordering in the packing.

Figure 10

(Color online) The probability distribution of local contact numbers as the cutoff used in defining neighbors is increased. Rattlers are clearly seen, and a relative maximum at $Z=6$ is seen. Note that only one particle with 11 neighbors is observed, and very few have as many as 10 neighbors. No particle with 12 contacting neighbors has been observed in any of our packings, indicating a lack of crystallinity.

Figure 11

(Color online) The first shell $g_2^{\left(\delta \right)}\left(l\right)$ for a collection of fcc crystal packings with a fraction $p$ of the spheres removed, starting with $N=13500$ particles. The inset shows the packing with most vacancies, where every fourth sphere is removed to form a cubic sublattice of vacancies (colored dark). Intermediate $p$’s are achieved by randomly adding back some of the spheres to the sublattice. The density has been reduced by $\delta =\sqrt{2}\bullet {10}^{-11}$ from close packing.

Figure 12

(Color online) The force probability distribution for the collection of fcc crystal packings from Fig. 11. For the pure crystal and the crystal with the most vacancies, all of the particle pairs are identical and therefore the probability distribution would be a $\delta$ function if forces are averaged over an infinite time horizon. For the intermediate $p$’s, multiple relatively broad peaks are observed. In contrast with the disordered case, very small forces are not observed.

Figure 13

(Color online) Compression of an initially liquid system with $\varphi =0.5$ to jamming with several different expansion rates $\gamma$ (the mean thermal velocity is 1, in comparison). The pressure is plotted on a reciprocal scale (the tickmarks being equally spaced in equal increments of $p^{-1}$, increasing in the usual direction), to highlight the expected linear relation (4) near jamming. The pressure-density curves for the perfect fcc crystal [23], the accepted fluid-solid coexistence region, and the widely known Carnahan-Starling equation of state for the fluid branch are also shown for comparison. Sufficiently fast compression suppresses crystallization and leads to densities around 0.64–0.65, and slower compression allows for partial crystallization, typically occurring around $\varphi \approx 0.55$, which is the end of the coexistence region (i.e., the density where the crystal necessarily becomes thermodynamically favored). This produces denser packings which exhibit more crystal ordering the denser they are.

Figure 14

(Color online) The evolution of the peaks in $g_2\left(r\right)$ as crystalline order is increased, for the packings from Fig. 13. The formation of peaks at distances typical of the fcc lattice, such as $r=\sqrt{2}$, is clearly seen. It is interesting to note that a peak is observed at $\sqrt{11∕3}\approx 1.91$, which is a fifth-neighbor distance in the hcp (but not the fcc) lattice (a similar hcp peak at $\sqrt{8∕3}\approx 1.63$ is barely visible). This is in agreement with numerous previous theoretical and numerical investigations of crystallization [33].

Figure 15

(Color online) Compression of an initially (metastable) liquid system with $\varphi =0.6$ to jamming with several different expansion rates, as in Fig. 13. For this range of expansion rates, crystallization is suppressed due to the large initial density and all final packings are apparently disordered and would be ordinarily identified as random; however, it is clear that slower compression leads to higher densities, and thus the final packings are not all identical, but rather some are more ordered than others, as can be verified by the slight increase in bond-orientational order metric $Q_6$ [18], for example.

Figure 16

(Color online) The evolution of $P_f\left(f\right)$ as crystalline order is increased, showing the disappearance of the peak at small forces. The inset shows a log-log view of the plot and is consistent with exponential decay for large forces.