Cuddling Ellipsoids: Densest local structures of uniaxial ellipsoids

Connecting the collective behavior of disordered systems with local structure on the particle scale is an important challenge, for example in granular and glassy systems. Compounding complexity, in many scientific and industrial applications, particles are polydisperse, aspherical or even of varying shape. Here, we investigate a generalization of the classical kissing problem in order to understand the local building blocks of packings of aspherical grains. We numerically determine the densest local structures of uniaxial ellipsoids by minimizing the Set Voronoi cell volume around a given particle. Depending on the particle aspect ratio, different local structures are observed and classified by symmetry and Voronoi coordination number. In extended disordered packings of frictionless particles, knowledge of the densest structures allows to rescale the Voronoi volume distributions onto the single-parameter family of $k$-Gamma distributions. Moreover, we find that approximate icosahedral clusters are found in random packings, while the optimal local structures for more aspherical particles are not formed.


I. INTRODUCTION
"How many candies can I fit in the jar?" Children and scientists alike have always been fascinated by packing problems [1]. Nevertheless, it remains a challenge to understand how collective properties of packings arise from microscopic mechanisms on the particle level. Microscopic interactions are encoded in local structural motifs comprising a particle and its immediate neighbors. These structural motifs serve as the building blocks for extended packings. A pivotal role in the analytic modeling of granular matter is played by the densest local structure [2,3]. This structure coincides with the solution of the kissing problem in mathematics [4,5], which asks to maximize the number of spherical particles simultaneously in contact with a central one. The icosahedral cluster, depicted in Fig. 1, top row, is the most symmetric way to arrange the twelve kissers and maximizes the local packing density. Ideal icosahedra are about 1% denser than the best space-tiling arrangements of congruent spheres, given by stacked hexagonal lattice planes [6]. By embedding distorted icosahedral clusters, disordered granular packings also locally exceed the density limit for space-tiling sphere packings.
Much like candies, the particles found in nature and in industrial applications are not always spherical and often randomly shaped. For example, pebbles and sand grains vary widely in size and shape, and the kissing and packing problems lack rigorous answers. Generalizations of the sphere-packing problem to congruent aspherical particles have been intensely studied [7][8][9][10][11], motivated both by applications in granular matter [12][13][14][15][16][17] and by advances in the synthesis of colloidal particles with prescribed shapes [18,19]. In general, aspherical shapes pack denser than spheres, and it could recently be shown that the sphere is a local pessimum for lattice packing [20,21]. A convenient shape for studying the effect of asphericity * sebastian.kapfer@fau.de are ellipsoids. Even though dense ellipsoid crystals are known [7], there is currently no mathematical proof of optimality.
The kissing problem has been generalized to tetrahedra [22], but we are not aware of results for ellipsoidal particles. Motivated by recent work on the packing properties of disordered granular ellipsoids [12][13][14][15], we here consider a modified kissing problem for uniaxial ellipsoids 1 : 1 : α, also known as spheroids, or ellipsoids of rotation. The number α is the aspect ratio of the particles, that is, α < 1 (α > 1) corresponds to oblate (prolate)  The solid line at the bottom of the graph represents the densest-known crystal structure, and "ico" is the icosahedral cluster. Usually, in the optimal structures, the contact number C is equal to the number of Voronoi neighbors N . One exception is the oblate N = 14 structure which loses contacts around α ≈ 0.85 (see kink). A detailed description of the structures is given in the text.
ellipsoids. The conventional kissing problem maximizes the number C of particles in contact with a given central particle. Here, we instead maximize the local packing density φ l = V α /V , defined via the volume V of the central particle's Voronoi cell, where V α = 4πα/3 is the particle volume. The appropriate generalization of the Voronoi diagram for aspherical particles is the Set Voronoi diagram [23], also known as navigational map [24,25]. The Set Voronoi diagram allocates space by the distance to the particle surfaces instead of the distance to particle centers. Typical Set Voronoi cells are non-convex and have curved surfaces. The solution of our modified kissing problem is the structure minimizing the Voronoi volume V . We define V min (α) as the minimal Voronoi cell volume for each aspect ratio, and φ max (α) = V α /V min (α) as the corresponding local packing fraction. As a first key result of the present article, we numerically determine the optimal local structures for uniaxial ellipsoids of aspect ratios 0.7 ≤ α ≤ 1.4.
The description of granular matter and other out-ofequilibrium particle systems by macroscopic variables in the spirit of thermodynamics currently is a nascent field. A convenient characterization of the particle system is given by Voronoi cell partitions. In particular, distributions of Voronoi cell volumes are sensitive to structural transitions in granular assemblies [15,[26][27][28][29][30][31][32]. The theoretical modeling of such distributions remains an open problem. For non-interacting particles (ideal gas, Poisson point process), the distribution of Voronoi volumes is almost perfectly fit by a three-parameter Gamma distribution [33,34]. No fitting curve of comparable accuracy exists for interacting particles, let alone dense packings such as considered here. Aste et al. [2,3] propose a simplified model for sphere packings. This analytic model yields the full distribution of the Voronoi volumina, but crucially depends on the minimal cell volume V min (α). Our results for V min (α) permit to extend this model also for random packings of ellipsoidal particles.
In the following section II, we discuss the numerical optimization procedure and the resulting densest local structures. Section III investigates the occurrence of these motifs as building blocks of extended random packings. Moreover, we test the analytic model of Aste et al. for ellipsoidal particles and discuss future directions.

II. DENSEST LOCAL STRUCTURES
Finding the densest local structure of ellipsoids amounts to locating the minimum in a complex energy landscape, with the energy given by the Voronoi cell volume V . The function V depends on the position and orientation of the N Voronoi neighbors of the central particle. In order to solve this optimization problem in 5N variables, we use a two-stage approach: First, candidate structures are explored by a simulated annealing scheme; second, the densest candidate structures found during the annealing phase are optimized by a downhill algorithm. The numerical optimization considers only a central (red) particle and the N first-shell neighbor particles that share a Voronoi facet with it. Since the optimal number N is not known a priori, we propose several values of N in separate annealing runs; during one run, the number N is fixed. While the red particle is immobile, the other particles can translate and rotate. The mobile particles have a tendency to drift away and detach from the central Voronoi cell. We thus restrict their centers to an ellipsoid-shaped arena around the red particle. Each trial move proposes to displace and rotate a single mobile particle. The move is rejected if any overlap between the particles is created. Otherwise, we compute the volume change ∆V of the Voronoi cell belonging to the central particle, and accept the move with probability exp(−κ∆V ). The dimensionless pressure κ is 50 initially, and increases to 1000 over the 2 × 10 6 Monte Carlo steps comprising a run. We keep track of the best (densest) configuration observed so far, which is restored at the end of the annealing run. Finally, we use a downhill algorithm to find a local optimum of the cell volume.
The key feature of the simulated annealing algorithm is its ability to escape from local optima, but it does not always return the global optimum. We thus perform ten independent annealing runs for each set of parameters (N, α). In cases where multiple structures are in close competition, the number of independent runs was increased to 100. The major computational expense of our optimization procedure is the computation of the Set Voronoi cell volume. To be able to perform sufficiently many trial moves in the optimization, the resolution of the Set Voronoi volume computation must be limited. A typical optimization run uses a discretization [35] of 126 Voronoi seeds per particle (320 for the downhill phase). One evaluation of the cell volume takes ≈ 10 ms (up to 50 ms in the downhill phase) on a Xeon E5-2630 CPU. For the final configurations at the end of the downhill phase, we also compute high-precision cell volumes, and estimate the discretization error from the difference (see error bars in Fig. 2).
As the result of our numerical optimization scheme, we identify the presumed optimal structure for each aspect ratio α and Voronoi coordination number N , see Fig. 2. In many of the optimal structures, neighbor particles are related by discrete symmetry operations. These groups of neighbors share the same color in all figures. The colored boxes indicate the number of particles in a ring, with the same color code. Black and white belts mark the equator of each ellipsoid. The distinguished axis of the central particle is indicated by the red line. To obtain more precise coordinates, we also perform constrained simulations which enforce the rotational and mirror symmetries inferred from some of the discovered structures (see below). The reduced degrees of freedom lead to improved convergence of the simulated annealing algorithm [36].
As the particles become more aspherical, a general trend towards denser structures and increased N exists (see Tab. I, row II). In this sense, the sphere is not only a worst case of the lattice packing problem [20,21] but also of the local packing problem. For practical reasons, we limit our study to moderately aspherical ellipsoids. First, the pointed tips of very aspherical particles necessitate a fine discretization in the Set Voronoi computation. Second, dense structures of very aspherical objects strongly depend on minute details of the particle shape [10] and are less generic than the results for moderate asphericity considered here.
In the densest local structure of spheres [37], the symmetric icosahedral cluster (ico; see Fig. 1, top row), the first-shell neighbors do not touch each other, leaving space to be filled. In fact, the twelve kissing spheres leave so much space that it is possible for them to switch places without ever losing contact with the central sphere. Thus, packing density can be increased in two ways: One can either expand the neighbors in a direction parallel to the central sphere's surface, leading to prolate ellipsoids, or squash the neighbors in the normal direction, producing oblates. Thus, for small asphericity, oblate and prolate ellipsoids prefer different orientations of neighbors with respect to the central particle's surface, and pack densest in two different structures, see Fig. 3. The two branches cross in the ico configuration at α = 1.
Moving away from the spherical case, the icosahedral symmetry is broken. Even though the neighbor particles twist, a threefold rotation axis is preserved both for oblates and prolates, and the densest structure consists of a succession of four three-rings (3-3-3-3 structures). While both oblates and prolates show three twofold axes, the oblates have dihedral mirror planes (point group D 3d ), and prolates preserve an improper sixfold axis (S 6 ). The evolution of these structures with α is illustrated in Fig. 4. The icosahedral cluster also features fivefold axes, which would imply a 1-5-5-1 structure, see Fig. 1 (top right). The 1-5-5-1 structure is degenerate with the 3-3-3-3 for sphere-like oblates, see gray and yellow curves in Fig. 2. For all other aspect ratios, the fivefold symmetry is not stable with respect to unconstrained global optimization.
In the bottom row of Fig. 1, we illustrate two ways to generalize the ico cluster towards higher N : to introduce additional rings of particles (left), or to increase the number of particles per ring (right). Most of the dense structures found in this study can be categorized in this fashion, see Fig. 2. For example, prolate ellipsoids with 1.26 α 1.42 pack densest in a structure with threefold symmetry, composed of five three-rings of particles (3-3-3-3-3, point group D 3 ). At higher α, a structure with fivefold symmetry and N = 17 is preferred (1-5-5-5-1, C 5h ). For some of the structures, the rotational symmetries of the particle rings are broken. Consider for example the prolate 1-2-4-4-2-1 structure with N = 14 neighbors, which is optimal in the range 1.19 α 1.26. In this structure, the rings are not stacked along the central particle's axis. Consequently, the rotational symmetry is lifted by the central particle, and only two mirror planes are preserved (C 2v ).
Oblate particles below α ≈ 0.86 form a 14-coordinated structure instead of the icosahedral cluster, see Fig. 2. The full sixfold symmetry is reached only for α 0.85 (1-6-6-1, D 6 ) where all neighbors touch the central particle (contact number C = 14). Above α ≈ 0.85, due to steric constraints, not all of the Voronoi neighbors can touch the central particle, hence C < 14. The sixfold symmetry is lost and the packing fraction of the N = 14 structures declines rapidly. Once the packing fraction falls below the N = 12 branch, the cluster tends to expell a particle (orange particle in Fig. 2). For very oblate ellipsoids, α 0.76, fifteen neighbors without any obvious symmetry pack densest around the central particle (labeled N = 15 in Fig. 2 and Fig. 4). We conjecture that for even larger asphericities than considered here, the densest structures are typically disordered without any symmetries.  5. a) Distributions of the rescaled variable (φ l − φg)/σ(φ l ). For reference, a Gaussian distribution is shown (black dotted curve). b) Behavior of the standard deviation of local packing fractions σ(φ l ). The gray points are reference data from Ref. [15].

III. DENSE RANDOM PACKINGS
Having established the densest local structures for each aspect ratio, we proceed to analyze the occurrence of these motifs in dense disordered packings. We generate, for each aspect ratio, at least 100 random packings of 5000 monodisperse ellipsoids each using the Lubachevsky-Stillinger-Donev protocol [38,39]. The packing generation protocol does not implement gravity, and consequently the packings are essentially free of orientational order [40]. We set the expansion rate of the particles to 3 × 10 −6 · · · 10 −5 times the thermal velocity, which produces dense packings, but avoids the eventual formation of crystalline domains. Thus, our packings are representative of the random close-packed (or maximally random jammed) state for the respective aspect ratio [12,13]. Our data is consistent with published tomography experiments and numerical sedimentation data for oblate ellipsoids [15]. Fig. 5a shows the distributions of local packing fraction. Due to high statistics, we can resolve the fringes of the distribution which are not yet accessible in experiment. The standard deviations σ(φ l ) are also in very good agreement with earlier data for oblate ellipsoids (gray data points in Fig. 5b). Moreover, Fig. 5b demonstrates a surprising symmetry between ellipsoids of reciprocal aspect ratios: Both global packing fractions φ g = 1/ 1/φ l and the standard deviations σ(φ l ) approximately coincide. This finding continues a series of unexplained correspondences between oblate and prolate ellipsoidal particles, for example the equilibrium phase diagram [9] and packing fractions φ cry of the best known crystals [7]. Fig. 6 displays the frequency of Voronoi cells in our dense disordered packings with a given coordination number N and packing fraction φ l . The bold vertical bars indicate, for each Voronoi coordination number, the maximal packing fraction found in Fig. 2. At small asphericity, the densest local structure is the icosahedral cluster. Random packings of spheres are known to con- tain distorted variations of icosahedral clusters [41], while the probability for perfect ico clusters vanishes [30]. Our data for α = 0.9, 1.0, 1.1 in Fig. 6 confirms these results for packings of moderate asphericity. More aspherical ellipsoids (α ≤ 0.8 and ≥ 1.25) could pack denser with fourteen or more neighbors. As Fig. 6 demonstrates, such structures are not formed in significant amounts. Instead, the densest clusters are again N = 12 cells. Visual inspection and a quantitative analysis using Minkowski structure metrics [42] (see appendix) confirm that very aspherical particles tend not to form their optimal structures, but do form distorted variations of the N = 12 optimal structures. The densest clusters in a random packing readily exceed the best known ellipsoid crystals of packing fraction φ cry . These supracrystalline clusters (φ l > φ cry ) are almost uniformly distributed in the packing. The presence of a supracrystalline structure implies reduced packing fraction in its vicinity as these motifs  cannot be periodically continued. The radial distribution function g(r) of supracrystalline clusters is reminiscent of a hard-core fluid with short-range repulsion and quickly decays to unity, see top right plot in Fig. 6. We return to this point in the conclusion. While Voronoi diagrams are a powerful tool for characterizing granular packings, there is currently no established theory for aspherical particles. Aste et al. propose an analytic model [2,3] which predicts, under simplifying assumptions, the full distribution of Voronoi cell volumina. They find a so-called k-Gamma distribution, with probability density Γ k (x) ∝ x k−1 exp(−kx), where the quantity x is the re-scaled volume For jammed spheres, the parameter k varies between 11 and 15, which agrees with typical neighbor counts [2]. The value of k has been linked to a granular temperature [3]. The shape of the k-Gamma distributions strongly depends on V min (α), which was previously unknown for ellipsoids. Combining our present results for V min (α) and the k-Gamma model, we can now predict the full distribution of Voronoi volumina, with only a single parameter k. Fig. 7a contrasts our data for dense random packings with the predictions from the k-Gamma model. Both for spheres and ellipsoids, the model curves satisfactorily reproduce the skewed distributions of the data. Upon closer inspection, however, the data exhibit a shoulder at small volumes which deviates from the k-Gamma curve. The position of the shoulder is lowest for spheres, and shifts into the main peak (larger x) for aspherical particles. A similar feature is also present in packings of bidisperse disks (see Fig. 2 and 3 in Ref. [28]) and randomclose-packed colloidal particles (see Fig. 4 in Ref. [43]).
In order to isolate the origin of the shoulder contribution, we fit an unnormalized k-Gamma probability density, m 1 Γ k (m 2 x), to the upper portion of the data (x > 0.6 or x > 0.8, depending on α). The majority of Voronoi cells, m 1 /m 2 ≥ 95%, can be attributed to the k-Gamma-distributed background. For our packings, the best-fit value for k is between 12 and 17, see Tab. I, row III. In general, k increases with asphericity, while the average Voronoi neighbor count stays almost constant, N ≈ 14. Fig. 7b shows the residual after subtracting the k-Gamma background from the Voronoi volume distributions. We find an excess of Voronoi cells at low volumes, indicating that there are 'attractor' motifs at these packing fractions which are preferentially formed. This excess mainly stems from N = 12 cells which show a bimodal φ l distribution (see Fig. 6, bottom row). Evidently, the attractors which cause our packing to deviate from the k-Gamma model are the N = 12 motifs for the relevant aspect ratio. Excluding such specific packing motifs, the k-Gamma model provides useful and accurate parametrization of dense frictionless ellipsoid packings.

IV. CONCLUSION
In the present article, we explored the influence of asphericity α on local packing motifs. Opening a new chapter in the story of the classical kissing problem, we propose a modification for arbitrarily-shaped particles which is tractable by numerics. We present the results of an extensive numerical search for optimal structures. Surprisingly, many of our densest structures exhibit a high degree of symmetry. With increasing asphericity, we find a trend towards higher coordination numbers and packing fractions. Even though we do not maximize the contact number, in most of our optimal structures, all neighbor particles touch the central one, i. e., C = N . One exception is the oblate N = 14 structure which loses contacts for particles too spherical (see kink in Fig. 2). Our optimal structures also are good candidates for solving the classical kissing problem for ellipsoidal particles. A numerical approach similar to ours was recently employed to study clusters of Platonic solids with interesting physical properties [44]. The continuous degree of freedom present in our ellipsoid clusters could be useful in realizing new kinds of symmetries, and tailoring the structure to applications.
Furthermore, we connect these densest motifs to the properties of disordered ellipsoid packings. Knowledge of the minimal cell volume V min (α) permitted us to rescale our data of aspherical particle packings onto the curves predicted by the k-Gamma model. The rescaled distributions of Voronoi volumina are well described by the model curves. The remaining small deviations can be explained by the excess formation of certain types of N = 12 clusters. The success of the k-Gamma model is unexpected as it includes correlations between neighboring particles only via the minimum packing volume. In particular, mechanical stability of a packing would also imply an upper cutoff V max for the Voronoi cell volume distribution, not included in the model. V max would also depend on additional microscopic parameters such as friction and hence would be even more difficult to establish than V min . We expect that V max will be essential for the description of loose packings; in the present random-close-packed limit, low compactivity effectively masks this effect. The local packing motifs studied here can be regarded as the building blocks of disordered packings, and are analogous to the atoms in a solid. The interaction between neighboring 'atoms' is of similar complexity as the energy landscape of glasses, and the treatment of disordered phases in granulates, glasses and complex fluids remains a challenging task. However, the frequency and spatial distribution of supracrystalline motifs may provide some insight into the hidden structure of packings. Interestingly, only variations of icosahedral clusters are preferentially formed in our random packings. The optimal structures of our more extreme ellipsoids are not realized in random packings, an effect which is increasing with asphericity. One could imagine that slight polydispersity might amplify the formation of optimal motifs. Supracrystalline clusters exhibit nontrivial correlations on length scales above the particle scale (see inset in Fig. 6) which require further investigation. It will be interesting to see whether and how the history of packing formation, such as sedimentation under gravity or compaction by shaking are reflected in these signatures, and what additional information they can reveal about the architecture of the packing.
As demonstrated here, ellipsoids are a useful testing ground for the effects of variations in grain shape. Our findings are immediately relevant for the realistic modeling of granular materials consisting of aspherical, polydisperse and randomly-shaped grains. Such materials are found in geological processes such as the dynamics of dunes or tali, and industrial applications. There is no trend towards the 1-6-6-1 structure (green bullet).