Two classes of events in sheared particulate matter

TL;DR

This paper identifies two distinct types of contact changes in sheared particulate matter, revealing that at finite pressure, not all contact changes lead to instabilities, unlike the zero-pressure case.

Contribution

It distinguishes between network and rearrangement contact events in sheared particles and analyzes their scaling behavior across different conditions.

Findings

At finite pressure, contact changes include both stable and unstable events.

The fraction of each event type remains constant regardless of system size or pressure.

Zero-pressure limit exhibits a highly singular nonlinear response.

Abstract

Under shear, a system of particles changes its contact network and becomes unstable as it transitions between mechanically stable states. For hard spheres at zero pressure, contact breaking events necessarily generate an instability, but this is not the case at finite pressure, where we identify two types of contact changes: network events that do not correspond to instabilities and rearrangement events that do. The relative fraction of such events is constant as a function of system size, pressure and interaction potential, consistent with our observation that both nonlinearities obey the same finite-size scaling. Thus, the zero-pressure limit of the nonlinear response is highly singular.

Figures

b). We refer to these are *rearrangements*. Second, we observe events where individual contacts are broken or created but the stress remains smooth (Fig. 1c). We refer to these as *network events*.

a inset). We can convert this to a relative fraction of rearrangements: $\eta=0.135\pm 0.015$ for Hertzian 2D disks, $\eta=0.145\pm 0.010$ for Hookean 2D disks, and $\eta=0.143\pm 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.

Equations

$$\Delta Z=\frac{1}{N}W(N^{2(\xi -1)}p)=\frac{1}{N}W(N^{\beta }p)$$

$$⟨f_{min}⟩\sim \frac{p^2}{N},$$

Full text

††preprint: APS/123-QED

Under shear, systems of interacting particles undergo a sequence of transitions between mechanically stable states. For supercooled liquids and glassy solids, localized excitations that allow the material to flow and fail under an applied strain correspond to low-energy saddle points of this landscape Brito and Wyart (2009); Manning and Liu (2011); Dasgupta et al. (2012). While this scenario is well-established for materials with attractive interactions or at significant finite pressures, it is not clear whether this scenario holds for pressures close to zero.

The physics of weakly compressed systems of repulsive soft particles is governed by the jamming transition. Linear response properties, such as the shear and bulk moduli, robustly follow finite size scaling as a function of the number of particles $N$ and pressure $p$ Goodrich et al. (2012); Dagois-Bohy et al. (2012); Corwin (2012); Goodrich et al. (2014, 2016). These and other results for the linear response suggest that jamming is a bona fide phase transition. The number of contacts between particles plays a crucial role: at zero pressure, systems of repulsive spheres are isostatic, so that the number of degrees of freedom equals the number of constraints. The mechanical response of such systems is highly anomalous, and the removal of any contact will destabilize the system Maxwell (1864); Calladine (1978); Lubensky et al. (2015). At finite pressure, the excess number of contacts above isostaticity $N\Delta Z$ also is governed by finite size scaling, revealing a small pressure regime with only one excess contact  Goodrich et al. (2012); Dagois-Bohy et al. (2012); Goodrich et al. (2014, 2016).

Here we investigate jammed systems under large shear deformations, focusing on contact changing events van Deen et al. (2014, 2016). Near the critical point, breaking a single contact might have a dramatic effect. In particular, at zero pressure, any contact breaking corresponds to a loss of stability, while at finite pressures, where the excess number of contacts is positive, it could be possible to break a contact without initiating an instability. It is however not yet clear how often this happens, or even whether it is possible to distinguish between true instabilities that correspond to a saddle point, and innocuous contact changes not associated with a saddle point.

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 both the shear stress and irreversibiliy, while network events are not. We study the frequency of network events and rearrangements as a function of $N$ and $p$ and find that the fraction of all contact changes that correspond to instabilities is remarkably constant as a function of $N$, contact potential, and $p$. This highlights that while all contact breaking events are instabilities at isostaticity, at all pressures above isostaticity the two are abruptly no longer one-to-one.

Furthermore, we find that near isostaticity the nonlinear properties, such as stress jumps, degree of irreversibility associated with rearrangements, and the mean strain steps between subsequent network events and subsequent rearrangements, exhibit the same finite-size scaling that has previously been reported for linear response Goodrich et al. (2016). Specifically, the scaling exponents extracted from simulations appear to be integer, dependent on the interaction potential (Hertzian or Hookean), and can be rationalized using simple arguments. This surprising correspondence between the linear and nonlinear scaling is similar to what is found in Random Field Ising Models Kent-Dobias and Sethna (2018).

Finally, Wyart has recently argued that if finite contact potentials at small pressure introduce small perturbations to the hard-sphere system, then one expects the minimum force $f_{min}$ across a Hertzian contact at the onset of instability to scale as $p^{3}$ 111Unpublished work by Matthieu Wyart. 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 instead scales as $f_{min}\sim p^{2}$.

Taken together, our results highlight that the zero pressure limit, which corresponds to the physics of hard spheres, is highly singular and cannot by itself describe the nonlinear finite-pressure behavior of particulate matter in low dimensions. In other words, the phenomenology of jamming at small but finite pressures is different from that at zero pressure.

Model.— We simulate systems of 2D bidisperse disks in a 50-50 mixture with size ratio 1:1.4 in a square box under athermal quasistatic shear Maloney and Lemaître (2006). We perform an infinite temperature quench to obtain an initial packing and use pairwise particle contact potentials to define the energy: $U=\sum_{ij}\Theta(\varepsilon_{ij})\varepsilon_{ij}^{\xi}$ where $\varepsilon_{ij}$ is the overlap between particles i and j, $\Theta$ is the heaviside function, and $\xi$ represents the power of the potential. The stress $\sigma$ is calculated via the Born-Huang approximation Born and Huang (1954). We study Hertzian ($\xi=5/2$) and Hookean ($\xi=2$) disks in 2D and 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 $\gamma$, we utilize standard athermal quasi-static shear with Lees-Edwards boundary conditions where the energy is minimized after each strain step via the FIRE algorithm Bitzek et al. (2006) 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 van Deen et al. (2014) with a resolution of $\Delta\gamma<10^{-13}$ 222There is a finite density of particles with less than $D+1$ contacts, which are termed rattlers. In our statistics, we only consider particles which are not rattlers either before and after an event.. As we frequently need to distinguish small numbers from those that are zero to numerical precision, we employ quad precision calculations and perform the minimizations on GPUs 333pyCudaPack available upon request at https://github.com/simonsglass/pycudapacking/ to access a broad range of parameters with high numerical precision.

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. For each case, we apply shear until there have been either 30 contact changes in 2D or 60 contact changes in 3D, which allows us to stay within the pre-yielding regime. Our first main 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. 1b). We refer to these are rearrangements. Second, we observe events where individual contacts are broken or created but the stress remains smooth (Fig. 1c). 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 $\Delta\sigma$ associated with a contact change is defined using a linear extrapolation of the shear modulus $G$: $\Delta\sigma=(\gamma_{+}-\gamma_{-})G_{-}-(\sigma_{+}-\sigma_{-})$, where the indices $-$ and $+$ represents the last point before the contact change and the first point after the contact change. 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 inter-particle forces before and after an infinitesimal strain loop: \Delta F(\gamma_{n})=\big{(}\Sigma_{k}|\textbf{f}^{(k)}_{\rightarrow}-\textbf{f}^{(k)}_{\leftarrow}|^{2}\big{)}^{1/2}, 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 $\Delta F$ vs $|\Delta\sigma|$ reveal that network events and rearrangements are clearly distinguished (Fig. 1d,e): Network events exhibit no stress drops to within numerical precision $\Delta\sigma<10^{-14}$, while stress drops associated with rearrangements are finite 444While we cannot rule out the existence of saddle points that are not associated with contact changes a priori, we have found no evidence of such events, and our data places an upper bound on the stress drops that are associated with such events of $N\Delta\sigma<10^{-10}$..

Our data unambiguously shows 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 $|\Delta\sigma|$ exhibits power-law scaling Salerno et al. (2012), implying that arbitrarily small stress drops are possible, more recent work  Franz and Spigler (2017) 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 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 will necessarily correspond to a saddle point or instability Wyart (2012), which is clearly not the case at any finite pressure.

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^{\beta}p$ Goodrich et al. (2012); Dagois-Bohy et al. (2012), implying that jamming has the finite size scaling properties of a real phase transition Widom (1965); Fisher and Barber (1972); Goodrich et al. (2014, 2016); Sastry (2016).

Therefore, we use our data to study the statistics of stress drops $\Delta\sigma$ and irreversibility in the force network $\Delta 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 $\Delta Z$ to generic contact potentials O’Hern et al. (2002, 2003):

[TABLE]

where $\beta=3$ for Hertzian systems and $\beta=2$ for Hookean systems. Systems with $N^{\beta}p<1$ have on average only one contact above isostaticity, and exhibit a different behavior than those with $N^{\beta}p>1$. Figures 2a and 2b illustrate that the scaling collapse of the nonlinear response – including both stress drops $\Delta\sigma$ and force network irreversibility $\Delta F$ across saddle points – is consistent with the previously reported scaling of linear response near jamming Goodrich et al. (2012).

In the region where $N^{\beta}p>1$, we find that the scaling of stress drop and the irreversibility metric before rearrangements scale with $N$ and $p$ as $\langle\Delta\sigma\rangle\sim\frac{p}{N}$ and $\langle\Delta F\rangle\sim p$ respectively. The linear scaling of both with pressure implies that the only force scale in the system is the one set by the pressure. And while the typical size of a 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.

Characteristics Strain Scales.— In order to measure the relative frequency of network events and rearrangements, we measure the average strain between network events $\Delta\gamma_{net}$ and the average strain between rearrangements $\Delta\gamma_{rear}$ and analyze the scaling of both as a function of $N^{\beta}p$ (Fig. 3(a,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 Maloney and Lemaître (2006) predicted the exponent for strain between rearrangements in steady-state flow by arguing that all accumulated energy must be dissipated, so $\Delta\gamma_{rear}\sim\frac{\Delta\sigma}{G}$. Given the numerical observation that the shear modulus scales with the pressure and interaction potential as $G\sim p^{(\xi-3/2)/(\xi-1)}$ O’Hern et al. (2003), we obtain $\Delta\gamma_{rear}\sim\frac{p^{1/2(\xi-1)}}{N}=\frac{p^{1/\beta}}{N}$.

Surprisingly, Figure 3 shows this applies even in the pre-yielding regime before steady-state flow, and not only to the distance between rearrangements $\Delta\gamma_{rear}$, but also to the distance between network events $\Delta\gamma_{net}$. This is in agreement with the results of previous work by van Deen van Deen et al. (2016), which we extend in the supplement (Supplemental section B.)

A simple regression analysis between $\Delta\gamma_{rear}$ and $\Delta\gamma_{net}$ for Hertzian 2D disks yields $\Delta\gamma_{rear}/\Delta\gamma_{net}=5.88\pm 0.35$ (Figure 3a inset). We can convert this to a relative fraction of rearrangements: $\eta=0.135\pm 0.015$ for Hertzian 2D disks, $\eta=0.145\pm 0.010$ for Hookean 2D disks, and $\eta=0.143\pm 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.

Wyart recently predicted the low-pressure scaling of the minimum force across rearranging contacts, by assuming that small overlaps at low pressures generate a small perturbation to the energy landscape of hard spheres at zero pressure Note (1).

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, which is consistent with our data shown in Supplemental Figure S4.

For Hertzian spheres in 2D and 3D we measure the minimum force before an event, and find that

[TABLE]

as shown in Fig 4. Our data are not consistent with the prediction from a small-overlap expansion around hard spheres, which for a Hertzian interaction potential is $f_{min}\sim p^{3}$. These results also suggests the hard sphere limit is singular.

*Discussion.—*We have identified two types of contact changes that exist in finite pressure systems: network events and rearrangements.

First, the statistics of each are governed by finite size scaling relations with integer exponents, corresponding to those of linear response quantities. It is far from obvious that the linear and nonlinear response should obey the same scaling ansatz, although similar relationships have been observed in other systems Kent-Dobias and Sethna (2018). In jammed particulate matter, it could be that the low-frequency localized excitations that contribute to the linear response in low dimensions also govern the transition states, as postulated previously Karmakar et al. (2010); Manning and Liu (2011); Gartner and Lerner (2016). Second, the relative proportion of rearrangements is small (about 14 %), and independent of interaction potential, system size, and pressure. Third, the minimum force at the onset of instability scales differently than predicted based on an small-overlap expansion around the hard-sphere zero-pressure state Note (1).

Together, these results highlight that the zero-pressure isostatic limit of purely repulsive particles is highly singular. At zero pressure, every contact breaking event should correspond to an instability, but Fig. 3 demonstrates that at any finite pressure and system size the fraction of contact changes that are instabilities is a fixed constant far from unity. This may come from the fact that in soft spheres, the limit of zero pressure always corresponds to one excess contact, whereas in hard spheres, there are zero excess contacts.

Looking forward, it will be interesting to study whether the vibrational spectrum of soft spheres at low pressure are different in important ways from those of hard spheres at zero pressure. The eigenvalue spectrum seems to be well-described by an expansion around zero pressure – Wyart et al Wyart et al. (2005) and Goodrich et al Goodrich et al. (2012) demonstrate that isostatic systems at zero pressure must possess soft, system-spanning modes associated with breaking a single contact, and that these extended modes give rise to a boson peak in the density of states that shifts to higher frequencies as the pressure increases away from zero.

The eigenvector spectrum remains more mysterious. A recent manuscript suggests that extended modes found at zero pressure also destabilize the system at small but finite pressures, and that instabilities continuously become more localized as the pressure increases Shimada et al. (2018). However, other work hints that there might be different types of instabilities in 2D and 3D systems. Charbonneau and collaborators Charbonneau et al. (2015) studied the statistics of small forces in jammed hard-sphere packings as a function of dimension, and found that they scale as a power law predicted by mean-field Wyart (2012); Lerner et al. (2013); DeGiuli et al. (2014) and infinite-dimensional Charbonneau et al. (2014) theories, provided that localized buckler excitations were removed. The existence of these localized excitations, which are much more prevalent in 2D and 3D, suggest that mean-field extrapolations from zero pressure may not fully explain the mechanical response. It would be interesting to study whether the finite-size scaling of $f_{min}$ we observe, which is not consistent with zero-pressure mean field predictions, might instead be consistent with a different, localized class of excitations that govern instabilities in low dimensions.