{'fraction_non_alphanumeric': 0.0442613948380011, 'fraction_numerical': 0.009976203551162365, 'mean_word_length': 4.072987278298821, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Similarities between GSH, Hypoplasticity and KCR

Accounting for elasto-plastic motion in granular media, hypoplasticity is a state-of-the-art constitutive model derived from data accumulated over many decades. In contrast, GSH, a hydrodynamic theory, is derived from general principles of physics, with comparatively few inputs from experiments, yet sporting an applicability ranging from static stress distribution via elasto-plastic motion to fast dense flow, including non-uniform ones such as a shear band. Comparing both theories, we find great similarities for uniform, slow, elasto-plastic motion. We also find that proportional paths and the Goldscheider rule used to construct barodesy, another, more recent constitutive model, are natural results of GSH's equations. This is useful as it gives these constitutive relations a solid foundation in physics, and in reverse, GSH a robust connection to reality. The same symbiotic relation exists between GSH and KCR, or Kamrin's non-local constitutive relation, a model that was successfully employed to account for a wide shear band in split bottom cells.


I. INTRODUCTION
Being a subject of practical importance, elasto-plastic deformation of dense granular media has been under the focus of engineering research for many decades if not a century [1][2][3][4][5][6].
The state of geotechnical theories, however, remains confusing for outsiders: Innumerable models compete, employing strikingly different expressions. In a recent book on soil mechanics, Gudehus employed metaphors such as morass of equations and jungle of data to describe the situation [6].
Granular dynamics is frequently modeled employing the strategy of rational mechanics, by postulating a function C ij , such that the constitutive relation, holds, whereby C ij is a function of the Cauchy stress σ ij , strain rate v kℓ (the symmetric part of the velocity gradient) and density ρ. (More generally, ∂ t ≡ ∂ ∂t is to be replaced by an objective derivative, such as Jauman -here and below.) Together with mass and momentum conservation, it forms a closed set of equations for σ ij , ρ and the velocity v i . Hypoplasticity [3,4] is such a theory, same as barodesy [7].
Both use a single constitutive equation, without the recourse to yield or potential surfaces.
Although elasto-plastic theories are somewhat different, Einav [8] showed that the latter can be seen as a special case of the former. All these models yield realistic accounts of the complex elasto-plastic motion, providing us with vast amount of highly condensed and intelligently organized empirical data. We concentrate on comparing GSH to hypoplasticity here, because elasto-plastic theories, with its case separations, are not easy to deal with analytically, and because a comparison to barodesy already exists [9].
Typically, the hypoplastic model starts from the rate-independent constitutive relation, with H ijkℓ , Λ ij and ǫ functions of stress and density. Though continually improving functional dependence is being proposed for them, some known and systematic shortcomings remain: 1. The applicability of hypoplasticity is confined to uniform systems. Lacking a characteristic length scale, it does not account for shear band or evanescent creep flow (ie. the exponentially decaying penetration of flows into the static region). There have been two main approaches to ameliorate this: introduction of gradient terms [10], or addition of variables to account for the Crosserat rotation and couple stress [11]. The latter works well, but leads to a rather more complex theory. And it provokes questions about the underlying physics: If couple stress is indeed important in the shear band, because it is fluid, why then is it not important in the uniformly fluid and gaseous state of granular media -or, coming to think of it, in nematic liquid crystals [12]?
2. In hypoplasticity, rate-independence is an input. Yet it is not a robust feature of granular behavior. Fast dense flow is rate-dependent, as is the critical stress in the presence of external perturbations (such as given by tapping). Rate-independence needs to be explained, understood, and the extent to where it is valid clarified.
3. Elastic waves exist in granular media. Yet if one tries to account for them employing Eq.
(2), one finds falsely that they should always be over-damped. More generally, granular media appear softer and more dissipative during elasto-plastic motion than in many other circumstances -foremost elastic waves, but also incremental stressstrain relations and ratcheting. Focusing on the former, hypoplasticity struggles to account for any of the latter. (See however works on intergranular strains [13].) 4. Fast dense flow and collisional flow are not accounted for. 5. Finally, hypoplasticity lacks energetics. One does not know whether a given state is stable because it has an energy minimum, or instable because the energy turns concave. Nor do we know how much energy is being dissipated say in the critical state.
Moreover, H ijkℓ has 36 elements, Λ ij six, all functions of a tensor and a scalar. These are a lot of functions to choose. So it would be helpful if the second derivative of the energy were related to H ijkℓ and/or Λ ij (contradicting in fact the name hypoplasticity).
Granular solid hydrodynamics (GSH) is derived from basic principles of physics, including conservation laws and the second law of thermodynamics, by employing the hydrodynamic procedure [14,15]. The result is a continuum-mechanical theory [16][17][18][19][20][21][22] that leaves relatively little leeway, with a handful of scalar parameters to fit experiments. The most important one is the energy w, chosen such, first of all, that elastic waves, incremental stress-strain relation and static stress distributions are well accounted for. Second, w contains the information that no stable static stress distribution exists if the density ρ is too small for enduring contacts, or if the shear stress is too large for the given pressure. Third, by reducing GSH to Eq.(2) for slow, uniform shear rates, w's second derivative is found to yield H ijkℓ and Λ ij .
The thermodynamic consistency of hypoplasticity has been demonstrated independently by Fang et al, see [23,24]. This was done generally, without specifying an expression for w.
The energy w is taken as a function of the density ρ, granular temperature T g and elastic strain u ij . The state variable u ij accounts for the slowly-varying, large-scaled elastic deformations of the grains, with the elastic stress given by the hyper-elastic relation, The state variable T g quantifies the fast fluctuating elastic and kinetic energy of the grains.
It is frequently proportional to the shear rate, T g ∝ v s (≡γ), and a useful indicator for the system's behavior, from the quasi-static via the dense flow to the collisional one.
For collisional and fast dense flow, with a shear rate of 10 3 /s and more, there is vigorous jiggling and aggitation, so T g is obviously important. For elasto-plastic motion (usually referred to as "quasi-static" but see below), with shear rates around 10 −3 /s, T g may seem irrelevant: Aside from an occasional slip, grains only participate in the motion dictated by the macroscopic shear rate, with no perceptible deviations that may contribute to T g .
Nevertheless, such slips lead to vibrations in the region that imply a T g still many orders of magnitude larger than the true temperature of the grains. More importantly, as we shall see, these vibrations are what lead to plasticity.
If T g ∝ v s is yet smaller, below 10 −5 /s, GSH predicts a rate-dependent transition, to irrelevance of T g and diminishing plasticity. This is supported by recent experiments [25].
In what follows, we shall briefly present GSH in Sec.II, first its three coupled relaxation equations and an expression for the elastic energy, followed by a discussion of how best to understand T g and u ij , and why a second T g for the fluctuating elastic energy alone is counter-productive. In Sec.III, we obtain the critical state as the steady-state solution of the three relaxation equations. The approach to the state, for both high and low initial densities, is given by a simple analytic solution. In addition, GSH accounts for the variously observed fact that disturbing the critical state with some tapping, its stress becomes highly rate-dependent. We then compare GSH in detail to hypoplasticity and neo-hypoplasticity, In Sec.V, we discuss the transition to the quasi-elastic regime for extremely slow rates.

A. The state variables and their dynamics
A complete set of state variables is one that uniquely determines a macroscopic state, including its energy and stress. GSH's state variables are the density ρ, velocity v i , the granular temperature T g and the elastic strain u ij . The latter two quantify, respectively, the temporally fluctuating elastic and kinetic energy of the grains and their slowly-varying large-scaled deformations. T g is increased by viscous heating (f v s ) 2 , decreased by relaxation, and its balance equation is with R T T g the relaxation rate, and R T ξ 2 ) T a is the "ambient temperature" maintained by tapping or a sound field. For a steady, uniform v s , we have implying T g = f v s if there is no external perturbation, and T g = T a if the shear rate is zero.
As is easy to see from Eqs. (6,8) below, it is T g = f v s that leads to rate-independence. The easiest way to destroy it is therefore to introduce some T a , see Sec.III A 4.
The importance of T g in the collisional state is easy to accept. Recent progresses show it exists equally in the turbulent flows of rapid dry granular avalanches [26]. That T g is relevant even in elasto-plastic motion may surprise, yet as already discussed in the introduction, this is a result of T g , many orders of magnitude larger than the true temperature, being the physics underlying plasticity, see the next few paragraphs.
The elastic strain u ij alone suffices to account for elasticity and plasticity. To understand this, one needs two steps: First, granular stresses relax with a variable relaxation time.
Second, this relaxation is closely related to plasticity.
The free surface of a granular system at rest is frequently tilted. When perturbed, when the grains jiggle and T g = 0, the tilted surface will decay and become horizontal. The stronger the grains jiggle and slide, the faster the decay is. We take this as indicative of a system that is elastic for T g = 0, transiently elastic for T g = 0, with a stress relaxation rate ∝ T g . A relaxing stress is typical of any viscous-elastic or transiently elastic systems, such as polymers [27][28][29][30]. In granular media, the relaxation rate is not a material constant, but a function of the state variable T g -a behavior that we name variable transient elasticity.
For given ρ, the stress π ij ≡ −∂w/∂u ij is a monotonic function of u ij (as long as the energy w is stable and convex, ∂ 2 w/∂u ij ∂u kℓ > 0). Therefore, there is no principle difference in considering either relaxation. Yet since strain is a geometric quantity, stress a physical one that includes material parameters such as the stiffness (stress dependent and highly anisotropic in granular media), considering strain relaxation is simpler.
We take the equation of motion for u ij such to reflect variable transient elasticity, finding that it suffices to capture elasto-plasticity. Starting from pure elasticity, ∂ t u ij − v ij = 0, adding a relaxation term, ∂ t u ij − v ij = −u ij /τ , we allow for two tensorial coefficients, λ ijkℓ and α ijkℓ , to parameterize the efficacy of relaxation and deformation rate, If stress anisotropy changes appreciably, convective terms (or the so-called objective time derivative) become important, then one needs the substitution (see the derivation in [27]) Since ∂ t π mn = (∂π mn /∂u ij )∂ t u ij = −(∂ 2 w/∂u mn ∂u ij )∂ t u ij , and T g = f v s , see Eq.
Given the equality of structure, GSH is as capable of accounting for elastoplasticity as hypoplasticity, if an appropriate energy expression can be found. (Allowing the density to vary yields slightly more complicated expressions.) Next, dividing u ij into its trace /v s , and specifying the two tensors with two elements each, α, α 1 , λ, λ 1 , we arrive at which may equivalently be rewritten as two relaxation equations, Eqs. (4,11,12) are the three relaxation equations of GSH. They contain three inverse time scales: R T T g , λT g , λ 1 T g , a length scale ξ T , and the parameters f, α, α 1 . In dense media, R T T g is large, of order 10 3 /s or more. With T g = f v s , the scales λT g , λ 1 T g are, for the shear rates typical of soil-mechanical experiments, orders of magnitude smaller, around 1/s for v s = 10 −2 /s. The length scale ξ T is a few to a few tens granular diameters in dense media.
The steady state solution: T g = f v s , u s = u c and ∆ = ∆ c , implying constant elastic strain and stress at given shear rate, is clearly the critical state, see Sec. III A below.
Conservation of momentum and mass, close the set of equations, where the Cauchy stress σ ij is see [17] for derivation. The (off-diagonal Onsager) coefficients α, α 1 , a necessary consequence of the same coefficients in Eqs. (9,10), soften and mix the stress components. The term preceded by α 1 is smaller by one order in the small elastic strain u ij and may usually be neglected (though not in Sec. III B 1). α, α 1 are functions of the density for elasto-plastic motion, though they vanish identically for static stresses, see Eq.(77). The pressure P T ∝ T 2 g is exerted by jiggling grains, and the last term is the viscous stress. Both are ∝ v 2 s for T g ∝ v s and relevant only for fast collisional and dense flows, giving rise to Bagnold scaling there.

B. The energy
An explicit expression for the stress via Eqs. (3,14), is delivered by the energy w that we take to be w = w ∆ (∆, u s , ρ) + w T (T g , ρ). The part relevant for elasto-plastic flow is w ∆ .
(w T , important for collisional and dense flow, is neglected here.) Denoting with u * ij , π * ij the respective traceless tensors, we take w ∆ as Given by the energy of linear elasticity multiplied by √ ∆, this form is clearly inspired by the Hertzian contact, though its connection to granular elasticity goes far beyond that. It includes both stress-induced anisotropy and the convexity transition (see below). In addition, the expression for the elastic stress π ij has been validated for: • Static stress distribution in three classic geometries: silo, sand pile, point load on a granular sheet, calculated employing ∇ i π ij = ρg i , see [31][32][33].
The GSH framework is general enough that one may employ any energy expression one deems appropriate. In fact, the above expression has been chosen for simplicity and manipulative ease, with a more realistic one given in [39]. In Sec. III B, we shall show and compare a number of different elastic energy expressions for realistic modeling.
Stress-induced anisotropy: In linear elasticity, w ∝ u 2 s , the velocity of elastic waves ∝ ∂ 2 w/∂u 2 s does not depend on u s , or equivalently, the stress. For any exponent other than two, the velocity depends on the stress, and is anisotropic if the stress is.

Convexity Transition:
There are values of stress and density in granular media for which elastic solutions do not exist: First, for densities smaller than the random loose packing density, ρ < ρ ℓp , at which the grains start to loose contact with one another, no static elastic states can be maintained. Second, static elastic states collapse when the shear stress is too large for the given pressure. On a tilted plane, the angle at which this happens is frequently referred to as the angle of stability, ϕ st . The actually measured angle varies with the spacial dependence of the stress components and does not have a unique value. But in an infinite geometry, for uniform stresses, elastic solutions fail at a well-defined value.
Note that ϕ st , accounting for collapses of static stresses at T g = 0, is unrelated to the critical state, which exists only in the presence of a finite granular agitation, T g = f v s = 0.
In Sec. IV C, we shall identify the angle of repose ϕ re with the critical angle. Again, the measured angle may vary because the stress is not uniform.
In the space spanned by the density and stress components, the surface where the second derivative of the elastic energy changes its sign (turning from convex to concave) is a divide: On one side elastic solutions are stable, with grains at rest, on the other elastic solutions are instable, infinitesimal perturbations suffice to destroy any solution, and the grains are always agitated. The elastic energy of Eq (16) is convex only for It turns concave when this condition is violated. (The second condition may be derived by On a plane inclined by the angle ϕ, with y the depth of the granular layer on the plane and x directed along the slope, we take the stress to be π xx = π yy = π zz = P ∆ , π xy = π s / √ 2, π yz , π xz = 0. Integrating ∇ j π ij = g i ρ assuming a variation only along y, we find π xy = g sin ϕ ρ(y)dy and π yy = π xy / tan ϕ. The angle of stability ϕ st (for infinite geometry) Modifications from the proximity to walls or floor are not included here.
For the parameterization of the coefficients A, B, we need a second characteristic density, the random-close one, ρ cp . It is the highest density at which grains may remain uncompressed. Assuming B/A is density-independent (typically 5/3), and denotinḡ with B 0 > 0 a constant. This expression accounts for three granular characteristics: • The elastic energy is concave for ρ < ρ ℓp , rendering all elastic solutions instable.
• The energy is convex between ρ ℓp and ρ cp , ensuring the stability of elastic solutions in this region. In addition, the density dependence of sound velocities (as measured by Harding and Richart [40]) is well rendered by an expression ∝ B(ρ). (The usual formula used to report the resuslts of [40] implies an unphysically concave energy.) • The elastic energy diverges, slowly, for ρ → ρ cp , approximating the observation that the system becomes quite a bit stiffer there.

C. Two Elastic Limits
Because there is no space for granular rearrangements at the closed-packed density ρ cp , plastic motions cannot take place and the system should be fully elastic. Taking α, α 1 , λ, λ 1 → 0 for ρ → ρ cp in Eqs. (9,10,14), reduces GSH to simple elasticity, Close to ρ cp , the dependence on ρ cp − ρ is more sensitive than any on ρ directly. Hence we take, with a 1 , a 2 , a 3 , λ 0 , λ 10 > 0 and r ≡ 1 − ρ/ρ cp : Granular media should also be fully elastic for T g → 0, see the discussion at the end of the introduction, and in Sec.V. In this limit, the stress increases steeply with the shear rate, until the yield stress is reached, at which point the system collapses. There is no stress relaxation, hence neither a critical state.
D. More on T g and u ij As T g and u ij are novel variables that entail conceptual subtleties, a discussion is useful for those who want to look beyond the present equations. It is also useful because Sun Qichen and coworkers [41] have variously proposed a "Twin-Granular-Temperature Theory" that we find necessary to comment on.
1. The granular temperature T g In a uniform medium (such as water or copper) there are two length scales, macroand microscopic. All degrees of freedom may be divided into one of the two groups. A hydrodynamic theory takes the degrees of freedom from the first group as explicit variables, including mass, momentum and energy density, also the strain field in the case of solids. All degrees of freedom from the second group are taken summarily, with their contribution to the energy lumped together as heat, and characterized by the temperature T .
In granular media, there is an intermediate group, consisting of momentum and deformation of each grain. In the kinetic theory or DEM, these are taken as independent variables.
But a summary inclusion again suffices here, with the associated energy lumped into granular heat, characterized by T g . As the kinetic energy changes rapidly into an elastic one, and back, during a collision, both types of energy interact, equilibrate, and must be lumped into one single macroscopic variable T g . On the other hand, distinguishing T g from T is necessary: T g is many orders of magnitude higher than T in elasto-plastic motion, and when this is the case, granular stress relaxes, giving rise, as discussed, to plasticity.
There are also grain-sized deformations, such as in force chains, that are static or slowly varying in time. Their contribution is not included in T g , but in that of u ij , see below.

The elastic strain u ij
Accounting for the macroscopic elastic stress that varies slowly in space, the elastic strain is also slowly varying. As energy is positive and conserved, the macroscopic energy density w is the sum of all (slowly varying) mesoscopic contributions, w mes . In other words, w is the coarse-grained mesoscopic one, w = w mes , and since the latter includes the contribution from the force chains, the former also does.
It is important to realize that the elastic strain cannot be obtained by coarse-graining, . This is because both the energy and the elastic stress are coarse-grained quan- 3. Rebuttal of "Twin-Granular-Temperature Theory" Sun and coworkers have written a number of papers [41] proposing a generalization of GSH by introducing two granular temperatures, T k and T c , with T k taking account of the kinetic, and T c the elastic energy (including contributions from force chains). It should be clear from the last two sections that neither is a well-behaved macroscopic variable. The temporally fast varying elastic energy from collisional deformation is necessarily part of T g , while the slowly varying contribution from force chains are part of u ij .

E. Summary
With the differential equations given and the energy density specified, GSH is a fairly well-defined theory. It contains clear ramifications and provides little leeway for retrospective adaptation to observations. As we shall see in the rest of the paper, a wide range of granular phenomena is well encoded in these equations.
The three relaxation equations, (4,11,12), possess a steady state solution: The physics of the critical state is: Deformed grains give rise to an elastic stress. When grains are sheared past one another, they rattle and jiggle. And when they do, the deformation and stress is slowly lost, because of the loosen contact to their neighbors. As a consequence, a shear rate not only increases the deformation, as in any elastic medium, but also decreases it. A steady state exists in which both processes balance, such that the elastic deformation and stress remain constant over time. This is the critical state. Moreover, the critical stress is rate-independent, because (1) the increase in deformation is ∝ v s , (2) the relaxation is ∝ T g , and (3) T g ∝ v s in T g 's steady state.

A Stationary Elastic Solution
Taking ∂ t T g , ∂ t u * ij , ∂ t ∆ = 0 in Eqs. (4,11,12) for v * ij = const, v ℓℓ = 0, we have With ∆ = ∆ c , u * ij = u cvij and ρ being constant, so is the associated elastic stress π ij (u kℓ , ρ). The stress value is obtained by inserting Eqs. (25,26) into Eqs. (17,18,19), The loci of the critical states thus calculated [9] (though employing the slightly more general energy of [39]) greatly resembles those calculated using either hypoplasticity or barodesy [7] Same as the Coulomb yield of Eq (21), the critical ratio σ c /P c = π c /P ∆ is also frequently associated with a friction angle. Since the former is relevant for vanishing T g , while the latter requires an moderately elevated T g ∝ v s , it is appropriate to identify one as the static friction angle, or yield angle, and the other as the dynamic one, or the critical angle. In Sec. IV C, we shall present reasons why the critical angle should be identified with the angle of repose ϕ re , the angle at which granular flows will come to a halt for uniform stress.
Textbooks on soil mechanics frequently state that the friction angle is independent of the density -although they do not usually distinguish between the dynamic and the static one.
We shall assume, for lack of better information, that both are.

Approach to the Critical State at Constant Density
Solving Eqs. (4,11,12) for u s (t), ∆(t), at constant ρ, v s , with the initial conditions: ∆ = ∆ 0 , u s = 0 (same as P = P 0 , σ s = 0), the relaxation into the critical state is given as (drained monotonic triaxial) experiment, and II to the simulation by Thornton and Antony [44], both in the plots as originally given.
This is an exponential decay for u s , and a sum of two decays for ∆. It is useful, and quite demystifying, that a simple, analytical solution exists.
Typically, we take λ/λ 1 ≈ 3, as the decay of u s and f 1 are faster than that of f 2 . The prefactor f 2 may be negative, and ∆(t) is then non-monotonic. The associated pressure and shear stress are those of Eqs (27,28,29). For a negative f 2 , neither the pressure nor the shear stress is monotonic. For the system to complete the approach to the critical state, the yield surface of Eq.(20) must not be breached during the non-monotonic path.

Approach to the Critical State at Constant Pressure
Frequently, the critical state is approached at constant pressure P , or a stress eigenvalue see Eqs. (11,12), implying a rate-dependence via Eq.(5), If there is no tapping, T a = 0, we retrieve the unperturbed values,ū c = u c ,∆ c = ∆ c , σ c = σ c . With tapping,ū c ,∆ c ,σ c decrease for increasing T a , and increase with increasing shear rate f v s ≡ f |v s |. see Fig 2. (Note we have only considered the critical state at given shear rate, not the approach to it. So the result holds both at given density and pressure.) The above consideration is the basic physics of the observation of [45]. It helps to put rate-independence, frequently deemed a fundamental property of granular media, into the proper context. A quantitative comparison is unfortunately made difficult by the highly nonuniform experimental geometry. Nevertheless, some comparison, even if unabashedly qualitative, may still be useful. In [45], the torque τ on the disk on top of a split-bottom shear cell is observed as a function of its rotation velocity Ω and the shaking acceleration Γ. Now, (τ, σ s ), (Ω, v s ), (Γ, T a ) are clearly related pairs. Assuming the lowest order terms suffice in an expansion, we take σ c ∝ τ , Γ ∝ T a . If v s were uniform, Ω ∝ v s would also hold.  [22] for details]. Inset is the experimental curve of [45], with the torque τ denoted as T , as in [45]. (The stress dip at large Ω, neglected here, is explained in [46].) Since it is not, Ω ∝ v n s , n > 1 seems plausible, because with additional degrees of freedom such as position and width of the shear band, the system has (for given Ω) more possibilities to decrease its strain rate v s . We take Ω = c 2 v 3 s with c 2 = 1rs 2 [replacing T a /v s with Γ/ 3 √ Ω in Eq.(31)] to fit Fig 2, emphasizing that qualitative agreement exists irrespective of the fit.

Load and Unload
Inserting T g = f v s , v ℓℓ = 0 into Eqs (11,12), keeping the shear rate along one direction and allowing reversal, η ≡v ijûij = ±1, we have and conclude that the reason for the difference in apparent stiffness: ∂ t u s /v s and ∂ t ∆/v s , between load and unload, is that η turns negative when the shear rate v * ij is reversed. The strain then relax towards the new stationary values, with a simple relaxation dynamics not at all history-dependent. That loading (η = 1) and unloading (η = −1) have different slopes is frequently referred to as incremental nonlinearity in soil mechanics. It is the reason why Barodesy, a recent model again proposed by Kolymbas [7], has a rather more complicated rate dependence than the one shared by hypoplasticity and GSH. But our comparison of GSH and barodesy did yield essentially quantitative agreement [9,47]. The idea of barodesy is to have a more modular and better organized structure than hypoplasticity, with different parts in the constitutive relation taking care of specific aspects of granular deformation, especially that of proportional paths as summed up by the Goldscheider rule. In Sec. III B 5, we shall focus on the latter showing how the Goldscheider rule and proportional paths naturally arise from GSH.

A Comparison of the Hypoplastic Coefficients
Writing the hypoplastic model as we have the dimensionless coefficients H 1−7 that are functions of ρ and the 3 stress invariants: In [48], the authors present the parameter values for more than forty different granular media, with characterized by three functions F, f 0 , f d and one constant a, given as e 0d = 0.6, e 0c = 0.98, e 0i = 1.15.
Note σ ij is here negative from that in [48]; note also that the stress invariant σ t in cos 3θ and F will be approximated below with cos 3θ = 1, ie. σ 3 t = σ 3 s / √ 6, valid for axial compressions.
Returning to GSH, see Eqs. (11,12) with T g = f v s , we have four dynamic coefficients, and two static coefficient B and B/A = 5/3.
Defining an operator that is valid for any function f (∆, u s , ρ) such as the stress σ ij , in the form of Eq.(34) to obtain We first parameterize the static coefficients. Using the expression for the energy, Eqs. GSH, especially its energy, has been chosen stressing manipulative ease and simplicity of expressions. GSH's robust and qualitative result is the fact that hypoplastic coefficients may be obtained from the second derivative of a material-dependent scalar potential.
Third, starting with the Gibbs potential and choosing n = 0.6, α = 0.1, Niemunis, Grandas Tavera and Wichtmann most recently use its second derivative to obtain the coefficient H ijkl of Eq.(2) -a model that they call neo-hypoplasticity [49]. This is clearly a step that we highly welcome. Calibrating the coefficients using incremental stress-strain relations, they found good agreement with observation. The next three sections are about comparing four more recent hyperelastic models, to help optimize the choice of the energy.

Elastic Energy versus the Gibbs Potential
Elatic constitutive relations that possess an explicit elastic energy or Gibbs potential is usually termed a hyperelastic theory. Engineers tend to look for appropriate Gibbs potentials, because they prefer the stress as variable. On the other hand, only the energy is conserved and useful for the hydrodynamic procedure. Although one can obtain one from the other, simple energy expressions typically possess complicated Gibbs potentials, and vice versa. In this section, we collect some general results connecting both, starting with the basic one of the Legendre transformation: Assuming the elastic energy is a homogeneous function of degree m, we have with B a constant. This leads immediately to (where prime ′ means derivative). The Legendre transformed Gibbs potential is then with u s , ∆ taken as functions of P ∆ , π s , especially ∆ = (BP ∆ F ′ / π s ) A remarkble property of the above energy expression is the simple relation, As a result, the Gibbs potential also has the factorized form: where the dimensionless shearing factor F = F ( u s ( π s )) is For the Hertz contact m = 5/2, we have w = B −2/3 P 5/3 ∆ F .

Comparing Four Hyperelastic Models
We compare four hyperelastic models, with ξ = 2B/5A a density independent constant. GSH's elastic energy (with a unspecified) was proposed and considered in 2003, see [51]. In our later works, we concentrated on a = 1/2, the Hertz contact value, to be definite, and because Hertz contact seems the appropriate granular picture. Three works on hyperelastic models appeared since, all of the form Eq.(45), by Houlsby, Amorosi and Rojas [52], Einav and Puzrin [53], and, as already mentioned, as neo-hypoplasticity, by Niemunis et al [49]. We first note that the power of the "spheric part" P

Houlsby-Amorosi-Rojas
The shearing factor F for these models are therefore, see Eq. (52), with F 0 ≡ F ( π s = 0) and The shearing factor of GSH explicitly shows the elastic instability that it contains: The dimensionless shear stress is restricted to 0 ≤ π s ≤ a (2 + a) ξ. The root becomes imaginary at the right edge and instability sets in. The other three models remain stable for all π s > 0.
To account for the instability of a strictly static slope that is too steep, since there is no dynamics, one has only the energy to work with. A convex-concave transition -very much in analogy to the van de Waals theory of the liquid-gas phase transition -is an appropriate candidate. Therefore, we prefer to use the elastic energy that contains these instabilities.
One may alternatively use brute force and put in the end of elasticity by hand. But then an abrupt transition is forced onto the system, and one cannot account for any precursors.
Circumstances are similar if one aims to account for the lack of elastic solutions at low densities, ρ < ρ lp , relevant for the transition from elasto-plastic motion to fast flows. occurs too early, as compared to the angle of stability (the angle up to which a static layer of sand will remain at rest). And he proposed adding a term given by the third invariant of the elastic strain, u 3 t = u ij u jk u ki , to increase the predicted angle of stability [50], Unfortunately, considering only negative values of χ, he came to the erroneous conclusion that this term does not yield any improvement. We subsequently found that it is possible to yield a realistic angle of stability of approximately 30 • with a positive χ. What is more, incremental stress-strain relations and granular acoustics also support the same χ [39].

Proportional Paths and Barodesy
Barodesy's starting point are the proportional paths as summed up by the Goldscheider rule (GR). Taking pεp and pσp to denote proportional strain and stress paths, they are • A pεp starting from zero shear stress σ s = 0 is associated with a pσp.
• The same pεp starting from an arbitrary σ ij = 0 leads asymptotically to the same pσp obtained when starting at σ s = 0.
Any constant strain rate v ij is a pεp: In the principal strain axes (ε 1 , ε 2 , ε 3 ), a constant v ij means the system moves with a constant rate along its direction, with ε 1 /ε 2 = v 1 /v 2 , ε 2 /ε 3 = v 2 /v 3 , independent of time. gr states there is an associated stress path that is also a straight line in the principal stress space, that there are pairs of strain and stress paths. And if the initial stress value is not on the right line, it will converge onto it.
To understand GR using GSH, we need the stationary solution for an arbitrary pεp, ie.
we need to generalize the stationary solution as given by Eq. (25,26) to include v ℓℓ = 0. Using the superscript ni to mean non-isochoric, we have where the approximate sign holds because v ℓℓ ≪ v s . We also note that the deviatoric stress, elastic strain and strain rate share the same direction,σ ij =π ij = −û ij = −v ij , see Eqs. (14,19,26). If the strain path is isochoric, v ℓℓ = 0, these are simply the critical state, with both the deviatoric elastic strain and stress being constant. They are hence dots that remain stationary in their respective space, given by u * ij | c = u c (ρ)v ij and σ * ij = −σ s (ρ)v ij . If however v ℓℓ = 0, with the density ρ[t] changing accordingly, u c (ρ) and σ s (ρ) will also change, but not the directionv ij , making the dots walking down a straight line alongv ij .
This would indeed be the first rule, except that GR states that it is the total (and not the deviatoric) stress that possesses a pσp. This is a slightly more involved point. With we again have a density dependent prefactor, [1 − α(ρ)]P ∆ (ρ), multiplied by This factor is a constant direction, giving rise to a pσp, if we assume that the system is already in the critical state -as discussed in Sec.III A 1), π s /P ∆ does not depending on the density. Otherwise, we need to rely on π s /P ∆ being only weakly density dependent.
The second rule is much easier to understand. Given an initial elastic strain deviating from that prescribed by Eq (58), u 0 ij = u cvij , ∆ 0 = ∆ c , Eqs (11,12) clearly state that any deviating component will relax toward zero, implying the elastic strain and the associated stress will converge onto the prescribed line.

IV. EXPERIMENTS AT GIVEN SHEAR STRESS
In this section, we examine the ramification of GSH if, instead of the shear rate, the stress or elastic strain is being held constant. As we shall see, taking ∂ t u * ij = 0 in Eq.(11) results in a creep motion of the magnitude f v s /T g = u s /u c , such that an initial T g will co-relax with v s , at the altered rate R T (1 − σ s /σ c ). Moreover, T g ∝ v s relax back to the equilibrium value T g , v s = 0 only for σ s < σ c , growing without bounds for σ s > σ c , though a shear band solution is stable for σ s > σ c . Note hypoplasticity is not applicable here. This behavior has implications for the angle of repose ϕ re (for uniform stresses): As long as the shear stress is held below the critical one, σ s < σ c , the T g ∝ v s -relaxation will run its course, and the system is in a static, mechanically stable state afterwards. For σ s ≥ σ c , however, the system does not come to a standstill. Therefore, defining ϕ re same as in Eq. (21), but employing the critical rather than yield value, we have tan ϕ re = σ c / √ 2 P, with ϕ re < ϕ st .
The inequality holds because the critical state is an elastic solution, while ϕ st of Eq. (21) is the angle at which all elastic solutions become unstable.
Stress-controlled experiments cannot be performed in a conventional triaxial apparatus possessing stiff steel walls, because the correcting rates employed by the feedback loop to keep the stress constant are usually too strong. As a result, too much T g is excited that distorts its relaxation. The situation is then more one of consecutive constant rates, less of constant stress. Instead, one may employ a soft spring to couple the granular system with its driving device, to enable small-amplitude stress corrections without exciting much T g .
Typical examples of stress-controlled experiments are creep (both uniform in space and exponentially decreasing), shear band, also flow on an inclined plane or in a rotating drum.
The flow either comes to a halt, becomes jammed at the angle of repose ϕ re , or starts flowing from that jammed state, becomes fluidized at the angle of stability ϕ st .

A. Uniform Creep
For T g = 0, the system is static, σ s = const, and there is no dynamics, only a stable, static elastic solution. But if T g = 0 initially, both T g and the elastic strain will relax according to Eqs. (4,11,12), and with them also the stress σ s . Maintaining a constant u s , or σ s , therefore requires a compensating shear rate v s , which is observed as creep. As we shall see, the characteristic time of T g is then ∝ (1 − σ s /σ c ) −1 , hence rather long close to the critical stress. So the accumulated shear strain ε s (t) = v s dt can be expected to be large.
In a recent experiment involving a fan submerged in sand, Nguyen et al. [54] pushed the system to a certain shear stress at a given (and fairly fast) rate, producing an elevated T g .
Then, switching to maintaining the shear stress, they observed the accumulation of a large total strain ε s (t) that appears to diverge logarithmically. They called it creeping.
In the experiment, a very soft spring was used to couple the fan and the motor. This fact we believe is essential why the experiment turned out as observed. Usually, triaxial apparatus with stiff steel walls and a feedback loop is used to keep the stress constant.
The correcting rates are strong, and much T g is excited that distorts its relaxation. The situation is hence more that of consecutive constant rates, less of constant stress. A soft spring coupling the granular system with its driving device has much lower correcting rates.
Because of the fan, the stress distribution in the experiment is rather nonuniform, ren-dering a quantitative comparison to GSH difficult. Nevertheless, the experiment's new and structurally robust results are worth closer scrutiny. And we shall employ GSH assuming uniformity to yield a qualitative understanding. Also, we shall assume it is the elastic shear strain u s that is being kept constant, not the shear stress σ s ∝ √ ∆ u s , as both cases will turn out to be rather similar. The relevant equations are still Eqs (4,11,12) (with T a = 0).
At the beginning, as the strain is being ramped up to u s employing the constant rate v 1 , granular temperature acquires the initial value T 0 = f v 1 . Starting at the time t = 0, the elastic strain u s is being held constant. Inserting ∂ t u s = 0 into Eq (11), we have with v s the rate needed to compensate the stress relaxation. Inserting this into Eqs (4,12), we arrive at with the effective T g -relaxation rate reduced from R T to Assuming uniformity, both equations may be solved analytically, if the coefficients are constant. Being functions of the density, they are if the density is. The stress P (t), σ s (t) will then change with time, as will ∆(t). The first equation accounts for the exponential decay of ∆, from both below and above, to the steady state value The relaxation is faster the more elevated T g is. Employing the initial condition T g = T 0 at t = 0, the solution to the second equation is Because of Eq (62), the same solution also holds for the shear rate, v s = v 0 /(1 + r v v 0 t), with v 0 ≡ (u s /u c )T 0 /f , r v ≡ (f u c /u s )r T , implying a logarithmically divergent total shear strain, However, we note that ε s does not actually diverge, because as T g → 0, it enters the quasielastic regime, see Eq.(78) below, where its relaxation becomes exponential. So even if creep is large close to u c , it comes to a halt eventually, and the system is mechanically stable.
Assuming a large T 0 , ∆ is quickly relaxed. Fixing u s is then equal to fixing the shear one may rewrite the factor in r T as The T g -relaxation is slower the closer σ s is to σ c , infinitely so for σ s = σ c . Then u s = u c , ∆ = ∆ c , with T g (t) = T 0 = f v 1 a constant that does not relax. This is indistinguishable from the rate-controlled critical state, which, clearly, may be maintained also at given stress.
For σ s > σ c , the relaxation rate is negative, and T g ∝ v s will grow, or explode, seemingly without bound, see the next two sections, each with a possibility of what happens next.
Only a T g sufficiently large will explode, destabilizing a static shear stress exceeding σ c , an infinitesimal T g will not. This is because the critical stress diverges for T g → 0 [as can be seen calculating the critical stress employing Eq. To summarize, as long as σ s < σ c , both the temperature and the shear flow will relax to zero, with the rate r T ≡ R T [1 − σ s /σ c ]. For σ s = σ c , we have the rate-independent critical state. For σ yield s > σ s > σ c , the system is meta-stable, and easily perturbed into developing a shear band as given in the next section, IV B. For σ s > σ yield s > σ c , the stress necessarily has the full form of Eq. (14), containing both the gaseuous pressure, viscous stress, and also the elastic part that is always in the critical state, see Sec.IV C.

B. Narrow Shear Bands
For π yield s > π s > π c , a stable, localized steady-state solution exists that we may identify as the shear band: Taking ∂ t T g = 0 in Eq.(64), the steady-state solution for 0 ≤ x ≤ ξ sb is The velocity difference across the band is ∆v = v s dx = v 0 s sin(πx/ξ sb )dx, hence This is to be combined with v s ∝ T g ≡ 0 for x ≤ 0 and x ≥ ξ sb . We note that although T g and v s are continuous at 0, ξ sb , neither is differentiable there. Also, the density ρ must be lower in the band, because only than can the shear stress σ s (that has to be uniform) be smaller than the critical stress σ c (ρ) in the quiescent region, but larger in the band.

C. Angle of Stability and Repose
If the system stays uniform at a shear stress exceeding the critical one, σ s > σ c , the variables T g ∝ v s will grow with r T , until the gaseous pressure and viscous stress become relevant, see Eq. (14), while the elastic stress stays in the critical state, This is the stress expression for fast dense flow, again a stable and uniform solution, though no longer rate-independent. We shall not dwell on it here, see [22,42] for results.
On a tilted plane with a ongoing granular flow and a stress as given by Eq.(72), reducing the inclination angle will reduce both the shear stress and the shear rate v s . And the flow will come to a stop, T g ∝ v s → 0, if the stress drops below the critical one, as discussed in Sec. IV A. Hence the angle of repose (for uniform stress) is tan ϕ re = σ c / √ 2 P , see Eq.(61.
Given a static elastic state on a plane inclined by the angle ϕ, the angle of stability ϕ st

D. Wide Shear Bands and KCR
Finally, we consider the coexistence between fast dense flow, T g = f v s , and one of static elastic stress, T g , v s = 0. Some T g will diffuse into the static region, giving rise to an exponentially decreasing creep, as observed by Komatsu et al [55], Crassous et al [56].
We take this "liquid-solid boundary" to be at x = 0, with x > 0 being solid. In the fluid phase x < 0, the shear rate is a constant, as is T g , providing the boundary conditions for the solid part. Pressure P and shear stress σ s are uniform, as is v ℓℓ ≡ 0, but ρ is discontinuous at x = 0, though constant otherwise. (The liquid density needs to be lower than that of the solid one, such that the same stress exceeds the critical one in the liquid, but is below it in the solid.) The variable T g = f v s varies alongx. These circumstances are similar to that of Sec IV A, though variation is in space rather than time. Given stationarity of Eqs (11,12), see Eqs.(62,66). With ∆, u s fixed, so are P, σ s , where especially P = P c if σ s = σ c . Note that the above relations imply ∆/u s = (∆ c /u c )(f v s /T g ). And since 1/µ ≡ P/σ s increases monotonically with ∆/u s , see Eq. (19), the friction µ decreases for increasing f v s /T g .
The balance equation for T g with ∂ t T g = 0, see Eq.(64), yields an evanescent creep, Kamrin et al [57,58] proposed a constitutive relation (KCR) to account for the steady flow in the split-bottom cell [59,60]. A key ingredient of the theory is a new variable fluidity: g ≡ v s /µ, with µ ≡ σ s /P and µ s ≡ σ c /P c . It is taken to obey the equation As g loc = 0 for µ < µ s , this relation is rather similar to Eq.(74), with g assuming the role of T g , and the decay length diverging at the critical stress. For µ ≥ µ s , the system is fluid, and g = g loc essentially constant. Again, T g is constant in the liquid region.
KCR is well capable of accounting for the steady flow in the split-bottom cell. This is fortunate, because again, there is a symbiotic relation relation between GSH and KCR, similar to that between GSH and hypoplasticity. But there are also three differences.
First, KCR assumes σ s > σ c in the liquid and σ s < σ c in the solid, with σ c a constant. This violates momentum conservation, as the stress is continuous across the solid-fluid interface (see also the discussion in Sec.IV B and IV D). Second, we have, for GSH, f v s /T g = 1 in liquid, and f v s /T g = σ s /σ c in solid; yet v s /g = µ ≡ σ s /P on both sides for KCR. Third, Kamrin and Bouchbinder constructed a "two-temperature continuous mechanics" in a recent paper [61], defining a configuration temperature θ c and one for the vibrational degrees of freedom θ v . (Microscopic degrees of freedom are not considered. There would be three temperatures if one did.) They hypothesize that the fluidity g may be related to θ c . We disagree here. As discussed above, g is essentially T g of GSH, which is closer to θ v .

V. THE QUASI-ELASTIC REGIME
Quasi-static is what many in soil mechanics call the rate-independent elasto-plastic motion. The implication is: Being the slowest one possible, elasto-plastic motion stays as is however slow the rate. We believe there is, at very small T g and very low rates -possibly v s ∼ 10 −6 /s, or I ∼ 10 −8 − 10 −9 (where I ∝ v s / √ P is the inertial number of the µ(I) rheology), a rate-dependent transition to another rate-independent regime. If true, only the latter should be called quasi-static, or quasi-elastic. The reasons for our believe are the following six points [see also Sec.II C and the list after Eq.(19)]: 1. We have T g = 0 in static stress distributions, which is well accounted for by the fully elastic Eqs. (23). This should remain so for very slow stress changes.
2. The fact that cyclic load or ratcheting are less dissipative than hypoplasticity predicts [13] may well be interpreted as the onset of a transition to quasi-elasticity. One does not need to invoke intergranular strain here.
3. Quasi-static motion -a consecutive visit of static, equilibrium states -should occur at T g = 0, and be reversible. Yet elasto-plastic motion occurs at elevated T g , is strongly dissipative and irreversible. Reactive and dissipative terms are comparable in size, being exactly equal in the critical state. (This is obvious only within the context of GSH, but one can identify the reactive and dissipative terms in hypoplasticity by comparing it to GSH. Then the same holds true there.) 4. Rate-independence is a hallmark of quasi-static motion. Yet the rate-independence of elasto-plastic motion is (as discussed in Sec.III A 4) easily destroyed. 5. T g is essentially zero in elastic waves, since there is not enough time to excite appreciable amount of it. Waves are well accounted for by the elastic Eqs.(23), but not by hypoplasticity -as they should if there were only the elasto-plastic regime. Employing hypoplasticity, one finds that waves are always over-damped: Starting from momentum conservation, ∂ t (ρv i ) = −∇ j σ ij , or ∂ 2 t (ρv i ) = −∇ j ∂ t σ ij , one inserts Eq.(2) to find the term ∝ H ijkℓ yielding the velocity and the one ∝ Λ ij damping. Both terms are of the same order in q, and given the values of the two tensors appropriate for elasto-plastic motion, they are also comparable in size, implying that elastic waves are always over-damped. (In solids, the damping term is an order in the wave vector q higher than the velocity term, rendering it much smaller for long enough wave lengths. This is the reason elastic waves propagate there.). In fact, elastic waves are a fast yet quasi-static phenomenon, same as incremental stress-strain relations.
6. Happily and perhaps most convincingly, Peng Zhen et al. have most recently observed indications of a rate-dependent transition away from the elasto-plastic regime [25].

A. Strain Increments versus Slow Rates
Faced with the quasi-elastic response of sound waves and ratcheting, a frequent suggestion is to take a small strain increment to be elastic and free of dissipation, but a large one as elasto-plastic and dissipative. Unfortunately, this is not compatible with the notion of quasistatic motion, a consecutive visit of equilibrium, static states: Starting from an arbitrary static state of given stress, applying a small incremental strain that is elastic, the system ends up in another static state, with a slightly higher stress. This second state, as static as the first, is again a valid starting point -any static state is. And the next small increment must again be elastic. Many consecutive small increments yield a large change in strain, and if the small ones are not dissipative, neither can their sum be.
In GSH , it is the strain rate, not the strain amplitude, that decides whether the response is elastic or elasto-plastic. Small strain increments produced with a brief puls of shear rate will produce an elastic response, if T g does not have time to get to a sufficiently high value.

B. How to Observe Quasi-Elasticity
There are at least two ways to observe the elusive quasi-elastic behavior, both by fixing the stress rate and keeping a low T g : Note that a given stress rate may be associated with either of two shear rates, a high, elasto-plastic one at elevated T g and a low, elastic one at vanishing T g . The first method is to slowly incline a plane supporting a layer of grains. In such a situation, the shear rate remains very small, and the system starts flowing only when the yield stress is breached.
A second method is to insert a soft spring between the granular medium and the driving device. If the spring is softer by a large factor a than the granular medium (which is itself rather soft), it will serve as a "stress reservoir," and absorb most of the displacement, leaving the granular medium deforming at a rate smaller by the same factor a. (Employing a feedback loop in a triaxial apparatus with stiff steel walls to maintain a stress rate will not usually work, because the correcting motion typically has strain rates that are so high, that the system slips into an elasto-plastic motion quickly.)
Instead of taking λ, λ 1 → 0 with T g , we choose to generalize T 2 g = f 2 v 2 s to because this was a natural result arising during the derivation of GSH, see [21,22]. For T g ≪ T 1 , we have T g = (f 2 /T 1 ) v 2 s quadratically small in v s , and with it also the relaxation rates ∝ T g . The ratio, in Eqs. (9,10), between the irreversible, plastic term (∝ v 2 s ) and the reversible, elastic term (∝ v s ) diminishes with v s → 0, as has been observed in DEM [35].

VI. CONCLUSION
Taking the stress increment as a function of density, strain rate and the stress itself, the Eq.(1) this paper starts from to set up a constitutive model, is an appropriate and physically sound way to come to terms with the complexity of triaxial results. A rather broader view of granular behavior are captured by introducing a pair of internal state variables, the elastic strain field u ij and the granular temperature T g , with the first accounting for the coarsegrained deformation of the grains, and the latter accounting for their quickly changing elastic and kinetic energy. The theory that does it is called GSH, and the phenomena accounted for include static stress distributions, creep, shear band, angle of repose and stability, also the critical state under tapping, all in addition to the trixial results.