# Static Structural Signatures of Nearly Jammed Disordered and Ordered Hard-Sphere Packings: Direct Correlation Function

###### Abstract

The nonequilibrium process by which hard-particle systems may be compressed into disordered, jammed states has received much attention because of its wide utility in describing a broad class of amorphous materials. While dynamical signatures are known to precede jamming, the task of identifying static structural signatures indicating the onset of jamming have proven more elusive. The observation that compressing hard-particle packings towards jamming is accompanied by an anomalous suppression of density fluctuations (termed “hyperuniformity”) has paved the way for the analysis of jamming as an “inverted critical point” in which the direct correlation function , rather than the total correlation function diverges. We expand on the notion that provides both universal and protocol-specific information as packings approach jamming. By considering the degree and position of singularities (discontinuities in the -th derivative) as well as how they are changed by the convolutions found in the Ornstein-Zernike equation, we establish quantitative statements about the structure of with regards to singularities it inherits from . These relations provide a concrete means of identifying features that must be expressed in if one hopes to reproduce various details in the pair correlation function accurately and provide stringent tests on the associated numerics. We also analyze the evolution of systems of three-dimensional monodisperse hard-spheres of diameter as they approach ordered and disordered jammed configurations. For the latter, we use the Lubachevsky-Stillinger (LS) molecular dynamics and Torquato-Jiao (TJ) sequential linear programming algorithms, which both generate disordered packings, but can show perceptible structural differences. We identify a short-ranged scaling as that accompanies the formation of the delta function at that indicates the formation of contacts in all cases, and show that this scaling behavior is, in this case, a consequence of the growing long-rangedness in , e.g., as for disordered packings. At densities in the vicinity of the freezing density, we find striking qualitative differences in the structure factor as well as between TJ- and LS-generated configurations, including the early formation of a delta function at in the TJ algorithm’s packings, indicating the early formation of clusters of particles in near-contact. Both algorithms yield structure factors that tend towards zero in the low-wavenumber limit as jamming is approached. Correspondingly, we observe the expected power-law decay in for large , in agreement with previous theoretical work. Our work advances the notion that static signatures are exhibited by hard-particle packings as they approach jamming and underscores the utility of the direct correlation function as a sensitive means of monitoring for the appearance of an incipient rigid network.

###### pacs:

61.50.Ah, 05.20.Jj## I Introduction

Packings of hard particles in -dimensional Euclidean space have been used ubiquitously as a powerful model to describe many-body systems such as liquids, glasses, colloids, granular materials, particulate composites, and biological systems, among others Bernal (1965); Zallen (1983); Cates et al. (1999); Torquato and Stillinger (2001); Torquato et al. (2003); Torquato and Stillinger (2003); Zachary and Torquato (2009); O’Hern et al. (2003); Manoharan et al. (2003); Makse et al. (2004); Xu et al. (2005); Zohdi (2006); Majmudar et al. (2007); Gevertz and Torquato (2008); Karayiannis et al. (2009); Parisi and Zamponi (2010); Jiang et al. (2010); Xu et al. (2011); Charbonneau et al. (2012); Dagois-Bohy et al. (2012); Gillman et al. (2013); Klatt and Torquato (2014); Yamchi et al. (2015); Royer and Chaikin (2015); Schaller et al. (2015); Ramola and Chakraborty (2016). In three dimensions, the venerable hard-sphere model is particularly useful, owing to its mathematical simplicity and the rich diversity of equilibrium and non-equilibrium behavior that it exhibits.

It has been shown that bringing hard-particle packings towards jamming (roughly speaking, mechanical stability) is accompanied by an anomalous suppression of long-range density fluctuations Silbert and Silbert (2009); Zachary et al. (2011a, b); Hopkins et al. (2012); Chieco et al. (2016); Atkinson et al. ()—a phenomenon known as “hyperuniformity” Torquato and Stillinger (2003); Zachary and Torquato (2009). A many-particle system is hyperuniform if the structure factor (trivially related to the Fourier transform of pair statistics in direct space) tends to zero in the limit that the wavenumber tends to zero. Hyperuniformity may be conceptualized as an “inverted critical point” in which the direct correlation function , which is defined through the Ornstein-Zernike integral equation for a system with number density :

(1) |

becomes long-ranged, i.e., its volume integral diverges Torquato and Stillinger (2003). This is to be contrasted with the usual thermal critical point (e.g., liquid-vapor or Curie critical points) in which the total correlation function (rather than ) becomes long-ranged. Accordingly, a static length scale, obtained from the Fourier transform of , , grows as a system approaches a hyperuniform state and ultimately diverges at this critical state. In hard-particle packings, this occurs as jamming is approached Hopkins et al. (2012), and a similar analysis can be used to obtain meaningful information about the nature of glassy states of particles with soft interaction potentials Marcotte et al. (2013). Because of this, it has been an intriguing prospect to investigate disordered hyperuniform systems by adapting standard tools used to investigate critical phenomena.

At the same time, the direct correlation function has proven to be a fruitful starting point for efforts to characterize the structure of disordered systems such as simple liquids since the pair statistics may be obtained through Eq. (1) Hansen and McDonald (1990). Physically, the equation suggests that pair statistics between any two particles may be decomposed into a “direct” contribution encoded in as well as an “indirect” contribution mediated through “chains” of particles, expressed mathematically through the convolution between and . For systems with suitably well-behaved interactions, one may equivalently think of as describing the linear response of a system to a perturbation in an externally applied potential field Hansen and McDonald (1990).

It has been observed that maximally random jammed (MRJ) hard-sphere packings, which constitute the most disordered configurations as measured by some scalar order metric subject to the constraint of jamming and isostaticity Torquato et al. (2000); Torquato and Stillinger (2010), are hyperuniform and exhibit a nonanalytic linear behavior in the structure factor for low , namely, Donev et al. (2005a). Hopkins et al. Hopkins et al. (2012) studied the behavior of very large () sphere packings produced by the Lubachevsky-Stillinger event-driven molecular dynamics algorithm Skoge et al. (2006) under rapid compression so as to study the approach to the MRJ-like states at densities well above the freezing density and close to jamming. They found evidence that the nonanalytic linear behavior in was evident considerably in advance of jamming, and that upon further compression, the extrapolated value at the origin tended towards zero, implying that corresponding long-ranged behavior in might be observable for this protocol. Their computations for the Fourier-transformed direct correlation function showed extending towards negative infinity near the origin as the packings were compressed, supporting this prediction.

In the current work, we further develop the view that one may find static, structural precursors to jamming in hard-particle systems. Because is known to generally possess a qualitatively simpler functional form while still encoding the complete pair statistics of the system, we will focus primarily on the signatures therein, paying particular attention to features that point towards the development of an incipient contact network and hyperuniform density fluctuations (i.e. long-rangedness in ). Moreover, is known to possess various singularities at jamming (e.g., a Dirac delta function at contact, discontinuity at a distance of two diameters), and we determine to what extent these features are inherited by .

However, one should also bear in mind that different packing protocols will tend to produce different ensembles of disordered jammed states Jiao et al. (2011). In this paper, we will also bring attention to qualitative differences in protocols’ approach to their jammed states. Additionally, it is nontrivial to ensure that standard protocols approach properly-jammed states and avoid becoming stuck in unstable mechanical equilibria. This has been found to be related to a “critical slowing down” that becomes of practical concern for large systems Atkinson et al. (). Therefore, we carry out our current investigation as if our systems are indeed jammed and hyperuniform with the important caveat that this is more difficult to do with high precision in practice than previously thought. If one had a protocol that were able to produce better-jammed packings, we expect that the packings would possess stronger structural signatures consistent with hyperuniformity.

To this end, we analyze computer-generated packings of monodisperse hard-spheres of diameter created by the LS algorithm as well as the Torquato-Jiao (TJ) sequential linear programming algorithm Torquato and Jiao (2010). We consider these two algorithms since (i) they are both known to generate highly disordered packings under suitable conditions, but also because (ii) the jammed states they produce possess considerable differences in their macroscopic properties, including density, rattler fraction, and degree of order as measured by various standard order metrics Atkinson et al. (2013). By investigating multiple protocols that differ considerably, we seek to discern what features are protocol-dependent, and which are in common to a diversity of MRJ-like states. We also briefly consider the behavior of the hard-sphere FCC crystal, conjectured to be the equilibrium phase Mau and Huse (1999); Pronk and Frenkel (2003) at packing fractions along the solid branch ending at close packing with a packing fraction of . This case provides valuable information about jamming under an arguably more well-behaved setting, where one need not worry about metastability, and hyperuniformity may be approached to arbitrary numerical precision with minimal practical issues.

We observe that the TJ and LS protocols exhibit markedly different qualitative behavior in , even at packing fractions far from jamming. Specifically, we find that the direct correlation of packings produced by the TJ algorithm exhibit signs of a delta function at at packing fractions below the freezing density , whereas features in exhibited by LS for are substantially more subtle up until much higher densities. With the development of the expected delta function at , we observe a concomitant development of a dominant scaling for as . We show using theoretical arguments that one can be predict this numerical observation. Interestingly, we observe the power-law scaling in the limit that is a consequence of the linear trend in for small . Observing this behavior is difficult in practice because it requires an accurate measurement of for small wavenumbers, which requires that one consider large packings. Our work advances the notion that static signatures are exhibited by hard-particle packings as they approach jamming and underscores the utility of the direct correlation function as a sensitive means of monitoring for the appearance of an incipient rigid network.

The remainder of the paper is organized as follows: in Sec. II, we discuss some relevant analytical results pertaining to the structure of disordered hard-sphere packings through and the Ornstein-Zernike equation. In Sec. III, we discuss quantitatively the manner by which inherits singularities from . In Sec. IV, we review some known facts regarding the critical scaling behavior expected for systems approaching hyperuniformity. In Sec. V, we review the protocols that we use to generate nearly-jammed hard-sphere packings. In Sec. VI, we present the structure factor and direct correlation functions of our ordered and disordered packings as they approach jamming and point out emergent static structural features exhibited by each. Conclusions and discussion are presented in Sec. VII.

## Ii Ornstein-Zernike Equation and Jamming

We begin by reviewing a number of relationships between the standard pair statistical descriptors of point processes with the goal of relating the direct correlation function to other familiar statistical descriptors. We then proceed by reviewing theoretical progress that has been made in obtaining an accurate description of disordered hard-sphere systems by use of the direct correlation function, including shortfalls that persist with the current state of the art which point to the necessity for our present numerical investigations.

The structure factor is defined for a translationally-invariant system in as

(2) |

where is the Fourier transform of the total correlation function. This is related to the scattering intensity , defined for a single system of particles within a fundamental cell under periodic boundary conditions as

(3) |

which includes forward scattering, i.e., . This is to be contrasted with the definition of Eq. (2), in which is related to the volume integral of . Apart from at , the scattering intensity is identical to the structure factor for a single configuration. For an ensemble of periodic point configurations such (e.g. derived from the particles centers of our packings), the ensemble average of is directly related to the structure factor via

(4) |

where is the Dirac delta function and the limit being taken on and the fundamental cell volume are such that the relevant physical system (e.g. unjammed packings at some constant density or at some given distance to the jamming density ) is preserved. In practice, we directly compute from our simulation data and average the data from a large number of packings with the same system size to approximate for the ensemble of packings generated by a given protocol, keeping in mind that finite-system artifacts are expected to persist to some degree. In all cases, the quantity must be inferred through extrapolation to the origin. Using this ensemble average, one may combine the Fourier transform of Eq. (1) with Eq. (2) in order to express the Fourier transform of the direct correlation function as

(5) |

from which may be computed by inverse Fourier transform.

A number of closure relations have been proposed to offer approximate solutions to Eq. (1).
In particular, the Percus-Yevick (PY) closure demands that for and for for monodisperse spheres of unit diameter ^{1}^{1}1To be mathematically precise, the statistical functions , , , and are all functions which map ; however, since these functions are often isotropic, we often express them in terms of or as appropriate..
Physically, these requirements specify that no two spheres may overlap, and that direct interactions (in the sense of Eq. (1)) are absent beyond the particles’ hard cores.
We would like the former criterion to apply in any solution for hard-spheres, but the latter assumption is increasingly violated as the packing fraction increases.
The cubic polynomial form of produced by the PY closure that solves Eq. (1) accurately describes much of the equilibrium liquid branch of the hard-sphere system with good quantitative agreement for and qualitative agreement for .
However, it possesses various shortcomings at higher densities still within the liquid branch that one might like to improve: it underpredicts (and thus the pressure), and oscillations in the pair correlation function are out of phase and decay too slowly with increasing Torquato (2002).
At higher densities, fails basic satisfiability criteria such as nonnegativity and hence ceases to be physical.

In order to address this shortcoming, a variety of adjustments have been made to the PY approximation to improve the range of densities over which it may apply. A classic approach is to introduce a Yukawa term beyond the core Waisman (1973); Høye and Stell (1976); Høye et al. (1976), i.e., . This improves the degree to which the approximation matches qualitative features in and provides additional fitting parameters to allow for a quantitative match of additional system properties. To this end, the recent work of Jadrich and Schweizer Jadrich and Schweizer (2013) used a two-Yukawa generalized mean spherical approximation which allowed them to match the system’s compressibility as well as and its first derivative in an attempt to accurately describe the behavior of the hard-sphere system along some metastable branch leading towards a disordered jammed state. By allowing and to approach infinity, this model may capture the appearance of a delta function at , and predicts a functional form for inside the core that departs from the solution to the Percus-Yevick approximation. For a single-Yukawa form with , one finds Waisman (1973)

(6) | |||||

for , where , , , satisfies

(7) | |||||

(8) |

and . Interestingly, upon taking and towards infinity so as to obtain a delta function at , we see that the fourth term in Eq. (6) scales as for as , where . Thus, as is taken to infinity, goes to zero. For , the term saturates to a constant. Note that, in order for this scaling behavior to imply that diverges at the origin, one must have . This is not necessarily implied by this analysis foot:Coulomb ().

While this approach improves upon these features, it still lacks the correct long- scaling behavior for MRJ-like packings, i.e., as . Accordingly, the small- behavior of the structure factor is not in qualitative agreement either. In addition to this, several salient features in the pair correlation function, including (i) the power-law divergence as , (ii) the cusp at , and (iii) the correct step discontinuity at , remain elusive.

## Iii Inheritance of features in from

Here, we present a derivation for the magnitude of the step discontinuity in in terms of the information in for the specific case of disordered, jammed hard-sphere packings in three dimensions. Critically, we assume throughout this analysis that (i) and are isotropic (radial) functions and that (ii) the packings are isostatic, meaning that the average contact number for backbone spheres in the packing is , where is the number of backbone spheres in a packing with rattler fraction . It is estimated that for TJ and for LS Atkinson et al. (2013). The particular form of the term depends upon whether one is considering collective or strict jamming, but becomes irrelevant in the infinite-system limit Donev et al. (2005b). The pivotal observation from which this analysis follows is that a singularity of a given order in may not contribute to a singularity in of lower order. For example, a step discontinuity in cannot cause a delta function to appear in . This assumption is justified by recursively inserting the form of into the Ornstein-Zernike relation to obtain

(9) | |||||

and noting that each successive convolution ought to increase the order of the differentiability class of the term by one.

We begin with the delta function at . By Eq. (9), we see immediately that the only contribution to this must come from a delta function in , and that the two must be of equal magnitude. The number of particles that are separated at a distance is given by . We substitute , where is the backbone fraction, and obtain the result that the strength of the delta function in is equal to .

We now proceed to identify the source of the jump discontinuity found at . We begin by decomposing into three parts:

(10) |

where is the delta function contribution from above; , where is the Heaviside theta function, captures the step discontinuity predicted in at , and captures the remainder of the direct correlation and is assumed to be at least continuous at . Substituting this into Eq. (9) gives

(11) | |||||

Substituting in and into this and subtracting the two equations from each other while letting gives:

(12) |

where we have used the shorthand and defined

(13) |

We notice that , meaning that the only surviving term of the first three is . Furthermore, by considering the decomposition of Eq. (10) as applied to Eq. (13), we can see that the term with the differentiability class of the lowest degree is given by

(14) |

and that should be nonzero so long as the delta function represented by has a nonzero amplitude (i.e. is not trivially zero everywhere). All other terms contributing to are continuous at , as are all higher-order convolutions. Thus, we are left with

(15) |

To evaluate the first term, we define , so that

(16) |

We invoke the convolution theorem to write

where we have invoked translational invariance to remove and the symmetry of the two terms within the integrand with respect to . Carrying out the inverse Fourier transform gives the result

(17) | |||||

Because represents the “sharpest” contribution from the single convolution term, we can conclude immediately that no other terms within the first convolution term (and no further convolution terms) will contribute to the quantity since they will be too smooth (i.e. they do not retain a step discontinuity following convolution). Thus, we arrive at the result that the magnitude of the step discontinuity in the total correlation function (and, equivalently, the pair correlation function) at is

(18) |

By extending this analysis, we can claim that even in the absence of any additional nonanalytic features in , is expected to have discontinuities in successively higher derivatives at further integer values of . For example, there ought to be a discontinuity in the first derivative of at , a discontinuity in the second derivative at , and so on.

## Iv Scaling relations for systems in the vicinity of hyperuniformity

In this section, we recall various scaling behaviors for various pair statistics in direct and Fourier spaces for ordered and disordered packings of hard spheres in the vicinity of jamming that are particularly germane to this paper.

Torquato and Stillinger have shown Torquato and Stillinger (2003) that a hyperuniform system with a structure factor that scales as

(19) |

may be thought of as an “inverted critical point.” At this point, the direct correlation becomes long-ranged, scaling as

(20) |

in dimension , where is a critical exponent such that . Additionally, the inverse of the structure factor at the origin exhibits critical scaling behavior in the vicinity of its critical density, i.e.,

(21) |

for densities close to, but below Torquato and Stillinger (2003).

In the case of the equilibrium crystal branch, we may exploit the compressibility relation relating the structure factor to the isothermal compressibility at temperature

(22) |

along with the free-volume equation of state Salsburg and Wood (1962) which predicts that the pressure behaves as

(23) |

to obtain the result that in the vicinity of jamming. From this, it follows that

(24) |

Note that this result is independent of dimension. One may also define a correlation length with the critical behavior

(25) |

which may be related to the previous critical exponents through

(26) |

As mentioned before, one such length scale can be defined by the volume integral of the direct correlation function: .

For packings along glassy metastable branches leading to MRJ-like states, one cannot use the compressibility relation because the states are nonequilibrium in nature Hopkins et al. (2012). However, one may reconcile the differing pictures presented by and by introducing a “nonequilibrium index”, defined as

(27) |

Hopkins et al. studied the behavior of under rapid compression toward jamming for MRJ-like packings prepared by the LS algorithm for system sizes up to and found that, as increased toward jamming, Hopkins et al. (2012).

By combining this result with the observation that the pressure in MRJ-like packings again diverges according to the free-volume equation of state, we get the result that

(28) |

By Eq. (26), this implies that for MRJ packings as well. This is a noteworthy result because it tells us that while both MRJ-like and equilibrium crystalline packings become increasingly hyperuniform as they are compressed, they are different universality classes with respect to the critical exponents and associated with hyperuniformity. While ensembles of packings are known to exist that interpolate between these extremes, the interpolation is not unique, owing to the large diversity of jammed packings that are known to exist. Thus, it is an interesting, outstanding question how these critical exponents might evolve between these two extreme states.

## V Simulation Methods

In order to study the behavior of the equilibrium hard-sphere FCC crystal for densities between and the close-packing density , we used standard event-driven molecular dynamics Skoge et al. (2006). Configurations of spheres with were placed on their lattice sites and allowed to equilibrate at fixed packing fraction within a cubic fundamental cell with periodic boundary conditions for collisions per sphere before taking statistics. Measurements of the structure factor were made every collisions per sphere to verify that equilibrium had been attained.

To generate disordered, strictly jammed sphere packings in three dimensions, we begin with initial conditions produced by random sequential addition (RSA) at an initial packing fraction of . We use the Torquato-Jiao (TJ) sequential linear programming method Torquato and Jiao (2010) with system sizes of up to using the same parameters as those used to study the MRJ state in Ref. Atkinson et al. (2013). The final mean density of these packings is for a system size of . After the algorithm terminates and a putatively jammed state is reached, the packing is equilibrated within its jamming basin using event-driven molecular dynamics at fixed density. We also compare these packings to those generated using the well-known event-driven Lubachevsky-Stillinger (LS) molecular dynamics algorithm Lubachevsky and Stillinger (1990). For LS, we use an initial dimensionless growth rate of until packings reach a dimensionless pressure of , at which point the expansion rate is slowed to , and packing continues until . The mean density of the final packings as prepared under this protocol is at a system size of .

## Vi Results

In this section, we present results pertaining to our computer-generated hard-sphere packings as they approach both ordered and disordered jammed states. We will present our analysis assuming that the particle diameter is unity unless otherwise specified.

### vi.1 Fcc

We begin by examining the behavior of the hard-sphere FCC crystal because the behavior exhibited by the crystal undergoing thermal motion away from jamming serves as a interesting starting point from which we can make several observations to guide our subsequent investigation of disordered jammed systems. We will investigate the crystal’s approach to close-packing with attention to the fluctuating component of the structure factor as well as the implications it has for the qualitative form of .

Figure 1 shows plots of the radially-averaged structure factor for the equilibrated FCC crystal for a variety of densities along the solid branch; our computation follows the collective coordinate formulation of Eq. (3). Curves are averaged over ensembles of 100 packings with . Figure 2 shows the corresponding radially-averaged direct correlation function evaluated numerically using discrete Fourier transform techniques following Eq. (5).

In addition to the Bragg peaks that arise from the crystal geometry from the “frozen-in” structure of the packings, the curves possess a background contribution derived from thermal fluctuations which scales as starting at a sufficiently-high wavenumber (on the order of unity) and saturates to . While the domains of interest in the scaling descriptions of Eqs. (19) and (20) as written are and , one might reasonably suspect that it is possible to invert the two limits so as to infer the scaling behavior of from . In the case of our FCC data, we see that this implies that the direct correlation function ought to scale as for . Looking at Fig. 2, this seems to be the case. We note in passing that this is also consistent with the prediction given by the analysis of the Yukawa form mentioned in Sec. IV.

At sufficiently low wavenumbers, converges to a constant, given to a good approximation as , confirming that the critical exponent as predicted in Eq. (24). Importantly, at the point of exact jamming, the FCC crystal is trivially stealthy (its structure factor is identically zero up to some positive wavenumber), meaning that the critical exponent may be thought of as being infinite. This discontinuous change from the equilibrium behavior at even vanishingly small distances to jamming highlights the singular nature of jamming and underscores the need to be careful when considering limiting behavior.

It is important to note that a structure factor that scales as can be obtained by applying random, uncorrelated displacements to each particle in the crystal Gabrielli and Torquato (2004). Given the inherent anharmonicity of the system owing to the singular nature of the hard-sphere interaction potential, we are motivated to ask whether the probability distribution of the pair separation of nearest-neighbor spheres in the crystal might be, to a good approximation, statistically independent from pair to pair on certain length scales. On larger length scales (lower wavenumbers), the exact form of is in excellent agreement with a normal mode analysis based on the fictitious interparticle potential derived in Brito and Wyart (2006), suggesting underlying correlated displacements on corresponding length scales. In the second paper in this series, we investigate the consequences of this by analyzing the percolation properties of an intimately-related “cherrypit model” Atkinson et al. (2016).

### vi.2 Disordered Packings and the MRJ state

We now turn our attention to characterizing the approach to jamming in disordered packings produced by the LS and TJ algorithms as described in Sec. V. Unlike the case for the FCC crystal, MRJ-like packings do not strongly indicate signs of an incipient jammed structure at densities far below that of jamming. Furthermore, there is a strong protocol dependence yielding qualitative differences in various packing protocols’ approaches towards disordered, jammed states.

Figure 3 shows the ensemble-averaged structure factor of our packings for a variety of densities and protocols. As jamming is approached, approaches zero, implying that the jammed state is hyperuniform, in agreement with previous investigations Donev et al. (2005a); Hopkins et al. (2012). Because the data is presented on a log-log scale, the vertical offset between the nearly-jammed TJ and LS configurations corresponds to a difference in slope in their respective linear behavior in the vicinity of the origin. Interestingly, the packings produced by the TJ algorithm display anomalous behavior well before jamming is reached, including a structure factor that increases as the wavenumber approaches zero for intermediate densities. By contrast, for LS-generated packings seems to monotonically decrease for at all packing fractions leading up to jamming. This suggests that the intermediate configurations that TJ creates on its way to jamming are far from equilibrium—even at packing fractions below the freezing density Hoover and Ree (1968). For configurations that are close to jamming at a density of , we group packings according to the quantity rather than .

The TJ algorithm seems to find configurations that are consistently more disordered than those visited by LS starting at intermediate densities and continuing up to jamming ^{2}^{2}2for very low densities, the RSA configurations that are given to TJ as initial conditions are rather similar to the equilibrium fluid..
One may consider the order metric

(29) |

which may be thought of as quantifying extent to which a given configuration differs from a Poisson process (for which for all ). This order metric was used in Torquato et al. (2015) with . The order metric is also reminiscent of the direct-space order metric that we have utilized before Atkinson et al. (2013), which measures deviations in the pair correlation function from unity, as well as the two-body excess entropy Truskett et al. (2000); Green (1952); Nettleton and Green (1958); Baranyai and Evans (1989). Here, we keep finite to prevent from diverging due to the contributions of Dirac delta functions in the corresponding direct-space statistics characteristic of jammed packings; a similar issue arises in the aforementioned order metrics as discussed in Truskett et al. (2000).

Looking at Fig. 3, one can see immediately that for TJ-generated packings remains much closer to unity at intermediate densities than for LS. This difference persists up to jamming. Figure 4 shows the quantity computed for ensembles of 1000 packings created by the TJ and LS algorithms corresponding to the densities used for Fig. 3 as a function of . Consistent with other order metrics including as well as the standard bond-orientational order metric , which primarily measure short-range order, demonstrates that the packings generated by TJ are also more disordered than those produced by LS on larger length scales Atkinson et al. (2013). Thus, provides complementary information to and .

Near jamming and at sufficiently high wavenumbers, the behavior of the integrand in Eq. (29) is dominated by the contribution arising from the delta function in at . This causes to diverge toward infinity at a rate that is linear in . The thick black line in Fig. 4 illustrates the slope associated with this behavior. If the delta function is not sharp (because of a spread in nearest-neighbor distances), then will saturate to a constant. This is seen for the LS and TJ packings for and as well as for LS at . This suggests that, while near-contacts accumulate starting at low densities as TJ densifies its packings, they remain spread out over a range of pair distances beyond contact.

Because of the noise in measuring the structure factor numerically (due to both a finite number of packings in the ensemble as well as finite system sizes), the decaying oscillations converging to will eventually saturate to white noise. Therefore, will begin to grow with increasing cutoff as . The beginning of this trend is visible in Fig. 4 for our ensembles at lower densities.

Figure 5 shows the corresponding direct correlation functions for these ensembles of packings. The first salient feature is that the direct correlation function for the packings produced by the TJ algorithm exhibits a prominent peak at that is clearly visible at packing fractions as low as . This is accompanied by a steep decrease in for that is dominated by a scaling as . As mentioned in Sec. II, the analysis of a Yukawa-like beyond the core as goes to infinity Waisman (1973); Henderson and Grundke (1975); Høye and Stell (1976); Høye et al. (1976); Jadrich and Schweizer (2013) provides the rationalization for the appearance of a scaling behavior of this form. However, it only becomes dominant if the volume integral of outside the core is sufficiently large. That is, the short-ranged behavior of seems to be, fascinatingly, communicating its growing long-rangedness in the sense of the theoretical considerations of Sec. IV.

The early appearance of a delta function in at leads us to suggest that cluster formation in TJ occurs long before the packing is confined to a jamming basin. This result is likely related to the observation of Shen et al. Shen et al. (2012) that athermal packings of spheres compressed from low densities in the presence of a viscous background exhibit a “contact percolation” which is accompanied by the emergence of a nontrivial mechanical response to applied stress at densities significantly below that of jamming. Because of the manner in which the sequential linear programming algorithm searches for local optimizations in packing fraction, which require little reconfiguration at low densities, there is reason to believe that the TJ algorithm explores available configuration space in a similar fashion to the procedure of Shen et al. for low to intermediate densities.

We also note in passing that the direct correlation functions exhibit both a cusp at and a mild step discontinuity at , mirroring the features found in the pair correlation function. We have noticed that the step discontinuity observed in of jammed packings is consistently larger than what is produced as an effect of the delta function at ; this is explained in light of the result contained in Eq. (18). Note also that the particular form of Eq. (18) relies on the assumption that the packing is isostatic, meaning that there are the minimum number of backbone contact pairs necessary to ensure jamming. In the infinite system limit, this corresponds to a mean backbone coordination number of , where the vanishing term reflects the difference between collective or strict jamming Donev et al. (2005b).

Figure 6 shows plotted on a log scale. As jamming is approached, we observe that for large , confirming numerically the prediction of Eq. (20) for in the case of disordered packings This scaling behavior is difficult to obtain numerically since one must accurately obtain data for low wavenumbers in order to extract the large- behavior of . In particular, one must necessarily extrapolate , as a direct computation using Eq. (3) contains a forward scattering contribution which must be omitted. Additionally, to improve numerical stability, our data were multiplied by the Fourier transform of a narrow triangular window so that the real-space data is smoothed accordingly via convolution. Details of the procedure are give in Appendix A.

The difference in the protocols’ approaches to jamming may be readily traced back to differences in the dynamics involved: on one hand, the constant thermal motion inherent to LS acts to equilibrate packing and avoid metastable branches terminating at low-density jammed states; an aggressive expansion rate works against this, though one must worry about the algorithm becoming trapped in an unstable mechanical equilibrium (which is, by definition, not a jammed state). The possible displacements obtained by TJ are highly degenerate since, in general, there are many different displacements which allow for the same increase in the packing fraction (which is limited to a small value so that the linear approximations in the LPs’ formulation remain reasonably accurate). Therefore, TJ tends to displace spheres in a “lazy” fashion, only moving what is necessary to increase the packing fraction and no more.

Evidence of this aforementioned qualitative difference may be observed directly; Figure 7 shows snapshots of two-dimensional packings of monodisperse disks created by TJ and LS (using a rapid compression rate) at a packing fraction of . In two dimensions just as in three, we can see that the structures of the packings are qualitatively different. In particular, the TJ algorithm exhibits clustering of particles that might be quickly dispersed through thermal motion; in the absence of this, the clusters continue to combine and aggregate as jamming is approached. We point out that particles within these clusters do not necessarily contact one another; some separation is expected to remain between particles owing to the nonlinearities that are not captured in the TJ algorithm’s linear approximation to the packing problem. In LS, on the other hand, particles tend to space themselves out more uniformly through their thermal motion. We note in passing that, while this difference does not prevent the LS algorithm from discovering MRJ-like states in 3D, the two-dimensional case was recently shown to be considerably more subtle—the difference in how TJ creates jammed packings has led to the first observations of MRJ-like packings of monodisperse disks in two dimensions, whereas the LS algorithm and other standard protocols are unable to observe them, finding significantly more ordered, polycrystalline structures even under rapid compression Atkinson et al. (2014).

## Vii Conclusions and Discussion

In this work, we have compared and contrasted the approach of both ordered and disordered hard-sphere packings towards jammed states through considering the behavior of their structure factors and direct correlation functions. By considering the degree and position of singularities in as well as how they are changed by the convolutions found in Eq. (1), we have established quantitative statements about the structure of the direct correlation function with regards to features it inherits from . These relations provide a concrete means of identifying what features must be expressed in if one hopes to reproduce various details in accurately.

Moreover, we found that the LS and TJ protocols approach their respective jammed states in markedly different manners, as shown by various pair statistics. Specifically, the structure factor of TJ-generated packings shows anomalous increasing behavior for small at intermediate densities, and generally remains closer to at all densities leading up to jamming when compared to LS. The order metric compares a configuration’s pair statistics to that of an uncorrelated (Poisson) point process, which may be thought of as a maximally disordered reference state. In this sense, may be thought of as a “disorder metric”. At low to intermediate densities, suggests that packings created by TJ are more disordered on large length scales, but more ordered on short length scales as evidenced by the crossover as the truncation in the integration domain is made increasingly large. This is consistent with the intuition that TJ does not disturb the packings as they are compressed as much as LS does, leaving the large-scale characteristics similar to the initial conditions obtained from low-density RSA. On the other hand, the formation of near-contacts well before jamming may be interpreted as a sort of “ordering” that LS avoids through equilibration; therefore, LS yields configurations that are more disordered locally at densities far from jamming.

TJ shows signs of particles in close proximity at surprisingly low densities as evidenced by the appearance of a clear precursor to the delta function at and corresponding scaling within the core as . We have shown that the latter numerical observation can be predicted from theoretical considerations using a Yukawa model for . By evaluating , we see that these near-contacts that cause the delta function to appear are distributed across a range of pair separations, and the delta function’s precursor is not “sharp” until higher densities. This is to be expected because of the linear approximations that TJ makes as the packing is compressed; the inaccuracies due to nonlinear contributions are largest when large changes in the system configuration (particle translations and box deformations) are made. This is the case at densities far from jamming, where the linear approximations to the packing problem still leave a large amount of configuration space accessible. Nonetheless, this feature in the intermediate-density structures produced by the TJ algorithm suggests that it possesses important qualitative commonalities with the physical process of compressing hard-spheres embedded in a dampening background, providing a conceptual physical analog to the algorithm as witnessed in practice.

It has been suggested previously Torquato and Stillinger (2003); Hopkins et al. (2012) that the hyperuniform, linear nonanalytic behavior of for MRJ-like packings ought to give rise to a long-ranged direct correlation function which exhibits a power-law decay of . We have confirmed this numerically using simulated packings of hard-spheres generated by two very different protocols, adding to the evidence to the conjectured link between jamming and hyperuniformity Donev et al. (2005a); Torquato and Stillinger (2010) and supporting the idea that the emergence of large- scaling behavior consistent with hyperuniformity may be regarded as a structural precursor to jamming. It would be interesting to consider new semi-empirical forms for incorporating this long-range behavior in order to gain understanding regarding its structural consequences.

It is interesting that the structure factor of the FCC crystal exhibits scaling that is constant for low , but gives way to beginning at wavenumbers on the order of unity, extending to an increasingly large maximum wavenumber as jamming is approached. We noted above that this would imply that the average pair distance between any given pair of nearest-neighbor particles might be spatially uncorrelated to a good approximation. It has been observed elsewhere Ikeda and Berthier (2015) that the fluctuating component of the structure factor in disordered packings of thermally-excited soft-spheres exhibits a similar quadratic scaling at densities slightly above the jamming transition density. We have noticed that this behavior is also exhibited for disordered hard-sphere packings at packing fractions below, but close to jamming, suggesting similarly that the pair separations between nearest-neighbor particles may fluctuate in an uncorrelated manner to a good approximation.

We suspect that the aforementioned differences between the LS and TJ algorithms should be evident in other ways. In a subsequent paper Atkinson et al. (2016), we will study our hard-sphere systems in the context of two different percolation problems. In the first, we decorate the hard cores with a perfectly penetrable shell (this is known as the “cherrypit model” Torquato (2002)) and tuning the size of this shell for various configurations to explore percolation criticality. Based on our findings here, one would expect for TJ that the critical shell thickness would become very small rather quickly, whereas it might decrease more steadily towards zero for LS as jamming is approached. Moreover, one might expect to find structural differences in the percolating clusters between these two algorithms. Second, we investigate percolation by measuring the time-averaged magnitudes of pair forces in nearly-jammed, structurally-arrested configurations as a function of a minimum “threshold” force. This approach serves to “average out” fluctuations in particle position, and thus provides insight into the role of fluctuations in the incipient structure of disordered, jammed systems. Both approaches provide insight into the jamming process using static structural features, expanding upon the work presented here.

###### Acknowledgements.

The authors thank Yang Jiao and Adam Hopkins for their careful reading of this manuscript and valuable comments. This work was supported in part by the National Science Foundation under Grant No. DMS-1211087.## Appendix A Numerical procedure for obtaining from data

In this appendix, we provide procedural details regarding our numerical computation of the direct correlation function used to ascertain the large- behavior.

We begin by measuring at wavevectors that are integer combinations of the columns of the reciprocal matrix defined as , where the columns of span the fundamental cell of our system in direct space. Measurements are then binned and averaged according to their wavenumber with a bin width of , so that we have reported measurements at wavenumbers given by for . Because the density of wavevectors scales as , we randomly select wavevectors with a probability proportional to the inverse of this density such that, for higher wavenumbers, the expected number of measurements per bin is . At smaller wavenumbers, the structure factor is measured at every available wavevector.

Once we have obtained data for all of our packings, we perform an ensemble average. If there are empty bins, then we linearly interpolate a value for those bins, expecting that this is representative of the large-system limit. We also linearly extrapolate down to zero if any bins are missing data; our results are qualitatively robust against small variations in this extrapolation.

Our data is then converted via Eq. (5) to give us . Because of the asymptotic behavior of , we find that it is necessary to apply a convolution in order to eliminate artifacts caused by difficulties associated with numerically integrating this high-frequency behavior. We do this by multiplying with the Fourier transform of a triangular window in direct space given by . The Fourier transform of this radial function is

(30) |

For a general three-dimensional, radial function , the Fourier transform and its inverse may be expressed Torquato (2002) as

(31) | |||||

(32) |

Note that both the forward and inverse transforms are equivalent up to their scaling coefficients. In order to take advantage of the usual one-dimensional fast Fourier transform algorithm for our three-dimensional, radial , we take our discrete for and compute

(33) |

We then compute the usual one-dimensional inverse FFT on this data, defined here as

(34) | |||||

where . We then apply the prefactor and a phase correction because to correct for the fact that the index corresponds to to obtain

(35) |

where and . Through analogy with Eq. (32), one expects that is completely imaginary, while one expects the Fourier transform to be completely real-valued. This is reconciled by dropping the imaginary unit from ; doing so is justified since the imaginary prefactor is expected if one applies Euler’s formula to the exponential term in Eq. (34), but is missing as a prefactor to the sine term in Eq. (32). Finally, the value for corresponding to is obtained by integrating numerically according to Eq. (32).

Once the direct space has been found (represented discretely through ), one must then be sure to truncate the data at where is the width of the simulation box; data beyond this point is subject to finite-size artifacts.

## References

- Bernal (1965) J. D. Bernal, “Liquids: Structure, properties, solid interactions,” (Elsevier, Amsterdam, 1965) pp. 25–50.
- Zallen (1983) R. Zallen, The Physics of Amorphous Solids (Wiley, New York, 1983).
- Cates et al. (1999) M. E. Cates, J. P. Wittmer, J. P. Bouchaud, and P. Claudin, Chaos 9, 511 (1999).
- Torquato and Stillinger (2001) S. Torquato and F. H. Stillinger, J. Phys. Chem. B 105, 11849 (2001).
- Torquato et al. (2003) S. Torquato, A. Donev, and F. H. Stillinger, Int. J. Solids Struct. 40, 7143 (2003).
- Torquato and Stillinger (2003) S. Torquato and F. H. Stillinger, Phys. Rev. E 68, 041113 (2003).
- Zachary and Torquato (2009) C. E. Zachary and S. Torquato, J. Stat. Mech.: Theory & Exp. 2009, P12015 (2009).
- O’Hern et al. (2003) C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev. E 68, 011306 (2003).
- Manoharan et al. (2003) V. N. Manoharan, M. T. Elsesser, and D. J. Pine, Science 301, 483 (2003).
- Makse et al. (2004) H. A. Makse, J. Brujic, and S. F. Edwards, The Physics of Granular Media , 45 (2004).
- Xu et al. (2005) N. Xu, J. Blawzdziewicz, and C. S. O’Hern, Phys. Rev. E 71, 061306 (2005).
- Zohdi (2006) T. I. Zohdi, Mech. Mater. 38, 969 (2006).
- Majmudar et al. (2007) T. Majmudar, M. Sperl, S. Luding, and R. P. Behringer, Phys. Rev. Lett. 98, 058001 (2007).
- Gevertz and Torquato (2008) J. L. Gevertz and S. Torquato, PLoS Comput. Biol. 4, e1000152 (2008).
- Karayiannis et al. (2009) N. C. Karayiannis, K. Foteinopoulou, and M. Laso, Phys. Rev. E 80, 011307 (2009).
- Parisi and Zamponi (2010) G. Parisi and F. Zamponi, Rev. Mod. Phys. 82, 789 (2010).
- Jiang et al. (2010) S. Jiang, Q. Chen, M. Tripathy, E. Luijten, K. S. Schweizer, and S. Granick, Adv. Mater. 22, 1060 (2010).
- Xu et al. (2011) N. Xu, D. Frenkel, and A. J. Liu, Phys. Rev. Lett. 106, 245502 (2011).
- Charbonneau et al. (2012) P. Charbonneau, E. I. Corwin, G. Parisi, and F. Zamponi, Phys. Rev. Lett. 109, 205501 (2012).
- Dagois-Bohy et al. (2012) S. Dagois-Bohy, B. P. Tighe, J. Simon, S. Henkes, and M. van Hecke, Phys. Rev. Lett. 109, 095703 (2012).
- Gillman et al. (2013) A. Gillman, K. Matouš, and S. Atkinson, Phys. Rev. E 87, 022208 (2013).
- Klatt and Torquato (2014) M. A. Klatt and S. Torquato, Phys. Rev. E 90, 052120 (2014).
- Yamchi et al. (2015) M. Z. Yamchi, S. S. Ashwin, and R. K. Bowles, Phys. Rev. E 91, 022301 (2015).
- Royer and Chaikin (2015) J. R. Royer and P. M. Chaikin, Proc. Natl. Acad. Sci. 112, 49 (2015).
- Schaller et al. (2015) F. M. Schaller, S. C. Kapfer, J. E. Hilton, P. W. Cleary, K. Mecke, C. D. Michele, T. Schilling, M. Saadatfar, M. SchrÃ¶ter, G. W. Delaney, and G. E. SchrÃ¶der-Turk, Europhys. Lett. 111, 24002 (2015).
- Ramola and Chakraborty (2016) K. Ramola and B. Chakraborty, arXiv preprint arXiv:1604.06148 (2016).
- Silbert and Silbert (2009) L. E. Silbert and M. Silbert, Phys. Rev. E 80, 041304 (2009).
- Zachary et al. (2011a) C. E. Zachary, Y. Jiao, and S. Torquato, Phys. Rev. E 83, 051308 (2011a).
- Zachary et al. (2011b) C. E. Zachary, Y. Jiao, and S. Torquato, Phys. Rev. E 83, 051309 (2011b).
- Hopkins et al. (2012) A. B. Hopkins, F. H. Stillinger, and S. Torquato, Phys. Rev. E 86, 021505 (2012).
- Chieco et al. (2016) A. Chieco, C. Goodrich, A. Liu, and D. Durian, in APS Meeting Abstracts (2016).
- (32) S. Atkinson, G. Zhang, A. B. Hopkins, and S. Torquato, Phys. Rev. E 94, 012902 (2016).
- Marcotte et al. (2013) É. Marcotte, F. H. Stillinger, and S. Torquato, J. Chem. Phys. 138, 12 (2013).
- Hansen and McDonald (1990) J.-P. Hansen and I. R. McDonald, Theory of simple liquids (Elsevier, 1990).
- Torquato et al. (2000) S. Torquato, T. M. Truskett, and P. G. Debenedetti, Phys. Rev. Lett. 84, 2064 (2000).
- Torquato and Stillinger (2010) S. Torquato and F. H. Stillinger, Rev. Mod. Phys. 82, 2633 (2010).
- Donev et al. (2005a) A. Donev, F. H. Stillinger, and S. Torquato, Phys. Rev. Lett. 95, 090604 (2005a).
- Skoge et al. (2006) M. Skoge, A. Donev, F. H. Stillinger, and S. Torquato, Phys. Rev. E 74, 041127 (2006).
- Jiao et al. (2011) Y. Jiao, F. H. Stillinger, and S. Torquato, J. Appl. Phys. 109, 013508 (2011).
- Torquato and Jiao (2010) S. Torquato and Y. Jiao, Phys. Rev. E 82, 061302 (2010).
- Atkinson et al. (2013) S. Atkinson, F. H. Stillinger, and S. Torquato, Phys. Rev. E 88, 062208 (2013).
- Mau and Huse (1999) S.-C. Mau and D. A. Huse, Phys. Rev. E 59, 4396 (1999).
- Pronk and Frenkel (2003) S. Pronk and D. Frenkel, Phys. Rev. Lett. 90, 255501 (2003).
- (44) To be mathematically precise, the statistical functions , , , and are all functions which map ; however, since these functions are often isotropic, we often express them in terms of or as appropriate.
- Torquato (2002) S. Torquato, Random Heterogeneous Materials: Microstructure and Macroscopic Properties (Springer-Verlag, New York, 2002).
- Waisman (1973) E. Waisman, Mol. Phys. 25, 45 (1973).
- Høye and Stell (1976) J. Høye and G. Stell, Mol. Phys. 32, 195 (1976).
- Høye et al. (1976) J. Høye, G. Stell, and E. Waisman, Mol. Phys. 32, 209 (1976).
- Jadrich and Schweizer (2013) R. Jadrich and K. S. Schweizer, J. Chem. Phys. 139, 054502 (2013).
- (50) In one dimension, we have a closed-form solution Percus (1964) for the direct correlation function of an equilibrium system of equal-sized hard rods. This solution posesses linear behavior in the vicinity of , which is the Coulombic potential for this dimension. Here, we have shown that three dimensions also exhibits Coulombic behavior for small . This leads us to suggest that mimics the Coulombic interaction in any dimension (e.g. in two dimensions and for ). .
- Donev et al. (2005b) A. Donev, S. Torquato, and F. H. Stillinger, Phys. Rev. E 71, 011105 (2005b).
- Salsburg and Wood (1962) Z. Salsburg and W. Wood, J. Chem. Phys. 37, 798 (1962).
- Lubachevsky and Stillinger (1990) B. D. Lubachevsky and F. H. Stillinger, J. Stat. Phys. 60, 561 (1990).
- Gabrielli and Torquato (2004) A. Gabrielli and S. Torquato, Phys. Rev. E 70, 041105 (2004).
- Brito and Wyart (2006) C. Brito and M. Wyart, Europhys. Lett. 76, 149 (2006).
- Atkinson et al. (2016) S. Atkinson, F. H. Stillinger, and S. Torquato, (2016), in progress.
- Hoover and Ree (1968) W. G. Hoover and F. H. Ree, J. Chem. Phys. 49, 3609 (1968).
- (58) For very low densities, the RSA configurations that are given to TJ as initial conditions are rather similar to the equilibrium fluid.
- Torquato et al. (2015) S. Torquato, G. Zhang, and F. Stillinger, Phys. Rev. X 5, 021020 (2015).
- Truskett et al. (2000) T. M. Truskett, S. Torquato, and P. G. Debenedetti, Phys. Rev. E 62, 993 (2000).
- Green (1952) H. S. Green, “The molecular theory of fluids,” (North-Holland Pub. Co., Amsterdam, 1952) pp. 70–73.
- Nettleton and Green (1958) R. Nettleton and M. Green, J. Chem. Phys. 29, 1365 (1958).
- Baranyai and Evans (1989) A. Baranyai and D. J. Evans, Phys. Rev. A 40, 3817 (1989).
- Henderson and Grundke (1975) D. Henderson and E. Grundke, J. Chem. Phys. 63, 601 (1975).
- Shen et al. (2012) T. Shen, C. S. O’Hern, and M. Shattuck, Phys. Rev. E 85, 011308 (2012).
- Atkinson et al. (2014) S. Atkinson, F. H. Stillinger, and S. Torquato, Proc. Natl. Acad. Sci. 111, 18436 (2014).
- Ikeda and Berthier (2015) A. Ikeda and L. Berthier, Phys. Rev. E 92, 012309 (2015).
- Percus (1964) J. Percus, in The equilibrium theory of classical fluids: a lecture note and reprint volume, edited by H. L. Frisch and J. L. Lebowitz (WA Benjamin, 1964).