Jammed Spheres: Minkowski Tensors Reveal Onset of Local Crystallinity

The local structure of disordered jammed packings of monodisperse spheres without friction, generated by the Lubachevsky-Stillinger algorithm, is studied for packing fractions above and below 64%. The structural similarity of the particle environments to fcc or hcp crystalline packings (local crystallinity) is quantified by order metrics based on rank-four Minkowski tensors. We find a critical packing fraction \phi_c \approx 0.649, distinctly higher than previously reported values for the contested random close packing limit. At \phi_c, the probability of finding local crystalline configurations first becomes finite and, for larger packing fractions, increases by several orders of magnitude. This provides quantitative evidence of an abrupt onset of local crystallinity at \phi_c. We demonstrate that the identification of local crystallinity by the frequently used local bond-orientational order metric q_6 produces false positives, and thus conceals the abrupt onset of local crystallinity. Since the critical packing fraction is significantly above results from mean-field analysis of the mechanical contacts for frictionless spheres, it is suggested that dynamic arrest due to isostaticity and the alleged geometric phase transition in the Edwards framework may be disconnected phenomena.

deduced from the fact that the volume fraction of polytetrahedra increases with packing fraction up to φ ≈ 0.646 and then decreases, as these structures transform into crystalline order [13,14].
In addition to these phase transition scenarios, where φ RCP is interpreted as the density of the disordered phase at coexistence, the notion of the maximally random jammed state (MRJ, [15]) has been proposed (as the maximally disordered state amongst all jammed packings); with respect to a number of common measures of order, the MRJ packing fraction has been estimated as φ MRJ ≈ 0.63 [7].
The resolution of the RCP problem relies on suitably defined order metrics to quantify packing structure. A common approach to local structure characterization is by analysis of nearest neighborhoods [16,17]. Often, the bond-orientational order metrics q l defined by Steinhardt et al. [16] are applied to sphere packings [7,18]. However, these and other neighborhood-based order metrics [19] have shortcomings. First, they suffer from the ambiguous definition of the nearest neighborhood [20]. Second, in their common use as single scalar order metrics [7,18], they are insufficient to conclusively distinguish order and disorder [21]. A non-negligible fraction of non-crystalline environments is often incorrectly identified as crystalline (false positives), since their q 6 are close or identical to the reference q 6 values in crystals.
Alternatively, the structure of monodisperse sphere packings can be characterized by analysis of the Voronoi cells of the particle centers, see Fig. 1. Suitable morphological descriptors, such as Minkowski tensors [22,23], can then be used to quantify the cell shape and hence the local structure. Here, we show that crystalline order metrics can be constructed from rank-four Minkowski tensors of the Voronoi cells that give stringent criteria for fcc or hcp crystalline order. For jammed sphere packings (Color online) (left) The Voronoi diagram of a packing associates convex cells with the individual particles. Each Voronoi cell K contains one particle k such that the distance of any point p ∈ K to particle k is smaller than the distance of p to any other particle. (right) For the evaluation of q6, we define the set of nearest neighbors of particle k as those 12 other particles (6 in 2D) which are closest to k. generated by the LS algorithm, these order metrics reveal an abrupt onset of crystallinity at a critical packing fraction φ c ≈ 0.649.
Eigenvalue ratios of rank-two Minkowski tensors quantify anisotropy of the particle environments in jammed bead packs [22]. The Voronoi cells of the fcc and hcp close packing are isotropic, in the terminology of Ref. [22], while cells found in disordered packings typically are not. However, rank-two tensors are insufficient to distinguish different types of isotropic cells. These cell types can be classified using the rank-four Minkowski tensor W 0,4 1 . The tensor W 0,4 1 of a Voronoi cell is given as the sum of tensor products of the facet normals, weighted by the facet areas, where n i := (n(f )) i with i = 1, 2, 3 are the cartesian components of the facet normal, and a(f ) is the surface area of facet f ; further, A := f a(f ) is the total surface area; all sums run over the facets of the Voronoi cell K. In close analogy to the stiffness tensor of continuum mechanics, symmetry under permutations of indices allows the reduction of W 0,4 1 to a 6 × 6 symmetric matrix [24]. The six eigenvalues (ς 1 , . . . , ς 6 ) of this matrix are dimensionless due to normalization by A −1 and are rotational invariants [25]. A concise quantitative measure for the similarity of a given Voronoi cell K to the Voronoi cell K fcc of a crystalline fcc packing is given by the fcc crystalline order metric An analogous order metric ∆ hcp is defined for hcp cells. Figure 2 supports our claim that ∆ fcc and ∆ hcp measure deviations of a Voronoi cell's shape from the ideal fcc or hcp cell. The vertices of an ideal fcc or hcp lattice are displaced by small random vectors. Figure 2 shows averages and standard deviations (as error bars) of ∆ fcc and ∆ hcp as function of the root mean square displacement (RMSD) from the ideal lattice points, demonstrating an approximately linear relationship between deviations from the ideal crystalline shapes and the crystalline order metrics ∆ fcc and ∆ hcp . The quantitative agreement between the functions for hcp and fcc cells justifies the use of the same threshold value for selecting both fccand hcp-like cells from a packing. The crystalline order metrics ∆ fcc , ∆ hcp are used to identify crystalline clusters in jammed sphere packings generated by the Lubachevsky-Stillinger protocol [5]. Figure 3 shows, as key result of this paper, that (a) crystalline fcc and hcp order is absent for packing fractions below a critical value φ c ≈ 0.649; and that (b) above φ c the fraction of crystalline fcc or hcp in LS simulations is non-zero and rapidly increases by several orders of magnitude. As expected for an athermal system, no systematic difference between the number of hcp and fcc cells is observed, in contrast to crystallization dynamics in equilibrium [26].
We measure the fraction of fcc cells as n fcc (φ) = [N fcc ] 0.5 /N , where N fcc is the number of cells with ∆ fcc < δ = 0.005, in a sample of N = 4 × 10 4 spheres. [N fcc ] 0.5 is the median of N fcc over M ≈ 20 packings of similar φ, see Fig. 3. (In general, for a random variable X with a probability density f (X), the symbol [X] p denotes the p-quantile, i.e. the value F −1 (p) with the cumulative distribution function F (X) = X −∞ f (ξ)dξ.) The accuracy of our results is, for small n fcc and n hcp , limited by the finite system size, preventing the measurement of probabilities smaller than 1/N . (We further note, that the number of isotropic cells, in the terminology of Ref. [22], vanishes below φ c , and becomes nonzero at φ c .) The values of n fcc and n hcp depend, of course, on the choice of the threshold δ. The increase in local crystallinity is, however, also evident in the lowest occurring values of ∆ fcc . Figure 4 shows the first percentile [∆ fcc ] 0.01 as a robust estimate for the most crystal-like cell, i.e. the lowest occuring value of ∆ fcc . A sharp drop of [∆ fcc ] 0.01 is observed for φ φ c ; in packings below φ c , the most fcc-like cells are substantially different from fcc cells, while above φ c the differences quickly decay to close to zero. The value of φ c (x) is estimated by the intersection of two straight lines fitted to the data for [∆ fcc ] x . The insert of Fig. 4 shows the φ c estimates extracted by this approach, giving φ c ∈ [0.6492, 0.6499] for x ∈ [0.001, 0.1]. These values of φ c are larger than published values for the RCP or the MRJ packing fraction. Importantly, the data of Fig. 4 demonstrate that the drastic increase in n fcc in Fig. 3 is a robust result that is not sensitive to the value of the threshold δ. The value of δ is, within bounds, an irrelevant parameter. We do not observe differences between packings of N = 10 4 and N = 4 × 10 4 particles besides improved statistics.
The observed abrupt appearance of crystalline cells at φ c is difficult to observe using the bond-orientational order metrics q l and w l developed in seminal work by Steinhardt et al. [16]. Most frequently, q 6 is considered, which is deemed particularly sensitive to formation of fcc; it is Song et al. [11] Anikeenko et al. [13] FIG  with the spherical harmonics Y lm , the polar angles θ jk and ϕ jk of the bond vector between particles j and k, and · denoting the average over the 12 closest neighbors j of k (Fig. 1). Figure 5 shows probability distributions of q 6 values observed in LS packings above φ c that demonstrate the principal deficiency of using only q 6 as an order metric. The frequency f (q 6 ) of q 6 values develops sharp peaks at the values corresponding to fcc (q 6 = 0.57452) and hcp (q 6 = 0.48476) for φ φ c , not present for samples with φ < φ c . These peaks, however, sit on top of a dominant background of non-crystalline cells. The data clearly show that cells (false positives) exist which are distinctly different from fcc but that are identified as fcc by q 6 , i.e. |q 6 − q fcc 6 | < 5 · 10 −4 . For example, the cell displayed in (d) has eleven facets, several of which are five-sided; analogous hcp examples exist. If cells that are identified as either fcc or hcp by W 0,4 1 are excluded from the q 6 distribution, these peaks vanish; the residual smooth distribution represents the non-crystalline background. Thus, for reliable detection of crystallinity, more information is required than contained in q 6 alone. This can be achieved by using multiple q l metrics and their distributional properties [21], or more specialized bond-orientational order metrics such as θ fcc , θ hcp [17]. Minkowski tensors, such as used in the present study, represent a more general approach; it is not necessary to choose a set of neighbors or bonds associated with a particle in order to evaluate the Minkowski tensors, and the Minkowski tensors are continuous functions of the particle coordinates. At the same time, they contain the information necessary to discriminate in a robust and specific way between disordered structure and different types of crystallinity. It can be shown that a relation exists linking the rotational invariants of higher-rank Minkowski tensors to variants of the bond-orientational order metrics q l , w l which have been amended by weighting factors proportional to the Voronoi facet areas [20].
In conclusion, we have demonstrated that local crystallinity in Lubachevsky-Stillinger sphere packings sets in when the packing is compactified beyond a critical packing fraction φ c ≈ 0.649. The packing fraction φ c marks the density below which LS configurations show no de-tectable degree of local crystallinity. Compactified above this limit, the system responds by the formation of local crystallinity.
The value φ c ≈ 0.649 is higher than experimental estimates for the RCP limit [1,2] and than the prediction based on mechanical contact numbers [11], but also than the MRJ packing fraction [15]. However, φ c is close to the packing fraction of ≈ 0.646 where polytetrahedral aggregates are most prevalent (Fig. 3 in [13]); the crystalline order metrics thus identify the conversion of polytetrahedral aggregates into crystalline structure, detected indirectly by Refs. [13,14]. The small but significant discrepancy between the critical packing fractions ≤ 0.64 of Ref. [1,2,8,11] on the one hand and ≈ 0.65 of Refs. [13,14] and of the present work on the other raise the caution that the mechanisms of dynamic arrest and isostaticity may be distinct from the alleged geometric order-disorder transition.
Given the nonequilibrium nature of jammed matter, one might be tempted to attribute the observed difference in packing fraction to details of the preparation protocol. The critical packing fraction of ≈ 0.65 is, however, not specific to the LS algorithm. For example, for packings generated using the force-balancing 'split algorithm', the geometric (rather than mechanical) contact number exhibits an anomaly close to φ ≈ 0.65 (Fig. 14 of Ref. [27]). Data by Bargiel and Tory for Jodrey-Tory packings can be successfully fitted with a critical packing fraction of ≈ 0.6495 [17]. Recent results by Klumov et al. for Jodrey-Tory and Lubachevsky-Stillinger packings, in terms of quantiles of the w 6 distribution, fix the geometric transition around φ ≈ 0.65 [21], but do not exclude crystallinity below φ c .
Future work needs to focus on the precise nature of the geometric transition occurring at φ c . Is the firstorder phase transition scenario viable, and if so, what are the coexistence densities? What is the signature of the transition in the ∆ fcc distribution? How does the local structure (quantified by Minkowski tensors) relate to the observed "Kauzmann" density ( Fig. 8 of Ref. [14])?