Contact Changes of Sheared Systems: Scaling, Correlations, and Mechanisms

We probe the onset and effect of contact changes in 2D soft harmonic particle packings which are sheared quasistatically under controlled strain. First, we show that in the majority of cases, the first contact changes correspond to the creation or breaking of contacts on a single particle, with contact breaking overwhelmingly likely for low pressures and/or small systems, and contact making and breaking equally likely for large pressures and in the thermodynamic limit. The statistics of the corresponding strains are near-Poissonian. The mean characteristic strains exhibit scaling with the number of particles N and pressure P, and reveal the existence of finite size effects akin to those seen for linear response quantities. Second, we show that linear response accurately predicts the strains of the first contact changes, which allows us to study the scaling of the characteristic strains of making and breaking contacts separately. Both of these show finite size scaling, and we formulate scaling arguments that are consistent with the observed behavior. Third, we probe the effect of the first contact change on the shear modulus G, and show in detail how the variation of G remains smooth and bounded in the large system size limit: even though contact changes occur then at vanishingly small strains, their cumulative effect, even at a fixed value of the strain, are limited, so that effectively, linear response remains well-defined. Fourth, we explore multiple contact changes under shear, and find strong and surprising correlations between alternating making and breaking events. Fifth, we show that by making a link with extremal statistics, our data is consistent with a very slow crossover to self averaging with system size, so that the thermodynamic limit is reached much more slowly than expected based on finite size scaling of elastic quantities or contact breaking strains.

PACS numbers: 83.80.Fg, 83. 10.Rs, 62.20.fg How does a jammed system fail? Failure of amorphous systems under increasing driving generally leads to a complex chain of events, where an initial linear response gets gradually eroded by local micro events that lead to plasticity and eventually organize in persistent flows [3][4][5][6][7][8][9][10][11][12]. For systems near the critical jamming point, the question of failure is even more vexing, as the characteristic strain for the first deviations from linear response is vanishing, both with the number of particles in the system N, but also when the confining pressure P is lowered towards the critical jamming point. Moreover, near the unjamming point disordered solids are extremely fragile, and the tiniest of perturbations can cause an intrinsically nonlinear response [9,10,[13][14][15][16][17]. Hence, one may question the validity of linear response for athermal amorphous solids, as the range of validity may vanish [7,14,18,19]. Finally, the unjamming transition at vanishing P bears hallmarks of a critical phase transition: properties such as the contact number and elastic moduli exhibit power law scaling [19][20][21][22][23][24][25][26][27][28], time and length scales diverge [19,24,25,[28][29][30][31], the material's response becomes singularly non-affine [31,32] and finite size scaling governs the behavior for small numbers of particles N and/or small P [1,2,33]. The question we want to address is how, near jamming, when linear response vanishes and criticality dominates, a jammed system reacts and fails under increasing driving.
Earlier work on contact changes has focused on vibrations [14,34,35] or hard particle systems [18,36]. We instead focus on soft particle systems, as they are descriptive for a wider range of experimentally relevant systems, use experimentally relevant simple shear deformations, and focus on the (a) The first contact change in a sheared packing (N = 64, P = 10 −6 ) occurs at a strain γ * = 9.003851(2) × 10 −7 , when the two marked particles lose their contact. (b) The corresponding stress-strain curve remains continuous but exhibits a sharp kink; we define G 0 as the shear modulus of the undeformed packing, and G 1 as the shear modulus of the packing just above γ * . first unambiguous deviation from strict linear response: contact changes under quasistatic shear (Fig. 1a) [37].
We address the following questions: (i) What is the nature of the first contact changes near jamming? In systems far from jamming, rearrangements organize into avalanches: collective, plastic events in which multiple contacts are broken and formed and the stresses exhibit discontinuous drops [4][5][6][7]38]. For hard particles which represent a singular case where motion always involves unjamming, even a single contact break may induce a complete loss of rigidity [7,13,18,39]. In contrast, we find that near jamming the first events are the making or breaking of a single contact, and that the stress remains continuous. The probabilities for contact making and breaking are governed by finite size scaling, with making and breaking equally likely for N 2 P 1, but contact breaking dominant for N 2 P 1.
(ii) What is the mean strain γ cc at which the first contact change arises? What are the mean strains of the first contact breaking γ bk or contact making γ mk events? We first show that we can use linear response calculations to accurately capture these strains, and then show that all these characteristic strains vanishes when either N diverges or P vanishes. All strains obey finite size scaling: γ cc ∼ γ bk ∼ P and γ mk ∼ 1/N 2 for small systems close to jamming (N 2 P 1), whereas and γ cc ∼ γ bk ∼ γ mk ∼ √ P/N for N 2 P 1. As log-corrections to scaling are expected for jamming in 2D [1,2,40], and in additional alternative corrections to scaling have been proposed [13,36], we carefully study our data from this perspective, and find that our data is consistent with both -in 2D, extremely large systems are needed to distinguish between these different corrections.
(iii) How do contact changes affect linear response? For finite systems close to jamming, even a single contact change can strongly affect the elastic response (Fig. 1b). Clearly, calculations based on the Hessian matrix of the undeformed packing are then no longer strictly valid. As a result, the relevance of the linear response scaling relations are currently under dispute for systems close to jamming, at finite temperature, or in the thermodynamic limit [14,34,35,[41][42][43]. By comparing the shear modulus before (G 0 ) and after (G 1 ) the first contact change, we find that their ratio again is governed by finite size scaling, and while the ratio G 1 /G 0 approaches 0.2 for small N 2 P, for large N 2 P, G 1 /G 0 → 1. We also study the statistics of G 1 /G 0 by the standard deviation σ of its distribution, and find three regimes: for small N 2 P, σ ≈ 0.3, for N 2 P ≈ 1, the fluctuations are strongest and values of G 1 < 0 are most likely, whereas for large N 2 P, σ scales roughly as [N 2 P log 10 (N) −0.7 ] 0. 35 . The latter scaling allows us to estimate the cumulative effect of a diverging number of contact changes that occur when the strain is fixed and N → ∞, and shows that this is limited: effective linear response, quantified by the shear modulus at finite strain, appears well-defined.
(iv) We explore sequences of multiple contact changes under shear, and find strong correlations between alternating making and breaking events. A surprising effect is that while initial contact breakings drive the system precariously close to catastrophic failure (too few contacts to maintain rigidity), the subsequent sequence of contact making and breaking extends the range before such failure sets in.
(v) Fifth, we show that by making a link with extremal statistics, our data is consistent with a very slow crossover to self averaging with system size, so that the thermodynamic limit is reached much more slowly than expected based on finite size scaling of elastic quantities or contact breaking strains.
Our work paints a clear and coherent picture of the role of contact changes near the critical jamming point. While the range of strict validity of linear response vanishes for small P and large N, macroscopic quantities such as the shear mod- ulus are relatively insensitive to contact changes as long as P 1/N 2 . Hence, linear response quantities remain relevant for finite P and large N, while for P 1/N 2 , a single contact change already changes the packing significantly. The qualitative differences in the nature of contact changes close to and far from jamming suggests that plasticity, creep, and flow near jamming are controlled by fundamentally different mechanisms than plastic flows in systems far from jamming [4,6,7,[9][10][11][12]38].

I. METHOD & PROTOCOLS
We simulate bidisperse packings of massless, frictionless soft spheres in two dimensions [21,23,24]. Recently, it was shown that such finite packings are not guaranteed to have positive shear moduli, nor have zero residual stress [2,33,44], which both could lead to problems when studying contact changes. We therefore focus on so-called ε + all packings that have positive moduli and zero residual shear stress as described in [2,33]. In Appendix A, we describe in detail how to create and shear such packings, which in particular necessitates the use of non-square unit cells [2,33,44]. Here, we will focus on our algorithm to detect contact changes.
To find contact changes, we apply a strain (Eq. A2) where we increase ζ = 0, 1, ... until we detect a change in the contact network (δ i j = 0 ↔ δ i j > 0 for any pair i, j). We then move back to the state before the contact change, and use bisection to determine the strain at the contact change γ * until ∆γ/γ * < 10 −6 . Rattlers require special attention: because they are free to move, their behavior is ill-defined. In our simulations, we encounter rattlers in two distinct types of events. First, rattlers may become part of the load-bearing network. As the rattlers' position is ill-defined, the strain at which this occurs is algorithm dependent. We therefore exclude such contact making events in our analysis of the first contact change. Second, a particle with three contacts can become a rattler, where force balance dictates that all three contacts go to zero overlap simultaneously. This is detected correctly in our simulations, and the event is recorded as a single break event. In linear response calculations (to be discussed below), the creation of a rattler is also well-defined. In Fig. 2b, we show the overlap δ ri of particle r with its neighbours A, B and C. In the simulations (symbols), we find the overlaps smoothly go to zero while approaching the contact change strain γ * . In linear response, calculated at γ = 0, we find a slightly different contact change strain for each contact, but they are within |∆γ/γ * | < 10 −4 .

II. NUMERICAL RESULTS
In this section, we discuss the results of direct numerical simulations to determine the properties of the strain γ * at which the first contact change occurs. We first discuss the relative prevalence of contact making and breaking events. We then study in detail the statistics of γ * at given P and N, and finally discuss how the ensemble averages γ cc = γ * scale with N and P.
A. The first contact change For each packing in an (N, P) ensemble, we determine the strain of the first contact change γ * , as described in Sec. A. In Fig. 3a we show the cumulative distribution function (cdf) of γ * for N = 256 ensembles at various pressures. We observe that, first, the typical scale of the strain γ * increases with pressure P, and secondly that their shape is mostly independent of P.
In Fig. 3b, we show a stacked probability graph of the different contact change types. We distinguish events where one or more contacts are broken (break), events where one or more contacts are created (make) and events where contacts are both broken and created (mixed). The number of mixed events in- creases with pressure, but is less than 5%, independent of N.
Within the make class, we can distinguish events where a particle which originally was a rattler now becomes part of the contact network (make (rattler)). Of all make events, 5 − 15% involve rattlers. This is consistent between ensembles, with no clear dependence on either N or P. At low pressures, we find that the vast majority of events consists of contacts being broken. At large pressures, we find that roughly half of the events create a new contact. In Sec. III C, we will show how these probabilities vary as a function of N 2 P. In the remainder of this paper, we will focus on the simple make and break cases.

B. Strain distributions
We now take a more detailed look at the distributions of γ * and show that contact changes can essentially be described as a Poisson process, as the cdf close resembles an exponential distribution, with Pr(γ * ≤ γ) = 1 − e −γ/β . In Fig. 4a we show Pr(γ * > k γ * ), i.e. the complimentary cdf of γ * , rescaled by the ensemble mean γ * . If γ * is exponentially distributed, the ccdf is a simple exponential: Pr(γ * > k · γ * ) = e −k (k ≥ 0), and as Fig. 4a shows, our distributions for N = 256 are close to exponential. This is consistent with a Poisson process, where contact changes are independent of each other.
To check conformance to an exponential distribution as a function of N and P, we use the Anderson-Darling test [45], with which we test the hypothesis "these values of γ * were drawn from an exponential distribution". We use a 5% confidence interval, i.e., there is a 5% probability we reject the hypothesis for samples that were drawn from an exponential distribution. In Fig. 4c, we show the results of this test. We observe deviations from exponential behavior for small systems and low pressures. The boundary between rejection and non-rejection corresponds with the transition between systems for which finite size effects dominate and large systems, at N 2 P ≈ 1 [1,2,37]. This suggests that for large systems (N 2 P 1), contact changes are uncorrelated, while for small systems, correlations build up.
How do distributions for systems in the finite size regime deviate from exponential? In Fig. 4b, we show rescaled ccdfs for N = 16 systems at various pressures. The most significant deviation is at low k, where we find Pr(γ * > k · γ * ) is larger than expected for an exponential distribution. As Pr(γ * > k · γ * ) is the survival probability, this indicates a lack of events at small strain, which means that, in small systems, events are antibunched. Notwithstanding this deviation from exponential behavior, the mean remains well-defined, as is further evidenced by recent work which shows the number of contact changes scales linearly with strain [19].

C. Scaling
We now discuss the variation of the mean contact change strain γ cc = γ * with N and P. As discussed in [37], we can obtain data collapse for γ cc when we plot N 2 γ cc as a function of N 2 P. As shown in Fig. 5a, this results in a good (but not great) data collapse. It has been suggested that the upper critical dimension for jamming is two, which implies logarithmic corrections to scaling [2]. Using the form suggested in [2], we find a very good data collapse (Fig. 5b).
How do we think about these strains? As we will show later, strictly linear response captures the deformations well up to the first contact change. It is thus useful to consider, on the one hand, the overlaps and "underlaps" between pairs of particles in (near) contact and, on the other hand, the relative motion of such pairs. The former are set by the packing, and in particular, the overlaps scale trivially with the pressure. The latter follow from the full linear response via the set of eigenmodes that characterize the system (Appendix C). This way of thinking strongly suggests that we should consider the behavior for N 2 P smaller or larger than one separately. For N 2 P 1, the number of contacts is constant, and the eigenmodes are essentially independent of P (Appendix C). Hence, here the main variation with P is in the overlaps, which vanish when P → 0. Therefore, we expect breaking to happen at much smaller strains than making, and hence that the contact change strain is simply linear in P -consistent with the data in Fig. 5. Moreover, this simple picture suggests that the amount of shear stress at the first contact change is proportional to P.
The situation for N 2 P 1 is more complex, because here the eigenmode spectrum changes with P, and indeed, it is known that the relative motions normal and transverse to a contact pair's center-to-center line scale as u ∼ P 1/4 γ and u ⊥ ∼ γ/P 1/4 , respectively [32]. As the transverse motion diverges near jamming (for large N 2 P), it dominates the change δ ∼ u 2 ⊥ / in the center-to-center distance . A naïve argument for the pressure dependence of the breaking strain can then be constructed by balancing δ with the typical overlap in the initial condition, δ (γ bk ) ∼ δ, yielding the prediction γ bk ∼ P 3/4 . Indeed, a strain proportional to P 3/4 also arises in a recent scaling theory of the jamming transition [46]. While this argument correctly predicts nontrivial P-dependence in the characteristic strains for N 2 P 1, the 3/4 exponent is inconsistent with our data shown in Fig. 5. We believe the essence of this discrepancy is that the assumption that the first broken contact is typical of all contacts is incorrect. We note, in passing, that studies that consider "typical" contacts, do find a scaling consistent with a 3/4 exponent [14].
After this introduction, we we now discuss two distinct arguments that both lead to a scaling relation which is consistent with our data: an argument for compressive strain, and a stress-based argument for shear strain.
Compression: We start with a compressional argument, based on estimating the strain scale for making and breaking a contact under compression. There is a clear relationship between compression and the number of contacts: we gain contacts if we compress the system and we lose contacts if we expand the system. The scaling relation that relates the excess contact number N exc = N∆z/2 to N and P is well-known from earlier work [1,2,33], and is shown in Fig. 6. There are two branches: a plateau N exc ∼ 1 at low pressures and a square root pressure dependence N exc ∼ √ N 2 P at higher pressures.
How far do we need to expand or compress a system at given N and P to induce a contact change? In the highpressure regime, the derivative ± ∂ ∂P √ N 2 P ∼ ±N/ √ P gives the number of contacts changed due to unit pressure change.
Its inverse δP ∼ ± √ P/N, then gives the pressure change needed for a single contact change. The compressional strain is the pressure change divided by the bulk modulus K: ε cc ∼ ±δP/K. As K is independent of N and P [24], we simply find ε cc ∼ ± √ P/N.
In the low-pressure finite size regime, the number of contacts is independent of pressure. Nevertheless, the plateau has a finite length. On the one hand, the plateau ends at P = 0, as we unjam our system and lose all contacts. On the other hand, the plateau ends when we enter the large system size regime at N 2 P ∼ 1 and gain one new contact.
The scales for making and breaking a contact are thus no longer the same in the finite size regime: To break a contact, we unjam the system by reducing the pressure with δP ∼ P, and we find ε bk ∼ −P. To create a contact, we increase pressure up to the beginning of the large system regime, at P target = 1/N 2 . As we are initially in the small system regime, the current pressure P 1/N 2 and can be neglected, and the pressure change δP = P target − P ≈ −1/N 2 . We thus need to apply a strain ε mk ∼ −1/N 2 . The contact change strain, independent of direction, will be given by the minimum of the absolute making and breaking strains. As P 1/N 2 , we thus expect ε cc ∼ P.
Summarized, this argument leads to these characteristic strains for contact changes under compression: As we will see, arguments based on shear, as well as our our results, find the same scaling for these strains. 5248860101717435 Shear: We can also formulate an argument for the scaling of γ cc under shear from dimensional analysis. Other than taking γ cc constant, there is no clear strain scale, so we will construct the argument using stress instead. We will start by determining the typical stress scale σ cc .
There are three stress scales in the system: the confining pressure P, the bulk modulus K and the shear modulus G. As we are describing shear, it seems unlikely that K is relevant. If the stress scale σ cc were to scale with G, we would end up with a constant strain, and we have already seen that γ cc is not constant. This suggests that the only relevant stress scale is the confining pressure P, which we already have argued to govern the behavior for N 2 P 1; we now assume it also to govern the large system limit, and take σ cc ∼ P. The stress scale must also depend on the system size. Say we have a packing with N particles, which has a contact change at σ = σ cc . If we duplicate this system, we will have 2N particles, and two contact changes will have happened at the same stress σ cc , so that we expect that σ cc ∼ 1/N. Combining these two scalings leads to the following suggested scaling: We determine the strain scale γ cc via the shear modulus G = σ/γ. From earlier work [1,21,33], we know G scales as which, combined with the stress scaling we derived, suggests the following scaling for γ cc : consistent with the scaling proposed in Eq. (2) Finally, we note that Eq. (3) suggests to plot the stress at the first contact change, σ cc as a function of P/N. As shown in Fig. 7, this gives a reasonable, but not excellent, data collapse. Nevertheless, the quality of the scaling collapse of γ cc shows that ultimately the proposed scaling is correct, despite the hand waving nature of the underlying arguments to derive it.

III. LINEAR RESPONSE
We now show and utilize that many properties of the first contact change can be deduced from the initial state at γ = 0 using linear response. The idea is to estimate the trajectories of (non-rattler) particles from their linear elastic response: is calculated at the initial state. From the linear trajectories, we extract the variation of all overlaps (contacts) and underlaps (gaps between particles) with strain. Contact changes then correspond to sign changes of the overlaps and underlaps. As we will see, this strategy not only allows us to accurately obtain the strain for the first contact change, but also gives us insight into the microscopic mechanisms. In particular, linear response allows us to probe the closing of contacts in detail, which is difficult in direct numerical simulations (DNS) since, at low N 2 P, it becomes exceedingly rare for the first contact change to be a closing event (Fig. 3b). In this picture, the contact changes stem from a combination of geometric and linear response properties not explicitly considered before.
In this section, we show that the response remains essentially linear up to the first contact change: the nonlinear behavior of jammed packings under deformations arises mainly due to the cumulative effects of many contact changes. First, the stess-strain response is essentially linear between contact changes (III A). Then, we show that linear response predicts the contact change strains with surprising accuracy (III B): Linear response predicts its own demise. Finally, we investigate the first breaking and first closing events according to linear response (III C).

A. Stress response
First, we will show that the stress-strain response of our systems is essentially linear in the DNS simulations up to the first contact change. From the simulations, we obtain the shear stress σ(γ) at various strains before the first contact change (Figs. 8a-b). We fit this response with a second-order polynomial σ = G 2 γ + λγ 2 , and quantify the relative contribution of the quadratic component as the ratio between the quadratic and linear contributions at γ * : For a given N and P, the fluctuations in Q are much larger than the mean. Hence, the relative importance of the nonlinearities is given by the width of the distribution P(Q). Our data indicates that these distributions exhibits fat tails, i.e. decays significantly slower than exponential, and that the second moment of Q is ill-defined. Therefore, we characterize the width of P(Q) by halve of the 16%-84% width that we denote S Q -for Gaussian distributions this corresponds to one standard deviation. We have checked that S Q allows to collapse the CDFs of Q (such integrals over P(Q) are more robust to small sample fluctuations than PDFs), such that S Q presents a robust measure of the fluctuations and magnitude of Q. In Fig. 9, we plot S Q as function of N and P. The most important observation is that S Q remains small in the vast majority of cases, and the strongly nonlinear case shown in Fig. 8b is truly exceptional. The two regions where S Q appears to be largest are for small N and large P, and for large N and small P. The origins for these deviations are different. For large systems at low pressure, the larger deviation is caused by inherent nonlinearities in the system. Small systems at high pressures exhibit also significant deviations from linear response, as the characteristic strains at the first contact change become large when N is small and P is large. Nevertheless, the quadratic contribution to the stress, λγ 2 is small compared to the linear contribution G 2 γ, and we therefore expect to be able to predict the response of the system directly from linear response.

B. Contact change strains
In this section we describe how to calculate the contact change strains from linear response, and compare these values to the results from direct numerical simulations. First, for for various system sizes as indicated at (a) P = 10 −2 and (b) P = 10 −6 . For each PDF, the standard deviation σ is indicated. each particle pair i, j, we determine the contact change strain γ i j , defined as the strain where the particles, assuming linear trajectories, break contact or make a new contact. By minimizing over all these strains, we calculate the strain at which the first new contact is made γ LR * ,mk , the strain at which the first contact breaks γ LR * ,mk , and their minimum gives the strain at the first contact change γ LR * . We then, for each packing, compare these values to their counterparts obtained by simulations.
Calculating γ LR * : For each particle pair i, j, we determine the center-to-center distance r i j , and use linear response at γ = 0 to determine u i = ∂ x i /∂γ. The inter-particle velocities are then given by u i j = u i − u j − n y,i j L yyx , where the last term incorporates the velocity between the copies of the periodic box. Combining these, we can solve | r i j + γ i j u i j | = R i + R j for γ i j to determine when the overlap δ i j = 0.
We determine the first broken and closed contact independently: γ LR * ,mk ≡ min i, j not in contact which allows us to study opening and closing events directly and independently, which is impossible in DNS simulations. The first contact change for the entire system is then determined by taking the minimum of the strain over all particle pairs i, j: Comparison with DNS simulations: We now show that linear response accurately predicts the contact change strain. For each individual system, we compare the linear response values γ LR * to the corresponding strain γ DNS * from the DNS simulations. In Fig. 10, we plot pdfs of γ LR * /γ DNS * to quantify the relative deviation from the simulation. We observe that γ LR * is a good predictor for γ DNS * . First, these distributions are peaked around 1, which shows the mean strain found in linear response matches that of the simulations very well. Secondly, the standard deviation of the distributions, σ is of the order of 5% for small systems and 1% for large systems. At P = 10 −2 , the largest packings have a standard deviation of 7 · 10 −3 , which increases to 5 · 10 −2 for small systems. The largest standard deviation is obtained for very small systems (N = 16) at high P (10 −2 ). We find a large dependency on pressure: for P = 10 −6 , the distributions become very narrow around 1. The standard deviation remains on the order of 10 −2 due to outliers. We conclude that for all parameters considered, the differences between the strains obtained by linear response and direct numerical simulation are small. In addition to determining the right contact change strain, we found that in over 90% of cases linear response also correctly identifies the contact i, j where the first contact change takes place.
In conclusion, linear response provides us with a powerful tool to predict the behavior of packings. It allows us to predict the correct first contact change, as well as determining microscopic properties unavailable in the DNS simulations. We note in passing that the correct prediction of contact changes suggests that shearing jammed packings might be modeled in terms of a discrete event simulation, where, instead of slowly stepping through strain space, we immediately jump from contact change to contact change.

C. Scaling of ensemble averages obtained in linear response
We now use linear response to study the strains at which contacts are broken or created in detail. Based on Eq. 2, we expect three scaling regimes for the contact change strains: for low N 2 P, γ mk ∼ 1/N 2 and γ bk ∼ P, while for high N 2 P, both γ bk and γ mk are expected to scale as P/N 2 . As before, these scalings suggest scaling collapse if we plot N 2 γ as a function of N 2 P: In Fig. 11a, we plot our linear response data using this rescaling. As in Sec. II C, applying log corrections [2] improves the collapse. For low N 2 P, we find that the data is well described by the expected power laws γ mk ∼ (N 2 P) 0 and γ bk ∼ (N 2 P) 1 . For high N 2 P, we find that neither branch cleanly scales as γ ∼ √ N 2 P. Nevertheless, the linear response data support the scaling arguments, and in particular reveal the plateau for γ mk that cannot be obtained in DNS. We expect that for larger systems the clean square root scaling will be recovered for both branches, as here both N 2 P can be large while P remains small. We have seen that linear response provides us with a powerful tool to understand what happens in the simulations. We not only predict the first contact change with surprising accuracy, we can also capture the prevalence of different types of events.

D. Log-corrections versus freely adjustable exponents
Here we investigate how accurately we can determine the power laws via scaling collapse of our data, and compare the log corrections we applied in Sec. II C to power law corrections. In Sec. II C, we provided three arguments that predict the following scaling for the first contact change strain γ cc : where F(x) ∼ x for small N 2 P and F(x) ∼ x 0.5 for large N 2 P.
In the same section, we have seen the results from the simulation collapse when plotted in this way. Furthermore, we have seen that by adding the log correction with the same F(x) the collapse improves. First, we investigate for which exponents in N the collapse, without the log correction, is satisfactory, i.e., for what values of q and r does give an acceptable collapse? To make this quantitative, we measure the running maximum (starting at low N r P) and the running minimum (starting at high N r P), and calculate the effective area between the curves ℵ = log 10 (M(N r P)) − log 10 (m(N r P)) d log 10 (N r P), In Fig. 12, we show collapse plots for q = 1.8 . . . 2.2 and r = 1.6 . . . 2.0. We observe that all plots with are reasonable (ℵ 1), and that N 2 γ ∼ F N 1.8 P has the best overall scaling collapse (ℵ = 0.32). Our log-corrected collapse is very close to this, with ℵ = 0.37.
Secondly, we can wonder about the correct asymptotical behavior of F(x). To find this behavior, we fit F(x) = C · x β separately for both the upper (N 1.8 P > 10) and lower (N 1.8 P < 0.1) branches (Fig. 13a). Here, we find which means that the best overall scaling of γ becomes The error bars are given by the variation of the parameters when the fit range is increased or decreased by a decade. When p and q are varied within the collapse region, the exponents vary by ∼ ±0.05. When we compare the power laws to our expected scaling, we find the scaling of γ with P is as expected, but note two differences from the expected scaling of γ with N. First, we observe γ decreases as N −0.17 for small systems, instead of the independence of N our scaling model predicted. Secondly, for large systems, we observe γ cc scales as N −1.1 instead of N −1 .
Comparison between power law and log corrections: We can interpret the 1.8 exponent in N as a correction to the predicted N 2 P scaling: N 1.8 P = N −0.2 (N 2 P). In Fig. 13c, we compare this correction to the log correction described in Sec. II C. We observe the corrections produce largely the same effect in the range of N that our simulations cover. When we plot the ratio of the two (Fig. 13d), we observe that the deviations between both corrections are less than 35%, over a range where N 2 changes by three orders of magnitude. To achieve a measurable difference of a factor three, systems of at least 60000 particles are required. Alternatively, simulations can be performed in three dimensions, in which case the log corrections disappear [2]. As the variation in the quality of the collapse is small, caution is warranted.
To conclude, we find our deviations from the expected scaling can be described by both a log correction and a power law correction. Much larger or three-dimensional simulations are required to fully distinguish the two corrections.

IV. MULTIPLE CONTACT CHANGES
In this section, we discuss the behavior of our systems when they are strained beyond their first contact change, focussing on the the implications of contact changes for continuum elasticity, and reveal intriguing patterns of subsequent make and break events.

A. Shear modulus
As we have seen, the first contact change happens at lower and lower strains as systems get larger. Schreck et al. [14] suggested that this implies that linear response is no longer valid for disordered systems at large N. It is clear that changing a single contact can have a large effect on small systems, but one would expect the effect to vanish in larger systems: in the thermodynamic limit, systems are expected to behave increasingly like an elastic solid, and this apparent paradox lead to a lively debate [34,35,43].
Here we show how the effect of a single contact change on the shear modulus becomes smaller and smaller when the system size is increased. We note that, as long as the shear modulus does not change significantly, we can consider the system to have an effective linear response, even though it is no longer strictly linear. To quantify the effect of a single contact change, we calculated the shear modulus before (G 0 ) and after (G 1 ) the first contact change using Eq. B16. For each value of N and P, we have calculated the probability distributions ρ(G 1 /G 0 ), and from these determine in particular ρ(G 1 < 0), and the width of these distributions (Fig. 14a). We find that the shape of these distributions varies strongly and that we can organize our data using the finite size parameter N 2 P log 10 (N) −0.7 , and as function of this parameter we distinguish three regimes.
(i) N 2 P log 10 (N) −0.7 1: In the small system size limit, we find that ρ(G 1 /G 0 ) is a strongly asymmetric distribution, with most weight around zero. We find that the mean G 1 /G 0 ≈ 0.2, and that 0 < G 1 < G 0 . To understand this, we note that in this regime, the first contact change is a breaking event, which weakens the system. We find that G 1 is significantly smaller than G 0 because, in this regime, there is typically only a single excess contact (N c −2N = 1). Surprisingly, the system does not unjam immediately, for reasons we will discuss in Sec. IV B.
(ii) N 2 P log 10 (N) −0.7 ≈ 1: In the intermediate regime, the number of excess contacts remains small, contact changes are predominantly contact breaking events, and we observe that G 1 < G 0 . However, the probability that G 1 < 0 becomes finite, inc contrast to the behavior in regime (i). This follows from the variation of prestress: without prestress, G has to be non-negative [2,29], but as P increases in regime (ii) there is sufficient prestress to allow for negative values of G 1 , in up to 35% of cases (Fig. 14b).
(iii) N 2 P log 10 (N) −0.7 1: For large systems, we enter the continuum regime, where the distribution ρ(G 1 /G 0 ) peaks around one and becomes increasingly symmetric and narrow. Hence G 1 ≈ G 0 , and this is the essence of the solution of the apparent paradox. The symmetry of the distribution is consistent with out observation that contact creation and contact breaking becomes equally likely in this regime.
A simple scaling argument for the width of this distribution can be obtained from combining the scaling of G with P, G ∼ ∆z ∼ √ P with the observation that making and breaking of contacts is equally likely. As a single contact change modifies ∆z by ±1/N, we thus expect G ± 1 ∼ ∆z 0 ± 1/N. The width of this distribution scales as We measured the width of this distribution using the standard deviation σ, and observe that it vanishes as (N 2 P log 10 (N) −0.7 ) −β with β = 0.35 ± 0.01 (Fig. 14c). We suggest that the contacts changed under a shear deformation have a relatively large impact on the shear modulus -a relatively small number of contacts contribute disproportionally to the elastic moduli [47].
Nevertheless, the observed diminishing of the width of the distribution ρ(G 1 /G 0 ) is sufficiently strong to be consistent with an effective linear response picture. We call a material effectively linear if, for a small fixed deformation γ t , the standard deviation of G(γ t ) vanishes for N → ∞. In terms of contact changes, we thus need to establish how the number of contact changes experienced up to γ t grows with N, and how the effect of single contact changes decreases with N. We estimate the number of contact changes between γ = 0 and the test strain γ t as n = γ t /γ cc = γ t /( √ P/N).
We then assume that all contact changes are independent of each other, and assume each contact change causes a change in G drawn from the distribution ρ(G 1 /G 0 ) with standard deviation σ ∼ (N 2 P) −β . The central limit theorem then states the standard deviation after n contact changes is given by Combining these, we find that the standard deviation after a strain γ t is given by which vanishes for large N as long as 1 2 − 2β < 0, or Clearly, 0.35 > 1/4, so, for N → ∞, our systems approach the continuum limit. Significant correlations between subsequent values of G i+1/G i could in principle lead to a more problematic approach to the continuum limit. However, recent work by Boschan et al. [19] found that the ensemble-averaged stress-strain curve is linear with a slope compatible with G 0 up to a strain of order P. Though not a definitive test, on the basis of these results we consider that strong correlations are unlikely to be present. Our data is thus consistent with the picture where, for large N, the effective value of G depends on the applied shear γ rather than the number of contact changes n [19,48].

B. Alternating contact changes
Here, we investigate correlations between consecutive contact changes, focussing on the N 2 P 1 regime. In Fig. 15, we show the number of contacts in the system, N c , as a function of the number of contact changes for systems with N = 16 particles, at P = 10 −6 . Before shearing, N c reflects the number of rattlers, with N c = 33, 31, 29 corresponding to zero, one and two rattlers, respectively. The presence of these rattlers accounts for the parallel tracks in the dominant transition pathways ( Fig. 15). For definiteness, let us focus on the case where the initial packing has no rattlers (N c = 33). When the system is sheared, rattlers occasionally form, and N c is then seen to drop by three. In roughly one in ten packings, shearing cause three successive breaking events, causing the system to unjam. In most cases, however, we find that there are first two breaking events, followed by a series of alternating making and breaking events; clearly, throughout this process the pressure remains finite and the system remains jammed. This The circle area represents the fraction of systems with a given number of contacts; the thickness of the lines represent transition probabilities. Initially, the systems start off with the minimum number of contacts 2N + 1 = 33 (31 or 29 when there are one or two rattlers, respectively). In the first and second contact change, the system loses one contact (three when a rattler is created). In the following events, the system alternately gains and loses a contact. alternating behavior stays apparent at least until the 10 th contact change. This evidences correlations between subsequent events. We note that for larger pressures (N 2 P > 1), such correlations are absent.
To interpret the values of N c , we recall that the initial condition of these simulations are ε + all packings that have positive moduli and zero residual shear stress; for these packings it is well known that the minimal number of contacts equals 2N + 1, consistent with the initial values of 33, 31, . . . observed here [2,33]. The reason the system under shear remains jammed for lower contact numbers, is that the boundary conditions during shear, and during initial equilibration are different. Once the system is equilibrated, the box shape parameters α and δ (see Appendix A) are fixed, the system has two degrees of freedom less, and can remain jammed down to N c = 2N − 1 [2,33]. The situation is somewhat subtle though. We have observed that whether we fix the simulation box volume (as shown in Fig. 15) or fix the pressure does not change the minimal contact number during shear. However, if we fix the deviatoric (pure shear) stress τ = (1/2) σ xx − σ yy instead of the pure shear strain δ = (L yy − L)/L = L yy /L xx − 1 , we find that the minimal contact number is 2N instead of 2N − 1.
We note that the same contact is often involved in multiple contact changes, although typically not in subsequent contact changes. This is an example of the intriguing correlations in the spatiotemporal patterns of contact changes that invite further studies. We already discussed one aspect of the boundary conditions. In the constant volume protocol, the pressure increases with shear due to dilatancy -for the example shown in Fig. 15, the pressure becomes of order 10 −4 in the strain interval leading to the first contact creation event. How simulations at constant pressure, and/or constant shear stress influence this phenomenology is an open question.  Fig. 11b. (c) The ratio γ dist bk /γ LR bk varies slowly with N 2 P log 10 (N) −0.7 , from γ dist bk /γ LR bk ≈ 0.5 to γ dist bk /γ LR bk ≈ 1.0.

V. EXTREMAL VALUE SCALING
Second, we will approach the problem from a statistical perspective. Starting from the distribution of γ * of all contacts in all packings, we apply extreme value analysis to find the expected mean first contact change. We find that this does not yield a good prediction for the measured value, and determine that this cannot be explained by a few weak contacts, but rather points to strong correlations involving the whole system -i.e., the statistics of the first n changes in the system are different from the statistics of the first contact change in n systems.
In this section we probe whether we can predict the scaling of γ cc and distribution of γ * based on the distribution of all contact change strains ρ(γ i j ) for a given ensemble (N, P). Note that before (Secs. II and III), we have determined the scaling of γ cc by determining γ * for each packing, and averaging over those values. We have found that the distribution of γ * is close to a exponential distribution. Assuming that large enough packings are statistically similar, it should be possible to predict γ cc from the distribution of ρ(γ i j ) using extremal statistics. In particular, one might expect that ρ(γ i j ) takes on a simple form for sufficiently large N, possibly even amenable to a theoretical description. Deviations from this picture may point to lack of self-averaging or other subtleties, and as such provide important information for developing a deeper theoretical understanding for the characteristic strains of the first contact change. Before starting, we note that for contact creation, it is difficult to establish which potential contacts should be considered, and we therefore focus on the breaking of contacts only, using γ LR * ,mk from linear response. We will also limit our discussion to contacts that break for shear in the positive direction, i.e., γ > 0. As a first probe of the usefulness of extremal value statistics for contact breaking, we compare the results of two distinct methods to calculate the mean contact breaking strain. First, we define γ LR bk = γ LR * ,mk , the mean of the contact breaking strains determined for an ensemble of pack-ings, as we have done in Sec. III. Second, we determine γ dist bk from the distribution of positive contact change strains ρ(γ † i j ) by solving To implement this, we first compute the numerical cdfPr(γ † i j < γ) based on the breaking strain γ i j for every contact in every packing in the ensemble and then solve where N bk ≈ 0.5 N c is the mean number of contacts that break under positive strain, for which we take the numerical ensemble average. This procedure is illustrated in Fig. 16a for the N = 1024, P = 10 −2 ensemble, where N bk = 1147. For this particular example we find that γ LR bk = 1.5 × 10 −4 whereas γ dist bk = 1.4 × 10 −4 . These values are close but distinct (γ dist bk /γ LR bk = 0.93) -as we will show below, there are systematic deviations between these numbers which provide insight into the statistics of contact breaking.
We can repeat this procedure for a synthetic ensemble of uncorrelated systems. From the frequentist distribution of contact breaking strains ρ(γ † i j ) of the N = 1024, P = 10 −2 ensemble, we draw N bk = 1147 contacts for each of N s = 1000 systems (bootstrapping). For each system, we calculate the minimum strain γ * . We then compare the mean breaking strain γ bk = γ * = 1.34(4) × 10 −4 to γ dist bk = 1.4 × 10 −4 . Here, we find γ dist bk /γ bk = 1.05 ± 0.04 > 1. Values below 1 thus indicate significant deviations from uncorrelated systems.
Distribution of strains: We now probe the distribution of strains of first contact breaks. Consider an ensemble of M packings of N particles, each with N bk (m) contacts for which we calculate the breaking strains γ † i j . This yields a total of Σ M m=1 N bk (m) ≡ M N bk samples (values of γ † i j ), as illustrated in Fig. 17 for a synthetic data set, as well as for two data sets at fixed P and N. First, we can collect all breaking strains in a distribution ρ(γ † i j ) (black curves in panels b,e,h). As illustrated in Fig. 17 there are now two operations we can perform. Equivalent to what we do to determine γ LR bk in linear response, we can determine the minimum breaking strain for each of the M packings, obtaining M breaking strains (red crosses in panels a,d,g) and the corresponding distribution ρ(γ LR * ,mk ) (shown as red curves in panels b,e,h, as a fraction of ρ(γ † i j )). Alternatively, we may also consider the M smallest values out of M N bk samples taken out of the distribution ρ(γ † i j ) (blue circles), which yields the distribution ρ(γ < ) := ρ(γ|γ ≤ γ dist bk ) (blue curve). The mean values considered above are related to these distributions as follows: γ LR bk is the mean of the ρ(γ LR * ,mk ), whereas γ dist bk is the maximum value of γ < in ρ(γ < ). Clearly, the distributions ρ(γ LR * ,mk ) and ρ(γ < ) in general will be different, but if the different packings are statistically indistinguishable and large enough to allow for self-averaging, so that γ LR bk ≈ γ dist bk , these distributions are directly related (see below), which yields a statistical test on the nature of the contact breaking strains.

Results:
We have determined γ bk and γ dist bk for all (N, P) ensembles. In Fig. 16b we plot N 2 γ dist bk vs N 2 P log 10 (N) −0.7 , and in Fig. 16c we plot the ratio γ dist bk /γ LR bk vs N 2 P log 10 (N) −0.7 . At low N 2 P log 10 (N) −0.7 , we find that γ dist bk and γ LR bk exhibit similar scaling with N 2 P log 10 (N) −0.7 , but that their ratio γ dist bk /γ LR bk ≈ 0.6 < 1.05 ± 0.05 points to deviations from self-averaging. At very high N 2 P log 10 (N) −0.7 , γ dist bk increases faster than γ LR bk and appears to reach equality for the highest values of N 2 P -we suggest that here the packings are large enough to be self-averaging.
To further characterize the origins of this breakdown of self averaging in small systems, we take a closer look at the distributions ρ(γ * ,mk ) and ρ(γ < ) in Figs. 17 and 18. In Fig. 17(ac), we plot each value of γ † i j for the first 100 systems in the synthetic ensemble described above. When we compare the pdfsof the per system ρ(γ * ,mk ) (red curves in panel b) and distribution minima ρ(γ < ) (blue curves in panel b), we note they are similar for small values of γ † i j , but different for larger values of γ † i j . In Fig. 18a we compare the cdfof the per system minima to the cdfof the whole distribution. In the synthetic data, we can deduce that the inverse cdfof minima Pr(γ * ≥ γ) relates to the cdfof the distribution Pr(γ i j < γ) as for large enough N bk for a given N bk Pr(γ i j < γ). In Fig. 18a, we plot Pr(γ * ≥ γ) as a function of N bk Pr(γ i j < γ) for both the synthetic distribution described above, as well as for a synthetic distribution with small N bk . We observe the exponential scaling predicted in Eq. 26 for both. Hence, one expects 63% of the N s per-system minima γ * to be present in the set of N s global minima γ < . In Fig. 17(d-f), we plot each value of γ † i j for the first 100 systems, taken from the N = 1024, P = 10 −2 ensemble. The relation between the pdfsof the per system ρ(γ LR bk ) (red curves in panel e) and distribution minima ρ(γ < ) (blue curves in panel e) are similar to those of the synthetic data, and γ dist bk = 1.4 × 10 −4 and γ LR bk = 1.5 × 10 −4 are quite similar. Consistent with this, a plot of Pr(γ * ≥ γ) as a function of N bk Pr(γ i j < γ) is approximately exponential, although slight deviations can be seen in the tails of these distributions (Fig. 18b).
In Fig. 17(g-i), we plot each value of γ i j for the first 100 systems, taken from the N = 16, P = 10 −6 ensemble. The differences between the pdfsof the per system ρ(γ bk ) (red curves in panel h) and distribution minima ρ < (γ dist bk ) (blue curves in panel h) are more significant, and γ bk = 1.6 × 10 −6 and γ dist bk = 1.1 × 10 −6 are quite distinct. Consistent with this, a plot of Pr(γ * ≥ γ) as a function of N bk Pr(γ i j < γ) deviates significantly from an exponential (Fig. 18b). This deviation points to a lack of self-averaging in small systems.
Interpretation: We now discuss two possible scenarios to explain the deviations for small N 2 P log 10 (N) −0.7 . First, each finite packing could have a different distribution of γ i j , but between packings these distributions are related by an overall scale factor. The data shown in Fig. 17g suggests that this is possible. To understand the effect of such 'overall scale factor' for the statistics, we draw an overall system scale from a uniform distribution U(0, 1) for each of the synthetic systems, and multiply the strains for each system with this scale factor. The resulting behavior is shown in Fig. 18a, where we see the decay is much slower than for uncorrelated systems. The reason for this is that packings with a low minimum will typically come from a system which contains other low strains. This saturates the low strain region of the overall distribution with strains that are not system minima. The data extracted from our direct simulations (Fig. 18b) show a similar decay, slower than exponential, with slower decays for lower pressures. To directly check whether a per-system scale can explain the behavior, we divide all strains by the mean strain for each system, and show the results in Fig. 18c. In the case of a simple scale incorporated in synthetic data, this brings the be-  ( N bk = 1147). For the same ensemble, data with a single value from a distribution with lower mean (dot-dashed blue) and for systems with an overall per-system scale (dashed purple) are also shown. Dotted red: Synthetic data, from ρ(γ † i j ) in the N = 16, P = 10 −6 ensemble ( N bk = 16). The gray line indicates Pr(γ * ≥ γ) = exp(− N bk Pr(γ i j < γ)). (b) Data from our simulations. We observe the curves decay slower than exponential, indicating correlations between contacts. Curves from top to bottom: N = 1024, P = 10 −6 ; N = 16, P = 10 −6 ; N = 16, P = 10 −2 ; N = 1024, P = 10 −2 ; (c) Data from (b), but with all strains rescaled to the mean of strains within one system. This reduces the effect of a per-system scale (dot-dashed red), but does not completely negate it. The behavior for the packing-derived data is unchanged as compared to (b).
havior closer to the simple exponential (dashed purple line). The behavior is still not purely exponential due to the subtle effects we induce with this normalization step. Nevertheless, we note that the rescaling has very little effect on the contact change strains shown in Fig. 18b. We therefore conclude the correlations cannot be simply explained by an overall system scale.
Second, inspired by Lerner et al. [49], we now investigate whether we can recover the behavior of γ LR bk using extremal value statistics by assuming that most contacts are drawn from a distribution with mean k, but a limited number of 'weak' contacts are drawn from a distribution with mean k k. In the case of one extraordinarily weak contact in each packing, we expect most of the k system minima to show up in the lowest k values of the entire set of strains. We have simulated this by dividing one strain in each of the synthetic packings by 10 3 . As we see in Fig. 18a, Pr(γ * ≥ γ) decreases much more rapidly than exponential, and drops to Pr(γ * ≥ γ) = 0 around N bk Pr(γ i j < γ) ≈ 3 -in other words, the k minima are all found in the lowest 3k values of the full set. The exact point of intersection depends on how weak the contact is, and on how many weak samples are in the packing. However, our data for actual packings shows a slower than exponential decay, thus discounting the 'weak contact' hypothesis as source for the correlations in our systems.
Hence, in conclusion: for sufficiently large systems, packings are self averaging, and extremal value statistics may be sufficient to determine the mean value and distribution for the first contact break strains.

VI. CONCLUSION
We have presented a systematic analysis of the first contact changes in soft spheres sheared quasistatically under controlled strain. There are several important conclusions. To begin, contact changes are strongly sensitive to both the system size and the distance to jamming, and finite size corrections play an important role. We find distinctly different scaling relations in the limit N 2 P 1, which is relevant as the system size shrinks or the confinement pressure drops, and in the limit N 2 P 1, which describes thermodynamically large ensembles of soft particles. The characteristic strains describing made and broken contacts can be rationalized via simple mean field-like scaling arguments and Poisson statistics, while log corrections or weak corrections to scaling can improve their accuracy.
Contact changes are also reflected in the mechanical response, including the shear modulus. We have shown that the ensemble-averaged differential shear modulus is unaltered by contact changes in thermodynamically large systems. This finding rationalizes the ability of effective linear response (i.e. Hooke's law) to describe bulk mechanical properties even at finite values of the strain [19,48]. However extreme value analysis suggests that the thermodynamic limit is reached more slowly than one would infer from the system size dependence of contact change strains and mechanical properties.
Finally, we have demonstrated surprising correlations in the spatiotemporal patterning of successive contact changes. These suggest the need for further study of particle scale dynamics on finite strain and time scales. Open questions include the interplay between jamming physics and microscopic irreversibility, as well as the role of viscous interactions. Both can be addressed, e.g., with simulations of oscillatory rheology [48,50].
Appendix A: Creating and shearing a packing Boundary conditions: We use periodic boundaries in a non-square box, where each particle has periodic copies at r = r i + n x · L x + n y · L y , where r i is the canonical position of the particle, n x and n y are integers and L x = (L xx , L xy ) and L y = (L yx , L yy ) describe the box. The area of the unit cell is L 2 = L xx · L yy , the Lees-Edwards shear strain is α = L yx /L and the pure shear strain in δ = L yy −L L = L yy /L xx − 1 . For square cells L x = (L xx , 0) and L y = (0, L yy ), and consequently α = δ = 0 [23,24,51]. In contrast, here we require that the energy is at a minimum with respect to α and δ for the initial condition, which guarantees that we obtain ε + all packings where the shear modulus is positive and the residual shear stresses are zero [2,33], as one expects for a physical system at rest. We keep L xy = 0 as allowed by rotational symmetry.
Interactions, energy and stress: Our system consists of a bi-disperse mix of soft disks with repulsive harmonic interactions, using N/2 small particles with R s = 1 and N/2 large particles with radius R l = 1.4. The interaction between particles is determined by their overlap δ i j = max(0, | r i j | − R i − R j ), where r i j is the center-to-center distance of the two particles: r i j = r i − r j − n x,i j L x − n y,i j L y , where r i and r j are the canonical particle positions, and n x,i j (n y,i j ) is 0 if the closest copy of j to i is the canonical copy, +1 if it is across the right (top) boundary and −1 if it is across the left (bottom) boundary. Contact forces have magnitude f i j = kδ i j , where k = 1 is the spring constant, and result in the harmonic potential U i j = (k/2)δ 2 i j . The internal energy is given by the sum of all inter-particle potentials, U = i, j U i j = i, j (k/2)δ 2 i j . Length scales, stresses and energies are expressed in units R s , k and kR 2 s respectively. The boundary stresses are the simple shear stress σ yx = σ xy , the deviatoric (pure shear) stress τ = 1 2 σ xx − σ yy , and the volumetric stress P int = 1 2 σ xx + σ yy , which are computed using the Born-Huang approximation [52,53] where a, b ∈ {x, y} and the sum is over all particle pairs i, j.
Preparing a packing: To create ε + all packings at given pressure P between 10 −7 and 10 −2 , we minimize the enthalpy H = U + PL 2 . We place our bidisperse N particles within a square box with size L 2 init = φ init N 2 πR 2 s + N 2 πR 2 l , where φ init ≡ 0.8 is chosen to be far below the jamming density φ J ≈ 0.84. We use a combination of the Conjugate Gradient method [54] and the Fast Inertial Relaxation Engine (fire) [55] algorithms. The latter is much faster, but is unstable when the overlaps between particles are large. We therefore initially relax the packing using standard Conjugate Gradient methods to resolve the largest overlaps with fixed boundaries, and minimize the energy until |∆E| ≤ 10 −2 · E. We then use the fire algorithm, allowing the boundaries (i.e., L xx , L yy , and L yx ) to deform, and relax the system until |∆H| ≤ 10 −17 · H, and |σ yx | ≤ 10 −15 .
As we will study changes in individual contacts, and in particular probe the strain at which the first contact change takes place, we anticipate the need to study finite size effects. Moreover, we anticipate that many quantities will rescale with N 2 P as has recently been found in [1,2,37]. We therefore prepared ensembles of sheared systems at a range of N and P. Most ensembles contain 100 systems, with some ensembles containing up to 5000 systems. To characterize the behavior at the first contact change, we created a set of ensembles having N and P on a log-spaced grid, with N = 16, 32, . . . , 1024 and P = 10 −7 , 10 −6 5 6 , . . . Simple shear, contact changes and rattlers: We perform quasistatic shear, so viscous damping is irrelevant and only the elastic interactions between particles are taken into account. We apply shear by distorting the unit cell as L y (γ) = L y (0) + γL ·x, i.e., we change α → α + γ, while keeping L 2 and δ constant. We then use the fire algorithm to relax the system (keeping the boundaries fixed) until |∆H| < 10 −13 · H, where we sacri-fice a small error in the particle positions for simulation speed. We found that this is accurate enough for the detection of contact changes and to determine the stress and energy at the contact change: The details of the relaxation do not influence the detection of contact changes, and the relative error in σ xy is typically less than 10 −6 . Note that in this strained state, we are now no longer in an enthalpy minimum with respect to the boundary conditions, so σ xy 0 and G can become negative.
Appendix B: Calculating the linear response In this appendix we will briefly review how, based on the initial particle positions, box size and box shape, we determine the linear response of the system. Given an applied deformation of the box, we can determine the resulting particle motion, forces and energy cost [2,24,31,33].
The state of the system can be described as a vector where (x n , y n ) is the position of particle n and the four parameters L i j describe the box size and shape. We only include particles that are part of the load bearing network (non-rattlers). We then prescribe a displacement |∆q . We determine the energy in the new state |q + ∆q by expanding U up to second order: is the Jacobian and the extended Hessian at q| [28]. Because the initial state is at an energy minimum, the Jacobian term is zero, and the leading contribution to the energy comes from the extended Hessian. For a given displacement, the energy cost is thus given by and the resulting forces on particles and boundaries by However, typically, we do not know the displacement of each particle. Instead, we wish to calculate the displacement of the particles given a change in the boundaries, i.e., find a state where, given the new boundaries, the sum of forces on each particle is zero. To find this state, we split the extended Hessian into four parts: where the ordinary Hessian H xx describes the particle-particle interactions, H bx the interactions between boundaries and particles, and H bb those between different boundaries. We can then rewrite Eq. B6 as follows: where |∆q x and |∆q b are the displacements of particles and boundaries, and |∆ f x and |∆ f b the corresponding forces. Setting the forces on the particles to zero, we find Solving for |∆q x gives us the particle displacement as a function of the deformation of the simulation box Unfortunately, H −1 xx cannot be calculated due to the two zeroenergy translational modes. Instead, we choose to use the Moore-Penrose pseudoinverse H + xx , which fixes the zeroenergy translational modes in place [56, §6.4]: To calculate the energy cost and the stress on the boundary, we use the full displacement vector and, again using Eq. B8, find The corresponding stress can be calculated as but in practice, it is more convenient to calculate the stress by using the Born-Huang approximation [52,53], on the new particle positions |q x = |q x + |∆q x . The stress also allows us to determine the elastic modulus corresponding to a given boundary deformation For the resulting energy change we use |∆ f x ≡ 0 to find We now have all ingredients in place to calculate, for a given boundary deformation, the particle displacements, stress response and energy change from linear response.  At low N 2 P, σ is independent of pressure and at high N 2 P we recover a scaling to σ ∼ δ 0.25 , consistent with [32]. (b) Same, for σ ⊥ . At low N 2 P, σ is independent of pressure. At high pressure we find a scaling σ ∼ δ −0.2...−0.15 , somewhat slower than the δ −0.25 found in [32].
Appendix C: Finite size scaling of ρ(u ,i j ) and ρ(u ⊥,i j ) In this appendix, we will discuss the the distributions of u ,i j and u ⊥,i j , which provide a continuum description of the interparticle motion. For each particle pair i, j, we split the interparticle velocity u i j = ∂ x i j /∂γ in components parallel and perpendicular to the contact: Using every contact in every packing in an ensemble, we then build the frequentist distributions ρ(u ,i j ) and ρ(u ⊥,i j ).
In the following, we will discuss the relationship between the shape and scale of these distributions and N and P. Earlier work [32] has focused on Hertzian systems at intermediate to high pressure (P 2/3 ∼ δ ≥ 3 · 10 −4 ). They find the shape of the distribution does not depend on P, and find a simple single scaling of the overall scale with P. We extend this with harmonic systems much closer to jamming (P ∼ δ ≥ 10 −7 ). At high pressures, we recover the same behavior, but close to jamming, we find (i) the shape of the distributions depends on the pressure P, and (ii) the widths of the distributions scale with N 2 P, with two distinct scaling regimes.
Shape of distributions: In Fig. 19, we plot the probability density functions of u ,i j and u ⊥,i j , rescaled by their standard deviations σ and σ ⊥ , for ensembles with different system sized and pressures. We note that, even though the different distributions cannot be collapsed with a single scale parameter, the majority of the behavior is captured in the standard deviation σ. For both distributions, we observe the distributions become increasingly peaked near 0, and, although neither pdfdiverges, this peak appears to develop a sharp kink for small pressures. We observe the shape changes with P, and, for large enough N, is largely independent of N -N 2 P is not the relevant scaling parameter here. Surprisingly, this means the overabundant low values are still present for large systems at P ≈ 10 −3 , which would normally not be considered 'close to jamming'.
Scaling of standard deviations: Ellenbroek et al. [32] find the width of the distributions scale as where δ is the mean overlap between pairs of particles in contact in the ensemble. If we assume (i) the standard deviations will scale with N 2 P and (ii) the distributions are independent of N for large N, Eq. C3 and Eq. C4 suggest plotting should collapse our data. We note that, because the shape of the distribution varies, the choice of the scaling parameter (e.g. a percentile rather than the standard deviation) can have a rather large effect on the collapse (which can reach ±0.2 in the scaling exponent), and we therefore do not expect a perfect match.
In Fig. 20a, we find the best scaling collapse for σ is close but not equal to the expected scaling: we find σ ∼ N −0.4 at low N 2 P rather than σ ∼ N −0.5 . Nonetheless, we suggest that the scaling is close enough to be consistent with the proposed scaling. At low pressures, we find that σ only depends on N, and no longer depends on P. For N 2 P 1, we find the expected σ ∼ P 0.25 power law.
For σ ⊥ , we find Eq. C6 provides a rather good collapse (Fig. 20b). At low N 2 P, we find σ ⊥ becomes independent of P, and at high N 2 P, we find behavior similar, but different from the expected σ ⊥ ∼ P −0.25 power law.
Surprisingly, we find both σ ⊥ and σ reach a pressureindependent plateau for low N 2 P. This has important implications for the behavior close to jamming -in contrast to what is generally assumed, σ ⊥ /σ does not diverge for low pressures, but reaches a plateau whose value diverges as σ ⊥ /σ ∼ N 0.9 in the thermodynamic limit.

Appendix D: Discussion
Finally we will discuss our findings in the light of alternative scaling models that have surfaced in the literature. Nonlinearities in jammed packings at finite temperature were studied in Schreck et al. [14], and these authors find a different scaling that we attribute to their averaging over modes. Moreover, Combe and Roux [18] and Lerner et al. [36] have approached the problem from a hard particle perspective, and find a scaling law very close to the behavior we find close to jamming.

Excited eigenmodes
Schreck et al. [14] investigated contact breaking in jammed sphere packings using excited eigenmodes. They displace particles along an eigenmode: where r 0 is the original state, r the excited state, N the system size,ê k the eigenvector for eigenmode k, and δ the excitation amplitude. The system is then allowed to evolve at fixed energy. For small excitations δ, the system oscillates around a base state, and most energy is contained in the initial eigenmode. However, for excitations larger than a critical excitation amplitude δ c (k) there is a sharp increase in how much energy spreads into the other eigenmodes of the system. Schreck et al. find that δ c is directly related to the first contact change in the system. Surprisingly, they find that contacts only break, even for large systems (N = 1920) at high densities (∆φ = 10 −2 ).
For each system, δ c (k) is calculated for every eigenmode k. The authors then measure the average energy where ω k is the eigenfrequency of eigenmode k and the mean is taken over all eigenmodes.
For the scaling of the energy per particle E/N with the density "where A(∆φ) is only weakly dependent on ∆φ and β ≈ 1.7" [14]. Close to jamming (N∆z = 0 . . . 2), they find A(∆φ) is constant and β = 1 . . . 2 [57]. Writing this in terms of E, taking A(∆φ) as constant and using ∆φ ∼ P: To compare this with our results, we note that To test whether this matches the data, we plot N (β+3)/2 γ as a function of N 2 P in Fig. 21a, using the published value β = 1.7. We find, firstly, that the collapse is not very good. Secondly, we find the 0.75 power law for the upper branch overestimates the actual strains. To a lesser extent, the lower branch also deviates from Eq. D7. This is also reflected in the residuals in Fig. 21b -neither branch collapses onto a constant value.
We expect these differences arise due to the averaging in Eq. D2, which means the energy is effectively an average over 2N modes within the same system. In Sec. V, we will see that averaging over all contacts loses many of the features we found for the first contact change.

Hard particle systems
The question of contact breaking and plasticity has also been studied in systems of hard particles. These systems are isostatic [59], which means a contact change will cause the system to unjam, and thus contact changes are directly connected to plastic events. Isostaticity also implies that the force distribution is unique, and can be derived directly from the particle positions [60]. On the other hand, because the systems are isostatic, the results can only describe the N 2 P 1 limit of soft particle systems.
Combe and Roux [18] investigated the prevalence of and distance between strain jumps in a system under uniaxial stress-controlled compression. They found that the spacing between events is described by a exponential distribution in δq(N/1024) 1.16 , where δq is the relative uniaxial stress increment ∆σ/P. This is consistent with modeling contact changes as a Poisson process.
To calculate the scaling of γ bk with N and P in this system, we first note that the mean stress required to break the first contact scales as ∆σ ∼ P ∆q ∼ P/N 1.16 . (D8) We can then calculate the γ bk using the uniaxial compression modulus E. Using that K ∼ 1 and G ∼ 1/N near jamming, E is given by [61] E = 4 and the expected mean strain to break the first contact is thus given by which is very close to the P/N 0.20 scaling we found by fitting our data to a pure power law (Eq. 18). A theoretical argument for this power law, based on the concept of "weak" contacts that connect to local motion, and "strong" contacts that are connected to global motion, was introduced in [36]. Wyart [58] uses this to predict that the strain for the first contact change should scale as which is close to the value found in [18]. In Fig. 21c, we show this scaling also provides a good match to our data -the 0.15 exponent can be seen as a power law correction to our initial γ ∼ P scaling near jamming, and is essentially indistinguishable from either log or 0.2 power law corrections.