1. Introduction

Over a broad range of length scales, many-particle systems exhibit a rich variety of structures with varying degrees of long-range order, spanning from crystals, quasicrystals, and polycrystals to amorphous solids and liquids. Consequently, it is natural to focus attention on the statistical mechanics of the arrangement of the particles. In the case of a macroscopic system containing a large number N of particles, a full configurational description of that system usually is neither feasible, desirable, nor necessary. For most practical purposes, it suffices to determine, or to describe, the distribution functions of low orders n N. Conventionally, this information is conveyed in the form of correlation functions. For statistically homogeneous systems consisting of identical spherical particles in a volume V, these correlation functions are defined so that ⁿg⁽ⁿ⁾(r₁,r₂,...,r_n) is proportional to the probability density for simultaneously finding n particles at locations r₁, r₂, ..., r_n within the system,¹ where = N/V is the number density. With this convention, each g⁽ⁿ⁾ approaches unity when all particle positions become widely separated within V.

The present study concerns the special circumstances for which the constituent particles are spherically symmetric and identical, and the system is statistically homogeneous and isotropic. These conditions can be satisfied if the system contains a single fluid or amorphous solid phase. The correlation function of primary interest is g⁽²⁾(r₁₂), depending configurationally just on scalar pair distance r₁₂, and thus specifying how many pair distances of a given length occur statistically within the system. The third-order function g⁽³⁾(r₁₂,r₁₃,r₂₃) reveals how these pair separations are linked into triangles. This additional information strictly speaking cannot be inferred from the knowledge of g g⁽²⁾ alone, although the Kirkwood superposition approximation^1,2 presumes to fill that knowledge gap. The fourth-order function g⁽⁴⁾ controls the assembly of triangles into tetrahedra, and is the lowest-order correlation function that is sensitive to chirality of the medium.

On account of their probability interpretation, all of the g⁽ⁿ⁾ must be non-negative; in particular, for all r 0 we must have

Two recent studies^4,5 have examined the theoretical possibility of controlling pair interactions in an isothermal many-particle system over a nonvanishing density range, starting at = 0, in such a way that g(r) remains invariant over that range. The cases examined have assigned forms to the invariant g that were the zero-density limits appropriate for rigid rods, disks, and spheres,⁴ and for the hard core plus square well pair potential.⁵ In each of these examples an upper terminal density * could be identified such that over

The principal objective of the present project is to apply the g-invariance technique to the still-challenging problem of random packings of spheres. It has now been well established that the old concepts of "random loose packing" and "random close packing" are ill-defined.⁷ Instead, a nontrivial density range exists over which irregular packings of various types (locally, collectively, or strictly jammed¹¹) can be formed, the preferred densities and local structures of which depend on the preparation algorithm. Our objective has been to study the logical connections between qualifying rigid-sphere g's and the maximum corresponding densities, emerging as the terminal *'s.

Section 2 introduces a parametric family of qualitatively reasonable, but functionally elementary, pair correlation functions for the amorphous-state sphere spatial arrangement problem. This family contains a set of five adjustable parameters, whose values must of course be consistent with the two basic inequalities (1) and (3). Section 3 describes a numerical search procedure over these parameters, and its results, the goal of which was to produce the largest *. We believe it is significant that even with such a simple parametric family of pair correlation functions, * can come close to that obtained in many experimental and simulational preparations of random sphere packings. Section 4 offers some interpretive remarks stimulated by the numerical results in section 3 and indicates the natural and useful directions for future investigation. Finally, in an appendix, we derive a lower bound on the maximum density * for random sphere packings in d dimensions and show that it is sharper than a well-known lower bound for regular lattice packings for d 3.

2. Model Family

Our interest in this paper is to study models in which long-range order is suppressed and short-range order is controlled. Information that is available from previous determinations of local order in amorphous sphere packings provides useful guidance in choosing a model family of functions for the present investigation. In particular, we note that a survey of several experimental and computer-simulation protocols,^8-10 using distinct packing preparation procedures, appear to agree on the presence of some qualitative features. (Figure 1 shows the pair correlation function for a dense random packing of spheres as generated by us employing the Lubachevsky-Stillinger "compression" protocol.¹⁰) Using the sphere-pair distance of closest approach as the natural length unit, these g(r) attributes are the following:

Figure 1 Pair correlation function g(r) vs r as obtained by averaging over 10 configurations at a packing fraction = 0.64. The "binned" peak value of g(1) (not shown) is approximately 24.

(i) Obviously, g(r) must vanish for all 0 r 1.

(ii) On account of the jamming, virtually all spheres (a few "rattler" spheres can be present as exceptions) are rigidly in contact with neighbors. The number of such contacts must average at least 4 to meet the definition of "local" jamming.¹¹

(iii) For r 1, g(r) displays finite-amplitude oscillations about unity, which decay to zero with increasing r. The length scale of these oscillations is roughly comparable to the sphere diameter.

(iv) A pair of distinctive g(r) peaks appear at distances approximately equal to and . These are often termed a "split second peak", and appear in modified form for amorphous deposits of soft-sphere and attracting particles.¹⁴

(v) As r increases through r = 2, g(r) experiences a discontinuous drop in magnitude.

It is not the objective of the present work to try to include all of these features slavishly. Instead, we have chosen as a first step to represent only attributes i-iii by a simple parametric function family, and to see how close the largest corresponding density would come to the approximate experimental range^12,13 of "random" packing densities:¹⁵

Consequently, we have elected to write g(r) as a linear combination of three portions, corresponding respectively to (i), (ii), and (iii) above:

The Fourier transforms required to evaluate the structure factor,

The five adjustable parameters Z, A, B, C, and D are subject to some obvious constraints. Clearly, we must demand that

3. Numerical Search Procedure and Results

The problem of finding the maximal packing fraction *(Z,A,B,C,D) can be posed as an optimization problem. Specifically, this can be posed as a two-level "min-max" problem: one wants to maximize (Z,A,B,C,D) over the parameters Z, A, B, C, and D with the restrictions (16)-(18) such that the minimum of S(k;Z,A,B,C,D) in the variable k and the minimum of g(r;Z,A,B,C,D) in the variable r are both non-negative, i.e.,

The search procedure is implemented for six different cases summarized in Table 1. In the first case (case I), no restrictions on the functions, other than the ones described above, are imposed. In the remaining cases, we impose additional restrictions. In particular, it is known that dense random packings are typically spatially uniform, implying that S(k=0) 0. Therefore, in some instances, we carry out the search subject to the condition that S(k=0) = 0.

Case I. Not surprisingly, the least restrictive case yields the largest value of the maximal packing fraction: * = 0.627. The corresponding values of the parameters are listed in Table 2. Figure 2 shows the structure factor and pair correlation function. These functions reveal structural features that are not characteristic of typical dense random packings obtained either experimentally or computationally. For example, the structure factor at k = 0 is unusually high, implying significant density fluctuations in the infinite-wavelength limit. Moreover, at three different finite wavelengths (k 3.1, k 6.3, and k 10) the structure factor is essentially zero or nearly zero, implying vanishingly small density fluctuations at these wavelengths. Atypically, the pair correlation function g(r) exhibits appreciable oscillations before attaining its long-range value at about r = 8. This feature could be the result of polycrystallinity, but the fact that g(r) vanishes at r 1.4 (another anomalous characteristic) evidently eliminates the possibility that such putative crystallites are face-centered cubic arrangements. Interestingly, the coordination number Z 5.8 is approximately equal to the value observed in experiments and simulations of typical dense random packings.^7-9 Theoretical arguments have been put forth predicting that Z = 6 for random packings of identical frictionless spheres in three dimensions.^18-20

Figure 2 Structure factor (a) and pair correlation function (b) for case I. Note the appearance of a vertical line at contact in (b), indicating a -function contribution there.

Case II. In the second case, we conduct the search under the condition that S(k=0) = 0. For this condition to hold, Z must be given by

Figure 3 Structure factor (a) and pair correlation function (b) for case II. Note the appearance of a vertical line at contact in (b), indicating a -function contribution there.

Case III. In the third case, we suppress the damped-oscillating component [g_III(r) = G_III(k) = 0]. Here we find that * = 0.41 and Z = 3.1504.

Case IV. In the fourth case, we suppress the damped-oscillating component [g_III(r) = G_III(k) = 0] and we also impose the condition S(k=0) = 0. This problem can be solved exactly. Here we find

Case V. In the fifth case, we suppress the -function component (Z = 0). Here we obtain * = 0.375 (see Tables 1 and 2).

Case VI. In the sixth case, we suppress the -function component (Z = 0) and we also impose the condition S(k=0) = 0. We find that * = 0.3535 (see Tables 1 and 2).

4. Concluding Remarks

For the family of pair correlation functions specified by (8), the optimal packing fraction is characterized by unusual structural features such as substantially large density fluctuations in the infinite-wavelength limit, vanishing or nearly vanishing fluctuations at several finite wavelength values, and an interparticle radial distance (r 1.4) at which particle centers are prohibited. Nonetheless, the maximal packing fraction and coordination number (* = 0.627 and Z = 5.8) are consistent with values for dense random packings generated experimentally and computationally. Clearly, however, properties iv and v ("split second peak" and the discontinuous drop at r = 2) that are characteristic of typical random packings are absent in the optimal solution. In future work, one may want to consider other families of functions that are not as smooth as (8) away from r = 1, e.g., piecewise continuous functions. Such extensions would quantify the significance of properties iv and v in raising * above the optimal value of 0.627, while presumably driving S(k=0) downward toward zero.

Assuming that the optimal packing can actually be realized, there remain many open questions. Do the spheres form a contacting percolating network? Given the high density that is achieved, we suspect that the answer is in the affirmative. If so, what is the geometry of the contact set? Are rattlers present in the optimal solution? Is the packing locally, collectively, or strictly jammed?¹¹ The answers to all of these questions would greatly be facilitated if we could determine whether there are packings that achieve the optimal solution. This can be accomplished using stochastic reconstruction techniques that enable one to obtain realizations of sphere packings that have a targeted pair correlation function or structure factor.²¹ We will attempt such a reconstruction in a future study. Another interesting extension of the present work is to generalize the family of pair correlation functions to the case of spheres of different sizes and to determine the maximal packing fraction.

Acknowledgment

We are grateful to Aleksandar Donev for his help with the interval-arithmetic calculations. S.T. was supported by the Petroleum Research Fund as administered by the American Chemical Society.

Appendix: d-Dimensional Generalization of Case IV

In this appendix, we obtain an exact expression for the optimal value of * for the d-dimensional generalization of case IV. A consequence of this result is a lower bound on * for random sphere packings in d dimensions, which we compare to a well-known lower bound for regular lattice packings.

Consider the evaluation of the structure factor S(k) of relation 12 in d dimensions but without the short-ranged (damped-oscillating) contribution, i.e.,

It immediately follows that

In very high dimensions (d ~ 1000), the densest known packings are nonregular lattices.²³ Thus, it is of interest to compare the lower bound (36) to the Minkowski-Hlawka theorem,²³ which gives a lower bound on the maximum packing fraction for d-dimensional regular lattices of identical spheres: