Table Of Content

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 arguments, 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îj 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 contributions 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 perturbations 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 repulsive 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 threshold, 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 mechanical 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 µî s − 2 ± 1±χ(λ)µî·µˆj X 1 b2 b2 i rîj µˆ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µê σ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 energy 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îj andψîj 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îj (1+χcos[θi−θj])rij ψ C = cos2(θ θ ). ψˆ ijτij μî 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. µî =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îstheangle 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 onset 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 calculating 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 aspect 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 contributions to the dynamical matrix are destabilizing with -2 1 2) −S < 0 for all frequencies near jamming, and H stabilizes 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î onto the stiffness and stress widerangeofδ inresponseto perturba≈tionsmalionngeigen- matrices,H=mˆ†iHmî andS =mˆ†iSmî,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