Memory of jamming–multiscale models for soft and granular matter

Soft, disordered, micro-structured materials are ubiquitous in nature and industry, and are different from ordinary fluids or solids, with unusual, interesting static and flow properties. The transition from fluid to solid—at the so-called jamming density—features a multitude of complex mechanisms, but there is no unified theoretical framework that explains them all. In this study, a simple yet quantitative and predictive model is presented, which allows for a changing jamming density, encompassing the memory of the deformation history and explaining a multitude of phenomena at and around jamming. The jamming density, now introduced as a new state-variable, changes due to the deformation history and relates the system’s macroscopic response to its micro-structure. The packing efficiency can increase logarithmically slow under gentle repeated (isotropic) compression, leading to an increase of the jamming density. In contrast, shear deformations cause anisotropy, changing the packing efficiency exponentially fast with either dilatancy or compactancy as result. The memory of the system near jamming can be explained by a micro-statistical model that involves a multiscale, fractal energy landscape and links the microscopic particle picture to the macroscopic continuum description, providing a unified explanation for the qualitatively different flow-behavior for different deformation modes. To complement our work, a recipe to extract the history-dependent jamming density from experimentally accessible data is proposed, and alternative state-variables are compared. The proposed simple macroscopic constitutive model is calibrated from particles simulation data, with the variable jamming density—resembling the memory of microstructure—as essential novel ingredient. This approach can help understanding predicting and mitigating failure of structures or geophysical hazards, and will bring forward industrial process design and optimization, and help solving scientific challenges in fundamental research.

For many years, scientists and researchers have considered the jamming transition in granular materials to occur at a particular volume fraction, φ J (35).In contrast, over the last decade, numerous experiments and computer simulations have suggested the existence of a broad range of φ J , even for a given material.It was shown that the critical density for the jamming transition depends on the preparation protocol (11,17,28,(36)(37)(38)(39)(40)(41)(42)(43), and that this state-variable can be used to describe and scale macroscopic properties of the system (22).For example, rheological studies have shown that φ J decreases with increasing compression rate (7,42,44,45) (or with increasing growth rate of the particles), with the critical scaling by the distance from the jamming point (φ − φ J ) being universal and independent of φ J (19,28,36,46,47) Recently, the notion of an a-thermal isotropic jamming "point" was challenged due to its protocol dependence, suggesting the extension of the jamming point, to become a J-segment (32,45,48).Furthermore, it was shown experimentally, that for a tapped, unjammed frictional 2D systems, shear can jam the system (known as "shear-jamming"), with force chain networks percolating throughout the system, making the assemblies jammed, rigid and stable (6,33,34,49), all highlighting a memory that makes the structure dependent on history H.But to the best of our knowledge, quantitative characterization of the varying/moving/changing transition points, based on H, remains a major open challenge.
Here, we consider frictionless sphere assemblies in a periodic system, which can help to elegantly probe the behavior of disordered bulk granular matter, allowing to focus on the structure, without being disturbed by other non-linearities (6,50) (e.g.friction, cohesion, walls).For frictionless assemblies, it is often assumed that the influence of memory is of little importance, maybe even negligible.However, we demonstrate its relevance and quantitatively explore its structural origin in systems where the structure and its re-arrangements are the only possible mechanisms leading to the range of jamming points.
In this study, we probe the jamming transition concept by two pure deformation modes: isotropic compression and deviatoric pure shear (volume conserving), which allow us to combine the J-segment concept with a history dependent jamming density.Assuming that all other deformations can be superimposed by these two pure modes, we coalesce the two concepts of isotropic and shear induced jamming, and provide the unified model picture, involving a multiscale, fractal-type energy landscape (17,51,52); in general, deformation (or the preparation procedure) modify the landscape and its population; considering only changes of the population already allows to establish new configurations and to predict their evolution.The observations of different φ J of a single material require an alternative interpretation of the classical "jamming diagram" (4).
Our results will provide a unified picture, including some answers to the open questions from literature: (i) What happens to the shear-jamming regime in 3D and is friction important to observe it?-as posed by Bi et al. (6); (ii) What lies in between the jammed and flowing (unjammed) regime?-as posed by Coniglio et al. (48); (iii) Is there an absolute minimum jamming density?-as posed by Coniglio et al. (48); (iv) What protocols can generate jammed states?-asposed by Torquato et al. (41).Eventually, accepting the fact that the jamming density is changing with deformation history, significant improvement of continuum models is expected, e.g., for anisotropic models (2), GSH rate type models (53), or continuum models with a length scale (54).For this purpose we provide a simple (usable) analytical model that can be used for modifications or generalization of continuum models.Only allowing φ J (H) to be history dependent, as key modification, explains the multitude of reported observations and can be applied for real-world applications in e.g.electronic industry related novel materials (55).Results

M
The minimum (lower bound) of all M φ J,i is defined as the shear-jamming limit point, φ SJ = 0.6567.The solid lines through the data are universal fits to a stretched exponential (56)(57)(58) with only one single variable parameter φ max J , i.e., the upper limit jamming point for M → ∞, which depends on φ max i (SI section S2).(b) Schematic jamming phase diagram in volume fraction φ -strain ε d space, which includes the shear unjammed (SU), fragile (F) and shearjammed (SJ) states.Below φ SJ , the states are unjammed and application of ε d does not lead to jamming.Over-compression to very high volume fraction over many cycles M → ∞, see (a), lets the jamming point φ J (H) move towards its upper bound of available jamming points, φ max J , above which the systems are always jammed.The double arrow indicates the movement of φ J (H), increasing for over-compression and decreasing for shear; the solid black line indicates the shear-jamming transition between (SU+F) and (SJ) states, due to changes in M φ J,i .
Many different isotropic jamming points can be found in real systems and -as shown herealso for the simplest model material in 3D.We define a jamming "point" as the volume fraction, where the pressure on the unloading branch of an isotropic deformation cycle drops to zero (a cycle means loading and unloading 0, as described in the section Methods (2) (SI Fig. S1.(a)).From a relaxed, stress free initial state with volume fraction, φ t = 0.64 < φ J , we compress it isotropically to a maximum volume fraction, φ max i , and decompress back to φ t , as illustrated in Supporting Movie 1, and repeat this for M = 100 cycles.Fig. 1(a) shows the evolution of the isotropic jamming points M φ J,i , which increase with increasing M and with over-compression φ max i ; for subsequent cycles M of over-compressions, the jamming density M φ J,i grows slower and slower and is best captured by a Kohlrausch-Williams-Watts (KWW) stretched exponential relation : with the three universal "material"-constants φ SJ = 0.6567 (SI section S6), µ i = 1, and β i = 0.3, the lower limit of possible φ J 's, the relaxation (cycle) scale and the stretched exponent parameters, respectively.Only ∞ φ J,i , the equilibrium (steady-state or shakedown (59)) jamming point limit (extrapolated for M → ∞), depends on the over-compressions φ max i .Very little overcompression, φ max i φ SJ , does not change/increase φ J with cycles M , and thus provides a lower limit, φ SJ , of the jamming point range.Thus, the isotropic jamming point φ J is not a unique point, not even for frictionless particle systems, and is dependent on the previous deformation history of the system (48,60), e.g.over-compression or tapping/driving (data not shown).Both (isotropic) modes of deformation lead to more compact, better packed configurations (6,33,58).Considering different system sizes, and different preparation procedures, we confirmed that the jamming regime is the same (within fluctuations) for all the cases considered, see SI Fig. S3.All our data are consistent with an important conclusion: Smooth/slow isotropic deformation from random, dilute, unjammed states leads to jamming at a lower-limit density φ SJ .Unfortunately this limit is not well defined, i.e. it is highly protocol dependent, not well reproducible, and thus hard to determine experimentally and numerically as well.Reason is that any slow deformation (e.g.compression from below jamming) also leads to perturbations (like tapping leads to granular temperature): the stronger the system is perturbed, the better it will pack, so that usually φ J > φ SJ is established.Repeated perturbations, see SI Fig. S2.(a), lead to a slow stretched exponential approach to an upper-limit jamming density φ J → φ max J that itself increases slowly with perturbation amplitude, see SI Fig. S2.(b).The observation of different φ J of a single material, was referred to as J-segment (48,60), and only requires an alternative interpretation of the classical "jamming diagram" (4,6) and to give-up the misconception of a single, constant jamming "point".The state variable φ J varies due to deformation, but possibly has a lower limit that we denote for now as φ SJ .Jammed states below φ SJ might be possible, but require different protocols (61), or different materials, and are thus not addressed here.The concept of shear jammed states (6) below φ J , as illustrated in Fig. 1(b), is discussed next.

Shear jamming below φ J (H)
To study shear-jamming, we choose several unjammed states with volume fractions φ below their jamming points 1 φ J,i , which were established after the first compression-decompression cycle, for different history, i.e., various previously applied over-compression to φ max i (SI Fig. S4).Each configuration is first relaxed and then subjected to four isochoric (volume conserving) pure shear cycles (see Methods).We confirm shear jamming, e.g., by a transition in the coordination number C * , from below to above its isostatic limit, C * 0 = 6, for frictionless grains (2,12,25,29).This was consistently (independently) reconfirmed by using percolation analysis (6,24), allowing us to distinguish the three different regimes namely, unjammed, fragile and shear jammed states during (and after) shear, as shown in SI Fig. S6.(a).For this, we study the k−cluster, defined as the largest network of strong forces, f ≥ k f (62,63).When the initially unjammed isotropic system is sheared, it is growing and percolating first in the compressive direction, then in the neutral, non-mobile direction, and last in the extension direction.The best criterion for identifying growing clusters, the largest of which will percolate in all the three directions, is k = 2.2, different from k = 1 for 2D frictional systems (6).During shear deformation, the fraction of non-rattlers, f NR , increases from initially zero to large values, still well below unity, due to the always existing rattlers.For f NR > 0.82 ± 0.01, we observe that the growing force network is percolated in all three directions (SI Fig. S6), which is astonishingly similar to the value reported for the 2D systems (6).The evolution of the force network with f NR and its percolation in different directions, as viewed from two planes, is illustrated in Supporting Movie 2. From this perspective, when an unjammed material is sheared at constant volume, and it jams after application of sufficient shear strain, clearly showing that the jamming point has moved to a lower value.Shearing the system also perturbs it, just like over-compression; however, in addition, finite shear strains enforce shape-and structure-changes and thus allow the system to explore new configurations; typically, the elevated jamming density φ J of a previously compacted system will rapidly decrease and exponentially approach its lower-limit, the shear jamming "point" φ SJ , below which no shear-jamming exists.These quantification of history dependent jamming densities φ J (H), due to shear complementing the slow changes by cyclic isotropic (over)compression in Eq. ( 1), is discussed next.

Jamming phase diagram with history H
We propose a jamming phase diagram with shear strain, and present a new, quantitative history dependent model that explains jamming and shear-jamming, but also predicts that shearjamming vanishes under some conditions, namely when the system is not tapped, tempered or over-compressed before shear is applied.Using ε d and φ as parameters, Fig. 2(a) shows that for one initial the history dependent jamming state at 1 φ J,i , there exist sheared states within the range φ SJ ≤ φ ≤ φ J (H), which are isotropically unjammed.After small shear strain they become fragile, and for larger shear strain jam and remain jammed, i.e., eventually showing the critical state flow regime, where pressure, shear stress ratio and structural anisotropy have reached their saturation levels and forgotten their initial state (SI Fig. S5).The transition to fragile states is accompanied by partial percolation of the strong force network, while percolation in all directions indicates the shear-jamming transition.Above jamming, the large fraction of non-rattlers provides a persistent mechanical stability to the structure, even after shear is stopped.showing the different states: unjammed, isotropic jammed, shear unjammed, fragile and shear jammed, for one particular case of φ J (φ max i = 0.82, M = 1) =: 1 φ J,i = 0.6652.(b) Plot of minimum strain needed to jam states prepared from the first over-compression cycle with different φ max i , as given in the legend.Inset shows collapse of the states using a scaled definition that includes distance from both isotropic jamming point M φ J,i and shear-jamming point φ SJ , using Eq. ( 2).We only show data for the states for φ < 1 φ J,i that after the first isotropic compression decompression cycle jam by applying shear.
For φ approaching φ SJ , the required shear strain to jam ε SJ d increases, i.e., there exists a divergence "point" φ SJ , where 'infinite' shear strain might jam the system, but below which no shear jamming was observed.The closer the (constant) volume fraction φ is to the initial 1 φ J,i , the smaller is ε SJ d .States with φ ≥ 1 φ J,i are isotropically jammed already before shear is applied.
Based on the study of many systems, prepared via isotropic over-compression to a wide range of volume fractions φ max i ≥ φ SJ , and subsequent shear deformation, Fig. 2(b) shows the strains required to jam these states by applying pure shear.A striking observation is that independent of the isotropic jamming point 1 φ J,i , all curves approach a unique shear jamming point at φ SJ ∼ 0.6567 (SI section S6).When all the curves are scaled with their original isotropic jamming point M φ J,i as φ sc = (φ − φ SJ ) / M φ J,i − φ SJ they collapse on a unique master curve shown in the inset of Fig. 2(b), with power α = 1.37 ± 0.01 and shear strain scale ε 0 d = 0.102 ± 0.001 as the fit parameters.Hence, if the initial jamming point M φ J,i or φ J (H) is known based on the past history of the sample, the shear-jamming strain ε SJ d can be predicted.From the measured shear-jamming strain, Eq. ( 2), knowing the initial and the limit value of φ J , we now postulate its evolution under isochoric pure shear strain: Inserting, , φ J = φ and φ J = φ SJ , respectively.This means the jamming point evolution due to shear strain ε d is faster than exponential (since α > 1) decreasing to its lower limit φ SJ .This is qualitatively different from the stretched exponential (slow) relaxation dynamics that leads to the increase of φ J due to over-compression or tapping, see Fig. 3(a) for both cases.

Slow dynamics model
The last challenge is to unify the observations in a model that accounts for the changes in the jamming densities for both isotropic and shear deformation modes.Over-compressing a soft granular assembly is analogous to tapping (33, 58) more rigid ones, in so far that both methods lead to more compact packing structures, i.e., both represent isotropic perturbations.These changes are shown in Fig. 1(a), where the originally reported logarithmically slow dynamics for tapping ( 57) is very similar to our results that are also very slow, with a stretched exponential behavior; such slow relaxation dynamics can be explained by a simple Sinai-Diffusion model of random walkers in a random, hierarchical, fractal-type free energy landscape (3,56) in the a-thermal limit, where the landscape does not change -for the sake of simplicity.The granular packing is represented in this picture by an ensemble of random walkers in (arbitrary) configuration space with (potential) energy according to the height of their position on the landscape.(Their average energy corresponds to the jamming density and a decrease in energy corresponds to an increase in φ J (H), thus representing the "memory" and history dependence.)Perturbations, such as tapping with some amplitude (corresponding to "temperature") allow the ensemble to find denser configurations, i.e., deeper valleys in the landscape, representing larger (jamming) densities.Similarly, over-compression is squeezing the ensemble "down-hill", also leading to an increase of φ J , as presented in Fig. 3(b).Larger amplitudes will allow the ensemble to overcome larger barriers and thus find even deeper valleys.Repetitions have a smaller chance to do so, which explains the slow dynamics in the hierarchical multi-scale structure of the energy landscape.
In contrast to the isotropic perturbations, where the random walkers follow the "down-hill" trend, shear is anisotropic and thus pushing parts of the system "up-hill".For example, under planar simple shear, one (eigen) direction is tensile (up) whereas an other is compressive (down).If the ensemble is random, shear will only re-shuffle the population.But if the material was previously forced or relaxed towards the (local) land-scape minima, shear can only lead to a net up-hill drift of the ensemble, i.e., to decreasing φ J , referred to as dilatancy.For ongoing perturbation, if volume is conserved, both coordination number and pressure slowly decrease (SI Fig. S5) whereas for fixed confining pressure (data not shown) the volume would decrease (compactancy).This process is much faster than relaxation, since it is driven by shear strain amplitude.For large enough strain the system will be sufficiently re-shuffled, randomized, or "re-juvenated" such that it can be close to its quenched, random state φ SJ .

Prediction: minimal model
Finally, we test the proposed history dependent jamming point φ J (H) model, by predicting p and C * , when a granular assembly is subjected to cyclic isotropic compression to φ max i = 0.73 for M = 1 and for M = 300 cycles, with ∞ φ J,i = 0.667, as shown in Fig. 4(a-b) (prediction details discussed in SI section S7).It is observed that using the history dependence of φ J (H),  , leading to an increase in φ J (H) by slow stretched exponential relaxation.Dashed lines represent the much faster decrease in φ J (H) due to shear strain ε d , using Eq. ( 3).(b) The sketch represents only a very small, exemplary part of the hierarchical, fractal-type landscape.The red horizontal line represents the (quenched) average, while the dotted horizontal line indicates the momentary average φ J (H).The blue solid arrows show relaxation due to perturbations, while the dashed arrows indicate re-arrangements (rejuvenation) due to large shear strain.The green dots represent with their size the population after some relaxation, in contrast to a random, quenched population where all equal size valleys are equally populated (52).
the hysteretic behavior of the isotropic quantities, p and C * , is very well predicted, qualitatively similar to isotropic compression and decompression of real 2D frictional granular assemblies, as shown in Fig. 7 by Bandi et al. (43).
In Fig. 4(c), we show the evolution of the deviatoric quantities shear stress ratio τ/p and deviatoric fabric F d , when a system with φ = 0.6584, close to φ SJ , and initial jamming point φ J (0) = 0.6652, is subjected to three shear cycles (lowest panel).The shear stress ratio τ/p is initially undefined, but soon establishes a maximum (not shown) and decays to its saturation level at large strain.After strain reversal, τ/p drops suddenly and attains the same saturation value, for each half-cycle, only with alternating sign.The behavior of the anisotropic fabric F d is similar to that of τ/p.During the first loading cycle, the system is unjammed for some strain, and hence F d is zero in the model (observations in simulations can be non-zero, when the data correspond to only few contacts, mostly coming from rattlers).However, the growth/decay rate and the saturation values attained are different from those of τ/p, implying a different, independent stress-and structure-evolution with strain -which is at the basis of recently proposed anisotropic constitutive models for quasi-static granular flow under various deformation modes (2).The simple model with φ J (H), is able to predict quantitatively the behavior the τ/p and F d after the first loading path, and is qualitatively close to the cyclic shear behavior of real 2D frictional granular assemblies, as shown in Supplementary Fig. 7 by Bi et al. (6).
At the same time, also the isotropic quantities are very well predicted by the model, using the simple equations from SI section S7, where only the jamming point is varying with shear strain, while all material parameters are kept constant.Some arbitrariness involves the sudden changes of φ J at reversal, as discussed in SI section S7.Therefore, using a history dependent φ J (H) gives hope to understand the hysteretic observations from realistic granular assemblies, and also provides a simple explanation of shear-jamming.Modifications of continuum models like anisotropic models (2), or GSH type models (53), by including a variable φ J , can this way explain also transitions around jamming..73(small symbols) and decompression (big symbols) back to φ t , with ∞ φ J,i = 0.667, for M = 1 (red '•') and for M = 300 (blue ' ').(c) Deviatoric stress ratio τ/p and deviatoric fabric F d , fraction of non-rattlers f NR , coordination number C * , pressure p and history dependent jamming point φ J (H) over three pure shear strain cycles (bottom panel) for φ = 0.6584 and initial jamming point φ J (φ max i = 0.82, M = 1) =: 1 φ J,i = 0.6652.Solid lines through the data are the model prediction, involving the history dependent jamming point φ J (H), using Eq.(1) for isotropic deformation and Eq.(3) for shear deformation, and others, as described in SI section S7.Horizontal red lines in f NR and C * represent transition from unjammed to shear jammed states, whereas in φ J (H) indicates the shear jamming point φ SJ .

Interpretation and
In summary, the questions posed in the introduction can now be answered: (i) Shear-jamming occurs in 3D without any friction; (ii) the transition between the jammed and flowing (unjammed) regimes is controlled by a single, isotropic, history dependent state variable, the jamming density φ J (H), which (iii) has an absolute minimum jamming density; so that (iv) the protocol dependence of jamming is completely explained by the new state variable.
Given an extremely simple model picture, starting from an isotropically unjammed point, shear-jamming is not anymore a new effect, but is just due to the shift of the state variable jamming density to lower values during shear.(The evolution equations with their material parameters are determined quantitatively from a set of different simulations).Shear jamming occurs when the variable φ J (H) crosses the (fixed) φ of the system.The model implies now a minimum φ J ≥ φ SJ that represents the critical (steady) state in the limit of vanishing confining stress, i.e., the lower limit of all jamming points.This is nothing else than the lowest stable random density a sheared system locally can reach due to continuously ongoing shear, for vanishing confining stress.
This lower limit might be difficult to access in experiments and simulations, since every shear also perturbs the system leading at the same time to (slow) relaxation and thus a competing increase in φ J (H).However, it can be obtained from the relaxed, critical state values of pressure, extrapolated to zero, i.e., from the envelope in SI Fig. S8.The history dependent jamming point φ J (H) is difficult to access directly, but can consistently be extracted from experimentally measurable quantities, e.g.pressure p, coordination number C * or fraction of non-rattlers f NR .We explain the methodology to extract φ J (H) experimentally, and confirm by indirect measurement, as details in SI section S7 -S8 that the jamming density is indeed increasing during isotropic deformation and decreasing during shear.
Experiments should be performed to calibrate our model given real materials.Over-compression is possible for soft materials, but not expected to lead to considerable relaxation due to the small possible compressive strain for harder materials.However, tapping or small-amplitude shear can take the role of over-compression, also leading to perturbations and increasing φ J ; in contrast, large-amplitude shear is decreasing φ J and can be calibrated indirectly from different isotropic quantities.The accessible range of φ J − φ SJ is expected to much increase for more realistic systems, e.g., with friction, or non-spherical particle shapes.A measurement of the landscape, e.g. the valley width, depth and shapes (52) should be done to verify our model-picture, as this remains qualitative so far.Finally, the link of our model to glassy dynamics and frequency dependent response, encompassing also stretched exponential relaxation, see for example Ref. (65), is another open challenge open for future research.

Methods
Discrete particle simulations are used to model the deformation behavior of systems with N = 9261 soft frictionless spherical particles with average radius r = 1 [mm], density ρ = 2000 [kg/m 3 ], and a uniform polydispersity width w = r max /r min = 3, using the linear visco-elastic contact model in a 3D box with periodic boundaries (2).The particle stiffness is k = 10 8 [kg/s 2 ], contact viscosity is γ = 1 [kg/s], and the smallest contact duration is t c = 0.2279 [µs] (2).Firstly, the particles are generated with random velocities at φ = 0.3 and are isotropically compressed to φ t = 0.64 < φ SJ , and later relaxed.The system is further isotropically compressed to different φ max i and decompressed back to φ t , and is repeated over M cycles, which provides the M φ J,i (SI section S1).Later, several isotropic configurations φ, such that φ t < φ < 1 φ J,i from the decompression branch are chosen as the initial configurations.We relax them and apply pure (volume conserving) shear with the strain-rate tensor Ė = ± ǫd (−1, 1, 0), for four cycles (SI Figs.S4 and S5).The x and y walls move, while the z wall is stationary (Supporting Movie 3).Note that the strain rate of the (quasi-static) deformation is small, ǫd t c < 3.10 −6 , to avoid transient behavior.
Pressure P = tr(σ)/3, where σ is the stress tensor (2,6).The fabric tensor is defined as F = (1/V ) P∈V V P c∈P n c ⊗ n c , where V P is the particle volume for all particles P lying in V , and n c is the normal unit branch-vector pointing from center of particle P to contact c.Isotropic part of fabric is F v = tr(F).The corrected coordination number (2, 6) is C * = M 4 /N 4 , where, M 4 is total contacts of the N 4 particles having at least 4 contacts, and the non-rattler fraction is f NR = N 4 /N .For any tensor Q, its deviatoric part is, Q d = sgn (q yy − q xx ) 3q ij q ij /2, where q ij are the components of the deviator of Q, and the sign function accounts for the shear direction.Both pressure P and shear stress Γ are non-dimensionalized by 2 r /k to give dimensionless pressure p and shear stress τ.The shear-jamming point is obtained using the relaxed, critical state data sets of p, as presented in SI section S6.

Supplemental Materials: Memory of jamming and shear-jamming (in soft and granular matter)
In this report, we present in detail the methodology and results from elaborate data analysis used to understand the slow relaxation/compaction experiments and their structural origin, and the relation between jamming and shear-jamming in three-dimensions for frictionless, polydisperse spherical particles in a periodic box.The report is organized as follows: In section S1, we present a procedure to identify the jamming-points and their range.In section S2, we show the effect of cyclic over-compression to different target volume fractions and present a model that captures this phenomena.Section S3 is devoted to show a statistical analysis and sensitivity of the jamming point with the system size (number of particles) and different initial configurations.Section S4 is devoted to the effect of cyclic pure shear experiments and their effect on the microscopic and macroscopic quantities.In section S5, we discuss the results on growing, percolating clusters during cyclic shear with focus on reversal and their link with the non-rattler fraction.In section S6, we explore the phase space of macroscopic quantities during shear and the effects of relaxation, compaction and creep in critical (steady) state, after large strain amplitude, when the initial state is forgotten.There, we also present the methodology used to obtain the shear jamming point.In section S7, we show the analytical model, especially the history dependent jamming point, that is used for the prediction of macroscopic quantities.Finally in section S8, we present a procedure to extract the history dependent jamming point from experimentally measurable quantities.

S1 Identification of the jamming point
In the Main text of this paper, we have presented ranges of the isotropic jamming points, where pressure vanishes, that are possible due to different sample histories.Before we explain the procedure for identification of jamming points, we want to point out here that the stress tensor σ, that is used to define the dimensionless pressure p and shear stress τ, has static and dynamic contributions, with the latter being four orders of magnitude smaller than the former (in the cases studied) and hence the dynamic contributions are neglected.
When a sample is over-compressed isotropically, the loading and unloading paths are different in pressure p.This difference is most pronounced near the jamming point φ J , and for the first cycle.It brings up the first question of how to identify a jamming point, φ J .The unloading branch of a cyclic isotropic over-compression along volume fraction φ is well described by a linear relation in volumetric strain, with a quadratic correction: where p 0 , γ p presented in Table S1, and the jamming point φ J are the fit parameters.−ε v = log(φ/φ J ) is the true or logarithmic volumetric strain of the system, defined relative to the reference where p → 0, i.e. jamming volume fraction.C is the ratio of total non-rattler contacts M 4 and total number of particles N , i.e., C = M 4 /N = (M 4 /N 4 ) (N 4 /N ) = C * f NR , with corrected coordination number C * and fraction of non-rattlers f NR .Using Eq. (S1), we can extrapolate p to zero, to get φ J .We apply the same procedure for different over-compressions, φ max i , and many subsequent cycles M to obtain M φ J,i , for which the results are discussed below.The material parameter p 0 is finite, almost constant, whereas γ p is small, sensitive to history and contributes mainly for large −ε v , with values ranging around 0 ± 0.1; in particular, it is dependent on the over-compression φ max i (data not shown).Unless strictly mentioned, we shall be using the values of p 0 and γ p given in Table S1.= 0.82, the jamming point, as realized after unloading, is increasing.Also, with each cycle, from M = 1 to M = 100, the jamming point moves to larger values.Note that the difference between the loading and the unloading curves becomes smaller for subsequent over-compressions.Fig. S1.(b) shows the scaled pressure, i.e., p normalized by φC/φ J , which removes its non-linear behavior.p represents the average deformation (overlap) of the particles at a given volume fraction, proportional to the distance from the jamming point φ J .In the small strain region, for all over-compression amplitude and cycles, the datasets collapse on a line with slope p 0 ∼ 0.04.Only for very strong over-compression, −ε v > 0.1, a small deviation (from linear) of the simulation data is observed due to the tiny quadratic correction in Eq. (S1).The inset is the zoomed in regime near the jamming point.Lines are connecting the datasets only (b) Scaled pressure pφ J /φC plotted against volumetric strain −ε v = log(φ/φ J ) for the same simulations in (a).The φ J are extracted using a fit to Eq. (S1).Lines passing through the data represents scaled pressure, when Eq. ( S1) is restructured.

S2 Isotropic cyclic over-compression
In this section, we present a simple model that captures the increasing M φ J,i behavior with increasing over-compression, φ max i , as well as cycles, M .In Fig. S2.(a) (the same data as in Main text Fig. 1(a) are shown on a logarithmic scale, with fits using Main text Eq. ( 1)), weak over-compression does not lead to a significant increase in φ J,i , giving us information about the lower limits of the isotropic jamming points, which is the shear jamming point φ SJ = 0.6567, see section S6.With each over-compression cycle, M φ J,i increases, but for large M it increases less and less.This is analogous to compaction by tapping, where the tapped density increases logarithmically slow with the number of taps.Several models are available in literature to depict this density relaxation behavior e.g.hyperbolic tangent, stretched exponential, inverse logarithmic, reciprocal linear laws etc..We tried these four models and found that they fit our data qualitatively well; however, the data is best predicted by a Kohlrausch-Williams-Watts (KWW) stretched exponential relation (Main text Eq. ( 1)).It is interesting to note that with increasing φ max i , the limit value ∞ φ J,i increases and gives an upper bound to the isotropic

S3.(b):
Figure S3: (a) Effect of system size (number of particles N ) on the isotropic jamming point φ J .The preparation procedure is the same for all cases.All the samples were over-compressed from φ t = 0.64 to φ max i and decompressed back to φ t .(b) Effect of different initial system configurations on the isotropic jamming point φ J , which were prepared by applying time t of compression from φ = 0.3 to φ t .The compression time t is scaled with the smallest contact time between two particles t c .Later, all the samples were compressed to various φ max i and then decompressed back to φ t with the same rate for the first cycle.The red line indicates the shear jamming limit φ SJ .configurations, and hence independence on the initial configurations, provided that the main cycle-experiments are performed in the quasi-static limit.

S4 Cyclic Shear results
Here, we will discuss the effect of cyclic pure shear on some microscopic and macroscopic quantities, namely the corrected coordination number C * , the pressure p, and stress τ, the shear stress ratio τ/p, and the anisotropic fabric F d and their evolution for different volume fractions below isotropic jamming point.Firstly, from the observations in section S2, it is easy to say that there exist φ such that φ SJ ≤ φ ≤ φ max J , where it is possible to generate isotropic unjammed states and jammed states corresponding to higher and lower over-compressions respectively, is shown in Fig. S4.Secondly, shearing at constant volume can generate shear unjammed, fragile and finally shear jammed states for different shear strain ε d values, as seen in Fig. S6.
In Fig. S5, we plot the evolution of microscopic and macroscopic quantities during volume conserving cyclic shear experiments for φ = 0.6593.Since the state is below jamming, C * grows slowly from 0 and crosses the isostatic coordination number C * 0 = 6 as shown in Fig. S5.(a).As we apply further shear strain, C * attains an asymptotic value with small fluctuations.At reversal, the particles start to loose contacts due to changes in the main force network direction, hence C * decreases sharply.Afterwards, with increase in opposite shear strain, C * starts to build up again and reaches a similar, fluctuating limit value as for initial strain.For p, similar behavior like for C * vs ε d , is seen in Fig. S5.(b), where p grows from 0 (up to finite strain ε d ≤ 0.08) to higher values until it reaches its asymptotic value, with relatively larger fluctuations.After the initial cycle, both C * and p vs ε d collapse on top of each other.Note that the plot of C * and p vs ε d is symmetric, not about the original box configuration but about the mean of the large extremes ε max d and ε min d , i.e., 0.5 * (0.28 − 0.16) = 0.06, so in a way providing history-independence, and deviatoric limit state symmetry.In Figs.S5.(c) and S5.(d), we plot the variation of C * and p vs ε d for different volume fractions, respectively, for the second cycle.For higher volume fractions, we observe higher C * and p limits, because of a higher probability for making contacts with the neighboring particles.Figs.S5.(e) and S5.(f) show the behavior of shear stress ratio, τ/p and deviatoric fabric F d with strain ε d respectively.For higher shear strains, both τ/p ≈ 0.22 and F d ≈ 0.16 saturate, with no visible effect of volume fraction, and similar characteristics for high opposite shear strains.However, their growth rates with applied strain are different.It is interesting to note that at reversal, the systems stay highly anisotropic in fabric and stress for some strain ε d , whereas the coordination number C * and pressure p drop down much more rapidly at reversal.After reversal, it also takes some finite strain for C * and p to recover to C * > C * 0 = 6 and p > 0.

S5 Percolation analysis
Here, we will discuss the percolation analysis, that allows to distinguish the three regimes namely, unjammed, fragile and shear jammed states during (and after) shear, as shown in Fig. S6.(a).We study how the k−cluster, defined as the largest force network, connecting strong forces, f ≥ kf avg , with k = 2.2, percolates during shear.More quantitatively, for an exemplary volume fraction φ (φ max i = 0.82, M = 1) = 0.6584, very close to φ SJ , Fig. S6.(b) shows that f NR increases from initially zero to large values well below unity due to the always existing rattlers.The compressive direction percolating network ξ y /L y grows faster than the tension direction network ξ x /L x , while the network in the non-mobile direction, ξ z /L z , lies in between them.For f NR > 0.82 ± 0.01, we observe that the growing force network is percolated in all three directions (Fig. S6.(a)).The jamming by shear of the material corresponds (independently) to the crossing of C * from the isostatic limit of C * 0 = 6, as presented in Fig. S6.(b).Next, we present the evolution of the strong force networks in each direction during cyclic     shear, as shown in Fig. S5, for the same initial system.After the first loading, at reversal f NR drops below the 0.82 threshold, which indicates the breakage/disappearance of strong clusters, i.e. the system unjams.The new tension direction ξ y /L y drops first with the network in the non-mobile directions, ξ z /L z , lying again in between the two mobile direction.With further applied strains, f NR increases and again, the cluster associated with the compression direction grows faster than in the tension direction.For f NR above the threshold, the cluster percolates the full system, leading to shear jammed states again.At each reversal, the strong force network breaks/fails in all directions, and the system gets "soft" or even un-jams temporarily.However, the network is rapidly is re-established in the perpendicular direction, i.e., the system jams and the strong, anisotropic force network again sustains the load.It is important to note that for some systems with volume fractions away from φ SJ can resist shear strain reversal (data not shown).

S6 Relaxation effects
In this section, we will discuss the system stability by looking at the macroscopic quantities in the saturation state (after large shear strain), by relaxing them sufficiently long to have non-fluctuating values in the microscopic and macroscopic quantities.Every shear cycle after defining e.g. the y−direction as the initial active loading direction, has two saturation states, one during loading and, after reversal, the other during unloading.In Fig. S8, we show values attained by the isotropic quantities pressure p, isotropic fabric F v and the deviatoric quantities shear stress τ, shear stress ratio τ/p, and deviatoric fabric F d for various φ given the same initial jamming point φ J (φ max   .82,M = 1) =: 1 φ J,i = 0.6652.Black 'x' symbols represent the initial loading cycle, while the green '+' and blue ' * ' represent states attained for φ < φ J and φ > φ J , respectively for the subsequent shear.Cyan '•' and the brown ' ' are states chosen after large strain during loading and unloading shear respectively, and are relaxed.The red and purple lines indicate the shear jamming point φ SJ = 0.6567 and the jamming point 1 φ J,i respectively.
as well as at the two relaxed saturation states (averaged over four cycles), leading to following observations: (i) With increasing volume fraction, p, F v and τ increase, while a weak decreasing trend in stress ratio τ/p and deviatoric fabric F d is observed.(ii) There is almost no difference in the relaxed states in isotropic quantities, p and F v for the two directions, whereas it is symmetric about zero for deviatoric quantities, τ, τ/p, and F d .The decrease in pressure during relaxation is associated with dissipation of kinetic energy and partial opening of the contacts to "dissipate" the related part of the contact potential energy.However, F v remains at its peak value during relaxation, showing that the contact structure is almost unchanged and the network remains stable during relaxation.(iii) For small volume fractions, close to φ SJ , the system becomes strongly anisotropic in stress ratio τ/p, and fabric F d rather quickly, during (slow) shear (envelope for low volume fractions in Figs.S8.(d) and S8.(e)), before it reaches the steady state.(iv) It is easy to obtain the shear jamming point φ SJ from the relaxed critical (steady) state pressure p, and shear stress τ, by extrapolation to zero, as the envelope of relaxed data in Figs.

S8.(a) and S8.(c).
We use the same methodology using Eq.(S1), to extract the shear jamming point φ SJ .When the relaxed p is normalized with the contact density φC, we obtain φ SJ = 0.6567 ± 0.0005 by linear extrapolation.A similar value of φ SJ is obtained from the extrapolation of the relaxed τ data set, and is consistent with other methods using the coordination number C * , or the energy (S1).

S7 Predictive power
In this section, we present the simplest model equations, as used for the predictions presented in Main text Fig. 4, involving a history dependent φ J (H), as given in Main text Eq. ( 1) for isotropic deformations and Main text Eq. ( 3) for shear deformations.The only difference to Ref. ), where these relations are taken from, based on purely isotropic unloading, is a variable φ J = φ J (H).During (cyclic) isotropic deformation, the evolution equation for the corrected coordination number C * is: with C 0 = 6 for the frictionless case and parameters C 1 and θ are presented in Table S1.The fraction of non-rattlers f NR is given as: with parameters ϕ c and ϕ v presented in Table S1.We modify Eq. (S1) for the evolution of p together with the history dependent φ J = φ J (H) so that, with parameters p 0 and γ p presented in Table S1, and the true or logarithmic volume change of the system is −ε v = log(φ/φ J (H)), relative to the momentary jamming density.The noncorrected coordination number is C = C * f NR , as can be computed using Eqs.(S4) and (S5).The isotropic fabric F v is given by the relation F v = g 3 φC, as taken from Imole et al. (S2), with g 3 ∼ = 1.22 for the polydispersity used in the present work.Also the parameters C 1 , θ for C * , ϕ c , ϕ v for f NR ; and p 0 , γ p for pressure p are similar to Ref. (S2), with the second order correction parameter γ p most sensitive to the details of previous deformations; however, not being very relevant since it is always a small correction due to the product γ p (−ε v ).
The above relations are used to predict the behavior of the isotropic quantities: dimensionless pressure p and coordination number C * , as shown in Main text Fig. 4(a-b) during isotropic compression, as well as for the fraction of non-rattlers in Main text Fig. 4(c) for cyclic shear, with corresponding parameters presented in Table S1.Note that during isotropic deformation, φ J (H) was changed only during the compression branch, using Main text Eq. ( 1) for fixed M = 1 using φ max i as variable, but is kept constant during unloading/expansion.During cyclic (pure) shear deformation, a simplified equation for the shear stress ratio Evolution parameters where Q a , Q c and Ψ are the fitting constants with values presented in Table S2.
For predictions during cyclic shear deformation, φ J (H) was changed with applied shear strain ε d using Main text Eq. (3).Furthermore, the jamming point is set to a larger value just after strain-reversal which is discussed next.Behavior of jamming point at strain reversal As mentioned in section S5, there are some states below φ J , where application of shear strain jams the systems.The denses of those can resist shear reversal, but below a certain φ cr < φ J , shear reversal unjams the system again.With this information, we postulate the following: (i) After the first large strain shear, the system should forget, where it was isotropically compressed to before i.e., M φ J,i if forgotten and φ J = φ SJ is realized.(ii) There exists a critical volume fraction φ cr , above which systems can just resist shear reversal and remain always jammed in both forward and reverse shear.(iii) Below this φ cr , reversal unjams the system.Therefore, more strain is needed to jam the system (when compared to the initial loading), first to forget its state before reversal, and then to re-jam it in opposite (perpendicular) direction.Hence, the strain necessary to jam in reversal direction should be higher than for the first shear cycle.(iv) As we approach φ SJ , the reverse strain needed to jam the system increases.
We use these ideas and measure the reversal shear strain ε SJ,R d , needed to re-jam the states below φ cr , as shown in Fig. S9.When they are scaled with φ cr as φ sc = (φ − φ SJ ) / (φ cr − φ SJ ) , they collapse on a unique master curve, very similar to Main text Eq. ( 2 2), that includes distance from both φ cr and shear-jamming point φ SJ , using Eq.(S10).
The above relations are used to predict the isotropic and the deviatoric quantities, during cyclic shear deformation, as shown in Main text Fig. 4(c), with the additional rule that all the quantities attain value zero for φ ≤ φ J (H).Moreover, for any state with φ ≤ φ cr , shear strain reversal moves the jamming point to φ cr , and the movement of jamming point follows Eq. (S10).
Any other deformation mode, can be written of as a unique superposition of pure isotropic and pure shear deformation modes, and hence the combination of the above can be easily used to describe any general deformation, e.g.uniaxial cyclic compression (data not presented).

S8 How to measure φ J from experiments
Last, but not the least, we conclude by plotting the history dependent jamming point φ J (H) from measurable quantities, indirectly obtained via Eqs.(3), (S5), (S6), and directly from Main text Eq. (3).There are two reasons to do so: (i) the jamming point φ J (H) is only accessible in the unloading limit p → 0, which requires an experiment or a simulation to "measure" it (however, during this measurement, it might change again); (ii) deducing the jamming point from other quantities that are defined for an instantaneous snapshot/configuration for p > 0 allows to indirectly obtain it -if, as shown next, these indirect "measurements" are compatible/consistent: Showing the equivalence of all the different φ J (H), proofs the consistency and completeness of the model and, even more important, provides a way to obtain φ J (H) indirectly from experimentally accessible quantities.same master-curve for first over-compression as seen in Figs.S10, independent of the maximum -all following deformation is dependent on the previous maximum density..For shear deformation Fig. S11 shows the evolution of φ J (H), measured from the two experimentally accessible quantities: coordination number C * and pressure p, using Eqs.(S4) and (S6) respectively during volume conserving shear with φ = 0.66, and the initial jamming point φ J (φ max i = 0.82, M = 1) =: 1 φ J,i = 0.6652 > φ and shows good agreement with the theoretical predictions using Main text Eq. (3) after shear jamming.Thus the indirect measurements of φ J (H) can be applied if φ J (H) < φ; the result deduced from pressure fits the best, i.e., it interpolates the two others and is smoother.φ J (H)   Figure S11: Evolution of the history dependent jamming point φ J (H) during pure shear, calculated back from the measured quantities: coordination number C * , fraction of non-rattlers f NR and pressure p, using Eqs.(S4), (S5), (S6) respectively, as indicated in the legend.The volume fraction is constant, φ = 0.66, and the initial jamming point φ J (φ max i = 0.82, M = 1) =: 1 φ J,i = 0.6652 is greater than φ (represented by horizontal cyan line).The solid black line represents Main text Eq. ( 3), and the dashed vertical line indicates the shear strain needed to jam the system, ε SJ d , from which on -for larger shear strain -the system is jammed.

Figure 1 :
Figure 1: Jamming phase diagram taking into account the changing isotropic jamming points.(a) Evolution of isotropic jamming points M φ J,i after performing M isotropic compressiondecompression cycles up to different maximum volume fractions φ max i , as given in the inset.With increasing φ max i , the range of the established jamming points M φ J,i = φ J (M, φ max

Figure 2 :
Figure 2: Phase diagram and scaling with φ SJ to replace the M φ J,i 's.(a) Phase diagram showing the different states: unjammed, isotropic jammed, shear unjammed, fragile and shear jammed, for one particular case of φ J (φ max

Figure 3 :
Figure 3: Relaxation dynamics and energy landscape due to memory effects.(a) Evolution of the jamming points φ J (H) due to history H. Solid lines represent many isotropic compression decompression cycles for three different φ max i

Figure 4 :
Figure 4: Model prediction of cyclic loading: (a) Dimensionless pressure p and (b) coordination number C * plotted against volume fraction φ for an isotropic compression starting from φ t = 0.64 to φ max i = 0.73 (small symbols) and decompression (big symbols) back to φ t , with ∞ φ J,i = 0.667, for M = 1 (red '•') and for M = 300 (blue ' ').(c) Deviatoric stress ratio τ/p and deviatoric fabric F d , fraction of non-rattlers f NR , coordination number C * , pressure p and history dependent jamming point φ J (H) over three pure shear strain cycles (bottom panel) for φ = 0.6584 and initial jamming point φ J (φ max

Fig. S1 .
(a)  shows the behavior of p with φ during one full over-compression cycle to display the dependence of the jamming point on the maximum over-compression volume fraction and the number of cycles.With increasing over-compression amplitude, e.g.comparing φ max i = 0.68 S1 and φ max i

Figure S1 :
Figure S1: (a) Dimensionless pressure p plotted against volume fraction φ and for an isotropic compression starting from φ t = 0.64 to φ max i

Figure S5 :
Figure S5: (a-b) Coordination number C * and pressure p for the four shear cycles at constant volume fraction φ = 0.6593 and jamming point φ J (φ max i = 0.82, M = 1) =: 1 φ J,i = 0.6652.(c-f) C * , p, shear stress ratio τ/p, and anisotropic fabric F d for the second cycle only, for three different volume fractions φ as given in the insets.

Figure S6 :
Figure S6: (a) Plot of C * and cluster sizes ξ/L in the three directions for tension in x− and compression in y− directions against the non-rattler fraction f NR , along the loading path for an isotropic unjammed initial state with volume fraction φ = 0.6584 and φ J (φ max i = 0.82, M = 1) =: 1 φ J,i = 0.6652.The upward arrow indicates the direction of loading shear strain.(b) Snapshots of unjammed, fragile and shear jammed states, when the force networks are percolated in none, one or two, and all the three directions, respectively.Only the largest force network, connecting strong forces, f ≥ k f , with k ≥ 2.2 are shown for the three states for clarity, and hence the white spaces in the background.

Figure S8 :
Figure S8: Scatter plots of isotropic quantities (a) pressure p, (b) isotropic fabric F v and deviatoric quantities (c) shear stress τ, (d) shear stress ratio τ/p, and (e) deviatoric fabric F d for various φ and jamming point φ J (φ max i = 0.82, M = 1) =: 1 φ J,i = 0.6652.Black 'x' symbols represent the initial loading cycle, while the green '+' and blue ' * ' represent states attained for φ < φ J and φ > φ J , respectively for the subsequent shear.Cyan '•' and the brown ' ' are states chosen after large strain during loading and unloading shear respectively, and are relaxed.The red and purple lines indicate the shear jamming point φ SJ = 0.6567 and the jamming point 1 φ J,i respectively.
with F d 0 and F d max the initial and maximum (saturation) values of the deviatoric fabric, respectively, and β F its growth rate.The four parameters (τ/p) max , β s for τ/p and F d max , β F for F d are dependent on the volume fraction φ and are well described by the general relation from Imole et al. (S2) as:

=
Figure S9: Phase diagram showing the minimum reversal shear strain ε SJ,R d needed to jam the states below φ cr , for states prepared from the first over-compression cycle with different φ max i , as given in the legend.Inset shows collapse of the states using a similar scaled definition as Main text Eq. (2), that includes distance from both φ cr and shear-jamming point φ SJ , using Eq.(S10).

Table S2 :
Parameters for Eqs.(S7) and (S8) using Eq.(S9), with slightly different values than from Ref.(S2), that are extracted using the similar procedure as in (S2), for states with volume fraction close to the jamming volume fraction.