Exact theory of dense amorphous hard spheres in high dimension. II. The high density regime and the Gardner transition

I Introduction

The Random First Order Transition (RFOT) scenario for glasses introduced by Kirkpatrick, Thirumalai and Wolynes Kirkpatrick and Wolynes 1987a; Kirkpatrick and Thirumalai 1987; Kirkpatrick et al. 1989; Wolynes and Lubchenko 2012 proposes that the glass transition is represented – at least at the mean-field level – by a freezing transition similar to the one of the Random Energy Model Derrida 1981, described mathematically by a one-step replica symmetry breaking (1RSB) ansatz Mézard et al. 1987. Although this scenario was first based on an analogy with spin glasses Kirkpatrick and Wolynes 1987a; Kirkpatrick and Thirumalai 1987, it was quickly realized in a pioneering work by Kirkpatrick and Wolynes Kirkpatrick and Wolynes 1987b that the glass transition of $d$-dimensional hard spheres in the limit $d\rightarrow\infty$ could be described within the same framework. Much later Parisi and Zamponi 2006; Parisi and Zamponi 2010, similar results were obtained using the replica method, which also allowed for a detailed description of the glass phase and in particular of the jamming point where the pressure of the glass becomes infinite, corresponding to its close packing – which was therefore called glass close packing (GCP) and is closely related to the random close packing concept introduced by Bernal much earlier Bernal and Mason 1960.

The main advantage of the replica method is that it allows for a unified treatment of both the glass and the jamming transitions, within a simple static RFOT scenario, and it also allows for a partial understanding of dynamical aspects. Moreover, there is hope that the replica method can provide an exact result in the limit $d\rightarrow\infty$. A first step in this direction was performed in the first paper of this series Kurchan et al. 2012, where we have shown that the thermodynamics of hard spheres in the limit of high dimensions may be exactly obtained from the knowledge of the distribution of the two-point correlation function between states, encoded in the Parisi parameter. In the same paper it was shown that once a 1RSB ansatz is made, one recovers exactly the Gaussian replica free energy that was used in Ref. Parisi and Zamponi 2006; Parisi and Zamponi 2010 to derive estimates of the various transitions that characterize the RFOT scenario at the 1RSB level.

In this paper we take a second important step, by investigating the stability of the 1RSB solution towards further levels of replica symmetry breaking. We find that for higher pressures, well above the RFOT one, there is a second transition (a so-called Gardner transition Gardner 1985) leading to a somewhat different physics in that limit, and in particular around the jamming point. We believe that this physics is intimately connected with the peculiar mechanical properties of jammed states of hard spheres, that have been recently characterized in much detail Liu et al. 2011; Van Hecke 2010.

The rest of the paper is organized as follows. We start our presentation by a general discussion of the RFOT scenario, of its connection with the physics of jamming, and of the main new features that are due to the presence of the Gardner transition. This discussion is reported in Sec. II and it is for the moment mostly speculative, although some parts of it have been previously studied in spin glass models. Next, we present our new results, which constitute a first important step to substantiate this picture for hard spheres in the $d\rightarrow\infty$ limit. In Sec. III we provide a proof of the correctness of the Gaussian ansatz for a generic form of the overlap matrix, extending the main result of Ref. Kurchan et al. 2012; in Sec. IV we recall a few important results of Ref. Parisi and Zamponi 2010; Kurchan et al. 2012 that are directly needed here; in Sec. V we present our main results for the Hessian matrix of the 1RSB solution and in particular its so-called replicon eigenvalue, that is responsible for the instability of this solution (i.e. the Gardner transition); in Sec. VI we discuss the cubic terms in the expansion around the 1RSB solution and from them we extract the dynamical exponents that characterize the glass transition; in Sec. VII we present an approximate calculation to obtain an order of magnitude for the Gardner transition pressure in finite $d$; in Sec. VIII we summarize and draw our conclusions.

II A general RFOT scenario for the glass and jamming transitions

As is by now well known, Kirkpatrick, Thirumalai and Wolynes’ scenario for the liquid-glass transition involves a first point at which the equilibrium state fractures into an exponential number of ergodic components: this is the dynamical temperature $T_{\rm d}$ (or pressure $P_{\rm d}$, see Fig. 1), also called Mode-Coupling temperature because in low dimensions it can be computed using Mode-Coupling theory. The ergodic components are only truly dynamically separated in the mean-field limit, while in a realistic short-range finite-dimensional situation the system is still ergodic, although the dynamics slows down. As the temperature is lowered, or the pressure increased, the number of metastable states contributing to equilibrium diminishes, until a point is reached where the equilibrium system is left with only the deepest amorphous states: this is the Kauzmann point, beyond which the thermodynamics stays dominated by (or “frozen in”) those states. From a purely equilibrium point of view, one may picture the situation at $P>P_{\rm K}$ (or $T<T_{\rm K}$) as in the sketch of Fig. 2, with widely separated states of “size” $q$, defined for example as:

with $a,b$ two copies (replicas) of the system and $k$ a vector of length comparable to the inverse of the inter-particle distance (alternative definitions of $q$ are possible, see Franz et al. 2013 for a review). A more precise way of stating the same thing is to introduce the effective potential $V(q)$ Franz and Parisi 1998, which counts the logarithm of the number of configurations having correlation exactly $q$ with a “reference” equilibrium configuration. One obtains a picture as in Fig. 2, where one sees that configurations are either close, with overlap $q=q_{\rm EA}$ (with probability $1-m$), or very far away – in other states, with overlap $q=0$ – with probability $m$, corresponding to the two minima in the effective potential and to the two peaks in the Parisi order parameter $P(q)$. Note that $P(q)$ may in general only be non-zero where $V(q)$ takes the minimal value, and that $q_{\rm EA}$ plays the role of the Edwards-Anderson order parameter Mézard et al. 1987.

This construction concerns equilibrium configurations, but may be generalized to describe metastable states Kirkpatrick and Thirumalai 1989; Monasson 1995 by choosing the “reference” configuration, rather than from equilibrium, from a system perturbed by a small “pinning field”, itself thermalized at a higher temperature $T^{\prime}=T/m$. Technically speaking, following Monasson Monasson 1995, this amounts to the following calculation: one considers $m$ weakly coupled replicas at temperature $T$, takes an equilibrium configurations of one of the replicas as the reference configuration, and couples to it an additional replica who is forced to stay at distance $q$ from it. Then one computes the free energy of the additional replica, and averages it over the other $m$. This amounts to fixing the Parisi parameter $m$, rather than choosing the value that maximizes the free energy. In this way, one obtains a bistable form for the effective potential, up to a threshold value $T^{\prime}_{th}=T/m_{th}$ at which the minimum close to $q_{\rm EA}$ disappears (Fig. 3), and at precisely the threshold level, the stability matrix corresponding to the minimum at $q_{\rm EA}$ develops zero modes, signaling the fact that the states close to the threshold level are marginal. This shows up within the replica scenario as the vanishing of the “replicon” eigenvalue De Dominicis and Kondor 1983, and within the Thouless-Anderson-Palmer approach Mézard et al. 1987 as the free-energy Hessian developing zero eigenvalues. A crucial result of Ref. Cugliandolo and Kurchan 1993 is that the out of equilibrium aging dynamics happens exactly at this threshold level, and it exploits the marginality of the threshold states to explore phase space.

Let us now consider higher pressures, or lower temperatures. In most systems, there is a second transition discovered by Gardner Gardner 1985; Montanari and Ricci-Tersenghi 2003; Montanari and Ricci-Tersenghi 2004 years ago, at which each state itself breaks into smaller substates. The sketch one usually makes is as in Fig. 4. To be more precise, we consider what happens with the effective potential and $P(q)$ for an equilibrium configuration beyond the Gardner point. The situation is depicted in Fig 4: there are many configurations at all distances between $q^{*}$ and $q_{\rm EA}$: the state of size $q^{*}$ has fractured into many subcomponents of smaller sizes. However, going away from a configuration, up to correlations smaller than $q^{*}$, one finds big barriers and no states, up until completely different states, having minimal overlap are reached. The Parisi parameter $m$ is now related to the probability of being in some state closer than $q^{*}$, i.e. within the “metabasin” Heuer 2008: this probability is given by $1-m$. This fracturing of a state into many smaller ones also happens at the level of metastable states Barrat et al. 1997; Montanari and Ricci-Tersenghi 2003: there is a line in the phase diagram where all states undergo a Gardner transition (Fig. 1). Metastable states may be found as above Montanari and Ricci-Tersenghi 2003, by considering $m$ as a free parameter. However, in the more complex regime beyond the Gardner transition it is not clear how to compute the threshold level that will dominate the dynamics Montanari and Ricci-Tersenghi 2003; Montanari and Ricci-Tersenghi 2004. It is possible that the threshold level could be identified by looking at the stability properties of the fluctuations at the level of $q^{*}$, but this need to be clarified. See Rizzo 2013 for some initial steps in this direction.

II.2 Vibrational modes and dynamics close to jamming

A system of hard spheres when compressed suddenly ends up in a configuration that is blocked, with the exception of a small percentage of “rattlers” that are free to move within a “cage” made by their neighbors. A mechanically stable system of hard constituents such as this may be hypostatic, hyperstatic or isostatic, depending on whether the number of contacts is less, more, or precisely just what it takes to guarantee mechanical stability. Because the system is prepared with a rapid compression, it seems unlikely that it will be hyperstatic, because if at some time during the compression it reaches stability, it is unable to move on to create further, redundant contacts. The hypothesis that is usually made is that, forgetting the rattlers, the rest of the system is precisely isostatic, the assumption being that there are no “rattling clusters” other than isolated rattling particles. For a detailed discussion of the fundamental role of isostaticity in jammed packings see Refs. Liu et al. 2011; Van Hecke 2010. An additional assumption that seems to be justified in practice is that isostaticity is “irreducible”, in the sense that there is no subset of particles that is separately isostatic: if in such a system a contact is broken, then by definition all the particles in the system eventually become mobile. Clearly, this is a critical situation. Indeed, it has been proposed in Ref. Wyart et al. 2005; Brito and Wyart 2009; Liu et al. 2011 that these packings are marginally stable from a mechanical point of view, and from this most of the anomalous scalings that are found numerically have been derived analytically. In particular, the criticality manifests itself in the spectrum of vibrations $D(\omega)$ at densities just below jamming Brito and Wyart 2006; Brito and Wyart 2007, which has a general shape as in Fig. 5, where one has to distinguish two features:

There is a branch of higher frequency modes, whose lowest frequency is $\omega^{*}$. The frequency $\omega^{*}$ goes to zero as the pressure goes to infinity.

Within the gap $0<\omega<\omega^{*}$ there are the acoustic modes, which exist even at finite pressures. Moreover, in Refs. Wyart et al. 2005; Brito and Wyart 2006; Brito and Wyart 2007; Lerner and Procaccia 2009; Xu et al. 2010; Manning and Liu 2011 it was shown that the softer modes do not look like plane waves, therefore acoustic modes are mixed with other kinds of soft modes.

In the rest of this section we will argue that the Gardner transition provides a natural explanation for the presence of soft modes at $\omega<\omega^{*}$. These modes should appear at all pressures beyond the Gardner transition. However, the connection between the soft modes observed in Refs. Wyart et al. 2005; Brito and Wyart 2006; Brito and Wyart 2007; Lerner and Procaccia 2009; Xu et al. 2010; Manning and Liu 2011 and the ones associated to the Gardner transition is not clear for the moment.

Consider the squared displacements $\widehat{\Delta}_{i}(t,t^{\prime})=|x_{i}(t)-x_{i}(t^{\prime})|^{2}$, where $i$ labels the $N$ particles of system, and its average over particles $\widehat{\Delta}(t,t^{\prime})=N^{-1}\sum_{i=1}^{N}\widehat{\Delta}_{i}(t,t^{\prime})$. In the following we will assume that the system has been prepared by some rapid compression at time $t=0$, in such a way that if $t>t^{\prime}>0$ and $t^{\prime}$ is large enough, the system is stuck into a glass state. The mean square displacement is given by the average of the squared displacement over the dynamical process,

and the variance of the squared displacement defines the so-called four-point susceptibility

These definitions can be made more precise to take into account the presence of rattlers, we refer the reader to Ref. Ikeda et al. 2013 for a detailed discussion. The “cage size” is the limit

where ‘$\infty$’ stands for times $t,t^{\prime}$ as large as the lifetime of the state. At pressure $P$, the natural scale of the displacements is $1/P$, therefore it is convenient to introduce a scaled cage size as

and the last relation is derived in Refs. Brito and Wyart 2009; Ikeda et al. 2013. For $d>2$, this quantity is finite for finite $P$, because $D(\omega)\sim\omega^{d-1}$ in the low frequency acoustic branch, however it diverges as $P\rightarrow\infty$ because $\omega^{*}$ goes to zero and the integral is dominated by $\int_{\omega^{*}}^{\infty}d\omega\;\frac{D(\omega)}{\omega^{2}}\sim\frac{D(\omega^{*})}{\omega^{*}}$ (while the integral in $0<\omega<\omega^{*}$ does not contribute to the divergence). The fluctuations of the cage size yield the four-point susceptibility (or “spin glass susceptibility”) Ikeda et al. 2013

On the one hand, from the theory of the Gardner transition, we expect $\chi_{4}(\infty)$ to diverge there, and to stay infinite up to infinite pressure. In fact, one may convince oneself that this is so just by considering the curvature of the effective potential above and below the Gardner transition, where $\frac{d^{2}V(q)}{dq^{2}}=0$. On the other hand, we may look at this from the point of view of normal modes: we split (6) in a contribution above, and one below $\omega^{*}$:

It has already been remarked in Ref. Ikeda et al. 2013 that in three dimensions even the acoustic modes will make the integral in $0<\omega<\omega^{*}$ diverge for finite $P$. However, the effect of acoustic modes shows up in $\chi_{4}(t,t^{\prime})$ only at very long time differences $t-t^{\prime}\gg 1/\omega^{*}$, so that in Ref. Ikeda et al. 2013 it was shown that the regime of $t-t^{\prime}\sim 1/\omega^{*}$ gives a good definition of the four-point susceptibility. In our large-dimensional case (actually, for all $d>4$), the density of acoustic modes is negligible, but we still expect that the first term in (7) diverges below the Gardner transition. The conclusion seems to be that there are other soft modes (below $\omega^{*}$) that do not contribute to the linear susceptibilities or to the short-time $t-t^{\prime}$ value of $\chi_{4}(t,t^{\prime})$, but dominate the limit of $\lim_{t-t^{\prime}\rightarrow\infty}\chi_{4}(t,t^{\prime})$. Although it is tempting to identify these modes with the ones observed in Refs. Wyart et al. 2005; Brito and Wyart 2006; Brito and Wyart 2007; Lerner and Procaccia 2009; Xu et al. 2010; Manning and Liu 2011, more work is needed to clarify the connection.

II.3 Out of equilibrium dynamics in the Gardner phase

The out of equilibrium dynamics of this system has not been solved, but from the structure of states one may already guess its main features. Below the Gardner transition line, the slow compaction (aging) dynamics should proceed close to the threshold level, defined as described above as the one where the stability at the level of $q^{*}$ is marginal. The relaxation process can be seen as a dynamical exploration of phase space starting from a completely correlated state ($q=1$) down to a completely decorrelated state ($q=0$). The relaxation should be fast from correlation $q=1$ down to $q_{\rm EA}$, and then proceed – in a progressively slower way as the system ages – down to $q^{*}$, and from there to zero. The fluctuation-dissipation properties may be studied by considering a system with hard spheres in a thermal bath of temperature unity, subjected to a pressure $P$ generated by either a piston or by coupling to the radii of all spheres. The response and correlation functions are as described in Refs. Parisi 1997; Berthier and Barrat 2002: the response is computed from the staggered displacement $R(t,t^{\prime})=\sum_{i}\xi_{i}\delta\langle x_{i}\rangle$ induced subjecting particles to random unit fields $\xi_{i}$ with an energy term $E_{field}=h\sum_{i}\xi_{i}\delta x_{i}$, per unit of $h$. The conjugate correlation may be taken to be the quadratic displacement $\Delta(t,t^{\prime})$ defined above. Response and correlations may be put together in a plot, which should look as in Fig. 6 for long times $t^{\prime}$, the time $t>t^{\prime}$ being used to produce a parametric plot of $R$ versus $\Delta$. Here, $\Delta^{*}$ and $\Delta_{\rm EA}$ are the values corresponding to the correlations $q^{*}$ and $q_{\rm EA}$. At every time, the first barriers encountered are the small ones close to $q_{\rm EA}$: beyond the Gardner transition, where small states are separated by relatively small excitations because there are states at all distances $q$ with $q^{*}<q<q_{\rm EA}$. This might help explain the paradox that the path between these small states is mainly along the flattest vibrational modes – as found in Ref. Berthier and Witten 2009 – while this is not what one expects for the large-scale relaxation within a supercooled liquid, at least within the RFOT scenario.

II.4 Low temperature excitations

A long standing problem in the physics of glassy and amorphous materials is the low temperature behavior of their specific heat and thermal conductivity, which turns out to be quite different to that observed in crystals. These features are all the more intriguing because they tend to be quite universal for all amorphous materials. The usual explanation for this phenomenon is to attribute it to localized quantum-mechanical two-level tunneling systems Anderson et al. 1972; Phillips 1981. These models assume that there are particles or groups of particles that evolve and tunnel in random local potentials. These potentials are usually proposed phenomenologically, although they are of course generated by the same interactions that produced the amorphous solid in the first place. Furthermore, these simple localized clusters or particles will be coupled, and their interactions might generate collective effects. A proposal to take these features into account Kuhn 2003 is to consider a system of strongly coupled localized deformations, which one may assimilate as “spin”-like excitations, and to assume that they have essentially random interactions – with long range, partly because elasticity is long range, and partly to make the system solvable. One obtains in this way a “spin-glass” of deformations, with elementary excitations which one may calculate and which tend to be universal because of their collective nature. Note that this way, phenomenology has been pushed one step up, to the effective interaction of excitations.

Quite clearly, the mechanisms that generate the coupling between low-temperature excitations, and the one responsible for the amorphous matrix on which they live, are one and the same. One would thus expect the same theory to explain both features. In the context of this paper, it is very tempting to interpret the large valleys (metabasins) of size $q^{*}$ as being the amorphous structure, and the excitations of all sizes between $q^{*}$ and $q_{\rm EA}$ as a “spin-glass” of small excitations within that amorphous structure. Formally this is clearly so, a fact that was already recognized by Gardner in her original paper, where she showed that the transition is essentially that of the Sherrington-Kirkpatrick model within each large state. More recently, the analogy between jammed packings and the Sherrington-Kirkpatrick model has also been underlined Wyart 2012, and the idea that there is a deep connection between the marginality of jammed packings and low-temperature anomalies in glasses was also proposed by S. Nagel (e.g. in his talk at the ACS meeting, Philadelphia, August 2012). A possible difference may be noted with respect to Ref. Kuhn 2003: here the spin-glass transition need not (and in general will not) coincide with the liquid-glass transition at which the amorphous matrix is formed.

III Replicated entropy of infinite-dimensional hard spheres

The above discussion provides several important motivations to look for an instability of the 1RSB solution in particle systems, akin to the Gardner transition of spin glasses Gardner 1985. Additional ones will be given by the more technical discussion that we now start, see Sec. IV.2. We will show that a Gardner instability indeed happens in hard sphere systems in the limit $d\rightarrow\infty$ (and probably also in finite dimensions within the mean-field RFOT scenario).

We will consider a system of $N$ hard $d$-dimensional spheres with unit diameter, enclosed in a volume $V$, hence at density $\rho=N/V$. The packing fraction is $\varphi=2^{-d}\rho V_{d}$, with $V_{d}=\pi^{d/2}/\Gamma(1+d/2)$ the volume of a sphere of unit radius. In the first paper of this series Kurchan et al. 2012, we derived an exact expression for the replicated entropy of this system in the limit of large dimension $d\rightarrow\infty$. The result is obtained by first writing the entropy in a manifestly rotationally and translationally invariant form, and then performing a saddle point evaluation in the limit $d\rightarrow\infty$. Within replica theory the resulting entropy is a function of the density function $\rho(\hat{q})$, where the matrix $\hat{q}$ is a $m\times m$ symmetric matrix that encodes the overlaps $q_{ab}$ between different replicas. The main result of Ref. Kurchan et al. 2012 was that a Gaussian assumption for $\rho(\hat{q})$ gives the exact result for the entropy, i.e. for all thermodynamic properties of the system, exclusively in terms of $q_{ab}$ (with $\sum_{a=1}^{m}q_{ab}=0$ for all $b$, because of translational invariance). In other words, no higher order parameters $q_{abc},q_{abcd}...$ are necessary in the large dimensional limit.

The proof of Ref. Kurchan et al. 2012 was restricted to the 1RSB form of $q_{ab}$. In this section, we will extend the results of Ref. Kurchan et al. 2012 to obtain the replicated entropy as a function of the overlap matrix, without making any assumption on the RSB structure. We will start by deriving the Gaussian replicated entropy for a generic overlap matrix, and then show that this form coincides with the exact result. Obviously, in this section we will often make reference to Ref. Kurchan et al. 2012, which we encourage the reader to consult before looking to the rest of the section. Another option is to skip this section and take the result – as expressed by Eqs. (15) and (16) – for granted. This will be the starting point to study the stability of the 1RSB solution (and much more) in the following.

We have to parametrize a generic Gaussian form of $\rho(\hat{q})$, or equivalently $\rho(\bar{u})$, where $q_{ab}=u_{a}\cdot u_{b}$ and $u_{a}$ are the $d$-dimensional vectors corresponding to replica displacements, with $\bar{u}=\{u_{1},\cdots,u_{m}\}$. We can choose a parametrization in terms of a $m\times m$ symmetric matrix $\hat{A}$ such that $\sum_{a=1}^{m}A_{ab}=0$ for all $b$. Calling $\hat{A}^{m,m}$ the $(m-1)\times(m-1)$ matrix obtained from $\hat{A}$ by removing the last line and column, the most general Gaussian form of $\rho(\bar{u})$ is

which is normalized according to $\rho=\int{\cal D}\bar{u}\rho(\bar{u})$ and ${\cal D}\bar{u}=m^{d}\delta(\sum_{a}u_{a})du_{1}\cdots du_{m}$. The parameters $A_{ab}$ are interpreted as

for $a,b\in[1,m-1]$, while $\left\langle u_{a}\cdot u_{m}\right\rangle=-\sum_{b=1}^{m-1}\left\langle u_{a}\cdot u_{b}\right\rangle=A_{am}$ and $\left\langle u_{m}\cdot u_{m}\right\rangle=\sum_{ab}^{1,m-1}\left\langle u_{a}\cdot u_{b}\right\rangle=A_{mm}$.

The saddle point value of $\hat{q}$, that dominates all the integrals over $\rho(\hat{q})$, is obtained as follows. We start from the normalization condition (note that a complete derivation of $J(\hat{q})$, that was not reported in Ref. Kurchan et al. 2012, is reported here in Appendix A)

and maximizing the exponent for $d\rightarrow\infty$ leads, for $a,b\in[1,m-1]$, to $(q^{m,m})^{-1}_{ab}=(\hat{A}^{m,m})^{-1}_{ab}/d$, hence $q^{sp}_{ab}=d\,A_{ab}$, consistently with Eq. (9).

To compute the replicated entropy using the general Gaussian ansatz we start from Eq. (45) of Ref. Kurchan et al. 2012, which gives the following expression for the replicated entropy:

where the function ${\cal F}$ is given in Eq. (37) of Ref. Kurchan et al. 2012. The ideal gas term is

In order to obtain a simple limit $d\rightarrow\infty$, it is convenient to define a matrix $\hat{\alpha}=\frac{d^{2}}{D^{2}}\hat{A}$ and a reduced packing fraction $\widehat{\varphi}=2^{d}\varphi/d$. With this choice we have

The matrix $\hat{\alpha}$ is a variational parameter and is therefore determined by maximization of the entropy. Defining ${\cal F}^{\prime}_{ab}(\hat{\upsilon})=\frac{d{\cal F}(\hat{\upsilon})}{d\upsilon_{ab}}$, the equation for $\hat{\alpha}$ is

Eqs. (15) and (16) provide the expression of the Gaussian replicated entropy and will be the starting point of all our calculations.

Note that the entropy of the equilibrium glass is obtained by optimizing $s[\hat{\alpha}]/m$, given in Eq. (15), with respect to the matrix $\hat{\alpha}$ and of $m$. Let us call $\hat{\alpha}^{*}$ and $m^{*}$ the optimal values, $\hat{\alpha}^{*}$ being the solution of Eq. (16). The reduced pressure $p=\beta P/\rho$ of the equilibrium glass is given by

This result shows that the pressure diverges whenever $m^{*}\rightarrow 0$ as $p\sim 1/m^{*}$. Hence, the density at which $m^{*}\rightarrow 0$ defines the jamming point Parisi and Zamponi 2010.

III.2 Exact computation

We now show that Eqs. (15) and (16) can be equivalently obtained by an exact evaluation of the saddle point equations derived in Eqs. (64) and (65) of Ref. Kurchan et al. 2012. In fact, we can make use of Eqs. (65) and (39) of Ref. Kurchan et al. 2012 to obtain a closed self-consistent equation for $\hat{q}^{sp}$, which as before is the point where the argument of the integral

is maximum (subleading terms for $d\rightarrow\infty$ have been neglected). Taking the derivative with respect to $\hat{q}$ and computing the result in $\hat{q}=\hat{q}^{sp}$ leads to the equation

Clearly, defining $\hat{\alpha}=\frac{d}{D^{2}}\hat{q}^{sp}$ this equation is equivalent to Eq. (16).

Evaluation of Eq. (18) at the saddle point gives the equation for $\lambda$

from which, using Eq. (78) of Ref. Kurchan et al. 2012, we obtain

Combining Eqs. (64) and (65) of Ref. Kurchan et al. 2012, using Eq. (21) and recalling the definition $\hat{\alpha}=\frac{d}{D^{2}}\hat{q}^{sp}$ we have

which coincides with Eq. (15). This completes the proof of the exactness of the Gaussian ansatz for the computation of the entropy. Note that as already observed in Ref. Kurchan et al. 2012 this does not imply that the Gaussian form (8) can be used to compute correlation functions (that encode structural properties), because the equivalence is only correct at the saddle point level for the entropy: this is consistent with the numerical observation of a non-Gaussian cage shape obtained in Ref. Charbonneau et al. 2012a. A computation of the cage shape is in progress and will be hopefully reported in future papers of this series.

IV 1RSB solution

in Ref. Kurchan et al. 2012; Parisi and Zamponi 2010; Parisi and Zamponi 2006 we studied the 1-step replica symmetry breaking (1RSB) ansatz, which amounts in this formalism to assuming that all replicas are equivalent Monasson 1995. For completeness let us recall here this result, which correponds to the simple choice

with $\widehat{A}=\frac{d^{2}}{D^{2}}A$. Within this ansatz, $A$ is the “cage radius”, as it is proportional to the long time limit of the mean square displacement in the glass Parisi and Zamponi 2010. Note that we use a small hat for matrices, while the wide hat just denotes reduced scalar variables. Hence $\hat{A}$ is a matrix while $\widehat{A}$ is a scalar, and they should not be confused.

Using the relation $\log\det\{[\widehat{A}(\delta_{ab}-1/m)]^{m,m}\}=(m-1)\log\widehat{A}-\log m$ and the results of Sec. VIIB of Ref. Kurchan et al. 2012, it is easy to check that Eq. (15) reduces to the result of Ref. Kurchan et al. 2012 for the 1RSB entropy, which is

The equation for $\widehat{A}$ is derived by optimizing the above results, leading to

The function ${\cal F}_{m}(\widehat{A})$ introduced here should not be confused with the function ${\cal F}(\hat{\alpha})$ introduced before. These are different functions, the first acts on a scalar while the second on a matrix.

The physical consequences of this expression for the entropy have been derived in Ref. Parisi and Zamponi 2006; Parisi and Zamponi 2010, where the expressions of the dynamical transition density, the Kauzmann transition density, and the GCP density have been derived, with the scalings sketched in Fig. 1. Furthermore, at the level of the 1RSB solution, we know that $\widehat{A}^{*}\sim m^{*}$ so we conclude that the cage radius vanishes as $\widehat{A}^{*}\sim 1/p$. As a consequence of this scaling, the scaled cage size $\Delta_{\infty}$ introduced in Eq. (5) is found to diverge as $\Delta_{\infty}\sim p^{2}\widehat{A}^{*}\sim p$, as noted in Ref. Ikeda et al. 2013.

IV.2 Inconsistencies of the 1RSB solution

The 1RSB predictions for physical quantities were carefully compared with numerical results, around both the glass and the jamming transitions Parisi and Zamponi 2010; Berthier et al. 2011; Charbonneau et al. 2011; Charbonneau et al. 2012a; Charbonneau et al. 2012b. Despite the good overall agreement with numerical data, one expects, as described above, that a Gardner transition to a full replica symmetry breaking scheme is generic. Furthermore, several inconsistencies have been found close to the jamming transition, at very high pressure:

The 1RSB solution predicts the existence of jammed packings with density $\varphi_{j}$ in the interval $\frac{\varphi_{j}}{2^{-d}d}=\widehat{\varphi}_{j}\in[6.26,\log d]$ Parisi and Zamponi 2010. However, only the packings with $\widehat{\varphi}_{j}\sim\log d$ are isostatic, with each particle in contact, on average, with $z=2d$ other particles. Instead, the packings with $\widehat{\varphi}_{j}$ of order 1 are found to be hyperstatic with $z>2d$, which is, as mentioned above, unexpected and inconsistent with numerical results.

In the glass phase, the exact relation between the pressure $p$ and the contact value $y(\varphi)$ of the pair correlation, $p=1+2^{d-1}\varphi y(\varphi)$ Hansen and McDonald 1986, is violated. In particular, it is found that when $\varphi\rightarrow\varphi_{j}$, $p\sim d\,\varphi_{j}/(\varphi_{j}-\varphi)$, consistently with numerical results, while

where the first term is the one that is consistent with the scaling of the pressure. Hence, the correct relation between $p$ and $y(\varphi)$ is recovered only if $2^{1-d}\,d/\varphi_{j}=2/\widehat{\varphi}_{j}\ll 1$ when $d\rightarrow\infty$, which again suggests that the 1RSB solution is inconsistent when $\widehat{\varphi}$ is of order 1 and might be stable only when $\widehat{\varphi}\gg 1$.

The scaling at large (reduced) pressure $p$ of the cage radius $A$ (the long time limit of the mean square displacement in the glass) predicted by the 1RSB solution is $A\sim p^{-1}$, while the marginal stability argument of Ref. Wyart et al. 2005; Brito and Wyart 2009 predicts that $A\sim p^{-3/2}$, which has been confirmed numerically in several studies, e.g. Refs. Wyart et al. 2005; Brito and Wyart 2009; Ikeda et al. 2013. This exponent controls all the other exponents that characterize criticality at the jamming transition Ikeda et al. 2013 and is directly related to the anomalous soft vibrational modes that appear at jamming Liu et al. 2011; Van Hecke 2010; Brito and Wyart 2009; Ikeda et al. 2013, hence reconciling the theoretical prediction with the numerical results is of extreme importance.

Other exponents that characterize the structure at jamming, for instance the famous (almost) square-root singularity in the pair correlation function O’Hern et al. 2002; Donev et al. 2005; Silbert et al. 2006; Charbonneau et al. 2012b, are not reproduced by the 1RSB solution, at least at the Gaussian level (a more detailed calculation of the structure functions based on the non-Gaussian theory developed in this series of papers is in progress and will hopefully be reported in a future paper).

All these considerations suggest strongly that the 1RSB solution is unstable, at least when pressure is large enough and $\widehat{\varphi}$ is of order $1$. They provide further motivations to study the stability of the 1RSB solution Gardner 1985, which is the subject of the next section.

V Second order expansion: the Hessian matrix

Here we obtain the expansion around the 1RSB solution at the quadratic order. For this we need to consider a more general ansatz or expand around the 1RSB solution, which is made possible by the general expression of the entropy obtained in Eqs. (15) and (16). The quadratic expansion of the entropy around the 1RSB solution provides a stability matrix whose eigenvalues allow one to determine the stability of the solution. We will find that, as it happens in a Gardner transition Gardner 1985; Montanari and Ricci-Tersenghi 2003; Montanari and Ricci-Tersenghi 2004, the 1RSB solution becomes unstable when pressure is large enough. In finite and arbitrarily large dimensions, this happens for all $\widehat{\varphi}_{j}$. However, for $d\rightarrow\infty$, the so-called glass close packing (GCP) Parisi and Zamponi 2010 which is the densest amorphous packing and has $\widehat{\varphi}\sim\log d$ becomes stable again, suggesting that the 1RSB predictions for jamming are still approximately useful as a starting point, but should be corrected to take into account its instability.

Because the matrix $\hat{\alpha}$ should have the sum of the elements of every column and every row equal to zero, we can say that the independent entries are the elements above the diagonal of the matrix, provided that the matrix is symmetric and the diagonal is fixed in such a way that the constraints on the sum over the elements in a row or in column is satisfied. Hence, in the following we denote as $\delta/\delta\alpha_{a<b}$ the derivative taken with respect of the element $\alpha_{ab}$ with $a<b$ which is assumed to be the only independent element (hence its variation induces a variation of $\alpha_{ba}$ and of the diagonal elements $\alpha_{aa}$ and $\alpha_{bb}$). Taking into account all this, we define the Hessian matrix as

where the replica structure is a consequence of the structure of $\hat{\alpha}^{\rm 1RSB}$. Although the matrix $M$ is defined by the above equation only for $a<b$ and $c<d$, we will define it for convenience also for $a>b$ and $c>d$ assuming that it is symmetric (hence the notation $M_{a\neq b,c\neq d}$). The prefactor $2\widehat{A}^{2}/d$ is chosen for later convenience and is positive, hence it does not affect the sign of the three different eigenvalues of the mass matrix, which are

Defining the “entropic” and “interaction” terms

We now compute these two terms separately.

V.2 The entropic term

To compute the entropic term it is convenient to introduce a shorthand notation $\hat{\beta}=\hat{\alpha}^{m,m}$. Let us also use indices $i,j,k,\cdots$ for $\beta$ to highlight that they run from $1$ to $m-1$. In the 1RSB solution $\beta^{\rm 1RSB}_{ij}=\widehat{A}(\delta_{ij}-1/m)$ has the same form of $\hat{\alpha}$ but on the reduced $(m-1)\times(m-1)$ space. Then we have

From the definition $\hat{\beta}\hat{\beta}^{-1}=I$ where $I$ is the identity matrix we have

V.3 The interaction term

We now consider the interaction term. Hence we need an expansion of the function ${\cal F}(\hat{\upsilon})$ around the 1RSB solution. Recall that $\hat{\upsilon}$ is a $m\times m$ symmetric matrix such that the sum of the elements in each row and in each column is zero. Starting from the results of Section V of Ref. Kurchan et al. 2012, we can write explicitly the function ${\cal F}(\hat{\upsilon})$, introducing $m$-dimensional vectors $x_{a}$ such that $x_{a}\cdot x_{b}=\upsilon_{ab}$, as

First let us compute ${\cal F}$ on the 1RSB solution where

Let us call $\delta_{a,{\rm min}}$ a function that is equal to one only if $a$ is such that $\lambda_{a}$ is the minimum among all the $\{\lambda_{a}\}$: or in other words $\min_{a}\lambda_{a}=\lambda_{a}\delta_{a,{\rm min}}$. Then we have

and so on. Eq. (40) provides the derivation of the interaction part of Eq. (24) Parisi and Zamponi 2010; Kurchan et al. 2012.

V.3.2 Monomials of nn

We now want to expand the quantity $\mathcal{F}$ around the 1RSB solution. The part of the Hessian matrix coming from the interaction term is:

where the functions $f$ are at most quadratic:

we have that the Hessian matrix is given by

Hence we want to compute averages of monomials of the $\{n_{a}\}$, which can be written as follows:

where the definition of the average has been changed to

The factor in parenthesis is a polynomial in $\{\lambda_{a}\}$, hence we now want to be able to write averages of monomials of $\lambda$. Using the replica symmetry of the average over $\{\lambda_{a}\}$, we need in particular the following objects:

V.3.3 Monomials of λ\lambda

We therefore need to compute several monomials of the $\{\lambda_{a}\}$, which are listed in the following. Calculations follow Eq. (40) and it will be convienient to define one more average over $\lambda$ as

V.3.4 The structure of the mass matrix

Thanks to replica symmetry the Hessian matrix has only three independent matrix elements. These are

It is however convenient to write the matrix in this form:

The above equation, together with Eqs. (50) and (53), give the complete expression of the interaction part of the Hessian matrix.

V.4 The replicon

The stability of the 1RSB solution depends crucially on the replicon eigenvalue, $\lambda_{R}=M_{1}=\widehat{A}^{2}M_{1}^{(E)}-4\widehat{\varphi}\,\widehat{A}^{2}M^{(I)}_{1}$. Collecting all the above results and simplifying some terms we get

We can compute this numerically on the 1RSB solution, where $\widehat{A}$ is the solution of Eq. (25), to find the point where $\lambda_{R}=0$ and the 1RSB solution becomes unstable. The instability curve in the $(m,\widehat{\varphi})$ plane, which corresponds to the line $\widehat{\varphi}_{G}(m)$ on which the replicon vanishes, is reported in Fig. 7.

Asymptotically we obtain $\widehat{\varphi}_{G}(m)\sim m^{-1/2}$ for $m\rightarrow 0$. To explain this we must investigate the asymptotics of the function $\Lambda(m,\widehat{A})$ when both $m$ and $\widehat{A}$ are small. It is convenient to use Eq. (25) to eliminate $\widehat{\varphi}$ instead of $\widehat{A}$. Doing this the equation for the instability becomes

which must be solved to obtain $\widehat{A}_{G}(m)$ and then $\widehat{\varphi}_{G}(m)$ using Eq. (25). We want to show that $\widehat{A}_{G}(m)\sim m^{2}$, and that this implies $\widehat{\varphi}_{G}(m)\sim m^{-1/2}$.

First of all let us examine the asymptotics of the different terms for large and positive $\lambda$. For $\lambda\rightarrow\infty$ we have

It will be convenient for the following to define

Now we expand Eq. (58) at small $\widehat{A}$. Let us recall that ${\cal G}_{m}(\widehat{A})=1-m\left\langle\Theta_{0}(\lambda)^{m-1}\right\rangle$. From this we obtain

and similarly (the fact that the horrible integral corresponding to $\Lambda_{m}(\widehat{A}=0)$ is exactly 0 can be proven by a series of integrations by parts):

Asymptotically for small $m$, the integrals in the above expressions can have different behaviors, depending on the behavior of the integrand for large $\lambda$ when $m\rightarrow 0$. In fact, if the integrand decays faster than $1/\lambda$, the integral is well defined and has a finite limit for $m\rightarrow 0$. In the opposite case, the integral is divergent and the divergence is dominated by the large $\lambda$ behavior: in this case one has to analyze the possibly divergent part to determine the behavior of the integral at $m\rightarrow 0$. In the case of $\Delta_{1}(m)$, thanks to a subtle cancellation, the large $\lambda$ contribution to the integral is

The first two terms give contributions that are not divergent when $m\rightarrow 0$, hence they are subleading with respect to the $1/\lambda^{4}$ term that gives a finite contribution. We conclude that $\Delta_{1}(m)$ has a finite limit given by

Instead, the leading large $\lambda$ behavior of the integrand of $\Delta_{2}(m)$ is, changing variable to $y=\sqrt{m}\lambda$:

We conclude that $\sqrt{\widehat{A}_{G}}\approx 0.8m$. Finally, we can show similarly that for small $m$

Both asymptotic results for $\widehat{A}_{G}$ and $\widehat{\varphi}_{G}$ are perfectly consistent with the numerical data.

V.5 The 2RSB solution

When the 1RSB solution becomes unstable, one must consider further RSB. We performed a 2RSB calculation. Then we can linearize the 2RSB solution close to the 1RSB one and obtain the line at which the 2RSB provides a better maximization of the entropy, hence becoming stable. This provides an independent calculation of the 1RSB instability, which we verified to be coherent with the one reported above. A complete characterization of the 2RSB solution (as well as the 3, 4, $\cdots$, $\infty$RSB ones) will be presented in future papers of this series.

VI Computation of the dynamic exponents from the cubic expansion

The same strategy allows to obtain the cubic terms in the expansion. From these, following the procedure of Ref. Caltagirone et al. 2012; Franz et al. 2013, one can compute the mean-field dynamical critical exponents at the dynamical glass transition, the so-called exponent parameter of Mode-Coupling theory, $\lambda_{\rm MCT}$. Although this calculation is not the main scope of this paper, we report it in this section.

Let us define, following the same notation as for the second order terms (hence for $a\neq b$, $c\neq d$, $e\neq f$ which we omit from now on)

Exploiting the replica symmetry, the two coefficients $w_{1}$ and $w_{2}$ can be written in the following form

and we then have $\lambda_{\rm MCT}=\frac{w_{2}}{w_{1}}$.

hence we have a similar relation for $w_{1}$ and $w_{2}$. We now compute these two terms separately.

Following the same strategy as in Sec. V.2, we obtain

Using Eq. (71), the results of Sec. V.2 and performing the traces we obtain

VI.2 The interaction term

Using the expression of $f$, expanding the products, and simplifying many monomials using the symmetries (e.g. $\left\langle n_{a}^{3}n_{d}n_{f}^{2}\right\rangle=\left\langle n_{a}^{3}n_{b}^{2}n_{c}\right\rangle$), we obtain

Now we convert the average over $n$ into an average over $\lambda$, using Eq. (48). Performing the derivatives and exploiting similar symmetries to simplify the result we obtain

This result, together with Eq. (53), allows for the explicit computation of these terms. After some simplifications, we obtain

which can be easily computed numerically.

VI.3 Numerical result

Collecting all the terms together we obtain the final result

When computed at the dynamical transition with $m=1$, $\widehat{\varphi}=\widehat{\varphi}_{d}=4.80677$ and $\widehat{A}=\widehat{A}_{d}=0.57668$ given by the solution of Eq. (25), we obtain

which implies that the MCT exponents are $a=0.324016$, $b=0.629148$ and $\gamma=2.33786$. The result for $\gamma$ is roughly consistent with the numerical estimate of Ref. Charbonneau et al. 2012a.

VII Phenomenological extension to finite dimensions

We can obtain quantitative results in finite $d$ by a phenomenological extension of Eq. (15). First we go back to non-rescaled density and we rearrange it as

We now recognize that $s_{liq}=1-\log\rho-2^{d-1}\varphi$. Furthermore, by comparison with the finite $d$ results obtained in the small cage expansion Parisi and Zamponi 2010, we know that the interaction term is renormalized by the contact value of the liquid correlation $y_{liq}(\varphi)$. We therefore can propose the following form for the entropy:

Although these equations are not obtained from a consistent derivation, recalling that for small $\widehat{A}$ we have ${\cal G}_{m}(\widehat{A})\sim\sqrt{\widehat{A}}G_{1}(m)$ and $G_{1}(m)=2Q_{0}(m)$, they reproduce the small cage expansion of Ref. Parisi and Zamponi 2010 at the leading order in $\widehat{A}$. Note that when expressed in terms of $\widehat{A}$ and $m$, the equation for the stability $\lambda_{R}=0$ is exactly the same as in $d\rightarrow\infty$, Eq. (58). Hence the result for $\widehat{A}_{G}(m)$ is independent of dimension.

We will check a posteriori that even in $d=3$ the Gardner transition happens at very large pressure, hence $m$ and $\widehat{A}$ are small. So we can use the asymptotic expansions to obtain quantitative estimates. The procedure is the following:

Recall that at small $m$ we have $\sqrt{\widehat{A}_{G}}\approx 0.8m$.

Now we obtain $\varphi_{G}(m)$ (or better $m_{G}(\varphi)$) by solving

We recall from the analysis of Ref. Parisi and Zamponi 2010 that

The Gardner transition happens when the two lines cross, hence $\varphi_{G}$ is the solution of $m_{G}(\varphi)=m^{*}(\varphi)$ which can be easily found numerically once an equation of state for the liquid has been chosen. Here we use the Carnahan-Starling equation already used in Ref. Parisi and Zamponi 2010.

Finally we use the result Parisi and Zamponi 2010

to estimate the Gardner pressure $p_{G}=p(\varphi_{G})$.

The numerical values of the Gardner pressure are reported in Tab. 1. Note that the non-monotonicity of $p_{G}$ at low dimension could be an artifact of the approximations used above. Note also that $p_{G}$ is always much larger than the pressure at the glass transition (reported in Ref. Parisi and Zamponi 2010; Charbonneau et al. 2011). Still, the reader should keep in mind that the value of $p_{G}$ corresponds to the Gardner instability of the equilibrium (“ideal”) glass. According to the phase diagram of Fig. 1, we expect that the Gardner instability for the metastable states that are reached out of equilibrium will happen at lower pressures. Unfortunately, quantifying this effect requires the use of “state following” techniques Krzakala and Zdeborová 2010 and goes far beyond the scope of this article.

The limit $d\rightarrow\infty$ is recovered as follows. Recall that $\widehat{\varphi}_{GCP}\sim\log d$ Parisi and Zamponi 2010. Moreover, $y_{liq}\rightarrow 1$. Hence $\mu\sim 2^{d-1}/d$ and the equation for $\widehat{\varphi}_{G}$ becomes

which shows that the distance between $\widehat{\varphi}_{G}$ and $\widehat{\varphi}_{GCP}$ shrinks as $(\log d)^{-2}$ and the Gardner pressure diverges as $p_{G}\sim d(\log d)^{3}$, as was sketched in Fig. 1.

Finally, at this level of approximation, it can be easily shown that $\lambda_{\rm MCT}$ does not depend on dimension. The small dependence of $\lambda_{\rm MCT}$ reported in Ref. Charbonneau et al. 2012a should be explained by corrections to this approximation, that have been neglected here.

VIII Conclusions

In this paper we were able to investigate the possibility of a Gardner transition for hard spheres in large spatial dimensions. Such a study was never been done before for particle systems, and was possible thanks to the expression of the entropy in terms of the overlap matrix obtained in the first paper of this series, and extended here to obtain Eqs. (15) and (16): because this expression has been shown to be exact, any discrepancy or instability is only attributable to the instability of the 1RSB ansatz.

The 1RSB solution is unstable in the equilibrium glass phase for (reduced) pressure higher than the Garner pressure $p_{G}$ Gardner 1985, and for metastable states in large region of the phase diagram, just as in $p$-spin glasses Montanari and Ricci-Tersenghi 2003; Montanari and Ricci-Tersenghi 2004. We provided an estimate of the Gardner pressure in finite dimensions, finding that it is quite high; moreover we showed that $p_{G}$ diverges (slowly) with increasing dimension. These pressures are directly accessible to numerical simulations, hence we expect that this transition should be quite easy to detect numerically. Estimating analytically the transition pressure for metastable glasses would be very useful to guide numerical simulations: this is however hard, as it requires a state following computation Krzakala and Zdeborová 2010. Although this is possible in principle, we leave it for future work. We expect that in any case the instability will happen at lower pressure for metastable glasses than for the ideal glass.

The physical consequences of this instability are very intriguing but for the moment not all its implications have been worked out. In fact, even for the simplest $p$-spin glasses, the impact of the Gardner instability on the out-of-equilibrium dynamics is not completely understood from a technical point of view Montanari and Ricci-Tersenghi 2004. The structure of the metastable states of the $p$-spin glass model and its impact on the out-of-equilibrium dynamics are being actively investigated Gardner 1985; Montanari and Ricci-Tersenghi 2003; Montanari and Ricci-Tersenghi 2004; Crisanti et al. 2005; Tonosaki et al. 2007; Krzakala and Zdeborová 2010; Zdeborová and Krzakala 2010; Rizzo 2013 and making progress on this simpler model will be crucial for understanding the technically more involved hard sphere case. We expect (hope) that the scenario we proposed in Sec. II will be confirmed by these studies.

Let us recall here some speculations on the possible impact of the Gardner instability on the physics of jamming that we discussed in this paper, leaving a more detailed investigation for future works.

It is reasonable to expect that at the Gardner transition, the 1RSB solution will transform continuously into a full RSB solution, although we have not yet constructed this explicitly. Such a solution describes a situation where glassy states are arranged in a complex and correlated pattern Mézard et al. 1987. More importantly, they are marginally stable Bray and Moore 1979; Mézard et al. 1987. This means that the spectrum of vibrations around a glassy state displays many soft modes. Hence, it is likely that a full RSB description of the problem will allow one to obtain information on the soft modes that are observed at the jamming transition Wyart 2005; Liu et al. 2011; Van Hecke 2010, especially those of frequency below the gap $\omega^{*}$. Some steps in this direction have been already performed in Ref. Wyart 2012, where the analogy with the full RSB physics of the Sherrington-Kirkpatrick model was noticed.

There should be several signatures of the Gardner transition. Suppose that the hard sphere system is prepared in a glass state in the region where the 1RSB solution is stable, and that pressure in slowly increased approaching the instability. As mentioned in Sec. II, the spin glass susceptibility Mézard et al. 1987 (which in this context is a four-point static susceptibility) diverges on approaching the instability. Moreover, even if the system was already equilibrated in the initial glass state, aging effects should appear below the instability when the state breaks down into many correlated sub-states.

The aging curves, and in particular the fluctuation-dissipation plots, should give a good indication of the transition: at pressures above the Gardner pressure these plots should crossover from two straight lines to two straight lines joined by a curved segment (although detecting the curved part might be numerically challenging).

There is probably a relation between the Gardner transition and the “dynamic criticality” defined in Ref. Ikeda et al. 2013. As noticed there, all the anomalous scalings at the jamming transition are related to the scaling of the cage radius with pressure, $A\sim p^{-3/2}$. Hopefully, this scaling, which is not found in the 1RSB solution, could be a property of the full RSB phase. In fact, it is well known that in the Sherrington-Kirkpatrick the presence of a full RSB phase changes the scaling of the overlap at low temperatures.

Brito and Wyart Brito and Wyart 2009 have demonstrated that close to jamming, the dynamics is characterized by sudden “cracks” at which the system leaves abruptly a locally stable structure to find a new one. At these cracks, the displacements of the particles are strongly correlated with the lowest frequency eigenvectors of the stability matrix of the structure that the system is leaving. As mentioned in the introduction, this is not what is expected in a 1RSB phase, where states are locally stable and one has to cross a barrier to jump from one state to the other – hence the vibrations at the bottom of the well should give no information on the shape of the barrier. However, in a full RSB phase the dynamics is much different and similar to the one found in Ref. Brito and Wyart 2009. It would be very nice to check whether the results of Brito and Wyart really fit into a full RSB picture, for example by calculating the spatial distribution of the displacements between two nearby states.

Finally, the response of full RSB magnetic systems to an external perturbation is very complex, being characterized by avalanches and intermittency, see e.g. Ref. Le Doussal et al. 2012. This is due to the existence of (relatively low) barriers separating nearby states – all this within a large basin. By analogy, we would expect the response of a hard sphere system to a mechanical perturbation in the full RSB phase to be similarly complex. Hence, the rheological properties in this phase could be very interesting and could explain some of the anomalous behavior found around the jamming transition. Analytical computations might be possible following the strategy introduced in Ref. Yoshino and Mézard 2010; Yoshino 2012a; Yoshino 2012b.

From the technical point of view, the next step to make progress is to investigate the $K$RSB solutions, with $K=2,3,4,\cdots$, eventually with $K\rightarrow\infty$ that corresponds to full RSB. Following Gardner’s example in the $p$-spin model, this can be done just below the Gardner transition. This investigation is in progress and will provide some answers to the above questions. In parallel, numerical simulations should be performed to detect the 1RSB instability. Also, an exact solution of the dynamics, along the lines of Ref. Mari and Kurchan 2011, could provide very useful complementary informations.

To conclude, let us mention that in this paper, from the expansion of the cubic terms around the 1RSB solution (see Section VI), we obtained an estimate of the mean-field dynamical critical exponents at the dynamical transition (the so-called Mode-Coupling Theory exponent parameter $\lambda_{\rm MCT}$) Götze 2009. We found that for $d\rightarrow\infty$, $\lambda_{\rm MCT}=0.70698$, which is consistent with numerical simulations Charbonneau et al. 2012a. This is important because a previously attempted calculation from the replicated HNC equations Franz et al. 2013 gives results that are quantitatively bad. The fact that in $d\rightarrow\infty$ we can obtain a good result implies that the negative result obtained in Ref. Franz et al. 2013 has to be attributed to the poor quantitative performances of the replicated HNC approximation, which indeed were already known Parisi and Zamponi 2010. Unfortunately, also the approach presented here gives poor quantitative results for the dynamical glass transition in low dimensions Parisi and Zamponi 2010. Obtaining an accurate theory of the dynamical glass transition in low dimension by improving the replicated HNC is therefore very important, see Jacquin and Zamponi 2013 for a preliminary step in this direction.

Appendix A Computation of the Jacobian J⁡(q^)J(\hat{q})

The Jacobian $J(\hat{q})$ is defined in the following way

Let us introduce the following notation. We define a $d\times(m-1)$ matrix $U$ whose first column contains the components of the $d$-dimensional vector $u_{1}$, the second column contain the components of $u_{2}$ and the $m-1$ column contains the components of $u_{m-1}$. Moreover we define the $m\times m$ matrix $\hat{q}$ and its reduced version that is the matrix $\hat{q}^{(m,m)}$ which is obtained from $\hat{q}$ by deleting the last column and the last row. Let us consider the integral in the expression (89):

where the last integral is with a flat measure over the entries of the matrix $U$. To perform this computation we start from a simple case. Let us consider the case where the matrix $\hat{q}^{(m,m)}$ has a diagonal structure

The second Dirac delta function says that the vectors $u$ have to be all orthogonal one to the other. The first one tells us about the length of these vectors. By going to polar coordinates we have that the previous expression is given by

The integral that appears in the expression for $C_{m,d}$ is very simple. In fact the Dirac deltas tell us that the unit vectors $\hat{u}$ must be all orthogonal. So we have to compute the phase space accessible to them. This can be done iteratively. Suppose that we have only one unit free vector. It has a phase space available given by the solid angle in $d$ dimensions which is $\Omega_{d}$. Then we add a second unit vector orthogonal to the first one. Clearly it has a phase space available that is the solid angle in the space orthogonal to the first vector that is $\Omega_{d-1}$. Going on by iteration we have

Now we want to generalize this to a matrix $\hat{q}^{(m,m)}$ that is not diagonal. Because the matrix $\hat{q}^{(m,m)}$ is symmetric we can diagonalize it and we can write

that has a Jacobian which is unitary due to the orthogonality property of the matrix $\Lambda$, then the previous expression becomes

that is the same that we discussed above. It follows that

which is the result that was used in Ref. Kurchan et al. 2012.

References