Compression Driven Jamming of Athermal Frictionless Spherocylinders in Two Dimensions

I Introduction

In a system of athermal ($T=0$) granular particles with only contact interactions, as the particle packing fraction $\phi$ increases, the system will undergo a sharp transition from a liquid-like state to a rigid but disordered solid state. In the liquid-like state, particles have sufficient room to avoid each other and so there are no contacts, and so no stress, in the system. As $\phi$ increases, particles come into mutual contact. In the disordered solid state, force chains percolate across the system giving it a finite elastic rigidity, and the system supports a finite stress. This transition from a stress-free liquid-like state to a stress-supporting solid state is known as the jamming transition O’Hern et al. (2003); Liu and Nagel (2010). For particles without intergranular friction, and for a given protocol for compacting the system, this jamming transition occurs at a well defined $\phi_{J}$ and the transition is continuous; stress increases continuously from zero as $\phi$ increases above $\phi_{J}$ O’Hern et al. (2003).

Much of the work that has been done to analyze behavior near the jamming transition has been for the simple case of perfectly spherical particles. It is therefore interesting to ask how the jamming transition may be modified if the particles have shapes with a lower rotational symmetry. Recently, several works have considered the cases of non-spherical particles, in particular monodisperse distributions of aspherical ellipsoids Donev et al. (2004a, b); Man et al. (2005); Donev et al. (2007), oblate ellipsoids Donev et al. (2004b); Man et al. (2005); Donev et al. (2007), and prolate ellipsoids Donev et al. (2004b); Man et al. (2005); Donev et al. (2007); Sacanna et al. (2007); Zeravcic et al. (2009); Schreck et al. (2012) in three dimensions (3D), and bidisperse distributions of ellipses Donev et al. (2007); Mailman et al. (2009); Schreck et al. (2012) in two dimensions (2D). Spherocylinders Williams and Philipse (2003); Wouterse et al. (2007); Azéma and Radjaï (2010, 2012) have been used to model rod-shaped particles, and other work has considered cut spheres Wouterse et al. (2007) in 3D. For a review, see Ref. Borzsonyi and Stannarius (2013).

In this work we consider in detail the compression driven jamming of athermal, frictionless, soft-core 2D spherocylinders. By “compression driven” we mean a protocol in which we start with a dilute system of non-overlapping particles, then isotropically shrink the system box to increase the density, passing through the point at which the system jams. A spherocylinder in 2D consists of a rectangle with two circular end caps. We will study the behavior of such spherocylinders as a function of the aspect ratio of rectangular length to end cap diameter. Spherocylinders are unlike ellipses and ellipsoids in that they have parallel flat sides that could in principle lead to configurations in which particles stack in parallel layers. In that sense they share a similarity with the cut spheres in 3D considered by Wouterse et al. Wouterse et al. (2007). We will pay particular attention to the effect of these parallel sides on the nature of elastic vibrational modes at jamming. We will focus on two issues that arise when considering particles with shape anisotropy: (i) do spherocylinders show any orientational ordering as they are compressed through the jamming transition; and (ii) are configurations at jamming isostatic, and if not, what are the characteristics of the unconstrained (to quadratic order in the energy) modes?

Orientational ordering: It has long been known, since the work of Onsager Onsager (1949), that hard-core (no overlaps allowed) rod shaped particles in thermal equilibrium will undergo a liquid to nematic phase transition as the density is increased. Despite the absence of any globally preferred direction in the system, there will be a spontaneous symmetry-breaking transition in which the rods will show a macroscopic alignment in a particular direction. Bolhuis and Frenkel Bolhuis and Frenkel (1997) mapped out the phase diagram for thermalized hard-core spherocylinders in 3D, finding the dependence of the nematic transition as a function of the spherocylinder aspect ratio $\phi_{c}(\alpha)$ (they also found smectic and crystalline transitions). Monte Carlo simulations for thermalized hard-core spherocylinders in 2D Bates and Frenkel (2000); Lagomarsino et al. (2003) similarly observed a transition, upon increasing density, from an isotropic liquid to a nematic phase with algebraically decaying (rather than long range) orientational correlations. Experimental studies of vibrated, but otherwise athermal, elongated grains in 2D have observed several different types of ordered states with both nematic and tetratic order Narayan et al. (2006); Galanis et al. (2006, 2010).

In contrast to the above, one can ask if a system of athermally compressed rod shaped particles, in which there is neither thermal nor mechanical random agitation of the particles, will show any orientational ordering. Will the elastic forces that act between particles as they are compressed into mutual contact cause a spontaneous alignment so that the particles can pack more densely, or will they jam into an orientationally disordered state as the packing fraction $\phi$ increases? Using a configurational statistical mechanics for athermal granular systems with elongated grains, Mounfield and Edwards Mounfield and Edwards (1994) argued that the grains need not order nematically to be minimally compact. Experimental studies of athermal 3D oblate ellipsoids Donev et al. (2004b) and aspherical ellipsoids Donev et al. (2007) at jamming reported only a small nematic order parameter that was interpreted as consistent with no orientational ordering. Experiments on prolate ellipsoids Sacanna et al. (2007) similarly reported no orientational ordering. Simulations on cut spheres in 3D Wouterse et al. (2007), however, did find nematic ordering when the aspect ratio was sufficiently large (i.e., thin disk shaped particles). In this work we will present a detailed investigation as to whether moderately elongated spherocylinders in 2D show any orientational ordering when athermally and isotropically compressed. We will find that they do not.

The remainder of the paper is organized as follows. In Sec. II we discuss the details of our model system and our procedure for slowly compressing the system through jamming. In Sec. III.1 we present our results on the lack of orientational ordering of moderately elongated spherocylinders. In Sec. III.2 we present our results for pressure as a function of packing fraction, and determine the packing fraction at jamming as a function of spherocylinder aspect ratio, $\phi_{J}(\alpha)$, for compression driven jamming in the quasistatic limit. In Sec. III.3 we describe our energy minimization method for constructing mechanically stable states, and show that it is necessary to treat side-to-side contacts, where two spherocylinders contact along their flat edges, carefully. Doing so, we find that spherocylinders are hypostatic for the entire range of aspect ratios we consider. We also note the strong propensity of spherocylinders to have contacts along their flat sides, even for very small $\alpha$, thus suggesting that the $\alpha\to 0$ limit is singular. In Sec. III.4 we analyze the eigenmodes of small vibrations for both nearly circular and moderately elongated spherocylinders near jamming, and relate these results to the hypostaticity of the system. In Sec. IV we summarize our conclusions.

II Model and Simulation Method

A two dimensional spherocylinder consists of a rectangle with two circular end caps. We will denote the half length of the rectangular part of spherocylinder $i$ as $A_{i}$. The radius of the end cap, which is also the half width of the rectangle, we denote as $R_{i}$, as illustrated in Fig. 1(a). We will refer to the “spine” of the spherocylinder as the axis of length $2A_{i}$ that goes down the center of the rectangle, as indicated by the solid lines in Fig. 1. For every point on the perimeter of the spherocylinder, the shortest distance from the spine is $R_{i}$. We define the aspect ratio of the spherocylinder as,

so that $\alpha_{i}\to 0$ describes a circular particle, and the ratio of the total tip-to-tip length to width is $1+\alpha_{i}$. In this work we consider only systems in which all particles have the same aspect ratio $\alpha$.

Our system consists of $N$ such spherocylinders confined within a square box of length $L$. We use periodic boundary conditions in both the $\mathbf{\hat{x}}$ and $\mathbf{\hat{y}}$ directions. The packing fraction is,

where $\mathcal{A}_{i}$ is the area of spherocylinder $i$. Unless otherwise stated, the results in this work are for a bidisperse mixture of spherocylinders, with equal numbers of big and small particles, with $R_{b}/R_{s}=1.4$. However we have also considered a monodisperse system.

We specify the position of a spherocylinder by the location of its center of mass $\mathbf{r}_{i}=(x_{i},y_{i})$, which lies at the center of the rectangle. The orientation of the spherocylinder is given by the angle $\theta_{i}$ that the spine makes with respect to the $\mathbf{\hat{x}}$ axis, as shown in Fig. 1(a). Two spherocylinders $i$ and $j$ come into contact when the shortest distance between their spines, $r_{ij}$, is less than the sum of their radii $d_{ij}=R_{i}+R_{j}$. When $r_{ij}<d_{ij}$ the contact between the spherocylinders may be one of three types, as illustrated in Figs. 1(b), (c), and (d), respectively: (i) tip-to-side, (ii) tip-to-tip, or (iii) side-to-side. In order to have a side-to-side contact (iii) rather than a tip-to-side contact (i), in principle it is necessary that the two spherocylinders be perfectly parallel, i.e. $\theta_{i}=\theta_{j}$; in practice, due to limitations in the numerical accuracy of our contact detection algorithm Pournin et al. (2005), we take two spherocylinders as parallel whenever $|\theta_{i}-\theta_{j}|<10^{-8}$. When this happens, we take the point of contact to be midway between the corresponding endpoints of the spines of $i$ and $j$, as indicated in Fig. 1(d).

To determine when two spherocylinders are in contact, and if so to then determine the value of $r_{ij}$ and the location of the contact point, we use the efficient algorithm defined in Ref. Pournin et al. (2005). In such a case we model the elastic contact force as a simple one-sided harmonic repulsion which acts at the point of contact only when $r_{ij}<d_{ij}$. The elastic force on spherocylinder $i$ due to contact with $j$ is thus given by,

where $k_{e}$ sets the energy scale, and $\mathbf{\hat{r}}_{ij}$ is the unit normal to the surface at the point of contact, pointing inward to spherocylinder $i$. The total elastic force on spherocylinder $i$ is therefore,

where the sum is over all spherocylinders $j$ in contact with $i$. Although the elastic force always acts normal to the surface, there can nevertheless be a torque exerted on the spherocylinder due to the non-circular shape. The total elastic torque on $i$ is,

where $\mathbf{s}_{ij}$ is the moment arm from the center of mass $\mathbf{r}_{i}$ of spherocylinder $i$ to the point of contact with spherocylinder $j$, as illustrated in Fig. 1(b).

In this work we are interested in the jamming of the system when it is uniformly compressed from a dilute state. To model a uniform compression we shrink the box length at a constant rate, $dL/dt=-\kappa L$. We can then regard the substrate area of the box as undergoing a similar affine contraction, with the local velocity of the substrate at position $\mathbf{r}$ being,

The shrinking box then interacts with the spherocylinders via a viscous drag force between spherocylinder and substrate. We denote the local velocity of a position $\mathbf{r}$ on spherocylinder $i$ by,

where the first piece is due to the center of mass motion and the second piece is due to the spherocylinder’s rotation. The total dissipative force on spherocylinder $i$ is then,

where the integral is over the area of spherocylinder $i$. There is similarly a dissipative torque on the spherocylinder,

Using $\int_{i}d^{2}r\,[\mathbf{r}-\mathbf{r}_{i}]=0$ by symmetry, taking the area of spherocylinder $i$ as $\int_{i}d^{2}r=\mathcal{A}_{i}$, and defining,

we can write the dissipative force and torque as,

We take the elastic and dissipative forces as the only forces acting on the spherocylinders; there is no interparticle frictional force or collisional dissipation. Taking an overdamped equation of motion,

Eqs. (13) and (14) can then be numerically integrated to find the motion of the center of mass $\mathbf{r}_{i}(t)$ and the orientation $\theta_{i}(t)$.

For determining mechanically stable configurations very close to the jamming transition, we will also use a conjugate gradient energy minimization applied to the configurations generated by the above compression protocol. We defer discussion of this minimization procedure to Sec. III.3.

III Results

To measure $n-$fold orientational order in two dimensions, the magnitude of the order parameter $S_{n}$ and its direction of orientation $\theta_{n}$, for any particular configuration, can be computed as Donev et al. (2006),

where the $\theta_{n}$ that maximizes the sum is the ordering direction. One can then show that,

Choosing $n=2$ measures the nematic order while $n=4$ measures tetratic order.

III.2 Jamming Transition

In this section we investigate the jamming transition of a bidisperse mixture of spherocylinders as a function of the aspect ratio $\alpha$. Here, and in subsequent sections, we use a system with $N=1024$ spherocylinders. At sufficiently small packing fraction $\phi$, the spherocylinders are dilute enough that they may avoid all contact with each other and the system is at zero pressure. As $\phi$ increases, the system will ultimately become so dense that spherocylinders will necessarily come into mutual contact, force chains will percolate across the system, and a finite pressure will develop. For frictionless particles, the pressure increases continuously from zero at a specific packing fraction $\phi_{J}$, known as the jamming transition.

We perform such compression runs for a bidisperse system with $N=1024$ spherocylinders, computing the pressure of configurations at regular time intervals and averaging over the $M_{s}$ independent samples. In Fig. 3 we plot the resulting average pressure $\langle p\rangle$ vs $\phi$ for the two specific cases of (a) nearly circular disks with $\alpha=0.01$, and (b) moderately elongated spherocylinders with $\alpha=4$. We use $M_{s}=6$ to $10$, depending on the compression rate $\kappa$.

We see that $p=0$ at low $\phi$ and then $p$ increases to finite values as $\phi$ increases above some $\phi_{J}(\kappa)$. As $\kappa$ decreases, $\phi_{J}(\kappa)$ increases, the curves sharpen up near $\phi_{J}(\kappa)$, and $\langle p\rangle$ increases linearly in $\phi$ sufficiently above $\phi_{J}(\kappa)$, as expected for our harmonic elastic force O’Hern et al. (2003). For $\kappa\leq 10^{-9}$ we see no change in the $\langle p\rangle$ vs $\phi$ curve, and we have reached the limit of quasistatic compression. The value of $\phi_{J}$ in this quasistatic limit is the critical packing fraction of the compression-driven jamming transition. The small tail that is seen near $\phi_{J}$ in this quasistatic limit is a finite size effect. For finite $N$, each sample $s$ has a slightly different, sample specific, value of $\phi_{Js}$, as has been observed previously for circular disks O’Hern et al. (2003) and as we confirm for spherocylinders below; as $N\to\infty$, this spread in $\phi_{Js}$ shrinks to zero.

To estimate the value of $\phi_{J}$ for each aspect ratio $\alpha$, we consider the runs at $\kappa=10^{-10}$, which are in the quasistatic limit. We look at each of the $M_{s}$ samples separately and fit the part of the $p$ vs $\phi$ curve where the pressure first develops a linear behavior upon increasing $\phi$, before there occurs any plastic rearrangements that may lead to discontinuous drops in pressure. Extrapolating this linear region to $p=0$ then determines $\phi_{Js}$ for this particular sample. In Fig. 4 we show two examples of such determinations for the case $\alpha=4.0$. We then average over these $\phi_{Js}$ to determine $\langle\phi_{J}\rangle$. In Fig. 5 we plot the resulting $\langle\phi_{J}\rangle$ vs aspect ratio $\alpha$. At $\alpha=0$ we find $\langle\phi_{J}\rangle=0.8412\pm 0.0005$, consistent with earlier results for circular disks O’Hern et al. (2003); Vågberg et al. (2011). As $\alpha$ increases, $\langle\phi_{J}\rangle$ increases to a maximum $\langle\phi_{J}\rangle\approx 0.8875$ around $\alpha\approx 1$, and then decreases. The results we see here for $\phi_{J}(\alpha)$ are qualitatively similar to those found in previous simulations of ellipsoids and spherocylinders in 3D Donev et al. (2004a, b); Man et al. (2005); Donev et al. (2007); Sacanna et al. (2007); Williams and Philipse (2003); Wouterse et al. (2007).

III.3 Mechanically Stable Configurations and Lack of Isostaticity

To investigate the question of isostaticity at jamming in spherocylinders, we will want to consider the density of states for vibrational modes at the jamming transition. The density of states is found by expanding the energy of the system about a mechanically stable state (i.e., a local energy minimum) to second order in small displacements of the degrees of freedom, and finding the eigenstates of the resulting dynamical matrix. Since the configurations we obtain from compressing are the result of a dynamical (albeit slow) process, they are not necessarily in exact mechanical equilibrium. We therefore wish to energy minimize the configurations we obtain from compression.

To get configurations close to jamming, we will consider configurations in which the total elastic energy $U$ (defined below) is fixed to the value $U_{0}/L^{2}=10^{-15}$. We choose configurations at a fixed value of $U$, rather than a fixed value of $\phi$, since the jamming point $\phi_{Js}$ varies slightly from sample $s$ to sample $s^{\prime}$; fixing $U$, rather than $\phi$, ensures that all our samples will be about the same distance from their sample specific jamming transition. For our harmonic elastic force we have $U/L^{2}\propto(\phi-\phi_{Js})^{2}$, and we find that $U_{0}/L^{2}=10^{-15}$ corresponds to $(\phi-\phi_{Js})\lesssim 10^{-7}$.

To locate configurations with the desired $U_{0}$, we start with a configuration with $U>U_{0}$ obtained from our continuous compression runs at a fixed rate $\kappa$, and energy minimize it using a conjugate gradient algorithm (see Appendix A for details). Depending on whether the resulting minimized energy $U$ is greater or less than $U_{0}$, we carry out an affine decompression or compression of the box,

and then energy minimize the resulting configuration. We continue such decompression or compression steps until $U$ crosses the value $U_{0}$. We then reduce $\lambda$ by half, and reverse direction, i.e., if we had been decompressing, we now compress, and vice versa. We continue in this fashion until we have narrowed in on the desired value $U=U_{0}$. We start this process with a value $\lambda=10^{-6}$ and stop when $\lambda<10^{-16}$, which we find gives and accuracy in the energy of $|U-U_{0}|/U_{0}\lesssim 10^{-7}$.

To implement the above minimization procedure we must define a global energy function consistent with the elastic forces of Eq. (3). We use,

where the sum is over all pairs of contacts between spherocylinders $i$ and $j$, $r_{ij}$ is the shortest distance between their two spines, and $d_{ij}=R_{i}+R_{j}$ is the sum of their end cap radii.

Once we have found energy minimized states sufficiently close to jamming, we will then wish to construct the dynamical matrix. We find it convenient to convert the orientation angle $\theta_{i}$ into a length, and thus we take the coordinates of a given spherocylinder $i$ to be written as $\bm{\zeta}_{i}=(x_{i},y_{i},A_{i}\theta_{i})$, where $\zeta_{i1}=x_{i}$, $\zeta_{i2}=y_{i}$ and $\zeta_{i3}=A_{i}\theta_{i}$. The dynamical matrix is then the $3N\times 3N$ matrix,

where $i,j=1,2\dots,N$, and $a,b=1,2,3$, and the derivatives are evaluated at the energy minimized configuration.

To evaluate $M_{ia,jb}$ we need to know how $r_{ij}$ depends on the coordinates $\bm{\zeta}_{i}$ and $\bm{\zeta}_{j}$ of the two spherocylinders in contact, since we have,

The dependence of $r_{ij}$ on the spherocylinder coordinates depends on which of the three types of contacts of Fig. 1(b), (c), (d) that one is considering. For the tip-to-tip contact of Fig. 1(c), a small displacement of any of the two spherocylinder’s coordinates will keep the contact tip-to-tip. We can therefore write,

with the appropriate signs taken so as to minimize $|\Delta x_{ij}|$ and $|\Delta y_{ij}|$.

For a tip-to-side contact, as in Fig. 1(b), a motion of either spherocylinder parallel to the side with the contact, or a rotation of the spherocylinder with the tip contact, will result in a sliding of the location of the side contact. The calculation of $r_{ij}$ must be done more carefully. If $i$ is the spherocylinder with the side contact and $j$ is the spherocylinder with the tip contact, then,

where the sign is taken so as to minimize $r_{ij}$.

For a side-to-side contact, if we take the location of the contact bond as illustrated in Fig. 1(d), then a small rotation of either spherocylinder changes the configuration from a side-to-side contact to a tip-to-side contact, with a resulting discontinuous jump in the location of the contact point, and hence in the torques on the spherocylinders. This discontinuity makes the derivatives needed for the dynamical matrix ill-defined, and moreover also causes difficulties carrying out the conjugate gradient minimization procedure. We therefore modify the contact energy for this case as illustrated in Fig. 6. Instead of a single contact located midway between the ends of the opposing spines (dotted line in Fig. 6) we now model the side-to-side contact as two bonds located at the corresponding ends of the spines (solid lines labeled ${r}_{ij}^{(a)}$ and ${r}_{ij}^{(b)}$ in Fig. 6). We use the same convention when doing our conjugate gradient minimization of the energy $U$, provided both $r_{ij}^{(a)}$ and $r_{ij}^{(b)}$ are points of spherocylinder overlap, i.e. $r_{ij}^{(a)},r_{ij}^{(b)}<d_{ij}$. Taking spherocylinder $j$ as the one whose tip comes closest to the spine of spherocylinder $i$, then the bond where the spherocylinder overlap is larger (i.e., $r_{ij}^{(a)}$ in Fig. 6), is given by the same relation as Eq. (28). The bond where the overlap is smaller (i.e., $r_{ij}^{(b)}$ in Fig. 6) is given by,

where the sign is taken so as to maximize $r_{ij}^{(b)}$.

Modeling a side-to-side contact by two contact bonds as described above, rather than one, is also physically reasonable since in the hard-core limit a side-to-side contact will constrain two degrees of freedom: translational motion perpendicular to the spherocylinder spine, as well as rotational motion.in Ref. Wouterse et al. (2007) a similar effect was noted for plane-to-plane contacts in cut spheres, where each such planar contact constrains three degrees of freedom. However the authors of that work continued to count such planar contacts as a single contact. Their result for the number of contacts as a function of aspect ratio is therefore an underestimate of the correct constraint counting that should be done to test for isostaticity. In contrast, tip-to-side and tip-to-tip contacts constrain only one degree of freedom. Hence when counting the number of contact bonds per particle $z$, we count each side-to-side contact as two bonds. A similar observation was made in Ref. Azéma et al. (2013) for polyhedral shaped particles.

Using the above procedure, we construct mechanically stable configurations with the desired $U_{0}/L^{2}=10^{-15}$, very close to jamming. Having obtained the mechanically stable configurations, we then remove any “rattler” particles. We take a rattler to be any particle which has only zero or one contact with another particle (here, and here only, a side-to-side contact is counted as one contact). For a particle with only two contacts, we also take it to be a rattler unless the two contacts are oriented on opposite sides parallel to the spherocylinder spine; in this case the spherocylinder may still be important for the stability of the contact network, even if it has a zero-energy sliding mode in the direction parallel to the spine. Passing through the configuration to remove such rattlers, we then iterate the process until no further rattlers are found.

Having such mechanically stable configurations, obtained by energy minimization as discussed above, will be essential for our analysis of the eigenmodes of the small elastic vibrations of the system, to be discussed in the next section. However, we find in practice (checking explicitly for $\alpha=0.01,1.0$ and 4.0) that $\langle z_{J}\rangle$ changes negligibly if we compare the value computed in these energy minimized configurations with the value computed in the quasistatically compressed configurations from which we start the minimization procedure. Rather than carry out energy minimization at all values of $\alpha$, we therefore use the values of $\langle z_{J}\rangle$ found from our quasistatically compressed configurations. In Fig. 8(a) we plot the resulting $\langle z_{J}\rangle$ vs $\alpha$. In this figure the circular data points give the values of $\langle z_{J}\rangle$ when we count each side-to-side contact as two bonds, as illustrated in Fig. 6; we believe this is the correct approach to properly count the number of constraints, and all cited numerical values for $\langle z_{J}\rangle$ represent values computed in this way. We see that $\langle z_{J}\rangle$ has a peak near the same value of $\alpha\approx 1$ that gives the peak in $\langle\phi_{J}\rangle$, and that it decreases as $\alpha$ increases further. Thus, unlike 2D ellipses and 3D ellipsoids Donev et al. (2004a, b); Man et al. (2005); Donev et al. (2007); Zeravcic et al. (2009); Mailman et al. (2009); Schreck et al. (2012), spherocylinders in 2D are not approaching the isostatic limit as they get increasingly elongated. At the peak value, $\langle z_{J}\rangle\approx 5.91\pm 0.01$, close to the isostatic value of 6, but still smaller, so the system is always hypostatic.

In Fig. 9 we show the fraction of contact bonds of each of the three different types (i.e., side-to-side, tip-to-side, tip-to-tip) as a function of the aspect ratio $\alpha$. As we do in computing $z$, each side-to-side contact is counted as two bonds. Not surprisingly, the fraction of side-to-side contacts increases as $\alpha$ increases. However, consistent with our preceding arguments concerning $\langle z_{J}\rangle$, we find that the fraction of side-to-side contacts remains finite as $\alpha\to 0$. The fraction of tip-to-side contacts similarly stays finite as $\alpha\to 0$. Indeed, for $\alpha=0.01$, we find that virtually all the particles ($96.3\%$ of them) have a contact on at least one of their two flat sides, even though the flat sides represent only $0.63\%$ of the perimeter length. This is readily seen in Fig. 7(a).

To examine this propensity for spherocylinders at small $\alpha$ to have contacts on their flat sides, we measure the probability for a spherocylinder to have a contact at a particular point on its surface. Defining $\varphi$ as the polar angle that a given point on the spherocylinder surface makes with respect to the spine (see inset to Fig. 10), we measure the probability density $P(\varphi)$ to have a contact at angle $\varphi$. We average over both big and small particles. In calculating this distribution, we will count side-to-side contacts as only a single bond, located along the flat side as in Fig. 1(b), since we are more interested in the geometry of the contacts rather than counting constraints. In Fig. 10 we plot $P(\varphi)$ vs $\varphi$ for values of $\alpha=1.0,0.12$ and $0.01$. We see clearly that as $\alpha$ decreases, a sharp peak grows at $\varphi=90^{\circ}$, i.e. along the flat side. In contrast, for a circular disk this distribution would be flat. The smaller, broader, side peaks observed near $\varphi=30^{\circ}$ and $150^{\circ}$ may be interpreted as a shadow effect; if a contact exists at an angle $\varphi$, then a neighboring contact is generally no closer than $\varphi\pm 60^{\circ}$.

The prevalence of contacts along the flat sides of the spherocylinders, even as $\alpha\to 0$ and the length of these flat sides becomes a negligible fraction of the total spherocylinder surface, suggests that the presence of flat sides makes the $\alpha\to 0$ limit in some sense singular. As the length $2L$ of the spherocylinder spine shrinks to zero, the system nevertheless seems to remember what direction that spine is in. This conclusion appears to be robust, as we have demonstrated by the following check. Rather than starting our energy minimization to obtain mechanically stable states from our quaistatically compressed configurations, we start from a jammed configuration of perfectly circular disks ($\alpha=0$). We then choose a random spine direction for each particle and distort it into a spherocylinder with $\alpha=0.01$. We then follow the procedure discussed above to vary the system box size, and energy minimize, so as to obtain a new mechanically stable state of the spherocylinders at $U/L^{2}=10^{-15}$, close to jamming. The resulting $P(\varphi)$ for these configurations is found to be the same as in Fig. 10.

In response to our above observation, Vanderwerf et al. Vanderwerf et al. (2017) have recently computed the analogous $P(\varphi)$ for a bidisperse distribution of 2D elliptical particles with minor to major axis ratio $b/a$. Although the effect is not as dramatic as we find for spherocylinders, they similarly find an increasing probability for contacts along the minor axis of the ellipse, as one takes the limit $b/a\to 1$. This suggests that the effect may hold generally for barely aspherical particles, rather than be specifically due to the flat sides of the spherocylinders.

III.4 Density of States and Eigenmodes

Having obtained mechanically stable configurations and eliminated rattler particles, in this section we analyze in detail the spectrum of the eigenmodes of the $3\mathcal{N}\times 3\mathcal{N}$ dynamical matrix $M_{ia,jb}$ of Eq. (23), determining the matrix eigenvalues $\lambda_{m}$ and the corresponding normalized eigenvectors,

gives the small displacement of spherocylinder $i$ in eigenmode $m$. The frequencies of elastic vibration are then $\omega_{m}=\sqrt{\lambda_{m}}$. We define the resulting density of vibrational states $D(\omega)$ by counting the number of modes within bins of equal relative widths $\Delta\omega/\omega$. We normalize the density of states so that $\int d\omega D(\omega)=1$.

We will also wish to characterize the nature of the eigenvectors, in particular whether they correspond to localized or extended modes, and the extent to which they involve translational or rotational motion of the particles. To measure the extent to which the modes are extended, involving the correlated motion of large groups of spherocylinders, or localized, involving only a few spherocylinders, we compute the participation ratio $P_{m}$, defined as Zeravcic et al. (2009),

where $\delta x_{im}$, $\delta y_{im}$, $A_{i}\delta\theta_{im}$ give the components of the eigenvector $\mathbf{\hat{u}}_{m}$ as in Eq. (31). For an extended mode in which each degree of freedom is excited equally, we have $P_{m}=3$; for a localized mode in which only a single degree of freedom is excited, we have $P_{m}=1/\mathcal{N}$. Taking the average of $P_{m}$ over all eigenmodes with frequencies $\omega_{m}$ within bins of equal width $\Delta\omega/\omega$, we then define the participation ratio $P(\omega)$ for modes at frequency $\omega$.

To measure the extent to which a given eigenvector involves translational motion parallel to the spherocylinder’s spine, translational motion perpendicular to the spine, or rotational motion about the center of mass, we define Zeravcic et al. (2009); Mailman et al. (2009); Schreck et al. (2012) quantities $u^{2}_{\parallel m}$, $u^{2}_{\perp m}$ and $u^{2}_{\theta m}$,

Because each eigenvector is normalized to unity, we have $u^{2}_{\parallel m}+u^{2}_{\perp m}+u^{2}_{\theta m}=1$. Taking the average of these quantities over all eigenmodes with frequencies $\omega_{m}$ within bins of equal width $\Delta\omega/\omega$, we then define $u^{2}_{\parallel}(\omega)$, $u^{2}_{\perp}(\omega)$ and $u^{2}_{\theta}(\omega)$ to describe the average behavior of modes at frequency $\omega$.

In this work we will treat in detail three specific cases: aspect ratio $\alpha=0.01$, corresponding to nearly circular particles, $\alpha=4.0$, corresponding to moderately elongated particles, and $\alpha=1.0$, corresponding to the peak value of the packing fraction and also the largest value of $\langle z_{J}\rangle\approx 5.91$. We start by considering $\alpha=0.01$. In Fig. 11(a) we plot our results for the density of states $D(\omega)$ vs $\omega$ for mechanically stable configurations at energy $U/L^{2}=10^{-15}$, very close to the jamming transition. For comparison, we also show results for $\alpha=0$, i.e., perfectly circular disks; for $\alpha=0$ there is also a delta function contribution (not shown) at $\omega=0$ that represents the $\mathcal{N}$ non-interacting rotational modes of the $\mathcal{N}$ disks. Our results here, and for other quantities in this section, are averaged over six independent samples.

For $\alpha=0.01$ we find behavior qualitatively similar to that found previously for elliptical and ellipsoidal particles Zeravcic et al. (2009); Mailman et al. (2009); Schreck et al. (2012). $D(\omega)$ shows two distinct bands of frequencies, separated by a clear gap. The upper frequency band consists of the finite energy modes usually associated with disordered granular solids near jamming. Rather than the $D(\omega)\sim\omega^{2}$ behavior found in uniform elastic solids, there is a proliferation of floppy modes (the “boson peak” O’Hern et al. (2003)) as $\omega$ decreases, causing $D(\omega)$ to drop sharply as one goes to the low frequency edge of this upper frequency band. We find that the total number of modes in this upper frequency band is precisely $\mathcal{N}\langle z\rangle/2$, corresponding to the $\mathcal{N}\langle z\rangle$ contacts that serve to constrain any large-length scale motion, so that the system is jammed and can support a finite pressure.

In Fig. 11(b) we show the participation ratio $P(\omega)$. We see that modes at the edges of either the high frequency band or the low frequency band are localized, with small values of $P(\omega)$. But modes in the center of either band are fairly extended. In Fig. 11(c) we plot the quantities $u^{2}_{\theta}(\omega)$, $u^{2}_{\parallel}(\omega)$ and $u^{2}_{\perp}(\omega)$. We see that modes in the high frequency band are mixed in nature, with roughly equal participation in each of the rotational and translational degrees of freedom. For the lower frequency band, however the situation is different. Modes in the upper part of this band involve primarily rotational motions, similar to what was found for the entire low frequency band for ellipses and ellipsoids Zeravcic et al. (2009); Mailman et al. (2009); Schreck et al. (2012). However, unlike with ellipses and ellipsoids, modes towards the lower edge of this band involve only translational motion, with motion parallel to the spherocylinder spine somewhat greater than motion perpendicular to the spine. That the lowest energy modes involve translational motion is presumably a reflection of the flat surfaces that exist on the sides of the spherocylinders.

In Fig. 12 we illustrate graphically two examples of eigenmodes in the low frequency band. In these figures, arrows on each spherocylinder are proportional to the translational displacement of the spherocylinder, while the color of each spherocylinder indicates the degree of rotation: blue is a counterclockwise rotation, while red is clockwise, with the darkness of the color proportional to the amount of the rotation. Fig. 12(a) is for a mode at $\omega_{m}=3\times 10^{-3}$, somewhat in the middle of the band. As expected from Fig. 11, one clearly sees that this mode is extended throughout the system and is mixed between translational and rotational motion. In contrast, the mode in Fig. 12(b) at $\omega_{m}=3\times 10^{-2}$, near the upper edge of the band, is clearly seen to be more localized and consists primarily of rotational motion.

When displacing the spherocylinders an amount $\delta$ in the direction of a given eigenmode $\mathbf{\hat{u}}_{m}$, the energy of the system will increase by $\Delta U=\omega_{m}^{2}\delta^{2}$, to lowest order in $\delta$. It is interesting to see how $\Delta U(\delta)$ behaves as one increases $\delta$ away from the small $\delta$ limit. Note, for a highly localized mode we expect that $\delta\sim 1$ corresponds to the displacement of a particle on the order of a particle diameter $2R_{s}$; for a highly delocalized mode, we expect that $\delta\sim 1$ corresponds to the displacement of particles on the order of $2R_{s}/\sqrt{3N}\sim 2R_{s}/50$. In Fig. 13 we plot $\Delta U(\delta)$ vs $\delta$ for several typical modes $m$ at the different frequencies $\omega_{m}$ as shown. For $\omega_{m}$ in the high frequency band, we see that $\Delta U\sim\delta^{2}$ for the entire range $0\leq\delta\leq 1$. For $\omega_{m}$ in the low frequency band, we see that $\Delta U$ crosses over from a form $A\delta^{2}$ at small $\delta$ to a form $B\delta^{2}$ at larger $\delta$, with $A\ll B$. The small $A$ corresponds to the small eigenvalue $\lambda_{m}=\omega_{m}^{2}$ of the lower band; this $\lambda_{m}$ would be zero if the system were exactly at $\phi_{J}$. The cross over region to the larger $B$ suggests the quartic nature of these low band modes (i.e., $\Delta U\sim\delta^{4}$) that has been predicted Donev et al. (2007) to hold exactly at the jamming $\phi_{J}$. As $\delta$ increases further we find that the contact network starts to change significantly; the “grazing” particle overlaps (overlap $\sim\delta^{2}$) that characterize the quartic nature of the low band modes at small $\delta$ start to break, and new “ordinary” contacts start to form as particles push into each other with overlaps similar to those found in the modes of the higher band (overlap $\sim\delta$). This results in the larger value $B$, which is comparable to similar values found in the high frequency band.

Finally we explore the dependence of the eigenmodes on the energy of the system $U$. In Fig. 14(a) we plot $D(\omega)$ vs $\omega$ for mechanically stable configurations at several different values of $U$. As was found previously for ellipses and ellipsoids Zeravcic et al. (2009); Mailman et al. (2009); Schreck et al. (2012), we see that as $U$ increases, the upper frequency band changes relatively little, while the lower frequency band increases, the gap between the two bands narrows, and ultimately the two bands merge. Defining $\bar{\omega}_{0}$ as the average frequency of the modes in the lower band, in Fig. 14(b) we plot $\bar{\omega}_{0}$ vs $U/L^{2}$. We see a perfect power law dependence, strongly suggesting that the lower band of modes collapses to $\omega\to 0$ as $U/L^{2}\to 0$ at jamming. In particular we find $\bar{\omega}_{0}\sim(U/L^{2})^{1/4}$. Since, for our harmonic elastic interaction $U/L^{2}\sim(\phi-\phi_{J})^{2}$, this gives, $\bar{\omega}_{0}\sim(\phi-\phi_{J})^{1/2}$, in agreement with results found previously for ellipses and ellipsoids Schreck et al. (2012).

III.4.2 Moderately elongated spherocylinders: α=4.0𝛼4.0\alpha=4.0

In Fig. 15(b) we plot the participation ratio $P(\omega)$, and in Fig. 15(c) we plot the quantities $u_{\theta}^{2}(\omega)$, $u_{\parallel}^{2}(\omega)$ and $u_{\perp}^{2}(\omega)$. We see that the high frequency band consists of a set of mostly extended modes with mixed rotational and translational motion, similar to what was found for $\alpha=0.01$. However we see that all the modes in the low and middle bands are strongly localized. In the low band these modes are entirely translational in the direction parallel to the spine of the spherocylinder. In the middle band these modes are primarily rotational.

In Fig. 16 we illustrate graphically two examples of eigenmodes, one in the low frequency band and one in the middle frequency band. Because these modes are highly localized, we show only a subregion of the system containing the spherocylinder that moves, rather than the entire system. Fig. 16(a) is for a mode at $\omega_{m}=10^{-8}$, in the low band. It is clear that this mode consists of only a single spherocylinder that slides parallel to its spine. Fig. 16(b) is for a mode at $\omega_{m}=10^{-4}$, in the middle band. Again this mode consists of only a single spherocylinder, but now the motion is primarily rotational. In both cases, although the isolated spherocylinder may move along one degree of freedom with low cost in energy, its presence is nevertheless clearly important for the global rigidity of the system.

In Fig. 17 we show the change in energy $\Delta U(\delta)$ as the spherocylinders are displaced a distance $\delta$ in the direction of several given eigenmodes $\mathbf{\hat{u}}_{m}$. For the modes in the high frequency band, we see $\Delta U\sim\delta^{2}$ for most of the range of $\delta$. For the two lowest modes shown, at $\omega_{m}=10^{-4}$ and $10^{-8}$ in the middle and low band respectively (these are the same two modes illustrated in Fig. 16), we do not have sufficient numerical accuracy to compute $\Delta U$ at the smallest values of $\delta$. For the middle band mode, which is mostly rotational, we see behavior similar to that found for the low band modes of $\alpha=0.01$. As $\delta$ increases, $\Delta U(\delta)$ transitions from a small $\delta$ behavior of $A\delta^{2}$, with $A=\omega_{m}^{2}$, to $B\delta^{2}$ with $A\ll B$; the transition region has a quartic dependence $\sim\delta^{4}$. For the translational mode in the low band, the small $\delta$ behavior $\Delta U=\omega_{m}^{2}\delta^{2}$ is too small to be computed accurately and so does not appear in our plot. We believe this very small energy is due to the fact that the spherocylinder involved in this mode is not exactly parallel with its neighbors, and so the energy in the side-to-side contact changes ever so slightly as the spherocylinder slides parallel to its spine. The sharp step upwards seen for this mode in Fig. 17 corresponds to a displacement large enough that the tip of the sliding spherocylinder starts to contact and overlap a neighbor that it previously did not touch. We suspect that exactly at $\phi_{J}$, such sliding modes may be strictly unconstrained (for small enough $\delta$), rather than quartic.

Finally in Fig. 18(a) we show how the density of states $D(\omega)$ changes as $U$ increases and the system moves further from jamming. As $U$ increases, the high frequency band changes little. For sufficiently large $U$, the middle band merges with the upper band. Defining the average frequency of the modes in the lowest band of localized, translational, modes as $\bar{\omega}_{0}^{t}$, and the average frequency of the modes in the middle band of localized, primarily rotational, modes as $\bar{\omega}_{0}^{r}$, we plot $\bar{\omega}_{0}^{t}$ and $\bar{\omega}_{0}^{r}$ vs $U/L^{2}$ in Fig. 18(b). We see that they both vary as a power law as $U/L^{2}$ decreases, strongly indicating that the frequencies of the modes in these two bands vanish exactly at $\phi_{J}$. However we see that they vanish with different power laws. We find for the middle band modes that $\bar{\omega}_{0}^{r}\sim(U/L^{2})^{1/4}\sim(\phi-\phi_{J})^{1/2}$, the same as was found for the low band modes for nearly circular spherocylinders with $\alpha=0.01$, and the same as was found for ellipses and ellipsoids Schreck et al. (2012). However for the low band modes we find that they vanish more rapidly as $U\to 0$, with $\bar{\omega}_{0}^{t}\sim(U/L^{2})^{1/2}\sim(\phi-\phi_{J})$.

III.4.3 Spherocylinders near the peak packing fraction: α=1.0𝛼1.0\alpha=1.0

Finally in this section we consider spherocylinders with aspect ratio $\alpha=1.0$, thus being a case between those considered in the two previous sections; $\alpha=1$ also corresponds to the aspect ratio that gives the peak packing fraction $\langle\phi_{J}\rangle\approx 0.8875$ and which is also closest to being isostatic, with $\langle z_{J}\rangle=5.91\pm 0.01$.

In Fig. 19(a) we plot our results for the density of states $D(\omega)$ vs $\omega$ for mechanically stable configurations at energy $U/L^{2}=10^{-15}$, very close to the jamming transition. In Fig. 15(b) we plot the participation ratio $P(\omega)$, and in Fig. 15(c) we plot the quantities $u_{\theta}^{2}(\omega)$, $u_{\parallel}^{2}(\omega)$ and $u_{\perp}^{2}(\omega)$. We see that the situation at $\alpha=1.0$ is a natural combination of the two previous cases. There are three distinct frequency bands, with the two upper bands looking essentially the same as was found for $\alpha=0.01$. States are localized near the edges of these bands but extended in the middle. The highest frequency band consists of modes that are mostly of mixed translational and rotational character. The middle frequency band is primarily rotational towards its upper edge, but primarily translational towards its lower edge. The lowest frequency band is like that found for $\alpha=4.0$, consisting of highly localized sliding modes that are purely translational parallel to the spherocylinder spine. We may thus speculate that this represents the generic case. As $\alpha$ decreases from unity, the low frequency band of sliding modes shrinks and disappears while the middle frequency band grows. As $\alpha$ increases from unity, the middle frequency band shrinks while the low frequency band grows. For $\alpha=1.0$ we find the fraction of modes in the low, middle, and high frequency bands to be 0.0039, 0.0105, and 0.986 respectively

It is interesting to note that three frequency bands were also reported for ellipses and ellipsoids at low aspect ratios Mailman et al. (2009); Schreck et al. (2012). In that case the authors argued that it was only their lowest frequency band that represented the quadratically unconstrained states. In our case of spherocylinders at $\alpha=1.0$, however, it is both the lowest two bands that represent quadratically unconstrained states. Our high frequency band is found to contain all the $\mathcal{N}\langle z\rangle/2$ modes expected for the quadratically constrained states, while the modes in the lowest two bands scale to zero as $U\to 0$ at the jamming transition. To demonstrate this, we compute the average frequency $\bar{\omega}_{0}^{r}$ of modes in the middle band, and the average frequency $\bar{\omega}_{0}^{t}$ of modes in the lower band, as a function of the system energy $U/L^{2}$. As we found in the preceding section for $\alpha=4.0$, we similarly find here that $\bar{\omega}_{0}^{r}\sim(U/L^{2})^{1/4}$ while $\bar{\omega}_{0}^{t}\sim(U/L^{2})^{1/2}$ (see Fig. 20).

IV Conclusions

We first considered the question of orientational ordering. Rod-shaped particles in thermal equilibrium are known to have a normal-liquid to nematic-liquid phase transition as the density increases, where the orientational ordering increases from zero as one goes above the transition. In contrast, we find that, for moderately elongated spherocylinders with aspect ratio $\alpha=4.0$, there is no orientational ordering as the system is athermally compressed to above jamming. We find this to be true for both monodisperse and bidisperse ensembles. Thus the fluctuations induced by athermal collisions under compression would seem to be qualitatively different from thermally induced fluctuations.

We then considered the limit $\alpha\to 0$, in which spherocylinders are approaching the limit of perfectly circular disks. Surprisingly, we found that this limit appears to be singular, with a strong probability developing for the spherocylinders to form contacts along their flat edges, even as those flat edges shrink to a negligible fraction of the particle surface. Similar results have recently been found for 2D ellipses Vanderwerf et al. (2017), suggesting that this may be a general behavior for non-spherical particles.

Our results confirm that small distortions of particles from a perfect circular shape result in hypostatic states at jamming, however we show that behavior for very elongated particles depends in detail on the particle shape: for ellipses and ellipsoids, particles approach isostaticity at jamming as the aspect ratio increases, while for spherocylinders they remain hypostatic. We believe that this hypostatic behavior for elongated spherocylinders is a consequence of the long flat sides of the particles, which have a strong effect on the nature of the quadratically unconstrained modes at jamming. As we were completing this work we learned of similar work being carried out by Vanderwerf et al. Vanderwerf et al. (2017).

Acknowledgements

This work was supported by National Science Foundation Grant No. CBET-1435861. Computations were carried out at the Center for Integrated Research Computing at the University of Rochester. We thank S. V. Franklin, P. Olsson and C. S. O’Hern for helpful discussions.

Appendix A

In this appendix we provide details of our energy minimization procedure, and tests that show how well our procedure results in mechanically stable jammed states. Starting from an initial state obtained by our slow compression algorithm, we energy minimize to obtain a mechanically stable configuration by using the Polak-Ribiere conjugate gradient method.

With $\bm{\zeta}=(\bm{\zeta}_{1},\bm{\zeta}_{2},\dots)$ a vector giving the initial position of our configuration in the $N$ particle coordinate space, we compute the steepest descent gradient $\mathbf{v}=-\partial U/\partial\bm{\zeta}$, take as our initial search direction the unit vector $\mathbf{\hat{v}}=\mathbf{v}/|\mathbf{v}|$, and perform a line search to find the approximate minimum in this direction. We begin the line search by choosing a small step size $\varepsilon$ and finding the energies of the current configuration $U(\bm{\zeta})$ as well as $U(\bm{\zeta}+\varepsilon\mathbf{\hat{v}})$ and $U(\bm{\zeta}+2\varepsilon\mathbf{\hat{v}})$. If the energy increases when moving to $\bm{\zeta}+\varepsilon\mathbf{\hat{v}}$, i.e., $U(\bm{\zeta})<U(\bm{\zeta}+\varepsilon\mathbf{\hat{v}})$, we take $\varepsilon^{\prime}=\varepsilon/2$ and find the new set of energies with $\varepsilon^{\prime}$. If the energy decreases monotonically to $\bm{\zeta}+2\varepsilon\mathbf{\hat{v}}$, i.e. $U(\bm{\zeta})>U(\bm{\zeta}+\varepsilon\mathbf{\hat{v}})>U(\bm{\zeta}+2\varepsilon\mathbf{\hat{v}})$, then we take $\varepsilon^{\prime}=2\varepsilon$ and find the new set of energies. Once we have a series of three points with the lowest energy at $\bm{\zeta}+\varepsilon\mathbf{\hat{v}}$ we make a quadratic fit to the points, and determine the location of the minimum $\bm{\zeta}_{0}$ of that quadratic fit. If $U(\bm{\zeta}_{0})<U(\bm{\zeta}+\varepsilon\mathbf{\hat{v}})$ we then move the configuration to $\bm{\zeta}_{0}$; otherwise we move it to $\bm{\zeta}+\varepsilon\mathbf{\hat{v}}$. We then use the Polak-Ribiere method to define the new, orthogonal, search direction. For our system size and packings near jamming, we find empirically that an initial value of $\varepsilon=10^{-4}$ is a good choice. Recall, in our units, the smaller spherocylinders have a diameter of unity.

As the algorithm narrows in on a local minimum of $U$, the step size $\varepsilon$ needed to complete a line search gets ever smaller. Once $\varepsilon<10^{-16}$ we no longer have sufficient machine precision in the particle coordinates to accurately compute the energy difference $U(\bm{\zeta})-U(\bm{\zeta}+\varepsilon\mathbf{\hat{v}})$, and so we stop the line search, and reinitialize the search using the steepest descent direction at the current configuration coordinates. When the search in the steepest descent direction similarly fails to find a new minimum with $\varepsilon>10^{-16}$ we terminate the search.

Here $\mathbf{s}_{ij}$ is the moment arm from the center of mass of spherocylinder $i$ to the point of contact with spherocylinder $j$, the second sum is over all spherocylinders $j$ in contact with spherocylinder $i$, and each contact gives rise to two terms, one for the torque on spherocylinder $i$ and one for the torque on spherocylinder $j$.

Appendix B

In this appendix we discuss our method to determine the the eigenmodes and density of states $D(\omega)$ for our system of spherocylinders with aspect ratios $\alpha=1.0$ and 4.0. Because the eigenmodes in the lowest frequency band have exceedingly small frequencies $\omega_{m}$, we find that a direct analysis of the dynamical matrix results in large errors in these smallest eigenvalues. We therefore follow Refs. Wyart et al. (2005a, b); Donev et al. (2007); Schreck et al. (2012) and split the dynamical matrix into two pieces,

From Eq. (22), for our harmonic elastic interaction $\partial^{2}V_{ij}/\partial r_{ij}^{2}=k_{e}/d_{ij}^{2}$ at each contact is $O(1)$ in our units (where $k_{e}=1$ and $d_{s}=1$), while $\partial V_{ij}/\partial r_{ij}=-k_{e}(1-r_{ij}/d_{ij})/d_{ij}$ is proportional to the particle overlap and hence very small close to jamming. Hence we can regard $S_{ia,jb}$ as a small perturbation of $H_{ia,jb}$. The eigenvectors $\mathbf{\hat{u}}_{m}^{(0)}$ and eigenvalues $\lambda_{m}^{(0)}$ of the stiffness matrix $H_{ia,jb}$ are therefore the zeroth order approximates to those of $M_{ia,jb}$. We thus compute these $\mathbf{\hat{u}}_{m}^{(0)}$ and find them to include a set of degenerate eigenvectors with $\lambda_{m}^{(0)}=0$. These would be the unconstrained (to quadratic order) modes present in a hypostatic system, were the system exactly at the jamming transition where $S_{ia,jb}=0$. At finite energy $U$ above jamming, we take these as the zeroth order approximates to the modes in the lower frequency bands, while the eigenvectors with $\lambda_{m}^{(0)}>0$ are the approximates to the modes in the upper frequency band.

We can then compute the first order corrections to the eigenvalues, due to the non-zero $S_{ia,jb}$ at finite $U$, in the usual way. For a mode in the upper frequency band we have $\delta\lambda_{m}=-\mathbf{\hat{u}}_{m}\cdot\mathbf{S}\cdot\mathbf{\hat{u}}_{m}$. For such a mode in the upper frequency band we find that $\mathbf{\hat{u}}_{m}^{(0)}$ and $\lambda_{m}^{(0)}+\delta\lambda_{m}$ differ negligibly from the $\mathbf{\hat{u}}_{m}$ and $\lambda_{m}$ we obtain from a direct analysis of $M_{ia,jb}$.

For the modes in the lower frequency bands, we project $S_{ia,jb}$ onto the subspace spanned by the set of degenerate eigenvectors $\{\mathbf{\hat{u}}_{m}^{(0)}\}$ with $\lambda_{m}^{(0)}=0$, and then diagonalize $-S_{ia,jb}$ on that subspace. The resulting eigenvectors $\mathbf{\hat{u}}_{m}$ and eigenvalues $\lambda_{m}$ are then the next level approximates to the eigenmodes of the lower frequency bands of the full dynamical matrix $M_{ia,jb}$, and these values are then used in the construction of the density of states $D(\omega)$.

References