path: root/frsb_kac-rice.tex

\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
\usepackage[style=phys,eprint=true]{biblatex}

\addbibresource{frsb_kac-rice.bib}

\begin{document}
\title{
  How to count in hierarchical landscapes: a `full' solution to mean-field complexity
}
\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 symmetry
breaking (1RSB) phase at a Kauzmann temperature, and, at a lower temperature,
another  transition to  a full RSB phase (see \cite{Gross_1985_Mean-field,
Gardner_1985_Spin}, the so-called `Gardner' phase
\cite{Charbonneau_2014_Fractal}).  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? This is a
question whose answers are significant to interpreting the results of myriad
experiments and simulations \cite{Xiao_2022_Probing, Hicks_2018_Gardner,
Liao_2019_Hierarchical, Dennis_2020_Jamming, Charbonneau_2015_Numerical,
Li_2021_Determining, Seguin_2016_Experimental, Geirhos_2018_Johari-Goldstein,
Hammond_2020_Experimental, Albert_2021_Searching} (see, for a review
\cite{Berthier_2019_Gardner}).  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'.
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.


Perhaps the most interesting application of this computation is in the context  of
optimization problems, see for example \cite{Gamarnik_2021_The, ElAlaoui_2022_Sampling, Huang_2021_Tight}. A question
that appears there is how to define a `threshold level'. This notion was introduced  \cite{Cugliandolo_1993_Analytical} in the context of the $p$-spin model, as the energy at which the patches of the same  energy in phase-space percolate  - hence
explaining why dynamics never go below that level.
The notion of a `threshold' for more complex landscapes has later been 
invoked several times, never to our knowledge in a clear and unambiguous 
way. One of the purposes of this paper is to give a sufficiently detailed
characterization of a general landscape so that a meaningful general notion
of threshold may  be introduced - if this is at all possible.


\section{The model}
\label{sec:model}

For definiteness, we consider the mixed $p$-spin spherical model, whose Hamiltonian
\begin{equation} \label{eq:hamiltonian}
  H(\mathbf s)=-\sum_p\frac1{p!}\sum_{i_1\cdots i_p}^NJ^{(p)}_{i_1\cdots i_p}s_{i_1}\cdots s_{i_p}
\end{equation}
is defined for vectors $\mathbf s\in\mathbb R^N$ confined to the sphere
$\|\mathbf s\|^2=N$.  The coupling coefficients $J$ are taken at random, with
zero mean and variance $\overline{(J^{(p)})^2}=a_pp!/2N^{p-1}$ chosen so that
the energy is typically extensive. The factors $a_p$ in the variances are
freely chosen constants that define the particular model. For instance, the
so-called `pure' models have $a_p=1$ for some $p$ and all others zero.

The variance of the couplings implies that the covariance of the energy with
itself depends only on the dot product (or overlap) between two configurations.
In particular, one has
\begin{equation}
  \overline{H(\mathbf s_1)H(\mathbf s_2)}=Nf\left(\frac{\mathbf s_1\cdot\mathbf s_2}N\right)
\end{equation}
where $f$ is defined by the series
\begin{equation}
  f(q)=\frac12\sum_pa_pq^p
\end{equation}
More generally, one need not start with a Hamiltonian like
\eqref{eq:hamiltonian} and instead invoke the covariance rule for arbitrary,
non-polynomial $f$, as in the `toy model' of M\'ezard and Parisi
\cite{Mezard_1992_Manifolds}.

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(\mathbf s)+\frac\mu2(\|\mathbf s\|^2-N)
\end{equation}
At any critical point, the gradient and Hessian are
\begin{align}
  \nabla H(\mathbf s,\mu)=\partial H(\mathbf s)+\mu\mathbf s &&
  \operatorname{Hess}H(\mathbf s,\mu)=\partial\partial H(\mathbf s)+\mu I
\end{align}
where $\partial=\frac\partial{\partial\mathbf s}$ always.  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 d\mathbf s \;\delta(\|\mathbf s\|^2-N)\, e^{-\beta H(\mathbf 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 $n\times n$ matrix $Q$, whose diagonal must
be one. When the solution is a Parisi matrix, the free energy 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 define $\chi(q)$ by
\begin{equation}
  \chi(q)=\int_q^1dq'\,\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$ and $\chi'(1)=-1$.
Using standar arguments, the free energy can be written as a functional over
$\chi$ as
\begin{equation}
  \beta F=-1-\log2\pi-\frac12\int_0^1dq\,\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{Landscape complexity}

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}. It gives the number of stationary points $\mathcal N$ as
\begin{equation}
  \mathcal N
    =\int d\mathbf s\, d\mu\,\delta\big(\tfrac12(\|\mathbf s\|^2-N)\big)\,\delta\big(\nabla H(\mathbf s,\mu)\big)\,\big|\det\operatorname{Hess}H(\mathbf s,\mu)\big|
\end{equation}
It is more interesting to count stationary points which share certain
properties, like energy density $E$ or index $\mathcal I$. These properties can
be fixed by inserting additional $\delta$-functions into the integral. Rather
than fix the index directly, we fix the trace of the hession, which we'll soon
show is equivalent to fixing the value $\mu$, and fixing $\mu$ fixes the index
to within order one. Inserting these $\delta$-functions, we arrive at
\begin{equation}
  \begin{aligned}
    \mathcal N(E, \mu^*)
    &=\int d\mathbf s\, d\mu\,\delta\big(\tfrac12(\|\mathbf s\|^2-N)\big)\,\delta\big(\nabla H(\mathbf s,\mu)\big)\,\big|\det\operatorname{Hess}H(\mathbf s,\mu)\big| \\
    &\hspace{10pc}\times\delta\big(NE-H(\mathbf s)\big)\delta\big(N\mu^*-\operatorname{Tr}\operatorname{Hess}H(\mathbf 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}.  The annealed complexity is known to equal the
actual (quenched) complexity in circumstances where there is at most one level
of replica symmetry breaking in the model's equilibrium. 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)$ (a
model we study in detail later), the annealed 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}

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, we must deal with the
logarithm. We use the standard replica trick to convert the logarithm into a
product, which gives
\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 d\mathbf s_a\,d\mu_a\,
    \delta\big(\tfrac12(\|\mathbf s_a\|^2-N)\big)\,\delta\big(\nabla H(\mathbf s_a,\mu_a)\big)\,\big|\det\operatorname{Hess}H(\mathbf s_a,\mu_a)\big| \\
    &\hspace{13pc}  \times\delta\big(NE-H(\mathbf s_a)\big)\delta\big(N\mu^*-\operatorname{Tr}\operatorname{Hess}H(\mathbf s_a,\mu_a)\big)
  \end{aligned}
\end{equation}
As discussed in \S\ref{sec:model}, it has been shown that to the largest order
in $N$, the Hessian of Gaussian random functions in independent from their
gradient, once both are conditioned on certain properties. Here, they are only
related by their shared value of $\mu$. Because of this statistical
independence, we may write
\begin{equation}
  \begin{aligned}
    \Sigma(E, \mu^*)
    &=\lim_{N\to\infty}\frac1N\lim_{n\to0}\frac\partial{\partial n}\int\left(\prod_a^nd\mathbf s_a\,d\mu_a\right)\,
    \overline{\prod_a^n \delta\big(\tfrac12(\|\mathbf s_a\|^2-N)\big)\,\delta\big(\nabla H(\mathbf s_a,\mu_a)\big)\delta(NE-H(\mathbf s_a))}\\
    &\hspace{10pc}
    \times
    \overline{\prod_a^n |\det\operatorname{Hess}(\mathbf s_a,\mu_a)|\,\delta\big(N\mu^*-\operatorname{Tr}\operatorname{Hess}H(\mathbf s_a,\mu_a)\big)}
  \end{aligned}
\end{equation}
which simplifies matters. The average of the two factors may now be treated separately.

\subsubsection{The Hessian factors}

The spectrum of the matrix $\partial\partial H(\mathbf s)$ is uncorrelated from the
gradient. In the large-$N$ limit, for almost every point and realization of
disorder it is a GOE matrix with variance
\begin{equation}
  \overline{(\partial_i\partial_jH(\mathbf s))^2}=\frac1Nf''(1)\delta_{ij}
\end{equation}
Therefore in that limit its spectrum is given by the Wigner semicircle with radius $\sqrt{4f''(1)}$, or
\begin{equation}
  \rho(\lambda)=\begin{cases}
    \frac1{2\pi f''(1)}\sqrt{4f''(1)-\lambda^2} & \lambda^2\leq 4f''(1) \\
    0 & \text{otherwise}
  \end{cases}
\end{equation}
The spectrum of the Hessian $\operatorname{Hess}H(\mathbf s,\mu)$ is the same
semicircle 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\mu_m\equiv\pm\sqrt{4f''(1)}$, the critical points have index density
\begin{equation}
  \mathcal I(\mu)=\int_0^\infty d\lambda\,\rho(\lambda+\mu)
  =\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]
\end{equation}
When $\mu>\mu_m$, the critical points are minima whose sloppiest eigenvalue is
$\mu-\mu_m$, and when $\mu=\mu_m$, the critical points are marginal minima.


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$ \cite{Ros_2019_Complex}. We therefore find
\begin{equation}
  \overline{\prod_a^n |\det\operatorname{Hess}(\mathbf s_a,\mu_a)|\,\delta\big(N\mu^*-\operatorname{Tr}\operatorname{Hess}H(\mathbf s_a,\mu_a)\big)}
 \rightarrow \prod_a^ne^{N{\cal D}(\mu_a)}\delta(\mu_a-\mu^*)
\end{equation}
where the function $\mathcal D$ is defined by
\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}
It follows that by fixing the trace of the Hessian, we have effectively fixed
the value of $\mu$ in all replicas to $\mu^*$, and therefore the index of
saddles in all replicas as well.

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:
\begin{itemize}
  \item For $|\mu|>\mu_m$ there is the possibility of a finite number of eigenvalues of
    the second derivative matrix
  \item The second
\end{itemize}

\subsubsection{The gradient factors}

The $\delta$-functions are treated by writing them in the Fourier basis.
Introducing auxiliary fields $\hat{\mathbf s}_a$ and $\hat\beta$, for each
replica replica one writes
\begin{equation}
  \begin{aligned}
    &\delta\big(\tfrac12(\|\mathbf s_a\|^2-N)\big)\,\delta\big(\nabla H(\mathbf s_a,\mu^*)\big)\delta(NE-H(\mathbf s_a)) \\
    &\hspace{12pc}=\int\frac{d\hat\mu}{2\pi}\,\frac{d\hat\beta}{2\pi}\,\frac{d\hat{\mathbf s}_a}{(2\pi)^N}
      e^{\frac12\hat\mu(\|\mathbf s_a\|^2-N)+\hat\beta(NE-H(\mathbf s_a))+i\hat{\mathbf s}_a\cdot(\partial H(\mathbf s_a)+\mu^*\mathbf s_a)}
  \end{aligned}
\end{equation}
Anticipating the Parisi-style solution, we don't label $\hat\mu$ or $\hat\beta$
with replica indices, since replica vectors won't be broken in the scheme.  The
average over disorder can now be taken for the pieces which depend explicitly
on the Hamiltonian, and since everything is Gaussian this gives
\begin{equation}
  \begin{aligned}
    \overline{
      \exp\left[
        \sum_a^n(i\hat {\mathbf s}_a\cdot\partial_a-\hat\beta)H(s_a)
      \right]
    }
    &=\exp\left[
        \frac12\sum_{ab}^n
        (i\hat{\mathbf s}_a\cdot\partial_a-\hat\beta)
        (i\hat{\mathbf s}_b\cdot\partial_b-\hat\beta)
        \overline{H(\mathbf s_a)H(\mathbf s_b)}
      \right] \\
    &=\exp\left[
        \frac N2\sum_{ab}^n
        (i\hat{\mathbf s}_a\cdot\partial_a-\hat\beta)
        (i\hat{\mathbf s}_b\cdot\partial_b-\hat\beta)
        f\left(\frac{\mathbf s_a\cdot\mathbf s_b}N\right)
      \right] \\
    &\hspace{-14em}=\exp\left\{
        \frac N2\sum_{ab}^n
        \left[
          \hat\beta^2f\left(\frac{{\mathbf s}_a\cdot {\mathbf s}_b}N\right)
          -2i\hat\beta\frac{\hat {\mathbf s}_a\cdot {\mathbf s}_b}Nf'\left(\frac{{\mathbf s}_a\cdot {\mathbf s}_b}N\right)
          -\frac{\hat {\mathbf s}_a\cdot \hat {\mathbf s}_b}Nf'\left(\frac{{\mathbf s}_a\cdot {\mathbf s}_b}N\right)
          +\left(i\frac{\hat {\mathbf s}_a\cdot {\mathbf s}_b}N\right)^2f''\left(\frac{{\mathbf s}_a\cdot {\mathbf s}_b}N\right)
        \right]
      \right\}
  \end{aligned}
\end{equation}

We introduce new matrix fields
\begin{align} \label{eq:fields}
  C_{ab}=\frac1N\mathbf s_a\cdot\mathbf s_b &&
  R_{ab}=-i\frac1N\hat{\mathbf s}_a\cdot{\mathbf s}_b &&
  D_{ab}=\frac1N\hat{\mathbf s}_a\cdot\hat{\mathbf 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 $\mathbf s_a$ and $\hat{\mathbf
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\right.\\
    &\hspace{2em}\left.+\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 the extrema
of this expression which minimize the complexity. Extremizing with respect to
$\hat\mu$ is not difficult, and results in setting the diagonal of $C$ to one,
fixing the spherical constraint. Maintaining $\hat\mu$ in the complexity is
useful for writing down the extremal conditions, but when convenient we will
drop the dependence.

The same information is contained but better expressed in the Legendre
transform
\begin{equation}
  \begin{aligned}
    &G(\hat \beta,\mu^*)
    =\mathcal D(\mu^*)+\\
    &\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}
Denoting $r_d \equiv \frac 1 n {\mbox Tr} R$, we can write down the double Legendre transform $K(\hat \beta, r_d)$:
\begin{equation}
  e^{N K(\hat \beta, r_d)} =\int\,dE\,d\mu^* e^{N\left\{\Sigma(E,\mu^*) -\hat\beta E+r_d\mu^* -\mathcal D(\mu^*)\right\}}
\end{equation}
given by
\begin{equation}
  \begin{aligned}
    &K(\hat \beta,r_d)
   = \lim_{n\to0}\frac1n\left(
      \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 the diagonal of $C$ is fixed to one and the diagonal of $R$ is fixed to
$r_d$.  The variable $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 Sherrington--Kirkpatrick 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, we often want to maximize the complexity as a
function of $\mu^*$, to find the 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$ and the critical points are minima, \eqref{eq:cond.mu} implies
\begin{equation} \label{eq:mu.minima}
  \mu^*=\frac1{r_d}+r_df''(1)
\end{equation}
When $\mu^*<\mu_m$ and the critical points are saddles, it implies
\begin{equation} \label{eq:mu.saddles}
  \mu^*=2f''(1)r_d
\end{equation}

It is 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, which has been studied in great detail in the complexity of the Thouless--Anderson--Palmer free energy \cite{Annibale_2003_The, Annibale_2003_Supersymmetric, Annibale_2004_Coexistence}. 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 counts the most common minima}
\cite{Annibale_2004_Coexistence}. When minima are not the most common type of
stationary point, the supersymmetric solution correctly counts minima that
satisfy \eqref{eq:mu.minima}, but these do not have any special significance.

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

The supersymmetric solution produces the correct complexity for the ground
state and for a class of minima. Moreover, it produces the correct parameters
for the fields $C$, $R$, and $D$ at those points. This is an important foothold
in the problem of computing the general complexity. The full saddle point
equations at $k$-RSB are not very numerically stable, and a `good' saddle point
has a typically small radius of convergence under methods like Newton's
algorithm. With the supersymmetric solution in hand, it is possible to take
small steps in the parameter space to find non-supersymmetric numeric
solutions, each time ensuring the initial conditions for the solver are
sufficiently close to the correct answer. This is the strategy we use in
\S\ref{sec:examples}.

\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\left(\hat\beta r_df'(1)+r_d^2f''(1)+\log r_d^2\right)
      +\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.

Evaluated at $\mu^*_\mathrm{ss}=r_d^{-1}+f''(1)r_d$, the complexity further simplies to
\begin{equation} \label{eq:functional.action.ss}
  \Sigma(E,\mu^*_\mathrm{ss})
    =
    \hat\beta E+\frac12\left(
      \hat\beta f'(1)r_d-f''(1)r_d^2+\frac1{f''(1)r_d^2}
    \right)
    +\log(f''(1)r_d^2)
      +\frac12\int_0^1dq\,\left(
        \hat\beta^2f''(q)\chi(q)+\frac1{\chi(q)+r_d/\hat\beta}
      \right)
\end{equation}
At the ground state, the solution $\chi$ is smooth for all values of $q$. In order to satisfy the boundary conditions, we must have $r_d=f''(1)^{-1/2}$ and $\hat\beta=\frac12f'''(1)/f''(1)^{3/2}$

\subsection{Expansion near the transition}
\label{subsec:expansion}

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}
The equations cannot be used to find the value of $\bar c$ without going to
higher order in $x_\mathrm{max}$, but the transition line can be determined by
examining the stability of the replica symmetric complexity. First, we expand the full form for the complexity about small $x_\textrm{max}$ in the same way as we expand the extremal conditions. To quadratic order, this gives
\begin{equation}
  \begin{aligned}
  \Sigma(E,\mu^*)
  =\mathcal D(\mu^*)+\hat\beta E-\mu r_d+\frac12\left[\hat\beta^2f(1)+(2\hat\beta r_d-d_d)f'(1)+r_d^2f''(1)\right]+\frac12\log(d_d+r_d^2) \\
  -\frac12\left[
    \frac12\hat\beta^2\bar c^2f''(0)+(2\hat\beta\bar r-\bar d)\bar cf''(0)+\bar r^2f''(0)
    -\frac{\bar d^2-2d_d\bar r^2+d_d^2\bar c^2+4r_d\bar r(\bar d+d_d\bar c)-2r_d^2(\bar c\bar d+\bar r^2)}{2(d_d+r_d^2)^2}
  \right]x_\textrm{max}^2
  \end{aligned}
\end{equation}
The spectrum Hessian of this expression evaluated at $x_\text{max}=0$ gives the
stability of the replica symmetric solution with respect to perturbations of
the type described above. When the values of $\bar r$ and $\bar d$ above are substituted in and everything is evaluated at the replica symmetric solution, the  eigenvalue of interest takes the form
\begin{equation}
  \lambda
  =-\bar c^2\frac{(f'(1)-2f(1))^2(f'(1)-f''(0))f''(0)}{2(f'(1)+f''(0))(f'(1)^2-f(1)(f'(1)+f''(1)))^2}(\mu^*-\mu^*_+(E))(\mu^*-\mu^*_-(E))
\end{equation}
where
\begin{equation}
  \mu^*_\pm(E)
  =\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)}E
\end{equation}
This eigenvalues changes sign when $\mu^*$ crosses $\mu^*_\pm(E)$.  We expect
that this is the line of stability for the replica symmetric saddle.


\section{General solution: examples}
\label{sec: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}
established to have a 2RSB ground state \cite{Crisanti_2011_Statistical}.
With this covariance, the model sees a replica symmetric to 1RSB transition at
$\beta_1=1.70615\ldots$ and a 1RSB to 2RSB transition at
$\beta_2=6.02198\ldots$. At these transitions, the average energies are
$\langle E\rangle_1=-0.906391\ldots$ and $\langle E\rangle_2=-1.19553\ldots$,
respectively, and the ground state energy is $E_0=-1.2876055305\ldots$.
Besides these typical equilibrium energies, an energy of special interest for
looking at the landscape topology is the \emph{algorithmic threshold}
$E_\mathrm{alg}$, defined by the lowest energy reached by local algorithms like
approximate message passing \cite{ElAlaoui_2020_Algorithmic,
ElAlaoui_2021_Optimization}. In the spherical models, this has been proven to
be
\begin{equation}
  E_{\mathrm{alg}}=-\int_0^1dq\,\sqrt{f''(q)}
\end{equation}
For full RSB systems, $E_\mathrm{alg}=E_0$ and the algorithm can reach the
ground state energy. For the pure $p$-spin models,
$E_\mathrm{alg}=E_\mathrm{th}$, where $E_\mathrm{th}$ is the energy at which
marginal minima are the most common stationary points. Something about the
topology of the energy function is relevant to where this algorithmic threshold
lies. For the $3+16$ model at hand, $E_\mathrm{alg}=1.275140128\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(\mathbf s_\sigma)
  =\frac1{\mathcal N}
  \int d\nu(\mathbf s)\,A(\mathbf s)
\end{equation}
with
\begin{equation}
  d\nu(\mathbf s)=d\mathbf s\,d\mu\,\delta\big(\tfrac12(\|\mathbf s\|^2-N)\big)\,\delta\big(\nabla H(\mathbf s,\mu)\big)\,\big|\det\operatorname{Hess}H(\mathbf s,\mu)\big|
    \delta\big(NE-H(\mathbf s)\big)\delta\big(N\mu^*-\operatorname{Tr}\operatorname{Hess}H(\mathbf s,\mu)\big)
\end{equation}
the Kac--Rice measure. Note that this definition of the angle brackets is not
the same as that used in \S\ref{subsec:expansion} 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)$.
Let $\mathcal S$ be the set of all stationary points with given energy density
and index. Then
\begin{equation}
P(q)=\frac1{\mathcal N^2}\sum_{\mathbf s_1\in\mathcal S}\sum_{\mathbf s_2\in\mathcal S}\delta\left(\frac{\mathbf s_1\cdot\mathbf s_2}N-q\right)
\end{equation}
{\em This is the probability that two stationary points randomly drawn from the ensemble
of stationary points happen to be at overlap $q$}

%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{equation}
  \begin{aligned}
    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
    =\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} \; \frac{1}{{\cal{N}}^2} \left\{\sum_{{\mathbf s}_1,{\mathbf s}_2} (\mathbf s_1\cdot\mathbf s_2)^p\right\}
    =\frac1{N^p} \lim_{n\to0} \left\{\sum_{{\mathbf s}_1,{\mathbf s}_2,\ldots,\mathbf s_n} (\mathbf s_1\cdot\mathbf s_2)^p\right\}
  \end{aligned}
\end{equation}
The $(n-2)$ extra replicas providing the normalization.
Replacing the sums over stationary points with integrals over the Kac--Rice measure, the average over disorder, and again, for given energy and index, gives
\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}{\int\overline{\prod_a^n d\nu(\mathbf s_a)}\,\left(\frac{\mathbf s_1\cdot\mathbf s_2}N\right)^p} \\
    &=\lim_{n\to0}{\int D[C,R,D] \,
    \left(C_{12}\right)^p\; } e^{nN\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^{nN\Sigma[C,R,D]}
  \end{aligned}
\end{equation}
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}

One can construct a schematic 2RSB model from two 1RSB models.
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 landscape of the pure models is much
simpler than that of the mixed because, in these models, fixing the stability
$\mu$ is equivalent to fixing the energy: $\mu=pE$. This implies that at each
energy level there is only one type of stationary point. Therefore, for the
pure models our formulas for the complexity and its Legendre transforms are
functions of one variable only, $E$, and each instance of $\mu^*$ inside mus be
replaced with $pE$.

In the joint model, we wish to fix the total energy, not the energies of the
individual two models. Therefore, we insert a $\delta$-function containing
$(E_1+E_2)-E$ and integrate over $E_1$ and $E_2$. This results in a joint
complexity (and Legendre transform)
\begin{eqnarray}
    e^{N\Sigma(E)}&=&\int dE_1\, dE_2\, d\lambda \, e^{N[ \Sigma_1(E_1) + \Sigma_2(E_2) + O(\varepsilon) -\lambda N [(E_1+E_2)-E]}\nonumber \\
    e^{NG(\hat \beta)}&=&\int dE\, dE_1\, dE_2\, d\lambda\, 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$. One adds to the Hamiltonian a random term $\varepsilon \tilde H_p = \varepsilon \sum_{i_1,...,i_p} \tilde J_{i_1,...,i_p} s_{i_1}...s_{i_p}$, where the $\tilde J$ are 
random Gaussian uncorrelated with the $J$'s. 
The response to these is:
\begin{equation}
  \begin{aligned}
  & \overline{  \frac{\partial \langle \tilde H_p \rangle_{\tilde J}  }  {\partial \varepsilon} } 
  %  \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}\frac1{N^p}\sum_{i_1\cdots i_p}\frac\partial{\partial J^{(p)}_{i_1\cdots i_p}}
  %  \int\left(\prod_a^nd\nu(\mathbf s_a)\right)\,s^1_{i_1}\cdots s^1_{i_p} \\
    & =\lim_{n\to0}\int\left(\prod_a^nd\nu(\mathbf s_a)\right)\sum_b^n\left[
      \hat\beta\left(\frac{\mathbf s_1\cdot\mathbf s_b}N\right)^p+
        p\left(-i\frac{\mathbf s_1\cdot\hat{\mathbf s}_b}N\right)\left(\frac{\mathbf s_1\cdot\mathbf s_b}N\right)^{p-1}
      \right]
  \end{aligned}
\end{equation}
Taking the average of this expression over disorder and averaging over the equivalent replicas in the integral gives, similar to before,
\begin{equation}
  \begin{aligned}
   \overline{  \frac{\partial \langle \tilde H_p \rangle_{\tilde J}  }  {\partial \varepsilon} } 
  %  \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}\int D[C,R,D]\,\frac1n\sum_{ab}^n(\hat\beta C_{ab}^p+pR_{ab}C_{ab}^{p-1})e^{nN\Sigma[C,R,D]}\\
    &=\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
perturbations to the variance of the tensors $J$. This can be found by taking
the expression for the complexity and inserting the dependence of $f$ on the
coefficients $a_p$, then differentiating:
\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}{\partial 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{Maimbourg_2016_Solution}. 

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

\printbibliography

\end{document}