\documentclass[fleqn]{article} \usepackage{fullpage,amsmath,amssymb,latexsym,graphicx} \usepackage{appendix} \usepackage[dvipsnames]{xcolor} \usepackage[ colorlinks=true, urlcolor=MidnightBlue, citecolor=MidnightBlue, filecolor=MidnightBlue, linkcolor=MidnightBlue ]{hyperref} % ref and cite links with pretty colors \begin{document} \title{Full solution of the Kac--Rice problem for mean-field models.\\ or Full solution for the counting of saddles of mean-field glass models} \author{Jaron Kent-Dobias \& Jorge Kurchan} \maketitle \begin{abstract} We derive the general solution for the computation of saddle points of mean-field complex landscapes. The solution incorporates Parisi's solution for the ground state, as it should. \end{abstract} \section{Introduction} The computation of the number of metastable states of mean field spin glasses goes back to the beginning of the field. Over forty years ago, Bray and Moore \cite{Bray_1980_Metastable} attempted the first calculation for the Sherrington--Kirkpatrick model, in a paper remarkable for being one of the first applications of a replica symmetry breaking scheme. As became clear when the actual ground-state of the model was computed by Parisi \cite{Parisi_1979_Infinite} with a different scheme, the Bray--Moore result was not exact, and in fact the problem has been open ever since. To this date the program of computing the number of saddles of a mean-field glass has been only carried out for a small subset of models, including most notably the (pure) $p$-spin model ($p>2$) \cite{Rieger_1992_The, Crisanti_1995_Thouless-Anderson-Palmer}. In a parallel development, it has evolved into an active field in probability theory \cite{Auffinger_2012_Random, Auffinger_2013_Complexity, BenArous_2019_Geometry} In this paper we present what we argue is the general replica ansatz for the computation of the number of saddles of generic mean-field models, which we expect to include the Sherrington--Kirkpatrick model. It reproduces the Parisi result in the limit of small temperature for the lowest states, as it should. To understand the importance of this computation, consider the following situation. When one solves the problem of spheres in large dimensions, one finds that there is a transition at a given temperature to a one-step one step symmetry breaking (1RSB) phase at a Kauzmann temperature, and, at a lower temperature, another transition to a full RSB phase (see \cite{gross1985mean,gardner1985spin}, the `Gardner ' phase \cite{charbonneau2014fractal}. Now, this transition involves the lowest, equilibrium states. Because they are obviously unreachable at any reasonable timescale, an often addressed question to ask is "what is the Gardner transition line for higher than equilibrium energy-densities"? (see, for a review \cite{berthier2019gardner}) For example, when studying `jamming' at zero temperature, the question is posed as to "on what side of the 1RSB-FRS transition are the high energy (or low density) states reachable dynamically. Posed in this way, such a question does not have a clear definition. In the present paper we give a concrete strategy to define unambiguously such an issue: we consider the local energy minima at a given energy and study their number and other properties: the solution involves a replica-symmetry breaking scheme that is well-defined, and corresponds directly to the topological characteristics of those minima. \section{The model} Here we consider, for definiteness, the mixed $p$-spin model, \begin{equation} H(s)=-\sum_p\frac1{p!}\sum_{i_1\cdots i_p}J^{(p)}_{i_1\cdots i_p}s_{i_1}\cdots s_{i_p} \end{equation} for $\overline{(J^{(p)})^2}=a_pp!/2N^{p-1}$ Then \begin{equation} \overline{H(s_1)H(s_2)}=Nf\left(\frac{s_1\cdot s_2}N\right) \end{equation} for \begin{equation} f(q)=\frac12\sum_pa_pq^p \end{equation} or, more generally, the `Toy Model' of M\'ezard and Parisi \cite{Mezard_1992_Manifolds} which involve non-polynomial forms for $f(q)$. These may be thought of as a model of generic Gaussian functions on the sphere. To constrain the model to the sphere, we use a Lagrange multiplier $\mu$, with the total energy being \begin{equation} H(s)+\frac\mu2(s\cdot s-N) \end{equation} At any critical point, the gradient and Hessian are \begin{align} \nabla H=\partial H+\mu s && \operatorname{Hess}H=\partial\partial H+\mu I \end{align} The important observation was made by Bray and Dean \cite{Bray_2007_Statistics} that gradient and Hessian are independent for random Gaussian disorder. The average over disorder breaks into a product of two independent averages, one for the gradient factor and one for any function of the Hessian, in particular its number of negative eigenvalues, the index $\mathcal I$ of the saddle (see Fyodorov \cite{Fyodorov_2007_Replica} for a detailed discussion).. \section{Equilibrium} Here we review the equilibrium solution. \cite{Crisanti_1992_The, Crisanti_1993_The, Crisanti_2004_Spherical, Crisanti_2006_Spherical} The free energy averaged over disorder is \begin{equation} \beta F = - \overline{\ln \int ds \; e^{-\beta H(s)}} \end{equation} Computing the logarithm as the limit of $n \rightarrow 0$ replicas is standard, and it take the form \begin{equation} \beta F=-1-\log2\pi-\frac12\lim_{n\to0}\frac1n\left(\beta^2\sum_{ab}^nf(Q_{ab})+\log\det Q\right) \end{equation} which must be extremized over the matrix $Q$. When the solution is a Parisi matrix, this can also be written in a functional form. If $P(q)$ is the probability distribution for elements $q$ in a row of the matrix, then \begin{equation} \chi(q)=\int_1^qdq'\,\int_0^{q'}dq''\,P(q'') \end{equation} Since it is the double integral of a probability distribution, $\chi$ must be concave, monotonically decreasing and have $\chi(1)=0$, $\chi'(1)=1$. The free energy can be written as a functional over $\chi$ as \begin{equation} \beta F=-1-\log2\pi-\frac12\int dq\,\left(\beta^2f''(q)\chi(q)+\frac1{\chi(q)}\right) \end{equation} We insert a $k$ step RSB ansatz (the steps are standard, and are reviewed in Appendix A), and obtain $q_0=0$ \begin{equation} \begin{aligned} \beta F= -1-\log2\pi -\frac12\left\{\beta^2\left(f(1)+\sum_{i=0}^k(x_i-x_{i+1})f(q_i)\right) +\frac1{x_1}\log\left[ 1+\sum_{i=1}^{k}(x_i-x_{i+1})q_i \right]\right.\\ \left.+\sum_{j=1}^k(x_{j+1}^{-1}-x_j^{-1})\log\left[ 1+\sum_{i=j+1}^{k}(x_i-x_{i+1})q_i-x_{j+1}q_j \right] \right\} \end{aligned} \end{equation} The zero temperature limit is most easily obtained by putting $x_i=\tilde x_ix_k$ and $x_k=\tilde\beta/\beta$, $q_k=1-z/\beta$, which ensures the $\tilde x_i$, $\tilde\beta$, and $z$ have nontrivial limits. Inserting the ansatz and taking the limit carefully treating the $k$th term in each sum separately from the rest, we get \begin{equation} \begin{aligned} \lim_{\beta\to\infty}\tilde\beta F= -\frac12z\tilde\beta f'(1)-\frac12\left\{\tilde\beta^2\left(f(1)+\sum_{i=0}^{k-1}(\tilde x_i-\tilde x_{i+1})f(q_i)\right) +\frac1{\tilde x_1}\log\left[ \tilde\beta z^{-1}\left(1+\sum_{i=1}^{k-1}(\tilde x_i-\tilde x_{i+1})q_i\right)+1 \right]\right.\\ \left.+\sum_{j=1}^{k-1}(\tilde x_{j+1}^{-1}-\tilde x_j^{-1})\log\left[ \tilde\beta z^{-1}\left(1+\sum_{i=j+1}^{k-1}(\tilde x_i-\tilde x_{i+1})q_i-\tilde x_{j+1}q_j\right)+1 \right] \right\} \end{aligned} \end{equation} This is a $k-1$ RSB ansatz with all eigenvalues scaled by $\tilde\beta z^{-1}$ and shifted by 1, with effective temperature $\tilde\beta$, and an extra term. This can be seen more clearly by rewriting the result in terms of the matrix $\tilde Q$, a $(k-1)$-RSB Parisi matrix built from the $\tilde x_1,\ldots,\tilde x_{k-1}$ and $q_1,\ldots,q_{k-1}$, which gives \begin{equation} \label{eq:ground.state.free.energy} \lim_{\beta\to\infty}\tilde\beta F =-\frac12z\tilde\beta f'(1)-\frac12\lim_{n\to0}\frac1n\left( \tilde\beta^2\sum_{ab}^nf(\tilde Q_{ab})+\log\det(\tilde\beta z^{-1}\tilde Q+I) \right) \end{equation} In the continuum case, this is \begin{equation} \label{eq:ground.state.free.energy.cont} \lim_{\beta\to\infty}\tilde\beta F =-\frac12z\tilde\beta f'(1)-\frac12\int dq\left( \tilde\beta^2f''(q)\tilde\chi(q)+\frac1{\tilde\chi(q)+\tilde\beta z^{-1}} \right) \end{equation} The zero temperature limit of the free energy loses one level of replica symmetry breaking. Physically, this is a result of the fact that in $k$-RSB, $q_k$ gives the overlap within a state, e.g., within the basin of a well inside the energy landscape. At zero temperature, the measure is completely localized on the bottom of the well, and therefore the overlap with each state becomes one. We will see that the complexity of low-energy stationary points in Kac--Rice computation is also given by a $(k-1)$-RSB anstaz. Heuristically, this is because each stationary point also has no width and therefore overlap one with itself. \section{Kac--Rice calculation} \cite{Auffinger_2012_Random, BenArous_2019_Geometry} The stationary points of a function can be counted using the Kac--Rice formula, which integrates a over the function's domain a $\delta$-function containing the gradient multiplied by the absolute value of the determinant \cite{Rice_1939_The, Kac_1943_On}. In addition, we insert a $\delta$-function fixing the energy density $E$, and another delta function to count the number of saddles with trace of the Hessian $=\mu^*$. The latter will give us everything we need to characterize the saddles, as we shall see later \begin{equation} \begin{aligned} \mathcal N(E, \mu^*) &=\int ds\, d\mu\,\delta\big(\tfrac12(\|s\|^2-N)\big)\,\delta\big(\nabla H(s,\mu)\big)\,\big|\det\operatorname{Hess}H(s,\mu)\big| \\ &\hspace{10pc}\times\delta\big(NE-H(s)\big)\delta\big(N\mu^*-\operatorname{Tr}\operatorname{Hess}H(s,\mu)\big) \end{aligned} \end{equation} This number will typically be exponential in $N$. In order to find typical counts when disorder is averaged, we will want to average its logarithm instead, which is known as the averaged complexity: \begin{equation} \Sigma(E,\mu^*)=\lim_{N\to\infty}\frac1N\overline{\log\mathcal N(E, \mu^*}) \end{equation} If one averages over $\mathcal N$ and afterward takes its logarithm, one arrives at the so-called {\em annealed} complexity \begin{equation} \Sigma_\mathrm a(E,\mu^*) =\lim_{N\to\infty}\frac1N\log\overline{\mathcal N(E,\mu^*)} \end{equation} This has been previously computed for the mixed $p$-spin models \cite{BenArous_2019_Geometry}, with the result %\begin{equation} % \begin{aligned} %\Sigma_\mathrm a(E,\mu) %=-\frac{E^2(f'(1)+f''(1))+2E\mu f'(1)+f(1)\mu^2}{2f(1)(f'(1)+f''(1))-2f'(1)^2}-\frac12\log f'(1)\\ % +\operatorname{Re}\left[\frac\mu{\mu+\sqrt{\mu^2-4f''(1)}} % -\log\left(\frac1{2f''(1)}\left(\mu-\sqrt{\mu^2-4f''(1)}\right)\right) % \right] % \end{aligned} %\end{equation} The annealed complexity is known to equal the actual (quenched) complexity in circumstances where there is at most one level of RSB. This is the case for the pure $p$-spin models, or for mixed models where $1/\sqrt{f''(q)}$ is a convex function. However, it fails dramatically for models with higher replica symmetry breaking. For instance, when $f(q)=\frac12(q^2+\frac1{16}q^4)$, the anneal complexity predicts that minima vanish well before the dominant saddles, a contradiction for any bounded function, as seen in Fig.~\ref{fig:frsb.complexity}. A sometimes more illuminating quantity to consider is the Legendre transform $G$ of the complexity: \begin{equation} e^{NG(\hat \beta, \mu^*)} = \int dE \; e^{ -\hat \beta E +\Sigma(\hat \beta, \mu^*)} \end{equation} There will be a critical value $\hat \beta_c$ beyond which the complexity is zero: above this value the measure is split between the lowest $O(1)$ energy states. We shall not study here this regime that interpolates between the dynamically relevant and the equilibrium states, but just mention that it is an interesting object of study. \subsection{The replicated problem} In order to average the complexity over disorder properly, the logarithm must be dealt with. We use the standard replica trick, writing \begin{equation} \begin{aligned} \log\mathcal N(E,\mu^*) &=\lim_{n\to0}\frac\partial{\partial n}\mathcal N^n(E,\mu^*) \\ &=\lim_{n\to0}\frac\partial{\partial n}\int\prod_a^n ds_a\,d\mu_a\, \delta\big(\tfrac12(\|s_a\|^2-N)\big)\,\delta\big(\nabla H(s_a,\mu_a)\big)\,\big|\det\operatorname{Hess}H(s_a,\mu_a)\big| \\ &\hspace{13pc} \times\delta\big(NE-H(s_a)\big)\delta\big(N\mu^*-\operatorname{Tr}\operatorname{Hess}H(s_a,\mu_a)\big) \end{aligned} \end{equation} The replicated Kac--Rice formula was introduced by Ros et al.~\cite{Ros_2019_Complex}, and its effective action for the mixed $p$-spin model has previously been computed by Folena et al.~\cite{Folena_2020_Rethinking}. Here we review the derivation. In practice, we are therefore able to write \begin{equation} \begin{aligned} \Sigma(E, \mu^*) &=\lim_{N\to\infty}\frac1N\lim_{n\to0}\frac\partial{\partial n}\int\left(\prod_a^nds_a\,d\mu_a\right)\, \overline{\prod_a^n \delta\big(\tfrac12(\|s_a\|^2-N)\big)\,\delta\big(\nabla H(s_a,\mu_a)\big)\delta(NE-H(s_a))}\\ &\hspace{10pc} \times \overline{\prod_a^n |\det\operatorname{Hess}(s_a,\mu_a)|\,\delta\big(N\mu^*-\operatorname{Tr}\operatorname{Hess}H(s_a,\mu_a)\big)} \end{aligned} \end{equation} \subsubsection{The Hessian factors} The spectrum of the Hessian matrix $\partial\partial H$ is uncorrelated from the gradient. In the large $N$ limit for almost every point and realization of disorder a GOE matrix with variance \begin{equation} \overline{(\partial_i\partial_jH)^2}=\frac1Nf''(1)\delta_{ij} \end{equation} and therefore its spectrum is given by the Wigner semicircle with radius $\sqrt{4f''(1)}$, or \begin{equation} \rho(\lambda)=\frac1{2\pi f''(1)}\sqrt{4f''(1)-\lambda^2} \end{equation} and the spectrum of $\operatorname{Hess}H$ is this shifted by $\mu$, or $\rho(\lambda+\mu)$. The parameter $\mu$ thus fixes the spectrum of the Hessian, when $\mu$ is taken to be within the range $\pm\sqrt{4f''(1)}=\pm\mu_m$, the critical (or in fact, any) points have index density \begin{equation} \mathcal I(\mu)=\int_0^\infty d\lambda\,\rho(\lambda+\mu) =N\left\{\frac12-\frac1\pi\left[ \arctan\left(\frac\mu{\sqrt{\mu_m^2-\mu^2}}\right) +\frac\mu{\mu_m^2}\sqrt{\mu_m^2-\mu^2} \right] \right\} \end{equation} When $\mu>\mu_m$, the critical points are minima whose sloppiest eigenvalue is $\mu-\mu_m$. To largest order in $N$, the average over the product of determinants factorizes into the product of averages, each of which is given by the same expression depending only on $\mu$: \begin{equation} \begin{aligned} \overline{\prod_a^n |\det\operatorname{Hess}(s_a,\mu_a)|\,\delta\big(N\mu^*-\operatorname{Tr}\operatorname{Hess}H(s_a,\mu_a)\big)} \rightarrow e^{nN{\cal D}(\mu^*)}\prod_a^n\delta(\mu_a-\mu^*) \end{aligned} \end{equation} with \begin{equation} \begin{aligned} \mathcal D(\mu) &=\frac1N\overline{\log|\det\operatorname{Hess}H(s,\mu)|} =\int d\lambda\,\rho(\lambda+\mu)\log|\lambda| \\ &=\operatorname{Re}\left\{ \frac12\left(1+\frac\mu{2f''(1)}\left(\mu-\sqrt{\mu^2-4f''(1)}\right)\right) -\log\left(\frac1{2f''(1)}\left(\mu-\sqrt{\mu^2-4f''(1)}\right)\right) \right\} \end{aligned} \end{equation} What we have described is the {\em typical} spectrum for given $\mu$. What about the deviations of the spectrum -- we are particularly interested in the number of negative eigenvalues -- at given $\mu$. The result is well known qualitatively: there are two possibilities:\\ $\bullet$ For $|\mu|>\mu_m$ there is the possibility of a finite number of eigenvalues of the second derivative matrix \subsubsection{The gradient factors} The $\delta$-functions are treated by writing them in the Fourier basis, introducing auxiliary fields $\hat s_a$ and $\hat\beta$, \begin{equation} \prod_a^n \delta\big(\tfrac12(\|s_a\|^2-N)\big)\,\delta\big(\nabla H(s_a,\mu^*)\big)\delta(NE-H(s_a)) =\int\frac{d\hat\mu}{2\pi}\,\frac{d\hat\beta}{2\pi}\prod_a^n\frac{d\hat s_a}{2\pi} e^{\frac12\hat\mu(\|s_a\|^2-N)+\hat\beta(NE-H(s_a))+i\hat s_a\cdot(\partial H(s_a)+\mu^*s_a)} \end{equation} $\hat \beta$ is a parameter conjugate to the state energies, i.e. playing the role of an inverse temperature for the metastable states. The average over disorder can now be taken, and since everything is Gaussian it gives \begin{equation} \begin{aligned} \overline{ \exp\left\{ \sum_a^n(i\hat s_a\cdot\partial_a-\hat\beta)H(s_a) \right\} } &=\exp\left\{ \frac12\sum_{ab}^n (i\hat s_a\cdot\partial_a-\hat\beta) (i\hat s_b\cdot\partial_b-\hat\beta) \overline{H(s_a)H(s_b)} \right\} \\ &=\exp\left\{ \frac N2\sum_{ab}^n (i\hat s_a\cdot\partial_a-\hat\beta) (i\hat s_b\cdot\partial_b-\hat\beta) f\left(\frac{s_a\cdot s_b}N\right) \right\} \\ &\hspace{-14em}=\exp\left\{ \frac N2\sum_{ab}^n \left[ \hat\beta^2f\left(\frac{s_a\cdot s_b}N\right) -2i\hat\beta\frac{\hat s_a\cdot s_b}Nf'\left(\frac{s_a\cdot s_b}N\right) -\frac{\hat s_a\cdot \hat s_b}Nf'\left(\frac{s_a\cdot s_b}N\right) +\left(i\frac{\hat s_a\cdot s_b}N\right)^2f''\left(\frac{s_a\cdot s_b}N\right) \right] \right\} \end{aligned} \end{equation} We introduce new fields \begin{align} \label{eq:fields} C_{ab}=\frac1Ns_a\cdot s_b && R_{ab}=-i\frac1N\hat s_a\cdot s_b && D_{ab}=\frac1N\hat s_a\cdot\hat s_b \end{align} Their physical meaning is explained in \S\ref{sec:interpretation}. By substituting these parameters into the expressions above and then making a change of variables in the integration from $s_a$ and $\hat s_a$ to these three matrices, we arrive at the form for the complexity \begin{equation} \begin{aligned} &\Sigma(E,\mu^*) =\mathcal D(\mu^*)+\hat\beta E-\frac12\hat\mu+\\ &\lim_{n\to0}\frac1n\left( \frac12\hat\mu\operatorname{Tr}C-\mu^*\operatorname{Tr}R +\frac12\sum_{ab}\left[ \hat\beta^2f(C_{ab})+(2\hat\beta R_{ab}-D_{ab})f'(C_{ab}) +R_{ab}^2f''(C_{ab}) \right] +\frac12\log\det\begin{bmatrix}C&iR\\iR&D\end{bmatrix} \right) \end{aligned} \end{equation} where $\hat\mu$, $\hat\beta$, $C$, $R$ and $D$ must be evaluated at extrema of this expression. The same information is contained, and better expressed in its Legendre transform \begin{equation} \begin{aligned} &G(\hat \beta,\mu^*) =\mathcal D(\mu^*)-\frac12\hat\mu\\ &\lim_{n\to0}\frac1n\left( \frac12\hat\mu\operatorname{Tr}C-\mu^*\operatorname{Tr}R +\frac12\sum_{ab}\left[ \hat\beta^2f(C_{ab})+(2\hat\beta R_{ab}-D_{ab})f'(C_{ab}) +R_{ab}^2f''(C_{ab}) \right] +\frac12\log\det\begin{bmatrix}C&iR\\iR&D\end{bmatrix} \right) \end{aligned} \end{equation} Denoting $r_d \equiv \frac 1 n {\mbox Tr} R$, we have the double Legendre transform $K(\hat \beta, r_d)$: \begin{equation} e^{N K(\hat \beta, r_d)} =\int \, d\mu^* \,dE \, e^{N\left\{\Sigma(E,\mu^*)+r_d\mu^* -\hat\beta E -\mathcal D(\mu^*)\right\}} \end{equation} given by \begin{equation} \begin{aligned} &K(\hat \beta,r_d) = \lim_{n\to0}\frac1n\left( \frac{\hat\mu}2\operatorname{Tr}(C-I) +\frac12\sum_{ab}\left[ \hat\beta^2f(C_{ab})+(2\hat\beta R_{ab}-D_{ab})f'(C_{ab}) +R_{ab}^2f''(C_{ab}) \right] +\frac12\log\det\begin{bmatrix}C&iR\\iR&D\end{bmatrix} \right) \end{aligned} \end{equation} $r_d$ is conjugate to $\mu^*$ and through it to the index density, while $\hat \beta$ plays the role of an inverse temperature conjugate to the complexity, that has been used since the beginning of the spin-glass field. In this way $K(\hat \beta,r_d)$ contains all the information about saddle densities. \section{Replica ansatz} Based on previous work on the SK model and the equilibrium solution of the spherical model, we expect $C$, and $R$ and $D$ to be hierarchical matrices, i.e., to follow Parisi's scheme. In the end, when the limit of $n\to0$ is taken, each can be represented in the canonical way by its diagonal and a continuous function on the domain $[0,1]$ which parameterizes each of its rows, with \begin{align} C\;\leftrightarrow\;[c_d, c(x)] && R\;\leftrightarrow\;[r_d, r(x)] && D\;\leftrightarrow \;d_d, d(x)] \end{align} This assumption immediately simplifies the extremal conditions, since hierarchical matrices commute and are closed under matrix products and Hadamard products. The extremal conditions are \begin{align} 0&=\frac{\partial\Sigma}{\partial\hat\mu} =\frac12(c_d-1) \\ 0&=\frac{\partial\Sigma}{\partial\hat\beta} =E+\lim_{n\to0}\frac1n\sum_{ab}\left[\hat\beta f(C_{ab})+R_{ab}f'(C_{ab})\right] \label{eq:cond.b} \\ 0&=\frac{\partial\Sigma}{\partial C} =\frac12\left[ \hat\mu I+\hat\beta^2f'(C)+(2\hat\beta R-D)\odot f''(C)+R\odot R\odot f'''(C) +(CD+R^2)^{-1}D \right] \label{eq:cond.q} \\ 0&=\frac{\partial\Sigma}{\partial R} =-\mu^* I+\hat\beta f'(C)+R\odot f''(C) +(CD+R^2)^{-1}R \label{eq:cond.r} \\ 0&=\frac{\partial\Sigma}{\partial D} =-\frac12f'(C) +\frac12(CD+R^2)^{-1}C \label{eq:cond.d} \end{align} where $\odot$ denotes the Hadamard product, or the componentwise product. Equation \eqref{eq:cond.d} implies that \begin{equation} \label{eq:D.solution} D=f'(C)^{-1}-RC^{-1}R \end{equation} In addition to these equations, one is also often interested in maximizing the complexity as a function of $\mu^*$, to find the dominant or most common type of stationary points. These are given by the condition \begin{equation} \label{eq:cond.mu} 0=\frac{\partial\Sigma}{\partial\mu^*} =\mathcal D'(\mu^*)-r_d \end{equation} Since $\mathcal D(\mu^*)$ is effectively a piecewise function, with different forms for $\mu^*$ greater or less than $\mu_m$, there are two regimes. When $\mu^*>\mu_m$, the critical points are minima, and \eqref{eq:cond.mu} implies \begin{equation} \label{eq:mu.minima} \mu^*=\frac1{r_d}+r_df''(1) \end{equation} When $\mu^*<\mu_m$, they are saddles, and \begin{equation} \label{eq:mu.saddles} \mu^*=2f''(1)r_d \end{equation} We will find it often useful to have the extremal conditions in a form without matrix inverses, both for numerics at finite $k$-RSB and for expanding in the continuous case. By simple manipulations, the matrix equations can be written as \begin{align} 0&=\left[\hat\beta^2f'(C)+(2\hat\beta R-D)\odot f''(C)+R\odot R\odot f'''(C)+\hat\mu I\right]C+f'(C)D \\ 0&=\left[\hat\beta f'(C)+R\odot f''(C)-\mu^*I\right]C+f'(C)R \\ 0&=C-f'(C)(CD+R^2) \end{align} The right-hand side of each of these equations is also a hierarchical matrix, since products, Hadamard products, and sums of hierarchical matrices are such. The equations for the continuous case are found by using the mappings to functions $c(x)$, $r(x)$ and $d(x)$, then carrying through the appropriate operations. \section{Supersymmetric solution} The Kac--Rice problem has an approximate supersymmetry, which is found when the absolute value of the determinant is neglected. When this is done, the determinant can be represented by an integral over Grassmann variables, which yields a complexity depending on `bosons' and `fermions' that share the supersymmetry. The Ward identities associated with the supersymmetry imply that $D=\hat\beta R$ \cite{Annibale_2003_The}. Under which conditions can this relationship be expected to hold? We find that their applicability is limited. Any result of supersymmetry can only be valid when the symmetry itself is valid, which means the determinant must be positive. This is only guaranteed for minima, which have $\mu^*>\mu_m$. Moreover, this identity heavily constrains the form that the rest of the solution can take. Assuming the supersymmetry holds, \eqref{eq:cond.q} implies \begin{equation} 0=\hat\mu I+\hat\beta^2f'(C)+\hat\beta R\odot f''(C)+R\odot R\odot f'''(C)+\hat\beta(CD+R^2)^{-1}R \end{equation} Substituting \eqref{eq:cond.r} for the factor $(CD+R^2)^{-1}R$, we find substantial cancellation, and finally \begin{equation} \label{eq:R.diagonal} 0=(\hat\mu+\mu^*)I+R\odot R\odot f'''(C) \end{equation} If $C$ has a nontrivial off-diagonal structure and supersymmetry holds, then the off-diagonal of $R$ must vanish, and therefore $R=r_dI$. Therefore, a supersymmetric ansatz is equivalent to a \emph{diagonal} ansatz. Supersymmetry has further implications. Equations \eqref{eq:cond.r} and \eqref{eq:cond.d} can be combined to find \begin{equation} I=R\left[\mu^* I-R\odot f''(C)\right]+(D-\hat\beta R)f'(C) \end{equation} Assuming the supersymmetry holds implies that \begin{equation} I=R\left[\mu^* I-R\odot f''(C)\right] \end{equation} Understanding that $R$ is diagonal, this implies \begin{equation} \mu^*=\frac1{r_d}+r_df''(1) \end{equation} which is precisely the condition \eqref{eq:mu.minima}. Therefore, \emph{the supersymmetric solution only counts dominant minima.} Inserting the supersymmetric ansatz $D=\hat\beta R$ and $R=r_dI$, one gets \begin{equation} \label{eq:diagonal.action} \begin{aligned} \Sigma(E,\mu^*) =\mathcal D(\mu^*) + \hat\beta E-\mu^* r_d +\frac12\hat\beta r_df'(1)+\frac12r_d^2f''(1)+\frac12\log r_d^2 \\ +\frac12\lim_{n\to0}\frac1n\left(\hat\beta^2\sum_{ab}f(C_{ab})+\log\det((\hat\beta/r_d)C+I)\right) \end{aligned} \end{equation} From here, it is straightforward to see that the complexity vanishes at the ground state energy. First, in the ground state minima will dominate (even if they are marginal), so we may assume \eqref{eq:mu.minima}. Then, taking $\Sigma(E_0,\mu^*)=0$, gives \begin{equation} \hat\beta E_0 =-\frac12r_d\hat\beta f'(1)-\frac12\lim_{n\to0}\frac1n\left( \hat\beta^2\sum_{ab}^nf(C_{ab}) +\log\det(\hat\beta r_d^{-1} C+I) \right) \end{equation} which is precisely \eqref{eq:ground.state.free.energy} with $r_d=z$, $\hat\beta=\tilde\beta$, and $C=\tilde Q$. {\em We arrive at one of the main results of our paper: a $(k-1)$-RSB ansatz in Kac--Rice will predict the correct ground state energy for a model whose equilibrium state at small temperatures is $k$-RSB } Moreover, there is an exact correspondance between the saddle parameters of each. If the equilibrium is given by a Parisi matrix with parameters $x_1,\ldots,x_k$ and $q_1,\ldots,q_k$, then the parameters $\hat\beta$, $r_d$, $d_d$, $\tilde x_1,\ldots,\tilde x_{k-1}$, and $\tilde c_1,\ldots,\tilde c_{k-1}$ for the complexity in the ground state are \begin{align} \hat\beta=\lim_{\beta\to\infty}\beta x_k && \tilde x_i=\lim_{\beta\to\infty}\frac{x_i}{x_k} && \tilde c_i=\lim_{\beta\to\infty}q_i && r_d=\lim_{\beta\to\infty}\beta(1-q_k) && d_d=\hat\beta r_d \end{align} \section{Full replica symmetry breaking} This reasoning applies equally well to FRSB systems. \begin{align} 0&=\hat\mu c(x)+[(\hat\beta^2f'(c)+(2\hat\beta r-d)f''(c)+r^2f'''(c))\ast c](x)+(f'(c)\ast d)(x) \label{eq:extremum.c} \\ 0&=-\mu^* c(x)+[(\hat\beta f'(c)+rf''(c))\ast c](x)+(f'(c)\ast r)(x) \label{eq:extremum.r} \\ 0&=c(x)-\big(f'(c)\ast(c\ast d+r\ast r)\big)(x) \label{eq:extremum.d} \end{align} where the product of two hierarchical matrix results in \begin{align} (a\ast b)_d&=a_db_d-\langle ab\rangle \\ (a\ast b)(x)&=(b_d-\langle b\rangle)a(x)+(a_d-\langle a\rangle)b(x) -\int_0^xdy\,\big( a(x)-a(y) \big)\big( b(x)-b(y) \big) \end{align} and \begin{equation} \langle a\rangle=\int_0^1dx\,a(x) \end{equation} \subsection{Supersymmetric complexity} Using standard manipulations (Appendix B), one finds also a continuous version of the supersymmetric complexity \begin{equation} \label{eq:functional.action} \Sigma(E,\mu^*) =\mathcal D(\mu^*) + \hat\beta E-\mu^* r_d +\frac12\hat\beta r_df'(1)+\frac12r_d^2f''(1)+\frac12\log r_d^2 +\frac12\int_0^1dq\,\left( \hat\beta^2f''(q)\chi(q)+\frac1{\chi(q)+r_d/\hat\beta} \right) \end{equation} where $\chi(q)=\int_1^qdq'\int_0^{q'}dq''\,P(q)$, as in the equilibrium case. Though the supersymmetric solution leads to a nice tractable expression, it turn out to be useful only at one point of interest: the ground state. Indeed, we know from the equilibrium that in the ground state $\chi$ is continuous in the whole range of $q$. Therefore, the saddle solution found by extremizing \begin{equation} 0=\frac{\delta\Sigma}{\delta\chi(q)}=\frac12\hat\beta^2f''(q)-\frac12\frac1{(\chi(q)+r_d/\hat\beta)^2} \end{equation} given by \begin{equation} \chi(q)=\frac1{\hat\beta}\left(f''(q)^{-1/2}-r_d\right) \end{equation} is correct. This is only correct if it satisfies the boundary condition $\chi(1)=0$, which requires $r_d=f''(1)^{-1/2}$. This in turn implies $\mu^*=\frac1{r_d}+f''(1)r_d=\sqrt{4f''(1)}=\mu_m$. Therefore, the FRSB ground state is exactly marginal! It is straightforward to check that these conditions are indeed a saddle of the complexity. This has several implications. First, other than the ground state, there are \emph{no} energies at which minima are most numerous; saddles always dominante. As we will see, stable minima are numerous at energies above the ground state, but these vanish at the ground state. \subsection{Expansion near the transition} Working with the general equations in their continuum form away from the supersymmetric solution is not generally tractable. However, there is another point where they can be treated analytically: near the onset of replica symmetry breaking. Here, the off-diagonal components of $C$, $R$, and $D$ are expected to be small. In particular, we expect the function $c$, $r$, and $d$ to approach zero at the transition, and moreover take the form \begin{equation} c(x)=\begin{cases}\bar cx&x\leq x_\mathrm{max}\\\bar cx_\mathrm{max}&\text{otherwise}\end{cases} \end{equation} with $x_\mathrm{max}$ vanishing at the transition, with the slopes $\bar c$, $\bar r$, and $\bar d$ remaining nonzero. This ansatz is informed both by the experience of the equilibrium solution, and by empirical observation within the numerics. Given this ansatz, we take the equations \eqref{eq:extremum.c}, \eqref{eq:extremum.r}, and \eqref{eq:extremum.d}, which are true for any $x$, and integrate them over $x$. We then expand the result about small $x_\mathrm{max}$ to quadratic order in $x_\mathrm{max}$. Equation \eqref{eq:extremum.r} depends linearly on $\bar r$ to all orders, and therefore $\bar r$ can be found in terms of $\bar c$, yielding \begin{equation} \begin{aligned} \frac{\bar r}{\bar c} &= -\hat\beta-\frac1{f'(1)+f''(0)}\left(r_d(f''(0)+f''(1))-\mu^*\right) +\frac12\frac{\bar c}{f'(1)+f''(0)}\left\{ -\hat\beta(4f''(0)-f'''(0))\right.\\ &\left.+\frac1{f'(1)+f''(0)}\left[ 8(\mu-r_d(f''(0)+f''(1)))f''(0) -(2\mu+r_d(f'(1)-2f''(1)-f''(0))f'''(0)) \right] \right\}x_\textrm{max}+O(x_\textrm{max}^2) \end{aligned} \end{equation} Likewise, \eqref{eq:extremum.d} depends linearly on $\bar d$ to all orders, and can be solved using this form for $\bar r$ to give \begin{equation} \begin{aligned} \frac{\bar d}{\bar c} &=-2r_d\frac{\bar r}{\bar c}-\frac1{f'(1)}(r_d^2f''(0)+d_d(f'(1)+f''(0))-1) +\left\{ 2\frac{\bar r^2}{\bar c}-2\bar c\bar b-4r_d\bar r\right.\\ &\left. -\frac{\bar c}{2f'(1)^2}\left[ 4(f'(1)+f''(0))((d_d+r_d^2)f''(0)-1)+(d_d+r_d^2)f'(1)f'''(0) \right] \right\}x_\mathrm{max}+O(x_\mathrm{max}^2) \end{aligned} \end{equation} Finally, the equations for $\bar c$ at first order imply that either $\bar c$ vanishes as $x_\mathrm{max}$ to zero, or that the linear coefficient vanishes as $x_\mathrm{max}$ to zero. Using $\hat\mu$, $\hat\beta$, $r_d$, and $d_d$ from the annealed solution, this coefficient vanishes when \begin{equation} \mu^* =\pm\frac{(f'(1)+f''(0))(f'(1)^2-f(1)(f'(1)+f''(1)))}{(2f(1)-f'(1))f'(1)f''(0)^{-1/2}} -\frac{f''(1)-f'(1)}{f'(1)-2f(1)}\epsilon \end{equation} We expect that this is the line of stability for the replica symmetric saddle. \section{General solution: examples} \subsection{1RSB complexity} It is known that by choosing a covariance $f$ as the sum of polynomials with well-separated powers, one develops 2RSB in equilibrium. This should correspond to 1RSB in Kac--Rice. For this example, we take \begin{equation} f(q)=\frac12\left(q^3+\frac1{16}q^{16}\right) \end{equation} With this covariance, the model sees a RS to 1RSB transition at $\beta_1=1.70615\ldots$ and a 1RSB to 2RSB transition at $\beta_2=6.02198\ldots$. At these points, the average energies are $\langle E\rangle_1=-0.906391\ldots$ and $\langle E\rangle_2=-1.19553\ldots$, and the ground state energy is $E_0=-1.2876055305\ldots$. In this model, the RS complexity gives an inconsistent answer for the complexity of the ground state, predicting that the complexity of minima vanishes at a higher energy than the complexity of saddles, with both at a lower energy than the equilibrium ground state. The 1RSB complexity resolves these problems, predicting the same ground state as equilibrium and that the complexity of marginal minima (and therefore all saddles) vanishes at $E_m=-1.2876055265\ldots$, which is very slightly greater than $E_0$. Saddles become dominant over minima at a higher energy $E_s=-1.287605716\ldots$. Finally, the 1RSB complexity transitions to a RS description at an energy $E_1=-1.27135996\ldots$. All these complexities can be seen plotted in Fig.~\ref{fig:2rsb.complexity}. All of the landmark energies associated with the complexity are a great deal smaller than their equilibrium counterparts, e.g., comparing $E_1$ and $\langle E\rangle_2$. \begin{figure} \includegraphics{figs/316_complexity.pdf} \caption{ Complexity of dominant saddles (blue), marginal minima (yellow), and dominant minima (green) of the $3+16$ model. Solid lines show the result of the 1RSB ansatz, while the dashed lines show that of a RS ansatz. The complexity of marginal minima is always below that of dominant critical points except at the red dot, where they are dominant. The inset shows a region around the ground state and the fate of the RS solution. } \label{fig:2rsb.complexity} \end{figure} \begin{figure} \centering \includegraphics{figs/316_complexity_contour_1.pdf} \hfill \includegraphics{figs/316_complexity_contour_2.pdf} \caption{ Complexity of the $3+16$ model in the energy $E$ and radial reaction $\mu$ plane. The right shows a detail of the left. The black line shows $\mu_m$, which separates minima above from saddles below. The white lines show the dominant stationary points at each energy, dashed when they are described by a RS solution and solid for 1RSB. The red line shows the transition between RS and 1RSB descriptions. The gray line shows where the RS description predicts zero complexity. } \end{figure} \begin{figure} \centering \begin{minipage}{0.7\textwidth} \includegraphics{figs/316_comparison_q.pdf} \hspace{1em} \includegraphics{figs/316_comparison_x.pdf} \vspace{1em}\\ \includegraphics{figs/316_comparison_b.pdf} \hspace{1em} \includegraphics{figs/316_comparison_R.pdf} \end{minipage} \includegraphics{figs/316_comparison_legend.pdf} \caption{ Comparisons between the saddle parameters of the equilibrium solution to the $3+16$ model (blue) and those of the complexity (yellow). Equilibrium parameters are plotted as functions of the average energy $\langle E\rangle=\partial_\beta(\beta F)$ and complexity parameters as functions of fixed energy $E$. Solid lines show the result of a 2RSB ansatz, dashed lines that of a 1RSB ansatz, and dotted lines that of a RS ansatz. All paired parameters coincide at the ground state energy, as expected. } \label{fig:2rsb.comparison} \end{figure} \subsection{Full RSB complexity} \begin{figure} \centering \includegraphics{figs/24_complexity.pdf} \caption{ The complexity $\Sigma$ of the mixed $2+4$ spin model as a function of distance $\Delta E=E-E_0$ of the ground state. The solid blue line shows the complexity of dominant saddles given by the FRSB ansatz, and the solid yellow line shows the complexity of marginal minima. The dashed lines show the same for the annealed complexity. The inset shows more detail around the ground state. } \label{fig:frsb.complexity} \end{figure} \begin{figure} \centering \includegraphics{figs/24_func.pdf} \hspace{1em} \includegraphics{figs/24_qmax.pdf} \caption{ \textbf{Left:} The spectrum $\chi$ of the replica matrix in the complexity of dominant saddles for the $2+4$ model at several energies. \textbf{Right:} The cutoff $q_{\mathrm{max}}$ for the nonlinear part of the spectrum as a function of energy $E$ for both dominant saddles and marginal minima. The colored vertical lines show the energies that correspond to the curves on the left. } \label{fig:24.func} \end{figure} \begin{figure} \centering \includegraphics{figs/24_comparison_b.pdf} \hspace{1em} \includegraphics{figs/24_comparison_Rd.pdf} \raisebox{3em}{\includegraphics{figs/24_comparison_legend.pdf}} \caption{ Comparisons between the saddle parameters of the equilibrium solution to the $3+4$ model (black) and those of the complexity (blue and yellow). Equilibrium parameters are plotted as functions of the average energy $\langle E\rangle=\partial_\beta(\beta F)$ and complexity parameters as functions of fixed energy $E$. Solid lines show the result of a FRSB ansatz and dashed lines that of a RS ansatz. All paired parameters coincide at the ground state energy, as expected. } \label{fig:2rsb.comparison} \end{figure} Here a picture of $\chi$ vs $C$ or $X$ vs $C$ showing limits $q_{max}$, $x_{max}$ for different energies and typical vs minima. \section{Interpretation} \label{sec:interpretation} Let $\langle A\rangle$ be the average of $A$ over stationary points with given $E$ and $\mu^*$, i.e., \begin{equation} \langle A\rangle =\frac1{\mathcal N}\sum_{\sigma}A(s_\sigma) =\frac1{\mathcal N} \int d\nu(s)\,A(s) \end{equation} with \begin{equation} d\nu(s)=ds\,d\mu\,\delta\big(\tfrac12(\|s\|^2-N)\big)\,\delta\big(\nabla H(s,\mu)\big)\,\big|\det\operatorname{Hess}H(s,\mu)\big| \delta\big(NE-H(s)\big)\delta\big(N\mu^*-\operatorname{Tr}\operatorname{Hess}H(s,\mu)\big) \end{equation} the Kac--Rice measure. The fields $C$, $R$, and $D$ defined in \eqref{eq:fields} can be related to certain averages of this type. \subsection{\textit{C}: distribution of overlaps} First, consider $C$, which has an interpretation nearly identical to that of Parisi's $Q$ matrix of overlaps. It can be shown that its off-diagonal corresponds to the probability distribution of the overlaps between stationary points $P(q)$. First, define this distribution as \begin{equation} P(q)=\frac1{\mathcal N^2}\sum_{\sigma,\sigma'}\delta\left(\frac{s_\sigma\cdot s_{\sigma'}}N-q\right) \end{equation} where the sum is twice over stationary points $\sigma$ and $\sigma'$. It is straightforward to show that moments of this distribution are related to certain averages of the form. These are evaluated for a given energy, index, etc, but we shall omit these subindices for simplicity. \begin{eqnarray} q^{(p)} &\equiv& \frac1{N^p}\sum_{i_1\cdots i_p}\langle s^1_{i_1}\cdots s^1_{i_p}\rangle\langle s^2_{i_1}\cdots s^2_{i_p}\rangle \nonumber \\ &=&\frac1{N^p} \; \frac{1}{{\cal{N}}^2} \left\{ \sum_{{\mathbf s}^1,{\mathbf s}^2}\; \sum_{i_1\cdots i_p} s^1_{i_1}\cdots s^1_{i_p} s^2_{i_1}\cdots s^2_{i_p}\right\} =\frac1{N^p} \lim_{n\to0} \left\{ \sum_{{\mathbf s}^1,{\mathbf s}^2...{\mathbf s}^n}\; \sum_{i_1\cdots i_p} s^1_{i_1}\cdots s^1_{i_p} s^2_{i_1}\cdots s^2_{i_p}\right\} \end{eqnarray} The $(n-2)$ extra replicas providing the normalization. The average over disorder, and again, for given enegrgy, index,etc reads: \begin{eqnarray} \overline{q^{(p)}} &=&\overline{\frac1{N^p}\sum_{i_1\cdots i_p}\langle s_{i_1}\cdots s_{i_p}\rangle\langle s_{i_1}\cdots s_{i_p}\rangle} =\lim_{n\to0}{\int\prod_\alpha^n d\nu(s_\alpha)\,\left(\frac{s_1 s_2}N\right)^p} \times \Big[ {\mbox {REPLICATED AVERAGED MEASURE}} \Big]\nonumber \\ &=&\lim_{n\to0}{\int D[C,R,D] \, \left(C_{12}\right)^p\; } e^{N\Sigma[C,R,D]} =\frac1{n(n-1)}\lim_{n\to0}{\int D[C,R,D] \,\sum_{a\neq b}\left(C_{ab}\right)^p\; } e^{N\Sigma[C,R,D]} \end{eqnarray} In the last line, we have used that there is nothing special about replicas one and two. Using the Parisi ansatz, evaluating by saddle point {\em summing over all the $n(n-1)$ saddles related by permutation} we then have \begin{equation} \overline{q^{(p)}}=\int_0^1 dx\,c^p(x) = \int dq\,q^p P(q) \qquad ; \qquad {\mbox{with}} \qquad P(q)=\frac{dx}{dq} \end{equation} The appeal of Parisi to properties of pure states is unnecessary here, since the stationary points are points. With this established, we now address what it means for $C$ to have a nontrivial replica-symmetry broken structure. When $C$ is replica symmetric, drawing two stationary points at random will always lead to the same overlap. In the case when there is no linear field, they will always have overlap zero, because the second point will almost certainly lie on the equator of the sphere with respect to the first. Though other stationary points exist nearby the first one, they are exponentially fewer and so will be picked with vanishing probability in the thermodynamic limit. When $C$ is replica-symmetry broken, there is a nonzero probability of picking a second stationary point at some other overlap. This can be interpreted by imagining the level sets of the Hamiltonian in this scenario. If the level sets are disconnected but there are exponentially many of them distributed on the sphere, one will still find zero average overlap. However, if the disconnected level sets are \emph{few}, i.e., less than order $N$, then it is possible to draw two stationary points from the same set. Therefore, the picture in this case is of few, large basins each containing exponentially many stationary points. \subsection{A concrete example} Consider two independent pure $p$ spin models $H_{p_1}({\mathbf s})$ and $H_{p_2}({\mathbf \sigma})$ of sizes $N$, and couple them weakly with $\varepsilon \; {\mathbf \sigma} \cdot {\mathbf s}$. The complexities are \begin{eqnarray} e^{N\Sigma(e)}&=&\int de_1 de_2 \; e^{N[ \Sigma_1(e_1) + \Sigma_2(e_2) + O(\varepsilon) -\lambda N [(e_1+e_2)-e]}\nonumber \\ e^{-G(\hat \beta)}&=&\int de de_1 de_2 \; e^{N[-\hat \beta e+ \Sigma_1(e_1) + \Sigma_2(e_2) + O(\varepsilon) -\lambda N [(e_1+e_2)-e]} \end{eqnarray} The maximum is given by $\Sigma_1'=\Sigma_2'=\hat \beta$, provided it occurs in the phase in which both $\Sigma_1$ and $\Sigma_2$ are non-zero. The two systems are `thermalized', and it is easy to see that, because many points contribute, the overlap between two global configurations $$\frac1 {2N}({\mathbf s^1},{\mathbf \sigma^1})\cdot ({\mathbf s^2},{\mathbf \sigma^2})=\frac1 {2N}[ {\mathbf s^1}\cdot {\mathbf s^2}+ {\mathbf \sigma^1}\cdot {\mathbf \sigma^2}] =0$$ This is the `annealed' phase of a Kac-Rice calculation. Now start going down in energy, or up in $\hat \beta$: there will be a point $e_c$, $\hat \beta_c$ at which one of the subsystems freezes at its lower energy density, say it is system one, while system two is not yet frozen. At an even higher $\hat \beta=\hat \beta_f$, both systems are frozen. For $\hat \beta_f > \hat \beta> \hat \beta_c$ one system is unfrozen, while the other is, because of coupling, frozen at inverse temperature $\hat \beta_c$. The overlap between two solutions is: $$\frac1 {2N}({\mathbf s^1},{\mathbf \sigma^1})\cdot ({\mathbf s^2},{\mathbf \sigma^2}) =\frac1 {2N}[ {\mathbf s^1}\cdot {\mathbf s^2}+ {\mathbf \sigma^1}\cdot {\mathbf \sigma^2}] = \frac1 {2N} {\mathbf s^1}\cdot {\mathbf s^2} $$ The distribution of this overlap is the one of a frozen spin-glass at temperature $\hat \beta$, a $1RSB$ system like the Random Energy Model. The value of $x$ corresponding to it depends on $\hat \beta$, starting at $x=1$ at $\hat \beta_c$ and decreasing wuth increasing $\hat \beta$. Globally, the joint Kac-Rice system is $1RSB$, but note that the global overlap between different states is at most $1/2$. At $\hat \beta>\hat \beta_f$ there is a further transition. \subsection{\textit{R} and \textit{D}: response functions} The matrix field $R$ is related to responses of the stationary points to perturbations of the tensors $J$: \begin{equation} \begin{aligned} \overline{\frac1{N^p}\sum_{i_1\cdots i_p}\frac{\partial\langle s_{i_1}\cdots s_{i_p}\rangle}{\partial J^{(p)}_{i_1\cdots i_p}}} =\lim_{n\to0}\overline{\int\prod_\alpha^nd\nu(s_\alpha)\,\frac1n\sum_{ab}\left[ \hat\beta\left(\frac{s_a\cdot s_b}N\right)^p+ p\left(-i\frac{\hat s_a\cdot s_b}N\right)\left(\frac{s_a\cdot s_b}N\right)^{p-1} \right]} \\ =\lim_{n\to0}\frac1n\sum_{ab}(\hat\beta C_{ab}^p+pR_{ab}C_{ab}^{p-1}) =\hat\beta+pr_d-\int_0^1dx\,c^{p-1}(x)(\hat\beta c(x)+pr(x)) \end{aligned} \end{equation} In particular, when the energy is unconstrained ($\hat\beta=0$) and there is replica symmetry, \begin{equation} \frac1N\sum_i\frac{\partial\langle s_i\rangle}{\partial J_i^{(1)}} =r_d \end{equation} i.e., adding a linear field causes a response in the average stationary point location proportional to $r_d$. If positive, for instance, stationary points tend to align with a field. The energy constraint has a significant contribution due to the perturbation causing stationary points to move up or down in energy. The matrix field $D$ is related to the response of the complexity to such perturbations: \begin{equation} \begin{aligned} \frac{\partial\Sigma}{\partial a_p} =\frac14\lim_{n\to0}\frac1n\sum_{ab}^n\left[ \hat\beta^2C_{ab}^p+p(2\hat\beta R_{ab}-D_{ab})C_{ab}^{p-1}+p(p-1)R_{ab}^2C_{ab}^{p-2} \right] \end{aligned} \end{equation} In particular, when the energy is unconstrained ($\hat\beta=0$) and there is no replica symmetry breaking, \begin{equation} \frac{\partial\Sigma}{\partial a_1}=-\frac14\lim_{n\to0}\frac1n\sum_{ab}D_{ab}=-\frac14d_d \end{equation} i.e., adding a random linear field decreases the complexity of solutions by an amount proportional to $d_d$ in the variance of the field. When the saddle point of the Kac--Rice problem is supersymmetric, \begin{equation} \frac{\partial\Sigma}{\partial a_p} =\frac{\hat\beta}4\overline{\frac1{N^p}\sum_{i_1\cdots i_p}\frac{\partial\langle s_{i_1}\cdots s_{i_p}\rangle}{\partial J^{(p)}_{i_1\cdots i_p}}}+\lim_{n\to0}\frac1n\sum_{ab}^np(p-1)R_{ab}^2C_{ab}^{p-2} \end{equation} and in particular for $p=1$ \begin{equation} \frac{\partial\Sigma}{a_1} =\frac{\hat\beta}4\overline{\frac1N\sum_i\frac{\partial\langle s_i\rangle}{\partial J_i^{(1)}}} \end{equation} i.e., the change in complexity due to a linear field is directly related to the resulting magnetization of the stationary points. \section{Conclusion} We have constructed a replica solution for the general problem of finding saddles of random mean-field landscapes, including systems with many steps of RSB. For systems with full RSB, we find that minima are, at all energy densities above the ground state, exponentially subdominant with respect to saddles. The solution contains valuable geometric information that has yet to be extracted in all detail. A first and very important application of the method here is to perform the calculation for high dimensional spheres, where it would give us a clear understanding of what happens in a low-temperature realistic jamming dynamics \cite{maimbourg2016solution}. \paragraph{Funding information} J K-D and J K are supported by the Simons Foundation Grant No. 454943. \begin{appendix} \section{RSB for the Gibbs-Boltzmann measure} \begin{equation} \beta F=-\frac12\lim_{n\to0}\frac1n\left(\beta^2\sum_{ab}f(Q_{ab})+\log\det Q\right)-\frac12\log S_\infty \end{equation} $\log S_\infty=1+\log2\pi$. \begin{align*} \beta F= -\frac12\log S_\infty -\frac12\lim_{n\to0}\frac1n\left(\beta^2nf(1)+\beta^2\sum_{i=0}^kn(x_i-x_{i+1})f(q_i) +\log\left[ \frac{ 1+\sum_{i=0}^k(x_i-x_{i+1})q_i }{ 1+\sum_{i=1}^k(x_i-x_{i+1})q_i-x_1q_0 } \right]\right.\\ +\frac n{x_1}\log\left[ 1+\sum_{i=1}^k(x_i-x_{i+1})q_i-x_1q_0 \right]\\ \left.+\sum_{j=1}^kn(x_{j+1}^{-1}-x_j^{-1})\log\left[ 1+\sum_{i=j+1}^k(x_i-x_{i+1})q_i-x_{j+1}q_j \right] \right) \end{align*} \begin{align*} \lim_{n\to0}\frac1n \log\left[ \frac{ 1+\sum_{i=0}^k(x_i-x_{i+1})q_i }{ 1+\sum_{i=1}^k(x_i-x_{i+1})q_i-x_1q_0 } \right] &= \lim_{n\to0}\frac1n \log\left[ \frac{ 1+\sum_{i=0}^k(x_i-x_{i+1})q_i }{ 1+\sum_{i=0}^k(x_i-x_{i+1})q_i-nq_0 } \right] \\ &=q_0\left(1+\sum_{i=0}^k(x_i-x_{i+1})q_i\right)^{-1} \end{align*} \begin{align*} \beta F= -\frac12\log S_\infty -\frac12\left(\beta^2f(1)+\beta^2\sum_{i=0}^k(x_i-x_{i+1})f(q_i) +q_0\left(1+\sum_{i=0}^k(x_i-x_{i+1})q_i\right)^{-1}\right. \\ +\frac1{x_1}\log\left[ 1+\sum_{i=1}^{k}(x_i-x_{i+1})q_i-x_1q_0 \right]\\ \left.+\sum_{j=1}^k(x_{j+1}^{-1}-x_j^{-1})\log\left[ 1+\sum_{i=j+1}^{k}(x_i-x_{i+1})q_i-x_{j+1}q_j \right] \right) \end{align*} $q_0=0$ \begin{align*} \beta F= -\frac12\log S_\infty -\frac12\left(\beta^2f(1)+\beta^2\sum_{i=0}^k(x_i-x_{i+1})f(q_i) +\frac1{x_1}\log\left[ 1+\sum_{i=1}^{k}(x_i-x_{i+1})q_i \right]\right.\\ \left.+\sum_{j=1}^k(x_{j+1}^{-1}-x_j^{-1})\log\left[ 1+\sum_{i=j+1}^{k}(x_i-x_{i+1})q_i-x_{j+1}q_j \right] \right) \end{align*} $x_i=\tilde x_ix_k$, $x_k=y/\beta$, $q_k=1-z/\beta$ \begin{align*} \beta F= -\frac12\log S_\infty- \frac12\left(\beta^2f(1)+\beta^2(y\beta^{-1}-1)f(1-z\beta^{-1})+y\beta\sum_{i=0}^{k-1}(\tilde x_i-\tilde x_{i+1})f(q_i)\right. \\ +\frac\beta{\tilde x_1 y}\log\left[ y\sum_{i=1}^{k-1}(\tilde x_i-\tilde x_{i+1})q_i+y+z-yz/\beta \right]\\ +\sum_{j=1}^{k-1}\frac\beta y(\tilde x_{j+1}^{-1}-\tilde x_j^{-1})\log\left[ y\sum_{i=j+1}^{k-1}(\tilde x_i-\tilde x_{i+1})q_i+y+z-yz/\beta-y\tilde x_{j+1}q_j \right]\\ \left.-\frac\beta{\tilde x_1 y}\log\beta-\sum_{j=1}^{k-1}\frac\beta y(\tilde x_{j+1}^{-1}-\tilde x_j^{-1})\log\beta+(1-\beta y^{-1})\log\left[ z/\beta \right] \right) \end{align*} \begin{align*} \lim_{\beta\to\infty}F= -\frac12\left(yf(1)+zf'(1)+y\sum_{i=0}^{k-1}(\tilde x_i-\tilde x_{i+1})f(q_i) +\frac1{\tilde x_1 y}\log\left[ y\sum_{i=1}^{k-1}(\tilde x_i-\tilde x_{i+1})q_i+y+z \right]\right.\\ \left.+\sum_{j=1}^{k-1}\frac1 y(\tilde x_{j+1}^{-1}-\tilde x_j^{-1})\log\left[ y\sum_{i=j+1}^{k-1}(\tilde x_i-\tilde x_{i+1})q_i+y+z-y\tilde x_{j+1}q_j \right] -\frac1y\log z \right) \end{align*} $F$ is a $k-1$ RSB ansatz with all eigenvalues scaled by $y$ and shifted by $z$. $\tilde x_0=0$ and $\tilde x_k=1$. \section{ RSB for the Kac--Rice integral} \subsection{Diagonal} \begin{align*} \lim_{n\to0}\frac1n\log\det(\hat\beta R_d^{-1} Q+I) =x_1^{-1}\log\left(\hat\beta R_d^{-1}(1-\bar q_k)+1\right)+\int_{q_0^+}^{q_{k-1}}dq\,\mu(q)\log\left[\hat\beta R_d^{-1}\lambda(q)+1\right] \end{align*} where \[ \mu(q)=\frac{\partial x^{-1}(q)}{\partial q} \] Integrating by parts, \begin{align*} \lim_{n\to0}\frac1n\log\det(\hat\beta R_d^{-1} Q+I) &=x_1^{-1}\log\left(\hat\beta R_d^{-1}(1-\bar q_k)+1\right)+\left[x^{-1}(q)\log[\hat\beta R_d^{-1}\lambda(q)+1]\right]_{q=q_0^+}^{q=q_{k-1}}-\frac{\hat\beta}{R_d}\int_{q_0^+}^{q_{k-1}}dq\,\frac{\lambda'(q)}{x(q)}\frac1{\hat\beta R_d^{-1}\lambda(q)+1}\\ &=\log[\hat\beta R_d^{-1}\lambda(q_{k-1})+1]+\frac{\hat\beta}{R_d}\int_{q_0^+}^{q_{k-1}}dq\,\frac1{\hat\beta R_d^{-1}\lambda(q)+1} \end{align*} \subsection{Non-diagonal} \end{appendix} \bibliographystyle{plain} \bibliography{frsb_kac-rice} \end{document}