Constraints and vibrations in static packings of ellipsoidal particles

0.85 MB·English

by  Carl F. Schreck

Constraints and vibrations in static packings of ellipsoidal particles
Carl F. Schreck1, Mitch Mailman2, Bulbul Chakraborty3, and Corey S. O’Hern4,1
1Department of Physics, Yale University, New Haven, Connecticut 06520-8120, USA
2Department of Physics, University of Maryland, College Park, Maryland 20742, USA
3Martin Fisher School of Physics, Brandeis University,
Mail Stop 057, Waltham, MA 02454-9110, USA and
4Department of Mechanical Engineering & Materials Science,
Yale University, New Haven, Connecticut 06520-8260, USA
We numerically investigate themechanical properties of static packings of ellipsoidal particles in
2Dand3Doverarangeofaspectratioandcompression∆φ. Whileamorphouspackingsofspherical
2
1 particles at jamming onset (∆φ = 0) are isostatic and possess the minimum contact number ziso
0 required for them to be collectively jammed, amorphous packings of ellipsoidal particles generally
2 possess fewer contacts than expected for collective jamming (z < ziso) from naive counting argu-
ments, which assume that all contacts give rise to linearly independent constraints on interparticle
n
separations. To understand this behavior, we decompose the dynamical matrix M = H −S for
a
static packings of ellipsoidal particles into two important components: the stiffness H and stress
J
S matrices. We find that the stiffness matrix possesses N(ziso−z) eigenmodes eˆ0 with zero eigen-
9
values even at finite compression, where N is the number of particles. In addition, these modes eˆ0
1
are nearly eigenvectors of the dynamical matrix with eigenvalues that scale as ∆φ, and thus finite

compressionstabilizespackingsofellipsoidalparticles. Atjammingonset,theharmonicresponseof
t] static packings of ellipsoidal particles vanishes, and the total potential energy scales as δ4 for per-
f
o turbations by amplitude δ along these ‘quartic’ modes, eˆ0. These findings illustrate the significant
s differences between static packings of spherical and ellipsoidal particles.
.
t
a PACSnumbers: 83.80.Fg61.43.-j,63.50.Lm,62.20.-x
m
-
d I. INTRODUCTION sume that all contacts give rise to linearly independent
n constraintsontheinterparticleseparations. Despitethis,
o these packings were found to be mechanically stable
c There have been many experimental [1–3], computa- (MS) [14, 15].
[ tional[4–6],andtheoretical[7,8]studiesofthestructural
In a recent manuscript [14] by Donev, et al., the au-
and mechanical properties of disordered static packings
1 thors explained this apparent contradiction—that static
v of frictionless disks in 2D and spheres in 3D. In these packings of ellipsoidal particles are mechanically stable,
5 systems,countingarguments,whichassumethatallpar-
yet possess z < z . The main points of the argument
6 ticlecontactsgiverisetolinearlyindependentimpenetra- iso
are included here. The set of N interparticle contacts
9 bility constraints on the particle positions, predict that c
imposes N constraints, f = 1 r /σ 0, where r
3 c ij ij ij ij
the minimum number of contacts required for the sys- − ≤
1. tem to be collectively jammed is N Nmin =N +1, is the center-to-center separation and σij is the contact
0 where N = Nd for fixed bouncd≥ary cconditiodnosf and distance along rˆij between particles i and j. In disor-
dof dered MS sphere packings with z = z , each of the N
2 iso c
N = Nd d for periodic boundary conditions [9, 10],
1 dof − interparticle contacts represents a linearly independent
where d is the spatial dimension and N is the number
v: of particles [11]. The additional contact is required be- constraint. In contrast, some of the Nc contacts for MS
packings of ellipsoidal particles give rise to linearly de-
i cause contacts between hard particles provide only in-
X pendent constraints. Linearly dependent constraints do
equality constraints on particle separations [9]. In the
r not block the degreesof freedom that appear in the con-
a large-system limit, this relation for the minimum num- straint equations for sphere packings, whereas they can
ber of contacts reduces to z z , where z = 2N /N
≥ iso c block multiple degrees of freedom in packings of ellip-
is the average contact number. Disordered packings of
soidal particles because they have convex particle shape
frictionlessspheresaretypicallyisostaticatjammingon-
with a varying radius of curvature [14].
set with z = z , and possess the minimal number of
iso For static packings of spherical particles, interparti-
contacts required to be collectively jammed [9]. Fur-
clecontactsgiverisetoonly“convex”constraints(Fig.1
ther, it has been shown in numerical simulations that
(a)),whilecontactscanyield“convex”or“concave”con-
collectively jammed hard-sphere packings correspond to
straintsinpackingsofellipsoidalparticles(Figs.1(a)and
mechanically stable soft-sphere packings in the limit of
(b)). Note that the distinction between convex and con-
vanishing particle overlaps [6, 12, 13].
caveconstraintsis different than the distinction between
In contrast, several numerical [14–17] and experimen- convex and concave particles. For example, ellipsoids
tal studies [18, 19] have found that disordered packings have a convex particle shape, but static ellipsoid pack-
of ellipsoidal particles possess fewer contacts (z < z ) ings can possess concave interparticle constraints.
iso
than predicted by naive counting arguments, which as- In jammed packings of ellipsoidal particles, there are
2
u u even at finite compression (∆φ > 0); and 3) The modes
2 2
eˆ are nearly eigenvectors of the dynamical matrix (and
0
thestressmatrix S)witheigenvaluesthatscaleasc∆φ,
−
with c > 0, and thus finite compression stabilizes pack-
u u ings of ellipsoidal particles [15]. In contrast, for static
1 1
packings of spherical particles, the stiffness matrix con-
tributions to the dynamical matrix stabilize all modes
near jamming onset. At jamming onset (∆φ = 0), the
(a) (b) harmonicresponseofpackingsofellipsoidalparticlesvan-
ishes, and the total potential energy scales as δ4 for per-
turbations by amplitude δ along these ‘quartic’ modes,
FIG. 1: Schematic diagram that illustrates locally (a) con- eˆ . Our findings illustrate the significant differences be-
0
vex (positive radius of curvature) and (b) concave (negative
tween amorphous packings of spherical and ellipsoidal
radius of curvature) interparticle constraints in packings of
particles.
hard ellipsoidal particles. Inaccessible regions (with f >0)
ij
are shaded blue. The axes labeled u1 and u2 indicate two The remainder of the manuscript will be organized as
representative orthogonal directions in configuration space.
follows. In Sec. II, we describe the numerical methods
thatwe employedto measureinterparticle overlaps,gen-
eratestatic packings,andassessthe mechanicalstability
of packings of ellipsoidal particles. In Sec. III, we de-
N(z z) special directions eˆ in configuration space
iso − 0 scriberesultsfrommeasurementsofthedensityofvibra-
along which perturbations give rise to interparticle over-
tionalmodesintheharmonicapproximation,thedecom-
laps that scale quadratically with the perturbation am-
positionofthe dynamicalmatrix eigenvaluesinto contri-
plitude, f δ2 [14, 15]. Displacements in all other di-
ij ∼ butions from the stiffness and stress matrices, and the
rectionsyield overlapsthatscalelinearly with δ, f δ,
ij ∼ relative contributions of the translational and rotational
as found for jammed sphere packings. This novel scal-
degrees of freedom to the vibrational modes as a func-
ing behavior for packings of ellipsoidal particles can be
tion of overcompression and aspect ratio using several
explained by decomposing the dynamical matrix M =
packing-generation protocols. In Sec. IV, we summarize
H S for these packings into two important compo-
− our conclusions and provide promising directions for fu-
nents: the stiffness matrix H that contains all second-
ture research. We also include two appendices. In Ap-
order derivatives of the total potential energy V with
pendix A, we show that the formation of new interparti-
respect to the configurational degrees of freedom, and
cle contacts affects the scaling behavior of the potential
the stress matrix S that includes all first-order deriva-
energy with the amplitude of small perturbations along
tives of V with respect to the particle coordinates. The
eigenmodes of the dynamical matrix. In Appendix B,
directions eˆ are the eigenvectors of the stiffness matrix
0 we provide analytical expressions for the elements of the
H with zero eigenvalues.
dynamical matrix for ellipse-shaped particles in 2D that
For static packings of ellipsoidal particles at the jam- interact via a purely repulsive linear spring potential.
ming threshold (∆φ = 0) that interact via purely repul-
sive linear spring potentials (i.e. V f2), we find that
∼ ij
the total potential energy increases quartically when the
systemisperturbedbyδalongtheeˆ directions,V cδ4,
0
∝
where the constant c > 0. Also, at the jamming thresh-
old, the stress matrix S =0 and zero modes of the stiff-
II. METHODS
ness matrix are zero modes of the dynamical matrix.
In this manuscript, we will investigate how the me-
chanical stability of packings of ellipsoidal particles is In this section, we describe the computational meth-
modified at finite compression (∆φ > 0). For example, ods employed to generate static packings of convex,
when a system at finite ∆φ is perturbed by amplitude δ anisotropic particles, i.e. ellipses in 2D and prolate el-
alongeˆ0,doquadratictermsinδ ariseinthetotalpoten- lipsoids in 3D with aspect ratio α=a/b of the major to
tialenergyordothe contributionsremainzerotosecond minor axes (Fig. 2), and analyze their mechanical prop-
order? If quadratic terms are present, do they stabilize erties. To inhibit ordering in 2D, we studied bidisperse
or destabilize the packings, and how do the lowest fre- mixtures (2-to-1 relative number density), where the ra-
quency modes of the dynamical matrix scale with ∆φ tio of the major (and minor) axes of the large and small
and aspect ratio? particles is a /a = b /b = 1.4. In 3D, we focused on a
l s l s
We find a number of key results for static packings of monodisperse size distribution of prolate ellipsoids. We
ellipsoidal particles at finite compression (∆φ > 0) in- employed periodic boundaries conditions in unit square
cluding: 1) Packings of ellipsoidal particles generically (2D)andcubic(3D)cellsandstudiedsystemssizesinthe
satisfy z < z [14–17]; 2) The stiffness matrix H pos- range from N =30 to 960 particles to address finite-size
iso
sesses N(z z) eigenmodes eˆ with zero eigenvalues effects.
iso 0
−
3
a (a) (b) (c) (d)
FIG.4: Ellipseswithα=2positionedattheGay-Bernecon-
b tact distance σa. For two ellipses with the same size, the
ij
(a) end-to-end configuration is exact, while the (b) side-to-
end configuration has a 5% relative error. For two ellipses
(a) with a /a = 1.4, the (c) end-to-end configuration has a rel-
j i
ativeerror of 1%, while the(d) side-to-end configuration has
a relative error of 10%.
a
scales will be expressed in units of ǫ, l = I/m, and
l m/ǫ, respectively, where m and I are the mass and
b p
moment of inertia of the ellipsoidal particles.
p
Perram and Wertheim developed an efficient method
for calculating the exact contact distance between ellip-
(b)
soidal particles with any aspect ratio and size distribu-
tionin2D and3D[28]. Intheir formulation,the contact
distance is obtained from
FIG. 2: Wefocus on (a) ellipses in 2D with aspect ratio α=
a/b defined as the ratio of the major to minor axis and (b)
prolate ellipsoids in 3D where α is the ratio of the polar to σij = minσij(λ), (2)
λ
equatorial lengths.
σ0(λ)
σ (λ) = ij ,
ij
(β(λ)rˆ µˆ β(λ)−1rˆ µˆ )2
1 χ(λ) ij · i± ij · j
µˆi s − 2 ± 1±χ(λ)µˆi·µˆj
X
1 b2 b2
i rˆij µˆj σi0j(λ) = 2sλi + 1 jλ,
σ −
ij 1/2
a2 b2 a2 b2
j χ(λ) = i − i j − j ,
 a2+(cid:0) 1−λb2(cid:1)(cid:0)a2+ λ(cid:1) b2 
j λ i i 1−λ j
(cid:0) (cid:1)(cid:16) 1(cid:17)/4
FpaIGrt.ic3le:sDieafinnditjiownitohfutnhietcvoencttoarcstµˆdisatnadncµˆe σtihjaftocrhealrliapcstoeirdizael β(λ) = a2i −b2i a2j + 1−λλb2i .
theorientationsoftheirmajoraxes.iσ isthjecenter-to-center (cid:0)a2 b2(cid:1)(cid:0)a2+ λ b2(cid:1)
ij j − j i 1−λ j
separationatwhichellipsoidalparticlesfirsttouchwhenthey (cid:0) (cid:1)(cid:16) (cid:17)
are brought together along rˆ at fixed orientation.
ij
The approximation σa = σ (λ = 1/2) is equivalent
ij ij
to the commonly used Gay-Berne approximation for the
contactdistance[27,29]. TheaccuracyoftheGay-Berne
A. Contact distance approximationdepends onthe relativeorientationofthe
two ellipsoidal particles, and in general is more accurate
for monodisperse systems. For example, in Fig. 4, we
In both 2D and 3D, we assume that particles interact
showσa forseveralrelativeorientationsofbothmonodis-
via the following pairwise, purely repulsive linear spring ij
perseandbidispersesystems. Therelativedeviationfrom
potential
the true contact distance can be as large as e 10% for
∼
2 a /a = 1.4 and α = 2. Thus, the Gay-Berne approxi-
ǫ 1 rij r σ j i
V(rij)=2 − σij ij ≤ ij (1) mationshouldbeusedwithcautionwhenstudyingpoly-
0 (cid:16) (cid:17) r >σ , disperse packings of ellipsoidal particles. For monodis-
 ij ij
perse ellipses with α = 2, 0% < e < 5%. We find sim-
whereǫisthecharacteristicenergyoftheinteraction,rij ilar results for 3D systems. Unless stated otherwise, we
isthecenter-to-centerseparationbetweenparticlesiand employ the exact expression for contact distance, and
j, and σ is the orientation-dependent center-to-center thus σ = σ (λ ), β = β(λ ), χ = χ(λ ), and
ij ij ij min min min
separation at which particles i and j come into contact σ0 = σ0(λ ), where λ is the minimum obtained
ij ij min min
as shown in Fig. 3. Below, energies, lengths, and time from Eq. 2.
4
B. Packing generation protocol axis of the ellipsoidal particles. In 2D, we solve
d2~r
m i = F~ b ~v (3)
dt2 ij − r i
We employ a frequently used isotropic compression i>j
X
method for soft, purely repulsive particles [30, 31] to d2θ
I i = T b θ˙ , (4)
generate static packings of ellipsoidal particles at jam- dt2 ij − θ i
ming onset (∆φ = 0). Static packings at jamming on- Xi>j
set are characterized by infinitesimal but nonzero total
where θ is the angle the long axis of ellipse i makes
potential energy and pressure. This isotropic compres- i
with the horizontal axis, ~v is the translational velocity
sion method consists of the following steps. We be- i
of particle i, θ˙ is the rotational speed of particle i, b
gin by randomly placing particles at low packing frac- i r
and b are the damping coefficients for the position and
tion (φ = 0.2) with random orientations and zero ve- θ
0 angle degrees of freedom, and the moment of inertia I =
locities in the simulation cell. We successively com-
m(a2+b2)/4. The force F~ on ellipse i arising from an
press the system by small packing fraction increments ij
δφ=10−3, with eachcompressionfollowedby conjugate overlapwith ellipse j is
gradient(CG)energyminimizationuntilthetotalpoten-
rˆ ∂lnσijψˆ
tial energy per particle drops below a small threshold, F~ = F ij − ∂ψij ij , (5)
Vpa/rNtic≤lebVettowl e=en1s0u−c1c6e,ssoirvethiteertaottiaolnspootfetnhteiamliennimerigzyatpioenr ij − ij 1+ ∂∂lnψσijij 2
q
routine is (V V )/V V . The algorithm switches (cid:0) (cid:1)
t+1− t t ≤ tol where
from compressionto decompression if the minimized en-
ergy is greater than 2V . Each time the algorithm tog-
tol ∂lnσ 2 ∂V(r )
glesfromcompressiontodecompressionorviceversa,the F = 1+ ij ij , (6)
ij
packing fraction increment is halved. s (cid:18) ∂ψij (cid:19) (cid:12) ∂rij (cid:12)
∂lnσij = χ σi2j (β2 χ(cid:12)(cid:12)(cid:12))sin[2(ψ(cid:12)(cid:12)(cid:12) θ )]+
The packing-generationalgorithm is terminated when ∂ψ −2 (σ0)2 − ij − i
ij ij
the total potential energy per particle satisfies Vtol < (β−2 χ)(cid:2)sin[2(ψ θ )]
V/N < 2Vtol. Thus, using this method we can locate − ij − j ×
thejammedpackingfractionφJ andparticlepositionsat 1 χ2cos2[θi θj] −1, (cid:3) (7)
jamming onset for each initial condition to within 10−8. − −
After jamming onset is identified, we also generate con- dV(r )/dr (cid:0)= σ−1(1 rij) for(cid:1)the purely repulsive
figurations at specified ∆φ = φ φ over six orders of − ij ij ij − σij
magnitude from 10−8 to 10−2 by−apJplying a prescribed linearspringpotentialinEq.1, andrˆij andψˆij areillus-
trated in Fig. 5.
set of compressions with each followed by energy mini-
To calculate the torque T = [r~c F~ ] zˆ in Eq. 4,
mization. ij ij × ij ·
we must identify the point of contactbetween particles i
and j,
Todetermine whetherthe accuracyofthe energymin-
imization algorithm affects our results (see Sec. IIID), b 1
~rc = i
we calculate the eigenvalues of the dynamical matrix as ij 2 α−2+tan2τ ×
ij
a function of the total kinetic energy (or deviation from
zeroinforceandtorquebalanceoneachparticle)ateach (cposθi−sinθitanτij)xˆ+
∆φ. To do this, we initialize the system with MS pack- ((cid:2)sinθitanτij +cosθi)yˆ , (8)
ings from the above packing-generation algorithm and
usemoleculardynamics(MD) simulationswithdamping where (cid:3)
termsproportionaltothetranslationalandrotationalve-
tan(ψ θ ) ∂lnσij
locitiesoftheellipsoidalparticlestoremoveexcesskinetic tanτ = α−2 ij − i − ∂ψij (9)
energy from the system [32]. The damped MD simula- ij 1+tan(ψ θ )∂lnσij
tions are terminated when the total kinetic energy per ij − i ∂ψij
particle is below K/N =Ktol, where Ktol is varied from and~rc , ψ , and τ are depicted in Fig. 5. From Eqs. 5
10−16 to 10−32. This provides accuracy in the particle ij ij ij
and 8, we find
positionsoftheenergyminimizedstatesintherangefrom
10−8 to 10−16. b F
i ij
T =
ij
− 2 ×
For the damped MD simulations, we solve Newton’s 1 α−2 tanτ
ij
− . (10)
equations of motion (using fifth-order Gear predictor-
1+α−4(cid:0)tan2τ (cid:1)1+α−2tan2τ
corrector methods [33]) for the center of mass position ij ij
and angles that characterize the orientation of the long q(cid:0) (cid:1)(cid:0) (cid:1)
5
j derivative involving angles, Fi = ∂V(r )/∂θ , where
θ − ij i
μˆj θj Fθi = 14χ σσi0j 2(2αA(B++B−)+χC(B+2 −B−2)), (12)
τji (cid:18) ij(cid:19)
y cosθ x sinθ
ij i ij i
r A = − ,
ij r
rc ij
ji Fˆji α(x cosθ +y sinθ ) α−1(x cosθ +y sinθ )
ij i ij i ij j ij j
rc B± = ± ,
ij Fˆij (1+χcos[θi−θj])rij
ψ C = cos2(θ θ ).
ψˆ ijτij μˆi i− j
ij θi Complete expressions for the matrix elements of the dy-
i namical matrix for ellipses in 2D are provided in Ap-
pendix B. In 3D, we calculated the first derivatives of V
with respectto the particle coordinatesanalytically,and
then evaluated the second derivatives for the dynamical
matrix numerically.
The vibrational frequencies in the harmonic approxi-
FIG. 5: Geometry of two ellipses in contact. The angles θ
i
mationcanbeobtainedfromtheNd dnontrivialeigen-
andθ characterizetheorientationofparticlesiandjrelative f
j −
tothehorizontalaxis,i.e. µˆi =cosθixˆ+sinθiyˆ. ψij givesthe values mi of the dynamical matrix, ωi = mi/ǫbs. d of
angle between the center-to-center separation vector ~rij and the eigenvaluesare zerodue to periodic boundary condi-
thehorizontalaxisandψˆ =−sinψ xˆ+cosψ yˆistheangle tions. For all static packings, we have veprified that the
ij ij ij
unit vector in polar coordinates. The unit vector Fˆ =−Fˆ smallest nontrivial eigenvalue satisfies m /N >10−10.
ij ji min
pointsinthedirectionoftheforceonparticleiduetoparticle Below, we will study the density of vibrational fre-
jatthepointofcontact. ~ricj pointsfromthecenterofparticle quencies D(ω) = (N(ω +∆ω) N(ω))/(Ndof∆ω) as a
i to the point of contact with particle j, and τij is the angle function of compression ∆φ an−d aspect ratio α, where
between µˆ and~rc .
i ij N(ω) is the number of vibrational frequencies less than
ω. We will also investigate the relative contributions
of the translational and rotational degrees of freedom
to the nontrivial eigenvectors of the dynamical matrix,
C. Dynamical matrix calculation mˆ = mj=1,mj=1,mj=1,...,mj=N,mj=N,mj=N
i { xi yi θi xi yi θi }
for ellipses in 2D and mˆ =
i
Toinvestigatethemechanicalpropertiesofstaticpack- mj=1,mj=1,mj=1,mj=1,mj=1,...,mj=N,mj=N,mj=N,
ings of ellipsoidal particles, we will calculate the eigen- { xi yi zi θi φi xi yi zi
mj=N,mj=N for prolateellipsoids in3D, wherei labels
values of the dynamicalmatrix and the resulting density θi φi }
the eigenvector and runs from 1 to Nd d. The
ofvibrationalmodesintheharmonicapproximation[34]. f −
eigenvectors are normalized such that mˆ2 =1.
The dynamical matrix is defined as i
∂2V
M = , (11) D. Dynamical matrix decomposition
kl
∂u ∂u
k l
The dynamical matrix (Eq. 11) can be decomposed
where u (with k = 1,...,d N) represent the d N
k f f into two component matrices M = H S: 1) the stiff-
degreesoffreedominthesystemanddf isthenumberof nessmatrixH thatincludesonlysecond−-orderderivatives
degrees of freedom per particle. In 2D d = 3 with ~u =
f of the total potential energy V with respect to the con-
x ,x ,...,x ,y ,y ,...,y ,l θ ,l θ ,...,l θ and in 3D
{ 1 2 N 1 2 N 2 1 2 2 2 N} figurational degrees of freedom and 2) the stress matrix
for prolate ellipsoids d = 5 with ~u = x ,x ,...,x ,y ,
y ,...,y ,z ,z ,...,z ,l1fθ ,l2θ ,...,lNθ ,{l φ1 ,l2φ ,..N., 1 S that includes only first-order derivatives of V. The kl
2 N 1 2 N θ 1 θ 2 θ N 3 1 3 2 elements of H and S are given by
l φ , where θ is the polar angle and φ is
3 N i i
}
the azimuthal angle in spherical coordinates, ∂2V ∂(r /σ )∂(r /σ )
ij ij ij ij
H = (13)
l2 = √a2+b2/2, l3 = a2+b2 /5, and kl ∂(rij/σij)2 ∂uk ∂ul
i>j
X
lθi = 2b2+(a2−b2)sin2φi /5. q(cid:0) (cid:1) S = ∂V ∂2(rij/σij), (14)
kl
Thqed(cid:0)ynamicalmatrixrequi(cid:1)rescalculationsofthe first −i>j ∂(rij/σij) ∂uk∂ul
and second derivatives of the total potential energy V X
withrespecttoallpositionalandangulardegreesoffree- wherethesumsareoverdistinctpairsofoverlappingpar-
dominthesystem. ThefirstderivativesofV withrespect ticles i and j. Since ∂2V/∂(r /σ )2 = ǫ for the purely
ij ij
to the positions of the centers of mass of the particles~r repulsivelinearspringpotential(Eq.1),thestiffnessma-
i
canbeobtainedfromEq.5. In2D,thereisonlyonefirst trix depends only on the geometry of the packing (i.e.
6
∂(r /σ )/∂u ). Also, at zero compression ∆φ = 0, 6.5
ij ij k
S =0, M =H,andonly the stiffness matrixcontributes
tothedynamicalmatrix. Thefrequenciesassociatedwith 6
the eigenvalues h of the stiffness matrix (at any ∆φ)
i
are denoted by ωhi = hi/ǫbs, and the stiffness matrix 5.5
eigenvectors are normalized such that ˆh2 =1.
p i
5
z

E. Contact number
4.5
When counting the number of interparticle contacts
4 (a)
N ,we removeallrattlerparticles (definedas those with
c
fewerthand+1contacts)anddonotincludethecontacts
thatrattlerparticlesmakewithnon-rattlerparticles[35]. 3.5
1 1.2 1.4 1.6 1.8
Removing these contacts may cause non-rattlerparticles
α
to become rattlers, and thus this process is performed
recursively. Note that for ellipsoidal particles with d+1
11
contacts,thelinesnormaltothe points(orplanesin3D)
of contactmust allintersect, otherwise the system is not
mechanicallystable. Thenumberofcontactsperparticle 10
is defined as z =N /(N N ), where N is the number
c r r
−
of rattlers. We find that the number of rattler particles 9
decreaseswithaspectratiofromapproximately5%ofthe
system at α=1 to zero for α>1.2 in both 2D and 3D. 8
z

7
III. RESULTS
6 (b)
Static packings of ellipsoidal particles at jamming on-
set typically possess fewer contacts than predicted by
isostatic counting arguments [14], z < z , over a wide 5
iso
1 1.2 1.4 1.6 1.8 2
range of aspect ratio as shown in Fig. 6. This finding
α
raises a number of important questions. For example,
are static packings of ellipsoidal particles mechanically
stable at finite ∆φ > 0, i.e. does the dynamical ma- FIG. 6: Average contact number z versus aspect ratio α for
trix for these systems possess nontrivial zero-frequency staticpackingsof(a)bidisperseellipsesin2Dand(b)prolate
modes at ∆φ > 0? In this section, we will show that ellipsoidsin3Datjammingonset. Theisostaticvaluesziso=
6 (2D) and 10 (3D) are indicated by dashed lines.
packings of ellipsoidal particles are indeed mechanically
stable (with no nontrivial zero-frequency modes) by cal-
culating the dynamical, stress, and stiffness matrices for
these systems as a function of compression ∆φ, aspect
ratio α, and packing-generation protocol. Further, we amorphouspackingsofellipsoidalparticlesandshowthat
will show that the density of vibrationalmodes for these the density of vibrationalmodes forthese systemsshows
systems possesses three characteristic frequency regimes significant qualitative differences from that for spherical
anddeterminethescalingofthesecharacteristicfrequen- particles.
cies with ∆φ and α. In Figs. 7 (a) and (b), we show D(ω) on linear and
log-log scales,respectively, for ellipse-shaped particles in
2D at ∆φ = 10−8 over a range of aspect ratios from
A. Density of vibrational frequencies D(ω) α = 1 to 2. We find several key features in D(ω): 1)
For low aspect ratios α<1.05, D(ω) collapses with that
A number of studies have shown that amorphous for disks (α = 1) at intermediate and large frequencies
sphere packings are fragile solids in the sense that the 0.25<ω <2.25;2)Forlargeaspectratiosα 2,D(ω)is
≥
density of vibrational frequencies (in the harmonic ap- qualitatively different for ellipses than for disks over the
proximation)D(ω) for these systems possesses anexcess entire frequency range;and 3) A strong peak near ω =0
of low-frequency modes over Debye solids near jamming andasmallersecondarypeakatintermediatefrequencies
onset, i.e. a plateauforms andextends tolowerfrequen- (evident onthe log-logscale in Fig.7 (b)) occur in D(ω)
cies as ∆φ 0 [6, 36, 37]. In this work, we will cal- for α > 1. Note that at finite compression ∆φ > 0, we
→
culate D(ω) as a function of ∆φ and aspect ratio α for do not find any nontrivial zero-frequency modes of the
7
2 (a) 0.5 (a)
0.4
1.5
)
ω
D( 1 ω)0.3
(
D

0.2
0.5
0.1
0
0 0.5 1 1.5 2 2.5
ω
0
0 0.5 1 1.5 2 2.5 3
4 1 2 3 ω
)
4 ω
( 2
D
3
0 ) 1 2 3
1 0 ω
g 2
o D(
) 2 l -2
ω 2 0
D( -6 -4 -2 0 g10
0 log ω lo
g1 10
lo 0 ω) 1 -2
(
D -4 -2 0
10 0 log10 ω
-2 (b) g
o
l
-6 -4 -2 0 -1
log ω
10
(b)
-2
FIG. 7: (a) The density of vibrational frequencies D(ω) for
N = 240 ellipse-shaped particles at ∆φ = 10−8 with aspect
-5 -4 -3 -2 -1 0
ratio α = 1.0 (solid), 1.001 (dotted), 1.05 (dashed), and 2.0
log ω
(dot-dashed). D(ω)forα=1hasbeenscaled by2/3relative
10
to those with α> 1 to achieve collapse at low aspect ratios.
(b) D(ω) for the same aspect ratios in (a) on a log-log scale.
The inset illustrates the three characteristic frequencies ω1, FIG. 8: (a) The density of vibrational frequencies D(ω) for
ω2, and ω3 in D(ω) for α=1.001. N =512 prolate ellipsoids at ∆φ=10−6 for α=1.0 (solid),
1.001 (dotted), 1.005 (dashed), and 1.2 (dot-dashed). D(ω)
forα=1hasbeenscaled by3/5relativetothosewithα>1
toachievecollapseatlowaspectratios. (b)D(ω)forthesame
dynamical matrix in static packings of ellipses and ellip-
aspect ratios in (a) on a log-log scale. The inset illustrates
soids. The only zero-frequency modes in these systems
the three characteristic frequencies ω1, ω2, and ω3 in D(ω)
correspondtothedconstanttranslationsthatarisefrom
for α=1.001.
periodic boundary conditions and zero-frequency modes
associated with ‘rattler’ particles with fewer than d+1
interparticle contacts.
at large frequencies. All three characteristic frequencies
To monitor the key features of D(ω) as a function of
increase with aspect ratio. Note that we only track ω
∆φ and α, we define three characteristic frequencies as 2
and ω for aspect ratios where ω < ω . For example,
shown in the inset to Fig. 7 (b). ω and ω identify the 3 2 3
1 2
theintermediateandhigh-frequencybandscharacterized
locations of the small and intermediate frequency peaks
by ω and ω merge for α 1.2.
in D(ω), and ω marks the onset of the high-frequency 2 3
3 ≥
plateauregimeinD(ω). Forouranalysis,wedefineω as As showninFig. 8, D(ω) for 3D prolateellipsoids dis-
3
the largest frequency (< 1) with D(ω) < 0.15, which is plays similar behavior to that for ellipses in 2D (Fig. 7)
approximately half of the height of the plateau in D(ω) for aspect ratios α < 1.5. For example, D(ω) for el-
8
3 (a) 0 (a)
2 -1
ω) 1 ω -2
D( 10
g
0 0 o -3
1 l
g
o
l
-1 -4
-2 -5
-3 -2 -1 0
-3 log (α−1)
-6 -5 -4 -3 -2 -1 0 10
log ω
10
(b)
0.5
2 (b)
)
φ1/2 0
Δ
1 ω /1
( -0.5
ω) 0
1
D( 0 og
l
0 -1
1
g
o -1
l
-1.5
-2
-4 -3 -2 -1 0 1
log (α−1)
10
-3
-6 -5 -4 -3 -2 -1 0
log ω FIG. 10: (a) Characteristic frequencies ω1 (circles), ω2
10 (squares), and ω3 (diamonds) from D(ω) as a function of as-
pect ratio α−1 for N = 240 ellipses in 2D at ∆φ = 10−8.
FIG. 9: The density of vibrational frequencies D(ω) for N = The solid (dashed) lines have slope 1/2 (1). (b) ω1/(∆φ)1/2
240 ellipses as a function of compression ∆φ = 10−7 (solid), for systems with N = 240 ellipses in 2D at ∆φ = 10−7 (cir-
10−5 (dashed), 10−3 (dotted), and 10−2 (dot-dashed) for (a) cles), 10−6 (squares), 10−5 (diamonds), 10−4 (upward trian-
α=1.05 and (b) 2. gles), 10−3, (downward triangles), and 10−2 (crosses). The
solid line has slope 1/2.
lipsoids possesses low, intermediate, and high frequency
regimes, whose characteristic frequencies ω , ω , and ω
1 2 3
packings as a function of compression∆φ for two aspect
increase with aspect ratio. Note that the intermediate
ratios, α = 1.05 and 2. We find that the low-frequency
andhigh-frequencybandsω andω mergefor α>1.02,
2 3
band(characterizedbyω )depends on∆φ,while thein-
which occurs at lower aspect ratio than the merging of 1
termediate and high frequency bands do not. The inter-
the bands in 2D. Another significantdifference is that in
mediate and high frequencies bands do not change sig-
3D D(ω) extends to higher frequencies at large aspect
nificantly until the low-frequency band centered at ω
ratios (α&1.2) than D(ω) for ellipses. 1
mergeswiththemat∆φ 10−3and 10−4forα=1.05
We note the qualitative similarity between the D(ω) ≈ ≈
and 2, respectively.
forα=1.005ellipsoidsshownin Fig.8(b) andD(ω)for
α=0.96presentedinFig. 1(c)ofRef.[16]forω >10−2. We plot the characteristic frequencies ω , ω , and ω
1 2 3
However,Zeravcic, et al. suggest that there is no weight versus aspect ratio α 1 for ellipse packings in Fig. 10
in D(ω) for ω < 10−2 except at ω = 0 for both oblate and ellipsoid packings−in Fig. 11. The characteristic fre-
andprolateellipsoids,incontrasttoourresultsinFig.8. quenciesobeythefollowingscalinglawsoveratleasttwo
In Fig. 9, we show the behavior of D(ω) for ellipse ordersofmagnitudeinα 1andfiveordersofmagnitude
−
9
0
(a) 6
-1
4
)
h
ω
ω
(
g10-2 D0 2
o 1
l g
o
l
0
-3
-2
-4
-3 -2 -1
-10 -8 -6 -4 -2 0
log (α-1)
10 log ω
10 h
0.6 (b) FIG. 12: The distribution of frequencies D(ωh) associated
with the eigenvalues of the stiffness matrix H for N = 240
ellipsepackingsasafunctionofcompression∆φ=10−5 (dot-
φ)1/2 0.3 ted),10−3(dashed),and10−2(dot-dashed)forα=1.05. The
Δ verticalsolidlineindicatesthe‘zero-frequency’toleranceωtol,
ω /1 whichisthelowestfrequencyobtainedforthedynamicalma-
trix for packings at α = 1.05 and the smallest compression
( 0
0 (∆φ=10−8) in Fig. 7.
1
g
o
l
-0.3
B. Dynamical Matrix Decomposition
-0.6
As showninFig.6,static packingsofellipsoidalparti-
-3 -2 -1 clescanpossessz <zisooverawiderangeofaspectratio,
log (α-1) yet as described in Sec. IIIA, the dynamical matrix M
10
contains a complete spectrum of Nd d nonzero eigen-
f
−
values m near jamming. To investigate this intriguing
i
FIG. 11: (a) Characteristic frequencies ω1 (circles), ω2 property, we first calculate the eigenvalues of the stiff-
(squares), and ω3 (diamonds) from D(ω) as a function of nessmatrixH,showthatitpossessesN ‘zero’-frequency
aspect ratio α−1 for N = 240 prolate ellipsoids in 3D at z
∆φ=10−6. The solid (dashed) lines have slope 1/2 (1). (b) modes whose number matches the deviation in the con-
ω1/(∆φ)1/2 for systems with N = 240 prolate ellipsoids at tact number from the isostatic value, and then identify
∆φ = 10−6 (circles), 10−5 (squares), and 10−4 (diamonds). the separate contributions from the stiffness and stress
The solid line has slope 1/2. matrices to the dynamical matrix eigenvalues.
In Fig. 12, we show the distribution of frequencies
D(ω )associatedwiththeeigenvaluesofthestiffnessma-
h
in ∆φ: trixforellipsepackingsatα=1.05asafunctionofcom-
pression∆φ. WefindthreestrikingfeaturesinFig.12: 1)
ω (∆φ)1/2(α 1)1/2, (15) Many modes of the stiffness matrix existnear and below
1
∼ − the zero-frequency threshold (determined by the vibra-
ω (α 1), (16)
2 ∼ − tional frequencies of the dynamical matrix at α = 1.05
ω3 (α 1)1/2. (17) and ∆φ = 10−8), 2) Frequencies that correspond to the
∼ −
low-frequency band characterized by ω are absent, and
1
Similar results for the scaling of ω and ω with α 1 3) The nonzero frequency modes (with ω > 10−2) do
2 3 h
−
were found in Ref. [16]. We will refer to the modes in not scale with ∆φ as pointed out for the dynamical ma-
the low-frequency band in D(ω) (with characteristic fre- trix eigenvalues in Eqs. 16 and 17. Further, we find that
quencyω )as‘quarticmodes’,andthesewillbediscussed the number of zero-frequency modes N of the stiffness
1 z
in detail Sec. IIIC. The scaling of the quartic mode fre- matrix matches the deviation in the number of contacts
quencieswithcompression,ω (∆φ)1/2,hasimportant from the isostatic value (N N ) for each ∆φ and as-
1 iso c
∼ −
consequences for the linear response behavior of ellip- pect ratio. Specifically, N = N N over the full
z iso c
soidal particles to applied stress [15]. range of ∆φ for 99.96% of the more−than 103 packings
10
1 bandscharacterizedbyω andω inFig.7),theeigenval-
2 3
ues of the stiffness and dynamical matrices are approxi-
0 mately the same, ω2. The deviation ω2 = ,
0 ∆ H≈ i i −H −S
-0.4- / Tshhouwsn, winetfihnedintsheatttothFeigin.t1e3rm(ae)d,isactaeleasnldinheaigrhlyfwreiqthue∆ncφy.
-0.8 modes for packings of ellipsoidal particles are stabilized
-1
by the stiffness matrix H.
-4 -2 0
g10 log ω2 In the main panel of Fig. 13 (b), we show that for fre-
o 10 quencies in the lowest frequency band (characterized by
l -2
ω ) the eigenvalues of the stress and dynamical matrices
1
are approximately the same, ω2. In the inset to
-3 Fig. 13 (b), we show that the−dSevi≈atioin ωi2−(−S) = H
(a) scales as (∆φ)2. Thus, we find that the lowestfrequency
modes for packings of ellipsoidal particles are stabilized
-4 bythestressmatrix S overawiderangeofcompression
-4 -3 -2 -1 0 1 −
∆φ. Similar results were found previously for packings
log ω2 of hard ellipsoidal particles [14]. In contrast, for static
10
packings of spherical particles, the stress matrix contri-
butions to the dynamical matrix are destabilizing with
-2
1 2) −S < 0 for all frequencies near jamming, and H stabi-
lizes the packings as shown in Fig. 14.
∆
-4 0 ( /0
-1g1
o C. Quartic modes
) l
-2
- -5 -4 -3
(
g 10 -6 log10 (ω2/∆ ) MWfoersphaocwkeidngisnoSfeecl.lipIIsIoAidatlhpaatrttihcelesdycnoanmtaiicnasl amcaotrmix-
o
l plete spectrum of Ndf d nonzero eigenvalues mi for
−
∆φ > 0 despite that fact that z < z . Further, we
-8 iso
showed that the modes in the lowest frequency band
(b) scale as ω1 (∆φ)1/2 in the ∆φ 0 limit. What hap-
∼ →
pens at jamming onset when ∆φ = 0, i.e. are these
-10 low-frequency modes that become true zero-frequency
-10 -8 -6 -4 -2
modes at ∆φ = 0 stabilized or destabilized by higher-
log ω2 order terms in the expansion of the potential energy in
10
powers of the perturbation amplitude?
FIG. 13: The two contributions to the dynamical matrix To investigate this question, we apply the following
eigenvalues, (a) H and (b) −S, plotted versus ω2 = H−S deformation to static packings of ellipsoidal particles:
for ellipse packings with α = 1.05 and ∆φ = 10−6 (circles),
10−5 (squares),10−4 (diamonds),and10−3 (triangles). In(a)
~u=~u +δmˆ , (18)
and (b), the solid lines correspond to H = ω2 and S = ω2. 0 i
In the main panel and inset of (a), only modes correspond-
where δ is the amplitude of the perturbation, mˆ is an
ing the intermediate and high frequency bands are included. i
eigenvectorof the dynamical matrix, and ~u is the point
In the main panel and inset of (b), only modes correspond- 0
ing the low frequency band are included. The insets to (a) inconfigurationspacecorrespondingtotheoriginalstatic
and(b),whichplot−S/(∆φ)2versusω2andH/(∆φ)2 versus packing,followedbyconjugategradientenergyminimiza-
ω2/∆φ,showthedeviationsω2−H=−S ∝∆φforhigh-and tion. We then measure the change in the total potential
intermediate-frequency modes and ω2−(−S) = H ∝ (∆φ)2 energy before and after the perturbation, ∆V.
for low-frequency modes. We plot ∆V/N versus δ in Fig. 15 for perturbations
along eigenvectors that correspond to the smallest non-
trivialeigenvaluem =ω2 ofthe dynamicalmatrixfor
1 min
for aspect ratio α < 1.1 and for 100% of the more than
static packings of (a) disks and ellipses and (b) spheres
103 packings for α 1.1. and prolate ellipsoids at ∆φ = 10−7. As expected, for
≥
In Fig. 13, we calculate the projection of the dynam- disksandspheres,we findthat∆V/N mω2 δ2 overa
ical matrix eigenvectors mˆi onto the stiffness and stress widerangeofδ inresponseto perturba≈tionsmalionngeigen-
matrices,H=mˆ†iHmˆi andS =mˆ†iSmˆi,wheremˆ†i isthe vectors that correspond to the smallest nontrivial eigen-
transpose of mˆ and ω2 = mˆ†Mmˆ = . Fig. 13 value. In contrast, we find novel behavior for ∆V/N
(a) shows that fior largeieigenvialuesiω2 oHf t−heSdynamical when we apply perturbations along the eigendirection
i
matrix (i.e. within the intermediate and high frequency that corresponds the lowest nonzero eigenvalue of the