Differences in plasticity between hard and soft spheres

Peter Morse ,1 Sven Wijtmans ,1 Merlijn van Deen ,2 Martin van Hecke,2,3 and M. Lisa Manning1,4 1Department of Physics, Syracuse University, Syracuse, New York 13244, USA 2Huygens-Kamerlingh Onnes Lab, Universiteit Leiden, P.O. Box 9504, NL-2300 RA Leiden, The Netherlands 3AMOLF, Science Park 104, 1098 XG Amsterdam, The Netherlands 4BioInspired Institute, Syracuse University, Syracuse, New York 13244, USA

Disordered solids such as granular materials, metallic glasses, colloids, and foams undergo a sequence of transitions between mechanically stable states when slowly sheared. Despite intense research effort, a complete picture and quantitative description of the phase space structures and mechanical properties and instabilities that govern these transitions is still lacking [1]. Two popular models to study this phenomenology are hard spheres and soft spheres, both in the absence of friction and at zero temperature [2][3][4][5][6][7][8][9][10][11][12][13][14][15]. In hard spheres, geometry completely dictates the transitions between stable states that occur whenever two spheres break a contact, and scaling predictions have been developed based on this framework [5,9]. Equilibrium states of soft spheres at zero pressure are equivalent to hard spheres,-in both cases the particles are undeformed. Therefore, it is tempting to speculate that nonlinear transitions between mechanically stable states in soft sphere systems are qualitatively similar to those in hard spheres, and hard spheres are often used as a limiting case of soft spheres [16][17][18][19].
However, soft spheres at finite pressure differ in important ways from soft spheres at zero pressure, and equivalently, truly hard spheres. At zero pressure, a system of repulsive spheres is isostatic-the number of degrees of freedom equals the number of constraints-while at any finite pressure, one excess contact above isostaticity emerges [6,7]. This subtle difference underlies the finite size scaling of the linear response of jammed soft spheres: the shear and bulk moduli exhibit finite size scaling as a function of the number of particles N and pressure p, with a plateau at small system sizes and low pressures where the system has one excess contact (see Supplemental Material Section B [44]), and power-law scaling at larger system sizes and pressures [2,4,6,7,10,13,20].
The consequences of this subtle difference between zero and small pressure for the nonlinear response is not clear. If the transitions between mechanically stable states of sheared soft spheres at finite pressure were similar to those of hard spheres (i.e., contact breaking leading to instabilities), the finite pressure case could be obtained by perturbing away from the hard sphere limit. This assumption informs prior work [16][17][18][19], and leads to concrete, testable scaling predictions for forces and pressures [21], (see Supplemental Material Section D). However, the zero pressure limit of soft spheres appears singular in linear response quantities. Hence, the crucial open question is how the nonlinear response of soft sphere packings evolves as a function of system size and pressure. Is there a smooth transition between the zero-pressure limit where all contact breaking events correspond to saddle points and higher pressures, or is the zero-pressure limit singular?
To address this question, we investigate jammed systems under finite shear deformations, focusing on contact changing events [11,14]. We find that at finite pressures, we can sharply distinguish two classes of contact changes: network events, which do not correspond to instabilities, and rearrangements, which correspond to saddles in the energy landscape. Rearrangements are associated with finite jumps in the shear stress, while network events are not. Moreover, we find that the fraction of all contact changes that correspond to saddle points is about 14% regardless of system size N, contact potential, dimension, and pressure p: at all finite pressures, contact changes are abruptly no longer one-to-one with saddle points, meaning that shear of hard spheres and soft spheres at finite pressure are qualitatively different.
Next we study the finite-size scaling of nonlinear observables associated with saddle points, including the magnitude of the stress drops, degree of irreversibility, and strain between saddle points, and find that they all exhibit the same finite-size scaling that has previously been reported for the linear response [13]. Specifically, the scaling exponents extracted from simulations appear to be integer, dependent on the interaction potential (Hertzian or Hookean), are independent of dimension, and can be rationalized using simple arguments. This surprising correspondence between the linear and nonlinear scaling is not yet understood, though it is similar to what is found in random field Ising models [22]. We also quantify the strain between network events, and find that it also exhibits the same finite-size scaling, suggesting that simple contact changes are also encoded in the potential energy landscape. Finally, with the correct finite-size scaling, we are able to resolve the statistics of the minimum force at very low pressures and find that it scales as f min ∼ p 2 , instead of as f min ∼ p 3 as extrapolated from hard spheres (see Supplemental Material Section D).
Taken together, our results highlight that the zero-pressure limit, corresponding to the physics of hard spheres, is highly singular. Moreover, the potential energy landscape of soft spheres is organized in a remarkable way to ensure that nonlinear contact changes, nonlinear saddle points, and the linear response all have the same finite-size scaling.

I. MODEL
We simulate systems of two-dimensional (2D) bidisperse disks in a 50-50 mixture with size ratio 1:1.4 in a square box under athermal quasistatic shear [3]. We perform an infinite temperature quench to obtain an initial packing and use pairwise particle contact potentials to define the energy: is the heaviside function, and ξ represents the power of the potential. The stress σ is calculated via the Born-Huang approximation [23]. We study Hertzian (ξ = 5/2) and Hookean (ξ = 2) disks in 2D and monodisperse Hertzian spheres in 3D, to understand whether the two different types of potentials, which generate distinct interparticle stiffnesses as a function of pressure, affect our results.
To apply a shear strain γ , we utilize standard athermal quasistatic shear with Lees-Edwards boundary conditions where the energy is minimized after each strain step via the FIRE algorithm [24] until the maximum unbalanced force on any particle is less than 10 −18 . As we apply shear, the system undergoes contact changes, and we use a bisection algorithm to bracket contact changes [11] with a resolution of γ < 10 −13 [25]. In order to distinguish small numbers from those that are zero to numerical precision, we employ quad precision calculations and perform the minimizations on GPUs [26] to access a broad range of parameters with high numerical precision.

II. TWO CLASSES OF EVENTS
We study contact changes in ensembles of 50 different initializations at pressures ranging from 10 −7 to 10 −2 and systems sizes from 32 to 4096. We shear the system using standard techniques (see Supplemental Material Section A) and quantify the material response. Here we focus on the preyielding regime where the average stress increases with increasing strain, stopping simulations before the postyielding regime where the average stress remains constant, as the nonlinear response can be quite different in those two regimes [27].
Our first result is that all contact changes are of two mutually distinct types (Fig. 1). First, we observe events where the stresses exhibits a finite jump and the particles undergo discontinuous motion [ Fig. 1(b)]. We refer to these as rearrangements. Second, we observe events where individual contacts are broken or created but the stress remains smooth [ Fig. 1(c)]. We refer to these as network events.
To investigate the nature and statistics of these events, we focus on infinitesimal strain loops associated with the smallest bisection interval around a contact change event. First, the stress drop σ associated with a contact change is defined using a linear extrapolation of the shear modulus where the indices − and + represent the last point before the contact change and the first point after the contact change respectively. Second, to investigate the (ir)reversibility of each contact change, we study the change in the force contact network, defined using the quadrature sum of the differences in interparticle forces before and after an infinitesimal strain loop: where k is a sum over all contact forces and the right (left) arrow indicates the force at the beginning (end) of the loop across the contact change event. Scatter plots of F vs | σ | reveal that network events and rearrangements are clearly distinguished [Figs. 1(d) and 1(e): Network events exhibit no stress drop to within numerical precision σ < 10 −14 , while stress drops associated with rearrangements are finite [28]. Our data unambiguously show that network events are perfectly reversible, while rearrangements are irreversible. Hence the two classes of events can be separated by stress drop and (ir)reversibility. Analogous results are true for 3D Hertzian spheres (Fig. S2).
We note that while previous work suggested that | σ | exhibits power-law scaling [8], implying that arbitrarily small stress drops are possible, more recent work [15] suggests a finite cutoff to the smallest possible stress drops. Our simulations are consistent with the latter scenario, and the existence of a finite minimum magnitude of stress drops associated with saddles (Supplemental Material Fig. S1) makes it possible to distinguish network events from rearrangements.
Together, our results indicate a major difference between the physics at finite albeit small pressures, and the hard sphere scenario at strictly zero pressure. For hard spheres at zero pressure, contact breaking events necessarily correspond to a saddle point or instability [5], which is clearly not the case at finite pressure.

III. FINITE SIZE SCALING
Previously, Goodrich et al. showed the finite-size scaling of linear elastic response-the shear and bulk moduli-collapses as a function of N β p [6,7], implying that jamming has the finite size scaling properties of a real phase transition [10,13,[29][30][31]. Therefore, we use our data to study the statistics of stress drops σ and irreversibility in the force network F across saddle points at very small pressures, and test the scaling ansatz developed initially for linear response.
We first extend the previous scaling ansatz for Z, the number of contacts above isostaticity plus one [6,7,32,33] (see Supplemental Material Sections B and C), to generic contact potentials [2,34]: where W is an as yet unknown scaling function and β = 3 in Hertzian systems and β = 2 in Hookean systems. Systems with N β p < 1 have on average only one contact above isostaticity [6], and exhibit different behavior than those with N β p > 1. Figures 2(a) and 2(b) illustrate that the scaling collapse of the nonlinear response-including both stress drops σ and force network irreversibility F across saddle pointsis consistent with the previously reported scaling of linear response near jamming [6]. In the N β p > 1 region, we find that the stress drop and the irreversibility metric across rearrangements scale as σ ∼ p N and F ∼ p, respectively, where the averages represent the geometric mean, as each quantity is log-normally distributed. The linear scaling of both with pressure implies that pressure sets the only force scale in the system. And while the typical stress drop goes to zero with increasing system size, the size of the irreversibility metric remains constant, implying that rearrangements involve a characteristic number of particles in the pre-yielding regime. This is different from recent observations in Lennard-Jones systems, where σ exhibits a very slight system size dependence [27], giving an avenue for further study.

IV. CHARACTERISTIC STRAIN SCALES
. In order to measure the relative frequency of network events and rearrangements, we measure the average strain between network events γ net and the average strain between rearrangements γ rear and analyze the scaling of both as a function of N β p (Figs. 3(a) and 3(b)]. We again find excellent scaling collapse. Moreover, the fraction of contact changes corresponding to rearrangements remains constant at all pressures and system sizes.
Maloney and Lemaitre [3] predicted that accumulated strain between rearrangements in steady-state flow is γ rear ∼ σ G by arguing that all accumulated energy must be dissipated. Given the numerical observation that the shear modulus scales with the pressure and interaction potential as G ∼ p (ξ −3/2)/(ξ −1) [2], we obtain γ rear ∼ p 1/2(ξ −1) Fig. 3 shows that this also applies in the preyielding regime before steady-state flow, and not only to the distance between rearrangements γ rear , but also the distance between network events γ net , in agreement with the results of previous work by van Deen [14] (see Supplemental Material Section C.) A simple regression analysis between γ rear and γ net for Hertzian 2D disks yields γ rear / γ net = 5.88 ± 0.35 [ Fig. 3(a) inset]. We can convert this to a relative fraction of rearrangements: η = 0.135 ± 0.015 for Hertzian 2D disks, η = 0.145 ± 0.010 for Hookean 2D disks, and η = 0.143 ± 0.015 for Hertzian 3D spheres. These fractions are statistically indistinguishable, implying that the fraction of rearrangements does not depend on the system size, pressure, or interaction potential, and does not seem to vary between 2D and 3D. As this fraction corresponds to the ratio of saddle and contact change boundaries crossed by a worldline in phase space, it points to a new, perhaps universal, statistical measure of the configuration space of generic jammed particle packings.
In addition, measurements of forces across contacts at low pressure can also provide information about whether the zeropressure limit is singular. Wyart has recently argued that if finite contact potentials at small pressure introduce small perturbations to the hard-sphere system, then the minimum force f min across a 2D Hertzian contact at the onset of instability scales as p 3 [21]. We outline this argument in Supplemental Section D.
We note that as the forces across contacts in Hookean packings do not depend on the overlap, we expect f min not to correlate with pressure in this case, consistent with our data shown in Supplemental Material Figure S4. For Hertzian spheres in 2D and 3D we measure the minimum force before an event, and find that as shown in Fig 4, which is not consistent with the smalloverlap expansion around hard spheres, again suggesting that the hard sphere limit is singular.

V. DISCUSSION
We have identified two types of contact changes that exist in finite pressure systems: network events and rearrangements. The statistics of each are governed by finite size scaling relations with integer exponents, corresponding to those of linear response, though it is far from obvious that the linear and nonlinear response should obey the same scaling ansatz. In jammed particulate matter, it could be that the lowfrequency localized excitations that contribute to the linear response in low dimensions also govern the transition states, as postulated previously [35][36][37]. In addition, the relative proportion of rearrangements is small and independent of interaction potential, system size, and pressure.
Together, these results highlight that the zero-pressure isostatic limit of purely repulsive particles is highly singular. We speculate that this difference between hard and soft spheres may be driven by the fact that soft sphere packings have at least one excess contact above isostaticity, and these contact networks may generically have different geometries from hard sphere packing at zero pressure where there are zero excess contacts.
Importantly, the universality we observe in the potential energy landscape may be a property of the generic, highly disordered minima found using a standard technique of quenching from infinite temperature. Other initialization techniques that generate more stable or ultra stable packings result in far fewer rearrangements in the linear response regime [38], and likely change the relative proportion of rearrangements to network events. Similarly, we focus here on the response of as-quenched systems with no shear stabilization, as minimizing the box shape as a degree of freedom enforces that the first contact change is a network event [11], and therefore biases the location of the minimum and alters the relative frequency of each class of contact change. Investigating the nonlinear response of these specialized locations in configuration space is an interesting avenue for future work.
Looking forward, it will also be interesting to study whether the vibrational spectrum of soft spheres at low pressure is different from those of hard spheres. While some work suggests that features of both eigenvalues and eigenvectors can be predicted from expansions around the zero pressure limit [39,40], it remains unclear if these theories describe all the localized excitations in the spectra [12]. Furthermore, recent advances suggest that this could be measured in experiments [41][42][43].