Shear hardening in frictionless amorphous solids near the jamming transition

Abstract The jamming transition, generally manifested by a rapid increase of rigidity under compression (i.e. compression hardening), is ubiquitous in amorphous materials. Here we study shear hardening in deeply annealed frictionless packings generated by numerical simulations, reporting critical scalings absent in compression hardening. We demonstrate that hardening is a natural consequence of shear-induced memory destruction. Based on an elasticity theory, we reveal two independent microscopic origins of shear hardening: (i) the increase of the interaction bond number and (ii) the emergence of anisotropy and long-range correlations in the orientations of bonds—the latter highlights the essential difference between compression and shear hardening. Through the establishment of physical laws specific to anisotropy, our work completes the criticality and universality of jamming transition, and the elasticity theory of amorphous solids.


INTRODUCTION
The jamming transition occurs in a wide range of soft materials, ranging from granular matter to colloidal suspensions to glasses.As a non-equilibrium, athermal phase transition, its criticality is specified by a set of scaling laws [1][2][3][4].In particular, the scaling relationship (under the harmonic approximation of the interparticle interaction) G ∼ ∆Z, between the shear modulus G of the jammed phase and the excess coordination number ∆Z, has been derived by the mean-field replica glass theory in infinite dimensions [5] and microscopic elasticity theories [6,7], with a promising agreement to data from isotropic compression simulations [2].Here we demonstrate that this relationship cannot fully account for the shear hardening behavior -the influence of fabric anisotropy (quantified by the macroscopic friction coefficient µ = σ/P , which is the ratio of shear stress σ to pressure P ) is significant.We report numerically a new scaling law unique to shear, µ ∼ σ β with β ≈ 0.25, between µ and σ, and derive theoretically an additional anisotropic contribution to the elasticity, G AI ∼ µ 2 , verified by our simulation data.
Shear hardening, meaning that G increases by shear, is closely related to several other interesting phenomena.The vestige of shear hardening in the unjammed phase is the phenomenon of shear jamming (the onset of rigidity under a fixed volume shear deformation) observed in many experiments [8][9][10][11][12] and simulations [13][14][15][16][17][18][19][20][21][22][23][24][25][26].Shear hardening is generally accompanied by increasing pressure under fixed-volume conditions, which implies that under fixed-pressure conditions a dilatancy effect would occur [27][28][29].Recent studies have shown the possibility of complete decoupling between friction and shear jamming/dilatancy [17,20,21,30].In this study we * yuliangjin@mail.itp.ac.cn show that friction is also not essential for shear hardening.Due to the lack of internal friction, the conventional sawtooth model [31] cannot be applied anymore, and frictionless mechanisms are demanded.To this end, here we provide a phenomenological explanation based on the generalized jamming phase diagram [21], and reveal microscopic origins of frictionless shear hardening from the elasticity theory of amorphous solids [6,7,[32][33][34][35].Our simulation and theoretical results demonstrate that fabric anisotropy makes important contributions to shear hardening, which resembles a similar connection between anisotropy and shear jamming [15].Note that the discussed shear hardening occurs in the solid phase (above jamming), which shall not be confused with shear thickening that presents in non-Newtonian fluids (below jamming) such as dense suspensions [36].
Shear hardening has been directly observed in previous studies [12,22,37,38].Simulations report a hardening regime on the stress-strain curve, following elastic and softening regimes [22,37].These previous studies [22,37] consider mechanically trained packings, which have effectively a moderate degree of annealing.The shear hardening scaling regime in such systems is limited (about one decade) due to the interruption of yielding (note that the yield stress σ Y vanishes in the rapid quench limit [39]), making a reliable estimation of the scaling exponents difficult.To overcome this issue, here we simulate ultra-stable packings annealed by an efficient swap algorithm [40][41][42][43][44].The yield stress of our ultra-stable packings is significantly larger than that of mechanically trained packings [38], and consequently the shear hardening scaling regime is extended up to about eight decades.This advantage allows us to accurately determine shear hardening exponents, which turn out to be independent of the degree of annealing.Plotted are stress-strain curves of (A) a rapidly quenched system near the jamming transition, (B) a deeply annealed system near the jamming transition, and (c) a deeply annealed system over-compressed well above the jamming transition.Data are obtained from constant volume AQS simulations of 3D frictionless soft spheres (SSs).The dashed and solid lines represent averaged and single-realization curves respectively.The inset in panel (B) is an enlarged view of the single-realization curve.

Conditions for shear hardening
As demonstrated by the linear plots in Fig. 1, the stress-strain curves, averaged over samples, typically display three kinds of behavior under constantvolume athermal quasi-static shearing (AQS).(i) Rapidly quenched systems with the minimum jamming density (see Materials and Methods), ϕ j ≈ ϕ J , behave like plastic flows in the zero pressure limit P 0 → 0, where P 0 ≡ P (γ = 0) is the pressure of unstrained configurations: σ(γ) ≈ 0 and G(γ) ≈ 0 (see Fig. 1A).On the stress-strain curve of an individual sample obtained in a single realization of simulation, the coexistence of solidlike (σ > 0 and G > 0) and liquid-like (σ = G = 0) states can be identified [39].(ii) Shear hardening, manifested in the increasing function of G(γ), appears in deeply annealed systems with ϕ j ϕ J , at small P 0 (see Fig. 1B).Note that in this case, there is a narrow shear softening regime (i.e., G(γ) is a descending function) before shear hardening, which can only be identified in log scales (see Fig. 3 and Supporting Information SI Fig. S1).(iii) If the systems are over-compressed well above the jam- Dependence of several characteristic stresses on P0 and ϕj.(A) The yielding stress σY, the onset stress σ h of shear hardening, and the onset stress σs of shear softening are plotted as functions of the unstrained pressure P0, with a fixed ϕj = 0.69.The dashed line is the fitting curve σY(P0) = 1.48 × 10 −3 + 0.317P0.The solid and dotted lines are the relations σ h (P0) = 2.0 × P 4/3 0 and σs(P0) = 0.02 × P 5/4 0 (see Fig. 4).The two lines σY(P0) and σ h (P0) intersect at P * 0 = 1.13 × 10 −2 , above which shear hardening is non-observable.(B) The ϕj-dependence of σY and σ h , with a fixed P0 = 10 −6 .The dashed line is the fitting curve σY(ϕj) = 0.0175(ϕj − ϕJ) + 0.736(ϕj − ϕJ) 2 , where ϕJ ≈ 0.655.It is known that σY vanishes linearly with ϕj − ϕJ in the vicinity of ϕJ [21]; here a quadratic correction is added to account for the deviation from this linear critical scaling in the large ϕj regime.
ming transition, then shear hardening disappears and only shear softening is observed independent of ϕ j (see Fig. 1C).
From the above analysis, one finds that shear hardening emerges near the jamming transition (P 0 → 0) in deeply annealed packings (ϕ j ϕ J ), under constantvolume AQS.These two conditions are demonstrated more quantitatively in Fig. 2. The shear hardening behavior can be only observed in a stress window σ h < σ < σ Y , where σ h is the onset stress of hardening (see below for how σ h is determined).For a fixed ϕ j with increasing P 0 , σ h increases faster than σ Y , as shown in Fig. 2A.Above P * 0 , where σ h (P * 0 ) ≈ σ Y (P * 0 ), the shear hardening regime disappears (The finite-size analysis is in SI, see Fig. S5).For a fixed P 0 , σ h is nearly independent of ϕ j (see Fig. 2B), while σ Y decreases with decreasing ϕ j and vanishes when ϕ j → ϕ J [39,45] (a more detailed study on σ Y (ϕ j ) can be found in Ref. [21]).Thus the scaling regime of shear hardening also vanishes as ϕ j → ϕ J .The finite-size effects are discussed in SI Sec. 5.In this study, we do not consider other factors, such as finite tempera-G ∼ σ 3. Shear hardening scalings.Simulation results in 3D of (A) the shear modulus G, (B) the excess coordination number ∆Z and (C) the macroscopic friction coefficient µ as functions of stress σ, for a few different P0.The solid lines represent scaling laws in the shear hardening regime, G ∼ σ 2/5 , ∆Z ∼ σ 2/5 and µ ∼ σ 1/4 .The crossovers at σs and σ h , which are determined in Fig. 4, are marked by black solid squares and pentagons respectively.(D) Universal relationship between µ and σ in 3D samples prepared by cyclic shear (ϕj = 0.662), cyclic compression (ϕj = 0.663), and swap thermal annealing (ϕj = 0.672, 0.677, 0.682, 0.690), as well as in two and four dimensions (P0 = 10 −6 in all cases).

Critical scaling laws
Below we focus on deeply annealed samples with a large jamming density ϕ j = 0.69, which are quenched from dense equilibrium liquids at ϕ eq = 0.643.Based on the log-log plots in Fig. 3, we identify consecutively (before yielding) the elastic, shear softening and shear hardening regimes, separated by two crossover points at σ s and σ h .
The data of shear modulus in the elastic and shear softening regimes can be described by a scaling function G(P 0 , σ)/G 0 = G(σ/σ s ), where G 0 = G(P 0 , 0) is the unstrained shear modulus, G(x → 0) = 1, G(x → ∞) = x −2 and σ s ∼ P 5/4 0 (see Fig. 4A).The scaling function G(x) was initially proposed in Ref. [47] to describe both linear and softening behavior in rapidly quenched packings.We find that G(x) agrees well with our data of deeply an-nealed systems in the linear and softening regimes, but deviate from the data starting from σ = σ h which is the onset point of shear hardening.The scaling function G(x) allows us to determine numerically the crossover stress σ s = cP 5/4 0 , where the prefactor c is chosen to be 0.02.The σ s estimated in this way (squares) reasonably separates linear and softening behavior, as can be seen in Fig. 3A.From G(x), the scaling behavior of linear and softening regimes can be extracted.In the linear regime (σ σ s or x = σ/σ s → 0), G(x → 0) = 1, which means that the shear modulus is simply a constant G = G 0 .In the softening regime (σ (where we have used the well know scaling G 0 ∼ P 1/2 0 for un-strained systems [2]).These results are fully consistent with the findings in Ref. [47].Our data show that the presence of a strong shear hardening effect does not change the scalings in linear and softening regimes.It is of interest to provide a theoretical explanation of the softening exponent

studies.
The deviation of the simulation data from the master curve G(x) in Fig. 4 is followed by the shear hardening behavior.The scaling of σ h is determined based on the rescaled plot of macroscopic friction µ in Fig. 4B, with a scaling function µ(P 0 , σ)/µ h = U(σ/σ h ), where 4B shows that elastic and softening data follow a universal scaling µ ∼ σ.Previously, Ref. [22] reports two scalings µ ∼ γ and µ ∼ γ 1/2 respectively in the elastic and softening regimes for a twodimensional model, which, together with the shear softening scaling σ ∼ γ 1/2 , give a consistent result µ ∼ σ.
In the shear hardening regime, σ h < σ < σ Y , the system's behavior is characterized by two scaling laws.Numerical fitting gives µ ∼ σ β with β = 0.248 ± 0.006 and ∆Z ∼ σ ν with ν = 0.411 ± 0.005 in the shear hardening regime (Fig. 3).The scaling of ∆Z is consistent with the ansatz ∆Z ∼ σ 2/5 proposed in Ref. [4].Note that this scaling is examined in [4] by looking at the stress variance σ 2 of isotropically compressed configurations with ϕ j = ϕ J , where the mean stress is zero.Here we provide direct verification of the ansatz, thanks to the emergence of shear hardening in deeply annealed packings (ϕ j ϕ J ).The scaling µ ∼ σ 1/4 (or equivalently P ∼ σ 3/4 ), is absent in the framework of Ref. [4], and to our knowledge, has never been reported previously.Based on the elasticity theory (see below for details), we derive that G is a linear combination of G I ∼ ∆Z and G AI ∼ µ 2 , which, together with the above two scalings, gives G ∼ σ 2/5 in the jamming limit where σ → 0 (see Fig. 3A for a comparison to the numerical data).While it is straightforward to translate these scalings, G ∼ ∆Z ∼ σ 2/5 and µ ∼ σ 1/4 , from the stress-dependent forms into strain-dependent forms based on the relation G = dσ/dγ, one should be cautious about the influence of plasticity in the measurement of stress and modulus [47].
It can be shown that the shear hardening scaling exponents are independent of the preparation method (thermal annealing or mechanical training), the degree of annealing represented by ϕ j , and the dimensionality in d = 2, 3, 4 dimensions (see SI Fig. S2, S3 and S4).In particular, the data of µ versus σ collapse remarkably without any shifting (see Fig. 3D).

Shear hardening as a consequence of memory loss by shear
To understand the origin of shear hardening, it is useful to examine the evolution of state on the generalized jamming phase diagram [21].The key point is that, for a general ϕ j , the isotropic unstrained packings and the asymptotic stationary states at large strains satisfy different sets of equations of states (EOSs); they coincide only when ϕ j → ϕ J as in the initial proposal of the phase diagram by Liu and Nagel [45].The pressure, excess coordination number and macroscopic friction of unstrained states are described by the following EOSs respectively, , and µ 0 = 0.The corresponding EOSs of stationary states are , where µ c ≈ 0.1 [21,48,49].
Under constant volume shear, it is easy to see from Fig. 5 that the system has to gain rigidity, because both the pressure and the number of contacts increase.The increase of anisotropy µ makes an additional contribution to the shear modulus, as demonstrated below.Shear hardening is thus a natural consequence of the memory erasing process caused by shear.Different unstrained states must evolve to the same stationary state that is independent of the initial condition.Initial states that are dense packings (ϕ j ϕ J ) created by deep annealing should also be brought back by shear to the generic pack- ings at ϕ J (corresponding to the rapid quench limit [2]) -this process is accompanied by the increase of pressure and shear modulus under the constant volume condition.

Microscopic origins of shear hardening
Previous studies have developed an elasticity theory for athermal disordered solids near the jamming transition, attributing compression hardening to the increase of interaction bonds via the relation [6,7,35], where c I = 1/30 (we have set both the bond stiffness κ and the unstressed bond length r 0 = D to one).While this theoretical result agrees well with the simulation data of sphere packings under isotropic compression [7], it cannot fully account for our shear hardening data (see Fig. 6A).
To understand the microscopic origin of this deviation, we consider a generalized formula [6] (see SI for details), for P → 0, where depends on the normalized forces fb and bond orientations (represented by unit vectors n b ), N b is the number of contacts and N NR is the number of non-rattler particles (rattlers have fewer than d + 1 contacts in d dimensions).In isotropic packings, the lack of longrange spatial correlations between bond orientations, under the effective media approximation [7], and all contact forces fb are positive).In sheared packings, the bonds tend to align in shear-preferred directions and are orientationally long-range correlated [16] , where r b is the position of bond b, decays to a finite value at large r (see Fig. 7).Consequently, G AI is non-zero in such a case.
Our theoretical analysis predicts a simple relation between G AI and the macroscopic friction µ, where with two constants c 0 and α.In the jamming limit (∆Z ≈ 0), the coefficient c 0 can be evaluated from the formula, c 0 = f 2 3 f 2 , where f and f 2 are the first two moments of the force distribution.In our system, we obtain from simulations f 2 f 2 ≈ 0.5, which gives c 0 ≈ 0.17 (see Fig. S6).Above jamming (∆Z > 0), G AI is lowered by a higher-order correction term ∼ ∆Zµ 2 .
To examine Eq. ( 3), we separate isotropic and anisotropic moduli in our simulation data as follows.Considering that ∆Z(γ) is increased during shear, at each γ we first decompress the configurations keeping γ unchanged, until ∆Z reaches ∆Z 0 ≡ ∆Z(γ = 0) that is fixed by the initial condition ∆Z 0 ∼ P 1/2 0 .The decompressed configurations at different γ thus now have the same ∆Z = ∆Z 0 , but different µ (we denote these configurations as a constant-Z ensemble).In this way, we remove the contributions from the change of ∆Z, left with G AI that changes with µ. Figure 6B shows that G AI is proportional to µ 2 for any constant ∆Z, as predicted by Eq. ( 3).The behavior of the slope c sim AI (∆Z) obtained from fitting the simulation data is also consistent with our theoretical prediction Eq. ( 4).The remaining shear modulus G I = G − G AI linearly depends on ∆Z (see Fig. 6A), with a P 0 -independent slope in agreement with the analytical result c I = 1/30 as in the case of pure compression [7].Note that here ∆Z is increased by shear instead of compression.We find in our simulations that the pressure P does not depend on µ, but solely on ∆Z.Thus the constant-Z ensemble considered here is in fact equivalent to a constant-P ensemble.
Our results show that, within the first-order approximation, the excess coordination number and anisotropy make additive contributions to the shear modulus.At a large strain, the contact network becomes strongly anisotropic, in which case the anisotropic correction to G is important.On the other hand, in the limit σ → 0 (or γ → 0), the anisotropic contribution is subleading: combining the scalings ∆Z ∼ σ 2/5 and µ ∼ σ 1/4 (Figs.3B-C) with Eqs.(1, 2, 3), one obtains G ∼ ∆Z ∼ σ 2/5 (Fig. 3A).
For comparison, we also derive a theoretical expression of the bulk modulus B near the jamming transition (see SI Sec.6), where the constant c 0 = f 2 3 f 2 is identical to the one appeared in Eq. ( 4) for the shear modulus.In the jamming limit ∆Z → 0, the bulk modulus is a constant B(∆Z = 0) = c 0 .In other words, B changes discontinuously at either compression or shear jamming, following a scaling B ∼ ∆Z 0 as previously reported [2,16].Equation ( 5) further gives the next order correction that depends linearly on ∆Z.Plugging c 0 ≈ 0.17 into Eq.( 5) gives B ≈ 0.17 + 0.084∆Z, which is close to the result B = 0.146 + 0.0897∆Z obtained from fitting the simulation data (see Fig. 8A).As expected and suggested by Eq. ( 5), the bulk modulus B is independent of the anisotropy µ (see Fig. 8B for the numerical verification).
In our theoretical analysis, the internal stresses is ignored.Because the considered interaction is purely repulsive, the non-zero internal stress should lower the energy and therefore the modulus.In other words, the modulus obtained from the current theory overestimates the actual value.Approaching the jamming transition where the internal stress vanishes, the current theoretical results are very close to the simulation data, suggesting that the effects of internal stress are not significant in this limit (see Figs. 6 and 8).
Interestingly, a previous study [50] provided an approximated expression for the bulk modulus using the so-called Reuss estimate [51].In Fig. 8, we compare B Reuss with our theory Eq. ( 5) and  5) and Reuss estimation with c0 ≈ 0.17.(B) Bulk modulus B of samples with a varying µ.Each curve has a constant ∆Z = ∆Z0, which is realized by decompression as in Fig. 6 (see the legend in (A) for the corresponding P0 values).simulation data.To make a further comparison, we expand B Reuss to the linear order of ∆Z, using Z ≈ 6+∆Z, which results in B Reuss ≈ c 0 + c0 6 ∆Z.Compared to our theoretical result Eq. ( 5), this final expression gives the same constant c 0 when ∆Z → 0, but a different coefficient of the linear term.

DISCUSSION
Our finding raises several challenges to the theoretical understanding of jamming transition and rheology of amorphous solids.In particular, the scaling µ ∼ σ β with β ≈ 0.25, which is universal in two to four dimensions (see SI), requires an quantitative theoretical explanation.The scaling ansatz proposed in Ref. [4], , N ∆Z , treats ∆Z as the only relevant order parameter for the jamming transition.Equation (2) suggests that the ansatz needs to be extended, since the anisotropic part G AI , which relies on µ instead of ∆Z, contributes independently to the total shear modulus.
The response of amorphous solids to compression and shear has been studied by a first-principle mean-field the-ory -the state-following glass theory [52,53].In the stable-glass phase, the theoretical outcome is σ ∼ Gγ and P ∼ P 0 + Rγ 2 , where R is the dilatancy parameter (note that the second relation reflects the symmetry P (γ) = P (−γ)).These relationships give a scaling µ ∼ σ in the limit γ → 0, which can only explain the first scaling (in the elastic and softening regimes) in Fig. 4

(B).
The issue is twofold: on the one hand, it remains technically difficult to perform the state-following computation in the marginally stable phase where jammed states belong to [53]; on the other hand, conceptually one should consider the possibility of "state rejuvenation" since simulations suggest ϕ j decreasing towards ϕ J in the shear hardening regime [22,38].
Finally, it would be interesting to examine the effect of friction on shear hardening.Recent progress in producing ultra-stable granular materials makes an experimental investigation of frictional shear hardening feasible [12,54].It is expected that the critical exponents are non-universal with different friction coefficients [54].In addition, the correction to the isotropic shear modulus should be even more significant, because friction amplifies spatial heterogeneity and orientational anisotropy of the contact network.

Simulation models
We consider a 3D granular model of frictionless polydisperse spheres with a continuous diameter distribution P (D) ∼ D −3 , for D min ≤ D ≤ D min /0.45 [20,43,44].The particles interact via a harmonic SS pair potential, , where r ij is the inter-particle distance and D ij = (D i + D j )/2 the mean diameter of particles i and j.
We set the average diameter D as the unit of length.The lowest jamming density, or the J-point density [2], of this model is ϕ J ≈ 0.655 [20,55].The reported results are obtained from systems of N = 2000 particles, averaged over 192 independent samples.In addition, we study 2D and 4D models with the same kind of interaction in SI.

Preparation of initial states
Two approaches are used to prepare ultra-stable jammed configurations.(i) Swap thermal annealing.The thermal annealing is realized by applying an efficient swap Monte Carlo algorithm [43,44], which generates equilibrium hard-sphere (HS) configurations at ϕ eq above the mode-coupling theory (MCT) crossover density ϕ MCT ≈ 0.594 [56].After that, we set the temperature to zero in the remaining simulations, and switch to the SS potential v ij (r ij ).By applying athermal quasi-static compressions (AQCs) with the FIRE energy minimization algorithm [57] as described in [21,55,58], we obtain jammed configurations at ϕ eq -dependent jamming densities ϕ j [20,55,59].In the main text, we focus on the case of ϕ j = 0.69 unless otherwise specified, corresponding to ϕ eq = 0.643.(ii) Mechanical training.The method is realized by cyclic AQS [21][22][23]60] or cyclic AQC [13].Random initial configurations are rapidly compressed over ϕ J , and become unjammed after a sufficient number of cyclic AQS or AQC under the constant volume condition [21], meaning that the jamming density is increased to ϕ j > ϕ J .In both approaches, ϕ j correlates to the stability of the jammed states, and reflects the degree of annealing.
The above procedures generate packings at ϕ = ϕ j and P = 6 × 10 −9 , which are extremely close to the jamming/unjamming transition.They are then compressed quasi-statically to over-jammed states at a series of pressures P 0 ≡ P (γ = 0).The pressure P 0 characterizes the distance to the jamming point, via the well-known scaling between P 0 and the density ϕ 0 for harmonic SSs [2], P 0 ∼ (ϕ 0 − ϕ j ).In order to remove residual stresses, the shear stabilization method [61] is applied, and these configurations with zero residual stresses are referred to as the unstrained states (σ = 0).According to this definition, in finite-size unstrained states, the stress σ is always zero, while the strain γ can fluctuate around zero from sample to sample.In our scaling analysis (see Fig. 3), it is better to use σ than γ as the independent variable, since the latter can introduce additional uncertainties.Of course, in the limit of large sample size and/or of larger number of samples, this sample-to-sample fluctuation disappears.In fact, a similar issue exists in the case of compression jamming: the jamming transition is defined at the zero pressure, while the jamming density can fluctuate around the mean value in finite-size simulations.Thus practically it is better to use the pressure than the density as the independent variable in scaling analyses.

Compression and shear protocols
Compression and decompression procedures are realized by rescaling particle diameters with energy minimization.Starting from the unstrained states, we apply constant-volume, simple AQS in the x-z plane [20,21].At each step, the shear strain γ is increased by δγ, followed by the FIRE energy minimization, where δγ is logarithmically increased from 10 −8 to 2 × 10 −4 .
The stress and pressure are calculated using the virial formula [62] where r b and f b are the contact vector and force on bond b, N NR the number of non-rattler particles (rattlers have fewer than d + 1 contacts in d dimensions), and N b the total number of bonds.The shear modulus is defined as G = δσ/δγ + P , averaged over 192 independent samples, where δγ is between 10 −10 and 10 −8 .The term δσ/δγ corresponds to the slope of a single-realization stress-strain curve, which is piecewise linear [63] (see the inset of Fig. 1B).The correction term P , which is typically three orders smaller than the first term, is added so that G is equivalent to the xzxz component of the stiffness matrix.We have checked that the value of G obtained in this way is identical to that directly calculated using the explicit expressions derived from the elastic theory [34].The bulk modulus is computed by the formula B = − δP/δ V − 1 3 P , where V is the volumetric strain and δ V ≤ 2 × 10 −9 .
The changes of two independent parameters Z and µ reflect local rearrangements of particles.The increase in the number of contacts is quantified by the excess coordination number, ∆Z = Z−Z iso , where Z iso = 2d−2d/N NR is the coordination number required to satisfy the isostatic condition at the jamming transition in d dimensions with a finite size correction [64,65], and rattlers are excluded in the computation of Z.The macroscopic friction coefficient µ = σ/P measures the degree of anisotropy caused by shear [66].
In our simulations, the stress σ has much smaller fluctuations than the strain γ (see Materials and Methods).Thus σ is chosen as the independent variable in our scaling analyses (see Fig. 3).For completeness, in Fig. S1, we plot G, ∆Z, µ and P as functions of γ.The data are essentially equivalent to those presented in Fig. 3, with γ treated as the independent variable.

S2. INDEPENDENCE OF SCALING EXPONENTS ON THE DEGREE OF ANNEALING
In the main text, we show that shear hardening presents in deeply annealed packings (ϕ j = 0.69).Here, we perform additional simulations on systems with different ϕ j (i.e., different degrees of annealing).The data in Fig. S2 show that the exponents in all three scalings, G ∼ σ 2/5 , ∆Z ∼ σ 2/5 and µ ∼ σ 1/4 , are independent of ϕ j .Although the curves of different ϕ j seem to collapse in the figure completely, one should note that the yielding stress vanishes (σ Y → 0) in the limit of ϕ j → ϕ J [38], and consequently so does the scaling regime of shear hardening.

S3. SHEAR HARDENING IN SYSTEMS PREPARED BY MECHANICAL TRAINING
Dense packings with a large ϕ j can be generated by either thermal or mechanical annealing.The latter refers to athermal training procedures such as cyclic shear or cyclic compression.The simulation data presented in the main text are obtained from well-annealed systems prepared by swap thermal annealing.Here, we present the results of samples annealed by mechanical training, and show that the scaling exponents are independent of the annealing protocol.
Two different training methods, cyclic shear and cyclic compression, are used.In the cyclic shear protocol [21], random configurations are firstly over-compressed to a density ϕ x > ϕ J .Then, multiple cycles of shear, 0 → γ max → −γ max → 0, are applied until the system becomes unjammed at ϕ x .In the cyclic compression protocol [13], the random configurations are overcompressed to a pressure P max , and then are decompressed to ϕ = 0.62 that is lower than the jamming density.Such cycles are stopped when the desired jamming density ϕ j = 0.663 is reached.Both mechanical training procedures are performed under the athermal quasi-static condition.The strain increment is δγ = ±10 −3 and the density increment is δϕ = ±10 −4 .
The jamming density ϕ j depends on the parameters γ max and ϕ x in the cyclic shear protocol, and P max in the cyclic compression protocol [13,21].In this study, we set γ max = 0.05 and ϕ x = 0.66 in the former, and P max = 10 −2 in the latter.With these parameters, we obtain ϕ j = 0.662 and ϕ j = 0.663 respectively, both of which are above the J-point density ϕ J ≈ 0.655.
These mechanically trained systems also show shear hardening behavior (Fig. S3), with scalings identical to those in Fig. 3, which are obtained from systems generated by swap thermal annealing.We thus conclude that the scaling laws of shear hardening are independent of the way of annealing.

S4. SHEAR HARDENING IN TWO AND FOUR DIMENSIONS
The systems are composed by N = 4000 (2D) and N = 2000 (4D) particles.The diameter distribution is p(D) ∼ D −d , for D min ≤ D ≤ D max /0.45.As in the 3D model, two particles interact with each other via a harmonic soft sphere potential v ij (r ij ), if they are in contact.The mean diameter of all particles is set as the unit length.All particles have the same unit mass.
In 2D, we use the swap thermal annealing protocol to prepare equilibrium hard disk liquids of density ϕ eq = 0.85.After that, we switch to the soft potential, and compress the system to a jammed state of ϕ j = 0.881 (the J-point density is ϕ J = 0.841).
We find that it is very easy to observe shear hardening in 4D.Deep annealing is unnecessary-it is enough to perform only one compression-decompression cycle.The jamming density of the obtained samples is ϕ j = 0.486, which is only slightly larger than the the J-point density ϕ J = 0.480.
The scaling behavior in 2D and 4D is reported in Fig. S4, together with the data obtained in 3D.The scaling exponents in the shear hardening regime are nearly independent of the dimensionality d.This observation is consistent with the previous conjecture of upper critical dimension d u = 2 for the jamming transition [64,67].

S5. FINITE-SIZE EFFECTS
To check finite size effects, we consider systems of different N , with the jamming density ϕ j = 0.69 and the unstrained initial pressure P 0 = 10 −6 fixed.The stressstrain curves with different N nearly coincide with each other before yielding (see Fig. S5A).The yielding stress σ Y dependents on N , satisfying σ Y (N ) − σ Y (N = ∞) ∼ N −0.5 (see the inset of Fig. S5A).In Fig. S5B, we plot µ as a function of σ, and find that there is no noticeable finite size effects on the crossover stress σ h .Thus the shear-hardening regime σ h < σ < σ Y becomes smaller in larger systems, but remains finite in the thermodynamic limit.

S6. ELASTICITY THEORY A. Setup
Model.The model consists of N soft spheres (rattlers are not considered in this theoretical model) interacting via a harmonic potential at zero temperature.The total energy is where (we set both the spring constant and particle diameter to be one) The vector r ij = r j − r i goes from particle i to j.We are interested in the response of the model to athermal quasistatic shear, meaning that the system is force balanced at each shear step.Without loss of generality, we only compute the xzxz component G = C xz,xz of the stiffness tensor C as the shear modulus.
The following approximations are used, which are nonessential for the problem under consideration.(i) In simulations, the interaction is purely repulsive, and therefore Eq. (S2) is truncated for r ij > 1.The theoretical potential is perfectly symmetric.Thus the amorphous solid is modeled by a spring network.(ii) The polydispersity is neglected.(iii) The pre-stress is ignored in the calculation of shear modulus.
Angular averages.The orientation of any interaction bond can be denoted by a unit vector in three dimensions, n = {sin(θ) cos(φ), cos(θ), sin(θ) sin(φ)}.For a simple shear deformation in the x-z plane, the bond angles are G ∼ σ Here the contact b is formed between particles m and l, and b determines the row index.The column index is determined by i and α.Plugging Eq. (S9) into Eq.(S8), one obtains, whose geometric interpretation is transparent: it projects the relative displacement between i and j on to the direction of bond.On the other hand, T (a 3N by N b matrix) converts forces f along bonds (a vector of size N b ) into the net forces F on particles (a vector of size 3N ): or equivalently, where c is the contact matrix (c ij = 1 if particles i and j are in contact, otherwise c ij = 0).The matrix T is the transpose of S, and their product is the Hessian matrix, where H is defined as, B. Decomposition of the shear modulus into affine and non-affine parts A microscopic elasticity theory for amorphous solids has been recently developed [7,32,33,35,69], which decomposes the shear modulus into affine and non-affine parts, The first term is the seminal Born-Huang formula derived for lattices [70], considering only affine displacements of particles.Here N b = N Z/2 and n b is the unit vector along the bond.The term G NA , which does not appear in lattices, originates from corrections caused by non-affine displacements that are necessary to satisfy the force balance conditions during shear, Here Ξ is the affine force (also called mismatch force) field acting on N particles; v k and λ k are respectively the k-th normalized eigenvector and eigenvalue of the Hessian matrix H.Both Ξ and v k are vectors of size 3N , and there are in total 3N eigenvectors of the Hessian matrix.The affine force Ξ can be computed explicitly for the given energy Eqs.(S1) and (S2): where is the strain matrix.

C. Decomposition of the shear modulus into isotropic and anisotropic parts
A formally different expression of shear modulus is derived by Wyart [6], based on the duality between force propagation and soft modes.The total shear modulus is decomposed into isotropic and anisotropic parts, where and Here fp is the p-th state of self stress, a normalized vector representing the set of interaction forces (along the bond directions) that can satisfy the force balance conditions on every particle.The 1 2 N ∆Z vectors of fp are orthogonal to each other.Among them, only f1 corresponds to the real forces f generated in simulations, which have to be positive for all contacts: The rest of the vectors fp with p ≥ 2 are "virtual" forces, having a roughly equal number of positive and negative components [6], i.e., where we have used the orthonormality of fp and Eq.(S7).Note that Eq. (S25) is independent of bond orientations (or the anisotropy parameter R A ), thus capturing only isotropic contributions.Indeed, Eq. (S25) agrees well with the shear modulus measured in isotropically compressed packings [7], but not in sheared packings (see Fig. 6A).In Eq. (S22), the term of positive real forces (p = 1) dominates; the positive and negative virtual forces are roughly equally probable and therefore the terms of virtual forces (p ≥ 2) make minor contributions.The anisotropic shear modulus can then be approximately expressed as, where we have used Eq.(S23).

D. Equivalence of two decompositions
Here we prove the equivalence between the above two decompositions, by deriving Eqs.(S21) and (S22) from Eqs. (S17) and (S18).Our derivation is based on the formalism developed in Ref. [6] The goal is to replace the set of eigenvectors v k , where k = 1, 2, . . .3N , in Eq. (S18) by the set of equilibrium forces fp that balance all particles, where p = 1, 2, . . . 1  2 N ∆Z.This can be done by imposing the conditions of mechanical equilibrium, which connects displacement fields to equilibrium forces.
The first step is to project the particle-based displacement fields v k onto the directions of bonds, by applying Eq. (S8): where the pre-factor 1 √ λ k ensures the normalization.Plugging Eq. (S9) into the above equation gives, Using this expression and Eq.(S19), we can rewrite Eq. (S18) as, Second, based on the equilibrium condition T fp = 0, it can be shown that, for any k and p, v k and fp have to be perpendicular to each other.This is because, where we have used Eqs.(S13) and (S27).Equation (S30) guarantees the minimization of energy at equilibrium.If we further require that fp 's are perpendicular to each other, then v k and fp form a complete orthonormal basis of the vector space of dimension N b (note that in the jammed phase, N b = 3N + 1 2 N ∆Z), which can be represented by an N b by N b orthonormal matrix, Because the columns of an orthonormal matrix are also orthonormal, we get, Finally, combining Eq. (S32) with Eqs.(S16), (S17) and (S29), we obtain Eqs.(S21) and (S22).

E. Anisotropic shear modulus
Next, we relate the anisotropic shear modulus Eq. (S26) to the macroscopic friction coefficient µ = σ/P .Using the virial expressions of σ and P (see Materials and Methods), we obtain, Note that the spatial correlation between forces is short-ranged [2], which means that the first term Our simulation results show that the force distribution p(f ) is independent of µ (see Fig. S6).From p(f ), it is easy to evaluate f 2 / f 2 ≈ 0.50.In addition, Z ≈ 2d = 6 near the jamming transition.With these values, when ∆Z → 0 we have G AI ≈ 0.17µ 2 .
It is also possible to consider the contributions from virtual forces fp for p ≥ 2. In analogy with Eq. (S33), we obtain,  f 2 and α is a constant that can be determined from the fit of simulation data (see the inset of Fig. 6B).Equation (S36) shows that the contribution of virtual forces is higher-order, which only appears in systems above jamming (∆Z > 0).

F. Bulk modulus
It is straightforward to generalize the theoretical analysis to other components of the stiffness tensor, C αβ,γδ .Similar to the above derivation, any component C αβ,γδ can be decomposed into isotropic and anisotropic parts.The isotropic part is where c αβ,γδ (∆Z = 0) = c 0 and the sub-leading correction is proportional to ∆Z.In particular, the bulk modulus B is, which is independent of the anisotropy.

FIG. 1 .
FIG.1.Stress-strain curves.Plotted are stress-strain curves of (A) a rapidly quenched system near the jamming transition, (B) a deeply annealed system near the jamming transition, and (c) a deeply annealed system over-compressed well above the jamming transition.Data are obtained from constant volume AQS simulations of 3D frictionless soft spheres (SSs).The dashed and solid lines represent averaged and single-realization curves respectively.The inset in panel (B) is an enlarged view of the single-realization curve.

FIG. 5 .
FIG. 5. Equations of states.Plotted are equations of isotropically jammed (γ = 0, solid lines) and steady states (γ = ∞, dashed lines), where the (A) pressure P , (B) excess coordination number ∆Z, and (C) macroscopic friction coefficient µ are plotted as functions of volume fraction ϕ.The vertical arrows present the constant-volume shear protocol.

7 P 4 P 2 G 1 ∆Z 0 = 4 . 1 × 10 FIG. 6 .
FIG.6.Decomposition of the shear modulus into isotropic and anisotropic parts.(A) Total shear modulus G (lines) and its isotropic part GI (points) as functions of ∆Z.The bold dashed line corresponds to the theoretical result GI = cI∆Z with cI = 1/30[7,35].(B) Anisotropic part of the shear modulus GAI as a function of µ 2 for a few fixed ∆Z = ∆Z0 (see the legend in (A) for corresponding P0 values).The bold dashed line indicts our theoretical result GAI = c0µ 2 for ∆Z → 0, with estimated c0 ≈ 0.17.The solid lines represent linear fitting to the data.The slopes are plotted in the inset, which are fitted to a linear function according to Eq. (4), giving c sim AI (∆Z) = 0.15 − 0.11∆Z (line).
FIG. S1.Simulation results of (A) the shear modulus G, (B) the excess coordination number ∆Z, (C) the macroscopic friction coefficient µ and (D) the pressure P as functions of strain γ.The crossovers γs and γ h are represented by solid squares and pentagons, respectively.

) σ Y
FIG. S5. (A)Stress-strain curves and (B) macroscopic friction coefficient µ as a function of stress σ for different system sizes N , with fixed P0 = 10 −6 and ϕj = 0.69.The system size dependence of yielding stress σY (i.e., the maximum stress) is shown in the inset of (A) and the solid line is the fitting curve σY(N ) = 0.0107N −0.5 + 0.00124.

p ≥ 2 .n x b n z b n x b n z b ≈
Under the effective media approximation, i.e., assuming independence between forces and bond orientations, fp,b fp,b fp,b fp,b n x b n z b n x b n z b , Eq. (S21) becomes fp,b2 = δ b1b2 .(S32)
The minus sign in Eq. (S35) comes from the property of virtual forces Eq. (S24), according to which, (