Translational and rotational velocities in shear-driven jamming of ellipsoidal particles

We study shear-driven jamming of ellipsoidal particles at zero temperature with a focus on the microscopic dynamics. We find that a change from spherical particles to ellipsoids with aspect ratio $\alpha=1.02$ gives dramatic changes of the microscopic dynamics with much lower translational velocities and a new role for the rotations. Whereas the velocity difference at contacts---and thereby the dissipation---in collections of spheres is dominated by the translational velocities and reduced by the rotations, the same quantity is in collections of ellipsoids instead totally dominated by the rotational velocities. By also examining the effect of different aspect ratios we find that the examined quantities show either a peak or a change in slope at $\alpha\approx1.2$, thus giving evidence for a crossover between different regions of low and high aspect ratio.

Background Dense collections of circular disks in two dimensions (2D) and spheres in 3D with contact interaction at zero temperature have been studied extensively during the last decades with the aim to understand the jamming transition. Since ordinary molecular dynamics that automatically explores phase space doesn't work at zero temperature one has mainly used two other computational methods. The first is what we call isotropic jamming, which is the generation of particle packings with energy-minimization methods in various ways [1][2][3][4][5][6] and the second is to model a shear flow by steadily shearing the system, either quasistatically [7][8][9] or at different constant shear rates [10][11][12][13][14][15][16].
Since real particles are seldom perfectly spherical an interesting generalization of this model with obvious experimental relevance is the change from spherical to aspherical particles. It has already been argued that the spherical limit is singular [17][18][19][20] and the purpose of the present Letter is to explore further consequences of this change to aspherical particles.
A key result from the study of jamming of spherical particles is that jamming occurs when the average number of contacts per particle is just enough to stabilize the degrees of freedom, according to the isostatic condition, z = 2d f . Since rotations become relevant for aspherical particles we find d f = 6 for general ellipsoids and d f = 5 for ellipsoids with a symmetry axis (spheroids), which is the kind of particles studied in the present work. If the isostatic conjecture were always valid one would expect a jump in z already for a slight deviation from the spherical limit. It has however long been known that there is no jump in z [21] and that the systems are thus hypostatic at jamming, i.e. have a smaller number of contacts than number of degrees of freedom.
To see how a hypostatic packing can be jammed it is helpful to consider a contact between two ellipses that barely touch at their respective waists. This contact will then not only hinder the motion of the ellipses toward each other but will also hinder their respective ro-tations, thus affecting two different degrees of freedom for each particle. Whereas the interaction energy will be quadratic in the position coordinate it will instead be quartic-to the fourth power-in the angle coordinates. To see this, consider an ellipse, x 2 + a 2 y 2 = r 2 0 with a < 1 and thus elongated along the y axis, as illustrated in the Supplemental Material (SM). The distance from the origin to the perimeter of the ellipse is then r 2 = r 2 0 +(1−a 2 )y 2 , which for small y, such that the angle is θ = y/r 0 , becomes r(θ)/r 0 − 1 ∼ θ 2 . With two ellipses barely touching at their respective waists, θ = 0, the effect of a rotation is an overlap ∼ θ 2 , and the fourth power follows since the energy is in turn quadratic in the overlap. For a jammed collection of particles this mechanism leads to a set of quartic vibrational modes [22][23][24], beside the ordinary modes for which the energy is quadratic in the displacements. (For a particle with several contacts, not at their waists, r(θ) for each contact also has a term linear in the rotation θ. In force-balanced jammed packings we expect these linear terms to cancel each other out, such that only the θ 2 -dependence remains.) From studies of static packings it has been found that the number of quartic modes exactly matches the deviation in contact number from the isostatic value [23].
In the present Letter we do shearing simulations of both spherical and ellipsoidal particles to examine how the asphericity affects the microscopic dynamics. We do find dramatic effects. Focusing first on aspect ratio α = 1.02 we find that the ellipsoidal particles have considerably lower translational velocities and that their rotational motion gets a different role. Since the quartic modes found in static packings are primarily rotational in character [23,24] we believe that these effects are consequences of the quartic modes on the shear-driven dynamics. Our work thus opens up a new window for the study of shear-driven elongated particles, which until now have mostly been analyzed with focus on possible ordering of the particles [17,20,25] and macroscopic quantities [19,26,27]. When doing simulations also with a set of larger aspect ratios we find several quantities to have features at α ≈ 1.2, which suggests the existence of different regions of high and low aspect ratio.
Model Our system is a bidisperse collection of particles with nominal diameters d b = 1.4d s , in equal proportion. The particles are prolate spheroids, i.e. ellipsoids with two equal minor axes with d (1) > d (2) = d (3) . The aspect ratio is α = d (1) /d (2) > 1. To make the particle volume independent of α we take the semi-axes to be given by d (1) = α 2/3 d and d (2) , d s for big and small particles. The method used to check for particle overlaps is described in Ref. [17]. For overlapping ellipsoids we define a scale factor δ < 1 such that the ellipsoids are just barely touching when they are rescaled by δ, keeping the center of mass position fixed. The elastic force between overlapping particles is then wheren ij is a unit vector pointing inwards to particle i at the point of contact with particle j.
Simulations We take a purely collisional dynamics where dissipation takes place at the particle contacts. With r C i,j ≡ r C j,i the position of the point of contact of particle i with particle j-a point on the rescaled ellipsoidsthe velocity of the surface of particle i at this point is where v i is the center of mass translational motion and v is the velocity at the point of contact due to ω i , the rotational velocity of particle i. The velocity difference at that point is v C ij = v i (r C i,j ) − v j (r C i,j ) and the dissipative force is given by The equations of motion are where I i is the moment of inertia tensor. This model is one of the simplest with a reasonable dynamics but it is unusual in that it has no coupling between the tangential dissipation and the strength of the elastic force. As shown in the SM the introduction of such a coupling however only gives minor changes, at the densities of interest, simply because the dissipating forces are there anyway typically much smaller than the elastic ones. We take k e = 1 and k d = 1/2 and simulate with N = 1024 particles. The density (filling factor) is . We take the unit of length to be equal to d s , the unit of energy to be k e and the unit of mass to be m s , and let the particle mass be proportional to the particle volume. The unit of time is Asphericity and shear rate dependence Before turning to the microscopic behavior a few words on the shear Measures of translational and rotational particle velocities for spheres, α = 1.00 and ellipsoids with α = 1.02. Panels (a) and (b) show v/ω vs φ for spheres and ellipsoids. (Legends for the differentγ are in panel (d).) For spheres this quantity only depends slowly on φ but for ellipsoids it decreases in a dramatic way, and even more so for the lower shear rates. The vertical dashed lines are the respective φJ . Panel (c) shows an extrapolation of v/ω at φJ (1.02) to thė γ → 0 limit, giving [v/ω]0 ≈ 0.11. Panel (d) shows W which is the ratio of the root-mean-square rotational velocities along the minor axes and the major axis, as described in the main text. The figure shows that the rotations are predominantly around the minor axes which are the ones that can help the particles fit in with their neighbors. rate dependence of slightly aspherical particles is in order. At high shear rates such particles behave mostly as spherical particles but as the shear rate is lowered the behavior crosses over to a low shear rate behavior that is controlled by the aspericity [28]. The reason for this is that typical particle overlaps at high shear rates may be bigger than the typical deviation of the particle surface from the spherical whereas they become much smaller at low shear rates. The particles then have time to rotate to pack more densely and the behavior is then controlled by the asphericity. Since the jamming density is expected to increase linearly with the distance from the spherical point, φ J (α) − φ J (1) ∼ |α − 1|, as found from isotropic jamming [21], one expects the main effect of the asphericity close to jamming to be to shift the curves to (somewhat) higher φ [19]. Fig. S3 in the SM illustrates this crossover for the shear viscosity.
Translational and rotational velocities A study of the microscopic dynamics however shows that the onset of particle asphericity also gives other effects than shifting the jamming transition. A simple and direct way to investigate the dynamics is to determine the ratio of the nonaffine translational velocity and the rotational velocity, v/ω. v and ω are here the root-mean-square non-affine translational velocity v na i = v i −γy ix and the root-meansquare rotational velocity show the ratio v/ω vs φ for spherical and ellipsoidal particles at several different shear rates. For spherical particles, panel (a), this ratio is to a good approximation independent ofγ and takes values from 0.4 through 0.9 in the density window of the figure. For ellipsoidal particles, panel (b), the behavior is similar to that of the spheres at low densities but becomes very different at higher densities, with a dramatic shear rate dependence and a big drop. The ratio v/ω changes by a factor of about five for our range of shear rates. As this is for an asphericity ≡ α − 1 of only 2%, the effect is surprisingly big. As shown in the SM the decrease in this ratio is largely an effect of a smaller v but there is also a contribution from an increase in ω.
To study the dependence on α we compare data at the respective jamming densities, φ J (α). For α = 1.00 we take φ J (1.00) ≈ 0.648 [29] and for α > 1 the φ J (α) are determined as the densities that give similar algebraic decays of p(γ). Details are given in the SM. For the slightly aspherical particles this gives φ J (1.02) = 0.652. We stress that these are approximate values since φ J is always difficult to determine, and this is even more the case for α close to the spherical limit. These (small) uncertainties in φ J (α) are however not of any big concern since the studied quantites depend only slowly on φ.
With the pronounced decrease of v/ω with decreasingγ in Fig. 1(b) it becomes interesting to try to extrapolate to the limit of vanishing shear rate. At φ J (1.02) = 0.652-which is also the density where v/ω has its minimum-we fit to constant plus algebraic be- Fig. 1(c) the fit is good and suggests that the limiting value is a finite constant, [v/ω] 0 ≈ 0.11. This quantity in also shown in Fig. 3 Rotations around different axes Slowly sheared spherical particles close to jamming display erratic translational displacements which are needed to fit into a constantly changing environment. The lower translational velocity of the ellipsoidal particles suggests a picture where this translational motion is largely replaced by rotations. Since only the rotations around the minor axes can help the particles fit in with their neighbors, one expects these rotations to be more prominent than the rotations around the major (symmetry) axis. With the rotation vector in the body frame given by ω = ( ω (1) , ω (2) , ω (3) ), where the major axis is direction "1", Fig. 1(d) shows W ≡ [ ω (2) rms + ω Contributions of the different terms of Eq. (4) to (v C ij ) 2 -and thus to the dissipation-for both spheres (open circles,γ = 10 −5 , no significantγ-dependence) and ellipsoids for many different shear rates,γ = 1 × 10 −8 , 2 × 10 −8 , 5 × 10 −8 ,. . . , 10 −4 . The figure shows the difference between spheres and ellipsoids in two ways: (i) For spheres the velocity difference is mainly due to the translational velocity difference and the effect of the rotations-panels (b) and (c)-is to reduce the velocity difference at contacts. (ii) For ellipsoids at densities around φJ (1.02) = 0.652 and in the limit of small shear rates (crosses forγ = 10 −8 ) the velocity difference is mainly due to the rotational velocity. The translational velocity-panels (a) and (b)-gives a very small contribution.
At low densities this quantity is close to unity, indicating that the rotations are equally strong around the different semi-axes, but it increases to higher values around the jamming density. This increase is small for highγ, reaches more than a factor of two for our lowestγ, and an extrapolation at φ J , again by fitting to constant plus algebraic behavior, suggests theγ → 0 limit W 0 = 2.25.
Velocity difference at particle contacts Whereas the quantities discussed above are single-particle quantities, we now turn to the velocity difference at particle contacts, and, more specifically, the role of rotations and translations for these velocity differences. Since the velocity difference is related to the dissipation, and thereby to the shear viscosity, this is a quantity with a more direct physical relevance than the single particle velocities, and it is also the quantity that most clearly shows the difference between spherical and aspherical particles.
For the analysis we use the notation in and after Eq. (1) and write the contact velocity difference as the sum of translational and rotational parts, v C ij = v ij + v R ij . The contact velocity difference squared then becomes The relative contributions of these three terms are shown in Fig. 2 both for spheres (where the behavior is independent ofγ) and for ellipsoids at several differentγ. We then compare the behavior of ellipsoids for the lowestγ with the behavior of spheres. For spheres v C ij is dominated by the translational velocities, v ij , and reduced by the rotations, v R ij , which is seen by panels (b) and (c) together giving a negative contribution. The relative contributions of the three terms in Eq. (4) at φ = 0.648 are 1.57, −1.13, and 0.56.
For the ellipsoids the main contribution is from the rotations, v R ij , and the contribution from v ij is very small. This is seen by the crosses, which are data for α = 1.02 andγ = 10 −8 , in the three panels at φ J being close to zero and unity, respectively. After extrapolating the three terms of Eq. (4) at φ J (1.02) = 0.652 toγ → 0 their relative contributions are found to be 0.04, −0.02, and 0.98, which shows that it is the rotations that totally dominate the velocity difference. The translations only contribute 4% of the dissipated power whereas the rotations, when including the negative mixed term, contribute the remaining 96%.
Different aspect ratios We now turn to the dependence on the aspect ratio, α, by comparing behaviors of several quantities at the different φ J (α), extraploated (as shown in the SM) to theγ → 0 limit. Fig. 3(a) shows that W 0 (α)-the relative rotation around the minor axesfirst increases with increasing α, reaches a maximum at α ≈ 1.2, then starts decreasing and is close to unity at α = 2.5. It is interesting to note that this peak in W 0 (α) not coincides with the maximum jamming density, which is at α ≈ 1.5 both in determinations from isotropically jammed packings [21,22] and in the present study. In panel (b) [v/ω] 0 shows a change from one linear region to another at α ≈ 1.2. In panel (c) which shows the relative contribution of the different terms in Eq. (4)excluding α = 1 which is altogether different-the same kind of change is also present, though not that sharp. That figure shows that the contribution to the dissipation from the rotations which is very big at α = 1.02, decreases monotonously with increasing α. Panel (d) shows that the rotational velocity around the z axis (of the lab frame) associated with the shearing, which is −[ω z /γ] 0 = 1/2 for spheres, first stays almost constant but then starts decreasing for α > 1.2.
Summary and discussion To summarize, we have shown that a change from spherical to slightly ellipsoidal particles with aspect ratio α = 1.02, gives an altogether different microscopic dynamics. Comparing spheres and ellipsoids close to jamming it is found that the translational velocity is reduced by 80%, and that the rotations get a different role for the ellipsoids, contributing to the particles ability to fit together. The dramatic difference  is also seen when separating the dissipation into contributions from the translational and the rotational velocities. For spheres the dissipation is dominated by the translational motion whereas the dissipation in the ellipsoids is instead dominated by the rotations. Several quantities show different behaviors above and below α ≈ 1.2, suggesting the existence of two different regions for ellipsoids with low and high aspect ratio.
What can now be concluded about the relation between this different dynamics and the quartic modes? To settle that question beyond possible doubt would require a study where both the quartic modes and the microscopic dynamics were determined simulataneously, which is well beyond the scope of the present work. When trying to find an explanation for the present findings from properties of isotopically jammed systems, it does however seem that the presence of quartic modes is a natural candidate simply because they are the outstanding difference between packings of ellipsoids and packings of spheres. What gives additional credibility to such a conclusion is the observation that this different dynamics very strongly affects the rotations, which is in agreement with the quartic modes being primarily rotational in character [23,24]. We therefore conclude that there is strong evidence for considering the very different dynamics in shear-driven ellipsoids to be effects of the quartic vibrational modes on the shear-driven systems.
We thank S. Teitel for many discussions. The simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at HPC2N. Supplemental Material to "Translational and rotational velocities in shear-driven jamming of ellipsoidal particles"

Yann-Edwin Keta and Peter Olsson
(Dated: June 11, 2020) Note that figures and equations in this Supplemental Material are prefixed with an "S", e.g. " Fig. S1". " Fig. 1", on the other hand, is a figure in the Letter.

I. GENERAL CONSIDERATIONS
A. Ellipses and the quartic modes Figure S1 shows a closeup of two ellipses, x 2 + a 2 y 2 = r 2 0 , with a < 1, touching at their waists. The dashed curve is from a circle with radius r 0 and we here detail the argument in the main text that the relative distance between the circle and the ellipsoid depends quadratically on the angle θ, r(θ)/r 0 − 1 ∼ θ 2 .
To show this in more detail we note that x 2 = r 2 0 −a 2 y 2 in r 2 (θ) = x 2 + y 2 becomes For a small θ such that y = r 0 θ this becomes tangential part of the dissipation at a given contact may jump discontinuously to zero when the particles lose contact. With a friction parameter µ that controls the maximum size of the dissipation, F dis t ≤ µF el , our model may be described as having µ = ∞ since any contact, however weak, will be able to sustain an arbitrarily big dissipative force, given by the velocity difference at the point of contact.
Since the model with µ = ∞ can be argued to be unphysical we have performed additional simulations with finite µ. It was then found necessary to take a value as low as µ = 0.1 to get any clearly visible changes. Fig. S2 is therefore a comparison of v/ω obtained with µ = 0.1 and the data in the paper, which are for µ = ∞. Even with such a small µ there are no significant differences at the density of interest in our studies, which is around φ = 0.65. The only notable difference is a dip in v/ω at φ = 0.62, which is well outside of our region of interest. Looking into the reason for this weak dependence on µ one finds that the tangential dissipative forces, are typically much smaller than the elastic forces, at the densities of interest, which means that taking µ = 0.1 only affects a small fraction of contacts with the weakest elastic forces. We therefore conclude that the results and the analyses of the present paper would not be significantly affected by instead using a model that couples the elastic and the dissipative forces. Figure S3 is shear viscosity for spheres and slightly aspherical particles, α = 1.02, for different shear strain rates,γ. The figure shows that these two cases behave essentially the same for highγ. This is also reasonable since the typical particle overlaps at contacts are then of about the same size as typical deviations of the ellipsoidal particles from the spherical shape, which means that the effect of the asphericity should be very limited. At loweṙ γ the particle overlaps are smaller and the behavior is also quite different. The shear viscosity of the ellipsoids at lowerγ is somewhat lower than for spheres, which is consistent with the jamming transition taking place at a higher density. Fig. S3(b) which displays this data vs φ/φ J (α), shows that a rescaling of the density can accomodate this change. Note however that this similarity is not the full story. The main message of the Letter is that there are big differences in the microscopic dynamics in spite of the small differences in certain macroscopic quantities. to be close to φJ (α).) This has also been found to be the case for spherocylinders in 2D [2].

II. TRANSLATIONAL AND ROTATIONAL VELOCITIES
To look into the origin of the dip in v/ω in Fig. 1(b) we now examine v and ω for spheres and ellipsoids with α = 1.02 at the densities φ J (1.00) = 0.648, and φ J (1.02) = 0.652. Fig. S4(a) and (b) show v/γ and ω/γ and we first note that these quantities-open circles in shows that the non-affine translational velocity at low shear strain rates is much lower for ellipsoids than for spheres. The rotational velocity, panel (b), is in contrast somewhat higher for ellipsoids than for spheres. These behaviors of v and ω together give a dramatic decrease in v/ω, shown in Figs. 1(b) and (c).  S5. Velocity difference at the point of contact. This is a quantity which is related to the dissipation and thereby to the shear viscosity. Panel (a) shows the behavior of spheres and ellipsoids at their respective φJ . Panel (b) which is data at the same φ for both cases is included to show that the ellipsoids, at a given density, have lower contact velocities than the spheres, and thereby also a lower shear viscosity. panels (a) and (b)-behave the same for α = 1.00, which is consistent with v/ω being essentially independent ofγ. In comparison, for α = 1.02-open squares in panels (a) and (b)-v/γ increases more slowly than ω/γ, which explains the decrease of v/ω. The figure shows clearly that it is the lower translational velocity which is the main reason for the small v/ω of the ellipsoids. The velocities in Figure S4 are single-particle quantities. To complement this, Fig. S5(a) shows the rootmean-square velocity difference at contact, v C ij , which can be argued to be a physically more relevant quantity than the single-particle velocities due to the relation to shear viscosity from power balance, N k d (z/2)(v C ij ) 2 = σγV [3]. Panel (a) shows these quantities at their respective φ J (α). v C ij , and thus also the shear viscosity, is found to be slightly bigger for ellipsoids than for spheres. (The ij to the velocity difference at contacts. It is found that vij and v R ij are almost perfectly anti-correlated at low densities, φ < 0.58, which means that the rotations almost entirely compensate for the translational velocity differences. For the spherical particles vij and v R ij remain strongly anti-correlated but for the aspherical particles at the lowest shear rates, these correlations almost vanish and the two contributions are almost independent. contact number z, which enters the relation to the shear viscosity, is also different for the two cases, but these differences are too small to significantly affect the comparison.) Fig. S5(b) is the same quantity at the same density, φ = 0.648, for comparison, which shows that the shear viscosity at the same density and smallγ, is indeed lower for the ellipsoids, consistent with expectations since φ J (1.02) > φ J (1.00). Fig. 2 shows an analysis of the dissipation at contacting particles which was done by splitting up the velocity difference on translational and rotational velocity. We here make use of the same data to extract information in a different way.

III. TRANSLATIONAL AND ROTATIONAL CONTRIBUTIONS TO THE DISSIPATION
The quantity in focus is the correlation coefficient, which is obtained from the second term in Eq. (4), but with a different normalization, Here v ij and v R ij are the root-mean-square values. Fig. S6 shows that both spheres and ellipsoids have a region with almost perfect anticorrelation, g RT ≈ −0.996, for φ < 0.58. (Note the different density range; φ in Fig. 2 only extends down to φ = 0.60.) In this region there is Bagnold scaling, p ∼γ 2 , and the particles have on average less than two contacts, which means that rotations can almost altogether compensate for the translational velocity differences. The behavior then jumps to a region with Newtonian behavior, p ∼γ and a larger number of contacts, z > 4. It is then no longer possible for the particles to eliminate the velocity differences by rotations, and the anti-correlation is weaker, down to g RT ≈ −0.6 around the jamming density. Just as in the other quantities discussed above, there is essentially nȯ γ-dependence for the spherical particles, which are shown forγ = 10 −5 , only.
In Ref. [1], Fig. 21, this correlation was explored in a 2D model by means of a scatter plot which suggested an essentially perfect anticorrelation at low densities (panel (a)) changing to a weaker anti-correlation at higher densities (panel (b)). (In the notation of Ref. [1] V ij,T is the tangential component of the translational velocity difference, which is almost the same as the full translational velocity difference, here denoted by v ij .) As jamming is approached the behavior of ellipsoids again differs considerably from the behavior of spheres; the correlation coefficient is dropping dramatically. An extrapolation to theγ → 0 limit by fitting g RT at φ J to an algebraic behavior plus a constant, gives −g RT 0 = 0.045 which implies that v ij and v R ij are almost entirely decoupled from each other, in sharp contrast to the essentially perfect anticorrelation at low densities. As mentioned in the Letter the determination of φ J (α) is a truly difficult task and this is especially so for small asphericities. For examining the dependence on α it is however necessary to compare behaviors at their respective φ J (α) and we have therefore had to resort to some approximate determinations. Since the quantities in focus do not depend very sensitively on φ, our conclusions do not depend on the precise values of these φ J (α).
In shear driven jamming the jamming density is the density where pressure and shear stress at smallγ decay algebraically withγ. This simple recipe is however difficult to use due to both corrections to scaling, finite size effects, and-especially for the case of small asphericities-a crossover from spherical behavior at larger shear rates to the true ellipsoid behavior, as illustrated in Sec. I C above. Reliable determinations of φ J need to take all these effects into account.
Since a precise determination of φ J (α) is not possible at the present time we chose to neglect all these complications and instead determine approximate values of φ J (α) as the densities where p(γ; α) behave similarly to p(γ; 1.00). This does not imply the assumption that the models are in the same universality class, but only that the exponent q in p ∼γ q for ellipsoids is not altogether different from the value for spheres. circles-and for ten different α > 1. We note that each p(γ, α) is algebraic and to a decent approximation behaves the same as p(γ, 1). Panel (b) is φ J (α).
From Fig. S7(a) one can extract the dependence of p on α at the respective φ J (α) and constantγ and it then appears that p first increases, then reaches a maximum at α = 1.2 or 1.3 and then eventually decreases as α increases further. This is thus in agreement with the conclusions of Fig. 3. It should however be mentioned that this is only a tentative result since pressure is a quantity that does depend strongly on φ, which means that we cannot exclude the possibility of a differently looking curve if the data had instead been taken at the true φ J (α).
The constant shearing of the system in the x-y plane leads to a rotational velocity around the z axis which for spheres is given by ω z = −γ/2. The determination of ω z is however difficult since the shear-driven rotation is a small signal compared to the erratic particle rotations. To handle this we have here taken data from some additional runs with N = 65536 particles that give better statistics. These runs have however only been done for a few different α, and that is the reason why there are only five curves in Fig. S8(c).