Softening and Yielding of Soft Glassy Materials

Solids deform and fluids flow, but soft glassy materials, such as emulsions, foams, suspensions, and pastes, exhibit an intricate mix of solid and liquid-like behavior. While much progress has been made to understand their elastic (small strain) and flow (infinite strain) properties, such understanding is lacking for the softening and yielding phenomena that connect these asymptotic regimes. Here we present a comprehensive framework for softening and yielding of soft glassy materials, based on extensive numerical simulations of oscillatory rheological tests, and show that two distinct scenarios unfold depending on the material's packing density. For dense systems, there is a single, pressure-independent strain where the elastic modulus drops and the particle motion becomes diffusive. In contrast, for weakly jammed systems, a two-step process arises: at an intermediate softening strain, the elastic and loss moduli both drop down and then reach a new plateau value, whereas the particle motion becomes diffusive at the distinctly larger yield strain. We show that softening is associated with an extensive number of microscopic contact changes leading to a non-analytic rheological signature. Moreover, the scaling of the softening strain with pressure suggest the existence of a novel pressure scale above which softening and yielding coincide, and we verify the existence of this crossover scale numerically. Our findings thus evidence the existence of two distinct classes of soft glassy materials -- jamming dominated and dense -- and show how these can be distinguished by their rheological fingerprint.

creased driving strains? Which aspects of these scenarios depend on material properties, and which aspects are universal? Here we use extensive numerical simulations of particle-based models of non-Brownian amorphous solids to disentangle how the preand post-yielding regimes are connected. We evidence that there are two qualitatively distinct scenarios for the strain dependent response of soft glassy materials, depending on the rigidity of the jammed state of the material at rest.
The jamming transition has in recent years been shown to play an important organizing role for the response of soft repulsive sphere packings, which are an effective model for foams, emulsions, suspensions and granular media [2][3][4][5] . First, there is overwhelming evidence that the distance to the critical jamming point, as measured by e.g. the confining pressure P, is the key control parameter governing a packing's quasistatic, linear elastic response 3,6-9 , and more recent work has extended these findings to the linear viscoelastic response at finite frequency 8,10-12 . Second, the distance to jamming has also been found to play a key role in organizing the steady state rheology of a wide range of soft materials [13][14][15][16][17] . As linear response and steady state rheology are connected to the limits of zero and infinite strain amplitude in oscillatory shear, the distance to jamming can be expected to play a crucial role at finite strains as well. Moreover, several characteristic scales that limit the range of linear behavior all scale with the distance to jamming: when the pressure is lowered towards the critical point, the material becomes dominated by nonlinear response 8,11,[18][19][20][21] , and both the yield stress [13][14][15] and the range of strains that are free of microscopic contact breaking events 8,[22][23][24][25][26] all vanish.
To fill in the gap between zero and infinite strain amplitudes and shed light on the rheological scenarios and microscopic mechanisms of yielding, we perform simulations of packings under oscillatory shear at varying strain amplitude γ 0 , pressure P, and number of particles N, focusing on the low frequency regime. Macroscopically, we identify not one but two distinct crossover strain scales. We find that linear response gives way to softening above a strain scale γ s , after which the material's response reaches a new plateau before the material yields at a scale γ y . Unlike prior observations of two-step yielding [27][28][29] , this scenario does not require interparticle attraction. Near jamming softening and yielding scales differ in their pressure dependence: while γ y is essentially constant, the softening strain γ s vanishes linearly with P. Far above jamming the scales merge, suggesting the existence of a characteristic pressure that distinguishes jamming-like and dense systems. The relation between these rheological phenomena and the microscopic particle motion is complex. Whereas yielding is associated with the onset of diffusive particle motion, the relation between softening and microscopic rearrangements is far more subtle. While rearrangements destroy strict reversibility in the particles' trajectories, bulk properties such as the storage and loss moduli remain linear long after the first rearrangements -leading to a regime where particles exhibit irregular, trapped motion but the bulk response appears linear. Bulk softening only becomes apparent after an extensive amount of contact breaking events; the sum of many of these singular events leads to non-trivial power law scaling of the elastic modulus with strain amplitude.
These findings provide a fresh perspective on the physics of yielding, evidence a characteristic pressure scale that distinguishes jamming-dominated systems (such as granular media and wet foams) and dense systems (such as dry foams), and suggest that these two qualitatively different classes of soft glassy materials can be distinguished by their experimentally accessible rheological fingerprint.

Macroscopic rheology: Softening and Yielding
We use oscillatory rheology to show that there are three distinct flow regimes, which we refer to as the linear response, softened, and yielding regimes. A time-varying shear strain γ = γ 0 sin (ωt) is applied to soft sphere packings consisting of N particles that are at pressure P when unsheared, and the resulting shear stress σ is measured (see Appendix A). In the linear response regime, σ is proportional to γ 0 , and the in-phase and out-of-phase components of σ /γ 0 are called the storage and loss moduli G and G , respectively; they characterize elastic stiffness and viscous damping. For finite strain amplitudes, the stress response ceases to be purely sinusoidal (see Appendix C), and the first harmonic terms in the Fourier expansion of σ are used to define the storage and loss moduli, which now depend on γ 0 . We focus on data taken at the low frequency ω = 10 −3 in units constructed from the mi- , ω = 10 −3 and P ranging from 10 −6 (black) to 10 −1 (light green) in two steps per decade. Each data point corresponds to an ensemble average of at least 25 packings at fixed P and γ. For each packing and value of P, G and G reach a clear plateau value for small γ that we use to define the linear response values G (0) and G (0) (For the variation of the linear response values with P and ω, see Appendix B). To focus on the variation of G and G with strain amplitude, we plot the rescaled elastic modulus G /G (0) and rescaled loss modulus G /G (0) as function of the strain amplitude γ 0 . (a-b) The elastic and loss moduli have a clear plateau at low strain amplitudes, before showing softening at a pressure dependent strain scale γ s , and yielding at a larger strain scale γ y . Dashed lines indicate power law decay of the moduli for large strains with slope −3/2 and −3/4 respectively. (c-d) The softening transitions in G and G collapse when plotted as function of the rescaled strain amplitude γ 0 /P. (e) Linear, softened and yielded regimes as function of the control parameters P and γ 0 . Squares and diamonds indicate yielding and softening obtained from the data for G /G (0) shown in panel (a) (see Appendix A). The dashed lines are guides to the eye, and indicate that γ y ≈ 10 −1 , γ s ≈ 10 −1 × P, leading to a crossover pressure scale P c ≈ 10 −2 .
croscopic stiffness and damping coefficients; in Appendix D we demonstrate that our conclusions are unchanged at the higher frequency 10 −2 . Moduli are reported for data recorded after the passage of an initial transient, which can be identified from the cycle-to-cycle diffusion statistics (see below).
Our simulations uncover surprisingly rich rheological scenarios, as shown in Figs. 1a-b, which display the variation of the rescaled storage modulus G (γ 0 , P)/G (0, P) and the loss modulus G (γ 0 , P)/G (0, P) as a function of strain -in Appendix B we demonstrate that the linear response quantities scale with pressure and frequency consistent with earlier predictions 10 . Our data evidences three different rheological regimes. (i) Linear response: For each pressure there is a finite strain range for which the moduli are essentially constant and equal to G (0) and G (0), indicative of linear response. (ii) Softening: Surprisingly, for low pressures and intermediate strains, both moduli fall below their linear response values but then reach new plateau values. We call this softening, to distinguish it from yielding, and associate a strain scale γ s with its onset. (iii) Yielding: At sufficiently large strain amplitudes, the elastic modulus rapidly decays, while the loss modulus peaks. We refer to this as yielding, with an associated yield strain γ y that does not strongly vary with P. For large strains the storage and loss moduli fall off rapidly as 1/γ ν 0 and 1/γ ν 0 , respectively, with ν ≈ 1.5 and ν ≈ ν /2 30,31 . Softening and yielding are distinct phenomena, each with their unique rheological signature and pressure dependence. First, the loss modulus goes down at softening, and up at yielding, signaling a qualitative difference between these two phenomena. Second, systems at lower pressures clearly soften at smaller strain amplitude, whereas the yielding strain appears pressure independent. To characterize the pressure dependence of γ s , we replot the data of Fig. 1a-b as a function of the rescaled strain γ 0 /P in Figs. 1cd. We find excellent collapse in the linear response and softening regimes of both the storage and the loss moduli for P up to 10 −2 .
Recently, several works have presented conflicting evidence for the scaling of the softening transition with pressure 8,19,32,33 , and both γ s ∼ P and γ s ∼ P 3/4 have been suggested -our data is inconsistent with an exponent 3/4.
We summarize our findings of a pressure-dependent softening transition, and a pressure-independent yielding transition in Fig. 1e. Here the red dot-dashed line indicates the softening crossover at γ s ∼ P, the red dashed line indicates the yielding crossover at γ y , and these lines meet at a characteristic pressure P c . We note that many scaling laws near jamming break down when P becomes of order 10 −2 34 , and consistent with this our data suggests that γ y ≈ 10 −1 and P c ≈ 10 −2 . Whereas such breakdown of scaling can be expected sufficiently far away from any critical point, we suggest that our data shows that there is a clear crossover pressure scale which separates the near-jamming and dense regimes; in the latter the physics is essentially pressure independent, and no longer controlled by the critical jamming point. Moreover, oscillatory rheology provides a specific experimental protocol to test which asymptotic regime is relevant for a given system: the jamming phenomenology is important when softening and yielding can be distinguished.

Reversible, trapped and diffusive dynamics
To shed light on the microscopic signatures of the rheological softening and yielding transitions, we now probe the microscopic particle trajectories. Recent years have seen considerable effort directed towards understanding the transition from reversible to irreversible particle trajectories under cyclic driving [35][36][37][38][39][40][41] . With few exceptions 42 , this work has been restricted to systems far from jamming, where the two step softening/yielding scenario identified here is absent. In that case the onset of irreversibility is found to correlate with macroscopic yielding. Here we characterize particle trajectories by stroboscopically sampling all particles at zero strain as a function of cycle number n. From these stroboscopic trajectories we compute the cycle-to-cycle squared dis-
Our data evidences three different dynamical regimes. (i): Reversible dynamics: Particle trajectories at sufficiently small strain amplitudes are reversible: ∆s 2 1 decays to the noise floor after an initial transient, and particles trace out ellipses in space consistent with a strictly linear response (Fig. 2c1). (ii): Trapped dynamics: At intermediate strains the particle trajectories are trapped: both ∆s 2 1 and ∆s 2 reach a finite plateau at large n. The particle trajectories do not form closed orbits, but remain bounded ( Fig. 2c2 and Fig. 2c3). (iii): Diffusive dynamics: For high strain amplitudes, ∆s 2 grows linearly with n and the particle motion becomes diffusive ( Fig. 2c4 and 2c5).
We have verified (via the total harmonic distortion, see below) that in the reversible regime no contact changes take place during the oscillatory driving, and that the transition from reversible to trapped dynamics occurs when contacts are broken and created during the strain oscillations. Trapped dynamics is reminiscent of caging in hard sphere glasses, although there is no ballistic motion at short times and particles generally remain in contact with multiple neighbors. We note that the loss of microscopic reversibility after a single contact change suggests that the tran-sition from reversible to caged dynamics is a finite size effect, as the strain needed to change one contact is O(1/N) in large systems (verified below). Directly obtaining the characteristic strain where the first contact changes from oscillatory rheology is numerically prohibitive, as it entails scans over N, P and γ 0 . However recent simulations of quasistatically sheared soft spheres determined the typical strain scale γ cc at which the first contact change occurs and found that it obeys finite size scaling in the linear response regime 8,24,25 : Our data for contact changes under oscillatory drive is consistent with this scaling, and we thus conjecture that the same scaling governs the strain scale where the transition from reversible to caged dynamics takes place.
The pressure dependence of the transition to diffusive dynamics can be deduced from the behavior of ∆s 2 1 and ∆s 2 as function of n; we find that the transition is essentially pressure-independent. In Fig. 3a we plot the large n plateau of ∆s 2 1 as a function of γ 0 for varying pressures. While the initial growth of ∆s 2 1 , associated with contact changes, depends on P -consistent with Eq. (1) and observations that contact breaking near jamming depends on P 8,23-26 -the asymptotic value of ∆s 2 1 becomes independent of P for large γ 0 . This increase signals a pressure-independent transition to diffusive motion, as is further evidenced by inspecting our data for ∆s 2 for all pressures, and fitting our data for ∆s 2 as n α . As Fig. 3b shows, the scaling exponent α sharply increases with γ and reaches a diffusive (α = 1) regime for large strains in an essentially pressure-independent manner.
We summarize our picture for the microscopic behaviors in Fig. 3c, and now discuss the relation with the rheological behaviors shown in Fig. 1e. Our data fully supports identifying the transition to diffusive motion with rheological yielding -both are pressure-independent and their characteristic strains are close. The correspondence between the onset of diffusion and yielding is consistent with recent experimental and numerical findings that focus on concentrated emulsions 40,42 . This link is reminiscent of the Lindemann melting criterion: once the relative particle motions reach a significant fraction of their separation at rest, structural information is lost, particles can freely diffuse, and macroscopic rigidity vanishes. A similar correspondence between diffusion and yielding also occurs in equilibrium systems, where the fluctuation-dissipation relation establishes the relation rigorously.
Surprisingly, the micro transition from reversible to trapped dynamics and the macroscopic softening are not directly linked. In fact, the transition from reversible to trapped dynamics is a finitesize artefact. This follows from the scaling of γ cc with N, which dictates that in the limit of large N, the onset strain for trapped dynamics is vanishingly small so that the reversible regime vanishes. Hence, for large systems, the transition from reversible to trapped dynamics is irrelevant, which is reasonable as strict microscopic reversibility should be absent in the thermodynamic limit. As the critical strain for softening is independent of N, for large systems there is a wide parameter range where the microscopic dynamics (1) for the transition from reversible to caged dynamics; here N 2 > N 1 . The regime P < 1/N 2 and γ 0 < 1/N 2 indicates the finite size regime where the first contact change arises at γ cc ≈ P. [43][44][45] . Outside the finite scaling regime, γ cc ∼ √ P/N and thus vanishes for large N; hence the reversible regime disappears in the thermodynamic limit, where only the trapped and diffusive regimes play a role.
is caged, but the rheology is still effectively linear.

Trapped motion and softening near jamming
In the remainder of this article, we will disentangle the relation between contact breaking and softening. Both are pressuredependent, and both have characteristic strains that vanish when P → 0 -they are thus connected to the jamming transition. The picture that will emerge is that contact changes have an O(1/N) effect on G , and occur at strains of order 1/N. Significant softening occurs after an extensive number of contact change events have happened, and this leads to a well-defined thermodynamic limit. Finally, the effect on the elastic modulus of contact breaks is cumulative, so that G decays linearly with |γ 0 |, thus signifying non-analytic behavior.
We first note that the contact change strain provides a way to rationalize the linear dependence of the softening strain on P if we postulate that contact changes are necessary for softening, so γ s ≥ γ cc . This is plausible because the shear response of packings can be mapped directly onto a spring network for infinitesimal strains. However networks stiffen rather than soften at finite strains 46 , which involve contact changes in packings but not in networks. We expect the lower bound on the softening strain will saturate for marginal packings that would unjam with the loss of just one contact. This marginal state is reached when P ∼ 1/N 2 43,45,47 . If we make the ansatz γ s ∼ P ν , the value of the exponent ν is then determined by requiring γ s ∼ γ cc when P ∼ 1/N 2 . The result is ν = 1, consistent with the observed scaling of γ s .
Returning to numerics, the link between contact changes and mechanical softening is illustrated in Fig. 4a, which focuses on a single packing. We detect contact changes via the sharp increase of anharmonic behavior of the time dependent stress signal, caused by abrupt changes in the contact stiffness when harmonic contacts open or close, and quantified by the total harmonic distortion T HD := Σ i>1 |σ 2 i |/|σ 2 1 |, where we have decomposed the stress signal in a Fourier series with coefficients σ 1 , σ 2 , . . . . At a characteristic strain γ * , Fig. 4a shows a sharp increase in the T HD over several orders of magnitude, accompanied by a small decrease in G . This corroborates our expectation that softening can only arise due to contact changes. Moreover, most of the softening of the storage modulus happens subsequent to, not at, the first contact break.
To uncover the link between contact changes and softening, we now explore the role of system size, and simulate oscillatory rheology for system sizes varying from N = 32 up to N = 2048 at fixed pressure P = 10 −3 , while varying the driving amplitude γ 0 over many decades. We found that G (γ 0 ), in particular for small systems, exhibits significant fluctuations, which requires a careful procedure to obtain meaningful averages. For each packing, we therefore define γ ∆ as the smallest strain where the deviation in G from linear response, 1 − G (γ)/G (0) reaches a value ∆, and then take ensemble averages to obtain γ ∆ (N, ∆).
In Fig. 4b we compare the mean values of γ * (closed symbols) -obtained from detecting jumps in the THD -and γ ∆ for a range of ∆ and N (open symbols). First, our data shows that γ * indeed decreases with system size, consistent with the 1/N of Eqs. 1. Second, for very small values of ∆, the data for γ ∆ closely approach γ * , consistent with the picture that any appreciable softening only arises after contact changes accumulate. Third, for large N or large ∆, γ ∆ becomes independent of N, evidencing a well defined continuum limit, where softening is due to an extensive number of contact changes.
As shown in Fig. 4c, all data for γ ∆ can be collapsed on a master curve by plotting γ ∆ /∆ as function of N∆ -data taken at different pressures shows the same trends (See Appendix D). This scaling collapse shows, first, that the typical effect of a single contact change on G is O(1/N). Second, the behavior for N∆ 1 is consistent with γ ∆ /∆ ∼ (N∆) −1 , or γ ∆ ∼ 1/N: hence, for small ∆ 1/N, γ ∆ scales as, and is close to, γ * . Third, the plateau for N∆ 1 confirms the existence of a well defined continuum limit, where γ ∆ , and hence the softening behavior, becomes independent of N and γ ∆ is linear in ∆. We conclude that softening first sets in once contacts start to break, only becomes significant when many contacts are changing, and is independent of system size for large N.
A final striking consequence of the softening being caused by  Fig. 4d we show data for ∆ for large systems and a range of pressures, which confirms this linear deviation of the plateau value of G with strain. As ∆(γ 0 ) needs to be an even function due to symmetry, this implies non-analytic behavior where ∆ ∼ |γ 0 |. We note that similar nonanalytic behavior has also been observed for strain stiffening in random spring networks at the rigidity transition 46 -consistent with the usual association of non-analyticity with a phase transition. In contrast, the non-analytic behavior found here occurs at a finite distance from the jamming point, and we suggest that it is inherited from the purely repulsive contact forces between particles, which are themselves non-analytic at the point of contact.

Discussion
Simulating large amplitude oscillatory shear, we have found evidence for two qualitatively distinct yielding scenarios for soft glassy solids. In dense systems, such as highly concentrated emulsions and Lennard-Jones glasses 48 , the macroscopic stress-strain response is linear up to the point of yielding, which occurs at a constant strain γ y . In marked contrast, weakly jammed solids such as wet foams and emulsions first soften at a pressure-dependent strain γ s , only to yield at a larger strain γ y . The ratio of γ s and γ y determines a characteristic pressure on the order of 10 −2 that marks the dividing line between these two material classes. Particle trajectories evidence an intricate link between microscopic and macroscopic behavior. The particle dynamics display just one clear transition that separates trapped and diffusive trajectories, and which corresponds to yielding. In contrast, softening has no sharp microscopic fingerprint, but results from the accumulation of an extensive number of contact changes leading to non-analytic rheological curves. Our measurements correspond well with recent experiments in emulsions 42 and pastes 49 . Softening is also observed in attractive glasses 28 and granular media 20,50 ; however the characteristic strains scale differently, likely due to attraction, friction, and/or non-harmonic contact force laws. Our data for cycle-to-cycle diffusion (Fig. 2a) are strikingly similar to data from emulsions, which also show a swift rise that sharpens with increasing P 42 .
So far, we have focused on the behavior at a single driving frequency (ω = 10 −3 ). However, our data for the characteristic strains and changes in diffusive behavior obtained at ω = 10 −2 are very close, as we detail in Appendix D. In both cases, the elastic contributions to the stress are dominant and well separated from the viscous and (yet smaller) inertial contributions. This suggests that our scenario describes a rate independent, low frequency regime. We note that recent simulations of transient rheology find similarly rate-independent characteristic strains at these frequencies 8 , and that the existence of rate independent characteristic strains is to be expected above jamming, when the system is solid.
Our findings clarify the notion of linear response 23,51 . In sufficiently large systems, vanishingly small strains lead to contact changes, perfectly reversible trajectories are not to be found, and linear response in the strict sense, in which the microscopic equations of motion can be linearized about the initial condition, is violated as soon as the first contact change occurs 8,[22][23][24][25][26]51,52 . However, we have provided conclusive evidence that the breaking of contacts does not significantly influence the macroscopic behavior 8,11,24,51 , leading to a well-defined effective linear response for macroscopic quantities.

Conflicts of interest
There are no conflicts of interest to declare.

Appendix A: Numerical Model
We perform MD simulations of 2D packings submitted to oscillatory shear. The initial packing configurations are so-called shear stabilized (SS) packings 44,45 . Unlike algorithms that relax particles within a fixed box, SS packings are equilibrated in a purely isotropic stress state; the shear modulus is guaranteed to be pos-itive, and there are no residual shear stresses. This is crucial, because packings with residual stress subjected to oscillatory driving are prone to long transients, making the use of SS packings crucial for numerical studies of oscillatory rheology.
The shear is imposed using Lees-Edwards periodic boundary conditions, and Newton's Laws are integrated using a velocity-Verlet algorithm modified for velocity dependent forces. Our packings are composed of N soft spheres of radii linearly spaced from 1 to 1.4 and mass density ρ = 1. The spheres interact through contact elastic repulsion (amplitude f el = kδ , where δ is the overlap between two particles and k = 1), as well as viscous damping ( where ∆ v is their velocity difference at contact and b = 1).
We measure the steady state stress tensor that develops in response to an imposed shear strain γ xy (t) = γ 0 sin(ωt). The first sum is over contacts, f i j is the sum of elastic and viscous contact forces between particles i and j, and r i j points between their centers. The second sum gives the kinetic stress, where m i is mass and v i is velocity. In the low-frequency parameter regimes we explore, the kinetic stress is orders of magnitude below the viscous and elastic stresses, and the latter dominates.
We determine G * via the complex ratio of the first harmonic in the stress and strain signals, carefully checking that we reach stationarity. As our numerical model is deterministic and well behaved, statistical error bars are set by numerical noise and generally very small. However, it is well known that ensembles of finite size jammed systems at fixed pressure exhibit a considerable spread in quantities such as the static shear modulus 25,[43][44][45] . Our data for G and G , both in linear response and at finite strain, show concomittantly strong ensemble fluctuations, and we report mean values G (γ)/G (0) and G (γ)/G (0) . Diffusivity: Our local dissipation law is numerically expensive but crucial to capture the correct physics 15 , and to properly resolve each oscillatory cycle, a large number of simulation steps are needed, in particular near jamming. Therefore, these simulations are numerically expensive, and pushing our numerical capacities, to obtain particle tracks we have simulated 33 realizations of systems of N = 128 particles for each value of P, Γ 0 and two values of ω; for ω = 10 −3 we have simulated 30 cycles, and for ω = 10 −2 , where transients are longer-lived, 300 cycles.
For each realization we calculate the second moment ds 2 (n, m) = (x i n − x i m ) 2 + (y i n − y i m ) 2 , where m, n are cycle numbers and · is the average over the non-rattler particles. These measures are highly sensitive to rattlers and drift in the center of mass, which is a Goldstone mode.
Before processing the position data, we first carefully identify the rattling particles at each time step and remove them from any analysis, then compute the Goldstone mode and remove it from the non-rattler positions. Nevertheless, some runs still exhibit residual drift or rattlers. As the reported second moments vary over 30 decades, drift makes the main trends less easy to observe. Therefore we focus on the median of the distributions of ds 2 (n, m), in particular the cycleto-cycle second moment ∆s 2 1 (n) := ds 2 (n − 1, n) and squared displacements ∆s 2 (n) := ds 2 (0, n) , where square brackets denote the median. In Appendix D we also show results for the mean, which show more fluctuations but do not lead to a different interpretation.
Determination of softening, yielding and diffusive onsets: We determine the data points for softening and yielding shown in Fig. 1e as follows. First, from our rheological data we estimate the strain where G /G (0) dips below 0.3 and 0.03 respectively. Second, to obtain an estimate of the critical strains that does not strongly depend on the choice of this cutoff, we assume that G /G (0) ∼ γ −3/2 , and determine the strain where G /G (0) = 1, which then gives the hypothetical intersection of the plateau at low strains and power law decay at larger strains. To determine the onset of diffusive motion shown in Fig. 3c, we take data for ∆s 2 as function of n, and perform linear fits to log(∆s 2 ) as α log(n), focusing on 3 ≤ n ≤ 30. We have checked that such fits are close to the data. The overshoot of α for intermediate strains is likely due to transients and is not expected to persist for larger n.

Appendix B: Scaling of G *
We report here our numerical results for the linear elastic and loss moduli, measured at γ 0 = 10 −10 . Here the stress signal is well described by a simple harmonic response of the form σ xy = G * γ 0 exp(iωt) + c.c.. In this regime, recent theoretical arguments predict precise scaling laws for G and G as function of P and ω 10,12 . In Fig. 5 we show our rheological data for a range of pressures and frequencies for a system of N = 1024 particles. Our data is in good agreement with the aforementioned scaling arguments: (1) all data collapses when plotted as function of ω/P; (2) for small values of ω/P, the elastic modulus G ∼ P 1/2 3,34 and the loss modulus vanishes linearly with ω as G ∼ ω/P 1/2 ; 9,10 at large ω/P, both G and G exhibit nontrivial scaling with √ ω 8-10 . We note that we can improve on the quality of the data collapse when we plot G * /P α vs ω/P β , and find the best collapse for β = 0.8, α = 0.45. We do not believe that this constitutes a significant deviation from the theoretical mean field exponents α = 1/2, β = 1 -even linear response calculations show slight deviations 10 , and data in 2D may suffer from log corrections, as 2 is believed to be the upper critical dimension for jamming 43 . Moreover, our value of β is strongly influenced by the data at the lowest ω, for which simulations are expensive and we only have a limited number of oscillation cycles. Furthermore, several effects limit our scaling range. First, for P > 0.1 we see substantial deviations from scaling, as is the case for many static properties near jamming, while for P < 10 −6 finite size effects start to dominate for our case of N = 1024 particles 8,12,[43][44][45]47 . Second, for our choice of microscopic parameters, inertial effects become detectable for ω 0.05 -if we limit our data to a smaller range of ω and P, the collapse becomes better but the scaling range shrinks.

Appendix C: Nonlinearity of response
To illustrate the nonlinearity of the large amplitude oscillatory shear response, in Fig. 6 we present Lissajous curves (stress-strain plots parameterized by time) for the same strain amplitudes presented in Fig. 2c. Anharmonic contributions are clearly visible for large strain amplitudes, but are by no means dominant, so that G and G remain meaningful.

Appendix D: Robustness of results
To show that the separation of softening and yielding is robust to changes in particle number and frequency, in Figs. 7 and 8 we show our data for G and G for N = 128 and both ω = 10 −3 and ω = 10 −2 -all features shown in Fig. 1 for N = 1024, ω = 10 −3 are also present here.