\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.
These include most notably the (pure) $p$-spin model ($p>2$) \cite{Rieger_1992_The, Crisanti_1995_Thouless-Anderson-Palmer}. 
The problem of studying the critical points of these landscapes
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,  including the Sherrington--Kirkpatrick model. It reproduces  the Parisi result in the limit
of small temperature for the lowest states, as it should.


\section{The model}

Here we consider, for definiteness, the  mixed $p$-spin model, itself a particular case 
of the `Toy Model' of M\'ezard and Parisi \cite{Mezard_1992_Manifolds}
\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}
Can 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 hessian is
\begin{equation}
  \operatorname{Hess}H=\partial\partial H+\mu I
\end{equation}
$\partial\partial H$ is 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$ fixes the spectrum of the hessian. By manipulating it, one
can decide to find the complexity of saddles of a certain macroscopic index, or
of minima with a certain harmonic stiffness. When $\mu$ is taken to be within
the range $\pm\sqrt{4f''(1)}=\pm\mu_m$, the critical points are constrained to have
index
\begin{equation}
  \mathcal I(\mu)=N\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$. Finally,
when $\mu=\mu_m$, the critical points are marginal minima.


The parameter $\mu$ fixes the spectrum of the hessian. When it is an integration variable,
and one restricts the domain of all integrations to compute saddles of a certain macroscopic index, or
of minima with a certain harmonic stiffness, its value is the `softest' mode that adapts to change the Hessian \cite{Fyodorov_2007_Replica}. When it is fixed, then the restriction of the index of saddles is `payed' by the realization of the eigenvalues of the Hessian, usually a
`harder' mode.

\section{Equilibrium}

Here we review the equilibrium solution. \cite{Crisanti_1992_The, Crisanti_1993_The, Crisanti_2004_Spherical, Crisanti_2006_Spherical}
The free energy is well known to 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 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}

\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$, giving the number of stationary points at
energy $E$ and radial reaction $\mu$ as
\begin{equation}
  \mathcal N(E, \mu)
  =\int ds\,\delta(NE-H(s))\delta(\partial H(s)+\mu s)|\det(\partial\partial H(s)+\mu I)|
\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 complexity:
\begin{equation}
  \Sigma(E,\mu)=\lim_{N\to\infty}\frac1N\overline{\log\mathcal N(E, \mu)}
\end{equation}
The radial reaction $\mu$, which acts like a kind of `mass' term, takes a
fixed value here, which means that the complexity is for a given energy
density and hessian spectrum. This will turn out to be important when we
discriminate between counting all solutions, or selecting those of a given
index, for example minima. The complexity of solutions without fixing $\mu$ is
given by maximizing the complexity as a function of $\mu$.

If one averages over $\mathcal N$ and afterward takes its logarithm, one arrives at the 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}.

\subsection{The replicated problem}

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 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\,\delta(NE-H(s_a))\delta(\partial H(s_a)+\mu s_a)|\det(\partial\partial H(s_a)+\mu I)|
  \end{aligned}
\end{equation}

As noted by Bray and Dean  \cite{Bray_2007_Statistics}, gradient and Hessian
are independent for a random Gaussian function, and the average over disorder
breaks into a product of two independent averages, one for the gradient factor
and one for the determinant. The integration of all variables, including the
disorder in the last factor, may be restricted to the domain such that the
matrix $\partial\partial H(s_a)-\mu I$ has  a specified number of negative
eigenvalues (the index $\mathcal I$ of the saddle), (see Fyodorov
\cite{Fyodorov_2007_Replica} for a detailed discussion). 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\right)\,\overline{\prod_a^n \delta(NE-H(s_a))\delta(\partial H(s_a)+\mu s_a)}
    \times
    \overline{\prod_a^n |\det(\partial\partial H(s_a)+\mu I)|}
  \end{aligned}
\end{equation}
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}
    \mathcal D(\mu)
    &=\frac1N\overline{\log|\det(\partial\partial H(s_a)+\mu I)|}
    =\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}
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(NE-H(s_a))\delta(\partial H(s_a)+\mu s_a)
  =\int \frac{d\hat\beta}{2\pi}\prod_a^n\frac{d\hat s_a}{2\pi}
    e^{\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+\\
    &\lim_{n\to0}\frac1n\left(
      -\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\beta$, $C$, $R$ and $D$ must be evaluated at extrema of this
expression.

\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. 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\beta}
  =E+\sum_{ab}\left[\hat\beta f(C_{ab})+R_{ab}f'(C_{ab})\right] \label{eq:cond.b} \\
  \frac{\tilde c}2I&=\frac{\partial\Sigma}{\partial C}
  =\frac12\left[
    \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. The
equation for \eqref{eq:cond.q} would not be true on the diagonal save for the
arbitrary factor $\tilde c$.  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}

\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?

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}
  \tilde cI=\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}
  (\tilde c-\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{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}
\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}
Using standard manipulations (Appendix B), one finds also a continuous version
\begin{equation} \label{eq:functional.action}
  \begin{aligned}
    \Sigma(E,\mu)
    =\mathcal D(\mu)
    +
      \hat\beta E-\mu R_d
      +\frac12(2\hat\beta R_d-D_d)f'(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^2/D_d}
      \right)
  \end{aligned}
\end{equation}
where $\chi(q)$ is defined with respect to $Q$ exactly as in the equilibrium case.

\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}

In the case where any FRSB is present, one must work with the functional form
of the complexity \eqref{eq:functional.action}, which must be extremized with
respect to $\chi$ under the conditions that $\chi$ is concave, monotonically
decreasing, and $\chi(1)=0$, $\chi'(1)=-1$. The annealed case is found by
taking $\chi(q)=1-q$, which satisfies all of these conditions. $k$-RSB is
produced by breaking $\chi$ into $k+1$ piecewise linear segments.

Forget for the moment these tricky requirements. The function would then be
extremized by satisfying
\begin{equation}
  0=\frac{\delta\Sigma}{\delta\lambda(q)}=\frac12\hat\beta^2f''(q)-\frac12\frac1{(\lambda(q)+R_d^2/D_d)^2}
\end{equation}
which implies the solution
\begin{equation}
  \lambda^*(q)=\frac1{\hat\beta}f''(q)^{-1/2}-\frac{R_d^2}{D_d}
\end{equation}
If $f''(q)^{-1/2}$ is not concave anywhere, there is little use of this
solution. However, if it is concave everywhere it may constitute a portion of
the full solution.

We suppose that solutions are given by
\begin{equation}
  \lambda(q)=\begin{cases}
    \lambda^*(q) & q<q_\textrm{max} \\
    1-q          & q\geq q_\textrm{max}
  \end{cases}
\end{equation}
Continuity requires that $1-q_\textrm{max}=\lambda^*(q_\textrm{max})$.

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\,\delta(NE-H(s))\delta(\partial H(s)+\mu s)|\det(\partial\partial H(s)+\mu I)|
\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
\begin{equation}
  \int dq\,q^p P(q)
  =q^{(p)}
  \equiv\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
\end{equation}
The appeal of Parisi to properties of pure states is unnecessary here, since
the stationary points are points. These moments are related to our $C$ by
computing their average over disorder:
\begin{equation}
  \begin{aligned}
    \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}\overline{\int\prod_\alpha^nd\nu(s_\alpha)\,\frac1{n(n-1)}\sum_{a\neq b}\left(\frac{s_a\cdot s_b}N\right)^p} \\
    =\lim_{n\to0}\frac1{n(n-1)}\sum_{a\neq b}^nC^p_{ab}
    =\int_0^1 dx\,c^p(x)
  \end{aligned}
\end{equation}
where $C$ is assumed to take its saddle point value. If we now change variables in the integral from $x$ to $c$, we find
\begin{equation}
  \overline{q^{(p)}}=\int dc\,c^p\frac{dx}{dc}
\end{equation}
from which we conclude $\overline{P(q)}=\frac{dx}{dc}\big|_{c=q}$.

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{\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{Ultrametricity rediscovered}

TENTATIVE BUT INTERESTING

The frozen phase for a given index ${\cal{I}}$ is the one for values of $\hat \beta> \hat \beta_{freeze}^{\cal{I}}$. 

[Jaron: does $\hat \beta^I_{freeze}$ have a relation to the largest $x$  of the ansatz? If so, it would give an interesting interpretation for everything]
 
 The complexity of that index is zero, and we are looking at the lowest saddles 
in the problem, a question that to the best of our knowledge has not been discussed
in the Kac--Rice context -- for good reason, since the complexity - the original motivation - is zero.
However, our ansatz tells us something of the actual organization of  the lowest saddles of each index in phase space. 



\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. 
The main results of this paper are the ansatz (\ref{ansatz}) and the check that the lowest energy is the correct one obtained with the usual Parisi ansatz. 
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.

\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{Solution}



\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*}

\end{appendix}


\bibliographystyle{plain}
\bibliography{frsb_kac-rice}




\end{document}