Jammed packings of 3D superellipsoids with tunable packing fraction, contact number, and ordering
Ye Yuan Kyle VanderWerf Mark D. Shattuck Corey S. O'Hern


fancy
plain
plain


[C]
    < g r a p h i c s >

[L]

    < g r a p h i c s >

[R]

    < g r a p h i c s >













 


iblabel[1]#1            
akefntext[1]
[0pt][r]thefnmark #1




1.125 



*

§ 
0pt4pt4pt
*

 §.§ 
0pt15pt1pt




[LO,RE]

    < g r a p h i c s >

[CO]

    < g r a p h i c s >

[CE]

    < g r a p h i c s >

[RO]1–LastPage  
[LE]  1–LastPage

 







 
 










[
  \begin@twocolumnfalse





    < g r a p h i c s >
     Jammed packings of 3D superellipsoids with tunable packing fraction, contact number, and ordering^† 


     
 


     Ye Yuan,^a,b Kyle VanderWerf,^c Mark D. Shattuck,^d and Corey S. O'Hern^∗^b,c,e,f 



    < g r a p h i c s >
     We carry out 
numerical studies of static packings of frictionless superellipsoidal 
particles in three spatial dimensions. We consider more than 200 different 
particle shapes by varying the three shape 
parameters that define superellipsoids. We characterize the structural and 
mechanical properties of both disordered and ordered packings using two 
packing-generation protocols. We perform athermal quasi-static compression 
simulations starting from either random, dilute configurations (Protocol 1) 
or thermalized, dense configurations (protocol 2), which allows us 
to tune the orientational order of the packings. In general, we find 
that the contact numbers at jamming onset for superellipsoid packings are hypostatic, 
with z_J < z_ iso, where z_ iso = 2d_f and d_f = 5 
or 6 depending on whether the particles are axi-symmetric or not. 
Over the full range of orientational order, we find that 
the number of quartic modes of the dynamical matrix for the packings always 
matches the number of missing contacts relative to the isostatic value. 
This result suggests that there are no 
mechanically redundant contacts for ordered, yet hypostatic packings 
of superellipsoidal particles. Additionally, we find that the packing 
fraction at jamming onset for diordered packings of superellipsoidal 
particles can 
be collapsed using two particle shape parameters, e.g. the asphericity 
𝒜 and 
reduced aspect ratio β of the particles. 



  



 \end@twocolumnfalse 


  ]





§ 







^a Department of Mechanics and Engineering Science, College of Engineering, Peking University, Beijing 100871, China; E-mail: yuanyepeking@pku.edu.cn
^b Department of Mechanical Engineering and Materials Science, Yale University, New Haven, Connecticut 06520, USA. 
^c Department of Physics, Yale University, New Haven, Connecticut 06520, USA. 
^d Benjamin Levich Institute and Physics Department, The City College of New York, New York, New York 10031, USA. 
^e Department of Applied Physics, Yale University, New Haven, Connecticut 06520, USA. 
^f Graduate Program in Computational Biology and Bioinformatics, Yale University, New Haven, Connecticut 06520, USA. 














§ INTRODUCTION


Athermal particulate materials, such as granular media, foams, and
emulsion droplets, typically jam, or become solid-like with a non-zero
static shear modulus when they are compressed to sufficiently large
packing fractions. <cit.>  Unwanted jamming occurs in many industrial
processes, such as clogging in hopper flows, <cit.> and controlled jamming
and unjamming has been used in robotics to grip soft, sharp, or fragile
objects. <cit.> Further, unjamming in geological systems, such as landslides 
and earthquakes, causes significant financial and human loss.  

Many prior studies have focused on jamming in model systems composed
of frictionless, spherical particles.  Disordered packings of
frictionless, monodisperse spherical particles are isostatic at
jamming onset with z_J = z_ iso contacts per particle, where
z_ iso = 2 d_f=6 and d_f = 3 is the number of translational
degrees of freedom for spheres, and with packing fraction at jamming onset
ϕ_J ≈ 0.64. <cit.> Previous work has characterized the critical
scaling of the structural and mechanical properties <cit.> and the anomalous
vibrational density of states <cit.> of jammed packings of spherical
particles. <cit.>

However, most athermal, particulate systems in industrial processes
and in nature are composed of highly non-spherical
particles. <cit.> In general,
disordered jammed packings of non-spherical particles are hypostatic
with z_J < z_ iso, where z_ iso=10 or 12 for
axisymmetric and non-axisymmetric particles,
respectively. <cit.> Thus, disordered jammed
packings can possess a range of contact numbers, 6 ≤ z_J ≤ 12,
and packing fractions at jamming onset that depend on the shape of the
constituent particles. In two spatial dimensions (2D), we showed
recently that disordered packings generated via athermal, quasistatic
compression for a wide variety of non-spherical shapes are
mechanically stable, despite the fact that z_J < z_
  iso. <cit.> We found that certain types
of contacts between nonsperical particles can constrain mulitple
degrees of freedom and that the number of missing contacts below the
isostatic value matches the number of quartic eigenmodes of the
dynamical matrix.  At jamming onset, perturbing the system along a
quartic eigenmode causes the total potential energy to increase as the
fourth power (not quadratically) in the perturbation
amplitude. <cit.>

Given that jammed packings of non-spherical particles can occur over a
wide range of contact numbers and packing fractions, is it possible to
a priori determine whether a system is jammed if we are only
given its z and ϕ?  For disordered packings of
monodisperse spheres, we know that if ϕ > 0.64 and z > 6, the
packing is jammed. For disordered packings of convex-shaped particles
in 2D, we found that the packing fraction at jamming onset can be
collapsed approximately onto a master curve that depends only on the
shape parameter A = p^2/4π a, where p is the perimeter and a
is the area of the particle. <cit.> In 2D, ϕ_J
≈ 0.84 for A=1, ϕ_J increases with A, reaching a peak
near A ≈ 1.1, and then decreases continuously with further
increases in A. Results for ϕ_J have also been reported for packings 
of nonspherical particles in 3D, but separately for each family of shapes, 
e.g., ellipsoids, <cit.>
spherocylinders, <cit.> and
spheropolyhedra. <cit.> Here, we will address the 
question of whether there is a  general relationship
between the packing fraction at jamming onset and one or more particle
shape parameters in packings of non-spherical particles in 3D.

Further, few studies have attempted to connect the contact number to
mechanical stability for ordered packings of non-spherical
particles, <cit.> despite the fact that packings of monodisperse particles
that deviate by less than 20% from perfect sphericity can possess
significant translational and orientational order. In particular, does
the relationship between the number of missing contacts below the
isostatic value and number of quartic modes hold for ordered or
partially ordered packings of non-spherical particles?  One might
expect that some of the “extra” contacts that occur in ordered
packings, may be mechanically redundant, <cit.> and
therefore will not contribute to the packing's stability, resulting in
a mismatch between the number of missing contacts and the number of
quartic modes.



We investigate these questions by generating static packings of monodisperse
frictionless, superellipsoidal-shaped particles in 3D using numerical
simulations. <cit.> We consider more than 200
different particle shapes by changing the shape parameters that define
superellipsoids. For each packing, we determine ϕ_J, z_J, the
orientational order, and the eigenvalues and eigenmodes of the
dynamical matrix. We carry out two packing-generation protocols. In
Protocol 1, we jam the packing via athermal quasistatic
compression,<cit.> starting from a random, dilute initial
configuration of particles. In Protocol 2, we thermalize an unjammed
configuration at an intermediate packing fraction before applying the
same athermal quasistatic compression protocol (Protocol 1). We find
that Protocol 1 generates globally disordered packings with a narrow
distribution of jammed packing fractions and contact numbers. Protocol
2, on the other hand, is able to generate packings of superellipsoidal
particles with a wide range of orientational order.

We find several key results. First, for disordered packings of
superellipsoidal particles in 3D generated via Protocol 1, we find
that the jammed packing fraction can be collapsed onto a master curve
using two shape parameters instead of only one as we found for
2D. <cit.>  In addition, we find that the
number of contacts, even in ordered packings of superellipsoids,
determines their mechanical stability. In particular, the number of
quartic eigenmodes of the dynamical matrix matches the number of
missing contacts relative to the isostatic value in ordered
superellipsoid packings, as well as in disordered packings.



The article is divided into several sections. In Sec. <ref>,
we review the definition of superellipsoids, describe the two
packing-generation protocols we implement, and define the
orientational order parameters we use to measure the degree of order
in jammed packings. In Sec. <ref>, we present our key
results. Finally, in Section 4, we summarize our results and discuss
directions for future research.  We also include three Appendices.  In
Appendix A, we show that we widely sample the two shape
parameters that characterize the shape of superellipsoids. In
Appendix B, we examine the local orientational order in 
superellipsoid packings. Finally, in Appendix C, we show the
correlation between the average curvature of the particles at
interparticle contacts and the average contact number for packings of
superellipsoids.




§ METHODS



In this section, we begin by defining the shape parameters for
superellipsoids, and explain the wide variation in particle shape that
is possible by tuning these parameters. Next, we describe our two
protocols, the athermal protocol 1, and the thermal protocol 2, which
we use to generate disordered and ordered jammed packings of these
shapes, respectively. We then discuss calculations of the eigenvalues
and eigenmodes of the dynamical matrix for superellipsoid packings to
measure their mechanical response. Finally, we define the two order
parameters that we use to quantify the orientational order in the
packings.



 §.§ Model of superellipsoidal particles


The surface of a superellipsodal particle located at the origin is defined by  

    |x/a|^2p + |y/b|^2p + |z/c|^2p = 1,

where a, b, and c (a ≤ b ≤ c) are the lengths of the
semi-major axes, and p is the deformation
parameter. <cit.> For superellipsoids,
there are three independent parameters that control the particle
shape, i.e., p and the two aspect ratios w_1=a/b and w_2=c/b. If a
= b, there is only one relevant asepct ratio w = c/a and if b =
c, w = a/c. Note that the particle shape reduces to a
superball when a = b = c. By tuning p, we can vary the
superellipsoid shape from ellipsoidal (p=1) to octahedral (p < 1)
and cuboidal (p > 1).  We focus our studies on five specific p-values: p =
0.75, 0.85, 1, 1.5, and 2.  


 

Instead of p and the aspect ratios, w_1 and w_2, the shape of
superellipsoids can also be characterized by p, the reduced aspect ratio
β=ac/b^2, and asphericity,

    𝒜 = 1 - (4π)^1/3(3V_p)^2/3/A_p,

where V_p and A_p give the particle volume and surface
area.<cit.> The shape parameter
β allows us to distinguish “flattened” (β < 1) versus
“elongated” (β > 1) shapes. The shape with β = 1 is termed
a self-dual ellipsoid, which shows anomalous properties in disordered<cit.> and dense<cit.>
packings. The asphericity satisfies 0 <  A < 1, and A=1 for
spheres. For the p values studied, the superellipsoidal particle
shape in the β-A plane is roughly bounded by the values
for prolate β_ max( A) and oblate ellipsoids
β_ min( A) as shown in Fig. <ref> in
Appendix A.  We focus on A-values from 0 to ∼
0.35 and sample β_ min( A) < β < β_
  max( A).


 



We consider pairwise, purely repulsive interactions between
superellipsoids using the Perram and Wertheim
formulation.<cit.>. For each pair of
overlapping superellipsoids i and j, we calculate the volume
scaling factor η_ij that brings the two superellipsoids to
exact tangency. The potential energy for particles i and j is then
defined by U_ij = ϵζ_ij^2/2, where ϵ is the
characteristic energy scale, ζ_ij = η_ij^2 - 1, and
η_ij≤ 1. The total potential energy is given by
U=∑_i>j U_ij. The repulsive force on particle i from j
f⃗_ij = ∇⃗_i U is given by

    f⃗_ij = 2ϵζ_ijη_ijn̂_ij/(r⃗_ij·n̂_ij)

where n̂_ij is unit normal of the tangent
plane between just-touching superellipsoids pointing toward i and
r⃗_ij is the center-to-center vector pointing from
superellipsoid j to i. The torque τ⃗_ij on particle i from
j is calculated using

    τ⃗_ij = l⃗_ij×f⃗_ij,

where l⃗_ij is the vector from the center of particle i to
the point of adjacency between superellipsoids i and j.  See
Fig. <ref> for an illustration of n̂_ij, r⃗_ij, and l⃗_ij for two contacting superellipsoids. We
will measure lengths, energies, and forces in terms of a,
ϵ, and ϵ/a.





 §.§ Packing-generation protocols


We generate jammed packings of N=400 frictionless, monodisperse
superellipsoidal particles in cubic simulation cells with periodic
boundary conditions using two compression protocols: 1) an athermal
protocol and 2) a thermal protocol.  For protocol 1, we first
initialize an overlap-free, dilute configuration of particles with
random positions and orientations. We then compress the configuration
in small increments of packing fraction, Δϕ=10^-3,
minimizing the total potential energy U using the L-BFGS
method<cit.> after each compression step. We terminate
the energy minimization procedure when the average normalized force on
a particle is below a small threshold, ⟨ |∑_j f⃗_ij|
⟩/⟨ f_ij⟩ < Δ, where Δ = 10^-4.
We stop compressing the system when the total potential energy per
particle first satisfies U/N > U_ tol, where U_ tol =
10^-10.  We then measure the packing fraction ϕ_J, contact
number z_J, and other quantities of the first jammed packing with
U/N > U_ tol that is closest to U_ tol.  We find that
the results presented here do not depend on the thresholds Δ
and U_ tol.  Examples of nine static packings of
superellipsoids with different shapes generated via protocol 1 are
shown in Fig. <ref>.

For protocol 2, we first thermalize unjammed configurations at
intermediate packing fractions ϕ_i ∼ 0.55, between the freezing
and melting packing fractions for hard superellipsoids,<cit.> using Monte
Carlo methods that do not allow particle overlaps for N_s steps. We then 
input these configurations into the compression and energy 
minimization procedure described in protocol 1. By varying 
N_s and ϕ_i, we can obtain jammed packings of superellipsoids with 
tunable ϕ_J, contact number z_J, and degree of orientational 
order. 

To calculate average quantities for the packing fraction, contact
number, and other quantities at jamming onset, we average over 5 to
10 independent initial conditions. We validated our methods for
generating jammed packings of superellipsoids by comparing our results
for ϕ_J from protocol 1 to those from recent studies of packings
of spheroids and superballs by Donev, et
  al.<cit.> (See
Fig. <ref>.)



 §.§ Dynamical matrix


The dynamical matrix, which provides all possible second derivatives
of the total potential energy with respect to the rotational and
translational degrees of freedom of the system, determines the linear
mechanical response of jammed particle packings. We define the
dynamical matrix as

    M_kl = ∂^2U/∂ξ_k∂ξ_l,

where ξ⃗ = {x_1, y_1, z_1, aθ_1, aϕ_1, aψ_1,…,
x_N, y_N, z_N, aθ_N, aϕ_N, aψ_N}, (x_i,y_i,z_i) is the
location of the center of particle i, and
(θ_i,ϕ_i,ψ_i) are the rotation angles about the x-,
y-, and z-axes used to define the orientation of particle i.
Thus, the dimension of the dynamical matrix is 6N× 6N.  For
jammed superellipsoid packings (in cubic simulation cells with
periodic boundary conditions), the dynamical matrix possesses 6N'-3
nonzero eigenvalues λ_i (with corresponding unit eigenvectors ê_i), where N'=N-N_r and N_r is the number of rattler
particles with unconstrained translational or rotational degrees of
freedom.

To determine M_kl, we calculated the first-order derivatives of
the dynamical matrix, ∂U/∂ξ_k analytically, and
calculated all of the the second-order derivatives numerically. We
find that the eigenvalues of the dynamical matrix do not depend
sensitively on the numerical derivatives for displacements <
10^-8.



 §.§ Order parameters


In packings of non-spherical particles, one can measure the degree of
order in the translational (i.e. positions of the particle centers)
and rotational (i.e. orientations of the particles) degrees of
freedom.  In the systems we study, when the particle orientations are
ordered, the particle positions also contain significant order.  Thus,
in these studies, we will focus on quantifying the orientational
order.
 
We measure the global nematic
S_2<cit.> and cubatic C_4 order
parameters.<cit.> S_2 is defined as the largest
eigenvalue of the 3× 3 matrix:

    S_αβ = 3/2⟨ŝ_α iŝ_β j⟩ - δ_αβ/2

where δ_αβ is the Kronecker delta, α,
β=x, y, and z, ŝ_α i is the α-compoent
of the unit vector that characterizes the orientation of particle i. 
and ⟨ . ⟩ indicates an average over all pairs of particles
i and j. 
ŝ_i is chosen as the shortest (longest) axis of the particle
when β < 1 (β > 1).  With this definition of ŝ_i, 
S_2 can capture stacking order that can occur in packings of flat shapes, 
as well as nematic order that can occur in packings of elongated shapes. 
S_2 = 0 for systems without orientational order and 1 for systems 
with complete particle alignment.  

The cubatic order parameter <cit.> C_4 is obtained 
by first calculating the fourth-order Legendre polynomial,  

    P_4(t̂,û_i) =1/8( 35 [t̂·û_i]^4 - 30[t̂·û_i]^2 + 3 ),

where t̂ is the unit vector aligned with one of the 3N
orientations of the semi-major axes of each of the particles and
û_i is a unit vector aligned with one of the three
orientations of the semi-major axes for particle i. For each
particle i in a given jammed packing, we select the û_i
that maximizes P_4(t̂, û_i) for a given t̂. We
then average P_4^ max(t̂) over all particles for a given
t̂ and define C_4 as the maximum over all 3N orientations
t̂. For C_4 ∼ 1, packings possess large cubatic order, which
can occur in packings of cube-like particles with p>1. In
Appendix B, we show results for the local nematic and cubatic order
in packings of superellipsoids.





§ RESULTS AND DISCUSSION



Our results are divided into two subsections. In
Sec. <ref>, we present our results for disordered packings
of superellipsoids generated via protocol 1. We show the global
nematic and cubatic order parameters for packings containing a wide
variety of superellipsoidal shapes. We find that the packing fraction
at jamming onset for disordered packings of superellipsoids can be
collapsed as a function of the two shape parameters, A and
β. In Sec. <ref>, we show that we can tune the packing
fraction and contact number at jamming onset by increasing the
orientational order of the packings generated via protocol 2. We
also show that, even for ordered packings, the number of quartic modes
of the dynamical matrix is equal to the isostatic number of contacts
minus the number of contacts in the packing.  Thus, we find
a direct link between the contact number and mechanical properties
even for ordered packings of superellipsoids.





 §.§ Disordered packings of superellipsoids



In this section, we focus on the structural propreties of
superellipsoid packings generated via protocol 1. In
Fig. <ref> (a), we show a scatter plot of the global nematic
S_2 and cubatic C_4 order parameters for all jammed packings
generated using protocol 1. We find that many of the packings are
disordered with S_2 and C_4 ∼ 1/√(N)∼ 0.07.  However,
as demonstrated in the inset to Fig. <ref> (a), S_2
increases as β decreases below 1 and the particle shape
flattens.  For more elongated shapes with β >1, S_2 is roughly
independent of β. We also find that the cubatic order increases
as the particles become more cube-shaped with p>1, even though the
packings were generated using the athermal protocol.  In
Fig. <ref> (b) and (c), we show example packings of flattened
and cube-like superellipsoids generated via protocol 1 wih elevated
values of S_2 and C_4. In (b), we show the local nematic order of
the particles for a packing of flattened superellipsoids with p=0.75
and w=0.3. In (c), we show the local cubatic order of the particles
for a packing of superballs with p=2. These packings possess local 
nematic and cubatic order. (See Appendix B.)

In Fig. <ref>, we show the packing fraction at jamming onset
ϕ_J as a function of the asphericity A for a variety of
superellipsoid shapes.  The relation between ϕ_J and
A is similar to that for packings of noncircular particles in
2D.<cit.> ϕ_J starts at a relatively low
value for spherical particles (i.e. random close packing for
monodisperse spheres with ϕ_J(0) ≈ 0.64), ϕ_J grows with
increasing asphericity, reaching a peak ϕ_J ∼ 0.70-0.74 near
A∼ 0.05, and then ϕ_J begins decreasing, falling
below ϕ_J(0) for A > 0.1. We also note that the data for
ϕ( A) does not collapse as well onto a single curve in 3D,
compared to the collapse of ϕ_J( A) for packings of 2D
noncircular particles.<cit.>

In Fig. <ref>, we show the contact number at jamming onset z_J
versus the asphericity A for (a) spheroids with an axis of
symmetry and z_ iso = 10 and for (b) all other particle shapes
with z_ iso = 12. z_J = 6 for isostatic packings of spherical
particles in the limit A→ 0. As found previously,
z_J for packings of nonspherical particles does not jump
discontinuously from 6 to z_ iso when A increases
above zero. Instead, z_J increases continuously with A.  z_J for some of the particle shapes reaches z_ iso for A < 0.35, e.g. oblate, prolate, self-dual, and general
elliposoids, but others, such as superellipsoids with p=0.75,
0.85, 1.5, and 2.0 do not. Note that z_ iso is smaller
for spheroids, compared to z_ iso for other non-axisymmetric
particle shapes, and thus the maximumum packing fraction for spheroids
is smaller than that for the other shapes we studied. We correlate 
values of z_J < z_ iso for superellipsoids with the curvature 
at interparticle contacts in Appendix C.

The packing fraction at jamming onset ϕ_J for packings of
superellipsoids does not completely collapse when plotted versus a
single shape parameter, e.g. the asphericity A.  (See
Fig. <ref>.) This result suggests that ϕ_J for packings of
nonspherical particles in 3D depends on two or more shape
parameters. In Fig. <ref> (a), we show ϕ_J versus the
reduced aspect ratio β for several values of the asphericity
A =0.1, 0.15, 0.2, and 0.25, excluding cube-like
superellipsoids with p>1. All of the curves ϕ_J(β) are
concave down for the different values of A. In
Fig. <ref> (b), we show a contour plot of ϕ_J as a
function of both β and A. We find that at small A, the largest ϕ_J, ϕ_J^ max, occurs near β =
1, however, ϕ_J^ max shifts to β >1 when A >
0.2. Thus, ϕ_J depends on both shape parameters A and
β.



 §.§ Tunable hypostaticity



In this section, we show that we can increase the nematic or cubatic
order in packings of superellipsoids using protocol 2 to generate
the packings. We compare the packing fraction and contact number at
jamming onset for packings generated via protocols 1 and 2. We
focus on packings of superballs with p = 1.25 and 1.5 and
packings of oblate ellipsoids with w = 0.3.

In Fig. <ref> (a), we show the global nematic S_2 and
cubatic C_4 order parameters versus the packing fraction at jamming
onset ϕ_J for single packings of oblate ellipsoids (with w=0.3)
and superballs (with p=1.25 and 1.5) generated via protocol
2. We also compare these results to those for packings of the same
shapes, but generated using protocol 1. Example packings are displayed
in Fig. <ref> (b).  We find that S_2 and C_4 < 0.1-0.2
for packings generated via protocol 1. However, S_2 and C_4 can
become larger than 0.7 for packings generated using protocol
2. For all shapes studied, ϕ_J increases with increasing
orientational order.



In Fig. <ref> (a), we show the eigenvalue spectrum of the
dynamical matrix (Eq.<ref>) sorted from smallest to largest for
packings of 6 different types of superellipsoids.  As found in
previous studies of packings of ellipsoids, the eigenvalue spectrum
has three distinct regimes. <cit.> For nearly spherical shapes, in regimes
2 and 3, the eigenmodes are purely rotational and translational,
respectively. Regimes 2 and 3 merge for systems with sufficiently
large asphericity 𝒜. “Quartic” modes occur in regime
1. When the system is perturbed along an eigenmode in this regime,
the change in the total potential energy Δ U between the unperturbed 
and pertrubed packings first increases quadratically
with the perturbation amplitude δ, but then scales as
δ^4 beyond a characteristic amplitude δ^* that scales to
zero with decreasing pressure. (See Fig. <ref> (a).) We found
in previous studies of nonspherical particles that the number of
quartic modes N_q matches the deviation in the number of contacts at
jamming onset from the isostatic value, i.e. N_c = N_c^ iso
-N_q, where N_c^ iso = d_f (N-N_r) -2.  We show this result
for packings of superellipsoids generated via protocol 1 in
Fig. <ref> (b).  This result shows that even though N_c <
N_c^ iso, disordered packings of superellipsoids generated via
protocol 1 are mechanically stable.



Is the relationship between the number of contacts and number of
quartic modes the same for packings of nonspherical particles with
significant orientational order?  For example, in ordered systems, it
is possible that some of the N_c contacts are redundant and
therefore do not provide independent constraints to block the degrees of
freedom in the packings.  In Fig. <ref> (a), we show the
contact number for packings of three types of superellipsoids
generated via protocol 2 that possess significant global nematic and
cubatic order (c.f. Fig. <ref> (a).)  The contact number
in these systems (z_J → 10) is much larger than that for
packings generated using protocol 1.

In Fig. <ref> (b), we show the eigenvalue spectrum of the
dynamical matrix for three packings of superballs with p = 1.5
generated using protocol 2. As shown previously, the spectrum
includes three regimes with a regime of quartic eigenmodes at the
lowest values.  Further, the crossover in behavior from Δ U ∼δ^2 to ∼δ^4 occurs at a similar value of δ^*
that scales to zero with decreasing pressure. (See Fig. <ref>
(b).) In the inset of Fig. <ref> (b), we show the number of
contacts N_c versus the number of quartic modes N_q for all of the
packings generated using protocol 2.  We find that even with
significant orientational order, the number of quartic modes matches
the deviation in the number of contacts from the isostatic value.
Thus, we find that there are no redundant contacts for hypostatic
packings of superellipsoids with z_J < z_ iso, and z_J determines 
their mechanical stability. 



 



§ CONCLUSIONS AND FUTURE DIRECTIONS



In this article, we carried out computational studies of jammed
packings of frictionless superellipsoids for more than 200 different
particles shapes in three spatial dimensions. We implemented two
protocols to generate static packings: protocol 1, which uses
athermal quasistatic compression, and protocol 2, which includes
thermal fluctuations and compression. Protocol 1 typically generates
packings with small values of the global nematic and cubatic
orientational order parameters, lower packing fraction ϕ_j and
contact number z_J at jamming onset.  In contrast, protocol 2
allows us to tune the orientational order (as well as ϕ_J and
z_J) in packings of superellipsoids over a much wider
range compared to those in protocol 1.

We found several important results. Prior studies of disordered jammed
packings of 2D nonspherical particles have found that the packing
fraction at jamming onset ϕ_J for a wide variety of shapes can be
collapsed onto a masterlike curve with respect to a single shape
parameter—the asphericity. <cit.> For disordered packings of
superellipsoids in 3D, we find that two shape parameters, e.g. the
asphericity A and reduced aspect ratio β, are required
to collapse ϕ_J. Additionally, prior studies have found that
packings of nonspherical particles are hypostatic with z_J < z_
  iso, and the number of missing contacts below the isostatic value
matches the number of quartic eigenmodes of the dynamical matrix. <cit.> Most
of these prior studies have considered disordered packings of
nonspherical particles with small values for global measures of
orientational order.  We find that for packings of superellipsoids
with significant orientational order generated via Protocol 2, the
number of missing contacts matches the number of quartic
eigenmodes. Thus, ordered packings of superellipsoids do not possess
any geometrically redundant contacts, and thus the contact number 
z_J < z_ iso directly determines their mechanical stability.

Our work opens up several new avenues of future research.  First, in
this work, we were able to generate packings of superellipsoids with
tunable orientational order, ϕ_J, and z_J. However, we only
considered packings with z_J < z_ iso. It will be interesting
to generate packings of nonspherical particles with even more order,
where z_J > z_ iso. In this case, do quartic modes still occur
and if so, what determines their number?  Another future research
direction involves packings of frictional non-spherical
particles. <cit.> Packings of frictional spherical particles can occur with
contact numbers that satisfy d_f +1 < z_J < 2d_f, where, d_f = 3
for spherical particles. <cit.> Prior studies have shown that packings of frictional
nonspherical particles can possess z_J < d_f + 1, <cit.> where for example d_f
= 5 for axisymmetric particles. Do these packings possess quartic modes, 
and if so, how many?  It is clear that much more work is needed to understand the number of contacts that are required to determine the mechanical stability 
of packings of frictional, nonspherical particles. 




§ CONFLICTS OF INTEREST

There are no conflicts to declare.



§ APPENDIX A: VARIATION OF THE SHAPE PARAMETERS Β AND A



In this Appendix, we show the range of reduced aspect ratio β
and asphericity A that can be achieved for superellipsoidal
particle shapes. For oblate and prolate ellipsoids, the asphericity
A(β) can be written explcity. For oblate ellipsoids, we
find

    A(β) = 2β^2/3/1+β^2/sinγln(1+sinγ/cosγ),

where γ = cos^-1β. For prolate ellipsoids, A(β) 
can be expressed as 

    A(β) = 2β^2/3/ 1+βα/sinα,

where α = cos^-1β^-1. In Fig. <ref>, we plot the 
relations between β and A for oblate and prolate ellipsoids, 
as well as β( A) for superellipsoids with p=2.  We find 
that these curves serve as upper and lower bounds for the shape 
parameters of all other shapes that we study for A < 0.35.   



§ APPENDIX B: LOCAL NEMATIC AND CUBATIC ORDER PARAMETERS



In the main text, for example in Figs. <ref> (a)
and <ref> (a), we showed results for the global nematic
S_2 and cubatic C_4 order parameters for packings of
superellipsoids. In these figures, we also show example packings from
the simulations with the particles colored according to the value of
the local nematic and cubatic order parameters. The local nematic
order parameter S_2^ loc is defined analogously to Eq. <ref>
as the largest eigenvalue of the 3 × 3 matrix:

    S^ loc_αβ = 3/2⟨ŝ_α i·ŝ_β j⟩_j
    - δ_αβ/2,

where ⟨ .⟩_j averages over particles j that overlap 
particle i. 

To define the local cubatic order parameter C_4^ loc for
particle i, we first calculate

    P_4(û_i, û_j) =1/8( 35 [û_j ·û_i]^4 - 30 [ û_j ·û_i]^2 + 3 ),

where û_j is a unit vector aligned with one of the three
orientations of the semi-major axes for particle j that overlaps
particle i. We first select the û_i orientation along one
of the three semi-major axes that maximizes P_4( u_i, 
  u_j) for a given u_j. We then average P_4^
  max(û_j) over all particles j that overlap i. The local
cubatic order parameter C_4^ loc is defined as the maximum over
the three orientations for u_j. We plot the global versus the
local orientational order parameters in
Fig. <ref>. For the nematic and cubatic order, the
global and local order grow proprotionately.
  




§ APPENDIX C: GAUSSIAN CURVATURE AT CONTACT POINTS



In this Appendix, we show that for packings of superellipsoids with
small contact numbers at jamming onset, the Gaussian curvature
K_G at the points of contact are typically small,
suggesting that two flat contacting surfaces can constrain multiple
rotational degrees of freedom.  In Fig. <ref>, we show the
probability distribution P(K_G), where K_G = K_G
(abc)^2/3, for superellipsoid packings generated via protocol
1. (Note that each contact point contributes two K_G values.) We
find that P(K_G) for packings of cube-like superellipsoids,
e.g. with p = 2 and small z_J, possess a wide tail that extends to
small values of K_G.  For other particle shapes, such as oblate and
prolate ellipsoids, P(K_G) is much narrower and does not extend to
small values of K_G.





§ ACKNOWLEDGEMENTS


We acknowledge support from NSF Grants Nos. CBET-1605178 (K. V. and
C. O.) and CMMI-1463455 (M. S.) and the China Scholarship Council Grant
No.  201806010289 (Y. Y.). This work was also supported by the High
Performance Computing facilities operated by Yale's Center for
Research Computing. In addition, we acknowledge the High-performance Computing
Platform of Peking University. We also thank S. Li at Peking University
for helpful discussions.  


















rsc 
