Mean field stability map of hard-sphere glasses

TL;DR

This paper maps the phase diagram of hard-sphere glasses in infinite dimensions, revealing the conditions for shear jamming and yielding, and how these boundaries evolve with different preparation densities.

Contribution

It provides a comprehensive theoretical analysis of the stability and phase transitions of hard-sphere glasses under shear in infinite dimensions, expanding previous models.

Findings

Identification of shear jamming and yielding lines in the phase diagram.

Characterization of how these lines evolve with glass preparation density.

Mapping of the stability boundaries for hard-sphere glasses under shear.

Abstract

The response of amorphous solids to an applied shear deformation is an important problem, both in fundamental and applied research. To tackle this problem, we focus on a system of hard spheres in infinite dimensions as a solvable model for colloidal systems and granular matter. The system is prepared above the dynamical glass transition density, and we discuss the phase diagram of the resulting glass under compression, decompression, and shear strain, expanding on previous results [P. Urbani and F. Zamponi, Phys.Rev.Lett. 118, 038001 (2017)]. We show that the solid region is bounded by a "shear jamming" line, at which the solid reaches close packing, and a "shear yielding" line, at which the solid undergoes a spinodal instability towards a liquid, flowing phase. Furthermore, we characterize the evolution of these lines upon varying the glass preparation density. This work aims to provide a general overview of yielding and jamming phenomena in hard-sphere systems by a systematic exploration of the phase diagram.

Figures

Equations

$$Z[{\Delta }_r,\hat{\phi },\gamma |Y,{\hat{\phi }}_g]=\int dX{e}^{-\beta V_{\gamma }[X;\hat{\phi }]}\delta \left({\Delta }_r-\Delta (X,Y)\right),$$

$$\begin{array}{rl}{s}_g& =\frac{{\Delta }_r-\Delta /2}{\Delta }+\frac{1}{2}\log \left(\frac{\pi e\Delta }{d^2}\right)+\\ & +\frac{{\hat{\phi }}_g}{2}{\int }_{-\infty }^{\infty }dh{e}^h{q}_{\gamma }\left(2{\Delta }_r-\Delta ;h\right)f(h).\end{array}$$

$$q(\Delta ;h)=\Theta \left(\frac{h+\Delta /2}{\sqrt{2\Delta }}\right),$$

$$q_{\gamma }(\Delta ;h)=\int \mathcal{D}\zeta q(\Delta +{\gamma }^2{\zeta }^2;h),$$

$$\begin{array}{rl}2{\Delta }_r=& \Delta +{\hat{\phi }}_g{\Delta }^2\times \\ & \times {\int }_{-\infty }^{\infty }dh{e}^h\frac{\partial }{\partial \Delta }\left[q_{\gamma }(2{\Delta }_r-\Delta ;h)f(h)\right],\\ \frac{2}{\Delta }=& -{\hat{\phi }}_g{\int }_{-\infty }^{\infty }dh{e}^h\left[\frac{\partial }{\partial {\Delta }_r}{q}_{\gamma }\left(2{\Delta }_r-\Delta ;h\right)\right]f(h),\end{array}$$

$$\frac{1}{\hat{\phi }}=-\Delta {\int }_{-\infty }^{\infty }dh{e}^h\log \left[q(\Delta ;h)\right]\frac{\partial q(\Delta ;h)}{\partial \Delta }\equiv \mathcal{F}(\Delta ).$$

$$\frac{1}{{\hat{\phi }}_{\mathrm{d}}}=\underset{\Delta }{\max }\mathcal{F}(\Delta ).$$

$${\hat{p}}_g=-\frac{{\hat{\phi }}_g}{2}{\int }_{-\infty }^{\infty }dh{e}^h{q}_{\gamma }\left(2{\Delta }_r-\Delta ;h\right)\frac{\partial f(h)}{\partial \eta }$$

$$\underset{\Delta \to 0}{\lim }\Delta \log q(\Delta ,h)=-\frac{h^2}{2}\theta (-h),$$

$$\begin{array}{rl}& {\Delta }_r=\frac{{\hat{\phi }}_g}{4}{\int }_{-\infty }^{\eta }dh{e}^h{q}_{\gamma }(2{\Delta }_r;h)(h-\eta {)}^2,\\ & 1=\frac{{\hat{\phi }}_g}{4}{\int }_{-\infty }^{\eta }dh{e}^h\left[\frac{\partial }{\partial {\Delta }_r}{q}_{\gamma }\left(2{\Delta }_r;h\right)\right](h-\eta {)}^2.\end{array}$$

$$H({\Delta }_r)={\Delta }_r-\frac{{\hat{\phi }}_g}{4}{\int }_{-\infty }^{\eta }dh{e}^h(h-\eta {)}^2{q}_{\gamma }\left(2{\Delta }_r;h\right)$$

$$M=\left(\begin{array}{cc}\frac{{\partial }^2{s}_g}{\partial \Delta \partial \Delta }& \frac{{\partial }^2{s}_g}{\partial \Delta \partial {\Delta }_r}\\ \frac{{\partial }^2{s}_g}{\partial {\Delta }_r\partial \Delta }& \frac{{\partial }^2{s}_g}{\partial {\Delta }_r\partial {\Delta }_r}\end{array}\right),$$

$$\begin{array}{rl}{\lambda }_{\text{R}}=1& -\frac{{\hat{\phi }}_g}{2}{\Delta }^2{\int }_{-\infty }^{\infty }dh{e}^h{q}_{\gamma }(2{\Delta }_r-\Delta ;h)f^{\prime \prime }(h{)}^2.\end{array}$$

$$\lambda =\frac{{\hat{\phi }}_g{\int }_{-\infty }^{\infty }dh{e}^h{q}_{\gamma }(2{\Delta }_r-\Delta ;h)f^{\prime \prime \prime }(h{)}^2}{\frac{4}{{\Delta }^3}+2{\hat{\phi }}_g{\int }_{-\infty }^{\infty }dh{e}^h{q}_{\gamma }(2{\Delta }_r-\Delta ;h)f^{\prime \prime }(h{)}^3}.$$

$$\dot{\Delta }(\lambda )$$

$$A(h)$$

Full text

I Introduction

A glass can be essentially described as an amorphous solid characterized by an enormous slowing down of the diffusive motion. This slowing down is accompanied by a spectacular increase of the relaxation time and the viscosity. The glass transition is empirically defined as the temperature below which (or density above which) the material becomes too viscous to undergo a diffusive motion on experimentally relevant timescales, and is then frozen in an out-of-equilibrium state.

Because the glass transition is not a sharp phase transition, but rather an anthropocentric concept, and the glass is an out-of-equilibrium state, protocol preparation plays a strong role in determining glass properties. As a result, glassy materials display peculiar features distinct from ordinary solids. For instance, the stress-strain curves of the glass strongly depend on its degree of annealing Rodney et al. (2011); Bonn et al. (2017); Ozawa et al. (2018). Furthermore, the response of glasses to an external weak stress is characterized by large deformations and very complex phenomena that can be generally explained by invoking marginal stability Liu and Nagel (2010); Müller and Wyart (2015). The marginal stability condition is responsible for unusual properties ranging from an abundance of zero-frequency modes Wyart et al. (2005) to the emergence of long-range correlations in vibrational dynamics Berthier et al. (2016).

In the following, we will focus on a thermodynamic approach to structural glasses, originally proposed and developed by Kirkpatrick, Thirumalai and Wolynes Kirkpatrick and Wolynes (1987); Kirkpatrick and Thirumalai (1987); Kirkpatrick et al. (1989), and known as the Random First Order Transition (RFOT) theory. Within this theory, glasses are considered as extremely long-lived metastable states, essentially described in terms of vibrations around an amorphous reference structure Singh et al. (1985), which encodes the protocol dependence. Because here we consider glasses obtained by slow cooling of a supercooled liquid, the reference structure is taken to be a typical equilibrium configuration at the state point where the liquid fell out-of-equilibrium Rainone et al. (2015); such a construction is called the Franz-Parisi potential method, originally developed in spin glass theory Franz and Parisi (1995, 1997); Barrat et al. (1997); Krzakala and Zdeborová (2010); Zdeborová and Krzakala (2010); Krzakala and Zdeborová (2013).

In the abstract limit of infinite spatial dimension $d$, the method is exact Charbonneau et al. (2017), but it can also be used in two or three dimensions in an approximate way Mézard and Parisi (1996); Cardenas et al. (1998); Yoshino and Mézard (2010); Parisi and Zamponi (2010). Results for infinite-dimensional hard spheres Rainone et al. (2015); Urbani and Zamponi (2017) show that the application of a strain $\gamma$ can result in very different and highly complex scenarios. Generically, at small $\gamma$, an elastic-like regime is observed, which is identified by a linear dependence between the response stress $\sigma$ and the applied strain, i.e. $\sigma\sim\mu\gamma$ with $\mu$ being the shear modulus Yoshino and Zamponi (2014). Upon increasing the strain, beyond a Gardner transition Gardner (1985); Kurchan et al. (2013); Charbonneau et al. (2014); Rainone et al. (2015); Rainone and Urbani (2016); Biroli and Urbani (2016, 2018), the response turns out to be characterized by intermittent drops or avalanches, related to collective rearrangements of the whole system, which is nevertheless still solid Rainone and Urbani (2016); Franz and Spigler (2017). Depending on the preparation and target densities, two scenarios are possible at larger strain. Either the system eventually looses its solidity and starts to flow, undergoing shear yielding, or it can reach a jamming (or close packing) point, where a rigid network of contacts is formed and the system can no longer be strained Urbani and Zamponi (2017). The shear jamming and the isotropic jamming (under compression without shear) transitions fall in the same universality class with the same critical exponents regulating the distributions of effective forces and gaps Urbani and Zamponi (2017). Furthermore, the theory has been extended to hard sphere with a short range attraction Sellitto and Zamponi (2013), to investigate a peculiar two-step yielding transition that characterizes colloidal systems Altieri et al. (2018). These predictions, obtained in the mean-field $d\rightarrow\infty$ limit, have been partially confirmed by extensive numerical simulations in $d=3$ Charbonneau et al. (2012, 2015); Baity-Jesi et al. (2017); Jin and Yoshino (2017); Jin et al. (2018) and experiments on a granular material in $d=2$ Seguin and Dauchot (2016).

In this paper, we will complete and extend the hard-sphere results obtained in Rainone et al. (2015); Urbani and Zamponi (2017). We will explore a wider regime of densities, in order to obtain a complete phase diagram. In particular, we will show that the shear yielding line terminates on the melting point at low density and vanishing shear strain, while on the other side it touches the shear jamming line at a critical point identified in Urbani and Zamponi (2017). We will study the preparation density evolution of these lines, and we will provide further details on the nature of the Gardner transition already studied in Rainone and Urbani (2016).

The article is organized as follows: in Sec. II, we introduce the mathematical formalism for a system of hard spheres in $d$ dimensions, with $d\rightarrow\infty$; in Sec. III, we describe the relevant phase transitions, including the jamming, yielding and Gardner transitions. Then, in Sec. IV we present our results, based on the determination of the glass phase diagrams, under compression or decompression and quasi-static shear strain. In Sec. V we conclude presenting some perspectives for future investigations.

II State following procedure

We consider a system of $N$ identical hard spheres in $d$ dimensions, which turns out to be a very convenient model for structural glasses and granular matter (monodisperse spheres are a good glass former as long as $d\geq 4$ Skoge et al. (2006); van Meel et al. (2009)). We then consider the limit $d\rightarrow\infty$, which makes the problem analytically tractable and provides an exact thermodynamical Kurchan et al. (2012); Rainone et al. (2015) and dynamical solution Maimbourg et al. (2016); Agoritsas et al. (2019). Note that, in order to obtain a finite $d\rightarrow\infty$ limit, one has to rescale the packing fraction $\varphi$, which is the fraction of volume occupied by the spheres, introducing $\widehat{\varphi}=2^{d}\varphi/d$. The other observables have to be rescaled as well, as detailed below. For hard spheres, thermal fluctuations are irrelevant, and the only control parameter in equilibrium is density.

To perform our analysis we employ a state following protocol Franz and Parisi (1995, 1997); Barrat et al. (1997); Krzakala and Zdeborová (2010); Zdeborová and Krzakala (2010); Krzakala and Zdeborová (2013), which describes glass formation during a very slow liquid cooling (see Rainone et al. (2015) for a detailed explanation in the context of structural glasses). The state following method works for very large relaxation times, which is actually the case in infinite dimensions where the equilibrium relaxation time diverges upon approaching the dynamical transition at packing fraction $\widehat{\varphi}_{\rm d}$ Maimbourg et al. (2016); Charbonneau et al. (2017) (or critical transition in Mode-Coupling Theory Götze (1999)). Several properties of glassy systems can be re-obtained in this framework, such as the stress overshoot and the shear modulus behaviors, as well as the dependence of the pressure and the specific volume on the cooling rate Rainone et al. (2015); Charbonneau et al. (2017).

Because the method has been discussed in details elsewhere Rainone et al. (2015); Charbonneau et al. (2017), we only provide here a brief summary of the main equations. We consider a system equilibrated at a value of density $\widehat{\varphi}_{g}>\widehat{\varphi}_{\rm d}$, where diffusion is frozen and an enormous number of metastable glassy states emerge, each with the same pressure as the equilibrium liquid prepared at $\widehat{\varphi}_{\rm d}$. In this region, the system is trapped in one such metastable state, whose evolution under compression/decompression and/or an applied shear strain is then studied Rainone et al. (2015).

To realize this construction in practice, we first extract an equilibrium configuration, $Y=\{\bf{y}_{i}\}$, from the Boltzmann-Gibbs measure at $\widehat{\varphi}_{g}$, and then we consider another configuration, $X=\{\bf{x}_{i}\}$, which is constrained to be close to the former. Their mutual mean-squared displacement is fixed by the condition $\Delta(X,Y)=\Delta_{r}$, where $\Delta(X,Y)=\frac{d}{N}\sum\limits_{i=1}^{N}|{\bf{x}}_{i}-{\bf{y}}_{i}|^{2}$. The resulting partition function then becomes:

[TABLE]

where $V_{\gamma=0}[X;\widehat{\varphi}]$ is a hard-core potential, at a density $\widehat{\varphi}$ that can differ from $\widehat{\varphi}_{g}$, thus realizing a compression or decompression. An applied shear strain $\gamma$ on configuration $X$ is described by a modified potential $V_{\gamma}[X;\widehat{\varphi}]$, in which a linear transformation is applied to configuration $X$ to describe the box straining, as derived in Yoshino and Mézard (2010); Rainone et al. (2015). In this way, an equilibrium theory for the rheology of amorphous solids is obtained, by keeping the reference configuration $Y$ unstrained and applying a strain on the other, constrained replica $X$. The physical idea behind this construction is that in the solid phase, after the application of the affine shear strain, the system relaxes via non-affine deformation to a stationary state which is close enough to the original state; the constrained replica $X$ describes precisely this long-time pseudo-equilibrium state. See Yoshino and Mézard (2010); Yoshino (2012) for a detailed discussion.

Unfortunately, the computation of the partition function in Eq. (1), as well as that of the associated free energy, is impossible because the constraint $\Delta_{r}=\Delta(X,Y)$, for fixed $Y$, explicitly breaks the translational invariance of the measure for $X$: the equilibrium configuration $Y$ plays the role of a quenched disorder. Yet, the free energy is self-averaging with respect to $Y$, which means that for given $\widehat{\varphi}_{g}$ and $(\widehat{\varphi},\gamma)$, the free energy for given $Y$ is equal to that averaged over $Y$ with probability one in the thermodynamic limit. One can then take the equilibrium average over $Y$, which can be done via the replica method Franz and Parisi (1995); Mézard et al. (1987). This procedure restores translational invariance and makes the problem tractable Rainone et al. (2015). One then introduces $s$ replicas of configuration $X$ and eventually considers the analytic continuation to $s\rightarrow 0$. In order to perform this analytic continuation, an ansatz has to be made on the behavior of the replicas under the permutation symmetry Mézard et al. (1987), the simplest case being the so-called replica symmetric (RS) ansatz, in which the replicas are symmetric under any permutation Mézard et al. (1987). A single order parameter, $\Delta=\Delta(X_{a},X_{b})$, has to be introduced to describe the mean-squared displacement between any pair $a,b$ of replicas of $X$. Once the analytic continuation is performed, the averaged free energy should be extremized with respect to the variational parameters $\Delta,\Delta_{r}$. Additional technical details can be found in Rainone et al. (2015); Rainone and Urbani (2016).

In the hard-sphere model considered thus far, energy vanishes and only entropy is relevant. The RS glass entropy (per particle and per spatial dimension, i.e. per degree of freedom) can be exactly computed in the limit $d\rightarrow\infty$, and is given by

[TABLE]

The function $q(\Delta;h)$ reads for $\gamma=0$

[TABLE]

with $\Theta(x)=\frac{1}{2}\left(1+\text{erf}(x)\right)$, and the shorthand notation $f(h)=\log q\left(\Delta;h-\eta\right)$ is used in Eq. (2) and in the following. The parameter $\eta=\log(\widehat{\varphi}/\widehat{\varphi}_{g})$ denotes the logarithm of the relative compression (if positive) or decompression (if negative). The function $q_{\gamma}(\Delta;h)$, which encodes the dependence on a shear strain $\gamma$, is given by

[TABLE]

where $\mathcal{D}\zeta$ denotes the Gaussian measure with zero mean and unit variance.

The differentiation of Eq. (2) with respect to the self mean-square displacement, $\Delta$, and the relative mean-square displacement, $\Delta_{r}$, provides two equations

[TABLE]

which can be solved iteratively to determine $\Delta,\Delta_{r}$. Because the existence of a unique solution is not guaranteed (see Altieri et al. (2018) for an explicit example), we start from the equilibrium solution and follow it adiabatically upon increasing $\gamma$.

At equilibrium $(\widehat{\varphi}=\widehat{\varphi}_{g},\gamma=0)$, we have $\Delta=\Delta_{r}$, and Eqs.(5) reduce to a single equation

[TABLE]

The function $\mathcal{F}(\Delta)\geq 0$ vanishes at $\Delta=0$ and $\Delta\rightarrow\infty$, and it therefore has a maximum in between, whole value provides the dynamical transition density $1/\widehat{\varphi}_{\rm d}=4.8067$:

[TABLE]

For $\widehat{\varphi}<\widehat{\varphi}_{\rm d}$, Eq. (7) has no solution, which translates into the impossibility of finding a stable glass phase: the system is liquid. Instead, for $\widehat{\varphi}>\widehat{\varphi}_{\rm d}$, a solution for $\Delta$ exists in equilibrium. We then need to select the initial glass state at $\widehat{\varphi}_{g}>\widehat{\varphi}_{\rm d}$.

To iterate Eqs.(5) within a given precision for $\widehat{\varphi}\neq\widehat{\varphi}_{g}$ and/or $\gamma>0$, we set as initial condition $\Delta=\Delta_{r}$ equal to the solution of Eq. (6) at $(\widehat{\varphi}=\widehat{\varphi}_{g},\gamma=0)$, we compute numerically the right hand side of Eqs.(5) at the new state point, and we iterate until convergence. Then we slowly change the density, the shear strain or both, and we use the values obtained at the previous step as new guess. Once we have checked the convergence of the above equations and obtained the resulting values of $\Delta$ and $\Delta_{r}$, we can compute the glass entropy, its reduced pressure per degree of freedom $\widehat{p}_{g}=\beta P_{g}/(d\rho)$,

[TABLE]

and all other thermodynamic quantities of interest.

III Phase transitions

We now aim at completely characterizing the phase diagram of a glass prepared at $\widehat{\varphi}_{g}$, and followed in compression/decompression and applied shear strain, $(\widehat{\varphi},\gamma)$. In the following, we briefly recap the main features of the phase transitions that are encountered in this hard-sphere system.

III.1 Jamming

Conventionally, the jamming transition line of hard spheres is identified as the locus of points where the pressure diverges and the mean-square displacement approaches zero Parisi and Zamponi (2010). We can then define a critical shear jamming line $\gamma_{J}(\eta)$, or $\eta_{J}(\gamma)$, in correspondence of which the jamming transition takes place. We take the limit $\Delta\rightarrow 0$ in Eqs. (5) and, using the relation

[TABLE]

we obtain two equations for the jamming line, one for $\Delta_{r}$ and the other for $\eta$:

[TABLE]

Note that these equations are equivalent to imposing that

[TABLE]

and its first derivative $H^{\prime}(\Delta_{r})$ both vanish. Their numerical solution is then quite easy.

Note that, according to the replica symmetric (RS) ansatz we use in this paper, in the jamming limit $\widehat{p}\propto\Delta^{-1}$. However, this prediction disagrees with the correct scaling solution, $\widehat{p}\simeq\Delta^{-\kappa}$, regulated by a universal exponent $\kappa$, whose value is exactly predictable only within a full replica symmetry breaking (RSB) approach Charbonneau et al. (2014). Yet, the RS approximation provides a decent approximation to the location of the jamming critical line. The employment of a full RSB solution might account for the identification of the exact shape of the shear jamming and shear yielding curves in the region where they merge into a critical point, and possibly give rise to a cusp-like trend. This peculiar behavior has been pointed out in recent numerical simulations Jin et al. (2018) reproducing the stability-reversibility map of hard-sphere glasses.

III.2 Yielding

Within mean field theory, a small enough strain applied on a solid gives rise to a static stress, up to a point where the applied stress destabilizes the solid: the solution for $(\Delta,\Delta_{r})$ undergoes a bifurcation and is lost. We call this point a yielding transition Rainone et al. (2015), because it is reasonable to conjecture that at larger shear strains the solid yields and starts to flow, reaching a stationary state in which dynamics is diffusive and the shear stress depends on the shear strain rate $\dot{\gamma}$. This should be, however, tested by solving the dynamical mean field equations Agoritsas et al. (2019).

More precisely, the stress $\sigma(\gamma)$ increases upon increasing $\gamma$, overshoots, and then decreases upon approaching the bifurcation (spinodal) point, where it displays a square root singularity, i.e. $\sigma(\gamma)-\sigma(\gamma_{Y})\propto\sqrt{\gamma_{Y}-\gamma}$. We determine the yielding transition line $\gamma_{Y}(\eta)$, or $\eta_{Y}(\gamma)$, by following the stress until the singularity is observed. Note that the appearance of a spinodal corresponds to having a zero longitudinal mode in the Hessian matrix of the free energy,

[TABLE]

which encodes the fluctuations within the replica symmetric sector. One could determine the yielding point by the condition $\det M=0$, but the calculation of $M$ is rather involved and we do not follow this route here.

Note that comparing this mean-field definition of yielding with experiments or simulations is non trivial, and multiple definitions of yielding are used in realistic systems Bonn et al. (2017). Sometimes the yielding point is identified with the maximum of the stress Koumakis et al. (2012), sometimes with the onset of energy dissipation and the appearance of plastic avalanches Dubey et al. (2016), or with the point where the highest stress drop is observed Ozawa et al. (2018). Whether some signature of the mean field spinodal survives in finite dimensions is a currently debated issue Jaiswal et al. (2016); Procaccia et al. (2017); Popović et al. (2018); Ozawa et al. (2018); Jin et al. (2018). Yet, mean field theory provides a qualitatively correct prediction for the shape of the yielding line in the phase diagram Jin et al. (2018). Also in this case, the RS ansatz is sufficient to obtain a correct estimate of the stress overshoot, and of the yielding spinodal point.

III.3 Gardner transition

Eqs. (5) are iterative equations for the self and the relative mean-square displacements, which have been obtained by assuming that the symmetry under permutation of replicas remains unbroken. This is not a general assumption and should instead be tested at each state point. In other words, we should check if the RS solution is actually a stable local minimum of the free energy.

We can then proceed with the computation of the stability matrix. Thanks to the symmetry of the replica indices, the matrix can be simply decomposed in terms of three independent sectors, the longitudinal one described above (two scalars), the anomalous (vector field) and the replicon (tensorial field). The longitudinal mode can generally be associated with the emergence of spinodal points, while the replicon provides information on instabilities that possibly lead to replica symmetry breaking solutions. We refer the interested reader to the specific literature for the computation of these eigenvalues De Dominicis and Kondor (1983); Temesvári et al. (2002); Rainone and Urbani (2016). We only report here the final expression for the replicon mode,

[TABLE]

A positive replicon, $\lambda_{R}>0$, defines a stable phase where the RS solution holds. Conversely, if the replicon changes sign, the solution becomes unstable and the RS ansatz is no longer correct. The Gardner transition Gardner (1985); Kurchan et al. (2013); Charbonneau et al. (2014) line $\gamma_{\rm G}(\eta)$, or $\eta_{\rm G}(\gamma)$, is then defined by the condition of vanishing replicon eigenvalue, $\lambda_{R}=0$. When $\lambda_{R}<0$, each amorphous state, now called basin, is fragmented into a hierarchy of sub-basins, hence determining a very rough energy landscape Mézard et al. (1987); Charbonneau et al. (2014); Rainone and Urbani (2016). Describing it requires a more complicated parameterization, via a function $\Delta(x)$ in the interval $x\in\left[0,1\right]$, called the Parisi order parameter Mézard et al. (1987). This function encodes the probability distribution of mean-square displacements between replicas, $\Delta=\Delta(X_{a},X_{b})$, via the relation $P(\Delta)=|dx/d\Delta|$. The RS solution corresponds to $\Delta(x)=\Delta$, $\forall x\in\left[0,1\right]$.

III.4 Nature of the Gardner transition

Note that, while the exact solution for $\Delta(x)$ above the Gardner transition density requires solving the RSB equations Rainone and Urbani (2016), which goes beyond the scope of this work, one can evaluate perturbatively its trend close to the Gardner transition, where the replicon is zero Sommers (1985). We can then compute the breaking point $\lambda$, in whose vicinity $\Delta(x)$ becomes non-trivial, and the slope $\dot{\Delta}(\lambda)$. Using these results, one can identify the universality class of the transition, using only information from the RS calculation. The breaking point turns out to be:

[TABLE]

while the slope $\dot{\Delta}$, evaluated at the breaking point, is:

[TABLE]

Based on the values of $\lambda$ and $\dot{\Delta}(\lambda)$, one can classify the Gardner transition in three universality classes, keeping in mind that $\Delta(x)$ must be a decreasing function of $x\in[0,1]$ Charbonneau et al. (2014); Rainone and Urbani (2016), as follows:

• •

If the parameter $\lambda$ takes values in $\left[0,1\right]$, and $\dot{\Delta}(\lambda)<0$, then one has a bona fide continuous transition towards a marginally stable fullRSB phase.

• •

If $\lambda\in\left[0,1\right]$, and $\dot{\Delta}(\lambda)>0$, then one has a bona fide continuous transition towards a non-marginal 1RSB phase.

• •

If, instead, $\lambda>1$, then the transition at $\lambda_{R}=0$ is unphysical, and must be preceded by another form of RSB transition, typically a discontinuous RFOT (exactly as it happens e.g. in the $p$-spin model at low enough external field Cavagna et al. (1999) or in the spherical perceptron model Franz et al. (2017)).

IV Results

Previous studies Rainone et al. (2015); Urbani and Zamponi (2017) used the state following method to analyze hard-sphere glasses in the shear jamming and shear yielding regimes. Ref. Urbani and Zamponi (2017) specifically performed a complete analysis of the phase diagram of a glass prepared at $\widehat{\varphi}_{g}=8$, under an applied compression to density $\widehat{\varphi}\geq\widehat{\varphi}_{g}$ and shear strain $\gamma$. The aim of the current work is to complete the exploration of the phase diagram at $\widehat{\varphi}_{g}=8$, in particular considering the previously unexplored regime in decompression, and to investigate systematically the dependence of the whole phase diagram on $\widehat{\varphi}_{g}$.

IV.1 Glass phase diagram

In Fig. 1 we show the mean-field state following phase diagram in the plane $(\widehat{\varphi},\gamma)$, considering both compression and decompression, for different values of the initial glass density $\widehat{\varphi}_{g}=5.5,8,12$, specifically chosen to be representative of the different physical regimes. Let us focus first on the phase diagram at $\widehat{\varphi}_{g}=8$ and later discuss two other instances, respectively at lower and higher densities.

When the glass is compressed at $\gamma=0$, it exists up to a jamming point at which pressure diverges and $\Delta\rightarrow 0$. When the glass is decompressed at $\gamma=0$, it exists up to a melting point where it undergoes a spinodal instability and melts into the liquid. These two transitions thus delimit the region of existence of the glass at $\gamma=0$. At $\gamma>0$, the jamming point extends into a shear jamming line (blue line in Fig. 1), defined as in section III.1, while the melting point extends into a shear yielding line (red line in Fig. 1), defined as in section III.2. The two lines merge into a critical point (black diamond in Fig. 1), whose location $(\widehat{\varphi}_{c},\gamma_{c})$ can be analytically computed considering that there the system is both jammed ($\Delta=0,\widehat{p}\rightarrow\infty$) and yielded (with a vanishing longitudinal eigenvalue). More precisely, it can be identified by imposing simultaneously that Eq. (11) as well as its first and second derivatives vanish, $H(\Delta_{r})=H^{\prime}\left(\Delta_{r}\right)=H^{\prime\prime}\left(\Delta_{r}\right)=0$. While it has been shown that the shear jamming and the isotropic (without shear) jamming transition lines fall in the same universality class Urbani and Zamponi (2017), we cannot make any prediction at this level on the critical point, for which two different kinds of instabilities should be taken into account simultaneously.

Finally, a Gardner transition line is also present (green line in Fig. 1), defined as in Sec. III.3, separating a stable RS phase (below it) from a RSB phase (above it). Note that the whole shear jamming line, as well as the critical point, fall into the RSB regime, and their RS value reported here is therefore only an approximation. In order to better understand the nature of the Gardner transition line, in Fig. 2 we report the values of the breaking point $\lambda$ and the slope $\dot{\Delta}(\lambda)$ computed along the line, as a function of $\gamma$. At small $\gamma$, we have $\lambda<1$ and $\dot{\Delta}(\lambda)<0$, indicating a continuous transition towards a fullRSB phase, according to the discussion of section III.4. Upon increasing $\gamma$, the slope $\dot{\Delta}(\lambda)$ diverges and changes sign, while still $\lambda<1$, indicating that the transition becomes a continuous transition towards a 1RSB phase. At even larger values of $\gamma$, $\lambda$ reaches one, indicating that the transition becomes a RFOT. This behavior is reminiscent of that of several spin glass models Cavagna et al. (1999); Franz et al. (2017).

We now discuss the two other phase diagrams in Fig. 1, for $\widehat{\varphi}_{g}=5.5,12$. The phase diagram for $\widehat{\varphi}_{g}=12$ mostly differs from the one at $\widehat{\varphi}_{g}=8$ for the relative location of the different transitions. In particular, the shear jamming and the Gardner transition lines appear much closer, with a visible shrinkage of the region delimited by the green and the blue lines as $\gamma\rightarrow 0$. Furthermore, the critical point $(\widehat{\varphi}_{c},\gamma_{c})$ has now moved at densities smaller than the preparation one, $\widehat{\varphi}_{c}<\widehat{\varphi}_{g}$, indicating that an applied shear on the initial equilibrium configuration would lead to shear jamming. The behavior of the breaking point $\lambda$ and slope $\dot{\Delta}(\lambda)$ are also qualitatively similar to the $\widehat{\varphi}_{g}=8$ case.

The phase diagram for $\widehat{\varphi}_{g}=5.5$ presents a more important qualitative difference. In fact, upon decreasing $\widehat{\varphi}_{g}$, the critical point moves to $\gamma_{c}=0$ and disappears. The shear jamming line correspondingly disappears, leaving a single shear yielding line that originates from the melting point in decompression at $\gamma=0$ and ends into an unphysical spinodal point in compression at $\gamma=0$ Rainone et al. (2015). This unphysical spinodal point has been shown to disappear once RSB is taken into account Rainone and Urbani (2016), and as a consequence the behavior beyond the Gardner transition line (shaded gray area in Fig. 1) might be qualitatively different once RSB is taken into account. In this case, the value of $\lambda$ (Fig. 2) remains always smaller than one, indicating that the Gardner transition is always continuous.

IV.2 Evolution of the glass phase diagram as a function of the preparation density $\widehat{\varphi}_{g}$

We now give a broader overview of the evolution of the phase diagram of Fig. 1 upon changing $\widehat{\varphi}_{g}$. In Fig. 3(left), we plot the critical values of density that define the phase diagram of Fig. 1 as a function of $\widehat{\varphi}_{g}$. This $(\widehat{\varphi}_{g},\widehat{\varphi})$ phase diagram then shows the evolution with density $\widehat{\varphi}$ and shear strain $\gamma$ of a glass prepared at $\widehat{\varphi}_{g}$.

For a given $\widehat{\varphi}_{g}\geq\widehat{\varphi}_{\rm d}=4.8067\ldots$, the glass exists in the region of $\widehat{\varphi}\in[\widehat{\varphi}_{m}(\widehat{\varphi}_{g}),\widehat{\varphi}_{J}(\widehat{\varphi}_{g})]$, delimited from above by the jamming (blue) line and from below by the melting (red dot-dashed) line. In between these two lines, we plot a few other important values of density: the Gardner transition point $\widehat{\varphi}_{G}$ at $\gamma=0$ (green line), the unphysical spinodal point at $\gamma=0$ (blue dashed line), and the value $\widehat{\varphi}_{c}$ (red line) of the density corresponding to the critical point located at $(\widehat{\varphi}_{c},\gamma_{c})$.

The interpretation of this phase diagram is the following. The glass state is prepared in equilibrium at $\widehat{\varphi}=\widehat{\varphi}_{g}$ (black dot-dashed line), and it can be compressed up to $\widehat{\varphi}_{J}$ or decompressed down to $\widehat{\varphi}_{m}$. Above $\widehat{\varphi}_{G}$, the compression brings the system into a RSB phase and in that case an unphysical spinodal can also be met before jamming. If the system is compressed to $\widehat{\varphi}<\widehat{\varphi}_{c}$ and then strained, shear yielding is observed, while if it is compressed to $\widehat{\varphi}>\widehat{\varphi}_{c}$, shear jamming is observed. Note that upon decreasing $\widehat{\varphi}_{g}$ towards $\approx 6.5$, the critical density $\widehat{\varphi}_{c}$ merges with $\widehat{\varphi}_{J}$, and the shear jamming regime disappears.

Note that, when the system yields, a stationary flow state with no memory of previous conditions should be reached at long times Berthier et al. (2000). Fig. 3(left) then suggests that for any density smaller than $\widehat{\varphi}_{c}(\widehat{\varphi}_{g})$, a stationary flow state exists. It is interesting to observe that the line $\widehat{\varphi}_{c}(\widehat{\varphi}_{g})$ first increases and then decreases, thus showing a single maximum. The value $\widehat{\varphi}^{\text{flow}}_{J}=\max_{\widehat{\varphi}_{g}}\widehat{\varphi}_{c}(\widehat{\varphi}_{g})\approx 11$ is then a lower bound to the value of density for which a hard-sphere system still admits a flowing state. This result should be confirmed by locating the jamming point under shear strain for hard sphere dynamics, using the dynamical equations derived in Agoritsas et al. (2019).

In Fig. 3(right), we present the same phase diagram in a slightly different representation, in the $(\widehat{\varphi},1/\widehat{p})$ plane. The glass is prepared on the equilibrium (dot-dashed black) line, at $\widehat{\varphi}=\widehat{\varphi}_{g}$ and $\widehat{p}=\widehat{\varphi}_{g}/2$. Adiabatic following of the glass state selected in equilibrium by changing its density $\widehat{\varphi}$ at $\gamma=0$ gives rise to a pressure-density equation of state for each $\widehat{\varphi}_{g}$ (black full line). The equation of state is delimited by the melting point $(\widehat{\varphi}_{m},1/\widehat{p}_{m})$ at low densities and by the jamming point $(\widehat{\varphi}_{J},0)$ at high densities. Note that in the vicinity of the melting point, we have $\widehat{p}-\widehat{p}_{m}\propto\sqrt{\widehat{\varphi}-\widehat{\varphi}_{m}}$, but the proportionality coefficient changes sign around $\widehat{\varphi}_{g}\sim 10$, from negative at lower $\widehat{\varphi}_{g}$ (hence leading to a pressure undershoot before melting) to positive at higher $\widehat{\varphi}_{g}$ (hence leading to a monotonic pressure-density curve). As in Fig. 3(left), we report the value of the Gardner point at $\gamma=0$ (green line), and the value of the critical point, here reported as $(\widehat{\varphi}_{c},1/\widehat{p}_{c})$ where $\widehat{p}_{c}$ is the pressure of the glass at $\widehat{\varphi}_{c}$ and $\gamma=0$.

V Conclusions and perspectives

In this work we have applied a microscopic mean-field approach to the simplest model of glass former: hard spheres in the $d\rightarrow\infty$ limit. We obtained the phase diagram of the glass states, followed adiabatically in compression or decompression, and under a quasi-static shear deformation. Our analysis completes and extends previous analytical studies of rheological properties of hard-sphere systems in high dimensions Rainone et al. (2015); Rainone and Urbani (2016); Urbani and Zamponi (2017).

Remarkably, the phase diagrams shown in Fig. 1 are qualitatively similar to those obtained in numerical simulations of three-dimensional systems Jin and Yoshino (2017); Jin et al. (2018), for which we also note that for increasing preparation density $\widehat{\varphi}_{g}$, the values of the yield strain $\gamma_{Y}$ and the yield stress increase as well. The main qualitative difference is the shape of the shear jamming and shear yielding curves in the vicinity of the critical point where they merge, which might be explained by RSB effects.

Indeed, while our analysis is restricted to the replica symmetric solution, replica symmetry is broken in important portions of the phase diagram, in particular in the region of high densities close to the shear jamming line and its associated critical point. Within the RS solution, at low preparation density, the shear jamming line disappears and the shear yielding line extends down to $\gamma=0$ giving rise to an unphysical spinodal of the glass under compression. Clearly, a major direction for future research would be to investigate RSB effects within the state following approach, or even better, to understand what is the dynamical behavior of glasses followed in the unstable region Folena et al. (2019).

Another very interesting problem is to identify the exact value of $\widehat{\varphi}^{\text{flow}}_{J}$, defined as the density above which no flowing state of hard spheres under constant shear rate exists. The configurations at $\widehat{\varphi}^{\text{flow}}_{J}$ are jammed and likely to be marginal, i.e. characterized by a non-trivial force distribution with universal exponents Lerner et al. (2012). It would be extremely interesting to better understand, within mean field theory, how the marginal stability condition affects the flow dynamics around the critical density.

Acknowledgments

We warmly thank P. Urbani for several insightful interactions on this subject. We also acknowledge Y. Jin and H. Yoshino for useful exchanges and discussions. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant agreement 723955 - GlassUniversality).