A Geometric-Structure Theory for Maximally Random Jammed Packings

Maximally random jammed (MRJ) particle packings can be viewed as prototypical glasses in that they are maximally disordered while simultaneously being mechanically rigid. The prediction of the MRJ packing density ϕMRJ, among other packing properties of frictionless particles, still poses many theoretical challenges, even for congruent spheres or disks. Using the geometric-structure approach, we derive for the first time a highly accurate formula for MRJ densities for a very wide class of two-dimensional frictionless packings, namely, binary convex superdisks, with shapes that continuously interpolate between circles and squares. By incorporating specific attributes of MRJ states and a novel organizing principle, our formula yields predictions of ϕMRJ that are in excellent agreement with corresponding computer-simulation estimates in almost the entire α-x plane with semi-axis ratio α and small-particle relative number concentration x. Importantly, in the monodisperse circle limit, the predicted ϕMRJ = 0.834 agrees very well with the very recently numerically discovered MRJ density of 0.827, which distinguishes it from high-density “random-close packing” polycrystalline states and hence provides a stringent test on the theory. Similarly, for non-circular monodisperse superdisks, we predict MRJ states with densities that are appreciably smaller than is conventionally thought to be achievable by standard packing protocols.

Predicting the packing density φ of jammed disordered packings as a function of particle characteristics, such as particle shape and size distributions, is a long-standing and notoriously difficult problem. This is because the final states depend on many factors, such as the preparation technique, degree of order and interparticle friction, to mention just a few 12,13,35 . The determination of the maximal packing density φ max for non-space-filling particles has been a source of fascination for centuries 36,37 . Recently, general organizing principles and conjectures have been proposed that enable one to predict dense packing arrangements of frictionless nonspherical particles based on the geometrical properties of the shapes [38][39][40][41] . This allowed the prediction of φ max of multiply connected nonconvex bodies such as tori 42 . The estimation of the MRJ density φ MRJ is even more challenging. Whereas predicting φ max only requires the consideration of relatively simple ordered arrangements of a small number of particles in a periodic fundamental cell, to predict φ MRJ it is necessary to capture complex structural features of very large disordered packings that simultaneously possess strong short-range and hyperuniform long-range correlations 13,18,[20][21][22] .
A number of schemes to predict RCP densities have been developed. Recently, a 3D "granocentric" model for random packings of colloidal spheres has been devised, which distills the complexity of global packing into local stochastic processes and yields accurate estimates of the overall packing density for a variety of particle size distributions 43 . Moreover, a mean-field theory for RCP based on Edwards' ensemble has been developed, which leads to a predictive framework to calculate packing density for a wide class of particles including both spherical and non-spherical shapes (e.g., ellipsoids, spherocylinders and polyhedra) based on the distribution of Voronoi volumes associated with a particle 44,45 .
An important family of two-dimensional nonspherical shapes are given by superdisks. A superdisk is a centrally symmetric body defined by where x 1 and x 2 are Cartesian coordinates and p ∈ [0, ∞] is the deformation parameter, which indicates to what extent the particle shape has deformed from that of a circular disk (p = 1). As p continuously increases from unity to infinity, a superdisk continuously changes its shape from a circle to a square. As p decreases from 1 to 1/2, another family of square-symmetric superdisks is obtained, but with the symmetry axes rotated 45 degrees with respect to that of the first family. When p < 1/2, the superdisk is concave and becomes a cross in the limit p→ 0. In this paper, using the geometric-structure approach 13 , we derive for the first time a highly accurate formula for the densities of MRJ packings of frictionless two-dimensional (2D) binary (two-component) convex superdisks, which constitutes a family of an uncountably infinite number of distinct types of packings (both in particle shape and relative concentrations, see Fig. 1). Our formula, which is based on a novel organizing principle for such MRJ packings derived here, explicitly incorporates specific attributes of MRJ states, including jamming induced local spatial correlations and hyperuniformity, and thus, is distinctly different from the aforementioned schemes [43][44][45] . Our formula for frictionless binary superdisks is in excellent agreement with computer simulations (see Methods) for the entire α-x plane with semi-axis ratio α and small-particle relative number concentration x, and is easily generalized to other smooth non-spherical shapes and other dimensions.
It is noteworthy that for the special limit of monodisperse circular disks (p = α = x = 1), our prediction φ MRJ = 0.834 is in very good agreement with the recently numerically discovered MRJ isostatic state with a density of φ MRJ = 0.827 17 . The latter is a significant development because it was previously thought that either RCP states for monodisperse disks did not exist (because standard numerical protocols typically produce highly ordered, polycrystalline packings due to a lack of geometrical frustration) or, by entropic measures, RCP states were indeed these highly ordered polycrystalline arrangements. The fact that our prediction for φ MRJ in this instance is consistent with this new data is a testament to the power of the approach. Moreover, for most monodisperse superdisks that are not circles, our predicted MRJ densities are significantly lower than packing densities of polycrystalline packings that would be achieved by standard numerical and experimental packing protocols.

Results
Novel organizing principle for MRJ packings. The numerically generated high-quality strictly jammed packing configurations allow us to obtain robust statistics of the interparticle contacts. Specifically, we compute the average number of contacts per particle as a function of interparticle gap tolerance, which possesses a plateau over a wide range of tolerance values (see Fig. 2a). This feature enables us to determine the contact neighbors of a particle and to locate the contact points to a heretofore unachieved accuracy. The distributions of the number of contacts per particle for different packings are shown in Fig. 2b. Note that in these packings, the total number of contacts (i.e., constraints) N c , which can be easily computed from the contact number distributions, is smaller than the total number of degrees of freedom (DOF) N f = 3N + 3, which includes the contribution from the 2 translational DOF and 1 rotational DOF-associated with the N particles in the packing as well as the 3 DOF of the deformable boundary 16 . Unlike MRJ packings of circular disks, which are isostatic (i.e., N c = N f ), MRJ superdisk packings are hypostatic (i.e., with N c < N f yet strictly jammed). This indicates that the orientations of the contacting superdisks are necessarily correlated to block the relative rotations of the particles, which play a trivial role in MRJ packings of circular disks 46 .
We now devise an organizing principle by analyzing hypostatic contacting particle configurations. We observe that in MRJ packings, the contact points on a particle associated with small surface curvatures are favored over those associated with large surface curvatures. This feature is explained by the fact that the small-curvature contacts more efficiently block particle rotations 46,47 . However, the value of surface curvature depends on both particle shape and size, and thus, is different for superdisks with different p and α values. Therefore, we examine contact angles instead of contact curvatures for our analysis. Specifically, a contact angle is defined as the angle between the line connecting the particle center to a contact point on particle surface and the horizontal principal axis when the particle is aligned with the axes of a local Cartesian coordinate system. For a specific p value, the local particle surface curvature κ at the contact point is related to the contact angle θ via: where θ π ∈ ( , / 0 2 ] . Because the particle possess 4-fold rotational symmetry, the curvature value associated with other angles can be obtained by translating the angles to the interval π ( , / 0 2 ] . Since small-curvature contacts are more efficient in blocking relative particle rotations, it is reasonable to assume that for any specific p value, the most probable contact angle θ * (i.e., the mean of the distribution) corresponds to the minimal curvature κ min associated with the particle shape. Due to the symmetry of the particle shape, we expect the distribution of contact angles f(θ) is also symmetric about θ * . We thus propose that f(θ) possesses a Gaussian form, i.e., where the mean θ * and standard deviation σ depend alone on the particle shape (i.e., the deformation parameter p). To validate (3), we compute the statistics of contact angles from the numerically generated packing configurations. For each particle in a MRJ packing, the locations of contact points on the particle surface are identified and the associated contact angles are computed. The contact angles for all particles in the packing are then binned to generate a histogram, or a distribution of the angles. We test the validity of the distribution using two different methods (i.e., Kolmogorov-Smirnov test and t-test) and find that for any given particle shape that we considered, the distribution of contact angles indeed possesses the Gaussian form (3). In addition, after simple re-scaling, the contact angle distributions for different particle shapes all collapse onto a universal standard Gaussian distribution curve with zero mean and unitary standard deviation, as shown in Fig. 3a. We note that the original mean value of contact angle for superdisks with p < 1 is π/4 and the mean for superdisks with p > 1 is π/2. Figure 3b shows the standard deviations of the un-scaled contact angle distributions as a function of p. This collapse strongly suggests the existence of a universal organizing principle for MRJ packings of frictionless superdisks independent of the detailed particle shape, which we elaborate on below.
As shown in Fig. 3b, the standard deviation σ is a monotonic decreasing function of − p 1 and possesses a singularity at p = 1 (i.e., the circle limit). This indicates that as the shape deviates from that of a circle (with uniform local curvature), the contacts are more likely to occur in low curvature regions on the shape, which results in significant orientational correlation in the packing. For nonspherical shapes, the requirement of jamming in the packing (i.e., mechanical stability) necessarily leads to correlations among the contacts, which block both the translational and rotational motions of any particles in the packing. The observed correlations among contacts, which are quantitatively characterized by the universal Gaussian distribution of contact angles, are sufficient and necessary to induce jamming and thus, are an attribute of MRJ states. Therefore, the capability of blocking relative particle motions via contacts at different local surface curvatures associated with different contact angles does not depend on the specific particle shape and spatial dimension. This is not only a new organizing principle for MRJ packings of frictionless superdisks but other frictionless nonspherical particles with a sufficiently smooth shape (e.g., ellipses). A quantitative consequence of this organizing principle is the universal Gaussian distribution of contact angles. Moreover, an accurate analytical formalism for predicting φ MRJ should quantitatively incorporate such jamming induced spatial correlations, as we do here.
Accurate prediction of MRJ packing densities. We now employ our organizing principle to estimate the MRJ packing density φ MRJ for frictionless superdisk systems for various values of p, x and α. The salient idea is to approximate the global packing density by a contact-configuration averaged local packing density. We first provide arguments why such a local analysis can lead to accurate predictions of φ MRJ . We note that the MRJ packings are hyperuniform with long-range correlations 13,18 . A unique feature of these hyperuniform systems is that the local-volume-fraction fluctuations decay faster than the inverse of the volume of a spherical observation window of radius R, i.e., faster than R −d 21 . This indicates that the local packing density τ(R) converges to the global packing density φ MRJ rapidly as the size of the window increases, i.e, effectively σ In other words, the long-range hyperuniformity property imposes strong suppression of density fluctuations at short scales. This rapid convergence property of disordered hyperuniform jammed packings is supported by previous theoretical studies 13,18,21 as well as experimental investigations 27 . Moreover, to further reinforce this point, we have computed here the variance σ ( ) τ R 2 of a MRJ packing of binary superdisks with R equal to the diameter of the large particles, and compared it to that of a corresponding saturated non-hyperuniform packing. This non-hyperuniform packing is generated by first randomly and sequentially placing particles with random orientations in a simulation domain without overlapping existing particles, i.e., via the random sequential addition (RSA) process, which leads to a nearly hyperuniform packing 48 . Then the RSA configuration of superdisks is compressed to the densest possible state while fixing the particle orientations. The variance of the MRJ packing in this small window is four times smaller than that of the compressed RSA packing. Therefore, although only first-neighbor shell is considered in our analysis, the organizing principle is sufficient to yield accurate prediction of φ MRJ .
Our density-prediction procedure works as follows: we first construct a local two-particle packing configuration, which is a region Ω enclosed by two line segments that are tangential to two contacting particles and portions of particle surfaces, as shown by the thick dash lines in Fig. 4. We note that this construction allows us to quantitatively incorporate the jamming-induced contact distribution, which encodes explicit two-body statistics. Our approach is to be contrasted with the "granocentric" approach 43 as well as the Edwards' ensemble-approach based on the Voronoi volume distribution 44,45 , which is a one-body statistic 49 . In the "granocentric" approach, which was originally devised for 3D colloidal spheres, the formation of an interparticle contact was considered independent of other contacts. However, extending this assumption to 2D is problematic because of the lack of intrinsic geometric frustrations (intrinsic to 3D packings), a consequence of which are substantial correlations between contacting particles. For the Edwards' ensemble-approach, although the Voronoi volume of a particle depends on the centroid positions of its contacting neighbors, unlike in our formalism, the spatial correlations between the central particles and its neighbors were not explicitly incorporated. We then compute the local packing density φ local = v particle /v local , where v particle is the total volume of the two particles and v local is the volume of the region Ω . In other words, we can consider that the 2D plane is approximately tessellated into polygons, each effectively containing two particles.
In general, φ local depends on the location of contact on the particles, which is characterized by the contact angle θ, the relative orientation of the particles Θ and size ratio α between the two contacting particles. For a binary packing (the focus of this letter), there are two classes of local contact configurations, i.e, the two particle are of the same size (whose local packing density is denoted by φ ( ) local 1 ) and two particles with different sizes (with local packing density denoted by φ ( ) local 2 ). Therefore, the global MRJ density φ MRJ can be estimated as follows: is the average local density over Θ at fixed α. For p < 1, θ 1 = 0 and θ 2 = π/2; for p > 1, θ 1 = π/4 and θ 2 = 3π/4. More generally, formula (14) is easily extended to the case of a polydisperse packing with particle size distribution P(R) (where R is the semi-axis of the superdisks) as follows: local is the average local density over Θ and α for a given superdisk with semi-axis R.
To verify its accuracy, Eq. (4) is employed to estimate the φ MRJ for a variety of binary frictionless superdisk packings with x ∈ [0.05, 0.95], α ∈ [0.2, 0.95] and p ∈ [0.85, 4.5]. We find that for all values of α and x considered the estimated φ MRJ always agrees very well with the corresponding simulation results, with the largest deviations smaller than 1.5%, as shown in Fig. 5 for selected values of x, α and p (see the Supplementary Information for the values of φ MRJ ). For the extreme case where p = 0.5 and ∞, the particle becomes a square with singular curvatures at the corners and our formalism does not hold. In addition, when α is very small (e.g., less than 0.2), even with x ∈ [0.05, 0.95], the large particles can form a jammed backbone with small particles moving freely within cases formed by large particles 50 , and thus, the organizing principle does not hold for these cases.

Discussion
When the particles are nearly identical in size (α ~ 1) or when the concentration of one species is overwhelmingly larger than the other (x ~ 0 or 1), the most "disordered" packings generated using typical packing simulation protocols, including the one employed here (i.e., DTS molecular dynamics method), possess polycrystalline regions and thus, a high degree of order. This calls into question whether our simulated packings, in these special cases, represent the actual MRJ state as well as the effectiveness of numerical packing protocols in generating true MRJ states. In ref. 16, high-fidelity MRJ isostatic packings of equal-sized circular disks were constructed for the first time using a novel sequential linear programming (SLP) packing protocol 51 . Such packings were shown to possess no crystalline regions and an average packing density φ MRJ = 0.827. This singular result provides a stringent benchmark for our predictions of φ MRJ for equal-sized superdisks. Specifically, in the limit α = 1 and p = 1, our prediction yields φ MRJ = 0.834, which agrees very well with the numerical result reported in ref. 16. This consistency strongly indicates the accuracy of our formula for the circle case and it is reasonable to expect it yields accurate predictions of φ MRJ for equal-sized superdisks with p ≠ 1. Indeed, as shown in Fig. 6, the predicted φ MRJ for these shapes are at least 5% lower than the corresponding polycrystalline packings produced using the DTS simulations. This provides strong evidence that such polycrystalline packings of superdisks with p ≠ 1 do not represent the actual MRJ states, and our predicted φ MRJ correspond to MRJ packing arrangements that are yet to be constructed via certain novel numerical and experimental protocols.
We expect that our formalism yields accurate predictions for φ MRJ for frictionless superdisks with x ∈ [0, 1], α ∈ [0.2, 1] and p ∈ (0.65, 5), which are beyond the parameter values studied in our numerical simulations. The reason is that it explicitly incorporates characteristics of MRJ states, including hyperuniformity-and jamming-induced local correlations between contacting neighbors, which are crucial ingredients in a novel local organizing principle we report here. Moreover, the derivation based on the novel organizing principle does not depend on specific particle shape (i.e., p values) and size distributions. It is notable that our analysis provides a clear systematic procedure to predict MRJ densities of other families of smooth particle shapes in both 2D and 3D, e.g., ellipses and ellipsoids. We note that after initial submission of our paper, we learned of a preprint that uses a local perturbation analysis to predict random close packing densities of the special case of smooth particle shapes that are very nearly spherical, which can continuously be deformed into a sphere 52 .

Methods
We generate via the molecular-dynamics-based Donev-Stillinger-Torquato (DST) method 53 high-quality strictly jammed packing configurations and then analyze the characteristics of individual packings in order to deduce a novel organizing principle. The packings of frictionless binary superdisks with x = n S /(n S + n L ) = 2/3 (where n L and n S are respectively the number of large and small particles) and α = R S /R L = 5/7 (where R L and R S are the corresponding semi-axis of large and small particles, Figure 6. Comparison of the MRJ packing density for monodisperse superdisks estimated using Eq. (4) (solid line) and obtained from DTS molecular dynamics (MD) simulations (red circles), which fail to represent MRJ states, as well as sequential linear programming (SLP) datum for the circle case (blue square) 17 . respectively). Specifically, small non-overlapping particles are initially placed in the periodic simulation box (fundamental cell) with random positions and orientations. The particles are then given random translational and rotational velocities and their motions follow Newtonian dynamics as they collide elastically and also expand uniformly with an expansion rate Γ , while the fundamental cell deforms to better accommodate the configuration. Large Γ values (~0.01) are used in the early stage of the simulation to suppress the formation of order; and very small Γ values (~10 −6 ) are used when the system is close to jamming in order to establish well-defined contact networks. The final packings are verified to be strictly jammed using the "infinitesimal shrinkage" method with deformable boundary 47 .