Jamming and irreversibility

We investigate irreversibility in soft frictionless disk packings on approach to the unjamming transition. Using simulations of shear reversal tests, we study the relationship between plastic work and irreversible rearrangements of the contact network. Infinitesimal strains are reversible, while any finite strain generates plastic work and contact changes in a sufficiently large packing. The number of irreversible contact changes grows with strain, and the stress–strain curve displays a crossover from linear to increasingly nonlinear response when the fraction of irreversible contact changes approaches unity.


Introduction
Among many other topics in the physics of granular matter, Bob Behringer's research has had outsized impact on the field of jamming [1][2][3]. His measurements of the jump and subsequent power law growth in the contact number above a critical packing fraction [4] represents the first, and still one of the few [5][6][7], measurements of a major hallmark of isotropic jamming. His work on shear jamming [8], dilatancy [9], and contact force statistics [10,11] upended the conventional view of the jamming phase diagram and illuminated how granular materials' rigidity encodes their loading history. Here, inspired by Bob's work, we ask how shearing can wipe out the memory of an initially isotropic state in a weakly jammed solid. In other words, how does irreversibility emerge near jamming?
Packings of soft spheres prepared at small but finite pressure are marginal solids-while their response to infinitesimal strains is elastic [1], a small shear stress suffices to instigate quasistatic plastic flow [12,13]. Recently there has been considerable interest in how the ensemble-averaged stress-strain curve for shear becomes nonlinear, and in particular on how the crossover from linear to nonlinear response depends on the distance to jamming [14][15][16][17][18][19][20][21][22][23]. The shear strain required to make or break a contact vanishes in the limit of large system sizes, so finite deformations necessarily involve topological changes to the contact network [24][25][26][27][28]. It is therefore natural to ask about the relationship between nonlinearity and plasticity, especially when one approaches (un)jamming. More precisely, we ask whether there is a correlation between the linear-to-nonlinear crossover and (ir)reversibile contact changes.
To probe nonlinearity and irreversibility near jamming, we study shear reversal in marginally jammed packings of athermal, frictionless, purely repulsive soft spheres. We begin from an isotropic state prepared at a targeted pressure p. We use this initial pressure (prior to shearing) to quantify the distance to unjamming at p = 0 . After preparation, the system is subjected to simple shear in small quasistatic steps to a maximum strain m . The shearing direction is then reversed, and the system is returned to zero strain. A load is reversible if the stress follows the loading curve back to its initial value at zero strain. Reversible and irreversible deformations are illustrated in Fig. 1 with data from our simulations. This complements similar irreversibility under volumetric strain as observed in [21] and interpreted in terms of a history-dependent critical packing fraction.
The present work builds on results from Boschan et al. [19,29], who studied the loading curve but did not consider shear reversal. The loading curve was found to be linear up to a strain scale † ∼ p . After † the stress continues to grow, albeit more slowly than an extrapolation of the initial linear trend. The crossover to steady plastic flow occurs later, at a distinct strain scale y ≃ 0.05 . Simulations of large amplitude oscillatory shear at finite rate also showed two distinct crossovers with identical scaling properties [22].
Boschan et al. [19] also studied contact changes, i.e. made and broken contacts during shearing. They found that the linear-to-nonlinear crossover at † is also evident in the contact change statistics, as detailed in Sect. 4. It is plausible that contact changes are a proxy for irreversible rearrangements, but this must be verified-while rearrangements involve contact changes, not all rearrangements are irreversible [22,[30][31][32][33][34][35].
Here we probe nonlinearity and and irreversibility during a loading-unloading cycle. We first monitor the plastic work performed during the cycle, and then correlate these results to the statistics of contact changes at the particle scale. We find, first, that there is finite plastic work even when the ensemble-averaged stress-strain curve is linear. Consistent with this observation, we also find that irreversible contact changes accrue prior to the loss of linearity. Second, prior to † , some fraction of the contact changes are reversible. After † , when the stress-strain curve is nonlinear, essentially all contact changes are plastic.

Model and methods
We perform two-dimensional simulations of athermal frictionless disks, a standard model with a jamming transition [2]. Particles experience a spring-like force proportional to their overlap ij = (R i + R j ) − r ij , where R i and R j denote the radii and r ij is the length of the vector ij pointing from the center of particle i to j. The contact force on particle i due to particle j is purely repulsive, and there is no interaction when the particles are not in contact, where a hat indicates a unit vector. We fix the units of stress by setting the spring constant k and mean particle size to unity. The stress tensor is where Greek indices denote Cartesian coordinates, and V is the total area of the unit cell.
Initial conditions are created by randomly populating the bi-periodic simulation box and then using a nonlinear conjugate gradient energy minimization protocol to quench instantaneously to a local minimum of the elastic potential energy at fixed volume [36]. The box is then allowed to undergo small changes in size and shape to achieve a target pressure p and zero shear stress-these are called "shear-stabilized" packings in the nomenclature of Dagois-Bohy et al. [37,38]. Packings are bidisperse to avoid crystallization; we use the standard [1,36] 50:50 mixture of small and large particles and a radius ratio of 1:1.4.
Once the initial state is prepared, we apply quasistatic simple shear using Lees-Edwards boundary conditions with small logarithmically-spaced steps ranging between = 10 −8 … 10 −3 . After each strain step the energy is reminimized [36] while holding the strain fixed, so particles follow quasistatic trajectories. Once a maximum strain m is reached, the direction of shear is reversed and the system is returned to zero strain, again via a series of small logarithmically-spaced steps.
In order to quantify irreversibility, we calculate the plastic work W p of the loading/unloading cycle, where upwards and downwards pointing arrows are used to indicate the loading and unloading curves, respectively. Clearly W p is zero when the response is reversible.
The phenomenology of a shear reversal test in weakly jammed soft spheres is illustrated in Fig. 1. In panel (a), the maximum shear strain m = 10 −5 is so small that no contact . 1 Sample output from a loading-unloading cycle in simulations. a If deformation is reversible, the loading curve ( ↑ ) and unloading curve ( ↓ ) coincide. b In an irreversible deformation there is hysteresis, and the enclosed area is equal to the plastic work changes occur [26,27]. The stress-strain curve is linear and the loading and unloading curves coincide. In panel (b), the maximum shear strain m = 10 −2 is substantially larger. On reversal the stress decreases but does not retrace the loading curve. The loading and unloading curves are both nonlinear. Because there is hysteresis, there is an associated plastic work. In addition to the plastic work, irreversibility can be quantified by the plastic strain p and a plastic stress p , corresponding to the intercepts of the unloading curve with the x-and y-axis, respectively.

Plastic work
We perform shear reversal tests for a range of preparation pressures p and varying maximum strain m . Figure 2 illustrates loading and unloading curves for p = 10 −4 and m ranging from 10 −5 to 10 −2 in half-decade steps. The result is representative of other pressures.
To quantify the appearance of irreversibility, we analyze the plastic work as a function of m and p, as shown in Fig. 3. We find nonzero W p for all investigated maximum strains, which are as small as 10 −5 . (As noted above, packings of finite size can be sheared reversibly if the contact network remains unchanged, but this strain interval vanishes in the large system size limit [26,27]). For each pressure W p has an approximately power law growth with m , with an apparent exponent that varies with pressure.
To better understand the pressure dependence of W p , we seek to collapse the data to a master curve. Anticipating a correlation with the onset of nonlinearity, we plot the rescaled variable x ≡ m ∕p ∼ m ∕ † . On the other axis we plot the rescaled work W ≡ W p ∕p for some exponent . To motivate , we note that for small values of m , the loading curve is associated with work W ↑ ∼ G 0 2 m , where G 0 ∼ p 1∕2 is the shear modulus for Hookean particles near jamming [1,39,40]. If we assume G 0 also sets the relevant scale for W p at small m , then we expect W p ∼ p 1∕2 2 m . Rearranging in favor of m ∕p gives W p ∕p 5∕2 ∼ m ∕p 2 , which requires = 5∕2 . This prediction is tested in the log-log plot of Fig. 3b, where we find data collapse to a curve with an initial slope of 2. When x ≫ x * ∼ O(1) the plastic work grows more slowly with m , with an exponent of roughly 3/2, Plasticity is indeed sensitive to † , because data for W p collapse with the rescaled variable x. But irreversibility does not "turn on" when the ensemble-averaged stress-strain curve becomes nonlinear, as indicated by measurable W p even when the curve is linear.

Contact changes
We now seek to relate irreversibility to the evolution of contact changes during loading and unloading.
As first shown in Ref. [19] and verified below, the scale † is apparent in the evolution of the number of made and broken contacts per particle, which we refer to as the contact change density n cc ( ) . We now monitor contact changes during unloading to see to what extent the original contact network is recovered (i.e. broken contacts are re-made and made contacts are re-broken). Contact changes are always identified with respect to the initial condition, even during unloading. The "plastic contact change density" n p cc , equal to n cc at the end of the unloading curve, is a measure of irreversible (i.e. plastic) contact changes. Figure 4 depicts loading and unloading curves for three values of m and three different initial pressures p = 10 −5 , 10 −4 and 10 −3 . For the lowest m , in panel (a), most contact changes are recovered at the end of the cycle and n cc has a nonzero slope. n p cc is nevertheless nonzero, and it increases as p tends to zero. Plastic contact changes also increase with increasing m (panels (b) and (c)). In the final panel a large fraction of the contact changes are unrecoverable, n cc hits the vertical axis with zero slope, and n p cc is nearly equal to n cc ( m ). Figure 5 depicts n cc during loading. The figure shows that data for different pressures can be collapsed to a master curve by plotting N ≡ n cc ∕( † ) 1∕2 ∼ n cc ∕p 1∕2 as a function of y = ∕p ∼ ∕ † . This collapse was first demonstrated in Ref. [19]; for completeness we present it in Fig. 5 using data from the present study. We find

Contact changes during loading
The crossover y * ∼ O(1) is compatible with x * from the plastic work. For later reference, we note that when m > † = y * p . We estimate a m ≈ 3.7 ± 0.1 by fitting Eq. (6) to N for y > 10.
We note that, by definition, n cc changes by an amount 1 / N when the system has undergone a strain cc sufficient to produce one contact change. Hence and the average strain interval between contact changes, cc can be read off from the slope of the curves in Fig. 4. (Alternatively, the probability of a contact change in the interval [ , + d ) is 1∕ cc ). In particular, when the loading curve is linear, there is a typical strain interval cc ∼ p 1∕2 ∕N between contact changes. Van Deen et al. [26,27] reached compatible results by directly resolving contact changes. As noted above, cc vanishes in the large system size limit.

Contact changes after reversal
To quantify to what extent the initial contact network can be recovered under reversal, we now monitor the plastic contact change density n p cc . Clearly n p cc = 0 if the initial contact network is fully recovered. Figure 6 plots n p cc as a function of m for three pressures and system sizes N = 128 , 512, and 1024. We find n p cc is an increasing function of m , and for a given m it is larger at smaller pressures. There is also dependence on N.
The system size-dependence in n p cc suggests that the contact change strain cc ∼ p 1∕2 ∕N plays a dominant role, as opposed to † ∼ p . To test this hypothesis, we attempt to collapse data to a master curve by plotting as a function of z ≡ m N∕p 1∕2 ∼ m ∕ cc . We find collapse plotting P ≡ n p cc N 1∕2 ∕p 1∕4 versus p 1∕2 ∕N , as shown in Fig. 6b. The master curve is The crossover value z * ∼ O(10 2 ) . Therefore after the system has undergone on the order of one hundred contact changes. The constant a p ≈ 3.5 ± 0.1.

Relating nonlinearity and irreversibility
We can use the above observations to interpret the strain scale † in terms of irreversibility. To this end, it is useful to introduce the "plastic fraction" where n m cc is the value of n cc at the end of loading. f p quantifies the extent to which marginal contact changes tend to be plastic. If f p ( m ) = 0 , then all marginal contact changes in an infinitesimal interval around m are reversible. If f p = 1 , all contact changes are plastic.
While a direct numerical evaluation of f p is noisy, we can infer its scaling properties by noting that

From contact changes to the stress-strain curve
A remaining challenge is to determine how plastic events impact stress buildup. Here we make a first attempt. We expect irreversible contact changes to have an associated stress drop p ∕N due to an eigenvalue of the Hessian matrix going to zero [41,42]. Then we assume that the infinitesimal stress d generated by a strain d has both an elastic contribution and an offsetting stress release due to irreversible events Using Eq. (9) and rewriting in dimensionless form gives It remains to determine the typical stress drop amplitude, p . The scaling relation p ∼ p suggests itself purely on dimensional grounds. Assuming this form then predicts that the right hand side of Eq. (15) depends on and p only via their ratio ∕p . Reassuringly, this is consistent with the linear-to-nonlinear crossover at † ∼ p , and with recent measurements of the secant modulus during shear startup 58 Page 6 of 7 [19] and the storage modulus in oscillatory shear [17,22]. We conclude that the typical stress drop is indeed linear in p. Eq. (15) can then be integrated to find This stress-strain curve is compatible with W p in Fig. 3b, including the 3∕2 scaling beyond ≈ † . The approach presented above is semi-empirical. A more fundamental motivation would require directly identifying plastic events to determine their frequency and associated stress drops. The necessary theoretical tools were recently developed in Refs. [42][43][44].

Discussion
We have investigated irreversibility at the macro and micro scale in systems near jamming, evidencing irreversibility in both the plastic work and the contact change statistics for small shear strains. Initially the average loading curve is linear and most contact changes are reversible. Increasing the maximum strain increases the number and fraction of plastic contact changes. For > † , the loading curve becomes nonlinear and nearly all contact changes are plastic. The onset of nonlinearity therefore corresponds not to the onset of irreversibility (as commonly assumed in continuum elasto-plastic theories), but to "fully developed" irreversibility, as reflected in the saturation of the plastic fraction f p . This crossover occurs earlier for smaller † ∼ p.
With hindsight, the above scenario is apparent in the contact change statistics. For small m , as in Fig. 4a, the plastic contact change density is much smaller than n m cc , and the unloading branch of the n cc curve ends with a nonzero slope-indicating that shearing the system "a little bit further" to ↓ < 0 would bring the system closer to its initial contact topology, i.e. fewer net contact changes. In contrast, for large m , as in Fig. 4c, n p cc is nearly equal to n m cc , and the unloading curve is flat-the system has effectively lost all memory of its initial condition.
Our work has correlated the onset of nonlinearity at the macro scale to a particle scale crossover from reversible to irreversible contact changes. Both of these crossovers are sensitive to the proximity to jamming. We have also suggested a phenomenological approach to relate irreversible rearrangements to the form of the loading curve, highlighting the need for a deeper understanding of the statistics of stress drops during loading.