Exact theory of dense amorphous hard spheres in high dimension. III. The full RSB solution

I Introduction

In the two previous papers of this series Kurchan et al. 2012; Kurchan et al. 2013 (on which the present work relies heavily), we have considered a system of identical hard spheres in spatial dimension $d\rightarrow\infty$, and we have obtained an exact expression of its replicated partition function. Within the Random First Order Transition (RFOT) general scenario for the glass transition Kirkpatrick and Wolynes 1987a; Kirkpatrick and Thirumalai 1988; Kirkpatrick et al. 1989; Wolynes and Lubchenko 2012, which relies on the conjecture that structural glasses behave in the same way as a certain class of mean field spin glass models, the replicated partition function describes the amorphous arrested states of the system, i.e. its glasses Kirkpatrick and Thirumalai 1989; Monasson 1995.

Previous work on structural glasses in the RFOT context always focused on the simplest replica scheme, the one-step replica symmetry breaking (1RSB). The 1RSB scheme indeed already predicts the most important phase transitions that happen upon approaching the glass phase, namely the dynamical and Kauzmann transitions. An approximate 1RSB treatment of finite-dimensional glasses was developed in Mézard and Parisi 1999; Mezard and Parisi 2012, and applied to hard spheres in Parisi and Zamponi 2010; these works have shown that the 1RSB approach gives quite accurate predictions of thermodynamic and structural quantities in three-dimensional glasses. Moreover, within the 1RSB scheme, the exact expression obtained in Kurchan et al. 2012; Kurchan et al. 2013 coincides with the infinite-dimensional limit of the finite-dimensional approach of Parisi and Zamponi 2010 while a series of numerical works Charbonneau et al. 2011; Charbonneau et al. 2012a in $d=3,\cdots,12$ have shown that the main qualitative properties of the system evolve smoothly with dimensions (although the investigated dimensions are quite far from the asymptotic $d\rightarrow\infty$ regime Charbonneau et al. 2011). One could therefore think that the 1RSB scheme is sufficient to describe the properties of glasses within RFOT theory in all the range of physical parameters.

However, when one tries to apply the 1RSB scheme to study the properties of hard or soft-sphere glasses at very high pressures and low temperatures, close to the jamming transition that marks the transition from a mechanically loose to a mechanically rigid glass state, one finds contradictory results. On the one hand, the behavior of the main thermodynamic quantities (pressure, energy, …) is quite well reproduced Parisi and Zamponi 2010; Berthier et al. 2011. On the other hand, the scaling properties of other observables are not. It is therefore natural to search for the origin of this discrepancy.

A first step in this direction was made through the investigation of nearly jammed sphere packings. It was shown that hard sphere glasses at high pressures are close to a mechanical instability, the so-called “isostatic” point where the number of mechanical constraints exactly equals the number of degrees of freedom Moukarzel 1998; Tkachenko and Witten 1999; Roux 2000. Due to this proximity, anomalous low-frequency modes appear in the vibrational spectrum O’Hern et al. 2002; O’Hern et al. 2003; Wyart et al. 2005; Brito and Wyart 2006; Brito and Wyart 2007. Based on this observation, in a series of papers Wyart and coworkers Wyart et al. 2005; Brito and Wyart 2009; Liu et al. 2011; Wyart 2012 assumed that hard sphere glasses at high pressure are marginally stable, and under this assumption they derived a scaling theory of the jamming transition that is able to describe most of its basic phenomenology. In particular, it was shown both analytically and numerically Wyart et al. 2005; Brito and Wyart 2009; Ikeda et al. 2013 that marginality is associated with a particular scaling of the mean square displacement $\Delta$ in the hard sphere glass, which vanishes when the pressure $p\rightarrow\infty$ as $\Delta\sim p^{-\kappa}$ with $\kappa\sim 3/2$ (naive free-volume considerations would suggest $\Delta\sim p^{-2}$). It was confirmed in Ikeda et al. 2013 that the exponent $\kappa$ plays an important role, and in fact controls all the criticality of the jamming transition.

This analysis suggests that a 1RSB description is incorrect in this regime, because 1RSB states are perfectly stable and do not show any sign of marginality. In fact, the 1RSB solution further wrongly Still, as noted in Ikeda et al. 2013, the 1RSB prediction shows a sign of an anomalous behavior of the mean square displacement with respect to the naive expectation. predicts $\Delta\sim p^{-1}$ Parisi and Zamponi 2010, and furthermore misses other critical exponents associated to the structure and the force distribution Charbonneau et al. 2012b; Wyart 2012; Lerner et al. 2013. In the context of spin glasses, it is well known that fullRSB phases are always associated with marginal stability and anomalous low-frequency modes Bray and Moore 1979; Mézard et al. 1987, which affect the low-temperature scaling of physical quantities Mézard et al. 1987. It appears therefore that a fullRSB phase is a natural candidate for explaining the marginal stability of low-temperature glasses.

Additional insights came from analyzing the out-of-equilibrium dynamical behavior of one of the simplest spin glass models that are at the basis of the RFOT scenario, the Ising $p$-spin model Montanari and Ricci-Tersenghi 2003; Montanari and Ricci-Tersenghi 2004; Rizzo 2013; Krzakala and Zdeborová 2013. It was shown that, although the 1RSB solution correctly describes the approach to the glass phase at equilibrium, when the system is instantaneously quenched out-of-equilibrium it evolves towards a region of phase space where the 1RSB solution is unstable, the so-called Gardner phase The existence of a transition at low temperature from a 1RSB to a fullRSB solution was discovered independently in Gardner 1985; Gross et al. 1985 and its precise location was found in Gardner 1985. Gardner 1985; Gross et al. 1985, in which a fullRSB solution is present Rizzo 2013. The behavior of the Ising $p$-spin model is not yet fully understood, but this model nonetheless strongly suggests that fullRSB effects may be observable in low-temperature glasses.

Many other observations hint to the presence of a fullRSB phase in low temperature or high pressure glasses (see e.g. Fullerton and Moore 2013), which we reviewed in Kurchan et al. 2013. Based on these observations, in the previous paper of this series Kurchan et al. 2013 we used our exact expression for the replicated partition function of infinite-dimensional hard spheres to investigate the stability of the 1RSB phase, a computation that could not be done previously because the finite-dimensional approximations of Mézard and Parisi 1999; Mezard and Parisi 2012; Parisi and Zamponi 2010 have been formulated only at the 1RSB level. Our analysis showed that the 1RSB phase is stable around the dynamical glass transition, but becomes unstable at high pressures close to the jamming transition. This result gives additional indications in favor of the presence of a different phase that could better describe the properties of high pressure (or low temperature) glasses.

The aim of this work is to write explicitly the replica equations in the $k$RSB scheme and send $k\rightarrow\infty$ to describe the fullRSB solution, in order to check if this solution predicts correctly the scaling properties of the jamming transition. In the first part of this work we derive the set of equations that give the replicated entropy at the level of $k$RSB. Interestingly, the equations we obtain are remarkably similar to those describing the Ising $p$-spin model with, however, some crucial differences. The proof of this formal similarity between the equations that describe a disordered mean field spin glass model and a system of interacting particles without any quenched disorder somehow completes the program initiated by Kirkpatrick, Thirumalai and Wolynes, who constructed the RFOT scenario by assuming that such a similarity existed and proved some arguments supporting it, see e.g. Kirkpatrick and Wolynes 1987b; Kirkpatrick and Thirumalai 1989.

In the second part of the paper, we extract physical quantities from the equations. Our main results are that (i) a fullRSB phase always exists at high pressures, hence the jamming transition always lies in the fullRSB region; (ii) as in the SK model Thouless et al. 1980; De Dominicis and Kondor 1983; Goltsev 1983; Kondor and de Dominicis 1986, the fullRSB phase is marginally stable, because one of the eigenvalues of the stability matrix of the replicated entropy is identically vanishing; (iii) jammed packings are predicted to be isostatic at all densities; (iv) the fullRSB solution predicts a different scaling of the mean square displacement in the glass, namely $\Delta\sim p^{-\kappa}$ with an exponent $\kappa$ very close to $3/2$; (v) the fullRSB solution predicts a power-law divergence of the pair correlation function of jammed packings at contact, $g(r)\sim(D-r)^{-\alpha}$; (vi) the force distribution vanishes at small forces as $P(f)\sim f^{\theta}$; (vii) we obtain analytical predictions for the exponents $\kappa,\alpha,\theta$.

In Wyart 2012; Lerner et al. 2013; DeGiuli et al. 2014, the exponents $\theta$ and $\alpha$ were argued to control the stability of packings and the presence of avalanches, and a scaling relation between them was derived assuming marginal stability. The value of one exponent was lacking however to have a complete scaling picture. Our predicted values of $\kappa,\alpha,\theta$ are perfectly compatible with the scaling relations derived in Wyart 2012; Lerner et al. 2013; DeGiuli et al. 2014 and with previous numerical investigations Donev et al. 2005; Silbert et al. 2006; Torquato and Stillinger 2010; Wyart 2012; Charbonneau et al. 2012b; Lerner et al. 2013 (except for $\theta$, where the situation remains somehow unclear). By means of additional numerical simulations, we show that the prediction for $\kappa$ is correct in all spatial dimensions.

These results open the way to many different studies. For example, one could now hope to compute the shear modulus Brito and Wyart 2006; Yoshino and Mézard 2010; Yoshino 2012; Yoshino 2013, the distribution of avalanche sizes Le Doussal et al. 2010, the complete scaling on both sides of the jamming transition Ikeda et al. 2013; Berthier et al. 2011, and so on. Moreover, the results could be extended to finite dimensional systems within the effective potential approximation scheme of Parisi and Zamponi 2010; Berthier et al. 2011, which would allow one to directly compare the predictions with numerical simulations and experiments. We discuss briefly these possible developments in the conclusions. An important open question is how the fullRSB structure affects the off-equilibrium dynamics Montanari and Ricci-Tersenghi 2003; Montanari and Ricci-Tersenghi 2004; Rizzo 2013.

Because the paper is quite long, we do not present a detailed plan here. At the beginning of each section, we explain what is the aim of the section and detail the structure of subsections. Reading the technical part of the paper requires familiarity with the general concepts of spin glass theory (see e.g. Mézard et al. 1987; Castellani and Cavagna 2005 for reviews), with its application to structural glasses Mezard and Parisi 2012 and jamming Parisi and Zamponi 2010, and with the previous papers of this series Kurchan et al. 2012; Kurchan et al. 2013. A short account of this work, where the main physical ideas behind it and the main results are presented in a more accessible way for the reader not interested in technical details, has been presented in Charbonneau et al. 2014.

Part I General equations

In the following, we consider a system of $N$ identical spheres of diameter $D$ in a $d$-dimensional volume $V$, in the thermodynamic limit $N=\rho V$ and $V\rightarrow\infty$. At sufficiently high density, this system exhibits a glass phase, i.e., its dynamics arrests and an amorphous solid phase forms due to self-induced frustration. A strategy to obtain a thermodynamic description of the glass states has been derived in Kirkpatrick and Thirumalai 1988; Kirkpatrick and Wolynes 1987b; Monasson 1995. The idea behind these works is the following. In the dense supercooled regime the liquid phase splits in a collection of glass basins that can be defined as amorphous minima of an appropriate density functional Kirkpatrick and Thirumalai 1988. Hence if one extracts one liquid configuration in equilibrium, this configuration falls into one of the multiple glass basins. To describe this glass basin thermodynamically, one has to consider a second configuration that is coupled to the first one. Its partition function gives the properties of the glass. The coupling to the first configuration acts as an external quenched disorder, and that disorder can be treated using replicas, in analogy with spin glasses Kirkpatrick and Thirumalai 1987. This construction was made precise through the introduction of the so-called Franz-Parisi potential Franz and Parisi 1995; Franz and Parisi 1997; Cardenas et al. 1999. In this way one can easily follow the adiabatic evolution with density and temperature of glassy states that originate from an equilibrium liquid configuration (also known as the “state following” procedure Krzakala and Zdeborová 2010). Dynamically, this corresponds to preparing glassy states through a slow annealing.

A slightly different (and computationally simpler) strategy was proposed by Monasson Monasson 1995. Here $m$ replicas are coupled to an external random field that selects glassy states. After averaging out the external field one is left with $m$ coupled replicas. Because one uses here a completely random field, this procedure gives the properties of the “typical” states that exist at a given temperature or density Monasson 1995; Parisi and Zamponi 2010. In this way one can compute the dynamical glass transition density, the Kauzmann density (when it exists), and the so-called dynamical line (or threshold line) that delimits the region of existence of glassy states. Dynamically, it is expected that after a fast quench the system becomes arrested in the glassy states that lie on this dynamical/threshold line Cugliandolo and Kurchan 1993; Castellani and Cavagna 2005. This strategy is very efficient and was used in many previous studies, see e.g. Mézard 1999; Mézard and Parisi 2000 for pedagogical introductions and Mezard and Parisi 2012; Parisi and Zamponi 2010 for reviews of previous results. In this paper, as in the previous papers of this series Kurchan et al. 2012; Kurchan et al. 2013, we follow this simpler approach and consider an $m$-times replicated system in order to describe glassy states. We will see that this treatment is sufficient to describe the marginal stability of jammed packings and to obtain the critical properties around the jamming transition. The state-following (or Franz-Parisi) computation is also possible, but is left for a future publication. In the second part of this paper, we will explain more precisely what kind of information on the phase diagram can be obtained from this procedure.

We consider $m$ identical copies of the original hard sphere system, assuming that the spheres are arranged in molecules, each molecule containing an atom of each replica. A molecule is thus described by $m$ vectors $\bar{x}=\{x_{1}\cdots x_{m}\}$, each $x_{a}$ being $d$-dimensional. The molecular liquid is translationally invariant, so each atom in a molecule fluctuates around the center of mass of the molecule $X=m^{-1}\sum_{a}x_{a}$. We call the displacement of one atom $u_{a}=x_{a}-X$. In the following, $m\times m$ matrices in replica spaces are indicated by a hat, e.g. the matrix $\hat{\alpha}$ has elements $\alpha_{ab}$, and we denote by $\hat{\alpha}^{a,b}$ the $m-1\times m-1$ matrix obtained from $\hat{\alpha}$ by removing line $a$ and column $b$. Wide hats are used to denote quantities that have been properly scaled to be finite in the limit $d\rightarrow\infty$; for example, $\widehat{\varphi}=2^{d}\varphi/d$ is the scaled packing fraction, with the usual packing fraction $\varphi=\rho V_{d}$, where $V_{d}$ is the volume of a sphere of unit radius Kurchan et al. 2012.

The aim of this section is to write the entropy of the replicated system in a convenient way. With the notations introduced above, we consider the general form of the replicated entropy that has been derived in paper II of this series Kurchan et al. 2013:

where the matrix $\hat{\alpha}$ with elements $\alpha_{ab}=d\left\langle u_{a}\cdot u_{b}\right\rangle/D^{2}$ encodes the fluctuations of the replica displacement vectors $u_{a}$. The term proportional to $\widehat{\varphi}$ comes from the density-density interaction and we refer to it below as the “interaction term”. The other terms encode the entropic contributions and we refer to them all as the “entropic term”. It is useful to define the scaled mean square displacement between two replicas

It has been shown in Kurchan et al. 2012; Kurchan et al. 2013 that the elements of the matrices $\hat{\Delta}$ and $\hat{\alpha}$ are finite in the limit $d\rightarrow\infty$. Because the vector $u_{a}$ has $d$ components, Eq. (2) implies that the variance of the norm of the vector $u_{a}-u_{b}$ is of order $1/d$, while the variance of each one of its components is of order $1/d^{2}$. In order to find the thermodynamic entropy of the replicated liquid, the replicated entropy (1) must be optimized For integer $m>1$, the optimum is a maximum, but for real $m<1$ analytic continuation changes the sign of some eigenvalues and the optimum becomes a saddle point Mézard et al. 1987. with respect to the matrices $\hat{\alpha}$ or $\hat{\Delta}$.

In Kurchan et al. 2013 a generic expression for the interaction term $\mathcal{F}(2\hat{\alpha})$ has been derived, but we now want to write it in a more convenient form. It reads:

We will now suppose that the diagonal part of the matrix $\alpha$ is constant and $\alpha_{aa}=\alpha_{d}$ for all $a$, because this is true for matrices with a general $k$RSB structure. This implies that

and therefore the function ${\cal F}$ can be rewritten in terms of the matrix $\hat{\Delta}$ as

Note that the diagonal elements of the displacement matrix $\hat{\Delta}$ are all zero. Moreover, we temporarily define a matrix

We assume that the parameter $\Lambda$ is positive. Using the identity

In the above expression, the integration measure of the $\lambda_{a}$ corresponds to independent random variables with a Gaussian probability distribution. Then, the probability distribution of the random variable $h=-\min_{a}\left(\sqrt{\Lambda}\lambda_{a}+h_{a}\right)=\max_{a}\left(-\sqrt{\Lambda}\lambda_{a}-h_{a}\right)$ is given by

because the product of integrals is the probability that $h\geq-\sqrt{\Lambda}\lambda_{a}-h_{a},\ \forall a$, hence it is the probability that $h\geq\max_{a}\left(-\sqrt{\Lambda}\lambda_{a}-h_{a}\right)$, which is nothing but the cumulative distribution of $h$. Here we introduced the function $\Theta(z)=(1+\text{erf}(z))/2$ and $\int_{-h}^{\infty}{\cal D}\lambda=\Theta(h/\sqrt{2})$. Therefore

By defining the Gaussian kernel $\gamma_{a}(z)=e^{-z^{2}/(2a)}/\sqrt{2\pi a}$, one obtains the relation

where we defined the $\gamma_{a}\star f$ operation as the convolution of the Gaussian kernel with the function $f(h)$. Then

Using this relation we obtain the compact result

The replicated entropy is therefore given by Eq. (1), with the function ${\cal F}$ given in Eq. (15). This expression is quite general and will be used to study $k$-step replica symmetry breaking schemes in the following. Note that it has been derived under the only assumption that the diagonal elements of $\alpha_{ab}$ are all equal.

III The replicated entropy for hierarchical RSB matrices

In spin glasses, it has been shown that a correct description of the system can be achieved by considering a special class of matrices, known as hierarchical $k$RSB matrices Mézard et al. 1987. Therefore, we want to specialize the general expression of the replicated entropy given by Eqs. (1) and (15) to these matrices. In Sec. III.1, we introduce hierarchical $k$RSB matrices and discuss some of their properties. In Sec. III.2, we compute the entropic term of the replicated entropy for a $k$RSB matrix, and in Sec. III.3 we compute the interaction term. In Sec. III.4, we present the final explicit expressions for the 1RSB, 2RSB, $k$RSB and fullRSB case.

The structure of hierarchical $k$RSB matrices is well-known and here we just summarize it briefly. Remember that with respect to the formalism of Mézard et al. 1987 we have an important difference, in that the diagonal elements of the matrices we consider are determined by the condition that $\sum_{b}\alpha_{ab}=0,\ \forall b$. The simplest class is that of 1RSB matrices In the standard notation of Mézard et al. 1987 this would be called a replica-symmetric matrix, but remember that here we are using the Monasson’s real replica scheme Monasson 1995 where we consider $m$ coupled replicas and treat $m$ as a parameter to select different metastable states. It is a standard convention to denote a RS matrix in the Monasson’s scheme as a “1RSB matrix”: the reason is that in models with quenched disorder the two schemes are indeed equivalent. , that has been studied in the first Kurchan et al. 2012 and second Kurchan et al. 2013 paper of this series. It corresponds to $\alpha_{ab}=-\widehat{\alpha}_{1},\ \forall a\neq b$:

At the 3RSB level, the blocks of $m_{1}$ replicas are each divided in $m_{1}/m_{2}$ sub-blocks of $m_{2}$ replicas, and the construction can be iterated to any desired level of RSB (which we denote $k$RSB). Note that diagonal elements of hierarchical matrices are all equal.

For each $k$, these “hierarchical” matrices form a closed algebra. Moreover, in the interesting case where one performs an analytic continuation to $m<1$, a generic hierarchical matrix $\hat{Q}$ can be parametrized by its diagonal element $q_{d}$ and a single function $q(x)$. This is done as follows. If we formally define $m_{0}=m$ and $m_{k}=1$, we can observe that for integer $m>1$, one has $m_{k}\equiv 1<m_{k-1}<m_{k-2}<\cdots<m_{1}<m_{0}\equiv m$. When $m<1$, these inequalities are reversed Mézard et al. 1987 and one has $1>m_{k-1}>m_{k-2}>\cdots>m_{1}>m>0$. Then we can define $q(x)$ to be a piecewise constant function for $x\in$, which in the interval $x\in[m_{i-1},m_{i}]$ takes the value of the elements $q_{i}$ in the corresponding sub-block (with $i=1\cdots k$), and $q(x)=0$ for $x\in[0,m_{0}]$. We write this parametrization as $\hat{Q}\leftrightarrow\{q_{d},q(x)\}$.

The matrices $\hat{\alpha}$ are therefore defined by $\hat{\alpha}\leftrightarrow\{\alpha_{d},-\alpha(x)\}$ where $\alpha(x)$ is a piecewise constant function given by the $\widehat{\alpha}_{i}$ in each block, while $\alpha_{d}$ is fixed by the condition $\sum_{b}\alpha_{ab}=0$ and is given by

For a given matrix $\alpha_{ab}$, one has a corresponding matrix $\Delta_{ab}=2\alpha_{aa}-2\alpha_{ab}$, which is therefore parametrized as $\hat{\Delta}\leftrightarrow\{0,\Delta(x)\}$ with

The seemingly strange conventions that we adopted for the signs of the above parametrization are justified by the fact that with this choice both $\alpha(x)$ and $\Delta(x)$ are positive functions. Clearly this is the case by definition for $\Delta(x)$, and from Eq. (20) we deduce that the same must be true for $\alpha(x)$. Moreover, $\Delta(x)$ must be a decreasing function of $x$ (contrary to the usual overlap function Mézard et al. 1987). This is because larger values of $x$ correspond to inner blocks of the matrix $\Delta_{ab}$, hence to replicas that are “closer” to each other in the usual interpretation, and therefore must have a smaller value of $\Delta$. An example is given in Fig. 1. Of course, the above construction can be generalized to the case where the function $\Delta(x)$ is allowed to have a continuous part: this can be thought as an appropriate limit of the $k$RSB construction when $k\rightarrow\infty$ and is therefore called “fullRSB” or “$\infty$RSB”.

III.2 The algebra of hierarchical matrices and the entropic term

In order to compute the entropic term, we need to compute $\log\det\hat{\alpha}^{m,m}$. For this we need to recall some standard results for the algebra of hierarchical matrices. Let us consider a generic hierarchical matrix $q_{ab}$, parametrized by the corresponding function $q(x)$ and with diagonal element $q_{d}$. Although our matrices $\hat{Q}_{m}$ are $m\times m$ matrices with fixed $m<1$, we can think to embed them in a $n\times n$ dimensional matrix $\hat{Q}_{n}$ with $n\rightarrow 0$, and just set to zero the element of the outermost block. This is consistent with the fact that $q(x)$ is defined in $x\in$, and it corresponds to the special choice that $q(x)$ vanishes for $0\leq x<m$. We can therefore use standard results for the algebra of $n\times n$ hierarchical matrices, in the limit $n\rightarrow 0$, see e.g. Mézard and Parisi 1991. The important result of Mézard and Parisi 1991 that we need is that, introducing the notations

we have for the determinant (Mézard and Parisi 1991, Eq.(AII.11))

and the inverse matrix is parametrized by (Mézard and Parisi 1991, Eq.(AII.7))

To adapt these results to our case, we note that if the outermost blocks vanish, we have $q(0)=0$. Furthermore, $q(x)=[q](x)=0$ for $x<m$. Finally, $\det\hat{Q}_{n}=(\det\hat{Q}_{m})^{n/m}$, hence $\log\det\hat{Q}_{m}=(m/n)\log\det\hat{Q}_{n}$, and the diagonal element of the inverse of $\hat{Q}_{n}$ and $\hat{Q}_{m}$ are identical. We conclude that in the $m\times m$ space

In order to compute $\log\det\hat{\alpha}^{m,m}$, we introduce a matrix $\hat{\beta}$ which is “regularized” in such a way that $\sum_{b}\beta_{ab}=\varepsilon$, and $\hat{\alpha}=\lim_{\varepsilon\rightarrow 0}\hat{\beta}$. For instance we can choose

for any $x_{0}\in[m,1]$. The matrix $\hat{\beta}$ is invertible, hence we can use the relation $\det\hat{\beta}^{m,m}=(\hat{\beta}^{-1})_{mm}\det\hat{\beta}$. Using Eq. (24) we get

Expressed in terms of the function $\Delta(x)$ by means of Eq. (20), we have that

When specialized to the 1RSB case Eq. (16), which corresponds to $\alpha(x)=\widehat{\alpha}_{1}$ and $\Delta(x)=\widehat{\Delta}_{1}=2m\widehat{\alpha}_{1}$, we get

which reproduces the results of Kurchan et al. 2012; Kurchan et al. 2013. The 2RSB solution is parametrized by a step function

Writing explicitly the $k$RSB expression requires the introduction of a function

In fact, if $\Delta(x)$ is a $k$RSB piecewise constant function and using the notations of Fig. 1, it is easy to check that $G(x)$ is also a piecewise constant function parametrized by $\widehat{G}_{i}$, with

Note that from this result we can recover the 1RSB and 2RSB results obtained above.

III.3 The interaction term

We now compute the interaction term of the replicated entropy for a generic hierarchical matrix $\Delta_{ab}$ parametrized by $\Delta(x)$. We start from Eq. (15), and we need to compute the function

Because $\Delta_{ab}$ is a hierarchical matrix, this computation can be done by taking the derivative with respect to the external fields in a hierarchical way Duplantier 1981. Let us define the matrix $I^{m_{i}}_{ab}$, which has elements equal to 1 in blocks of size $m_{i}$ around the diagonal, and zero otherwise. In other words, the matrix $I^{m_{i}}_{ab}$ is parametrized by a function $I^{m_{i}}(x)=1$ for $m_{i}\leq x\leq 1$ and zero otherwise. Then, recalling the notations of Fig. 1, and noting that $I^{m_{k}=1}_{ab}=\delta_{ab}$, one can easily check that

Inserting this form in Eq. (40), one obtains a sequence of differential operators acting on the product of theta functions, each of them being the sum of partial derivatives inside a block. This sequence of operations can be written as a recursion; since the procedure is very well explained in Duplantier 1981, we only report the main results here. When acting with the term containing $\widehat{\Delta}_{k}$ we obtain, recalling Eq. (13):

Then, the action of each of the terms $i=k-1,\cdots,0$ induces a recursion of the form

Note that the last of these iterations can be written explicitly as

This recursive procedure allows us to compute easily $g(m,h)$ for any $k$, and, according to Eq. (15), the function ${\cal F}(\Delta)$.

The last iteration might seem problematic, because the kernel $\gamma_{a}$ has a negative parameter $a=-\widehat{\Delta}_{1}$ and therefore cannot be represented as in Eq. (13). Luckily enough, this last iteration can be eliminated. In fact, we have

Using Eq. (42) we obtain, at the 1RSB level,

and so on. At the 1RSB level, we reproduce the result of Kurchan et al. 2012; Kurchan et al. 2013.

that has to be solved with the initial condition (42). Note that the partial differential equation is well defined because $\dot{\Delta}(x)\leq 0$. This equation has to be integrated from $x=1$ down to $x=m_{1}$. The resulting $g(m_{1},h)$ can be inserted in Eq. (46) to obtain ${\cal F}(\Delta)$.

The partial differential equation written above can be also put in a more convenient form by introducing the function

Remarkably enough, the equations written above are identical to those of the SK model Mézard et al. 1987; Duplantier 1981. The only difference is in the initial condition for the Parisi equation.

III.4 1RSB, 2RSB, kkRSB, fullRSB expressions of the replicated entropy

We therefore obtained the expression of the replicated entropy, under the assumption that the matrices $\hat{\alpha}$ and $\hat{\Delta}$ are hierarchical $k$RSB matrices, and we now summarize the results that we obtained at different levels of RSB for the replicated entropy (1). At any level of $k$RSB, the replicated entropy has the form

where ${\cal S}_{k{\rm RSB}}$ has been defined in such a way that it contains the non-trivial dependence on $\hat{\alpha}$ and it has a good limit for $m\rightarrow 0$ (as we will discuss later).

At the 1RSB level, using Eqs. (30), (47) and (46), we obtain

which coincides with the previously derived results Parisi and Zamponi 2010; Kurchan et al. 2012; Kurchan et al. 2013. At the 2RSB level, we have, using Eqs. (34), (48) and (46),

At the generic $k$RSB level we obtain from Eq. (39), (46), (42), (43), (44):

Finally, the fullRSB expression is, from Eqs. (29), (46), (51), (52), (53):

IV Variational equations

The replicated entropy depends on the function $\Delta(x)$ that parametrizes the hierarchical matrix $\hat{\Delta}$. This function is determined by optimization of the replicated entropy through a variational principle Mézard et al. 1987. In this section we derive the variational equations for the function $\Delta(x)$ that are obtained by optimization of the free energy, i.e. by imposing the equation $\frac{\partial s[\hat{\alpha}]}{\partial\Delta_{ab}}=0$. It is well known in the context of spin glasses Mézard et al. 1987 and structural glasses Mezard and Parisi 2012; Parisi and Zamponi 2010 that, for integer $m>1$, this corresponds to the usual maximization of the entropy or minimization of the free energy with respect to the variational parameters $\alpha_{ab}$, but when the problem is analytically continued to real $m<1$, the solution of these equations does not correspond to a maximum of the entropy, but to a saddle point. This not very important because to really characterize the stability of the saddle point one has first to compute the matrix of second derivatives $\frac{\partial^{2}s[\hat{\alpha}]}{\partial\Delta_{ab}\partial\Delta_{cd}}$, and then perform the analytic continuation of its eigenvalues to $m<1$. The continued eigenvalues are required to be positive.

In this section we derive the variational equations for $\Delta(x)$, and we postpone a partial analysis of the stability matrix to the following sections. In Sec. IV.1 we derive the equation for $\widehat{\Delta}_{i}$ in the case of a $k$RSB structure. In Sec. IV.2 we derive the fullRSB equations in two equivalent ways: first by using Lagrange multipliers, and then by taking the $k\rightarrow\infty$ limit of the $k$RSB equations.

We consider first the $k$RSB solution for fixed $k$. We start from Eq. (57) and we want to impose the condition $\frac{\partial{\cal S}_{k{\rm RSB}}}{\partial\widehat{\Delta}_{i}}=0$. To do this, we consider Eq. (40) and (46). We have, without taking into account that the matrix $\Delta_{ab}$ is symmetric and for $a\neq b$:

The next step is to take the derivative with respect to $\widehat{G}_{i}$ of the $k$RSB free energy, which we write in the form

and make use of the relation (38). We get

We observe from Eq. (61) that $N_{i}(m,h)$ behaves like a Gaussian for $h\rightarrow\pm\infty$, therefore we can safely integrate by parts, and the last expression can be written as

using the same trick as in Eq. (45). Finally, we can show that

provided $P(m_{i},h)$ satisfies the following recurrence equations:

Eqs. (73), (71)-(72), (42)-(43), (38) constitute a set of closed equations for the $\widehat{G}_{i}$, or equivalently the $\widehat{\Delta}_{i}$. They can be solved by the following iteration: starting from a guess for $\widehat{\Delta}_{i}$, one can solve first the recurrence (42)-(43) and then the recurrence (71)-(72). From the solutions one can compute the new $\widehat{G}_{i}$ using Eq. (73) and from these the new $\widehat{\Delta}_{i}$ using Eq. (38).

IV.2 Variational equations for the fullRSB solution

over $\Delta(x)$, $f(x,h)$, $P(x,h)$, $f(m,h)$ and $P(1,h)$. The first two equations can be obtained by taking the variation with respect to $P(x,h)$ and $f(x,h)$

Taking the variation over $P(1,h)$ and $f(m,h)$ we obtain

Finally, taking the variation of $G(x)$ (for $x\neq 1$ and $x\neq m$) we obtain the following equation

The system of Eqs. (75)-(79) can be in principle solved numerically, with the following procedure:

from this one solves Eq. (75) with boundary condition (77) to get $f(x,h)$;

then one can solve Eq. (76) with boundary condition (78) to obtain $P(x,h)$;

We now take the continuum limit of the discrete $k$RSB equations following the strategy of section III.3.3. It was already shown in that section that in this limit the recurrence equations for $g(x,h)$, Eqs. (42) and (43), become the Parisi equation (75) with boundary condition (77). The boundary condition for $P(x,h)$ in the discrete, Eq. (71), is clearly equivalent to the one in the continuum, Eq. (78). It is quite simple, following the lines as in section III.3.3, to show that Eq. (72) becomes, in the continuum limit, Eq. (76), so we do not report the derivation.

It remains to derive Eq. (79). We start from Eq. (73) which becomes in the continuum limit

We now show that Eq. (79) and Eq. (80) are equivalent. This amounts to showing that

which shows that Eq. (81) holds at $x=m$. Next, we compute the derivative with respect to $x$ of the arguments of the integrals that appear in Eq. (81). We have, using Eqs. (75) and (76), that

This proves that the derivatives of the two sides of Eq. (81) with respect to $x$ coincide and therefore completes the proof of Eq. (81), and of the equivalence of Eq. (79) and Eq. (80). We have therefore derived the set of fullRSB equations Eqs. (75)-(79) in two independent ways.

V Derivation within the Gaussian ansatz

Although the above results have been derived from an exact evaluation of the replicated entropy following Kurchan et al. 2012; Kurchan et al. 2013, they could be equivalently obtained from a suitable Gaussian ansatz in replica space Kurchan et al. 2012. Here we discuss the appropriate form of this ansatz. This approach is interesting for two reasons: it sheds some light on the physical interpretation of both the $k$RSB ansatz and the function $P(m_{i},h)$ introduced in Sec. IV, and it opens the way to extend the result (in an approximate way) to finite dimensions, following the approach of Parisi and Zamponi 2010; Berthier et al. 2011.

In general, the replicated entropy is a functional of the single molecule density $\rho(\bar{x})$, where $\bar{x}=\{x_{1},\cdots,x_{m}\}$ and $x_{a}$ are the $d$-dimensional vectors corresponding to the positions of particles in the molecule. In terms of this object, the replicated entropy for $d\rightarrow\infty$ is given in (Kurchan et al. 2012, Eqs.(2)):

In Eq. (85), we introduced a generic interparticle potential $v(r)$ at inverse temperature $\beta$. In this paper we restrict ourselves to the hard sphere potential, where

but since the Gaussian derivation allows one to consider a generic potential, it will be useful to write the expressions for a generic $v(r)$ because this will be surely useful for future applications, e.g. to the soft sphere case following Berthier et al. 2011. Note that in finite dimensions, the replicated entropy can be expressed as an infinite sum of diagrams Parisi and Zamponi 2010, but for $d\rightarrow\infty$, and for potentials that have a hard core or a properly scaled soft core, one can truncate the series at the lowest order Frisch and Percus 1999; Parisi and Zamponi 2010, hence obtaining Eq. (85).

The Gaussian ansatz consists in making an appropriate Gaussian assumption on the function $\rho(\bar{x})$, thus introducing a set of variational parameters that are related to the matrix $\hat{\alpha}$ considered above. In Sec. V.1 we discuss the proper Gaussian parametrization of the density function, we compute the entropic term and we show that it has the same form as the one we found before. Finally, in Sec. V.2 we compute the interaction term, we discuss the connection with the correlation function, and we show that in the limit $d\rightarrow\infty$ we recover the results obtained above.

For the entropic term, we just need to recall some results already discussed in Kurchan et al. 2012; Kurchan et al. 2013, to which we refer for details. Thanks to translational invariance, we can choose a parametrization of $\rho(\bar{x})$ in terms of vectors $\bar{u}=\{u_{1},\cdots,u_{m}\}$ such that $\sum_{a=1}^{m}u_{a}=0$, and 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}$ give the average of the squared replica displacements

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}$. Hence, the relative replica displacements $\operatorname{\mathsf{D}}_{ab}$ are given by

These corresponds to the physical “overlaps”, i.e. the mean square displacements between different states. By comparison with Eq. (2), we see that $\Delta_{ab}=d^{2}\operatorname{\mathsf{D}}_{ab}/D^{2}$ and $\alpha_{ab}=d^{2}A_{ab}/D^{2}$.

which corresponds indeed to the entropic term in Eq. (1), with the rescaling $\alpha_{ab}=d^{2}A_{ab}/D^{2}$. Apart from this rescaling, all the derivation of section III.2 can therefore be repeated within the Gaussian ansatz and we arrive to exactly the same results for the entropic term, which therefore has the same form both in infinite dimensions and in finite dimensions within the Gaussian ansatz.

V.2 kkRSB Gaussian parametrization and the interaction term

To compute the interaction term, we need to find a simpler parametrization of the $k$RSB form of the Gaussian single molecule density. Recall first that in the 1RSB case one has $A_{ab}=A_{1}(\delta_{ab}-1/m)$, hence $\operatorname{\mathsf{D}}_{aa}=0$ and for $a\neq b$, $\operatorname{\mathsf{D}}_{ab}=\operatorname{\mathsf{D}}_{1}=2A_{1}$. In that case Eq. (87) can be written as Mézard and Parisi 1999; Mezard and Parisi 2012; Parisi and Zamponi 2010:

is a $d$-dimensional centered Gaussian. The simplest proof of the equivalence of Eq. (87) and Eq. (91) for 1RSB matrices is obtained by computing for $a\neq b$

This shows that Eqs. (87) and (91) have the same (vanishing) first and (finite) second moments, and therefore they coincide because a Gaussian is only specified by its first two moments.

Let us consider next a matrix $A_{ab}$ with a 2RSB structure. In this case we have two physical overlaps, the Edwards-Anderson $\operatorname{\mathsf{D}}_{2}$ and the inter-state $\operatorname{\mathsf{D}}_{1}>\operatorname{\mathsf{D}}_{2}$. Let us call $B_{i}=\{1+(i-1)m_{1},\cdots,im_{1}\}$, with $i=1\cdots m/m_{1}$ the $i$-th block of the 2RSB matrix. Then we can write

It is easy to check, through a computation very similar to Eq. (93), that $\left\langle(u_{a}-u_{b})^{2}\right\rangle=\left\langle(x_{a}-x_{b})^{2}\right\rangle$ is equal to $\operatorname{\mathsf{D}}_{2}$ if $a,b$ belong to the same block $B_{i}$, and $\operatorname{\mathsf{D}}_{1}$ otherwise (and of course zero if $a=b$), hence we conclude that Eq. (94) is identical to Eq. (87) for a 2RSB matrix.

Generalizing this construction we easily obtain the structure of the $k$RSB Gaussian parametrization. We do not write it explicitly because this would require introducing a heavy notation for block indices, however it is clear that one should couple each group of replicas in the innermost blocks to reference points $X^{k}_{i}$, then group these reference points in blocks, each block being coupled to a reference $X^{k-1}_{i}$, and so on until the most external reference points are coupled to a single $X^{1}$. This tree structure, and how it enters in the computation of the interaction term, is illustrated in Fig. 2.

We now show how the computation of the interaction term is performed recursively. We note that the term $-1$ in the Mayer function is obviously independent of the integration variables and we neglect it for the moment. The procedure is illustrated in the left panel of Fig. 2. Consider the innermost blue dots, which represent the coordinates of two atoms. They interact through the “bare” interaction $e^{-\beta v(x-y)}$. When we integrate over $x,y$, we generate a “bare” interaction between the reference points $X^{k},Y^{k}$, also marked by blue dots. This interaction is

where we denoted by $\otimes$ the $d$-dimensional convolution. Now, each of the $k$-level reference positions are coupled to $m_{k-1}/m_{k}$ atoms, therefore the total bare interaction between each pair $X^{k},Y^{k}$ is ${\cal G}(1,X^{k}-Y^{k})^{m_{k-1}/m_{k}}$. Now we can integrate over these variables to obtain the bare interaction between $(k-1)$-level reference positions, and so on. It is easy to see that at each step of the iteration we have (we now use a similar notations as that of the previous sections):

The bare interaction at level $i$ is ${\cal G}(m_{i},r)^{m_{i-1}/m_{i}}$. Iterating this procedure, we finally obtain the bare interaction of the most external reference points which is ${\cal G}(m_{1},r)^{m/m_{1}}$. Therefore the interaction term per particle is

Before we show that these expressions give exactly the same result obtained before in the limit $d\rightarrow\infty$, let us consider the total “dressed” interaction between two $i$-level reference points. By dressed interaction we mean that we do not only integrate over the points “down the tree”, that are at level $j>i$, as we did before; we also integrate over all the other points. For the outermost reference points $X^{1},Y^{1}$, the dressed interaction is clearly given by

as we discussed before. Now, if we consider the points at level $2$, we have:

The interpretation of this relation should be straightforward by looking at right panel of Fig. 2. In fact, the first term is the bare interaction that comes from points down the tree, while the second term is the contribution that comes from points up the tree, in which we divide ${\cal P}$ by ${\cal G}$ because we need to remove the contribution of the branch of the tree under consideration. Iterating, at level $i$ we have

At level $k$ this procedure gives the dressed interaction between the innermost reference points $X^{k},Y^{k}$. The last step allows one to obtain the dressed correlations of two atoms $x,y$ with $r=x-y$. This is the two-body effective potential $\phi_{\rm eff}(r)$ of Parisi and Zamponi 2010; Berthier et al. 2011 and we obtain

One can easily check that in the 1RSB case this result gives back the effective potential used in Parisi and Zamponi 2010; Berthier et al. 2011. This result is particularly interesting for two reasons: on the one hand, in the limit $d\rightarrow\infty$ (or, in finite $d$, under the low-temperature approximation of Ref. Berthier et al. 2011), the effective potential coincides with the pair correlation function of the glass, $g_{\rm g}(r)$. On the other hand, in finite dimension one could plug this effective potential into some liquid theory integral equations to compute an approximation of the replicated entropy, following the strategy of Parisi and Zamponi 2010. We do not discuss this second issue here, leaving it for future work, and we keep our focus instead on the limit $d\rightarrow\infty$.

To take the $d\rightarrow\infty$ limit, we assume that the potential has the form $e^{-\beta v(r)}=e^{-\widehat{v}[d(|r|-D)/D]}$, and $e^{-\widehat{v}(h)}$ is a finite function Remember that in general we use a wide hat to denote quantities that are properly rescaled to be finite when $d\rightarrow\infty$. when $d\rightarrow\infty$. The hard sphere potential, the only one that we consider in this paper, has this form with $e^{-\widehat{v}(h)}=\theta(h)$. Note, however, that the interaction potential enters only in the initial condition for the function ${\cal G}(1,r)$.

We now show briefly that in the limit $d\rightarrow\infty$ the equations we just derived give back Eq. (57). First of all, we note that all the interaction functions ${\cal G}(m_{i},r)$ and ${\cal P}(m_{i},r)$ are rotationally invariant and therefore depend only on $|r|$. Moreover, all these functions tend to 1 for $|r|\rightarrow\infty$ and they decrease fast to zero when $|r|\ll D$. Actually, as we show below, for $d\rightarrow\infty$ the growth of these functions from 0 to 1 happens on a scale $\sim 1/d$ around $|r|=D$. We therefore define $|r|=D(1+h/d)$ and we consider ${\cal G}(m_{i},h)$ and ${\cal P}(m_{i},h)$ to be functions of $h$. We will see that with this convention, $h$ is the same variable that enters in the equations of the previous sections.

Next, to obtain a non-trivial limit we scale the parameters introducing $\widehat{\Delta}_{i}=d^{2}\operatorname{\mathsf{D}}_{i}/D^{2}$. Then, following exactly the strategy of (Parisi and Zamponi 2010, Appendix C.2.d) through the use of bipolar coordinates, one can see that in the limit $d\rightarrow\infty$, with the scaling $r=D(1+h/d)$ and $u=D(1+z/d)$, we have

With a similar reasoning, Eq. (96) becomes

This equation coincides with Eq. (43) if one makes for all $i$ the identification

and recalling that for hard spheres $e^{-\widehat{v}(h)}=\theta(h)$. Hence the interaction term is (using that the solid angle is $d\,V_{d}$ with $V_{d}$ the volume of a unit sphere)

which provides the exact result of Eq. (57).

Next, we analyze the behavior of ${\cal P}$ for $d\rightarrow\infty$. The initial condition from Eqs. (98) and (71) is

With a similar reasoning as before, one can show that identifying

for all $i$, Eq. (100) becomes in fact identical to Eq. (72) in the limit $d\rightarrow\infty$. We obtain from this identification a deep physical interpretation of the function $P(m_{i},h)$, which turns out to be related to the dressed interaction of the reference positions at level $i$, ${\cal P}(m_{i},h)$, by Eq. (107).

Finally, we can write Eq. (101) for $d\rightarrow\infty$. We get

This is a very important result because it allows one to obtain structural information about the pair correlation from the knowledge of the functions $P(m_{i},h)$ and $g(m_{i},h)$.

Part II Extraction of the results from the equations

In the following sections, we will investigate the phase diagram that one obtains from the study of the $k$RSB equations, and we derive the scaling properties at large pressure. The way in which one has to extract physical information from the replicated entropy has been explained in many reviews Mézard et al. 1987; Monasson 1995; Mezard and Parisi 2012; Parisi and Zamponi 2010. Although we will give additional details along the way, we assume that the reader is familiar with this kind of computations.

Before extracting the physics, let us summarize here the $k$RSB equations (57) together with the variational equations (71)-(72) and (73). These are

Because (as we will see below) we are mostly going to work at small $m$, and the jamming limit corresponds to $m\rightarrow 0$, it is convenient to write the equations in scaled variables that remain finite when $m\rightarrow 0$. These are $y_{i}=m_{i}/m$ (keeping in mind that $y_{0}=1$ and that $y_{k}=1/m$ diverges with $m$ and will play the role usually played by temperature), $\widehat{f}(y_{i},h)=mf(m_{i},h)$, $\widehat{\gamma}_{i}=\widehat{G}_{i}/m$, from which it follows that $\widehat{\Delta}_{k}=m\widehat{\gamma}_{k}$ while all the other $\widehat{\Delta}_{i}$ remain finite for $m\rightarrow 0$. It will also be convenient for numerical reasons to introduce $\widehat{P}(y_{i},h)=e^{-\widehat{\Delta}_{1}/2}e^{-h}P(m_{i},h)$. Note also that $\widehat{\Delta}_{i}-\widehat{\Delta}_{i+1}=(\widehat{\gamma}_{i}-\widehat{\gamma}_{i+1})/y_{i}$. In terms of these variables, and introducing auxiliary variables $\widehat{\kappa}_{i}$ (not to be confused with the exponent $\kappa$ discussed above) we have:

To solve numerically these equations it is convenient to have some control on the asymptotic behavior of the functions when $h\rightarrow\pm\infty$. We start by the function $\widehat{f}(y_{i},h)$. From the initial condition we see that $\widehat{f}(1/m,h\rightarrow\infty)=0$ and $\widehat{f}(1/m,h\rightarrow-\infty)\sim-h^{2}/(2\widehat{\gamma}_{k})$. Inserting these asymptotes in the evolution equation for $\widehat{f}(y_{i},h)$, one can show that

From this, we obtain that $\widehat{P}(y_{1},h\rightarrow-\infty)=0$ while $\widehat{P}(y_{1},h\rightarrow\infty)=e^{-\widehat{\Delta}_{1}/2}$. As a consequence, from the recurrence equation for $\widehat{P}(y_{i},h)$ one can show that

To obtain a simpler asymptotic behavior it is convenient to make a change of variable:

in such a way that the leading asymptotic term of $\widehat{f}(y_{i},h)$ in Eq. (110) is subtracted from $\widehat{j}(y_{i},h)$. Then we have

is not a symmetric function of $h$ and $z$, nor a function of $h-z$. However, the advantage of this formulation is that the kernel $K$ is an almost Gaussian function which is well behaved, and all the other functions that appear in the integrals are smooth. This allows for a stable numerical evaluation of the integrals.

Note that Eqs. (113) admit a perfectly smooth $m\rightarrow 0$ limit. First of all one has to set $1/y_{k}=m$, $\widehat{\Delta}_{k}=0$. Then, using the large $\lambda$ development of $\Theta(-\lambda/\sqrt{2})$, one can easily show that

and therefore $\widehat{j}(y_{k},h)=0$. All the other equations remain identical to the case $m>0$.

VI.2 The continuum limit

It will be convenient for later purposes to write explicitly the continuum limit of the equations in terms of scaled variables. These are

Having formulated the $k$RSB equations in a convenient way, we now proceed to extract the physical results from them. However, because the numerical solution of these equations is not trivial, we first investigate a certain number of asymptotic limits in which analytical results can be obtained.

VII Perturbative 2RSB solution around the Gardner line

The Gardner transition Gardner 1985 separates the region where the 1RSB solution is stable from the one where it is unstable. In our problem, the 1RSB solution is stable in a certain region of the phase diagram, and it becomes unstable on a line that has been computed and characterized in the previous paper of this series Kurchan et al. 2013. The aim of this section, following the analysis of Rizzo 2013, is to perform a perturbative computation around the Gardner line, in the region where the 1RSB solution is unstable, and discuss the existence of a 2RSB (or fullRSB) solution. In Sec. VII.1 we perform the perturbative computation and show that a 2RSB solution only exists in a certain region of the Gardner line; in Sec. VII.2 we discuss the behavior of the perturbative solution at large densities and pressures.

We start our analysis of the $k$RSB equations by examining what happens around the Gardner transition line that has been computed in Kurchan et al. 2013. We will follow closely the analysis of Refs. Montanari and Ricci-Tersenghi 2003; Rizzo 2013. Fig. 3 reports a schematic phase diagram in the $\widehat{\varphi}$, $m$ plane. The reader should keep in mind that $m\propto 1/p$ where $p$ is the reduced pressure Kurchan et al. 2013. Let us summarize briefly the phase diagram. A 1RSB solution exists above the 1RSB dynamical line, i.e. for $\widehat{\varphi}>\widehat{\varphi}_{\rm d}^{\rm 1RSB}(m)$ or equivalently $m>m_{\rm d}^{\rm 1RSB}(\widehat{\varphi})$. Above the Gardner line, i.e. for $m\geq m_{\rm G}(\widehat{\varphi})$ or $\widehat{\varphi}>\widehat{\varphi}_{\rm G}(m)$, this 1RSB solution is stable Kurchan et al. 2013 and one can extract the physical results from it Parisi and Zamponi 2010. Below this line, the 1RSB solution is unstable and we look for a $k$RSB solution with $k>1$. Since the instability is due to a vanishing mode Kurchan et al. 2013, we expect that the 1RSB solution will transform continuously in a $k$RSB solution with $k>1$, and we therefore start by looking at the new solution by doing perturbation theory around the 1RSB solution in the vicinity of the Gardner line. Based on the analogy with spin glass models Rizzo 2013, we expect the new solution to be a fullRSB one, but in perturbation theory it is enough to consider a 2RSB solution because the breaking is small and the two perturbative computations give identical results Caltagirone et al. 2011.

We have found in Kurchan et al. 2013 that the instability of the 1RSB solution is due to the vanishing of the “replicon” eigenvalue $\lambda_{R}(m)$. It is therefore natural to consider a perturbation of the 1RSB matrix that is proportional to the subspace associated with the vanishing mode. We need to consider the cubic expansion in this direction and check if the cubic terms can stabilize the negative quadratic part. Hence we consider a perturbative 2RSB matrix of the form:

(remember that the matrices $I^{m_{i}}_{ab}$ are defined in Sec. III.3). The matrix $\delta\hat{\alpha}$ belongs to the replicon subspace Temesvári et al. 2002. It has been shown in Temesvári et al. 2002 that the cubic terms of the expansion of the 2RSB entropy around the 1RSB solution are eight, but if the perturbation matrix $\delta\hat{\alpha}$ is in the replicon subspace, only the two terms proportional to the coefficients $w_{1}$ and $w_{2}$ that have been analyzed in Kurchan et al. 2013 survive. In practice we have

where $\lambda_{R}(m)$ is the replicon eigenvalue that has been computed in Kurchan et al. 2013. The relation between $\lambda_{R}(m)$ and $\widehat{\lambda}_{R}(m,m_{1})$ can be exploited to obtain the replicon eigenvalue from the 2RSB entropy and it gives

It expresses the fact that the replicon eigenvalue can be obtained from the second derivative with respect to $\delta\widehat{\alpha}_{1}$ of the replicated entropy computed on a matrix of the form (117). We want now to express $W$ explicitly in terms of $w_{1}$ and $w_{2}$ that we have already computed in Kurchan et al. 2013. Note that $w_{1}(m)$ and $w_{2}(m)$ depend only on $m$ and not on $m_{1}$ as a consequence of the fact that they are computed on a 1RSB matrix that does not depend on $m_{1}$. To do this one can use the following relations between matrices $\hat{I}^{m_{i}}$:

to obtain that the matrix $\hat{r}$ defined in (117) satisfies

To obtain the perturbative 2RSB solution in the $(m,\widehat{\varphi})$ plane we search for a non trivial stationary point solution for the expression (118). The trivial 1RSB solution $\delta\widehat{\alpha}_{1}=0$ can always be found but it is unstable below the Gardner line where $\lambda_{R}(m)<0$. To find the non-trivial solution we first optimize over $\delta\widehat{\alpha}_{1}$ and then we optimize over the breaking point $m_{1}$. The saddle point equation for $\delta\widehat{\alpha}_{1}$ gives

The entropy as a function of $m_{1}$ is obtained by plugging the above expression in (118) and it gives

Now we should search for the extremum in $m_{1}$. We obtain that the breaking point is

where $\lambda(m)$ is the Mode-Coupling theory exponent parameter that has been discussed in Kurchan et al. 2013. However, as it is usual in replica computations, the saddle point solution for the breaking point should satisfy $m<m_{1}<1$. This implies that a 2RSB solution exists only if

Note that, because we are perturbing around the 1RSB solution close to the Gardner line where the replicon mode vanishes, all the quantities $\widehat{\lambda}_{R}(m)$, $w_{1}(m)$ and $w_{2}(m)$ are computed on the Gardner line.

We now follow the discussion of Rizzo 2013. We know that at $m=1$ the Gardner line corresponds to the dynamical transition point and at that point $\lambda(m=1)=\lambda_{\rm MCT}=0.70698<1=m$ Kurchan et al. 2013. Hence, by continuity, close to $m=1$ we have $\lambda(m)<m$ and there is no non-trivial 2RSB solution. We conclude, following Rizzo 2013, that in this region only the 1RSB solution exists down to the Gardner line, and there is no other solution below the Gardner line. However, the condition (126) tells us that there might exist a point in the $(m,\widehat{\varphi})$ plane where

so that below this point a perturbative 2RSB solution can be found. The point $m^{*}$ can be computed using the expressions for $w_{1}(m)$ and $w_{2}(m)$ that we computed in Kurchan et al. 2013. The numerical solution of equation (127) gives

We therefore obtain the schematic phase diagram represented in Fig. 3, which is strongly similar to the one found in Rizzo 2013 for the Ising $p$-spin glass model. For $m>m^{*}$ or $\widehat{\varphi}<\widehat{\varphi}^{*}$, no solution exist below the Gardner line, which therefore delimits the region where the only non-trivial 1RSB solution exists. Instead, for $m<m^{*}$ or $\widehat{\varphi}>\widehat{\varphi}^{*}$ a non-trivial $k$RSB solution with $k>1$ exists below the Gardner line. It is reasonable to expect that this solution will exist in a finite region below the Gardner line. The region of existence of the 2RSB solution should be delimited by a dynamical line $\widehat{\varphi}_{\rm d}^{\rm 2RSB}(m)$ or $m_{\rm d}^{\rm 2RSB}(\widehat{\varphi})$, shifted with respect to the 1RSB dynamical line (see Fig. 3 for a schematic drawing of this line). We will show in Sec. VIII.3 how this line can be defined. Note however that the instability of the replicon mode suggests that the 2RSB solution is also unstable towards 3RSB and so on, until the correct fullRSB solution is found.

VII.2 Asymptotic expression for λ⁡(m)\lambda(m) at large densities on the Gardner line

In this section we compute the asymptotic behavior of $\lambda(m)$ for $m\rightarrow 0$ or $\widehat{\varphi}\rightarrow\infty$ on the Gardner line. The importance of this computation is twofold. First, we want to check that the condition $m<\lambda(m)$ holds for all $m<m^{*}$ on the Gardner line, which implies that a non-trivial 2RSB solution exists for all $\widehat{\varphi}>\widehat{\varphi}^{*}$. Second, we have shown that the breaking point of the 2RSB solution is $m_{1}=\lambda(m)$ on the Gardner line. Actually this is true also if we perform a perturbative fullRSB calculation around the instability line Caltagirone et al. 2011. It follows that the asymptotic computation of $\lambda(m)$ tells us what is the behavior of the breaking point at large densities, which will be useful to investigate the general properties of the $k$RSB solutions.

In the following we rely heavily on the notations and results of (Kurchan et al. 2013, Sec.V and VI). Let us first recall the expression for $\lambda(m)$, which can be written as

where $\widehat{\varphi}_{\rm G}(m)$ is the Gardner line, $\widehat{A}_{\rm G}(m)$ is the 1RSB cage radius on the Gardner line, and

Following the analysis and the notations of (Kurchan et al. 2013, Sec.V and VI), at the instability line we have $\widehat{\varphi}_{\rm G}(m)=1/\mathcal{F}_{m}(\widehat{A}_{\rm G}(m))$, where $\widehat{A}_{\rm G}(m)$ satisfies the equation $2\mathcal{F}_{m}(\widehat{A}_{\rm G}(m))=-\Lambda_{m}(\widehat{A}_{\rm G}(m))$ and

It follows that the expression for $\lambda$ can be put in the form

In (Kurchan et al. 2013, Sec. V D) it was shown that in the limit $m\rightarrow 0$, $\sqrt{\widehat{A}_{\rm G}(m)}\simeq 0.8m$. This means that we can hope to expand the numerator and the denominator in powers of $\sqrt{A_{\rm G}}$.

Let us study the numerator first. It happens that $w_{2}^{(0)}(m)=0$ with very good numerical accuracy, and the equality can be probably demonstrated by a series of integrations by parts (see Kurchan et al. 2013 for a similar computation). The first order term (multiplied by 2 for convenience) can be written as

As it was discussed in Kurchan et al. 2013, the behavior of the integral depends on how the function inside behaves as $\lambda\rightarrow\infty$ for small $m$. We have two possibilities: the integral decays as $\lambda^{\alpha}$ with $\alpha>1$ and in that case the integral is convergent; on the other case, we have a divergent contribution that has to be studied looking at the limit $m\rightarrow 0$. To see which of the two behaviors happens we develop asymptotically the integrand

from which it follows that the integral is finite at $m=0$ and it is given by

Let us now study the denominator. Also in this case the zeroth order term is zero with very good accuracy. The first order term of the denominator is

Also in this case the integral is finite at $m\rightarrow 0$ and it is given by

which shows that $\lambda(m)$ has a finite limit on the Gardner line and therefore the condition $m<\lambda(m)$ holds for all $m<m^{*}$.

VIII The jamming limit of the 2RSB solution

The opposite limit, with respect to the perturbative computation of Sec. VII, in which the problem simplifies a lot is the jamming limit $m\rightarrow 0$. In this limit, pressure diverges and one approaches the jamming point where particles are in contact. This has been investigated in full details at the 1RSB level Parisi and Zamponi 2010; Berthier et al. 2011: at the 1RSB level, in the limit $m\rightarrow 0$ the parameter $\widehat{\gamma}_{1}=2\widehat{\alpha}_{1}$ remains finite, in such a way that the mean square displacement in the glass, $\widehat{\Delta}_{1}=m\widehat{\gamma}_{1}$ vanishes proportionally to $m$, $\widehat{\Delta}_{1}\sim m\sim 1/p$.

In this section we discuss what happens at the 2RSB level. This is interesting because it allows us to determine the endpoint of the 2RSB dynamical line, which corresponds to the 2RSB threshold $\widehat{\varphi}_{\rm th}^{\rm 2RSB}$, see Fig. 3. In general, we can expect (based on the experience accumulated on the spin glass models Mézard et al. 1987) that the 2RSB computation is an extremely good quantitative approximation for the fullRSB result as far as thermodynamic quantities are concerned, therefore the 2RSB threshold should be a very good approximation of the fullRSB one. Moreover, we will encounter here most of the numerical difficulties that will also be relevant for the study of the fullRSB solution.

In Sec. VIII.1 we obtain the expression of the 2RSB entropy and the associated variational equations in the limit $m\rightarrow 0$. In Sec. VIII.2 we show that the variational equations can be solved analytically in a systematic high density expansion. In Sec. VIII.3 we discuss the results of a numerical solution of the variational equations; we show that the numerical results are consistent with the high density expansion, and we discuss how the 2RSB threshold is computed. Finally, in Sec. VIII.4 we summarize the phase diagram obtained from the 2RSB solution.

In order to discuss the behavior of the 2RSB entropy at small $m$, it is convenient to make a change of variables as follows. We eliminate $\widehat{\gamma}_{2}$ and $m_{1}$ in favor of

In terms of the variables $\widehat{\gamma}_{1}$, $\eta$, $\nu$, and using Eqs. (33) we can reconstruct the other parameters as follows:

Furthermore, from the condition that $\widehat{\Delta}_{1}\geq\widehat{\Delta}_{2}\geq 0$ we see that $\eta\in$, while from $1\geq m_{1}\geq m$ we have that $\nu\in[m,1]$. Writing the entropy (56) in terms of these variables, and writing explicitly the convolutions with some simple changes of variables, we obtain

As we will see later, this expression is also convenient to perform a numerical computation of the optimal values of the parameters $\widehat{\gamma}_{1}$, $\eta$, $\nu$.

We now want to check that Eq. (144) has a finite limit $m\rightarrow 0$ if all the other parameters $\widehat{\gamma}_{1},\eta,\nu$ are fixed (i.e. they do not scale with $m$). Note that in this limit both $\eta,\nu\in$ and moreover, according to Eq. (143), $\widehat{\Delta}_{2}\rightarrow 0$ while $\widehat{\Delta}_{1}$ remains finite. The limit $m\rightarrow 0$ of Eq. (144) can be taken easily. The only non-trivial part is the function $I_{2}^{m}(x)$. Using Eq. (115), we have

We see that we obtain a finite limit that corresponds to the complexity at $m=0$ Parisi and Zamponi 2010, and we also conclude that the three parameters $\widehat{\gamma}_{1}$, $\eta$, $\nu$ have finite values at $m=0$. From this, we reach the important physical conclusion that the mean square displacement inside a glass, $\widehat{\Delta}_{2}$, vanishes at jamming, as it should, but at the same time, the mean square displacement $\widehat{\Delta}_{1}$ between different sub-glasses inside a meta-glass remains finite. Hence, sub-glasses inside a meta-glass are microscopically distinct. Finally, we note that $m_{1}$ vanishes proportionally to $m$ because $\nu$ is finite.

VIII.2 The high density limit of the 2RSB solution at m=0m=0

The expression (146) is still quite difficult to handle numerically. Therefore, before discussing the numerical optimization it is convenient to obtain some asymptotic results for large density. We first observe that in the 1RSB case $\widehat{\gamma}_{1}\sim\widehat{\varphi}^{-2}$ and we therefore expect the same scaling also in the 2RSB case. Furthermore, on the Gardner transition line we showed in section VII that $m_{1}=\lambda(m)\rightarrow 0.124$ and that $m\sim 1.98\widehat{\varphi}^{-2}$, hence $\nu=m/m_{1}\sim 16.0\widehat{\varphi}^{-2}$. We therefore seek for a small $\widehat{\gamma}_{1}$ and small $\nu$ expansion of Eq. (146). Because both $\widehat{\gamma}_{1}$ and $\nu$ are of the same order of magnitude, we use $1/\widehat{\varphi}$ as the small expansion parameter. Note that we write $\widehat{\gamma}_{1},\nu\sim{\cal O}(2)$ to indicate that these quantities are of order 2 in $1/\widehat{\varphi}$, and similarly for other quantities.

According to the definition in Eq. (54), the part of the entropy that has to be optimized is

When $\nu$ is small, it is very useful to separate two contributions to this integral as follows:

The term ${\cal I}^{(a)}$ can be easily handled and expanded in a power series of $\nu$ and $\widehat{\gamma}_{1}$. Moreover, the integrand of ${\cal I}^{(n)}$ is well behaved at large $|x|$ (it decays as a Gaussian) for all values of $\nu$, therefore we can expand the integrand.

We can expand ${\cal I}^{(n)}$ as follows

The functions ${\cal I}_{l,k}(\eta)$ are defined by well convergent integrals, hence they are well behaved functions of $\eta$ and they can be differentiated as many times as one wishes by exchanging the derivative with the integration over $x$. In this way, expanding all the variables in powers of $1/\widehat{\varphi}$ as follows,

one can perform a systematic expansion of the entropy in powers of $1/\widehat{\varphi}$ and determine the coefficients $\widehat{\gamma}_{1,k}$, $\nu_{k}$, $\eta_{k}$ by optimizing order by order in inverse density. This computation can be very easily performed with the help of some algebraic manipulation software (we used Mathematica). Below we just give an example of the lowest orders in the computation.

At the lowest order we have, recalling that ${\cal O}(n)$ denotes an order $n$ in $1/\widehat{\varphi}$:

Collecting all together and rearranging the terms in increasing order in the expansion we have

We have now to optimize Eq. (154) order by order. At the leading order we obtain an equation for $\widehat{\gamma}_{1}$ which coincides with the 1RSB one:

We therefore have to look for a solution of the form $\widehat{\gamma}_{1}=\frac{8}{\pi}\widehat{\varphi}^{-2}+\widehat{\gamma}_{1,3}\widehat{\varphi}^{-3}+{\cal O}(4)$. Plugging this in Eq. (154) and expanding we get

which allows to determine $\widehat{\gamma}_{1,3}=128/\pi^{2}$. Finally we look for $\widehat{\gamma}_{1}=\frac{8}{\pi}\widehat{\varphi}^{-2}+(128/\pi^{2})\widehat{\varphi}^{-3}+\widehat{\gamma}_{1,4}\widehat{\varphi}^{-4}+{\cal O}(5)$ and $\nu=\nu_{2}\widehat{\varphi}^{-2}+{\cal O}(3)$, plug this in Eq. (154) and expand to get

This function must be optimized numerically and we obtain $\nu_{2}=5.4226$ and $\eta_{0}=0.6752$. We therefore obtain the asymptotic 2RSB solution for $m=0$ and $\widehat{\varphi}\rightarrow\infty$:

This last result shows that at the 2RSB level the value of $\widehat{\varphi}_{\rm GCP}$, which corresponds to the point where the complexity (equal to the replicated entropy) at $m=0$ vanishes, is slightly reduced with respect to the 1RSB level. However, this happens only at subdominant orders in large $d$, because the dominant orders are defined by a term $\log d$ that comes from the ideal gas term Parisi and Zamponi 2010.

Higher orders in the calculation can be easily obtained by iterating the above procedure. This requires adding more terms in the expansion in Eq. (154) which can be easily automatized with Mathematica. In Fig. 4 we report the results of the calculation done at order 11 in density, which allows to obtain $\widehat{\gamma}_{1}$ to order $\widehat{\varphi}^{-7}$, $\nu$ to order $\widehat{\varphi}^{-6}$ and $\eta$ to order $\widehat{\varphi}^{-4}$. The results are in perfect agreement with a numerical optimization of the 2RSB entropy (146) at $m=0$, that we describe in Sec. VIII.3.

We expect that this strategy to construct a high density expansion (which corresponds to the small cage expansion of Parisi and Zamponi 2010 in the 1RSB case) could be generalized to $m>0$ and to $k$RSB solutions with $k>2$ with a little bit of additional work. However, having tested the accuracy of our numerical optimization code, we do not pursue this strategy further and we turn to the discussion of the numerical results.

VIII.3 Numerical solution of the 2RSB equations at m=0m=0: the 2RSB threshold

We report here results from the full numerical optimization of the 2RSB entropy at $m=0$, given in Eq. (146). The code we used makes explicit use of the decomposition (149), in such a way that ${\cal I}^{(a)}$ is computed easily and ${\cal I}^{(n)}$ is a numerically stable integral (some care should be taken to write the error functions in a numerically stable way). Taking derivatives with respect to $\eta$ and $\widehat{\gamma}_{1}$, we obtain recurrence equations for these quantities, that we do not report because they are the specialization of Eqs. (113) to the case $k=2$ and $m=0$. For each fixed value of density $\widehat{\varphi}$ and breaking point $\nu=1/y_{1}$, these two equations can be solved by iteration to obtain $\eta$, $\widehat{\gamma}_{1}$, and the entropy $s_{\rm 2RSB}$.

We make now a few remarks on the 2RSB equation at $m=0$:

When $\nu=1$ (corresponding to $m_{1}=m$), the 2RSB entropy reduces to the 1RSB one, function of $\widehat{\gamma}_{2}=\widehat{\gamma}_{1}(1-\eta)$.

Similarly, when $\nu\rightarrow 0$ (corresponding to $m_{1}=1$), the 2RSB entropy reduces to the 1RSB one, function of $\widehat{\gamma}_{1}$.

Finally, when $\eta\rightarrow 0$ (corresponding to $\widehat{\gamma}_{1}=\widehat{\gamma}_{2}$), one again recovers the 1RSB entropy.

Let us call $f(\nu,\eta;\widehat{\varphi})=\min_{\widehat{\gamma}_{1}}s_{2RSB}$. This function has to be minimized with respect to $\{\nu,\eta\}\in^{2}$. Based on the considerations above, $f(\nu,\eta;\widehat{\varphi})$ is a constant equal to the 1RSB value of the free energy when $\nu=0$, $\nu=1$ or $\eta=0$, i.e. on three sides of the box $^{2}$. Moreover its derivative in $\eta=0$ is always strictly negative (except for $\nu=0$ and $\nu=1$) as a consequence of the instability of the replicon mode (it is easy to see that the replicon is related to the derivative with respect to $\eta$ of $s_{\rm 2RSB}$. We conclude that a minimum must exist at non trivial values of $\nu$ and $\eta$ (remember that as usual in replica computations, the entropy should be minimized and not maximized Mézard et al. 1987). However, it is also important to remark that for some values of the parameters $\{\nu,\eta\}$ the solution for $\widehat{\gamma}_{1}$ might not exist formally corresponding to $\widehat{\gamma}_{1}=\infty$, which corresponds to losing the non-trivial 2RSB solution to a trivial solution where all replicas are uncorrelated.