Nonlinear elasticity, yielding, and entropy in amorphous solids.

The holographic duality has proven successful in linking seemingly unrelated problems in physics. Recently, intriguing correspondences between the physics of soft matter and gravity are emerging, including strong similarities between the rheology of amorphous solids, effective field theories for elasticity, and the physics of black holes. However, direct comparisons between theoretical predictions and experimental/simulation observations remain limited. Here, we study the effects of nonlinear elasticity on the mechanical and thermodynamic properties of amorphous materials responding to shear, using effective field and gravitational theories. The predicted correlations among the nonlinear elastic exponent, the yielding strain/stress, and the entropy change due to shear are supported qualitatively by simulations of granular matter models. Our approach opens a path toward understanding the complex mechanical responses of amorphous solids, such as mixed effects of shear softening and shear hardening, and offers the possibility to study the rheology of solid states and black holes in a unified framework.

INTRODUCTION

Amorphous solids including, but not limited to, glasses, granular matter, colloidal suspensions, foams, and polymers yield to external shear (, ). The yielding point, typically characterized by a maximum on the stress-strain curve, is essentially the end point of the solid regime, as the material transitions to plastic flow thereafter. Understanding the nature of yielding in amorphous solids, in the context of statistical mechanics, has been becoming an active research area. Recent theories attempt to explain yielding as a depinning transition (), a first-order nonequilibrium phase transition (, ), a spinodal point of the glass state of equation (, ), or phase transitions in the universality class of a random field Ising model (). However, treating nonlinear effects is intellectually challenging (, ), and consequently, their impact on yielding remains unclear to a large extent.

Nonlinear elastic responses, particularly shear hardening, as evidenced by a rapid increase in the shear modulus before yielding, are not uncommon in amorphous solids. Shear hardening has been widely observed in polymers () and, more recently, in dense hard sphere (HS) colloidal glasses close to jamming (, ). In both cases, hardening is accompanied by entropy vanishing caused by structural constraints: The number of allowed configurations tends to become zero, approaching the maximum stretch limit in polymer chains and the jamming limit in HSs. In this study, we find a previously unidentified type of shear hardening in amorphous solids, where the entropy increases with strain, caused by shear-induced rejuvenation. This nonlinearity also leads to a negative correlation between the yielding strain and the degree of annealing, in sharp contradiction to previous results (, , ).

Let us sketch out our strategy. (i) First, a zero-temperature effective field theory (EFT) for nonlinear elasticity is constructed, on the basis of shift-symmetric Goldstone fields, the phonons, as fundamental building blocks (). The theory also takes advantage of the fact that amorphous solids are typically rotationally invariant (isotropic) and assumes that the solid is homogeneous at low energy, which is true for systems at large spatial scales and/or slow dynamics, such as granular matter under quasi-static shear. The nonlinear elasticity is implemented in the form of an effective potential (cf. strain energy function), from which the stress-strain curve, as well as the onset of instability (yielding), can be calculated.

(ii) To overcome the difficulty in treating dissipative (finite-temperature) effects in the EFT framework, we map the EFT, defined in a (2 + 1)–dimensional flat spacetime with coordinates (x, y, t), onto a gravitational dual model () in asymptotically anti–de Sitter (AdS) spacetime (x, y, t, u), according to the holographic duality (or gauge-gravity duality) (–), where the extra dimension u represents the energy scale (see Fig. 1). The computation of the entropy becomes particularly simple now, since it boils down to the estimation of the black hole entropy given by the famous Bekenstein-Hawking area law (). Such a holographic approach has been already successfully applied to study the dynamics of strongly coupled fluids (e.g., quark-gluon plasma) () and novel strongly correlated materials (e.g., cuprates) ().

Fig. 1.

A schematic representation of the holographic duality.

The gravitational models live in (3 + 1) dimensions, while effective field theories/amorphous solid simulations are in (2 + 1) dimensions.

(iii) To test the predictions from the EFT and the gravitational theory, we perform numerical simulations on soft sphere (SS) models of granular matter. Consistent with the theoretical setup, shear is applied quasi-statically, thermal motions are neglected (for macroscopic grains, the thermal energy is much smaller than the interparticle contact energy), and plastic effects are irrelevant. We consider a common type of shear, planar shear (such as simple or pure shear), where the material is unaffected along the dimensions perpendicular to the plane, and the effective rheology is two-dimensional (2D) (see section S6). The spatial dimensionality is thus reduced to d = 2 in the theories. The simulations are performed in both 3D (see “Simulations of a granular matter model”) and 2D (sections S7 and S8). Below, we discuss our results in detail.

RESULTS

Correlation between yielding and nonlinear elasticity predicted by an EFT

The EFT builds on two scalar fields ϕ(t, x), whose fluctuations δϕ(t, x) represent the displacement fields in the solid (I = 1,2). The stress-strain curves are obtained from a power-law scalar potential V, which plays the same role of the nonlinear strain-energy function in standard elasticity theories (see Materials and Methods and sections S1 and S2 for details). The assumptions used in the theory include (i) rotational invariance (amorphicity), (ii) homogeneity at large scales (larger than the granularity scale of the solid), and (iii) the absence of plastic effects. The validity of these assumptions is tested by a direct comparison between theoretical and simulation results.

The only input to the theory is a power-law form of the stress-strain relation (see Fig. 2A)which corresponds to a concrete choice of the potential V. Here, ν = 2A is an exponent characterizing the nonlinear elasticity, which cannot be determined directly from the effective theory but rather should be considered as a phenomenological parameter. A complete microscopic model, able, for example, to describe the dependence on initial conditions, would be needed to determine directly such an exponent from theory. To continue, an important quantity defined in the theory is the Poisson ratio (the negative ratio of transverse to axial strains)which is related to the unstrained bulk modulus K0 and the shear modulus G0 via in two dimensions (notations with subscript 0 are defined at zero strain). Because there are only two independent parameters A and B (see Materials and Methods), the model is completely fixed by a given combination of ν and .

Fig. 2.

Predictions from the EFT and the gravity theory.

(A) Stress-strain curves extracted from the EFT. Poisson’s ratio is fixed to , nearly the limit of incompressibility. The nonlinear shear exponent v is varied from 1.3 to 3.1 (from purple to blue) insteps of 0.2. The empty circles locate the breaking points. (B) The relation between the nonlinear exponent ν and the breaking point strain εY (red) or the breaking point stress σY (blue) for . (C) The entropy-stress curves for v = 1.00,1.50,3.00,3.75 (from purple to red), from the gravity theory at T = 0.1, with the large stress scaling Eq. 3 indicated (dashed lines). Inset: The theoretical relation .

An interesting prediction is the dependence of breaking point {εY, σY}, defined by a global instability where the speed of sound vanishes, on the nonlinear elasticity. Since the theory is constructed in the hydrodynamic (long-wavelength) limit, this instability is naturally associated with the yielding of the whole system instead of local plastic rearrangements. Analytic formulas for εY and σY are derived and presented in section S2. For a fixed, large (i.e., K0 ≫ G0), corresponding to a nearly incompressible material, ν is negatively correlated with εY and positively correlated with σY (see Fig. 2B and section S2 for more details). Moreover, in the case of a power-law form of the stress-strain curve as in Eq. 1, the EFT predicts a power-law correlation between the breaking point strain εY and the nonlinear exponent ν of the type εY ∼ ν−κ, where κ > 0 depends on the specific details of the EFT potential (see section S2). Thus, nonlinear elasticity increases the strength (maximum stress) of the material while making it more brittle and having an earlier breaking point.

Scaling of entropy predicted by a gravitational theory

In the holographic description (), gravity is coupled to a scalar potential W, which depends on two bulk fields φ (see Materials and Methods and section S3 for details). For a given shear strain ϵ, the stress σ is read off from the black hole geometry using the holographic dictionary, and a power-law scaling as in Eq. 1 is recovered (see section S4). The breaking point can also be obtained by looking at the gradient instabilities of the gravitational modes, and, at least in the decoupling limit, its behavior as a function of the nonlinear exponent ν is identical to the EFT results reported above (). Let us emphasize that the power-law exponent κ and even the power-law functional form are not expected to be universal but rather dependent on the details of the concrete model. Nevertheless, in all cases investigated, εY (or σY) is a monotonically decreasing (or increasing) function of the nonlinear exponent ν. Therefore, we conclude that the correlations between the breaking point and the nonlinear exponent are robust at the qualitative level.

The zero-temperature entropy sath = s(T → 0), computed from the theory, remains nonzero (), which unexpectedly resembles a key feature of amorphous solids. It is well known that the total entropy of an amorphous solid can be decomposed into configurational and vibrational parts, s = sconf + svib. When T → 0, the vibrational entropy svib → 0, while the configurational entropy sconf and the total entropy s remain finite, in contrast to crystals where s → 0. The finite sath indicates the possibility of multiple metastable states, revealing the glassy nature of black hole systems (–). As an essential result from our theory, the zero-temperature entropy, sath ∼ sconf, scales with the shear stress σ as (for large σ)where the exponent ξ = 2/(3 + ν2/3). The prefactor vanishes at T = 0, but the scaling is robust at sufficiently low temperatures (see section S5). Since ξ is positive, the entropy always increases under shear, implying a rejuvenation effect—the system becomes more entropic and less stable under shear. This effect is overlooked by state-following (SF) calculations in mean-field glass theories (), where, by construction, the entropy is kept constant during shear. The increase in entropy under shear deformation is in stark contrast with the behavior observed in ordered crystals, e.g., high-density face-centered cubic HS crystals jam, and therefore, their (vibrational) entropy vanishes (the configurational entropy sconf = 0 in crystals), under shear (). It suggests that our holographic models might share more commonalities with amorphous rather than crystalline systems, in agreement with recent related considerations ().

Simulations of a granular matter model

We perform computer simulations of a 3D granular matter model, which consists of polydisperse spheres interacting via frictionless short-range repulsive forces (see Materials and Methods). The system is compression quenched from initial configurations that are equilibrated at φg using an efficient swap algorithm (see Materials and Methods) to zero temperature where it jams randomly at φj. Thus, φg can be understood as a glass transition density, quantifying the degree of annealing (the larger the φg value, the deeper the annealing). Thanks to the swap algorithm, we are able to prepare ultrastable states corresponding to deep annealing. This is the key reason to observe substantially stronger nonlinear elasticity, compared to previous numerical studies (). Additional data are provided in the Supplementary Materials for 2D models without (section S7) and with (section S8) friction and systems mechanically annealed by cyclic shear (section S6), confirming that the reported behavior is qualitatively insensitive to friction, dimensionality, and preparation protocols.

The above compression quenching procedure generates isotropic configurations at zero temperature and a finite pressure P0 > 0, where P0 = P(ε = 0) characterizes the distance to isotropic jamming, P0 ∼ φ0 − φj. These isotropic configurations serve as unstrained (ε = 0) reference states to athermal quasi-static shear with simple strain deformations under constant volume conditions (see Materials and Methods). For a small P0 = 10−2, the stress-strain curves in Fig. 3A display clear nonlinear elasticity, following the scaling law in Eq. 1. The exponent ν > 1 reveals shear-hardening behavior (G increases with ε), which is equivalent to a dilatancy effect under a constant pressure condition. The stress-strain curve is reversible before yielding in cyclic shear (fig. S13), confirming that the observed nonlinearity has a dominating elastic origin.

Fig. 3.

Correlation between nonlinear elasticity and yielding obtained from granular simulations.

(A) Stress-strain curves for P0 = 10−2 and a few different φg values. The yielding point {εY, σY} is estimated at σY = cσmax, where cσmax is the maximum stress and c = 0.98 (see section S6 for other choices of c and the discussion therein). Fitting the data (solid lines) according to Eq. 1 gives the exponent v (see section S6 for a discussion on the fitting). (B) εY and σY as functions of v. (C) Dependence of on the rescaled degree of annealing , where φMCT is the mode-coupling theory (MCT) transition density. Besides the EFT and simulation results of athermal SSs (P0 = 10 and 10−2) obtained in this study, we also present simulation data of 3D thermal HSs () and the SF theoretical result (divided by 5) (). In addition, we plot the theoretical result from the elastoplastic model (EPM), where the x axis represents (Ac − A)/A, with A being the degree of annealing and Ac being the critical point (). Theories are indicated by lines and simulations by line points. It is clear that εY decreases with the degree of annealing due to shear-hardening nonlinearity (red); in all other cases, shear hardening is absent, and εY increases or remains nearly constant (blue).

The yield stress σY increases, and the yield strain εY decreases, with ν (see Fig. 3B), qualitatively consistent with our theoretical predictions in Fig. 2B. Note that near jamming, the unstrained shear modulus G0 ≈ 0, and the unstrained bulk modulus K0 is finite (); thus, we have set Poisson’s ratio close to 1.0 in the theory (see Fig. 2B).The negative correlation between the yielding strain εY and the degree of annealing φg is in contradiction with previous simulation (, ) and theoretical results (, ), where εY increases or is nearly unchanged with the degree of annealing (see Fig. 3C). Note that in all those previous cases, the nonlinear shear-hardening effect is absent. To confirm this point, additional simulations are performed at larger P0, where shear hardening is compensated by strong plasticity, and εY correspondingly becomes either independent of φg (fig. S14) or slightly increases with φg (fig. S16B). We thus conclude that the left-shifted yielding peak is caused by nonlinear corrections to the elasticity and therefore cannot be captured by linear elasticity theories (). Shear hardening disappears in poorly annealed (small φg; see Fig. 3A) or overcompressed (large P0; see figs. S14 and S16B) systems, which explains why it was not observed in many previous simulations (, ). Deep annealing and isostaticity (the coordination number Z = 2d) could be two key ingredients to this effect.

Next, we investigate the change in entropy during shear in simulations. Figure 4A shows that the jamming density φj(ε), at which P vanishes upon decompression for the given ε, decreases monotonically with ε. According to the mean-field glass theory, φj is positively correlated with the configurational entropy sconf of glass states (). Thus, Fig. 4A already suggests an increase in entropy induced by shear, qualitatively consistent with our theoretical prediction, Eq. 3.

Fig. 4.

Entropy evolution of simulated granular matter under shear.

(A) Jamming density φj versus strain ε for a few different φg values. (B) Configurational entropy sconf from () and (C) Edwards entropy ΔsEd as a function of stress σ, fitted to Eq. 3 (solid lines; see section S6 for a discussion on the fitting). The fitting parameters ξ and ξEd as a function of the nonlinear exponent ν are plotted in the inset of (C). Error bars represent the SE of the fitting exponents.

Because the direct computation of the zero-temperature entropy is impractical for systems of thousands of particles (), we estimate it indirectly by the following two independent ways. (i) In the first approach, we assume that the zero-temperature configurational entropy at φj is proportional to the finite-temperature configurational entropy sconf(φg) of the corresponding parent liquid state at φg. Under this assumption, we collect the data of sconf(φg) for the same model from () (estimated by taking the difference between the total entropy s and the vibrational entropy svib) and φj(φg) from (). Combining them with φj(ε) in Fig. 4A and σ(ε) in Fig. 3A gives sconf(σ) in Fig. 4B. The numerical value of ξ is estimated by fitting the data to Eq. 3. (ii) In the second approach, we consider the Edwards entropy, constructed using the framework of Edwards statistical mechanics of granular matter, which is a generalization of Boltzmann statistical mechanics to nonequilibrium, athermal systems (, ). The Edwards entropy is computed on the basis of the fluctuations of the local Voronoi volumes () (see section S6 for details) and has been applied to unstrained, isotropic granular systems in both simulations (–) and experiments (). Here, we apply the method to anisotropic systems under simple shear, using the Lees-Edwards boundary conditions (). The change in Edwards entropy ΔsEd per particle under shear [where we have taken the jammed states at the lowest jamming density, or the J-point density (), φJ ≈ 0.655 (), as the reference] is plotted in Fig. 4C and fitted to Eq. 3 to obtain ξEd. The two approaches give close results on the exponents (see Fig. 4C, inset). The numerical values of ξ and ξEd are of the same order as the theoretical prediction (Fig. 2C, inset) and decay similarly with the nonlinear exponent ν. We point out that the agreement remains mainly qualitative because of the phenomenological nature of our theories.

DISCUSSION

Our results shed lights on the correct and concrete physical interpretation of the theoretical models at hand (, ). In particular, the behavior of the entropy under shear suggests that the holographic models considered () are phenomenologically closer to amorphous solids rather than crystalline systems. On the other hand, the current versions of theories do not incorporate marginal stability and isostaticity and therefore cannot properly describe the jamming transition. The theories assume that the system is always stable, instead of marginally stable, before the breaking point and thus are applicable only to the regime away from jamming where the comparison to simulations is made. It would be extremely interesting to incorporate the necessary microscopic information and extend our theoretical formulation toward the jamming transition to explore its critical properties.

Our simulations show that shear hardening universally exists in 2D/3D, frictionless/frictional, and thermally/mechanically annealed granular models, which suggests that the phenomenon could be directly relevant to a number of experimental systems (see section S9) (–). In particular, ultrastable shear-jammed granular materials were realized in a recent experiment (), making a direct test of the discussed correlations possible in the laboratory.

The approach presented here can be generalized to study more complex nonlinear behaviors. For example, very close to jamming, the stress-strain curve of granular matter typically displays three consecutive elastic regimes with shear (): linear (ν = 1), shear softening (ν < 1), and shear hardening (ν > 1). In this study, we mainly focus on the third regime (shear hardening) at large strains. We expect the entire stress-strain curve to be captured by a generalized two-potential gravity theory. Further extensions can be made by taking into account the effects of finite temperature, finite shear rates, and viscoelasticity, as well as nonequilibrium relaxational dynamics.

MATERIALS AND METHODS

Effective field theory

Following (), the EFT description is implemented in terms of d = 2 scalar fields (I = 1…d), ϕ(t, x) = 〈ϕ〉 + δϕ(t, x), which play the role of co-moving coordinates. The variations from their equilibrium positions δϕ(t, x) coincide with the displacement fields used in standard elasticity theory. The EFT action is built only in terms of the derivatives of the fields, reflecting the invariance under the shift symmetry ϕ → ϕ + a with a constants. This symmetry follows from identifying the effective fields with phonons, which can be understood as Goldstone bosons for the spontaneously broken translations ()—fluctuations around the nontrivial vacuum expectation values 〈ϕ〉 = x, with zero energy.

The only independent scalar objects that can be constructed out of the derivatives of the displacement fields (in d = 2) are X ≡ Tr [∂μϕμϕ] and Z ≡ Det [∂μϕμϕ], where the index μ = (t, x) collectively describes the set of spacetime coordinates. The scalar potential then becomes a generic function V(X, Z), and it is the only ingredient that must be provided in the theory. In most of our discussion, we will consider a power-law potential, which is reminiscent of the so-called hyperelastic models and takes the form V(X, Z) = X(. Here, A and B are two phenomenological model parameters that cannot be fixed without first-principle calculations. The stress tensor of the theory can be derived using the standard quantum field theory variational prescription and the dispersion relation of the low-energy excitations by computing the action for the fluctuations around equilibrium at second order. More technical details are presented in sections S1 and S2.

Gravity theory

The gravitational description takes advantage of the so-called holographic duality in which a gravity system in a (d + 2)–dimensional spacetime is mapped to a many-body system “defined” on its (d + 1)–dimensional boundary. The boundary system is referred to as a “hologram” of the bulk. Thanks to the duality, introducing dissipative mechanisms and the effects of temperature becomes easy. More precisely, we use the bottom-up holographic duality in the large N limit. Here, the parameter N is interpreted as the number of effective degrees of freedom in the dual-field theory (). Large N in the gravity description corresponds to L ≫ l, with L as the AdS length scale and l as the Planck scale. This framework is more general than the original, string theory–inspired AdS/conformal field theory correspondence (), and it applies also to systems that are not critical (i.e., without conformal invariance). Concretely, the infrared physics of our gravitational system is governed by a nonrelativistic geometry with an AdS2 fixed point, which does not enjoy the conformal group. This type of geometry shares some features with amorphous systems and spin glasses and, in particular, has a finite entropy at zero temperature (). Last, in the large N limit (or equivalently L ≫ l), all quantum loops in the boundary field theory are suppressed by factors of 1/N, and the corresponding physics is effectively classical. All these assumptions justify the validity of our framework to describe the nonlinear rheology of classical particles in amorphous systems. The framework thus allows us to perform simple and robust computations of important physical observables such as the entropy. Our computations are based on the nonlinear generalization (, ) of the known holographic axion model ().

The benchmark model uses a bulk potential form, , where and Z = Det[J] with J ≡ ∂μφμφ. The index μ here spans a 4D spacetime (t, x, u) with x ≡ (x, y). While the bulk potential W(X, Z) has a similar form as the one in EFT, the connection between the bulk potential and the dual EFT potential is very nonlocal and subtle. To avoid any clutter, we will always use different symbols for bulk quantities and EFT ones. The stress tensor of the dual-field theory can be extracted from the gravitational action by using the standard holographic dictionary, while that of the entropy can be extracted using the famous Bekenstein-Hawking area law. See sections S3 to S5 for more details.

Granular model of 3D frictionless SSs

The model (, , ) is composed of N = 8000 SSs, with a diameter distribution P(D) ∼ D−3, where Dmin ≤ D ≤ Dmin/0.45. Two spheres interact via a potential , only if their separation r is less than their mean diameter D = (D + D)/2; otherwise, V = 0. The unit of length is the average diameter of all particles, the unit of energy is 103 × k, and all particles have the same unit mass. Simulation data are averaged over 96 independent samples.

Swap algorithm

The SS granular configurations are quenched from equilibrium states at φg (, ). To prepare these equilibrium configurations, we use HS potential and a very efficient swap Monte Carlo (MC) algorithm (). The HS configurations are generated by the hybrid of two different kinds of moves, the standard moves and the swap MC moves of exchanging the diameters of two randomly picked spheres. The swap moves are accepted only if the resulting configuration does not violate the HS constraint. Such moves help particles break out of cages formed by their neighbors and diffuse freely and hence facilitate the equilibration procedure. With the aid of this algorithm, we obtain equilibrium HS configurations over a wide range of φg.

Simulation protocol of compression quench

Once an equilibrium HS configuration is obtained by the swap algorithm, we switch off the temperature and switch to the SS potential. The system is then quenched to jamming density by a series of athermal quasi-static compression and decompression (). If the system is jammed, i.e., the energy per particle is larger than 10−13, then the system is decompressed; otherwise, it is compressed. During each compression (decompression) step, we instantaneously inflate (deflate) the spheres to increase (decrease) the packing density by δφ and then minimize the total potential energy using the fast inertial relaxation engine (FIRE) algorithm (). The energy minimization stops when the averaged residual force per particle is less than 10−11, which means that the configuration has reached a mechanically stable state. The initial δφ = 10−4; it is then reduced by a factor of 2, whenever the state alters from jammed to unjammed (in the meanwhile, we switch from decompression to compression), or vice versa. This procedure stops when δφ < 10−6 and the system is jammed. The jammed configurations are quasi-statically compressed by δφ = 10−5 to the target pressure P0 to obtain unstrained configurations at ε = 0.

Simulation protocol of athermal quasi-static shear

Shear is performed under the athermal quasi-static and Lees-Edwards boundary conditions (). Each shear step (δε = 10−4) involves an affine transformation of coordinates and then an energy minimization using the FIRE algorithm (). The energy minimization stops when the average residual force per particle is less than 10−11. The pressure P and stress σ (shear is applied in the xy plane) are calculated from the virial formulawhere r and f are the center-to-center vector and force between particles k and l (r and f are the x and y components), V is the volume of simulation box, and 〈kl〉 stands for all contacting pairs.