Un-jamming due to energetic instability: statics to dynamics

Jamming/un-jamming, the transition between solid- and fluid-like behavior in granular matter, is an ubiquitous phenomenon in need of a sound understanding. As argued here, in addition to the usual un-jamming by vanishing pressure due to a decrease of density, there is also yield (plastic rearrangements and un-jamming that occur) if, e.g., for given pressure, the shear stress becomes too large. Similar to the van der Waals transition between vapor and water, or the critical current in superconductors, we believe that one mechanism causing yield is by the loss of the energy’s convexity (causing irreversible re-arrangements of the micro-structure, either locally or globally). We focus on this mechanism in the context of granular solid hydrodynamics (GSH), generalized for very soft materials, i.e., large elastic deformations, employing it in an over-simplified (bottom-up) fashion by setting as many parameters as possible to constant. Also, we complemented/completed GSH by using various insights/observations from particle simulations and calibrating some of the theoretical parameters—both continuum and particle points of view are reviewed in the context of the research developments during the last few years. Any other energy-based elastic-plastic theory that is properly calibrated (top-down), by experimental or numerical data, would describe granular solids. But only if it would cover granular gas, fluid, and solid states simultaneously (as GSH does) could it follow the system transitions and evolution through all states into un-jammed, possibly dynamic/collisional states—and back to elastically stable ones. We show how the un-jamming dynamics starts off, unfolds, develops, and ends. We follow the system through various deformation modes: transitions, yielding, un-jamming and jamming, both analytically and numerically and bring together the material point continuum model with particle simulations, quantitatively.


Open challenges
Some open questions are: How can we understand phenomena that originate from the particle-or meso-scale, which is intermediate between atoms and the macroscopic, hydrodynamic scale? And how can we formulate a theoretical framework that takes the place of the Navier-Stokes equations?
A universal theory must involve all states granular matter can take, i.e., granular gases, fluids, and solids, as well as the transitions between those states. What are the state-variables needed for such a theory? And what are the parameters (that we call transport coefficients) and how do they depend on the state-variables?
Main goal of this paper is to propose a minimalist candidate for such an universal theory, able to capture granular solid, fluid, and gas, as well as various modes of transitions between these states. The model, remarkably, involves only four state-variables, density, momentum density (vector), elastic strain (tensor), and granular temperature. It is a boiled down, simplified case of the more complete theory GSH [47,48,[109][110][111][112], complemented by insights based on DEM, see Ref. [49] and references therein, modified such that it works below, above and during transitions. For the sake of transparency and treatability, whenever possible, we reduce most transport coefficients and parameters to constants-without loss of generality.
Each transport coefficient is related to the propagation or evolution of one (or more) of the state-variables that encompass the present state of the system. For simple fluids [1,113], it is possible to bridge between the (macroscopic) hydrodynamic and the (microscopic) atomistic scales; as an example, the diffusion coefficient quantifies mass-transport mediated by microscopic fluctuations.
In the case of low density gases, the macroscopic equations and the transport coefficients can be obtained using the Boltzmann kinetic equation as a starting point. For moderate densities, the Enskog equations provide a good, quite accurate description of dense gases (or fluids) of hard atoms [1], or of particles including the effects of dissipation, which results in what is nowadays referred to as standard kinetic theory (SKT) [7,114]. Beyond SKT one can only reach out (empirically) towards realistic systems [8,115,116], and beyond, see, e.g., [4,12]. The limit of granular fluids is where transport coefficients, like the viscosity, deviate from SKT. This involves a divergence [12,117,118] when the granular fluid becomes denser [8,116,118], approaching the jamming density from well below, i.e., the state that we could call a granular solid, as described by classical solid mechanics [120]. Recent research also considers soft particles [74,118,121] for which jamming changes from a sharp to a rather smooth transition. One objective of this paper is to bring together fundamental theoretical concepts of continuum mechanics [21,32,35,48,107,122,123] with observations made from particle simulations for simple granular systems in the gas, fluid, and solid states, including also the transitions between those states [8,49,116,118,[124][125][126]].

About states of granular matter
When exposed to external stresses, grains are elastically deformed at their contacts. In static situations, there is only elastic energy; in flowing states, some of the elastic energy is transferred to kinetic energy and back 1 , as sketched in Fig. 1 for the example of slow isotropic jamming, and-after an overcompression cycle-eventually un-jamming. Note that the jamming transition at small (yet finite) compression rate appears smooth/continuous, whereas the un-jamming transition is rather sharp/discontinuous [127].
Grains yield differently for vanishing or finite T g . In motion, for T g ≠ 0 , yield is a continuous phenomenon, i.e., state-variables vary continuously. If the grains are at rest initially, T g = 0 , yield is discontinuous-as evidenced by a layer of grains on a tilted plane. Discontinuity is mainly in the equilibrium value of the elastic stress. It is finite in the convex region and zero in the concave one, as it always relaxes away there. Any discontinuity of a phase transition is always in the equilibrium values of some quantities.
Particle systems, compressed and decompressed with different rates, are discussed in Ref. [127]. Here, a rather slow rate is chosen and the system response is plotted in Fig. 1. Even such little dynamics allows for a rich phenomenology below and above jamming, where the yellow area corresponds to a dense fluid with solid features-just below jamming, while the cyan area corresponds to a solid with fluid features-jammed, but strongly unstable (see the steps and wiggles in pressure and coordination number) [88].
The capability of granular solids to remain quiescent, in mechanical equilibrium, under a given finite stress is precarious. For small perturbations they will return to their original state. If pressure or shear stress become too large, the grains will, suddenly, start moving-either locally or globally [17,127]-with a decaying elastic stress. This qualitative change in behavior is an unambiguous phase transition. We shall refer to the region capable of maintaining the global, overall equilibrium of static grains as elastic, and its boundary (in the space spanned by the state-variables) as the yield surface. For local loss of elasticity we rather use the term plastic, irreversible events, see Refs. [19,61]. 1 As definition of states, flowing states range from dilute granular gases via inertial, collisional granular fluids, to quasi-static flows and plastically (irreversibly) deforming granular solids, excluding only perfectly static, elastic granular solids (e.g., probed by elastic waves). The most interesting regime is around quasi-static flows where both solid and fluid features are important [88], with considerable permanent and fluctuating energy densities, w e and w T , respectively, summing up to the total w = w e + w T . Note that w T not only contains all kinetic energy, E kin , but also the fluctuating part of the potential energy, E f pot , i.e., w e = (E pot − E f pot )∕V , as discussed next. The density (equivalent to the volume fraction, ∼ ) alone is not sufficient to characterize the state of a particle system. Even though "magic" densities like random-close or -loose, RCP and RLP , respectively, are often used, but being highly material dependent they are not unique state descriptors. In addition to density, the ratio of kinetic to potential, elastic energy in the system, K = E kin ∕E pot , is one more possibility to characterize its state: gas ( ≪ 1 , K ≫ 1 ), fluid ( < RLP , K > 1 ) dense collisional flow ( ∼ RLP , K ∼ 1 ), quasi-static flow ( ∼ RCP > RLP , K ≪ 1 ), granular solid ( K ≈ 0 ), static ( K = 0 ), and the extreme, athermal case ( K ≡ 0 , maintained at all times), as can be realized by energy minimization, e.g., see Ref. [126] and references therein. The contribution of potential energy to the total energy is 1∕(1 + K) , but using the fluctuating fraction of total energy defines the states as: gas ( w T ∕w = 2∕(1 + K) ≪ 1 , due to E pot = E f pot ), collisional ( w T ∕w ∼ 1 ), intermediate quasi-static and solid-like ( w T ∕w = 2K∕(1 + K) ≈ 2K ≪ 1 , due to equipartition E kin = E f pot ), static, solid ( w T ∕w = 0 ), and crystalline ordered, possibly at ≫ RCP . Main message to be made in this paper is that, besides density, only two additional scalar state-variables are sufficient to encompass all possible (isotropic) states and transitions of a system, namely the isotropic elastic strain-related to the jamming density, J , itself [49]-and a granular temperature T g ∝ v ∝ √ w T , involving the velocity fluctuations, v [18]. For anisotropic (sheared) states at least one more state-variable is needed. Granular systems will also un-jam for vanishing pressure and a continuous reduction of density, though we reserve the term yield for the (sudden) loss of elastic stability: Grains un-jam in either case, they yield only when the elastic stress, in particular the pressure, is finite.
Starting from the elastic region, decompression (tension) reduces the density and the elastic deformations of the grains-until the latter vanish and the system un-jams. Decompressing further just reduces the density accordingly. The system is now un-jammed in the sense that one can change the density without any restoring force, i.e., the elastic energy remains zero. In reverse, compression only increases the density, as long as it is smaller than the jamming density. At jamming both the elastic deformations and the associated energy start to increase with density. In contrast, there is a discontinuity leaving the elastic regime at finite values of elastic stress. It is a sudden transition from quiescent, enduringly deformed grains, to moving ones oscillatorily deformed due to "jiggling" particle motions. This transition needs to be explained, to have a model for.
And it is clear that the transition must be encoded in the elastic energy-the only quantity characterizing the quiescent state-not in the dynamic/fluctuating contributions to energy.
In the elastic region, grains appear solid when at rest, but they will flow if subject to an imposed finite shear rate, and appear liquid. Such a continuous change in appearance is well accounted for by any competent dynamic theory or rheology, it is not a transition 2 . Moreover, flowing grains in the elastic region do feature a macroscopic elastic shear stress, with an associated elastic energy (even though granular contacts switch continually), something no Newtonian liquid is capable of. Also, the shear stress remains finite when the grains stop flowing, which is not the case in Newtonian fluids.
So there are two different flowing states, either with finite elastic stress/strain, or with vanishing ones, which includes granular gases, as accounted for by the kinetic theory, see Refs. [8,116] and references therein. There is also a transition between them, as possibly related to (dry) liquefaction [37], but not to be confused with liquefaction due to a fluid between the particles, which is completely disregarded in this study of dry granular matter (even though the fluid stress can be considerable in wet system). We take both transitions, either leaving the quiescent state, or the flowing one, as the same transition, with the same underlying physics. In fact, encoding the transition in the elastic energy certainly affects the flowing state as well. The mechanism for yield is here related to elastic energy (irrespective whether the pressure or the shear stress is too large, or the density too small), as traditionally encompassed by concepts like plastic potentials, yield functions, or flow rules [21,30,32,34,106,122], see Fig. 2 and textbooks like Ref. [122].

Relation to other systems in physics
We do not think that the transition is due to spontaneously broken translational symmetry-the usual mechanism giving rise to static shear stresses, as in any fluid-solid transition. T g = 0 , as a function of the pressure P, shear stress s , and void ratio e, as rendered by an energy expression in [112]. Panel c is the 3D combination of a and b; with b depicting how the straight Coulomb yield line bends over, depending on the void ratio e-a behavior usually accounted for by cap models in elasto-plastic theories; while a depicts the maximal void ratio e = 1∕ − 1 (equivalent to the inverse density) plotted against pressure P, or the so-called virgin consolidation line (VCL). In panel (a), the dotted line is an empirical relation, e = e 1 − e 2 log(P∕P 0 ) , with P 0 = 0.5 MPa, e 1 = 0.679 and e 2 = 0.097 , approximating the VCL, but not valid for P → 0 . The thick solid line cuts the e-axis at e 0 , with the intersection being the lowest possible, random loosest packing value, see Ref. [112] for details, where also the thin solid line is discussed. Thus e 0 also defines the lowest possible jamming volume fraction, , see Ref. [49], with static, elastic states possible only below the VCL, as will be shown in Sect. 5 and 6 The quick argument is: Consisting of solid, grains already break translational symmetry. More importantly, the loss of equilibrium and granular statics is caused by the shear stress or pressure being too strong. This is an indication of an over-tightening phenomenon, of which the (pair-breaking) critical current is a prime example. If a superconductor conducts electricity without dissipation, it is in a current-carrying equilibrium state. If, however, the imposed current exceeds a maximal value, the system leaves equilibrium and enters a dissipative, resistive state. The superfluid velocity, v sf ∼ ∇ q , given by the gradient of a quantum mechanical phase, q , is the analogue of the strain. The dissipationless current, j sf = w∕ v sf , given by the derivative of the energy with respect to v sf , is the analogue of the elastic stress. The over-tightening transition in superconductivity is well accounted for by an inflection point, at which the energy turns from stably convex to concave, see the classic paper by Bardeen [130]. The close analogy between the two systems is a good reason to employ the same approach here, to postulate that the surface of the cone in Fig 2 can be related to an inflection surface of the elastic energy.

About elastic granular matter
The granular solid state is contingent on granular matter capable of being elastic, for which there is ample evidence, see e.g. Refs. [6,52,84,95,125,[131][132][133][134] and references therein. In addition to the material stiffness, many other material properties (including cohesion, friction, surfaceroughness, particle-shape) determine the elastic response of granular matter. For soft and stiff materials the deformations are, respectively, considerable and slight, but never zero. Because of their Hertz-like non-linear contacts, grains are infinitely soft in the limit of vanishing contact area (deformation). Therefore, at any given finite force, deformations are always sufficiently large to display the full spectrum of elastic behavior, including a considerable static shear stress (enabling a tilted surface), and elastic waves. Even the simplest model material, consisting of perfectly smooth spheres of isotropic, linearly elastic material, displays non-linearity due to their Hertz-type contacts, on-top of the contact network (fabric) and its re-structuring. Only in computer simulations is it possible to remove the first and focus on the second, see e.g. Ref. [49].

Yield: About the limits of elasticity
To envision the yield surface, we consider the space spanned by three parameters: pressure P, shear stress s , and void ratio e = (1 − )∕ (where = p , with material density p and volume fraction ), ignoring the granular temperature (i.e., fluctuations of kinetic and potential energy), as discussed in Ref. [144] and so many papers following. Based on the observation of the Coulomb yield and the virgin consolidation line, we assume that the yield surface is as rendered in Fig. 2. Elastic, jammed states, maintained by deformed grains, are stable and static only inside it 3 .
The Coulomb yield line, see Fig. 2b, can be reached by increasing the shear stress at given confining pressure. When the shear stress exceeds a certain level, the system yields, un-jams and becomes dynamic. No static, stable elastic state exists above the Coulomb yield line, as evidenced by a sand pile's steepest slope.
It is imperative to realize that (what we call) the Coulomb yield line is conceptually different from the peak shear stress achieved during the approach to the critical state at much larger strains. Coulomb yield is the collapse of static states-such as when one slowly tilts a plate carrying grains until they start to flow (max. angle of stability). Its behavior is necessarily encoded in the system's energy, because this phenomenon does not at all involve the system's dynamics. The critical state, including the peak shear stress-though referred to as "quasi-static"-is a fully dynamic and irreversible effect. It is accounted for by the stationary solution at given strain rates in GSH. The angle of repose (always smaller than the max. angle of stability) is in GSH given by the critical friction angle [47,48].
In the absence of shear stresses, the maximally sustainable pressure depends on the void ratio, e, as rendered in Fig. 2a. Starting from a given e, slowly increasing P, the grain-structure will collapse and yield at this pressure, to a smaller value of e, such that the final state is stable, static, and below the curve of Fig. 2a. This is because when applying a slowly increasing pressure, the point of collapse is (ever so) slightly above the curve; and the end point below it is typically also close. This evolution resembles a stair-case, with the granular medium increasing its density by hugging this curve, which frequently referred to as the virgin/primary consolidation line, or simply the pressure yield line. The line cuts the e-axis at the random loosest void ratio, e 0 , above which no elastic stable states exists.
Because of the pressure yield line, the Coulomb yield curve cannot persist for arbitrarily large P at given e. Rather, it bends over to form a "cap", as rendered in Fig. 2b, since an additional shear stress close to the pressure yield line will also cause the packing to collapse. (The shape of the cap depends on the interplay of isotropic and deviatoric deformations as well as the probability for irreversible, possibly large-scale re-structuring events of the micro-structure, i.e., the contact network, including also the sliding of contacts, but also breakage of particles, which is, however, excluded from this study. Whether the picture sketched here is sound without breakage remains an open question [21].) Merging Figs. 2a and 2b yields the elastic region below the yield surface, as given in Fig. 2c. Although the e-axis, for P, s = 0 , is also referred to as the loci of (isotropic) unjamming, the elastic stress goes continuously to zero here, because the grains are successively less deformed. There is, as already discussed above in Sect. 1.3, no discontinuous phase transition or yield here, except for the coordination number, see Fig. 1. Point is, concerning only this one point of the plot, if the elastic stresses vanish, both in the convex and concave region, nothing much resembling a discontinuous transition happens there. Isotropic jamming and un-jamming, as well as the discontinuity in the coordination number on the isotropic e-axis is discussed in detail at various spots in this paper, see Sects. 1.3, 3.2.3, and 4.1.1.
Next, all different symbols and nomenclatures are summarised.

Notation and symbols
This paper is a cooperation of co-authors, whose notational baggage from past publications clash with one another. In the dire need to compromise, we ask the readers to sufferwith us-using varying symbols and notations. Our statevariables are: density, , momentum density, v i , granular temperature, T g , and the elastic (true) strain, as summarized here.
1. The bulk density, , is related to the volume fraction, = ∕ p (with p the particles' material density), the porosity 1 − , and the void ratio e = (1 − )∕ . (Later, we shall choose units such that p = 1 , so that volume fraction and bulk density are identical 4 .) 2. The conserved momentum density defines the velocity v i = ( v i )∕ . The symmetric part of the velocity gradient is where differences between finite and linear strain theory are detailed, e.g., in Refs. [145][146][147]. Eigenvalues of the total strain rate ̇i j are positive for compression and negative for tension. The symbol v ij is usual in condensed matter physics [113,120,145,148]; it is also the one employed in most previous GSH-publications. The notation D ij is common in theoretical mechanics [32,107,146], while ̇i j , or ̇ , are used, e.g., in soil mechanics and related literature [106,122]. 3. Subscripts, such as i,j,k,l, refer to components of tensors in the usual index notation, with double-indices implying summation, the comma indicating a partial derivative, as in v (i,j) ; the superscript * denotes the respective traceless (deviatoric) tensor. Using the summation convention, the volumetric strain rate is abbreviated as: ̇v =̇l l = −v ll = −D ll = −tr , where the last term is in symbolic tensor notation. The deviatoric strain rate is thus is the second deviatoric invariant, insensitive to the sign convention. 4. The elastic (true) strain, e ij ≡ −u ij , as properly defined in theoretical mechanics, e.g., see Refs. [145,147], even for large deformations (implying the logarithmic definition due to its additivity, reversibility, and splitability into isotropic and deviatoric contributions), is the tensorial state-variable on which the elastic (potential) energy depends. It is always well-defined and unique, in contrast to the total or plastic strains, which are not, and thus will not be used as state-variables for (constitutive) modeling. The respective strain rates, however, are well-defined and thus are used. The strain rate was already given (see item 2.), ̇i j = −v ij , so that the plastic strain rate is defined as: ̇p ij =̇i j − d dt e ij (see also item 7.). 5. The isotropic elastic strain is positive for compression 5 . It may be seen as the true (finite) strain relative to a stress-free reference configuration-for finite > 0 . Arriving at = 0 , the system un-jams at u J = and the jamming density remains the actual one. 6 6. The norm of the deviatoric elastic strain is, in accordance to the general scheme, u s = √ u * ij u * ij = (2J u 2 ) 1∕2 . 7. In general, we take t as the partial time derivative, and d dt as the total one, including all convective terms. Hence, with the vorticity tensor given as , one has (as example) the total time derivative of the elastic strain Being off the focus here, the convective and vorticity terms are neglected, so that d dt ≡ t . The dots in ̇p ij and ̇i j are only a (convention preserving) indication of rates, but do not represent the mathematical operation above. 8. The total stress is not an independent state-variable, but rather given by the energy density and entropy production, as discussed in the classical GSH literature.
In the simplified version, it may be written as ij , with elastic, kinetic/granular temperature and viscous contributions. The isotropic stress is referred to as pressure, P = 1 3 kk , the elastic pressure is P = 1 3 kk , for three dimensions D = 3 , and the deviatoric elastic stress is denoted as Note that there are alternative definitions for shear/ deviatoric stress 7 .
9. The symbols B and G are used in the definitions of, respectively, isotropic and deviatoric (shear) elastic energy density, w e = w − w T , defined as the total minus the thermal energy density. In previous GSH-papers [47,48,112], the symbol A was used instead of G, but since A is here referred to as anisotropy, see Ref. [49], we stick to G 8 . 10. The granular temperature used in GSH is T g ∝ √ w T , encompassing both kinetic and potential fluctuating energy contributions. The granular temperature used in kinetic theory and DEM is different, denoted as T G = T K = 2E kin ∕MD , with total mass of all particles, M, in dimension D , ignoring the potential part. Comparing GSH-formulas in the gas limit to those of kinetic theory [7,8,114,116,119], one should remember In the following, we will use T g 9 , with units of velocity, i.e., if scaled by the particle diameter, that of an inverse time, or a rate, very similar to the fluidity, g, studied and discussed in Ref. [18], and references therein 10, 11 .

Overview
In what follows, we shall, in Sect. 2, consider the significance of an inflection surface, of a convex-concave transition in the energy, as relevant for classical systems, transiently 6 Generalizing GSH, we also allow negative "elastic" strains, = e v , interpreting those as the separation between particles-or their mean free path-in order to catch both jammed and un-jammed situations. Note that the elastic energy of a negative is identically zero, and that a negative is not independent of the density . Compressing from an un-jammed state, the system jams at = 0 , towards > 0 and > J . In isochoric situations (constant density), an evolution of the state-variable, , the isotropic elastic strain, implies an evolution of the (enslaved, dependent) jamming density, J = exp(− ) , as proposed and studied in detail in Ref. [49]. The physics clearly changes between positive (jammed) and negative (un-jammed) states, but for the sake of brevity, below jamming, we limit J ≥ J0 and thus Δ ≤ Δ 0 = log( ∕ J0 ) , in cases where it would drop below its absolute limit, J0 , which can be seen as the random loosest packing density. 7 Elastic stress and strain are chosen as energy conjugates such that w e = ij e ij = P + s u s = p e v + q e q , where the last term implies the geo-mechanical definition for elastic shear stress, q = √ 3∕2 s , and strain, e q = √ 2∕3 u s , equally energy conjugate, but different in the pre-factors. Another example, shear stress in simple shear, ignoring out-of-plane anisotropy, = xy ≈ √ 2 s , is conjugate to shear strain, xy ≈ √ 1∕2 s , but not to strain, , integrated in time over shear rate, ̇= 2̇x y , as discussed in Ref. [75]. 8 Note that (calligraphic) symbols B ≠ B , G ≠ G , and A = A ij , in general, are the (tangent) moduli, representing the second derivatives of the elastic energy density with respect to isotropic and deviatoric elastic strains, or mixed, respectively; symbols B , G are again different and are the secant moduli; for more details see Sects. 3.2.2 and 3.2.3. 9 The two temperatures T g and T G are different in the following sense. In thermal equilibrium of a static granular solid, one has T G = 0 , but T g = T , since it is defined as equal to the true temperature in equilibrium, see Eq. (23) and Refs. [47,48,112]. In granular gases, if thermal equilibrium could ever be reached, one would have T G = T , a relevant situation if one starts to consider dissipation and the consequent heating of the grains. By not claiming that T G is a "temperature" of the granular degrees of freedom, taking it only as a measure of the velocity fluctuation squared, T G ∼ | v i | 2 , one may go on using T G in denser ensembles too, ignoring the fluctuations of potential energy it has lost the meaning of a temperature. Conversely, one may use T g in granular gases, anyway, taking it as While T g = T does hold in granular static equilibrium, T G = T can never be reached, as any finite T G , for finite sized particles, translated into temperature, leads to unreasonably large values of order of the inner temperature of the sun. Only in the atomic/molecular limit of "particles" one has T G analogous to k B T . It is therefore more sensible to employ T g throughout. 10 The words elastic/plastic can mostly be interpreted as synonyms for reversible/irreversible in case of vanishing/finite fluidity.. 11 Referring to the states of granular matter, there are many words such as: solid-like, elastic, reversible behavior and fluid-like, plastic, irreversible behavior, where granular matter can also take states in between these extremes, as quantified by g. elastic systems and granular matter. We then present a review of GSH and new constitutive relations based on particle simulations, as well as a minimalist version, in Sect. 3, allowing for analytic solutions in Sect. 4, and numeric calculations to catch some transitions in Sect. 5. A quantitative comparison and calibration with particle simulations is carried out in Sect. 6, before we conclude in Sect. 7.

Equilibrium conditions and dissipative terms
In this section, we first revisit the reason for thermodynamic energy's convexity, and derive the equilibrium conditions for three systems: elastic, transiently elastic and granular media. There is one equilibrium condition for each state-variable, that maximizes its contribution to entropy or, equivalently, minimizes its contribution to energy. Examples for equilibrium conditions are uniform temperatures and stress gradients proportional to density. As these conditions represent extremal points, the energy needs to be convex to be minimal, for the system to be stable. Then we make the general point that every equilibrium condition, if not satisfied, is a dissipative channel that gives rise to a negative/dissipative term in the evolution equation of the associated state variable. As a result, the state-variable relaxes, towards satisfying the condition. In a closed system, all variables will eventually satisfy all their respective conditions, which is the state we called equilibrium.
If the energy is concave, equilibrium conditions represent maxima of the energy with respect to variation of a state-variable 12 . The dissipative terms can thus drive the system away from equilibrium, producing, e.g., non-uniformity in temperature and non-equilibrium stress fields. When this happens, what micromechanical mechanisms it originates from, is necessarily more specific. How the dynamics further evolves depends on the system one considers. In the classical van der Waals theory of the gas-liquid transition, droplet formation is the basic mechanism. In granular media, we propose the following mechanism.
In the stable region, within the cone of Fig 2, the dissipative term in the equation for the elastic strain serves to maintain stress equilibrium. It remains inconspicuous as long as one studies the evolution of stresses close to equilibrium.
Outside the cone, beyond stability, it can force the system to leave equilibrium. Non-equilibrium stresses accelerate grains in varying directions, producing jiggling and thus granular temperature which, in turn, allows the stress to relax, pushing the system back into the convex region. This is what we believe happens in grains at yield and beyond the transition. Setting up a dynamical model for following the system through the transition to different states is the main purpose of this paper.

Classical view on equilibrium states
Consider a system characterized by the state-variables density, , entropy density, s, and elastic strain, with elastic displacement vector (field) U i , relative to the stress-free reference configuration, and a thermodynamic energy density that is a function of all state-variables, w = w( , s, u ij ) [113,120].
Note that, for convenience, we spell out only the smallstrain approximation in Eq. (3) [145][146][147], however, throughout, we work with the (logarithmic) Hencky finite (true) strain definition, due to its more favorable properties: additivity, reversibility, split-ability into isotropic and deviatoric parts, which allows the interpretation of J as the stress-free reference density, and thus allows us to work not only with hard but also with very soft materials.
A textbook proof of energy convexity considers only the entropy as a variable, and involves that the system is connected to a heat bath. A temperature fluctuation (associated to entropy fluctuations) vanishes only if the energy is larger with it than without, which is shown to imply convexity [149].
In a more general consideration, we start with the assumption that the system is stable and has an equilibrium for given values of , s and u ij . Since the elastic stress, ij ≡ − w∕ u ij is symmetric, ij = ji , we may write the total differential of the energy density as: with gravitational potential , chemical potential = (w − )∕ 13 , and temperature T = w∕ s. Varying this energy by with Lagrange parameter L = const.; (iii) forbidding external work, i.e., assuming a closed system: ∮ ij U i dA j = 0 ; and dw = Tds − ij du ij +( + )d , 12 Note that non-local terms in the sense of diffusion of granular temperature are very similar to the "non-local" diffusive evolution equation for fluidity, see the discussion in Sect. 3.1. An essential difference here is that our particles are (strongly) deformable, which is not contained/considered in many other works; we do exclude breakage, however. 13 For standard situations, single species, the gravitational potential energy density is = g i x i + 0 (with arbitrary reference, 0 , e.g., such that = 0 at a free surface), position vector x i and gravitational acceleration vector, g i , that, e.g., has the components g 1 = g 2 = 0 , g 3 = −g [113], for gravity acting in negative vertical direction, or g 1 = sin( ) , g 2 = 0 , g 3 = − cos( ) , for the Cartesian system along a plane tilted by an angle from the horizontal [76].
(iv) using Gauss' theorem 14 , yields with either , s or U i varying independently, while keeping the others constant. The split-up equilibrium (extremal) conditions may now be written as: with the total equilibrium stress eq ij 15 . They represent an energy minimum-and stable equilibrium-only if deviations from them yield an energy increase.

The limit of solid elasticity
Elasticity (reversibility) corresponds to the unique dependence between density and isotropic elastic strain increment, u ii = − ∕ , or, integrated, = −u ii = log( ∕ 0 ) , the logarithmic (true) isotropic elastic strain, where 0 is the density of the stress free reference configuration 16 . This allows to rewrite 17 the integral which now implies, in conjunction with Eq. (5): (or ∇ j ij = g i , for P T = 0 ), with the thermal contribution to stress, P T = − ((w − w e )∕ )∕ (1∕ ) , using volume, i.e., inverse density 1∕ ∝ V , as state-variable, as discussed in more detail in Sect. 3.1.1.
Next, inserting T = T eq + T , ij = eq ij + ij , with ∇ i T eq = 0 and ∇ j eq ij = g i , we require Assuming first u ij ≡ 0 , we may write 2 w = T s = ( T∕ s)( s) 2 > 0 , implying or that the energy w is a convex function of s. As a result, temperature fluctuations will diminish, and the state characterized by a uniform temperature is a stable equilibrium. Conversely, if the energy is concave, 2 w∕ s 2 < 0 , the condition ∇ i T = 0 represents a maximum of energy, and the system is unstable. Any fluctuations in entropy will move it away from uniform temperature. In the case of the van der Waals transition between gas and liquid, a uniform singlephase system is moved to the coexistence of two phases, with different entropy densities, but the same temperature.
Next, assuming s ≡ 0 , as used explicitly below, in Sect. 3.2.2 and 3.2.3, we order the six components of ij and u ij each as a 6-tuple vector, denoted by Greek letters, and require This implies that the 6 × 6 Hessian matrix implying that the elastic energy w e is a convex function of the elastic strain u ij . If there is at least one negative eigenvalue, (11) 2 w e u u = − u has only positive eigenvalues, 15 The total stress ij has more contributions than only the elastic stress ij , i.e., a so-called thermal stress, P T ij = ij − ij . The unrealistic assumption of complete independence between and U i (equivalently, u ij ), would imply ∇ j ij = 0 , a situation where the stressgradient (due to gravity) could be carried completely by the thermal stress (i.e., by the chemical potential contribution to the total energy), as for example in an ideal gas or a fluid without elastic stress. The fully dependent case of elastic solids is discussed next. 16 For elastic solids, the incremental approach can be replaced by another (primed) of the many large strain definitions [145,146], so that − det(u � ij ) = ∕ 0 . Main reasons to use the true (logarithmic) strain include its additivity, reversibilty, and the clean, mechanically reasonable split into isotropic and deviatoric contributions (not straightforward for many of the other definitions, according to Ref. [147], but feasible if properly carried out [146]). 17 Using the incremental definition, the integral becomes: 14 According to Gauss' theorem, the surface integral transforms as: Using the definition of the stress or traction vector, t i = ijnj , the surface integral can be rephrased, ∮ ij U i dA j = ∮ t i U i dA , allowing to add tractions (or point/contact forces) at the surface of V, which would pop up on the right hand side of Eq. (5) but are not used here.
and mass-balance, removes the density gradient, Using the large strain definition improperly, −u � ii = ∕ 0 , would result in g i → 0 g i , while using the true strain definition, −u ii = − det(u ij ) = = log( ∕ 0 ) , yields the same result as the incremental strain, due to = ∕ . A more detailed analysis is in progress, to be published elsewhere.
Footnote 17 (continued) the condition ∇ j ij = g i no longer represents a stable state, because along the associated eigenvector, the energy is a maximum. The system can and will depart from its previously stable elastic state, initially by violating ∇ j ij = g i , typically rendering the stress in non-equilibrium.
To obtain stable, static elastic solutions, one has to solve ∇ j ij = g i for given boundary conditions. This is equivalent to looking for minima of the elastic energy. The solutions are stable if the elastic energy is convex and unstable otherwise.
The more general consideration, including both s and u ij , leads to a 7 × 7 matrix that, for stable equilibria, must possess seven positive eigenvalues.
A complete consideration for elasticity requires also the inclusion of density, , and momentum density v i (conjugate to v i ) as the energy's variables. This, being somewhat more lengthy, would distract from the present concern. The associated equilibrium conditions, with the chemical potential, = w∕ (conjugate to , see Refs. [123,150]), are: All three equations, together, express minimal energy, or maximal entropy.
If any of the equilibrium conditions 18 are not satisfied, dissipative currents appear to counteract: heat diffusion ∼ ∇ i T in the evolution equation for s, viscous stress ∼ v ij in the evolution equation for v i , and a term ∼ ∇ k ik − g i = ∇ k ik − ∇ i , in the equation for the "dissipative displacement rate": Analogous to heat-conductivity, quantifies the strength of the dissipation, carrying units of time divided by density, a mere abbreviation; taking it as a scalar is a simplification. All these terms serve the sole purpose of restoring the respective equilibrium conditions: The dissipative displacement rate, as a necessary result of thermodynamics, has been first recognized in the classical 1972-paper: "The unified hydrodynamic theory for crystals, liquid crystals, and normal fluids", by Martin et al. [148]. It drives the system, boundary conditions permitting, towards a constant stress gradient, ∇ k ik = g i . If on top of this equilibrium state, the stress varies, such as in elastic waves, the dissipative displacement rate, contributes to wave damping. If one concentrates on the evolution of stresses in equilibrium, this term vanishes and is irrelevant. However, if the energy is concave, this term can drive the system away from equilibrium and even can result in instabilities (numerical as well as physical). Writing the elastic stress gradient in the notation of the 6 × 6 matrix, see Eq. (11), as: since ∕ u = ∕ u = 0 , we see that, if the matrix, ∕ u , has a negative eigenvalue, the corresponding term will flip its sign. Instead of keeping maintaining stress equilibrium, it can drive the stress away from equilibrium. This, in turn, accelerates mass points, possibly leading to nonuniform velocities, v i , and thus finite strain rates, v ij ≡ −̇i j . Initially, the stress perturbation will grow along the direction associated with the negative eigenvalue, but for finite times, this is by no means true, as the system will try to move towards a new stable equilibrium state, whatever that is. See the next Sect. 4 and 5 about what happens in granular matter without gradients. Discussing the possibility of a relation to gradient plasticity is beyond the scope of this paper.
Inserting Eq. (15) in the definition of the elastic strain, Eq. (3), reads which seems to suggest that the "dissipative current", ij , is simply the plastic strain rate, ij =̇p ij , which apparently exists even in classical solids if the stress is not in equilibrium. This could be confusing, as it is not related to typical plasticity formulations, and neither is it connected to concepts of plastic potentials or flow functions (see Refs. [32,34,122]). The term plastic strain rate is more appropriate for the other dissipative contributions discussed in the next two sections, on transient elasticity and granular media.
Note that heat diffusion and viscous stress exist in any system, in which entropy and momentum are state-variables: liquids, solids, granular media, irrespective of the microscopic interaction. Same holds for the dissipative term ij , which exists in any system in which the elastic strain is a variable. This is the reason it also exists in granular media. Generally speaking, every dissipative term strives to satisfy its equilibrium condition by changing the value or distribution of the associated state-variable. Equilibrium is achieved if all equilibrium conditions are satisfied, as entropy is then maximal.

Transient elasticity and plasticity
There are many transiently elastic systems in nature. If quickly deformed, they are elastic and capable of restoring their original shape. But this does not happen if the deformation is kept longer; then the deformation is irreversible, plastic. One example are polymeric melts that consist of entangled elastic strands, which elastically deform, but disentangle if given enough time. This leads to a reduction, and eventually vanishing, of the elastic stress. For such systems, the equilibrium condition is: Consequently, Eq. (17) takes the form: with the plastic strain rate now a relaxation term, with a positive coefficient e . Employing essentially this equation, including the convective terms of Eq. (1), a wide range of polymer behavior including shear thinning/thickening and the Weissenberg or rod-climbing effect were reproduced [151,152]. It is noteworthy that the plastic strain rate in the form ̇p ij = − e u ij is a diagonal Onsager term, quantify dissipative currents that drive the system towards the respective equilibrium, in this case of the elastic strain. Hence, off-diagonal ones such as are also permitted (in this case strain-rate causes dissipative currents that drive elastic strain to equilibrium). They will turn out to be useful in granular physics, see below, Eq. (26), and the discussion below it.
The close link, even identity, between transient elasticity and strain relaxation on one hand, and plastic behavior of irreversible shape change on the other, is a useful insight. Similarly useful is the understanding of the difference between elasticity and transient elasticity. For the latter to be in equilibrium, the elastic stress has to vanish, while a constant stress suffices for the former. For verbal clarity, we denote where "plastic equilibrium" is short for "transiently elastic, long-term equilibrium".
There is a further subtlety that we must address here. If the polymer energy depends on both the density and the elastic strain, there are two contributions in the stress: the pressure as given by Eq. (14) and the elastic stress. Then the system may possess an equilibrium pressure even when Eq. (21) holds. However, if the density is not an independent state-variable, implying P ≡ 0 , an equilibrium pressure needs a finite ≡ −u ll to be sustained, and u ij = 0 cannot be the equilibrium condition. Rather, it is given as the vanishing of the deviatoric part, while the trace , not independent from the density, simply follows the dynamics of the density. It does not relax. Note that the relaxation time of and u s need not be the same. If that of is especially long, it may be neglected for certain phenomena, for which the dynamics is governed by ̇p ij = − e u * ij alone. When the system is crossing an inflection surface, the term − e u ij , in Eq. (19) is not affected, and continues to push the elastic strain towards u ij = 0.

Granular matter
GSH was set up in compliance with thermodynamics and conservation laws. Here, we discuss its structural part, necessary if one is to be consistent with the general principles of physics. In Sect. 3, a reduced complete version of GSH is presented, including a few, as simple as possible constitutive choices, which will be employed later to study the jamming and un-jamming dynamics.
Two basic pieces of physics characterize granular media: (1) They have two entropies: s g for the granular degrees of freedom and s for the much more numerous microscopic ones. (2) Depending on circumstances, granular media may be elastic or transiently elastic. Both elastic and plastic equilibria of Eqs. (21) are therefore relevant. However, note that the equilibrium (limit) state is not necessarily ever reached, neither under permanent deformation, nor under free relaxation. In the former case, the system is permanently pulled away from the equilibrium (steady state is not equal to equilibrium), while in the latter, if T g relaxes fast enough, the equilibrium cannot be realized by the other state-variables either.
Including s g as an extra state-variable, with T g ≡ w∕ s g , the equilibrium condition is T = T g , obtained by maximizing ∫ (s + s g )dV ≈ ∫ s dV , where s g ≪ s may be ignored. The equilibrium condition implies that all degrees of freedom, microscopic as well as granular ones, will eventually equilibrate with one another. Furthermore, since for particles of grain size well above molecular size, already for tiny velocity fluctuations (jiggling), one typically has T g ≫ T by many orders of magnitude, ∼ 10 10 , see item 10. in Sect. 1.7, we may set the equilibrium granular temperature to zero, In analogy to the relaxation terms discussed above, the evolution equation for s g must therefore possess a relaxation term ∼ T g , pushing s g towards s g ∝ T g = 0 . This dissipation/ relaxation takes place due to collisions, with rate ∼ T g , in the collisional gas-and fluid-like regime. In addition, analogous to the viscous heating term in the hydrodynamic theory of Newtonian fluids, which transfers kinetic energy into heat, there is a term that transfers kinetic energy into "granular heat", g v 2 s → T g t s g . Therefore, assuming ∇ i T g = 0 , and ignoring other gradients, the evolution equation for granular energy reads with coefficient = (T g ) dependent on T g , and the compressional viscosity neglected, like convective and diffusive terms, for the sake of brevity. To be used in the following, after division by T g and some re-writing, the evolution equation for granular temperature reads: The effective rate of dissipation T * g = T g + T e is discussed in more detail 19

below in Sects. 3.1 and 4.
For given deviatoric (shear) strain rate, v s = |v * ij | = | −̇ * ij | , the steady state solution is given and discussed in Sect. 4.6 in the limit cases 0 ≪ 1 T g and T e ≪ T g : a result known to hold in granular gases 20 , up to moderate densities [7,8,116]. In this case, the system is in the rate-independent elasto-plastic regime, where the granular temperature is proportional to the strain rate. For diminishing T g ≪ T e and 0 ≫ 1 T g , we have an exponential and much faster decay, t T g ∝ −T g , however, also here the steady state granular temperature persists and remains relevant, as T (e) g ≈ (T (ss) g ) 2 ∕T e , see Sect. 4.6. Returning to the elastic strain u ij , we note that granular media are elastic for quiescent grains, T g = 0 , as slopes of sand-piles demonstrate. If the particles "jiggle", T g ≠ 0 , the elastic shear strain and stress will diminish, and eventually vanish: Tapping a vessel of grains (with a finite number) long (and strong) enough results in a flattened granular surface, like in transient elasticity. Combining both conditions of Eqs. (21), the evolution equation for the elastic strain contains different types of plastic strain rates, see also Ref. [21,107] and Eqs. (17), and (20), where the first term on the right, pushing u ij towards the plastic minimum u ij = 0 , operates only for T g ≠ 0 , representing the fluctuation driven plastic strain rate. The second term represents strain-driven plastic deformations. Note that the probabilities p ijkl occur well within the macroscopic, elastically stable regime-involving possibly local events, on the particle scale-and will be simplified 21 using only the respective purely isotropic and deviatoric (shear) plastic deformation probabilities, p v and p s , see Sect. 4.1 22 . The micro-mechanical origins of these probabilities, are not addressed here, rather see Refs. [14, 18-20, 28, 34, 49, 129] and many more references therein. There-among other considerations-it is shown that (finite) granular systems can remain elastic for tiny strain, then have localized plastic events at larger strain with probability increasing, before (global) yield takes place with particular probabilities as cast into a meso-scale, stochastic master-equation approach, in Refs. [50,153,154].
The third term, the dissipative current, ij , depends on the equilibrium condition for the gradient of the elastic stress, see Refs. [47,112], and thus vanishes for stress in equilibrium-or constant stress, for homogeneous/uniform systems in the absence of gravity-as relevant in the following sections. As discussed above, around Eq. (17), the dissipative current, ij , pushes u ij towards the elastic equilibrium in the energetically convex region, and away from it in the concave one, since it changes sign at the transition.

Dynamics at constant shear rate or stress
Equation (26), in addition to the dynamics of T g , Eq. (24), render granular behavior rather more complex than the 19 Preempting the discussion in Sect. 3, to write down the final evolution equation for T g , for reasons detailed in [47,48,112], and partially in Sect. 3, we use: in order to work with parameters that do not depend on T g anymore, and mostly ignoring the Newtonian type viscosity 0 in the following. When inserting s g into Eq. (24) for energy, the time derivative of b is assumed to be small and thus neglected. superposition of behavior from polymers and elastic media. Imposing either a constant shear rate or a constant elastic stress in a polymer melt, Eq. (19), the steady state result is the same, v s = e u s , in either case. This symmetry does not hold for granular media-not even for the simplest case with T e = 0 , and p = 0 . A constant shear rate, v s , where the hat indicates the fact that this quantity is fixed/controlled, with the stationary solution (26), ignoring the p-terms on the r.h.s., leads to a rate-independent evolution equation for u ij that possesses the hypoplastic structure [155], since T g is taking a value proportional to the absolute value of the strain rate. The steady state elastic shear strain is thus: It accounts well for elasto-plastic motion [156], including the approach to the critical state and shear jamming [47,48,157,158].
On the other hand, controlling the stress or, equivalently, holding the elastic strain, û s , constant and inserting the stationary limit of Eq. (26), v (u) s = T gûs , into Eq. (24), yields the relaxation rate: the case when we find T g to relax to zero, pushing the system into a static state, v (u) s → 0. The relaxation rate vanishes (i.e., the relaxation time diverges) as the stress (or elastic strain) approaches the critical value, u c s . With a further increase of u s , the rate flips sign to positive above the critical value, see [47,48], creating an ever increasing strain rate v (u) s . Accordingly, switching from an imposed shear rate (say during an approach to the critical state) to an imposed sub-critical stress will render the system static due to the relaxation of T g , whereas a critical or super-critical stress will create T g and thus accelerate the flow, since v s ∝ T g .

Dynamics in the concave region
Within the cone of Fig. 2, in the energetically convex region, as long as one considers only the evolution of stress close to equilibrium, the dissipative current, , remains close to zero, see also Eq. (26) and the discussion below it. Serving to maintain stress equilibrium, it may simply be neglected. Perturbing the system by a (local) stress, ij , from a static situation, in the convex, stable region, results in a relaxation of the elastic strain, due to the sign of ij . In contrast, in the concave region, because of Eq. (16), the relaxation turns into an explosion, and can drive the stress towards further, stronger non-uniformity.
This accelerates the grains, locally, leading to non-uniform velocities v i and finite strain rates, v ij ≡ −̇i j ≠ 0 . The latter serve as a source for granular heat, see Eq. (24), and create considerable T g , which activates the first plastic term of Eq. (26), which relaxes the stress back into the stable, convex region. Hence, although the imposed perturbation creates a local stress response along the direction associated with the negative eigenvalue initially, it is the stress relaxation back to the convex region that dominates for finite times. If not strong/fast enough, the system will yield or un-jam dynamically. This is one way how GSH accounts for stability and un-jamming dynamics by instability, both mediated by the granular temperature Unfortunately, including the elastic stress-gradient driven plastic strain rate renders Eq. (26) an unstable partial differential equation, the solution of which requires increased technical efforts. This is undesirable in a first, qualitative study, and an approximation scheme may prove useful. We suggest to go on neglecting the elastic dissipative terms, and to add a stress term to Eq. (24), such that T g is directly produced by an elastic stress. The balance equations for s, s g , for the energetically convex region, are given as The equally permissible alternative, was not adopted, because any static ij would then produce T g , leading to its decay. This is not observed in static granular media that can sustain finite stresses indefinitely (see sand-piles)-if not perturbed externally. Yet the reasoning is not valid outside the cone, where static stresses are not stable. Hence we combine Eq. (27) with (30), noting the elastically stable cone. The explicit forms for ijkl and ̄i jkl , which carry the units of inverse stress (compliance) rate, are constitutive choices that will be discussed in the next section in some detail for ijkl , while a very simple model for ̄i jkl will be presented in Sect. 4.9 and studied in Sect. 5.

Second law of thermodynamics
The second law of thermodynamics (the balance of thermal and granular entropy) from Eqs. (27) and (28) can be summarized as R > 0 and R g + T 2 g > 0 (since both represent dissipative, irreversible processes), or combined: noting that the dissipation of granular energy − T 2 g irreversibly enters the thermal energy, thus cancelling itself in ijkl = 0 inside, and ijkl = 0 outside, More general, the total production of entropy can be re-phrased [48,107,109,110] as: where the dissipative displacement rate, Y i , as defined above in Eq. (15), being proportional to ∇ k ik − g i itself, renders the last term positive. Next, more complex expressions are needed for the evolution of the elastic strain, e ij =̇i j −̇p ij , and the a-priori unknown viscous/dissipative stress, D ij , for both of which the isotropic and deviatoric parts can-and are assumed to-evolve independently from each other. As example, where the dissipative stress and the dissipative current are not shown for the sake of brevity, inserting the constitutive relation from Eq. (26), leads to the following split-up: where the dots remind us of ignored terms (all ≥ 0).
The first and second term in Eq. (34) represent entropy production by fluctuation driven relaxation and are always positive. In general, the isotropic term can have a different coefficient, 1 ≠ . In contrast, the third and fourth term are due to plastic (re-arrangements) driven by isotropic and deviatoric strain, respectively. They can be either positive or negative, dependent on the direction of the strain rate 23 . The third term is negative for extension ( ̇v < 0 ), while the fourth term is negative, in particular, at strain reversal. If negative, these contributions must be compensated by positive production terms, as derived next 24 , or the probabilities could be set to zero, as will be discussed in more detail in Sect. 4.
The general approach to construct the viscous/dissipative stress from an inserted plastic strain rate is using the Onsager matrix (also to establish time-inversion symmetry, which is not discussed here), from the appendix in Ref. [107]. After ignoring gradients of temperature and thus heat fluxes, as well as all associated terms, one can transform all tensors into eigen-value form 25 , which also implies a signed variable s 26 . This results in only two independent invariants per state-variable tensor, ignoring the third for the sake of simplicity, yielding a 4 × 4 matrix form to determine plastic strains and dissipative stresses: with 4-tuple vectors of the invariants of the state-variables or their conjugates (pressure P = B and shear stress norm s = ±| ij | = sign( )| ij | = G e s ), with secant moduli assumed constant in the momentary configuration, above jamming (cases below jamming, < 0 , are briefly discussed at the end of this section).
One can interpret the diagonal (off-diagonal) terms as quantifying the dissipative currents that drive the system towards its equilibrium-of the respective (cross-coupled) state-variable. For example, e ss quantifies the dissipative current caused by e s , driving itself to equilibrium, while e vs quantifies the dissipative current caused by e s , driving the volumetric elastic strain, e v , to equilibrium. Other terms, like vv ( ss ), quantify the volumetric (deviatoric) dissipative stresses, i.e., momentum current density, caused by isotropic (shear) strain, driving themselves to For example, with constant volume boundary condition, ̇v = p v = 0 , the term strictly vanishes, while for constant stress (pressure) boundary conditions it can be active due to strain (volume) fluctuations. 24 For more details about possible choices of non-diagonal Onsager coefficients, see Appendix A in Ref. [107]. For example, the term p ijkl̇kl , coupling the plastic strain rate with the strain rate, requires a similar term of opposite sign, coupling to the elastic stress and contributing to dissipation via the viscous-dissipative stress: 25 Any rank-two tensor can be decomposed into (v = volumetric = isotropic, and s = shear = deviatoric) T ij = T v ij + T ŝ0ij + … , with unit tensor ij and (reference) unit deviator in diagonal form: , using the isotropic, first invariant T v = T ii ∕D , ignoring its third invariant, and a signed version of the second inva- [125], where the deviatoric strainrate unit-direction was chosen as reference. In the eigensystem, sign(T) switches sign whenever ̂T ij changes direction relative to the (constant) reference ̂0 ij . Translating back from the eigensystem (subscript s) to the Cartesian (subscript ij) removes the ambiguous sign (±), but is not carried out here explicitly, for the sake of brevity. 26 Note that s = ± (u s ∕ ) = (̇e * ij ∕ )̂0 ij , carries along the magnitude and sign (direction) of the deviatoric elastic strain. In some later examples, e.g., in critical state, the shortened syntax = s will be used since s is always positive in such situations. equilibrium, while cross-terms are ignored here 27 . The offdiagonal terms are non-symmetric, e.g., h vs cross-couples plastic strain (dissipative current) caused by shear strain that drives volumetric elastic strain to equilibrium and, at the same time, cross-couples dissipative shear stress with isotropic elastic stress, involving the modulus B 28 . The diagonal term − T g driving granular temperature to equilibrium, represents a source term for temperature, T, and is not shown here, like all related off-diagonal terms (set to zero).
The plastic strain rates (placeholder for strain-rate minus elastic strain, which is the state-variable) are thus: and while the dissipative stresses are and where the probabilities p v and p s vanish if the strain-rate is zero, i.e., they represent strain-activated dissipative stresses, different from the strain-rate activated viscous stresses, as in the first terms.
For the sake of brevity, some of the above terms are not used further on, i.e., while the 1 -term will be used, only a placeholder is given for sv = ṡv (1 − p v ) ; similarly, for the dissipative stresses, only the first terms are used in some (numerical) solutions (even though small), while the other terms are subject to ongoing studies. The only non-classical term used is the granular energy creation, active outside the elastically stable regime, abbreviated as f g =̄i jkl ij kl , defined in subsection 2.4.2, and used below in Sect. 4.
Inserting the above expressions into Eq. (33) results in an always positive total entropy production: which is true by construction, since all terms are quadratic in either elastic strains (stresses) or strain rates. The first two terms are a simple constitutive choice for the more general forms in Eq. (32), with purely isotropic, v = 1 T g ∕B , and deviatoric, s = T g ∕G , compliance rates.
Below jamming, one has (per definition) no elastic stress with purely plastic deformations, with probabilities, p v = p s = 1 , and thus only the viscous terms survive in Eq. (36). On the other hand, in the ideal elastic limit, above jamming, one has no plastic deformations, p v = p s = 0 , so that the classical GSH shows up with only the -preceded relaxation terms surviving (and the -preceded terms in plastic strains and dissipative stresses cancelling each other).
The question how to split up R and R g , as well as related questions about the discrete nature of plastic (re-arrangement) events, and the principle of time-reversal symmetry are subject of ongoing research [21,159].

Granular solid hydrodynamics (GSH)
As review, GSH is a continuum mechanical theory for granular media, set up in compliance with thermodynamics and conservation laws. GSH possesses the state-variables: g , and (vi) temperature T, not used in the following, conventions/ nomenclature are given in Sect. 1.7.
The GSH used here reduces to various different, more classical theories, in the respective limits-when set appropriately, as was shown in: Refs. [156] for hypoplasticity, [160] for the (I)-rheology, and [47] for fluidity, etc. The question is now if it is possible to catch the complex phenomenology at yielding, jamming, un-jamming, elasticity and loss of elasticity with a simple model that only knows about four state-variables: , , u s , and T g .
For the sake of completeness, we first recollect the more complex, more complete classical GSH, as published in the previous years, in Sect. 3.1, before we reduce GSH to an over-simplified minimal model in Sect. 3.2, which will allow for a better understanding of the structure of GSH. Note that the nomenclature of classical GSH is applied in Sect. 3.1, whereas we switch to the positive compressive strain convention and nomenclature in Sect. 3.2.

About classical GSH
The complete equations of GSH may be found in Refs. [47,112], a simplified version in Ref. [48], from which we boil down to a minimalistic version in Sect. 3.2, ignoring not only 27 The momentum density v i is a state variable, its conjugate variable is the velocity, v i . Vanishing total strain rate is the respective equilibrium condition, hence, total strain rate is the cause, a thermodynamic force like ∇ i T that we ignore here. 28 The appearance of a modulus in some matrix entries in Eq. (35) is just a consequence of using the elastic strain on the right hand side but the dissipative stress on the left; this could be avoided by using elastic stress on the right, which would render the h-terms anti-symmetric. Since we use elastic strain throughout this paper, we stick to this less convenient choice also here.
80 Page 16 of 41 momentum density and gradients, but also the density dependence of most transport coefficients and parameters, since those represent top-down constitutive assumptions, rather than basic (qualitative, bottom-up) theory. First, we discuss a few complications in the classical GSH nomenclature, that are not necessary for our present focus, but will become important if a more quantitative model is the goal, so that we keep them as reference for the sake of completeness.

The classical GSH constitutive model
The energy density has a granular thermal and an elastic part: This represents the first constitutive assumption at the core of classical GSH. In the following, we drop the explicit -dependence of B and G for convenience, but keep in mind and used it whenever needed. (In previous GSH-publications, G was denoted as A.). The elastic stresses are defined as the derivatives of w with respect to the elastic strain u ij : with P ≡ ∕3, which represent no constitutive assumptions, but are just a consequence of Eq. (37). Like the elastic stress, being conjugate to the elastic strain, the granular temperature is conjugate to the granular entropy, which allows to define the thermal pressure, P T , as the derivative of the granular thermal free energy with respect to volume, at constant s g or T g , as: where we note that the granular entropy is not needed, replaced by the density dependent (positive) function b = b( ) , decaying with density, b∕ < 0.
The elastic energy w has been tested for: (1) static stress distributions in silos, sand piles, point loads on a granular sheet [161]; (2) incremental stress-strain relations from varying static stresses [162]; (3) propagation of elastic waves at varying stresses [163]. As already observed in Ref. [112], w is convex if: For more details see Sect. 3.4. Because the macroscopic friction, or yield limit, 0 ∼ √ 2G∕B , is observed to be not (or only weakly) density dependent, in steady state, at least for cohesionless granular media, the next constitutive model assumption used is: G∕B = const. , and where B 0 > 0 is a constant, and ̄≡ 1 9 (20 p − 11 cp ) , with cp − p ≈ p −̄ . ( cp is the random-close packing density, the highest one at which grains may remain uncompressed, p is the random-loose packing density, the lowest one at which grains may stay static.) The expression for B was empirically constructed to account for three granular characteristics: (1) It provides concavity, for any density smaller than < p , and convexity between p and cp , ensuring the stability of elastic solutions in this region.
(2) The density dependence of sound velocities, c (as measured by Hardin and Richart [164]), is well approximated The slow divergence at cp mimics the fact that the system is much stiffer for = cp than at loose packing B( = p ) . Comparing these constitutive assumptions for G and B with particle simulations is subject of ongoing work 29 .
Finally, in 3D, the function b was chosen [47] as: with another small power law, a ≈ 0.1 , such that P T ≈ w T for → 0 , and P T ≈ w T ∕(1 − ∕ cp ) for → cp , limits which reduce b to first or second term, respectively.
The thermal pressure, explicitly given as: 29 To account for the un-jamming transition at the random loose density, p , a density dependence of B was seen as necessary in the classical GSH literature. To account for the virgin consolidation curve, higher order elastic strain terms in the elastic energy were proposed, with density dependent coefficients, see [112,165]. The Coulomb yield could be accounted for with no density dependence, as in Eq. (43). Since our illustrative examples are focused on the latter, hence B is set to constant in Sects. 4 and 5. A quantitative comparison with particle simulation data will show which assumptions or terms are really needed.
defines the (positive, dimensionless) abbreviation for the term in brackets, as set to constant in the following sections, a good approximation only for low densities, where the ideal gas pressure, In the regime of standard kinetic theory being valid, and with T g carrying units of diameter, d, and rate, i.e., inverse time, one has a different, analytically known G SKT = 1 + (D − 1)(1 + r) g 2 ( ) , with restitution coefficient r and pair correlation at contact, g 2 ( ). 30 For more details, in particular for the other transport coefficients, see Refs. [8,116,118,166].
Adding some speculative connection to other works, the function b is qualitatively similar to the density dependence, F = gd∕ v , of the scaled fluidity, g, as reported in Ref. [18]. However, note that the fluidity is based on shear stress and shear strain only, whereas the thermal stress, P T ∕G p = T 2 g = ( v∕d) 2 ∼ (g∕F) 2 is also defined for isotropic deformations, i.e., non-sheared systems. Whether g and F( ) are truly related with T g and G p ( ) , and how exactly, is subject of ongoing research and goes beyond the scope of this paper.

The evolution equations
For completeness, we specify the evolution equations in the classical GSH nomenclature, where we note the sign conventions = e v , u ij = − e ij and v ij = −̇i j , see Sect. 1.7. For the elastic strain one has: with 1 as an off-diagonal Onsager coefficient, see Sect. 2.5, accounting for Reynolds dilatancy.
Mass and momentum conservation read: with the total stress: ij = ij + P T ij − 1 T g v * ij , with viscosity, g = 1 T g .
Finally, the evolution equation for T g , with b as given by Eq. (41) and T * g ≡ T g + 0 ∕ 1 =∶ T g + T e , is given by Eq. (25).
The coefficients 1 , 0 , 1 , 1 , and b are all functions of the state variables, especially the density, which would require many more constitutive assumptions, so that they are over-simplified and taken as constants in the following sections. Alternative energy densities are compared next.

Minimal GSH type model for a material point
At the core of GSH, assuming a homogeneous representative volume, without convection, v i = 0 and gradients, ∇ i (...) = 0 , one has a postulated energy density, with an elastic and a dynamic, kinetic/granular contribution. The total stress is thus not an independent (state) variable, but can be abbreviated as where the five terms represent isotropic and deviatoric elastic stresses, kinetic/granular stress (with an over-simplified G p = 1 , which should depend-at least-on density, see Eq. (46)), and isotropic (v=volumetric) and deviatoric (s=shear) viscous stresses, with viscosities = v and = 1 = s , respectively, where the subscript 1 was used above. Additional dissipative stresses are derived and shown in Sect. 2.5, but are mostly ignored further on. Now, a few versions of the energy density are discussed, before elastic energy stability is considered in the next Sect. 3.4. The isotropic elastic pressure (defined in D dimensions) is: 30 Analytical integration of b SKT = −2 ∫ (G SKT ∕ )d is usually not possible, only for low densities it results in b 0 SKT = −2 log( ∕ 0 ) ; using the jamming density as reference, 0 = J , results in b SKT = 2 log( J ∕ )-below jamming-note the analogy to above jamming. and the deviatoric elastic stress is:

The linear elastic energy
Further differentiation yields the (constant) moduli: B lin = B lin , G lin = 2G lin , and no anisotropy A lin = 0 , which is surely too simple for granular matter.
The anisotropic linear elastic energy density: if > 0 with ̂e ij = e * ij ∕| e * ij | , as proposed in Refs. [91,125,167]. yields: and the deviatoric elastic stress: Further differentiation yields the (constant) moduli: B A = B lin , G A = 2G lin , and anisotropy A A = A lin̂e ij , which is the simplest possible anisotropic elastic model, with crosscoupling between isotropic and deviatoric elastic strains and stresses, as compared to particle simulations and discussed in detail in Refs. [91,125,167]. However, the anisotropy modulus is not constant and thus requires an evolution or state equation by itself, e.g., A lin ∕B lin = F dev , with deviatoric fabric, F dev , as observed from 3D particle simulations in Refs. [31,125].

The non-linear (Hertzian) elastic energy
One can derive the elastic stress ij = w∕ u ij , from the simplest (non-linear) elastic energy density: and w e = 0 if ≤ 0 , with u 2 s = e * ij e * ij , and B = B( ) , G = G( ) carrying the units of stress, while their possible functional dependencies on other state-variables (like density) are not carried along in the rest of this study. Two choices (out of many more) of the density dependence of the energy density (and its coefficients) are discussed below, where appropriate, but in other cases the density dependence is avoided completely in order to learn what the effect of this simplifaction would be. The isotropic elastic pressure (defined in D dimensions) is: and the deviatoric elastic stress is: implicitly defining the ( -dependent) bulk and shear secant moduli B and G , which mimic a linear -or e * ij -dependence of isotropic or deviatoric stress, respectively, not to be confused with the (true) tangent moduli: where the tensor nature of A is often dropped.
In absence of deviatoric elastic strain, u s = 0 , using the proportionality of pressure, P ∝ 3∕2 , this results in a pressure dependence of the moduli The notation details and alternative definitions of the statevariables e v = and u s = | e * ij | = | − u * ij | are given in Sect. 1.7.

The granular linear elastic energy
From particle simulations, using the linear contact model, see Refs. [49,125,168,169] and references therein, the (linear) elastic energy density is complemented by a prefactor, dependent (non-linearly) on the coordination number with positive constants 31 , see Sect. 6, so that: and w C = 0, if ≤ 0 , with u 2 s = e * ij e * ij , and B C , G C carrying the units of stress, while their possible dependencies on other state-variables are ignored in the rest of this study; density-dependent coefficients are discussed in the context of the Hertzian energy density. The nonlinear function is empirically chosen-to leading order-such that elastic pressure is proportional to C and dimensionless functions proportional to powers of elastic strain, 0, 1, and higher, i.e., [1 + O( ) + …] , constructed to remove the lowest order 31 The derivatives of C are: C � = C C 1 C −1 = C (C − C 0 ) −1 , and non-linearity ∝ C , for C = 1∕2 32 . The shear contribution to elastic energy is inspired by reasonable agreement with experimental data [21], and involves an additional coefficient C 2 . Assuming an ideally elastic system, i.e., ∕ = , the isotropic elastic pressure becomes: and the deviatoric elastic granular stress is: qualitatively different from the Hertzian type energy density. The tangent moduli are thus: where the tensor nature of A C , the unit tensor, ̂e ij , is mostly ignored in the following. In absence of deviatoric elastic strain, u s = 0 , and ≪ 1 , using the pressure, P ∕( CB C ) ≈ , results in a pressure dependence of the moduli B C ∝ P 0 , and G C ∝ P , or G C ∕B C ∝ P . For small this results in G C ∕B C ∝ (1 − C 1 2C 0 1∕2 ) , i.e., the shear modulus decreases with increasing compression, relative to the bulk modulus, as is also observed from particle simulations with a linear contact model [49,125], see Sect. 3.5.

Granular Hertzian energy density
The combination (superposition) of the two previous subsections, i.e., non-linear contacts combined with the granular (coordination number dependent) energy, is quite similar (different only in its dependence on and C) to the energy density proposed in Ref. [21], but goes beyond the scope of this study and will be discussed elsewhere.

Simplest GSH equations and discussion
For a material point, in absence of gradients, using t ∼ ∕ t ∼ d∕dt , the evolution of density with strain rate: has no free parameters. Here, positive strain rate corresponds to compression and negative to extension, i.e., density increase and decrease, respectively; density can also be seen as the volume fraction, related to each other by the (constant) material density, i.e., = ∕ p . Later, in Sect. 5, units will be chosen explicitly, such that p = 1 , so that = , as implied from now on.
In the evolution equation for the isotropic elastic strain: the first term couples elastic and total strain together, while the second term is relaxing towards zero 33 -in case of finite T g , with rate 1 T g . The third term can be positive (or negative, e.g., at strain reversal 34 ), and thus works against (or with) the relaxation term. The third equation defines the evolution of the deviatoric (shear) elastic strain where the first term creates deviatoric elastic strain, co-linearly with the strain rate, while the second term relaxes the deviatoric elastic strain, with rate T g . A third (dilatancy) term, sv , analogous to the third term in Eq. (63), is permitted by the Onsager relation, and was previously added for symmetry in Refs. [47,48,125], as discussed in Sect. 2.5, and used in Sect. 4, but not discussed here. 32 Equivalently, to first order, c( ) ≈ exp(− 1 5 C 1 C 0 C ) , asymptotically, more conveniently tends to zero for large . 33 Relaxation of → 0 , at fixed density, , implies that the granular temperature (jiggling) causes the jamming density to relax as J → , in both jammed and un-jammed states, increasing and decreasing, respectively. A decrease (an increase) of the elastic strain, , at fixed density, , corresponds to an increase (a decrease) of the jamming density, J , see Ref. [49]. On the other hand, at fixed confining pressure, P, a jammed system, at finite, but small T g (tapping) will develop to a state such that the elastic pressure, P = P − P T ≈ P , remains constant; relaxation of then corresponds to an increase of density, i.e., compaction. 34 After large strain, one has a positive product, e * ij̇ * ij > 0 , but at strain reversal the same term will be negative, for a while, until the elastic deviatoric strain reverts direction.
The fourth equation represents the evolution of the granular temperature with T * g = T g + T e , as specified in Sect. 4, and the abbreviation for the dissipation rate, R T = 1 ∕( b) = R T0 (1 − r 2 ) , proportional to the energy dissipation factor (1 − r 2 ) , where r is the (effective) restitution coefficient. The energy creation terms are condensed into the tensor function f T (̇i j ) , independent on r, so that one could separate them into two energy creation rates, R T0 f 2 , for shear and volumetric strain rates, respectively, with further terms from Sect. 2.5 ignored here.

Minimal elastic model with two variables
One could decompose the elastic stress and strain tensors into invariants (and their orientations). Under the assumption of fixed and co-linear tensor-eigensystems, and ignoring the third invariant for the sake of brevity, what remains are the isotropic and deviatoric stresses, = {P , s = * s } , and elastic strains, u = { , u s = u * s } , each as 2-tuple vectors, denoted by Greek indices. This provides the criteria for energy minima: Using the (positive) invariants yields the simple 2 × 2 Hessian matrix (for second order elastic work): If it has only positive eigenvalues, the (elastic) energy w e is a convex function of the elastic strain-invariants and u s . With other words, the elastic stability criterion is det(C) = BG − A 2 > 0 , where A is a scalar in the rest of this subsection.

Eigen-values and -vectors at elastic instability
First, we compute the eigen-values and -vectors from the matrix , before we introduce constitutive assumptions about the energy density and discuss those, separately, in the next sub-subsections.

GSH with Hertzian type elastic instability
In the case of a Hertzian type elastic energy density, w e , see Eq. (55), as typically used in the GSH literature [47], one has: B = (3∕2)B 1∕2 − (1∕4)Gu 2 s −3∕2 ≠ B , G = 2G 1∕2 , and A = G −1∕2 u s , i.e., the stability condition, BG − A 2 > 0 , translates to as previously shown in Eq. (12) in Ref. [48], and in Eq. (43) above, for elastic, static systems above jamming, for > 0 . Below jamming, for ≤ 0 , one has w e = 0 and thus trivially det(C) = 0 , while at the point of instability: Using w e in Eq. (55), the non-zero eigenvalue can be re-written as: 2B∕G , while the zero eigenvalue will be more (68) (u s ∕ ) 2 < g 2 e ∶= 2B∕G , 35 Note that eigenvectors are normalized, come with unspecified direction (±), are associated to an eigenvalue (superscript (0) or (1) In other words, considering the shear vs. normal stress space, one could see the limit of elasticity as one possible definition of the maximal (elastic) macroscopic (bulk) friction, defined by the ratio: e ∶= * s ∕P = G u s ∕(B ) , with the limit value taken at the loss of elastic stability: 0 e = √ 2G∕B = 2∕g e . Incremental changes of elastic strain along the eigenvector n (1) , with norm e , result in stress increments, (1) = C 1 en(1) , parallel to the slope 0 e = √ G∕B.

GSH with granular elastic energy instability
In the case of a granular (coordination number dependent) elastic energy density, w C , see Eq. (57), and the respective moduli, see Eq. (60), the stability condition, or a condition for the elastic shear strain: which simplifies in leading order to u e s ∝ √ , very close to jamming, ≪ 1 . This renders w C a choice qualitatively different to the Hertzian energy density, with respect to the elastic instability; note that the term ∝ 1∕2 kicks in quite strongly already for small . Some of the higher order terms are negative and drag down the limit of stability at higher density, i.e., at higher isotropic elastic strain 36 .
In other words, considering the shear vs. normal stress space, one could see the limit of elasticity as one possible definition of the maximal (elastic) macroscopic (bulk) friction, defined by the ratio: e C ∶= * s C ∕P C , with the limit value taken at the loss of elastic stability. e0 C = G C u e s ∕(P C ) . Incremental changes of elastic strain along the eigenvector n (1) , with norm e , result in stress increments, (1) = C 1 en(1) , parallel to the slope. Rather than detailing the elastic instability further, refer to the following subsection, which relates to particle simulations and is qualitatively similar to the present energy density for small , see Fig. 3. Finding an energy density in even better agreement with particle simulations, in particular improving the empirical c( ) , is work in progress.

Anisotropic, elastic-plastic moduli from DEM
In Refs. [49,125,167], an incremental (athermal) elastoplastic evolution model for the isotropic and deviatoric stresses was proposed (in their eigen-system, p and ) as with additional evolution equations for the anisotropy modulus A s (with subscript s indicating its tensor nature) [125,167], with eigen-directions of A s (fabric) and not necessarily identical, and the jamming density [49]-not detailed further here. In order to relate this model to the present GSH based evolution equations, assume (overly simplified, for the sake of clarity) that t p ≈ B t and t ≈ G t e s , to arrive at the evolution equations of the elastic strains: where the first terms represent their elastic and plastic responses, with probabilities for plastic deformations p v (see Sect. 4.1.1) and p s (see Sect. 4.1.2), while the second terms are anisotropy terms, cross-coupling isotropic and deviatoric strain actions and reactions.
In Ref. [125], the elements of the constitutive moduli matrix, , were directly deduced from particle simulations,  36 The condition for finding stability is: and took a form (slightly simplified here by implying that the fabric and the elastic strain are proportional): B = p 0 Z (with the product of volume fraction, , and coordination number, Z, which is a non-linear function of ), G∕B = G 0 ( )(1 − g s (u s ∕ ) 2 ) , and A s ∕B = u s ∕ , with G 0 ( ) = g 0 (1 − exp(− ∕ g )) and constants p 0 , g 0 , g s and g . The (elastic) shear modulus is plotted in Fig. 3, as obtained from the (elastic response only) particle simulation data (no details shown), with the (critical state) calibrated material parameters, p 0 ≈ 0.065 (from above B ), maximal (scaled) shear modulus, g 0 p 0 = 0.036 ± 0.001 , i.e., g 0 ≈ 1∕2 , characteristic elastic isotropic strain, c g = 0.070 ± 0.005 , and the (extrapolated) unjamming density, c J ≈ 0.655 ± 0.001 , in c = log( ∕ c J ) , where the superscript, c, implies the critical state. Note that for this model no energy density is available so far-work in progress.
In the elastic limit case one has p v = p s = 0 , and can identify the cross-term in Eq. (73) with the last term in Eq. From this, the condition for elastic instability, G∕B − (A∕B) 2 = 0 , with A = |A s | , translates to: , which implies a very narrow elastic regime for small , since G 0 ( ) → g 0 ∕ g , for vanishing ∕ g ≪ 1 , so that g e ∝ √ ∕ g . For large The "direction" (in elastic strain invariants) of maximal st ability becomes: n (1) , and with the perpendicular "direction" of maximal in-stability: , a f t e r u s i n g √ G∕B = A∕B = u g s ∕ . This model, derived from frictionless particle simulations [125], thus would result in a non-linear u g s = g e in Fig. 4-different from the other models presented before in this section. However, all models have in common that the stress response to a strain-increment in the direction of the unit-vector n (1) , is parallel to the slope 0 e = √ G∕B in stress space. Further consideration of this particle simulation based constitutive model for stress and fabric-and the question if an additional fabric state variable (tensor) is needed at allgo beyond the scope of the present study, but are subject of ongoing research. Nevertheless, anisotropy cross terms, as those discussed above, will be considered in the next section in some situations.

Special cases
In order to reduce the model complexity, and to understand what the eigen-vectors from the last subsection mean, it is instructive to consider a few simple special cases. Some of these cases will be later studied analytically and numerically.

Fig. 4
Sketch of (strain-driven) deformation cases in the space of the elastic strain invariants, i.e., u s plotted against . The slope g e indicates the elastic instability limit. The numbers at the black arrows indicate the case-number, where dashed, thin lines are continuing the trends outwards into the concave zone (elastically instable). The blue and red arrows give the eigen-vectors of stability n (1) , and instability, n (0) , respectively, where their unspecified direction (±) maybe directed outwards or into the convex zone (elastically stable), dependent on the situation and the boundary conditions. Note that in elastic stress space, the blue and red arrows are parallel and perpendicular to the slope 0 e = * s ∕P = √ G∕B , respectively, for the models presented above.
They represent simplifications that boil down a complicated theoretical framework to a simpler, possibly even transparent form that allows for better understanding and sometimes even for analytical solutions. We propose to apply those special cases to any new theory before one really applies the whole framework. Furthermore, the special cases allow to isolate a few of the terms and possibly calibrate the model parameters one by one. One traditional work on more complex, so-called proportional loading paths is Ref. [37], however, we reduce ourselves to the simplest cases only.
For the rest of this section, we use the stability results from the Hertz-like elastic energy density, as discussed in Sect. 3.4.2. Most of the cases are illustrated schematically in Fig. 4. Except for a few cases, most start from a jammed, elastically stable state with finite initial elastic strains (0) > 0 and u s (0) > 0.
(case 0) Assume the system unjammed, (0) < 0 , and apply a constant compressive strain rate, ̇v = −v ll > 0 . The density and the elastic strain, = log( ∕ J ) , will grow together until the system jams at J , from which on its evolution equation kicks in. It was shown in Refs. [170,171], and earlier works cited therein, that already below jamming, the jamming density (and thus ) depends on the procedure of preparation, in particular on the strain rate and on the granular temperature, however, this fact goes beyond the present focus and is thus not studied further.
(case 1) Assuming a purely isotropic de-compression, ̇v = −v ll < 0 , from a jammed state, one expects the elastic isotropic strain, , to decrease faster than its deviatoric (shear) counterpart, u s , until at u 2 s = (2B∕G) 2 , or u s = g e , the looser system cannot sustain the shear-stress anymore, so that un-jamming due to instability with respect to shear occurs. In order to remain at least marginally stable, one needs a decrease of u s → u 0 s = g e , a situation that could be referred to as shear-yielding [10,49,100].
(case 1b) In the situation without initial elastic shear strain, u s (0) = 0 , the stability criterion is always true and the system remains stable until isotropic un-jamming takes place at = 0.
(case 2) In the case of isotropic compression, the model remains stable, unless the virgin consolidation line is reached, where the system restructures to be able to carry the increasing stress. This situation is not detailed further, but the (case 2b) of loading without shear is studied in Sect. 4.1 to display the role of plastic deformations due to isotropic compression.
(case 3) Assuming a purely deviatoric (volume conserving) shear strain rate, ̇ * ij = −v * ij , from a jammed state, one expects the elastic deviatoric (shear) strain, u s , to increase faster than its isotropic counterpart, , could build up, until at → 0 (1 + ) = u s ∕g e , the system cannot sustain pressure (isotropic stress) anymore, so that an instability with respect to volume change occurs, and one has a consequent decrease of J , due to a possible further increase of , i.e., one origin of dilatancy. The evolution of inside the cone and at the limit of elastic stability are qualitatively different, as will be studied numerically later on.
(case 3b) Under the same purely deviatoric deformation, the isotropic elastic strain could also decrease, corresponding to J increasing. This leads to elastic instability at smaller elastic strains, → 0 (1 − ) = u s ∕g e , not much changing the considerations in case 3, but rather leading to compactancy instead of dilatancy, as to be expected for very loose packings. Furthermore, this could lead to isotropic unjamming, if J drops below .
Several of the cases discussed above will be next studied analytically (as far as possible) and numerically in Sect. 5.

Analytical results for special cases
After a summary of the equations that will be used in this section, we take several special cases, starting from the athermal limit, T g = 0 to highlight the probabilities for strain-induced irreversible deformations in Sect. 4.1. Various versions and limits of the (over-simplified) classical GSH model from Sects. 3.2.2 and 3.4.2, are discussed and analytically treated (in some special cases where this is possible), with focus on the effect of T g , while also several new terms and regularization schemes are proposed, to be later used in the numerical solutions. The set of model equations is summarized here for reference, with extensions from classical GSH (model 0): before some meaningful special cases (purely isotropic and deviatoric loading) are discussed below, for which analytical solutions are provided, if possible.
The elastic strains, for the sake of brevity, are to be read as deviations from the equilibrium, i.e., as: → − eq and e * ij → e * ij − e * ij eq , accounting for the stress equilibrium condition, ∇ j eq ij = g i , and other boundary conditions. Avoiding the terms not present in the original Eqs. (62)- (65), one arrives at what we refer to model 0 (having thus no valid athermal limit), which is used as starting point to study transients in Sect. 5. The terms p g and sv = ṡv (1 − p v ) , as defined in Sect. 2.5, are introduced here as place-holders for elements discussed below, in Sect. 4.7, or to be added in future, as introduced in Refs. [49,125,167]. The rate of cooling is modified in the elastic, jammed state ( > 0 ) by adding an "elastic dissipation rate" T e , referred to as model e, as T * g ∕T g = 1 + T e ∕T g = 1 + T e0 h ∕T g where only the special case h = 0 , i.e., T e = T e0 , will be treated below 37 . The presence of T e does not affect the dynamics too much for finite T g ≫ T e (Haff's free homogeneous cooling state (HCS), well below jamming), but in the limit of very small T g → 0 , for elastic, jammed systems, this (phonon/wave-driven) dissipation becomes important, generalizing HCS, providing an exponential decay of T g → 0 in absence of other driving mechanisms (and > 0 ), see Sect. 4.4. The dissipated granular energy is balanced by an equal positive term creating thermal energy or entropy, i.e., this term does not contribute to the second law of thermodynamics, see Sect. 2.5.
The new term f g (g * ) , in Eq. (79), is only active if the system is outside of the elastically stable regime, where g * = u s ∕ − g e > 0 , with the limit of elastic stability g e . It is continuous, inactive in the convex region, active outside. This term generates more granular temperature, jiggling, the more the system gets elastically unstable, due to concavity of the elastic energy.
The terms 1 − p v and 1 − p s represent the new probabilities for elastic deformations, with p v and p s the new probabilities for isotropic/deviatoric plastic deformations, respectively, see Ref. [49], as specified in Sect. 2.4, Eq. (26), and discussed next, in Sect. 4.1.

The granular athermal limit T g = 0
Enforcing the athermal case, T g = 0 , the system of equations reduces to: see Eqs. (73) and (74), where the cross-term sv is usually neglected, and the off-diagonal Onsager coefficients p v and p s were previously introduced in Ref. [112], but taken equal to 1 , while here they are, alternatively, interpreted as the probabilities for (isotropic and deviatoric) plastic (restructuring) events in the packing, as in Sect. 2.4, Eq. (74), and in Refs. [49,167]. Note that in Eqs. (80) and (81), the probabilites for isotropic and deviatoric plastic deformations are systematically attached to isotropic and deviatoric strain rates, respectively.
In the few plots in this section, dimensionless units are used, such that = , and stress is in units of 1 MPa, as discussed in detail at the beginning of Sect. 5. The quantitative calibration of the GSH based theory by particle simulation data, as well as an alternative choice of units, are discussed in Sect. 6.

Athermal isotropic loading
For isotropic loading ( ̇ * ij = 0 ), the system of equations reduces even further to The elastic limit, with probability p v = 0 , translates to constant J , whereas the fully plastic limit, p v = 1 , translates to ̇p v =̇∕ =̇v.
Relating this to the classical cam-clay model [122], where ̇p v =̇v ( − )∕ , allows to identify p v = 1 − ∕ as constant. However, in the following we will derive a more complex model, where p v , the probability for plastic deformations, is a function of the state-variables and the sign of deformation rate (i.e., compression or tension).
A simple constitutive assumption, p v̇v = − 1 T e , could be directly merged into the relaxation term as − 1 T * g , with T * g = T g + T e , and solved analytically 39 . This model displays the transient elastic behavior of polymer melts or glasses for which (in absence of any isotropic strain rate, for finite, constant T e ) → 0 . However, since the reality of granular matter, as measured from particle simulations in Ref. [49], is somewhat more complex, already for frictionless spheres-and even more for realistic frictional non-spherical particles-we have to come up with a better relation for the probability for isotropic plastic rearrangements.
The (un-)jamming density was reported, see Eq. (5) in Ref. [49], to reach after infinitely many isotropic loading/ unloading cycles the limit density: 37 For a Hertzian type bulk modulus, the time-scale of momentum (wave) propagation, for u s = 0 , can be estimated as , an exponent h = 1∕4 . This estimate, together with a Hertzian elastic pressure, P ∝ 3∕2 , yields an estimated wave speed v e ∝ P 1∕6 or bulk modulus B ∝ P 1∕3 . 38 The chain rule yields an identity between the plastic strain rate and the time-evolution of the jamming density: 39 Inserting the expression from above, this yields the athermal evolution of the elastic strain: t =̇v − 1 T e0 1+h . with the half-sided linear function [x > 0] + = x , and [x ≤ 0] + = 0 , otherwise, ∞ = J0 for < J0 40 . Using Eq. (82), and noting that the first loading is characterized by , with p v0 = 1 − e −1 , see Eq. (6) in Ref. [49], one can deduce the probability for plastic deformations from: yielding: expressed in terms of the difference between t h e l i m i t a n d t h e a c t u a l j a m m i n g d e n s i t y Cast into a pseudo-code, as implemented for the numerical solutions, see Figs. 5 and 6, the isotropic plastic strainrate is computed from Eq. (84): with a regularizing step-function smoothed around = 0 , with = 2 × 10 −5 . In words, the different terms/factors mean: First, the step function takes care that below jamming, one has no strain-driven plastic deformations.  = 0.08 . The dashed green curve represents the initial loading, up to max (green dots), with six un-/re-loading cycles, ending at the magenta dots. Note that the lowest max = 0.62 is un-jamming and re-jamming during the cycles. The upper (blue) curve represents the elastic limit case, with p v = 0 , i.e., without plastic rearrangements and the analytical pressure state-line: P = B 3∕2 , with B = 1 . The lowermost curves represent cyclic un-/re-loading from max = 0.68 with large amplitude, = 0.08 , down to min = 0.60 , well below the jamming-point, un-jamming and re-jamming during every cycle. The inset represents the void fraction, e, plotted against (logarithmic) pressure, P, similar to Fig. 2a 40 Thus, while J0 ∼ p , corresponds to random loose packing, ∞ ∼ cp takes the place of random close packing, cp , continuously growing with density. The higher densities could be achieved by overcompression of soft particles (rubber, gel, etc.), whereas hard particles (metal, glass, etc.) would break (not considered here). For hard/ rigid particles, one could replace Eq. (82) with an interpolation function between p and cp , with pre-factor b ∞ = cp − p , for ≥ 0 , more similar to ideas in Ref. [21]. (ii) Just above jamming, the denominator in Eq. (84) leads to a divergence, but is limited to p v ≤ 1. (iii) Second, in cases where the probability would become negative, typically for large J > ∞ , the system is set to be reversible, with positive probability, p v = −p v , and limited by a small constant slope, p − v ≥ 0. (iv) For positive, p v , the system is partly plastic during loading, but set to be elastic for unloading. (v) If the jamming density is below its isotropic lower limit, J0 , the system is set to be perfectly plastic. (vi) The last line makes sure that there is no strain-driven contribution to dissipative stress if static.
For some more details, see Sect. 6.2, and for the statistical analysis of the probability, see Refs. [18,19,159], and references therein. For experimental calibration of the model, and its relation to classical models, see Refs. [21,122].

Athermal deviatoric (pure shear) loading
Another special case that allows for analytical treatment is pure deviatoric (isochoric) shear, ̇v = 0 , the elastic strains develop as t = 1 e * ij̇ * ij 1 − p s and t e * ij =̇ * ij 1 − p s or, equivalently, for the plastic strain rate Postulating the existence of a constant "critical" steady state for the (elastic) macroscopic friction, ≈ e = s ∕P , i.e., the quasi static limit stress ratio, 41 ,42 , this allows to express the probability for plastic (shear) events as: for > 0 , and p s = 1 otherwise, where the hats denote unit-tensors, the ratio of the tangent moduli depends on the constitutive choice of the energy density, and the last approximation is only valid after sufficiently long steady shear, close to the critical state, and/or if ≈ c vanishes from Eq. (86), but not for strain reversal, similar to Eqs. (28) and (31) in Ref. [47]. The term in brackets limits p s ≥ 0 , as to keep it positive, i.e., [x > 0] + = x , and [x ≤ 0] + = 0 , with p s = 0 in the absence of shear, v s = 0 , and thus valid also for strain-reversal and during early transients, for which negative argument values result in perfectly elastic response, p s = 0 , as done similarly in Refs. [49,125,167] and references therein-based on, and in agreement with, DEM simulations 43 . The probability for plastic events in Eq. (81), specified above in Eq. (86), can be very small at the beginning of shear, but increases due to the build-up of elastic shear strain, u s , before it asymptotically approaches p s = 1 , for large strain in the perfectly plastic, critical state. At reversal of shear, the argument of the bracket-function becomes negative, i.e., the system is elastic with p s = 0 , until the shear strain adjusts to the new direction 44 . When stopping shear, v s = 0 , the system also becomes elastic, p s = 0. In cases where density and thus is reduced in magnitude after the system has reached the critical state (not studied here), also u c s ∝ c will reduce, and terms with (1 − p s ) could become negative, resulting in the decay of elastic strains, however, this is skipped here for the sake of brevity.
Analytical treatment is possible close to steady state, for ≈ c assumed constant, where the deviatoric elastic strain evolves as: t u s = v s 1 − p s ≈ v s (1 − u s ∕u c s ) , with analytical solution: and critical state elastic shear strain: as plotted in the inset of Fig. 7 as dashed lines, for the initial shear stress evolution dev = * s = 2G 1∕2 u s . This analytical solution is very similar to the solutions presented in Refs. [49,125,167], however, further discussion is beyond the scope of this paper 45 .
What remains is to also consider the variation/evolution of during pure (volume conserving) shear, due to variations in the jamming density. Noting the similarity between p s and the 1 -term in Eq. (80), in order to solve the problem analytically, one can rewrite the evolution equation for the isotropic elastic strain as: 41 For the Hertzian energy density, see Eq. (55), using ∶= * ∕P , the ratio of moduli, , between shear and bulk modulus, and allows to determine from the quadratic equation: c 0 (u c s ∕ c ) 2 − 4u c s ∕ c + c 0 g 2 e = 0 the shear to isotropic elastic strain ratio u c with real solutions for c 0 ≤ 2∕g e , as realized in cases modeled here (data not shown). 42 For the granular energy density, see Eq. (57), the ratio of tangent moduli and the solution for c 0 are not spelled out here, for the sake of brevity. 43 If one can assume: ≈ c , i.e., that the isotropic elastic strain is almost constant, close to its critical state limit already, Eq. (81) can be solved analytically, yielding an exponential approach of u s to its critical state limit, see Ref. [49]. 44 Like for p v , this could be merged into the relaxation term − T * g u * ij , if one would assume: −v * ij p s = − T e u * ij , the discussion of which goes far beyond this paper. 45 Note that since u c s depends (weakly) on , the system of equations is still coupled and the analytical solution is only approximate.
s v s p s 1 − p s ≈ 1 u c s p s t u s ≈ 1 u s t u s , for constant v s (not valid for strain-reversal). This equation has a critical state solution, = c , due to the term 1 − p s ≈ 0 , as well as a stable elastic solution with t ≈ 0 for p s ∝ u s ≈ 0 , in an initial isotropic state, see the infinite slopes in Fig. 6, for small shear strain and thus small shear stress. Due to the quadratic proportionality to u s , the variation in is much smaller than the variation in u s itselfat least as long as u s ∕ < 1 . During shear starting from an isotropic state, the variation in is strictly positive, in the athermal limit, corresponding to dilatancy, i.e., the jamming density decreases, while at shear strain-reversal the evolution is opposite, changing sign during evolution, allowing for "butterfly-shape" loading-unloading cycles of or P , see the inset in Fig. 8, as consistent with particle simulations [125]. Note that the above assumptions are reasonable well above jamming, > 0 , but not close to jamming, where the un-jamming can happen during every branch of the cycles, as shown in Ref. [49].

The granular thermal limit
Assuming that one could maintain a constant granular temperature steady state, e.g., by homogeneous driving/tapping, see Refs. [172], this would result in the set of equations: For vanishing strain rate ̇i j = 0 , the equations decouple and only the relaxation terms survive, This corresponds to the "plastic equilibrium" limit case = 0 , e * ij = 0 , which is approached exponentially fast, with rates 1 T g and T g . The term p g = 1 allows to choose the plastic equilibrium of transiently elastic systems, for which → 0 , or is needed in a form p g = 1 − ∞ ∕ , or p g = 1 − p v , so that a granular type plastic limit with > 0 can be achieved, see Sect. 4.7.
For finite ̇i j , the system will establish thermal, elastoplastic dynamic states that are not discussed further for the sake of brevity.
Strictly controlling density, i.e., fixing e, the situation is interesting again for granular matter. Any perturbation, as tapping or small-amplitude cyclic shear, will typically result in a decrease of both the elastic strain, , and consequently the pressure, P = B . In this situation, the pressure curve shifts to smaller densities (larger e), and changes slope, both moving it away further from the elastic state-line, like shown in Fig. 6 (which represents a zoom into the previous Fig. 5, , is chosen smaller than the elastic stability limit, 2∕g e = 1 (dashed blue line), such that the latter is never reached. The inset represents the shear stress evolution with strain, during the cyclic forward-backward shearing, where the higher density cases reach larger stress levels; the thick dashed lines represent the analytical solutions from Eq. (87) during initial shear from the isotropic states Fig. 8 Pressure plotted against density, from the same model solution as in Fig. 6. The lower curve represents the initial loading, up to max (green dots), with six cyclic forward-backward shear cycles, ending at the magenta dots, displaying the pressure-dilatancy caused by isochoric shear. The upper (dashed blue) curve represents the elastic limit compression, with p v = 0 . The inset displays the pressure during the shear cycles for the largest two densities but with different parameter J0 = 0.60 , which would more resemble frictional granular material. On the other hand, large strain shear results in (pressure) dilatancy, shifting the state-line to higher densities (or the void fraction to the right, towards the VCL, but not beyond), defining the critical state line (CS)-see magenta points in Fig. 8. The amplitude of pressure dilatancy increases with density (pressure), and so does the characteristic shear strain at which the system transits into a new state.

Isotropic jamming/un-jamming in minimal GSH
The model equations for isotropic compression/tension, with strain rates t =̇v ≠ 0 , and ̇ * ij = 0 , reduce to: The density is coupled to strain rate directly, while the second equation (91) is decoupled (for sv = 0 ) just relaxing an existing elastic shear strain to zero. The coupled evolution equations (90) and (92) could be (quantitatively) calibrated to particle simulations like in Ref. [49], in a future study, however, in Sect. 6 they are calibrated to the athermal case, that isolates the evolution of , as is relevant also for extremely small compression rates, ̇v , and thus T g ≈ 0 , if p v ≠ 0 . For finite positive (compressive) strain rate, the inhomogeneous solution leads to a divergent increase of T g with time due to the continuous energy input. The energy production term due to elastic instability in Eq. (92) would become active for finite u s , when < g e u s , but is ignored here, assuming u s = 0 (which is not strictly true in real systems, where there can be some small, local, random elastic deviatoric strain).
For finite positive (compressive) strain rate, one has a continuous energy input due to the viscous source term f T , that can lead to increase or decrease of T g , and thus affects also the evolution of . For negative (expansive) strain rate, the same is true, however, as soon as the system approaches unjamming, the behavior qualitatively changes due to T e → 0 , which is qualitatively, not quantitatively accounted for in the present version with constant parameters, in particular f v and R T0 ; more details are beyond the scope of this study.

Homogeneous cooling below and above jamming
In the absence of any strain-rate mode, or other means of energy input [172], and assuming that T g is so small that is practically constant, the evolution equation for T g , abbreviating = R T = R T0 (1 − r 2 ) , and assuming T e = 0 , results in an algebraic evolution: in the free, homogeneous cooling state, as relevant for systems below jamming in the granular gas state 46 . On the other hand, assuming the simplest model for T * g ≈ T e , with h = 0 (or for constant ), for a small perturbation from an elastic base state, one has as relevant for elastically stable systems, well above the jamming density, for which small perturbations decay exponentially fast.

Pure shear transients from an isotropic state
This case was studied in detail by particle simulations in Refs. [49,125], and should be studied analytically too with respect to questions about the build-up of anisotropy, and the degradation of the (shear) modulus, but is skipped for the sake of brevity.

Steady state pure shear (model 0 and e)
In case of deviatoric pure shear, the density equation vanishes, since v ll = 0 the density is conserved, t = 0 , and the terms with isotropic strain rate in the equations drop out. The remaining equations yield the steady state solution for the granular temperature: with T * g = T e + T g , so that (for T e = 0): 46 Remember that the granular temperature in the standard kinetic theory literature is T G ∝ T 2 g , but we do not consider all those details here and rather refer to the relevant literature, e.g., Refs. [7,8,115,116] and references therein.
where only the positive solution is reasonable.
In the "collisional" limit T g ≫ T e , one has the dynamic steady state: T (ss) g ≈ T (ss) g0 ∝ v s , while for T g ≪ T e , the steady state temperature in the "elastic" steady state is: , it vanishes quadratically for v s → 0.
For the deviatoric elastic strain one has: so that: while for the isotropic elastic strain one has: so that inserting Eqs. (95) and (97) yields the isotropic elastic strain in steady state: the former valid for model e, the latter for the simplest model 0, where the subscript 0 indicates T e = 0 ; model e is not indicated since it represents the default case. In the "elastic" limit T g ≪ T e , for v s → 0 , the other two state-variables, in model e, behave as: u (ss) s , and thus g (ss) = u s ∕ → v s , i.e., a leading order linear increase with (shear) strain rate.

Steady state pure shear (model 1)
In model 1, only the evolution equation of the isotropic elastic strain has to be modified: so that inserting Eqs. (95) and (97) yields the isotropic elastic strain in steady state: for model 1 for constant or -independent p g .
In some of the numerical implementations, we used p g = − ∞ , in order to make relax towards a finite value, with ∞ = log( ∞ ∕ ) , as defined in Eq. (82). This allows to re-write p g = log( J ∕ ∞ ) , which makes the relaxation term vanish for J = ∞ , negative for larger values and (96) increasingly positive for smaller jamming densities. Unfortunately, it also requires to solve a quadratic equation, resulting in i.e., an increased steady state elastic strain, representing strain-dilatancy. Note that this approach to achieve finite under steady state shear, increasing with density-as to be expected-is different in philosophy than making the bulk modulus factor B density dependent.
In other of the numerical implementations, we used p g = 1 − ∞ ∕ , in order to make relax towards a finite value, resulting in the simpler steady state expression: with ∞ < 0 for > ∞ , which even can change sign dependent on the relative magnitudes of (ss) and ∞ . Note that this approach to achieve finite under steady state shear, increasing with density-as to be expected-is different in philosophy than making the bulk modulus factor B density dependent.

Discussion of the steady state rheology
Dividing Eq. (97) by (98) yields the deviatoric to elastic strain ratio in steady state (in order to evaluate whether the system is elastically stable or not): If the ratio of elastic strains in Eq. (100) is smaller than the elastic stability limit g (ss) ≤ g e = √ 2B∕G the system remains in a possibly stable (elastic, jammed) state, while it looses stability if the ratio reaches and/or exceeds the limit value.
Solving numerically the system of equations, including the transient evolution, confirms that the steady state is independent of the density, for model 0, see Sect. 5, as does not appear in the steady state solutions above.
The elastic strain ratio, Eq. (100), which determines whether the system becomes elastically instable in steady state, is not the same as the macroscopic friction at which the material flows plastically. Dividing the steady state shear stress by pressure defines the macroscopic (bulk) "friction": = * ij ∕P , which results in the steady state bulk friction: In the slow strain rate limit, ̇i j → 0 , of Eq. (101), above jamming, > 0 , the second terms in nominator and denominator vanish, linearly and quadratically with T g → 0 , respectively, and one has For the special case g (ss) = g e , when the elastic limit of stability and the steady state ratio of elastic strains coincide, this translates to: (ss) 0 = 2∕g e .

Temperature unjam-regularization (model g)
In order to regularize the elastic unjamming instability, we introduce a measure for the distance from the elastic limit g s = (g − g e ) = (u s ∕ − √ 2B∕G) , which influences the temperature evolution with f g (g * ) = f g (g s )g s , and the step-function (g s > 0) = 1 , and 0 else, so that one has for steady-state pure shear (with model 0): i.e., just an elevated granular temperature that affects, in turn, the other state-variables (elastic strains) via their respective relaxation terms, as will be shown in the next Sect. 5.

Numerical solutions of GSH
In order to better understand GSH, including transients and transitions, and to validate the analytical solutions in the previous Sect. 4, we solve the system of equations numerically (with matlab, using ode45) and-focusing on a few terms only-discuss the features of the simplest GSH type model with mostly constant coefficients, see table 1, and the energy density from Eqs. (55). Symbols with a prime are dimensional, whereas all presented results are dimensionless (without prime), as explained next.
volume, V ′ p , of a single particle, so that its mass is: The protocol of the numerical solutions consists of three stages: The initial preparation by isotropic compression (green) is followed by the testing mode (various colors for different parameters), and finally by a relaxation phase without any strain rate (magenta). The testing mode is in the following examples pure deviatoric (volume-conserving) shear, for large strains, to approach the critical state.

Effect of elastic dissipation and unjamming
Next goal is to understand the behavior of the simplest version of the classical GSH model, see Sects. 3.2.2 and 3.4.2, and the effects of both elastic dissipation parameter, T e , and temperature regularization, f g , that controls the dynamics at elastic unjamming.
The initial preparation starts from an un-jammed state at (0) = 0.58 , with isotropic jamming taking place at density J0 = 0.60 , up to different target densities = 0.61 , 0.62, 0.63, 0.68, 0.74, and 0.80 during t p = 1000 . From this point on, pure shear is applied for t s = 5000 and the final relaxation is applied for t r = 4000.
First, the effect of T e on the system is studied in Figs. 9 and 10, by plotting shear stress against pressure and the ratio of the deviatoric-to-isotropic elastic strains against time. Due to the density-independent parameters, in particular B, all the different density configurations approach to the same steady state, as analytically predicted (solid point in upper panel). The overshoot in the transient decreases with increasing density, before the steady state ratio of u s ∕ is reached, and the relaxation kicks in after shear is stopped. In the former, Fig. 9, T e = 10 −6 (case A) is practically zero and has no effect, whereas in the latter, Fig. 10, the finite T e (case B) causes a reduced T g in steady state, as well as a much more rapid relaxation (exponential due to T e , instead of algebraic, like in the free cooling granular gas) to the static state (shorter magenta lines). Due the decreased T g , in steady state, the other state-variables and u s are increased, whereas their ratio is slightly decreased, see Eq. (100).
The effect of the new temperature production term with f g = 4.10 −4 is then tested in Fig. 11 (case C), with otherwise the same settings as in case B. Only those cases that overshoot g e are affected. One of them, the lowermost density case, is completely destabilized by the increase in T g in the unstable regime, reaching a completely different steady state (far out of plot range), and returning rapidly to elastic stability as soon as the shear strain is stopped. Another case (second lowest density) remains above, but moves closer to g e and remains there for some longer time before it reaches the analytically predicted steady state. This proofs that the production term of T g , due to the elastic instability, allows to regularize the systems behavior by dynamic means: Counterintuitively, an increased generation of T g can keep the system closer to the elastic instability, however, if too much T g is produced, this destabilizes the system and allows it to explore the plastic, collisional steady state with very large T g and-at the same time-comparatively small u s and (see the lower left corner in the upper panel and the out-ofbounds data in the lower panel).

Effect of dilatancy and dynamics
Next goal is to understand the behavior of the model at constant density, with different dilatancy parameters, 1 , and the effects of the elastic dissipation parameter T e and the temperature regularization f g .
The initial preparation starts from an un-jammed state at (0) = 0.58 , and is applied up to target density = 0.65 , during t p = 1000 . From this point on, pure shear is applied for t s = 5000 and the final relaxation is applied for t r = 4000 , like before.
The values of 1 are chosen such that a few of the data remain within the elastic instability limit u s ∕ < g e , but a few overshoot, as can be seen in the lower panels of Figs. 12, 13, and 14.
First, the effect of f g on the system is studied in Figs. 12, 13, and later the effect of T e in Fig. 14. Again, shear stress is plotted against pressure and the ratio of the deviatoric-toisotropic elastic strains is plotted against time. In the former, Fig. 12, T e and f g are practically zero and have no effect at all, but an increasing dilatancy parameter, 1 causes the system into decreasing levels of g = u s ∕ during shear steady state (ss). The two lowermost curves remain within the elastic instability limit, the intermediate value 1 = 1.25 displays a slight overshoot but hits g e = 2 in steady state, wheras the upper two curves are clearly beyond the elastically stable regime g > g e . In the shear stress to normal stress plot, the different 1 values lead to different steady states (thick dots) and a slow relaxation (magenta lines). and deviatoric-to-isotropic elastic strain ratio plotted against time (bottom). The green lines (on the horizontal axis) represent the isotropic preparation, the magenta lines (overlapping) the final relaxation, and the big solid (cyan) dot, or dashed cyan line in the lower panel, show the theoretically predicted steady state dev = (ss) 0 p , Eq. (102), being density-independent for the over-simplified models. The dash-dotted (blue) lines represent the elastic stability limit g e , Eqs. (43) and (68) When temperature regularization is active in Fig. 13, the curves in the stable regime are not affected, the intermediate case is slightly modified and the upper two curves (smaller two 1 ) are, again, considerably affected by the generation of T g , i.e., the much larger T g causes both elastic strains to relax towards the plastic limit-see the curves in the lower left corner of the shear to normal stress plot.
In the last Fig. 14, the finite T e causes a reduced T g in steady state, which results also in smaller u (ss) s ∕ , see Eq. (100). During final relaxation, T e is also causing a much more rapid (exponential) relaxation to the static state (shorter magenta lines).
Note that the elastic dissipation term, with finite T e , is reducing granular temperature within and outside, whereas the thermal activation, f g , increases T g , but only outside the elastically stable regime.

DEM particle simulations
The particle simulations to be compared to the GSH solutions are the simplest possible element tests in a periodic cubical cell, with only diagonal components of the strainrate active (isotropic compression/tension and pure shear). The N = 4913 frictionless particles ( p = 0 ), with particle diameters drawn from a random homogeneous size distribution with maximum to minimum width, d max p ∕d min p = 3 , are the same as used in Ref. [49], even though the simulations were re-run slower for the first compression and de-compression cycle, see Fig. 1.

Non-dimensionalization of DEM
The parameters given in the following with a prime, e.g.,  Fig. 9, with difference the elastic dissipation and the granular temperature creation both active. Meaning of colors, lines is the same as in Fig. 9 47 Units could be, e.g., length, x

Calibration of GSH with DEM
The energy density has the same units as stress, with the particle volume, V � p = V � ∕N = ( ∕6)⟨(d � p ) 3 ⟩ , and the contact number per particle C = 2M∕N , given the total number of contacts, M. The fraction of rattlers, f r , that relates C to the coordination number Z = C∕(1 − f r ) , is not studied here, as it was discussed in detail in Ref. [49] and references therein.
The dimensionless energy density is thus:  Fig. 12, where all parameters are the same, except for f g = 5.10 −5 , which determines the granular temperature production in the instable regime, see Eq. (104) (not shown explicitly), which causes the different behavior of the upper curves. Meaning of colors, lines is the same as in Fig. 12 48 The alternative dimensionless stress: , with the unit of time set by the shear rate, t � u =̇� , is more useful for collisional shear flows, see Ref. [118], and is thus not adopted here fore the sake of brevity. with all lengths non-dimensionalized by d ′ p , and B 0 = 3 ⟨d 3 p ⟩ ≈ 0.766 , accounting for the third non-dimensional moment of the size distribution ⟨d 3 p ⟩ ≈ 1.245 (equals unity for monodisperse particles), for more details see Refs. [31,168,169,171,174].
Various energy densities are plotted in Fig. 15. Note that up to the second-last term in Eq. (106), one has exact identities, with w * = w C ∕(B 0 C) = ⟨ 2 ⟩∕2 (small circles and magenta line), whereas the transformation to 2 is one ambiguous choice out of many, which only holds approximately (large circles and blue lines).
How to compute is actually based on the analysis of the incremental (tangent) stress-strain evolution in the elastic regime (details not shown-work in progress). For observing the (smoothed) non-dimensional bulk modulus: ⟨B⟩ = P * ∕ v ≈ const. , very slow isotropic compression simulation data are analyzed at finite, small strain-steps, v ≈ 10 −4 , with P * = P ∕(B 0 Ch 21 ) = ⟨B⟩ , ignoring data below jamming or during plastic, dynamic events with K > 10 −5 , which also assures that P ≈ P , due to negligible kinetic and dissipative stresses. The elastic strain, , based on the re-scaled pressure, P * ∕⟨B⟩ , is plotted in Fig. 16, showing that the linear contact model particle simulation data agree perfectly well with the linear granular model prediction from Sect. 3.2.3, with one adjustable parameter p 0 ≈ 0.0645 (and the correction term h 21 ). Note that the Hertzian model predicts a qualitatively different non-linear behavior of the re-scaled pressure, proportional to 3∕2 . Crucial here is the overlap-correction for stress, h 21 = 1 − ⟨ 2 ⟩∕(ad p ⟨ ⟩) with a = 2 , accounting for the reduced particle center-distances, see Ref. [168], for very soft, deformable particles that do not break/ fracture under strong compression 49 .  Fig. 13, where all parameters are the same, except for T e = 2.10 −4 , which leads to slightly lower steady states, and a much more rapid (exponential) dissipation of energy. Meaning of colors, lines is the same as in Fig. 12  Both w * e and w * C are based on the same solution of , for several largeamplitude, isotropic, loading-unloading cycles, see Sect. 4.1.1, where w * C is almost in agreement with w P sim = w * sim ∕h 2 21 (large circles), see main text 49 From the present simulation data (of first loading), one can identify (not shown) the (almost) constant ratios (dependent on the particle size and force probability distributions), between overlap and elastic strain: squared overlap and elastic strain: as well as overlap and squared overlap: using the elastic strain deduced from the scaled, dimensionless pressure, = * ∕h 21 = P∕(p 0 Ch 21 ) , with p 0 = B 0 ⟨B⟩ ≈ 0.0645 , and the large-overlap abbreviations h 1 = 1 − ⟨ ∕d p ⟩ and h 21 , with a = 1, see main text, and Ref. [168]. The ratios d 1 , d 2 , and d 3 , close to jamming, change slightly for further unloading, which reflects a minute change of the probability density function of overlaps (contact forces)-not detailed here, see Refs. [50,153,154]. The coordination number that is an essential ingredient of the linear granular model, see Sect. 3.2.3, is plotted in Fig. 17, requiring three parameters C 0 , C 1 , and C , as given in the inset. The particle simulation data are in perfect agreement with the analytical form of C( ) , from Eq. (56), during both loading and unloading.
This also allows to deduce the jamming density from the elastic strain, , or from the average contact overlap ⟨ ∕d p ⟩∕d 1 , representing the (virtual) stress-free reference configuration, directly measurable only at unjamming, as plotted in Fig. 18: which eventually allows to compare its evolution with the theoretical predictions for arbitrary strain paths, see also Fig. 5. In Eq. (107), the approximation deteriorates for large overlaps, where there are various similar, not exactly identical choices to (directly) deduce J from elastic strain, Δ , that can be approximated (indirectly) from overlap, .

Conclusion and Outlook
The focus of this paper was on yielding and un-jamming/ jamming of granular matter, a study inspired by the late Bob Behringer, to whom this work is dedicated. In an attempt to combine theoretical considerations with numerical/ experimental observations on granular matter, the authors propose a minimalist macroscopic model to capture qualitatively all states of granular matter, and which even can be (107) J = ∕ exp( ) ≈ exp(−⟨ ∕d p ⟩∕d 1 ) , solved analytically in various special cases. Furthermore, the paper contains a review of literature on GSH as well as on particle simulations, which are compared in relation to each other and eventually used to quantitatively calibrate the GSH theoretical model with existing numerical results from a simple, frictionless, soft particle model.
The system considered was a representative volume element (RVE) of granular matter, homogeneous, i.e., without gradients and with no walls. The granular material was considered in fluid-like and solid-like states, as well as during continuous changes between these states. Particular focus was on the transition from elastically stable to instable, which is a novel contribution since the latter states can be highly dynamic, a situation that is not treatable by, e.g., standard elasto-plastic approaches or critical state theory.  Fig. 18 Jamming density, J = J , with deduced from the pressure, in Fig. 16, from the same data as in Fig. 15, with simulation data for two loading-unloading cycles, and model data for three cycles. The parameters in the constitutive model introduced in Sect. 4.1.1, which lead to this good quantitative match with simulations are: J0 = 0.6567 , b ∞ = max = 0.021 , ∞ = 0.30 , as taken from Ref. [49] (where b ∞ = max ), and slope of the reversible branch, p − v = 0.02 ≈ b ∞ . Unjamming takes place when the jamming density hits the density (green line), and due to T g ≈ 0 , here, there is no evolution of J below jamming Based on the rather complex, but versatile granular solid hydrodynamics (GSH), a much simplified qualitative model that includes un-jammed, gas-or fluid-like states as well as jammed solid-like states (elastically stable) was proposed and studied both analytically as well as numerically. Furthermore, various transitions and intermediate states could be identified and better understood in the framework of this simple GSH type model, which has only four state-variables, density, elastic strain (isotropic and deviatoric) and granular temperature, unifying all the states and transitions of granular matter that we could imagine. In order to keep this universal model attempt transparent, the model equations were often over-simplified by setting most parameters constant, so that the structure of the model equations rather than the consequences of additional constitutive assumptions could be tested. Analytical solutions of the model were possible for cases where either one or more state-variables was fixed or set to zero, while other cases involved either purely isotropic or pure shear modes of deformation, which typically removes considerable fractions from the equations to render them solvable.
The model was generalized to include soft particle phenomenology, and even quantitative rheology, elastic and dissipative responses, as inspired by recent soft particle simulations. Also a strictly non-thermal limit (removing the granular temperature) was considered, as well as perfectly plastic, elastic or intermediate states-possibly related to the critical state and the elastic instability, which was actually the main focus and reason to start this research. A major open question about the size of the volume in which such instabilities occur, cannot be answered in this study, since we assumed homogeneity inside the RVE on the continuum theory level.
Even though rather simple, the minimal universal model is capable of following the granular system from very low (dilute granular gas) to very high densities (dense jammed granular solid), including various transitions and all the transients. In order to limit complexity, the model was considered for a homogeneous (gradient-free) system that could be either seen as a RVE, or as material point of a full continuum model. However, it is not clear which size this material element should have. From particle simulations with a few 1000 particles, it is clear (data not shown here) that the system is never really homogeneous, and that zones of plastic deformations can range from a few particles up to system spanning events. This inhomogeneity within the RVE was enclosed in the probabilities for plastic deformations that are an extension to the classical GSH.
Considering jamming, we report a very slow, "halfhearted" transition to the jammed state, as observed in both particle simulations and GSH, where the true jamming density J , is established above, not at the minimal possible jamming density J0 < J , even in the absence of perturbations due to granular temperature, T g , just due to the occurrence of plastic (irreversible) deformations that lead to better, more efficient packing during compression. A more detailed study of the rate dependence and thus dependence on T g , was beyond the scope of this study.

Modes of un-jamming
Once jammed, the first mode of isotropic un-jamming appears trivial: decompression of the system makes the density decrease and un-jamming takes place when the elastic strain vanishes. However, the density at which un-jamming takes place is not the same as the jamming density, it rather depends on the history of the packing. Perturbations by tapping or over-compression both can result in un-and re-jamming densities considerably larger than the lowest possible, the random loose packing density. The longer/stronger the system is perturbed, the larger the jamming density will be, but the approach to this upper limit is realized very slowly, so slow that it requires very many cycles to be reached. Whether there are well defined random loose and random close packing densities, below/above which the system cannot jam/un-jam anymore is an important open question. Both limit densities are very sensible to the protocol one uses to approach and realize/measure them, especially in the absence of friction as relevant for soft, gel-like particles that resemble many of the simulations referred to in this study.
The second mode of un-jamming is by plastic yielding, which involves irreversible deformations/re-structuring of the solid granular matter, but does not involve dynamics or granular temperature-at least not in the classical picture. Plastic events occur with a certain probability, see Ref. [49], which is larger the closer the system is to un-jamming or the larger the elastic shear strain (and/or stress) is, which was previously accumulated. This mode involves the more classical world of elasto-plastic continuum mechanics and rheology for example see Refs. [4,30,175]. The evident lack of a dynamic state-variable is at the origin of many difficulties with those elasto-plastic concepts, in particular when the deformation rates become larger and larger. Modern concepts like fluidity or non-local models have been proposed during the last years to overcome this problem [18,70,99,116,175,176], however, the proper account for the granular temperature in the elastic regime, and for unjamming, is still an open issue that is at least partly solved now.
The third mode of un-jamming is a transition occuring via an elastic instability, i.e., the loss of convexity, and then involves deformations of the solid granular matter that can occur without penalty (work), at the onset of concavity (elastic instability) or, are even activated/pushed by the external stresses (in the concave regime, or closeby). This mode is seemingly different from plastic yielding, since it allows for dynamics (granular temperature) to build up, grow, and eventually push back the system into a mechanically stable elastic state before/while it is dissipated.
How much different-if at all-plastic and elastic yielding really are has to be seen, and is subject of ongoing research.

Outlook and open questions
Besides extending the theory to general, inhomogeneous systems with gradients, further research is also needed to: • Generalize the present version to arbitrary tensorial form in three dimensions, involving also the so far lacking third invariant of the tensors. • Connect quantitatively the granular gas and fluid (standard) kinetic theory (SKT) with GSH. • Determine the proper shape of the energy density for granular solids from particle simulations, and for more realistic materials with friction, cohesion, polydispersity, etc. • Sort out if the present version of GSH without a microstructural (fabric) tensor is sufficient or needs to be improved-is the elastic strain enough? • Find the different mechanisms of relaxation, creation and destruction of energy in the elastic strain degrees of freedom as well as the dynamic, kinetic, granular ones. • Identify the relaxation/evolution dynamics and the interplay of the multiple mechanisms, represented by the various different terms, of the state-variables below, above and during un-jamming/jamming? • Find out how to bridge the gap between the discrete (in time) local and global plastic events (re-arrangements of the micro-structure) and the macroscopic/continuum picture presented here. • Discussion of the fundamental principles of time-reversal symmetry and other items related to relaxation, entropy production, etc., goes beyond the scope of this paper, but is subject of ongoing discussions/collaborations. Present research is aimed to address the remaining challenging questions: What are the differences and similarities of the driving forces/mechanisms? And, can they indeed all be combined in a single universal model as attempted in this study?