Hypostatic jammed packings of frictionless nonspherical particles

I Introduction

There have been a significant number of computational studies aimed at elucidating the jamming transition in static packings of frictionless spherical particles O’Hern et al. (2003); Liu and Nagel (2010); van Hecke (2009). Key findings from these studies include: i) sphere packings at jamming onset at packing fraction $\phi_{J}$ are isostatic (where the number of contacts matches the number of degrees of freedom, as shown in Fig. 1), ii) the coordination number, shear modulus, and other structural and mechanical quantities display power-law scaling as a function of the system’s pressure $P$ as packings are compressed above jamming onset at $P=0$, and iii) the density of vibrational modes develops a plateau at low frequencies $\omega$ that extends toward $\omega\rightarrow 0$ as the system approaches jamming onset. Many of these results are robust with respect to changes in the particle size polydispersity and different forms for the purely repulsive interparticle potential.

Most studies of jamming to date have been performed on packings of disks in 2D or spheres in 3D. More recently, both computational and experimental studies have begun focusing on packings of nonspherical shapes, such as ellipsoids Donev et al. (2007); Mailman et al. (2009); Schreck et al. (2012); Zeravcic et al. (2009); Basavaraj et al. (2006); Schaller et al. (2015); Donev et al. (2004); Man et al. (2005), spherocylinders Zhao et al. (2012); Blouwolff and Fraden (2006); Williams and Philipse (2003); Meng et al. (2016); Wouterse et al. (2007, 2009), polyhedra Jiao and Torquato (2011); Chen et al. (2014), and composite particles Schreck et al. (2010); Gaines et al. (2017); Miskin and Jaeger (2013); Baule et al. (2013). In particular, there is a well-established set of results on packings of frictionless ellipses (or ellipsoids in 3D). In general, static packings of frictionless ellipses are hypostatic with fewer contacts than the number of degrees of freedom using naïve contact counting. For amorphous mechanically stable (MS) packings of disks, the coordination number in the large-system limit is $z=2d_{f}$ (where $d_{f}=2$ is the number of degrees of freedom per particle) Tkachenko and Witten (1999). Thus, one might expect that the coordination number for ellipses in 2D would jump from $z=4$ to $z=6$ (with $d_{f}=3$) for any aspect ratio $\alpha>1$. However, $z(\alpha)$ increases continuously from $4$ at $\alpha=1$ and remains less than $6$ for all $\alpha$. We have shown that the number of missing contacts exactly matches the number of “quartic” eigenmodes from the dynamical matrix for which the potential energy increases as the fourth power of the displacement for perturbations along the corresponding eigenmode Schreck et al. (2012). In addition, the packing fraction at jamming onset $\phi_{J}(\alpha)$ possesses a peak near $\alpha\approx 1.5$, and then decreases for increasing $\alpha$.

Are these results for ellipses similar to those for all other nonspherical or elongated particle shapes? Prior results for packings of spherocylinders have shown that they are hypostatic Wouterse et al. (2007). However, packings of composite particles formed from collections of disks (2D) Schreck et al. (2010); Papanikolaou et al. (2013) or spheres (3D) Gaines et al. (2017) are isostatic at jamming onset. Unfortunately, very few studies explicitly check whether hypostatic packings are mechanically stable. The goal of this article is to determine which particle shapes can form mechanically stable (i.e. jammed), hypostatic packings, identify a key shape parameter that controls the forms of the coordination number $z$ and packing fraction $\phi_{J}$, and gain a fundamental understanding for why hypostatic packings are mechanically stable, even at jamming onset $P=0$.

To address these questions, we generate static packings using a compression and decompression scheme coupled with energy minimization for nine different convex particle shapes (ellipses, circulo-lines, circulo-triangles, circulo-pentagons, circulo-octogons, circulo-decagons Wang et al. (2015), dimers Schreck et al. (2010), dumbbells Han and Kim (2012), and Reuleaux triangles Atkinson et al. (2012)) in 2D. To study such a wide range of particle shapes, we developed a fully continuous and differentiable interparticle potential for circulo-lines and circulo-polygons, which allows us to generate packings with extremely accurate force and torque balance at very low pressure. We find several important results. First, we show that the jammed packing fraction for the nine particle shapes collapses onto a master-like curve when plotted versus the asphericity parameter ${\cal A}=p^{2}/4\pi a$, where $p$ is the perimeter and $a$ is the area of the particle. We also show that the coordination number $z({\cal A})$ follows a master-like curve when contacts between nearly parallel circulo-lines or nearly parallel sides of circulo-polygons are treated properly. In addition, for packings of circulo-lines, we calculate the fourth derivatives of the total potential energy along the quartic modes of the dynamical matrix Schreck et al. (2012) and show that the fourth derivatives are nonzero as $P\rightarrow 0$, which proves that these hypostatic packings are mechanically stable. Finally, we calculate the principal curvatures of the constraint surfaces in configuration space defined by each contact to identify which types of contacts in packings of circulo-lines allow them to be mechanically stable, while hypostatic.

This article is organized as follows. In Sec. II, we describe the compression and decompression plus minimization method we use to generate static packings of convex, nonspherical particles. In Sec. III.1, we present examples of static packings of several different convex, nonspherical particle shapes to put forward a conjecture concerning which nonspherical particle shapes form hypostatic packings and which always form isostatic packings. In Sec. III.2, we show that the packing fraction at jamming onset $\phi({\cal A})$ and coordination number $z({\cal A})$ display master-like forms when plotted versus the particle asphericity ${\cal A}$, for nine different nonspherical particle shapes. In this section, we also show results for the calculations of the fourth derivatives of the total potential of the static packings in the directions of dynamical matrix eigenmodes. Finally, in Sec. III.3, we calculate and analyze the principal curvatures of the constraint surfaces given by the interparticle contacts to understand the grain-scale mechanisms that enable hypostatic packings to be mechanically stable. We also include three Appendices. Appendix A describes the development of a continuous interparticle repulsive potential between circulo-lines and circulo-polygons, which allows us to generate extremely accurate force- and torque-balanced jammed packings near zero pressure. Appendix B describes how we generate different circulo-polygon shapes at constant asphericity ${\cal A}$. Finally, in Appendix C, we provide expressions for the elements of the dynamical matrix for packings of circulo-polygons.

II Methods

Using computer simulations, we generate static packings of frictionless, nonspherical, convex particles in 2D. The particles are nearly hard in the sense that we consider mechanically stable packings in the zero-pressure limit. We study nine different particle shapes: circulo-lines, circulo-triangles, circulo-pentagons, circulo-octagons, circulo-decagons, Reuleaux triangles, ellipses, dumbbells, and dimers. (See Table 1.) We focus on bidisperse mixtures in which half of the particles are large and half are small to prevent crystallization Speedy (1998); O’Hern et al. (2003). The large particles have areas that satisfy $a_{L}=1.4^{2}a_{S}$, where $a_{L,S}$ is the area of the large and small particles, respectively. Both particles have the same mass, $m$. We generated static packings at fixed asphericity ${\cal A}$ for the large and small particles over a wide range of ${\cal A}$. We employ periodic boundary conditions in square domains with edge length $L=1$ and system sizes that vary from $N=24$ to $480$ particles. Note that the term “convex particle shapes” stands for “shapes whose accessible contact surface is nowhere locally concave.” Our studies include dimers (which possess two points on the surface that are concave), circulo-lines (which contain regions of zero curvature), ellipses, and other explicitly convex particles.

We assume that particles $i$ and $j$ interact via the purely repulsive, pairwise linear spring potential,

where $k$ is the spring constant of the interaction and $\Theta(x)$ is the Heaviside step function. Below, lengths, energies, and pressures will be expressed in units of $L$, $kL^{2}$, and $k$, respectively. For disks, $r_{ij}$ is the separation between the centers of disks $i$ and $j$, and $\sigma_{ij}=R_{i}+R_{j}$ is the sum of the radii of disks $i$ and $j$.

For dimers, i.e. composite particles formed from two circular monomers, $r_{ij}$ is the center-to-center separation between each pair of interacting monomers, and $\sigma_{ij}$ is the sum of the radii of those monomers. A Reuleaux triangle is a shape that is constructed by joining three circular arcs of equal radius such that their intersection points (vertices of the Reuleaux triangle) are the centers of each circle. For this shape, we first identify whether two arcs are overlapping or whether a vertex is overlapping an arc. We then set $r_{ij}$ in Eq. 1 to be the distance between the center points of the overlapping arcs (in the case of two overlapping arcs) or the distance between the center point of the arc and the vertex (in the case of a vertex overlapping an arc). We set $\sigma_{ij}$ to be the sum of the radii of the two overlapping arcs, or the radius of the single arc when a vertex is overlapping an arc. (See Fig. 2.)

For ellipses, we take $r_{ij}$ to be the distance between the centers of the particles, and $\sigma_{ij}$ to be the center-center distance that would bring the particles exactly into contact at their current orientations Schreck et al. (2012). For dumbbell-shaped particles, we have multiple possible cases for $r_{ij}$, depending on their orientations. We calculate $r_{ij}$ either as the distance between each pair of circular ends, or between each circular end and the other particle’s shaft, with $\sigma_{ij}$ chosen to be the sum of the relevant radii in each case. The repulsive contact interactions between circulo-lines and -polygons are calculated in a similar fashion to dumbbells. However, because the regions of changing curvature in the case of circulo-lines and -polygons are accessible, unlike in the dumbbell case, additional constraints in the potential are necessary to prevent discontinuities in the pairwise torques and forces. For a thorough explanation of how we define a continuous, repulsive linear spring potential between circulo-lines and polygons, see Appendix A.

To generate static packings, we successively compress and decompress the system with each compression or decompression step followed by the conjugate gradient method to minimize the total potential energy $U=\sum_{i>j}U(r_{ij})$. We use a binary search algorithm to push the system to a target pressure $P=P_{0}$. If $P>P_{0}$, the system is decompressed isotropically, and if $P<P_{0}$, the system is compressed isotropically. Subsequently, we perform minimization of the enthalpy Smith et al. (2014) $H=U+P_{0}A$, where $A$ is the area of the system, the pressure $P=-dU/dA$, and $P_{0}=10^{-9}$ is the target pressure, with the particle positions and the box edge length as the degrees of freedom. Using this algorithm, we achieve accurate force and torque balance such that the squared forces $f^{2}_{i}$ and torques $\tau_{i}^{2}$ on a given particle $i$ do not exceed $10^{-25}$.

After generating each static packing, we calculate its dynamical matrix $M$, which is the Hessian matrix of second derivatives of the total potential energy $U$ with respect to the particle coordinates:

where $\xi_{i}=x_{i}$, $y_{i}$, and $\theta_{i}$, $x_{i}$ and $y_{i}$ are the coordinates of the geometric center of particle $i$, and $\theta_{i}$ characterizes the rotation angle of particle $i$. We then calculate the $3N$ eigenvalues $\lambda_{i}$ of $M$, and the corresponding eigenvectors ${\vec{\lambda}}_{i}$ with ${\vec{\lambda}}^{2}_{i}=1$. For more details on the calculation of the dynamical matrix elements, see Appendix C.

III Results

Our results are organized into three subsections. In Sec. III.1, we discuss which nonspherical particle shapes give rise to hypostatic packings, and then propose specific criteria that nonspherical particle shapes must satisfy to yield hypostatic packings. In Sec. III.2, we show the variation of the packing fraction $\phi$ and coordination number $z$ at jamming onset with particle asphericity ${\cal A}$ for packings of circulo-lines, circulo-polygons, and ellipses. Finally, in Sec. III.3, we calculate the principal curvatures of the inequality constraints in configuration space arising from interparticle contacts for hypostatic packings of circulo-lines to identify the specific types of contacts that allow static packings to be hypostatic, yet mechanically stable.

In this section, we discuss results for the contact number of static packings containing a variety of nonspherical particle shapes. Based on these results, we propose that frictionless convex particles will form hypostatic packings if both of the following two criteria are satisfied: (i) the particle has one or more nontrivial rotational degrees of freedom, and (ii) the particle cannot be defined as a union of a finite number of disks without changing its accessible contact surface. Below, we show several examples of systems that satisfy and do not satisfy these criteria.

First, disks do not satisfy (i) or (ii), and hence our conjecture predicts that disks will form isostatic, not hypostatic, packings. Next, we consider packings of circulo-lines that are prevented from rotating, and thus the particle’s orientation remains the same over the course of the packing simulations. (See Fig. 3 (a).) These particles obey criterion (ii), as a circulo-line can only be expressed as an infinite union of disks, but fail to meet criterion (i). Hence, the above conjecture predicts that these particles will form isostatic, not hypostatic packings. We also generated packings of bidisperse asymmetric dimers (Fig. 3 (b)). These particles meet criterion (i), since we allow them to rotate, but fail to meet criterion (ii), since dimers are made up of a union of two disks. Thus, our conjecture predicts that these particles will form isostatic, not hypostatic packings, as shown in Fig. 3 (b).

Finally, we generated packings of rotating circulo-lines, as well as Reuleaux triangles, examples of which are pictured in Fig. 3 (c) and (d), respectively. Both particles meet criterion (i), since they are allowed to rotate. Circulo-lines meet criterion (ii) as stated earlier. Reuleaux triangles also meet criterion (ii). Despite being comprised of a finite number of circular arcs, it is impossible to define them as a finite number of complete disks. Therefore, since both particle shapes meet both criteria, our conjecture predicts that they will form hypostatic, not isostatic packings. Ellipses also meet criteria (i) and (ii) and form hypostatic packings Schreck et al. (2012).

The importance of specifying “accessible contact surface” in criterion (ii) can be demonstrated by the two packings of dumbbells in Fig. 4. In both cases, the particles are allowed to rotate, so criterion (i) is satisfied. The packing in (a) also satisfies criterion (ii) because the shaft is part of the accessible contact surface of the constituent particles, and the shaft cannot be defined as a finite union of disks. Thus, we expect hypostatic packings for the dumbbells in Fig. 4 (a). In contrast, in Fig. 4 (b), the shaft is not part of the accessible contact surface, because it is too short to allow the end disks of other particles to come into contact with it. Thus, the particles in Fig. 4 (b) do not satisfy criterion (ii), because the accessible contact surface is a union of two disks. We expect packings generated using the dumbbells in Fig. 4 (b) to be isostatic.

III.2 Packing Fraction, Coordination Number, and Eigenvalues of the Dynamical Matrix

In this section, we describe studies of the packing fraction and coordination number of packings of nonspherical particles at jamming onset as a function of the particle asphericity ${\cal A}$. We also calculate the eigenvalues of the dynamical matrix for packings of circulo-lines and circulo-polygons and show the eigenvalue spectrum as a function of decreasing pressure. We find that hypostatic packings possess a band of eigenvalues, i.e. the ‘quartic modes’, for which the energy increases as the fourth power in amplitude when we perturb the system along their eigendirections. These quartic modes are not observed in isostatic packings. We further show that the fourth derivative of the total potential energy in the direction of these quartic modes does not vanish at zero pressure, proving that packings possessing quartic modes are mechanically stable, despite being hypostatic.

In Fig. 5, we plot the average packing fraction at jamming onset versus ${\cal A}$ for all of the nonspherical particles we considered. The data for $\langle\phi\rangle$ nearly collapses onto a master curve, which tends to $\langle\phi\rangle\approx 0.84$ for small ${\cal A}-1$, as found for packings of bidisperse disks Xu et al. (2005), forms a peak near ${\cal A}-1\approx 10^{-1}$, and decreases strongly for ${\cal A}-1>10^{-1}$. This result suggests that the asphericity can serve as common descriptor of the structural and mechanical properties of packings of nonspherical particles, i.e. jammed packings with similar ${\cal A}$ will possess similar properties.

In Fig. 6, we plot the average coordination number

where $N_{c}$ is the number of contacts in the packing. The $+1$ in the factor of $N_{c}+1$ is included to account for the $-1$ in the expression for the number of contacts $N_{c}=N_{c}^{0}=3N-1$ in isostatic packings of nonspherical particles in 2D, where $N_{c}^{0}$ is the isostatic number of contacts. $N_{r}$ is the number of rattler particles that have unconstrained translational and rotational degrees of freedom. $N_{s}$ is the number of ‘slider’ particles with a single unconstrained translational degree of freedom. An example of a slider particle is the yellow particle in the packing of circulo-lines in Fig. 3 (c), which can translate along its long axis without energy cost. Defining the coordination number as in Eq. 3 ensures that an isostatic packing of circulo-lines, circulo-polygons, or other nonspherical particles will have $\langle z\rangle=6$. If $\langle z\rangle<6$, the packing is hypostatic.

In Fig. 6, we show the coordination number $\langle z\rangle$ in Eq. 3 versus ${\cal A}-1$ for packings of ellipses and circulo-lines for two ways of defining a contact between two nearly parallel circulo-lines. At low asphericities, where the particle shape approaches a disk, a nearly parallel contact is only able to apply a small torque to the two contacting particles, making it unlikely to constrain a rotational degree of freedom in addition to a translational degree of freedom. Thus, at low asphericities, nearly parallel contacts should only be counted as a single constraint. In Fig. 6, we show that $\langle z\rangle$ for ellipses and circulo-lines (counting nearly parallel contacts once) both approach $4$ in the limit ${\cal A}-1$ tends to zero.

In contrast, at large asphericities, nearly parallel contacts between two circulo-lines prevent the particles from rotating and translating (in a direction perpendicular to their shafts). Thus, for large ${\cal A}-1$, nearly parallel contacts should be counted as two constraints. In Fig. 6, we show that the coordination number $\langle z\rangle$ for packings of circulo-lines approaches $6$ in the large ${\cal A}-1$ limit when nearly parallel contacts are counted twice. These results suggest that we must interpolate between counting parallel contacts once at low asphericities, and counting them twice as the asphericity increases.

One way to resolve the question of whether to count a nearly parallel contact between nonspherical particles as one or two constraints is to calculate the dynamical matrix (all second derivatives of the total potential energy with respect to the particle coordinates) of the static packings, and examine the spectrum of the dynamical matrix eigenvalues, which in the harmonic approximation give the vibrational frequencies of the packing Tanguy et al. (2002). For details on the calculation of the entries of the dynamical matrix for circulo-lines and -polygons, see Appendix C.

In Fig. 7, we show the eigenvalue spectrum (sorted from smallest to largest) for static packings of circulo-lines over a wide range of aspect ratios ${\cal A}-1$ from $\approx 10^{-6}$ to $1$ (decreasing from top to bottom). As found in Ref. Mailman et al. (2009) for ellipse packings, the eigenvalue spectra for packings of circulo-lines possess several distinct regions. Region ($\lambda_{i}\lesssim 10^{-14}$, which is set by numerical precision) corresponds to unconstrained degrees of freedom, such as overall translations from periodic boundary conditions, and rattler and slider particles. Region $1$ ($10^{-14}\lesssim\lambda_{i}\lesssim 4\times 10^{-8}$) corresponds to “quartic modes,” whose number is determined by the number of missing contacts relative to the isostatic contact number. For the asphericities we consider, regions $2$ and $3$ correspond to eigenmodes with predominantly rotational and translational motion, respectively.

If we focus on all but the three smallest asphericities (i.e. the three rightmost curves in Fig. 7), we can define a cutoff value $\lambda_{c}$ that clearly separates regions $1$ and $2$. For packings of $N=100$ circulo-lines at pressure $P_{0}=10^{-9}$, $\lambda_{c}\approx 4\times 10^{-8}$. For asphericities where $\lambda_{c}$ distinguishes regions $1$ and $2$, the number of contacts in packings of circulo-lines satisfies $N_{c}=N^{0}_{c}-N_{1}$, where $N_{1}$ is the number of eigenvalues in region $1$. A key observation is that defining the number of contacts in this way for intermediate and high asphericities is the same as if $N_{c}$ is determined by the number of particle contacts, with nearly parallel contacts counted twice. For asphericities where the difference between regions $1$ and $2$ is more ambiguous, we still use $\lambda_{c}$ to determine whether a given eigenvalue belongs to region $1$ or $2$. For the lowest asphericities, we find that defining $N_{c}=N^{0}_{c}-N_{1}$ corresponds to counting one constraint for each nearly parallel contact.

In Fig. 8, we plot the average coordination number $\langle z\rangle$ from Eq. 3 using $N_{c}=N^{0}_{c}-N_{1}$ versus ${\cal A}-1$ for $N=100$ packings of circulo-lines and circulo-polygons. At low asphericities ${\cal A}-1$, the coordination number for packings of circulo-lines and circulo-polygons, as well as ellipses, approaches $\langle z\rangle=4$, which is expected for bidisperse disk packings. At large asphericities, $\langle z\rangle=6$ for packings of circulo-lines and circulo-polygons as expected for isostatic packings with $2$ translational and $1$ rotational degree of freedom per particle. $\langle z\rangle$ for ellipse packings plateaus for large ${\cal A}-1$. However, the current data suggests that the plateau value is less than $6$, indicating that ellipse packings are hypostatic for all ${\cal A}-1$.

An interesting feature in $\langle z\rangle({\cal A})$ for static packings of circulo-lines and -polygons is the plateau in $\langle z\rangle$ that occurs near ${\cal A}-1\approx 10^{-5}$ in Fig. 8. Our results suggest that the plateau is likely an artifact of the small, but nonzero pressure of the static packings. If the particles are overcompressed, even slightly, nearly parallel contacts will be able to exert larger torques than they would at zero pressure, which causes more eigenvalues to be above the eigenvalue threshold $\lambda_{c}$, and contacts to be counted as two constraints instead of one. Thus, as we decrease the pressure to zero, we expect to count fewer of these nearly parallel contacts as two constraints and the plateau in $\langle z\rangle$ near ${\cal A}-1\approx 10^{-5}$ will decrease. As ${\cal A}-1$ decreases below $10^{-5}$, the effects from overcompression are less important, and the nearly parallel contacts are only counted once.

In Fig. 9, we plot the eigenvalues $\lambda_{i}$ sorted from smallest to largest for static packings of five different particle shapes (circulo-lines, -triangles, -pentagons, -octagons, and decagons) at the same asphericity, ${\cal A}=1.1$. We find that the eigenvalue spectra for all of these shapes are nearly identical. This behavior differs markedly from that in Fig. 7, where we show the eigenvalue spectra for packings with the same particle shape (circulo-lines), but at different values of the asphericity. Circulo-polygons with $n$ sides possess $2n-3$ parameters that specify their shape (not counting uniform scaling of lengths). Our results suggest that asphericity is a key parameter in determining the structure, geometry, and physical properties of hypostatic packings.

The reason why the eigenvalues in region $1$ (c.f. Fig. 7) are referred to as “quartic modes” is that, for perturbations along the corresponding eigenvectors, the total potential energy scales quartically with the amplitude of the perturbation, rather than quadratically, as one would expect for mechanically stable packings Mailman et al. (2009); Schreck et al. (2012). In Fig. 10, we plot eigenvalues from regions $1$ and $2$ as a function of pressure $P_{0}$ for a static packing of $N=32$ circulo-lines at asphericity ${\cal A}=1.03$. The eigenvalues from region $2$ are independent of pressure, whereas the eigenvalues from region $1$ scale linearly with pressure. Thus, for packings of circulo-lines and other particle shapes that yield hypostatic packings, the eigenvalues corresponding to the quartic modes are zero at jamming onset ($P_{0}=0$). This result agrees with prior studies of hypostatic packings of ellipses and ellipsoids Schreck et al. (2012).

Perturbations along the quartic modes are constrained to fourth order. In Fig. 11, we show the fourth derivatives of the total potential energy $d^{4}U/d{{\vec{\lambda}}}^{4}$ in the directions of the nine eigenmodes in region $1$ (that are depicted near the bottom of Fig. 10). We find that the fourth derivatives along eigenmodes in region $1$ do not depend on pressure, and thus remain nonzero at zero pressure. These findings demonstrate that hypostatic packings are fully constrained at zero pressure—in some directions by quadratic potentials and in other directions by quartic potentials.

III.3 Convex versus concave constraints

Why are hypostatic packings of circulo-lines and other nonspherical particles mechanically stable when they possess fewer contacts than the isostatic number, $N_{c}<N_{c}^{0}$? We have already shown that the number of missing contacts $N_{c}^{0}-N_{c}$ matches the number of quartic modes along which the energy increases quartically, not quadratically, with the perturbation amplitude. In the other $N_{c}$ eigendirections of the dynamical matrix, the energy increases quadratically with the perturbation amplitude. As a result, there are no directions in configuration space for which these hypostatic packings can be perturbed without energy cost, and thus they are mechanically stable.

To more fully address the question of how hypostatic packings of nonspherical particles can be mechanically stable, we consider the so-called “feasible region” of configuration space near each static packing for packing fractions slightly below jamming onset Donev et al. (2007). The feasible region near a given static packing includes all configurations for which there are no particle overlaps. The boundaries of this region are determined by all of the interparticle contacts, each of which corresponds to an inequality among the particle coordinates specifying when pairs of particles do not overlap. Points in configuration space that satisfy all of the inequalities are inside the feasible region. For mechanically stable packings, as the packing fraction is increased, the feasible region shrinks and becomes bounded and compact, preventing particle rearrangements that would allow the system to transition to a different packing. A static packing is mechanically stable if the feasible region of accessible configurations shrinks to a single point at jamming onset.

The number of constraints required to bound the feasible region depends on the curvature of the inequality constraints in configuration space, i.e. whether the constraints are concave or convex Donev et al. (2007). The inequality constraints that arise in disk packings are always concave. In particular, in disk packings, the curvature of each constraint is equal to minus the reciprocal of the sum of the radii of the two disks in contact. As a result, the number of contacts required to bound the feasible region for a mechanically stable packing of $N$ disks is $2N+1$ (minus $2$ from overall translations in periodic boundary conditions). Thus, hypostatic packings of nonspherical particles must possess contacts that give rise to bounding surfaces with convex curvature, which allows packings to be mechanically stable with fewer than the isostatic number of contacts.

In Fig. 12, we show a simple configuration involving a circulo-line that gives rise to a convex constraint. We consider three points at fixed positions. These points represent less strict constraints than contacts with other circulo-lines, and thus, if these three points can constrain a circulo-line, three contacting circulo-lines will constrain an interior circulo-line as well. We initialize a circulo-line at several locations between the three points, and then increase the size of the interior circulo-line until it is constrained by the three points. After the circulo-line is constrained, we shrink its diameter by $10^{-7}$ so that it no longer overlaps the bounding points. The feasible region of the slightly undercompressed circulo-line is shown in Fig. 12 (a).

For an isostatic system, four contacts are required to constrain a circulo-line. However, we find configurations in which a circulo-line is constrained by only three contacts. Fig. 12 (a) illustrates the reason that only three contacts are necessary: one of the contacts (open circle on the top shaft) gives rise to a constraint with convex curvature in configuration space. In contrast, the other two contacts (filled circles), which are on the end caps of the circulo-line, give rise to constraints with concave curvature. This example suggests that only certain types of contacts between circulo-lines generate constraints with convex curvature, and thus the number of contacts required for mechanical stability is less than the isostatic number when these types of contacts are present.

To verify that the circulo-line “packing” in Fig. 12 (b) is mechanically stable, we numerically calculated the volume $V$ and surface area $S$ of the feasible region as a function of the degree of undercompression, $\Delta R=R_{j}-R$, where $R_{j}$ is the radius of the interior circulo-line at which the system is jammed. In Fig. 13, we show that both $V$ and $S$ display power-law scaling with $\Delta R$, emphasizing that the feasible region for hypostatic packings shrinks to a point, and thus these packings are mechanically stable.

To further investigate the effect of convex and concave constraints on a hypostatic jammed packing, we measured the curvatures of the inequality constraints for each contact in a static packing with $N=24$ bidisperse circulo-lines with asphericity ${\cal A}=1.04$. We classified the contacts into five types as defined in Appendix A. Parallel contacts can involve the shaft of one circulo-line (middle) and the end cap of another (end). This arrangement gives rise to two types of contacts, one for the circulo-line with a contact on its end and another for the circulo-line with a contact on its middle. Similarly, the shaft (middle) of one circulo-line can be in contact with the end cap (end) of another, but the long axes are not parallel. This arrangement again gives rise to two types of contacts, one for the circulo-line with a contact on its end and another for the circulo-line with a contact on its middle. In addition, the ends of two circulo-lines can be in contact.

The average curvatures of the bounding surfaces for each contact type in a static packing of $N=24$ bidisperse circulo-lines are compiled in Table 2. (We find similar average values for other $N=24$ packings of bidisperse circulo-lines.) From this data, we can draw several conclusions about the contribution of each type of contact to the stability of circulo-line packings. First, end-end contacts yield concave constraints in configuration space, and thus on their own do not give rise to mechanically stable hypostatic packings. In contrast, end-middle contacts have a positive principal curvature for the circulo-line whose middle is in contact, and thus serve to stabilize hypostatic packings. Parallel contacts also possess a positive curvature associated with the circulo-line whose middle is in contact. However, note that the concave curvature for circulo-lines whose end is in parallel contact is much smaller than the concave curvature of the end circulo-line for end-middle contacts. This means that for circulo-lines with end contacts, the parallel contacts are more “stabilizing” than the end-middle contacts, and therefore they are more frequent in mechanically stable hypostatic circulo-line packings than other end-middle contacts.

The above observations about the curvatures of the inequality constraints in configuration space can help explain the distribution of contact angles $P(\psi)$ in static packings of elongated particles Tian et al. (2015); Marschall and Teitel shown in Fig. 14. This figure shows that, even for packings of circulo-lines at very small asphericities, parallel contacts are highly probable, despite the fact that the range of angles for parallel contacts at low asphericities is small. This behavior for $P(\psi)$ can be explained by the fact that end-middle and parallel contacts can contribute to making a hypostatic packing mechanically stable, whereas end-end contacts cannot. (See Table 2.) Thus, end-middle and parallel contacts (whose contact angles are close to $90^{\circ}$ at low asphericities) must be present to stabilize hypostatic packings of low-asphericity circulo-lines. As shown in Fig. 14, $P(\psi)$ is similar for both ellipse and circulo-line packings.

IV Conclusions and Future Directions

In this article, we carried out computational studies of static packings of frictionless nonspherical particles in 2D. We developed an interparticle potential for circulo-lines and -polygons that generates continuous pair forces and torques as a function of the particle coordinates. As a result, we are able to compare the structural and mechanical properties of mechanically stable packings of nine different nonspherical particle shapes: circulo-lines, -triangles, -pentagons, -octagons, -decagons, asymmetric dimers, dumbbells, Reuleaux triangles, and ellipses. Our studies place a particular emphasis on the question of which particle shapes give rise to hypostatic mechanically stable packings with fewer contacts than the isostatic number.

We conjecture that to form hypostatic mechanically stable packings, frictionless, convex particles must satisfy the following two criteria: (i) the particle has one or more nontrivial rotational degrees of freedom, and (ii) the particle cannot be defined as a union of a finite number of complete disks without changing its accessible contact surface. If the particle does not satisfy both criteria, we expect it to form isostatic packings. Packings of the nine particle shapes we considered are consistent with this conjecture. Future research can investigate methods to analytically prove this conjecture Roux (2000).

We then studied the packing fraction $\phi$ and coordination number $z$ at jamming onset for packings of a number of different types of nonspherical shapes in 2D as a function of asphericity ${\cal A}$. To do this, we resolved the ambiguity in the constraint counting of nearly parallel contacts of circulo-lines and -polygons using the branched structure of the eigenvalue spectra of the dynamical matrix. In future research, we will study the coordination number of packings of sphero-cylinders and -polygons in 3D, and compare the results to those in 2D, since it is extremely unlikely for sphero-cylinders and -polygons to form nearly parallel contacts.

We find that the packing fraction and coordination number obey approximate master curves when plotted versus the asphericity. Further, the eigenvalue spectra for different particle shapes, at the same ${\cal A}$, collapse. These results suggest that asphericity is a key parameter in determining the structure, geometry, and mechanical properties of hypostatic packings. For $n$-sided circulo-polygons, there are $2n-3$ parameters that specify their shape. In future studies, we will investigate additional shape parameters, such as the ratios of the area moments and others Schröder-Turk et al. (2013), to better understand the coupling between the shape parameter space and the properties of hypostatic packings of nonspherical particles.

We also demonstrated that hypostatic packings of circulo-lines (and by analogy circulo-polygons) are mechanically stable by showing that even though the eigenvalues of the dynamical matrix for the quartic modes tend to zero at zero pressure, the fourth derivatives of the total potential energy in the directions of the quartic modes do not. Thus, hypostatic packings of nonspherical particles are stable to perturbations in all directions in configuration space. Perturbations in some directions give rise to quadratic potentials, whereas other directions give rise to quartic potentials. In the directions with quartic potentials, we expect large anharmonic contributions to the vibrational and mechanical response Schreck et al. (2014).

In addition, we measured the curvatures of the inequality constraints that arise from interparticle contacts in hypostatic packings of circulo-lines to better understand the grain-scale mechanisms that allow hypostatic packings to be mechanically stable. The contacts in isostatic disk packings give rise to inequality constraints with only concave (negative) curvatures. In contrast, hypostatic packings of circulo-lines (and other nonspherical particles) possess different types of contacts (e.g. end-end and end-middle). Some types yield inequality constraints with concave curvatures and others yield inequality constraints with convex curvatures. We find that contacts with convex inequality constraints are present even at small asphericities. The contacts with convex inequality constraints allow the feasible region of slightly undercompressed hypostatic packings to be compact, bounded, and shrink to zero in the limit that the free volume tends to zero.

Acknowledgments

The authors acknowledge financial support from NSF Grant Nos. CMMI-1462439 (C.O.), CMMI-1463455 (M.S.), and CBET-1605178 (C.O. and K.V.), NIH Training Grant, Grant No. 1T32EB019941 (K.V.), and the Raymond and Beverly Sackler Institute for Biological, Physical, and Engineering Sciences (C. O. and K. V.). We also acknowledge the China Scholarship Council that supported Weiwei Jin’s visit to Yale University. In addition, this work was supported by the High Performance Computing facilities operated by, and the staff of, the Yale Center for Research Computing. We thank T. Marschall and S. Teitel for helpful conversations.

Appendix A Continuous Potential between Circulo-lines and -Polygons

The repulsive potential between two circulo-lines is given by Eq. (1), where $r_{ij}$ is the magnitude of ${\vec{r}}_{ij}$, which points from the location where the force is applied on circulo-line $j$ to the location where the force is applied on circulo-line $i$. These points of contact can be located on the ends or the shaft (middle) of a circulo-line. In this Appendix, we define the overlap distance $\delta=\sigma_{ij}-r_{ij}$, which will depend on the type of contact that occurs between two circulo-lines.

There are three types of interparticle contacts that occur in packings of circulo-lines: 1) the end of one circulo-line is in contact with the middle of another (Fig. 15), 2) the shafts of two circulo-lines are in contact and the circulo-lines are nearly parallel (Figs. 16 and 17), and 3) the ends of two circulo-lines are in contact (Fig. 18). Below, we define the overlap distance $\delta$ in the circulo-line potential (Eq. 1) for each type of contact.

End-middle contacts occur when the endcap of one circulo-line makes contact with the middle of another circulo-line, but does not overlap with either of the other circulo-line’s endcaps. (See Fig. 15.) In this case, we assume that the separation vector ${\vec{r}}_{ij}$ between circulo-lines points from the end of the shaft of the circulo-line with the end contact to the shaft of the other circulo-line. ${\vec{r}}_{ij}$ is perpendicular to the shaft of the circulo-line with the middle contact. The overlap between circulo-lines with an end-middle contact is $\delta=\sigma_{ij}-r_{ij}$, as shown in Fig. 15.

A.1.2 Parallel and Nearly Parallel Contacts

For parallel and nearly parallel contacts, an endcap of both circulo-lines overlaps the shaft of the other circulo-line. In this case, the spring potential in Eq. 1 for both overlaps is calculated as for end-middle contacts. If the circulo-lines are parallel, as in Fig. 16, $\delta=\sigma_{ij}-r_{ij}$ is the same for both overlaps. However, for nearly parallel contacts, as in Fig. 17, the separations are different for the two end-middle overlaps. This method for treating end-middle, parallel, and nearly parallel contacts ensures continuity of the potential as a function of the particle coordinates. If the circulo-lines in Fig. 17 rotate until their orientations match Fig. 15, the potential, force, and torque must all change continuously. Using our method, $\delta_{2}$ decreases continuously to zero as the contact evolves from that in Fig. 17 to that in Fig. 15. In addition, $\delta_{1}$ decreases continuously to zero as the circulo-lines in Fig. 17 rotate until $\delta_{2}$ is the only overlap.

A.1.3 End-End Contact

Also, suppose that we slide the two circulo-lines in Fig. 16 away from each other until they are similar to the configuration in Fig. 18 and form an end-end contact. In this case, we assume that the two overlap potentials add together as soon as the two relevant ends of the circulo-line shafts slide past each other. Therefore, to ensure continuity, the interaction potential for an end-end contact must be twice as large as that for an end-middle contact. Hence, we use $U=k\delta^{2}$ for end-end contacts.

However, this treatment of end-end contacts creates a discontinuity for the configuration in Fig. 19. If we imagine sliding the circulo-lines past each other until the overlap $\delta_{1}$ is associated with an end-end contact, the potential will suffer a discontinuous jump from $\tfrac{1}{2}k\delta_{1}^{2}$ to $k\delta_{1}^{2}$ since the end-end contact potential is twice as large as an end-middle potential. To remedy this discontinuity, we add the end-end contact potential between the two relevant endpoints as soon as they become close enough to overlap. However, we do not make the end-end potential twice as large in this case. Hence, the potential in this case is given by $U=\tfrac{1}{2}k\left(\delta_{1}^{2}+\delta_{2}^{2}\right)$. Thus, when we perform that same sliding transformation, the potential will grow continuously from $\tfrac{1}{2}k\delta_{1}^{2}$ to $k\delta_{1}^{2}$ as $\delta_{2}$ grows continuously from 0 to $\delta_{1}$. Note that we do not add this end-end overlap potential if two end-middle contacts are present, as in Fig. 20, because in that case, the potential will already change continuously as described in the previous subsection, and hence there is no discontinuity to remedy.

A.2 Generalizing the Circulo-line Potential to Circulo-Polygons

Generalizing our continuous circulo-line potential to circulo-polygons, such as those pictured in Fig. 21, is straightforward. We simply calculate the potential between all pairs of circulo-lines that comprise each circulo-polygon. For example, in Fig. 21, since the vertex of the top circulo-triangle is shared by two circulo-lines, we count the end-middle contact twice, and hence the overlap potential is $U=k\delta^{2}$.

Appendix B Generation of Circulo-Polygons

A circulo-polygon is formed through a Minkowski sum of a polygon and a disk with a radius $r$ Oks and Sharir (2006), which is equivalent to the sweeping of the disk around the profile of the polygon as in Fig. 22 (a). The shape of a circulo-polygon with $n$ edges is fully specified by $2n-3$ independent parameters. In this work, we focus on the asphericity shape parameter ${\cal A}$, which measures the deviation of a given shape from a circle in 2D.

We study bidisperse packings of circulo-polygons with asphericity ${\cal A}$ for which half of the circulo-polygons are large and half are small. The large circulo-polygons have areas that satisfy $a_{L}=1.4^{2}a_{S}$, where $a_{L,S}$ is the area of the large and small circulo-polygons, respectively. The large circulo-polygons (and small ones) have different shapes at the same ${\cal A}$. We generate different circulo-polygons at the same ${\cal A}$ using the following two-step approach: 1) We first randomly select $n$ points on a unit disk as the vertices of an $n$-sided polygon. The radius $r$ of the circulo-polygon is set to be $n$ percent of the perimeter of the polygon. 2) If the asphericity of the current circulo-polygon is smaller than the target $\cal A$, a vertex $J$ is randomly chosen and then stretched or shortened along the direction between the vertex $J$ and the center $O$ of the unit disk, by a distance randomly chosen between 0 and the distance between $J$ and the intersection of $JO$ with the line segment connecting the two neighboring vertices, as shown in Fig. 22 (b). This deformation is accepted only if the asphericity of the new shape is closer to $\cal A$ than the original and the new shape is still convex. If the asphericity of the current circulo-polygon exceeds ${\cal A}$, the radius $r$ is increased to match the target $\cal A$. We repeat step $2$ until a circulo-polygon with ${\cal A}$ is obtained.

Appendix C Dynamical Matrix Elements of Circulo-Polygon Packings

In this Appendix, we provide explicit expressions for the dynamical matrix elements for static packings of circulo-polygons that interact via the purely repulsive linear spring potential in Eq. 1. In this expression, $R_{i}$ is radius that forms the edge of circulo-polygon $i$ and ${\vec{r}}_{ji}$ is the separation vector from circulo-polygon $i$ to $j$, which is is given by

where $\vec{c}_{ji}=(x_{j}-x_{i},y_{j}-y_{i})\equiv(p,q)$ is the center-to-center separation between circulo-polygons, {\cal R}_{i}=\left[\begin{array}[]{cccc}\cos\theta_{i}&-\sin\theta_{i}\\ \sin\theta_{i}&\cos\theta_{i}\end{array}\right] is the rotation matrix in 2D, $\theta_{i}$ is the orientation of particle $i$ relative to the $x$-axis, $\vec{v}_{m}$ is the vector from the center of particle $i$ to the center of its corresponding edge $m$ when particle $i$ at zero rotation, ${\hat{u}}_{m}$ is the unit vector along edge $m$ at zero rotation, and $l_{m}$ indicates the distance between the contact point and the center of edge $m$.

The dynamical matrix requires the calculation of the second derivatives of the total potential energy $U$, and can be expressed in terms of the first and second derivatives of the contact distance $r_{ji}$ with respect to the particle coordinates:

where $\xi_{i}=x_{i}$, $y_{i}$, or $\theta_{i}$,

There are two types of contacts among circulo-polygons. The first type is a vertex-to-edge contact. Assuming the end point of edge $n$ on particle $j$ is in contact with edge $m$ on particle $i$, $l_{n}$ is half of the length of edge $n$, and $l_{m}$ can be written as

In Eq. 8, $\vec{v}_{m}$, $\vec{v}_{n}$, $\hat{u}_{m}$, and $\hat{u}_{n}$ are defined as

where $m$ and $n$ are the magnitudes of $\vec{v}_{m}$ and $\vec{v}_{n}$, respectively, and the angles $a$, $b$, $d$, and $e$ are the orientations of $\vec{v}_{m}$, $\vec{v}_{n}$, $\hat{u}_{m}$, and $\hat{u}_{n}$, respectively. The contact distance $r_{ji}$ is

The nonzero first and second derivatives can be expressed as:

In the expressions in Eqs. 16-28, $\Xi$ is defined as:

All of the other first and second derivatives are zero.

The second type of contact between circulo-polygons is a a contact between two vertices. In this case, $l_{m}$ and $l_{n}$ are each half the lengths of edges $m$ and $n$, respectively. The $x$- and $y$-components of separation vector $\vec{r}_{ji}$ are

For vertex-vertex contacts, the nonzero first and second derivatives are:

The other derivatives, $\frac{\partial\vec{r}_{ji}}{\partial\xi_{i}}\cdot\frac{\partial\vec{r}_{ji}}{\partial\xi_{j}}$ and $\vec{r}_{ji}\cdot\frac{\partial^{2}\vec{r}_{ji}}{\partial\xi_{i}\partial\xi_{j}}$, that are not listed above are zero.

References