## Abstract

Packing problems, such as how densely objects can fill a volume, are among the most ancient and persistent problems in mathematics and science. For equal spheres, it has only recently been proved that the face-centered cubic lattice has the highest possible packing fraction . It is also well known that certain random (amorphous) jammed packings have φ ≈ 0.64. Here, we show experimentally and with a new simulation algorithm that ellipsoids can randomly pack more densely—up to φ= 0.68 to 0.71for spheroids with an aspect ratio close to that of M&M's Candies—and even approach φ ≈ 0.74 for ellipsoids with other aspect ratios. We suggest that the higher density is directly related to the higher number of degrees of freedom per particle and thus the larger number of particle contacts required to mechanically stabilize the packing. We measured the number of contacts per particle *Z* ≈ 10 for our spheroids, as compared to *Z* ≈ 6 for spheres. Our results have implications for a broad range of scientific disciplines, including the properties of granular media and ceramics, glass formation, and discrete geometry.

The structure of liquids, crystals, and glasses is intimately related to volume fractions of ordered and disordered (random) hard-sphere packings, as are the transitions between these phases (*1*). Packing problems (*2*) are of current interest in dimensions higher than three for insulating stored data from noise (*3*), and in two and three dimensions in relation to flow and jamming of granular materials (*4*–*6*) and glasses (*7*). Of particular interest is random packing, which relates to the ancient (economically important) problem of how much grain a barrel can hold.

Many experimental and computational algorithms produce a relatively robust packing fraction (relative density) φ ≈ 0.64 for randomly packed monodisperse spheres as they proceed to their limiting density (*8*). This number, widely designated as the random close packing (RCP) density, is not universal but generally depends on the packing protocol (*9*). RCP is an ill-defined concept because higher packing fractions are obtained as the system becomes ordered, and a definition for randomness has been lacking. A more recent concept is that of the maximally random jammed (MRJ) state, corresponding to the least ordered among all jammed packings (*9*). For a variety of order metrics, it appears that the MRJ state has a density of φ ≈ 0.637 and is consistent with what has traditionally been thought of as RCP (*10*). Henceforth, we refer to this random form of packing as the MRJ state.

We report on the density of the MRJ state of ellipsoid packings as asphericity is introduced. For both oblate and prolate spheroids, φ and *Z* (the average number of touching neighbors per particle) increase rapidly, in a cusp-like manner, as the particles deviate from perfect spheres. Both reach high densities such as φ ≈ 0.71, and general ellipsoids pack randomly to a remarkable φ ≈ 0.735, approaching the density of the crystal with the highest possible density for spheres (*11*) . The rapid increases are unrelated to any observable increase in order in these systems that develop neither crystalline (periodic) nor liquid crystalline (nematic or orientational) order.

Our experiments used two varieties of M&M's Milk Chocolate Candies: regular and baking (“mini”) candies (*12*). Both are oblate spheroids with small deviations from true ellipsoids, Δ*r*/*r* < 0.01. Additionally, M&M's Candies have a very low degree of polydispersity (principal axes 2*a* = 1.34 ± 0.02 cm, 2*b* = 0.693 ± 0.018 cm, *a*/*b* = 1.93 ± 0.05 for regular; 2*a* = 0.925 ± 0.011 cm, 2*b* = 0.493 ± 0.018 cm, *a*/*b* = 1.88 ± 0.06 for minis). Several sets of experiments were performed to determine the packing fraction. A square box, 8.8 cm by 8.8 cm, was filled to a height of 2.5 cm while shaking and tapping the container. The actual measurements were performed by adding 9.0 cm to the height and excluding the contribution from the possibly layered bottom. After measuring the average mass, density, and volume of the individual candies, the number of candies in the container and their volume fraction could be simply determined by weighing. These experiments yielded φ = 0.665 ± 0.01 for regulars and φ = 0.695 ± 0.01 for minis. The same technique was used for 3.175 = mm ball bearings (spheres) and yielded φ = 0.625 ± 0.01. A second set of experiments was performed by filling 0.5-, 1-, and 5-liter round flasks (to minimize ordering due to wall effects) with candies by pouring them into the flasks while tapping (5 liters corresponds to about 23,000 minis or 7500 regulars) (Fig. 1A). The volume fractions found in these more reliable studies were φ = 0.685 ± 0.01 for both the minis and regulars (*13*). The same procedure for 30,000 ball bearings in the 0.5-liter flask yielded φ = 0.635 ± 0.01, which is close to the accepted MRJ density.

A 5-liter sample of regular candies similar to that shown in Fig. 1A was scanned in a medical magnetic resonance imaging device at Princeton Hospital. For several planar slices, the direction Θ (with respect to an arbitrary axis) of the major elliptical axis was manually measured and the two-dimensional nematic order parameter *S*_{2} = 〈2 cos^{2} Θ – 1〉 was computed, yielding *S*_{2} ≈ 0.05. This is consistent with the absence of orientational order in the packing (*14*).

Our simulation technique generalizes the Lubachevsky-Stillinger (LS) sphere-packing algorithm (*15*, *16*) to the case of ellipsoids. The method is a hard-particle molecular dynamics (MD) algorithm for producing dense disordered packings. Initially, small ellipsoids are randomly distributed and randomly oriented in a box with periodic boundary conditions and without any overlap. The ellipsoids are given velocities and their motion followed as they collide elastically and also expand uniformly. After some time a jammed state with a diverging collision rate is reached and the density reaches a maximal value. A novel event-driven MD algorithm (*17*) was used to implement this process efficiently, based on the algorithm used in (*15*) for spheres and similar to the algorithm used for needles in (*18*). A typical configuration of 1000 oblate ellipsoids (aspect ratio α = *b*/*a* = 1.9^{–1} ≈ 0.526) is shown in Fig. 1B, with density of φ ≈ 0.70 and nematic order parameter *S* ≈ 0.02 to 0.05.

We have verified that the sphere packings produced by the LS algorithm are jammed according to the rigorous hierarchical definitions of local, collective, and strict jamming (*19*, *20*). Roughly speaking, these definitions are based on mechanical stability conditions that require that there be no feasible local or collective particle displacements and/or boundary deformations. On the basis of our experience with spheres (*10*), we believe that our algorithm (with rapid particle expansion) produces final states that represent the MRJ state well. The algorithm closely reproduces the packing fraction measured experimentally.

The density of simulated packings of 1000 particles is shown in Fig. 2A. Note the two clear maxima with φ ≈ 0.71, already close to the 0.74 for the ordered face-centered cubic (fcc)/hexagonal close-packed (hcp) packing, and the cusp-like minimum near α = 1 (spheres). Previous simulations for random sequential addition (RSA) (*21*), as well as gravitational deposition (*22*), produce a similarly shaped curve, with a maximum at nearly the same aspect ratios α ≈ 1.5 (prolate) or α ≈ 0.67 (oblate), but with substantially lower volume fractions (such as φ ≈ 0.48 for RSA).

Why does the packing fraction initially increase as we deviate from spheres? The rapid increase in packing fraction is attributable to the expected increase in the number of contacts resulting from the additional rotational degrees of freedom of the ellipsoids. More contacts per particle are needed to eliminate all local and collective degrees of freedom and ensure jamming, and forming more contacts requires a denser packing of the particles. In the inset in Fig. 2B, the central circle is locally jammed. A uniform vertical compression preserves φ, but the central ellipsoid can rotate and free itself and the packing can densify. The decrease in the density for very aspherical particles could be explained by strong exclusion-volume effects in orientationally disordered packings (*23*). Results resembling those shown in Fig. 2A are also obtained for isotropic random packings of spherocylinders (*23*, *24*), but an argument based on “caging” (not jamming) of the particles was given to explain the increase in density as asphericity is introduced. Spherocylinders have a very different behavior for ordered packings from ellipsoids (the conjectured maximal density is , which is significantly higher than for ellipsoids), and also cannot be oblate and are always axisymmetric. The similar positioning of the maximal density peak for different packing algorithms and particle shapes indicates the relevance of a simple geometrical explanation.

By introducing orientational and translational order, it is expected that the density of the packings can be further increased, at least up to 0.74. As shown in Fig. 3 for two dimensions, an affine deformation (stretch) of the densest disk packing produces an ellipse packing with the same volume fraction. However, this packing, although the densest possible, is not strictly jammed (i.e., it is not rigid under shear transformations). The figure shows through a sequence of frames how one can distort this collectively jammed packing (*20*), traversing a whole family of densest configurations. This mechanical instability of the ellipse packing as well as the three-dimensional ellipsoid packing arises from the additional rotational degrees of freedom and does not exist for the disk or sphere packing.

There have been conjectures (*25*, *26*) that frictionless random packings have just enough constraints to completely statically define the system (*27*), *Z* = 2*f* (i.e., that the system is isostatic), where *f* is the number of degrees of freedom per particle (*f* = 3 for spheres, *f* = 5 for spheroids, and *f* = 6 for general ellipsoids) (*28*). If friction is strong, then fewer contacts are needed, *Z* = *f* + 1 (*29*). Experimentally, *Z* for spheres was determined by Bernal and Mason by coating a system of ball bearings with paint, draining the paint, letting it dry, and counting the number of paint spots per particle when the system was disassembled (*30*). Their results gave *Z* ≈ 6.4, surprisingly close to isostaticity for frictionless spheres (*31*).

We performed the same experiments with the M&M's, counting the number of true contacts between the particles (*32*). A histogram of the number of touching neighbors per particle for the regular candies is shown in Fig. 4. The average number is *Z* = 9.82. In simulations a contact is typically defined by a cutoff on the gap between the particles. Fortunately, over a wide range (10^{–9} to 10^{–4}) of contact tolerances, *Z* is reasonably constant. Superposed in Fig. 4 is the histogram of contact numbers obtained for simulated packings of oblate ellipsoids for α = 0.526, from which we found *Z* ≈ 9.80. In Fig. 2B we show *Z* as a function of aspect ratio α (*33*). As with the volume fraction, the contact number appears singular at the sphere value and rises sharply for small deviations. Unlike φ, however, *Z* does not decrease for large aspect ratios, but rather appears to remain constant.

We expect that fully aspherical ellipsoids, which have *f* = 6, will require even more contacts for jamming (*Z* = 12 according to the isostatic conjecture) and larger φ. Results from simulations of ellipsoids with axes *a* = α^{–1}, *b* = 1, and *c* = α (where α measures the asphericity) are included in Fig. 2A. At α ≈ 1.3 we obtain a surprisingly high density of φ ≈ 0.735, with no significant orientational ordering. The maximum contact number observed in Fig. 2B is *Z* ≈ 11.4. It is interesting that for both spheroids and general ellipsoids, *Z* reaches a constant value at approximately the aspect ratio for which the density has a maximum. This supports the claim that the decrease in density for large α is due to exclusion volume effects.

The putative nonanalytic behavior of *Z* and φ at α = 1 is striking and is evidently related to the randomness of the jammed state. Crystal close packings of spheres and ellipsoids show no such singular behavior, and in fact φ and *Z* are independent of α for small deviations from unity. On the other hand, for random packings, the behavior is not discontinuous, whereas the number of degrees of freedom jumps from three to five (or six) as soon as α deviates from 1. In several industrial processes such as sintering and ceramic formation, interest exists in increasing the density and number of contacts of powder particles to be fused. If ellipsoidal instead of spherical particles are used, we may increase the density of a randomly poured and compacted powder to a value approaching that of the densest (fcc) lattice packing.