Lots of rewriting for making it nice and clean.

7 files changed, 447 insertions, 439 deletions

diff --git a/figs/24_opt_d(x).pdf b/figs/24_opt_d(x).pdf
index 886fc65..724656f 100644
--- a/figs/24_opt_d(x).pdf
+++ b/figs/24_opt_d(x).pdf
diff --git a/figs/24_opt_q(x).pdf b/figs/24_opt_q(x).pdf
index da265b6..85e6a64 100644
--- a/figs/24_opt_q(x).pdf
+++ b/figs/24_opt_q(x).pdf
diff --git a/figs/24_opt_r(x).pdf b/figs/24_opt_r(x).pdf
index eef9866..8a5c648 100644
--- a/figs/24_opt_r(x).pdf
+++ b/figs/24_opt_r(x).pdf
diff --git a/figs/24_opt_xMax.pdf b/figs/24_opt_xMax.pdf
index 1035b6b..023c864 100644
--- a/figs/24_opt_xMax.pdf
+++ b/figs/24_opt_xMax.pdf
diff --git a/figs/316_detail.pdf b/figs/316_detail.pdf
index 360efdf..6372c9c 100644
--- a/figs/316_detail.pdf
+++ b/figs/316_detail.pdf
diff --git a/frsb_kac-rice.bib b/frsb_kac-rice.bib
index e09dc8a..b93932f 100644
--- a/frsb_kac-rice.bib
+++ b/frsb_kac-rice.bib
@@ -12,6 +12,20 @@
  doi = {10.1103/physrevlett.126.028001}
 }
 
+@article{Altieri_2021_Properties,
+ author = {Altieri, Ada and Roy, Felix and Cammarota, Chiara and Biroli, Giulio},
+ title = {Properties of Equilibria and Glassy Phases of the Random {Lotka}-{Volterra} Model with Demographic Noise},
+ journal = {Physical Review Letters},
+ publisher = {American Physical Society (APS)},
+ year = {2021},
+ month = {6},
+ number = {25},
+ volume = {126},
+ pages = {258301},
+ url = {https://doi.org/10.1103%2Fphysrevlett.126.258301},
+ doi = {10.1103/physrevlett.126.258301}
+}
+
 @article{Annibale_2003_Supersymmetric,
  author = {Annibale, Alessia and Cavagna, Andrea and Giardina, Irene and Parisi, Giorgio},
  title = {Supersymmetric complexity in the {Sherrington}-{Kirkpatrick} model},
@@ -496,6 +510,21 @@ stochastic localization},
  doi = {10.1103/physrevx.9.011049}
 }
 
+@inproceedings{Maillard_2020_Landscape,
+ author = {Maillard, Antoine and Ben Arous, Gérard and Biroli, Giulio},
+ title = {Landscape Complexity for the Empirical Risk of Generalized Linear Models},
+ publisher = {PMLR},
+ year = {2020},
+ month = {7},
+ volume = {107},
+ pages = {287--327},
+ url = {https://proceedings.mlr.press/v107/maillard20a.html},
+ booktitle = {Proceedings of The First Mathematical and Scientific Machine Learning Conference},
+ editor = {Lu, Jianfeng and Ward, Rachel},
+ pdf = {http://proceedings.mlr.press/v107/maillard20a/maillard20a.pdf},
+ series = {Proceedings of Machine Learning Research}
+}
+
 @article{Maimbourg_2016_Solution,
  author = {Maimbourg, Thibaud and Kurchan, Jorge and Zamponi, Francesco},
  title = {Solution of the Dynamics of Liquids in the Large-Dimensional Limit},
diff --git a/frsb_kac-rice.tex b/frsb_kac-rice.tex
index 7a3d798..05eccb3 100644
--- a/frsb_kac-rice.tex
+++ b/frsb_kac-rice.tex
@@ -18,43 +18,50 @@
 \addbibresource{frsb_kac-rice.bib}
 
 \begin{document}
+
 \title{
-  How to count in hierarchical landscapes: a `full' solution to mean-field complexity
+  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.
+    We derive the general solution for the computation of stationary 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
+first applications of a replica symmetry breaking (RSB) 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.
+program of computing the number of stationary points---minima, saddle points,
+and maxima---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} and for similar
+energy functions inspired by molecular biology, evolution, and machine learning
+\cite{Maillard_2020_Landscape, Ros_2019_Complex, Altieri_2021_Properties}.  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 stationary points 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,
+another  transition to  a full RSB (FRSB) 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,
+\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
@@ -64,7 +71,7 @@ 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
+temperature, the question is posed as to`on what side of the 1RSB-FRSB
 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
@@ -73,24 +80,27 @@ 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.
+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,' the lowest energy level that good algorithms can
+expect to reach. This notion was introduced  \cite{Cugliandolo_1993_Analytical}
+in the context of the pure $p$-spin models, 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 complicated 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.
 
 The format of this paper is as follows. In \S\ref{sec:model}, we introduce the
-mean-field model we study, the mixed $p$-spin spherical model. In
-\S\ref{sec:equilibrium} we review details of the equilibrium solution that will
-be relevant in our study of the landscape complexity. In \S\ref{sec:complexity}
-we derive a generic form for the complexity of the model. In \S\ref{sec:ansatz}
-we make and review the hierarchical replica symmetry breaking ansatz used to
-solve the complexity. In \S\ref{sec:supersymmetric} we write down the solution
-in a specific and limited regime, which is nonetheless helpful as it gives a
+mean-field model of study, the mixed $p$-spin spherical model. In
+\S\ref{sec:equilibrium} we review details of the equilibrium solution that are
+relevant to our study of the landscape complexity. In \S\ref{sec:complexity} we
+derive a generic form for the complexity. In \S\ref{sec:ansatz} we make and
+review the hierarchical replica symmetry breaking ansatz used to solve the
+complexity. In \S\ref{sec:supersymmetric} we write down the solution in a
+specific and limited regime, which is nonetheless helpful as it gives a
 foothold for numerically computing the complexity everywhere else.
 \S\ref{sec:frsb} explains aspects of the solution specific to the case of full
 RSB, and derives the RS--FRSB transition line. \S\ref{sec:examples} details the
@@ -108,123 +118,138 @@ For definiteness, we consider the mixed $p$-spin spherical model, whose Hamilton
 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 energy is typically extensive. The overbar will always denote an average
+over the coefficients $J$. 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}
+In particular, one finds
+\begin{equation} \label{eq:covariance}
   \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}.
+One needn't start with a Hamiltonian like
+\eqref{eq:hamiltonian}, defined as a series: instead, the covariance rule
+\eqref{eq:covariance} can be specified 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
+The family of mixed $p$-spin models may be thought of as the most general
+models 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
+For reasons that will become clear in \S\ref{subsec:hess}, we refer to $\mu$
+as the \emph{stability parameter}.  At any stationary point, the gradient and
+Hessian are given by
 \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
+where $\partial=\frac\partial{\partial\mathbf s}$ always. An important
 observation was made by Bray and Dean  \cite{Bray_2007_Statistics} that
-gradient and Hessian are independent for  random Gaussian disorder.  The
+gradient and Hessian are independent for Gaussian random functions.  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).
+for any function of the gradient and one for any function of the Hessian, and
+in particular the number of negative eigenvalues, the index $\mathcal I$ of the
+saddle (see Fyodorov \cite{Fyodorov_2007_Replica} for a detailed discussion).
 
 \section{Equilibrium}
 \label{sec: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
+Here we review the equilibrium solution, which has been studied in detail
+\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$
+Once $n$ replicas are introduced to treat the logarithm, the fields $\mathbf
+s_a$ can be replaced with the new $n\times n$ matrix field $Q_{ab}\equiv(\mathbf
+s_a\cdot\mathbf s_b)/N$. This yields for the free energy
 \begin{equation}
+  \beta F=-1-\ln2\pi-\frac12\lim_{n\to0}\frac1n\left(\beta^2\sum_{ab}^nf(Q_{ab})+\ln\det Q\right)
+\end{equation}
+which must be evaluated at the $Q$ which maximizes this expression and whose
+diagonal is one. The solution is generally a hierarchical matrix \emph{à la}
+Parisi. Each row of the matrix is the same up to permutation of their elements.
+The so-called $k$RSB solution has $k+2$ different values in each row: $n-x_1$
+of those entries are $q_0$, $x_1-x_2$ of those entries are $q_1$, and so on
+until $x_k-1$ entries of $q_k$, and of course one entry of $1$ (the diagonal).
+Given such a matrix, there are standard ways of producing the sum and
+determinant that appear in the free energy. These formulas are, for an
+arbitrary $k$RSB matrix $A$ with $a_d$ on its diagonal (recall $q_d=1$),
+\begin{equation} \label{eq:replica.sum}
+  \lim_{n\to0}\frac1n\sum_{ab}^nA_{ab}
+  =a_d-\sum_{i=0}^k(x_{i+1}-x_i)a_i
+\end{equation}
+\begin{equation} \label{eq:replica.logdet}
   \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)
+    \lim_{n\to0}\frac1n\ln\det A
+    &=
+    \frac{a_0}{a_d-\sum_{i=0}^k(x_{i+1}-x_i)a_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
+      a_d-\sum_{i=0}^{k}(x_{i+1}-x_i)a_i
+    \right]\\
+    &\hspace{10pc}-\sum_{j=1}^k(x_j^{-1}-x_{j+1}^{-1})\log\left[
+        a_d-\sum_{i=j}^{k}(x_{i+1}-x_i)a_i-x_ja_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
+where $x_0=0$ and $x_{k+1}=1$. Once these are substituted into the free energy,
+optimizing the anstaz is equivalent to maximizing $F$ with respect to the
+$q_0,\ldots,q_k$ and $x_1,\ldots,x_k$.
+
+The free energy can also be written in a functional form, which is necessary
+for working with the solution in the limit $k\to\infty$, the so-called full
+replica symmetry breaking (FRSB). If $P(q)$ is the probability distribution for
+elements $q$ in a row of the matrix, then define $\chi(q)$ by
 \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}
+  \chi(q)=\int_q^1dq'\,\int_0^{q'}dq''\,P(q'')
 \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
+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$. The
+function $\chi$ turns out to have an interpretation as the spectrum of the
+hierarchical matrix $Q$.  Using standard arguments, the free energy can be
+written as a functional over $\chi$ as
+\begin{equation}
+  \beta F=-1-\ln2\pi-\frac12\int_0^1dq\,\left(\beta^2f''(q)\chi(q)+\frac1{\chi(q)}\right)
+\end{equation}
+which must be maximized with respect to $\chi$ given the constraints outlined
+above.
+
+In our study of the landscape, the free energy will not be directly relevant
+anywhere except at the ground state, when the temperature is zero or
+$\beta\to\infty$. Here, the measure will be concentrated in the lowest minima,
+and the average energy $\langle
+E\rangle_0=\lim_{\beta\to\infty}\frac\partial{\partial\beta}\beta F$ will
+correspond to the ground state energy $E_0$.  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, one can show after some algebra that
 \begin{equation} \label{eq:ground.state.free.energy}
-  \lim_{\beta\to\infty}\tilde\beta F
+  \tilde\beta\langle E\rangle_0=\tilde\beta\lim_{\beta\to\infty}\frac{\partial(\beta F)}{\partial\beta}
   =-\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)
+    \tilde\beta^2\sum_{ab}^nf(\tilde Q_{ab})+\ln\det(\tilde\beta z^{-1}\tilde Q+I)
   \right)
 \end{equation}
-In the continuum case, this is
+where $\tilde Q$ is a $(k-1)$RSB matrix with entries $\tilde q_1=\lim_{\beta\to\infty}q_1$, \dots, $\tilde q_{k-1}=\lim_{\beta\to\infty}q_{k-1}$ parameterized by $\tilde x_1,\ldots,\tilde x_{k-1}$.
+This is a $k-1$ RSB ansatz whose spectrum in the determinant is scaled by
+$\tilde\beta z^{-1}$ and shifted by 1, with effective temperature
+$\tilde\beta$, and an extra term. 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 \langle E\rangle_0
+  =-\frac12z\tilde\beta f'(1)-\frac12\int_0^1 dq\left(
     \tilde\beta^2f''(q)\tilde\chi(q)+\frac1{\tilde\chi(q)+\tilde\beta z^{-1}}
   \right)
 \end{equation}
+where $\tilde\chi$ is bound by the same constraints as before.
 
 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,
@@ -239,7 +264,7 @@ each stationary point also has no width and therefore overlap one with itself.
 \label{sec: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
+which integrates 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}
@@ -249,7 +274,7 @@ the gradient multiplied by the absolute value of the determinant
 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
+than fix the index directly, we fix the trace of the Hessian, 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}
@@ -259,9 +284,9 @@ to within order one. Inserting these $\delta$-functions, we arrive at
     &\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:
+This number will typically be exponential in $N$. In order to find the typical
+count when disorder is averaged, we 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}
@@ -270,7 +295,7 @@ If one averages over $\mathcal N$ and afterward takes its logarithm, one arrives
   \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
+The annealed complexity 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
@@ -279,7 +304,7 @@ 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}.
+function.
 
 A sometimes more illuminating quantity to consider is the Legendre transform $G$ of the complexity:
 \begin{equation}
@@ -326,6 +351,7 @@ independence, we may write
 which simplifies matters. The average of the two factors may now be treated separately.
 
 \subsubsection{The Hessian factors}
+\label{subsec:hess}
 
 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
@@ -341,9 +367,10 @@ Therefore in that limit its spectrum is given by the Wigner semicircle with radi
   \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
+semicircle shifted by $\mu$, or $\rho(\lambda+\mu)$.  The stability parameter
+$\mu$ thus fixes the center of the spectrum of the Hessian. The value
+$\mu_m=\sqrt{4f''(1)}$ is a kind of threshold. When $\mu$ is taken to be within
+the range $\pm\mu_m$, the critical points have index density
 \begin{equation}
   \mathcal I(\mu)=\int_0^\infty d\lambda\,\rho(\lambda+\mu)
   =\frac12-\frac1\pi\left[
@@ -352,8 +379,10 @@ $\pm\mu_m\equiv\pm\sqrt{4f''(1)}$, the critical points have index density
   \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.
-
+$\mu-\mu_m$. When $\mu=\mu_m$, the critical points are marginal minima, with
+flat directions in their spectrum. This property of $\mu$ is why we've named it
+the stability parameter: it governs the stability of stationary points, and for
+unstable ones it governs their index.
 
 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
@@ -366,33 +395,23 @@ 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| \\
+    &=\frac1N\overline{\ln|\det\operatorname{Hess}H(s,\mu)|}
+    =\int d\lambda\,\rho(\lambda+\mu)\ln|\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)
+    -\ln\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}
+the value of the stability $\mu$ in all replicas to the value $\mu^*$, and
+therefore the index of saddles in all replicas as well.
 
 \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
+The $\delta$-functions in the remaining factor 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)) \\
@@ -434,7 +453,6 @@ on the Hamiltonian, and since everything is Gaussian this gives
       \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 &&
@@ -455,16 +473,18 @@ s}_a$ to these three matrices, we arrive at the form for the complexity
         \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}
+    +\frac12\ln\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.
+of this expression which minimize the complexity. Note that one cannot
+\emph{minimize} the complexity with respect to these parameters: there is no
+pure variational problem here. 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
@@ -478,7 +498,7 @@ transform
         \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}
+    +\frac12\ln\det\begin{bmatrix}C&iR\\iR&D\end{bmatrix}
     \right)
   \end{aligned}
 \end{equation}
@@ -495,7 +515,7 @@ given by
         \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}
+    +\frac12\ln\det\begin{bmatrix}C&iR\\iR&D\end{bmatrix}
     \right)
   \end{aligned}
 \end{equation}
@@ -506,30 +526,30 @@ 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}
 \label{sec: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
+be hierarchical matrices in 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
+matrix products and Hadamard products. In particular, the determinant of the block matrix can be written as a determinant of a product,
+\begin{equation}
+  \ln\det\begin{bmatrix}C&iR\\iR&D\end{bmatrix}=\ln\det(CD+R^2)
+\end{equation}
+The is straightforward to write down at $k$RSB, since the product and sum of
+the hierarchical matrices is still a hierarchical matrix. Recall the
+standard formulas for the product $C=AB$,
+\begin{align} \label{eq:replica.prod}
+  c_d&=a_db_d-\sum_{j=0}^k(x_{j+1}-x_j)a_jb_j \\
+  c_i&=b_da_i+a_db_i-\sum_{j=0}^{i-1}(x_{j+1}-x_j)a_jb_j+(2x_{i+1}-x_i)a_ib_i
+  -\sum_{j=i+1}^k(x_{j+1}-x_j)(a_ib_j+a_jb_i)
+\end{align}
+and for the sum $C=A+B$, $c_d=a_d+b_d$ and $c_i=a_i+b_i$. Using these, one can
+write down the hierarchical matrix $CD+R^2$, and then compute the $\ln\det$
+using the previous formula \eqref{eq:replica.logdet}.
+
+The extremal conditions are given by differentiating the complexity with respect to its parameters, yielding
 \begin{align}
   0&=\frac{\partial\Sigma}{\partial\hat\mu}
   =\frac12(c_d-1) \\
@@ -551,6 +571,13 @@ where $\odot$ denotes the Hadamard product, or the componentwise product. Equati
 \begin{equation} \label{eq:D.solution}
   D=f'(C)^{-1}-RC^{-1}R
 \end{equation}
+To these conditions must be added the addition condition that $\Sigma$ is extremal with respect to $x_1,\ldots, x_k$. There is no better way to enforce this condition than to directly differentiate $\Sigma$ with respect to the $x$s, and we have
+\begin{equation} \label{eq:cond.x}
+  0=\frac{\partial\Sigma}{\partial x_i}\qquad 1\leq i\leq k
+\end{equation}
+These stationary conditions are the most numerically taxing, because when the
+formulas above are chained together and substituted into the complexity it
+results in a complicated formula.
 
 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
@@ -571,9 +598,9 @@ When $\mu^*<\mu_m$ and the critical points are saddles, it implies
 \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
+inverses, so that the saddle conditions can be expressed using
+\eqref{eq:replica.prod} alone. 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 \\
@@ -581,25 +608,25 @@ as
 \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}
 \label{sec:supersymmetric}
 
 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
+absolute value of the determinant is neglected and the trace of the Hessian is
+not fixed. This supersymmetry has been studied in great detail in the
+complexity of the Thouless--Anderson--Palmer (TAP) free energy
+\cite{Annibale_2003_The, Annibale_2003_Supersymmetric,
+Annibale_2004_Coexistence}. When the absolute value is dropped, 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 to a specific
+line in the energy and stability plane.
+
+The identity $D=\hat\beta R$ 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}
@@ -609,7 +636,7 @@ Substituting \eqref{eq:cond.r} for the factor $(CD+R^2)^{-1}R$, we find substant
 \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.
+supersymmetric ansatz is equivalent to a \emph{diagonal} ansatz for both $R$ and $D$.
 
 Supersymmetry has further implications.
 Equations \eqref{eq:cond.r} and \eqref{eq:cond.d} can be combined to find
@@ -620,26 +647,26 @@ 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
+Understanding that $R$ is diagonal, we find
 \begin{equation}
   \mu^*=\frac1{r_d}+r_df''(1)
 \end{equation}
-which is precisely the condition \eqref{eq:mu.minima}. Therefore, \emph{the
+which is precisely the condition \eqref{eq:mu.minima} for dominant 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
+Inserting the supersymmetric ansatz $D=\hat\beta R$ and $R=r_dI$, one gets for the complexity
 \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\hat\beta r_df'(1)+\frac12r_d^2f''(1)+\frac12\ln 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)
+      +\frac12\lim_{n\to0}\frac1n\left(\hat\beta^2\sum_{ab}f(C_{ab})+\ln\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
@@ -650,31 +677,36 @@ $\Sigma(E_0,\mu^*)=0$, gives
   \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)
+      +\ln\det(\hat\beta r_d^{-1} C+I)
   \right)
 \end{equation}
-which is precisely \eqref{eq:ground.state.free.energy} with $r_d=z$,
+which is precisely the ground state energy predicted by the equilibrium
+solution\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
+{\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
+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
+x_1,\ldots,\tilde x_{k-1}$, and $c_1,\ldots,c_{k-1}$ for the
 complexity in the ground state are
 \begin{align}\label{eq:equilibrium.complexity.map}
   \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
+  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}
+Unlike the case for the TAP complexity, this correspondence between complexity
+and equilibrium solutions only exists at the ground state. We will see in our
+examples in \S\ref{sec:examples} that there appears to be little correspondence
+between these parameters away from the ground state.
 
 The supersymmetric solution produces the correct complexity for the ground
 state and for a class of minima. Moreover, it produces the correct parameters
@@ -691,13 +723,35 @@ sufficiently close to the correct answer. This is the strategy we use in
 \section{Full replica symmetry breaking}
 \label{sec:frsb}
 
-This reasoning applies equally well to FRSB systems.
+This reasoning applies equally well to FRSB systems. 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}
-  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}
+  C\;\leftrightarrow\;[c_d, c(x)]
+  &&
+  R\;\leftrightarrow\;[r_d, r(x)]
+  &&
+  D\;\leftrightarrow \;d_d, d(x)]
 \end{align}
-where the product of two hierarchical matrix results in
+the complexity becomes
+\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\int_0^1dx\left[
+      \hat\beta^2f(c(x))+(2\hat\beta r(x)-d(x))f'(c(x))
+      +r(x)^2f''(c(x))
+      \right]
+      +\frac12\lim_{n\to0}\frac1n\ln\det(CD+R^2)
+  \end{aligned}
+\end{equation}
+The determinant takes more work to treat.
+Define for any function of $x$
+\begin{equation}
+  \langle a\rangle=\int_0^1dx\,a(x)
+\end{equation}
+In general, the product of hierarchical matrices becomes
 \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)
@@ -707,62 +761,76 @@ where the product of two hierarchical matrix results in
     b(x)-b(y)
   \big)
 \end{align}
-and
-\begin{equation}
-  \langle a\rangle=\int_0^1dx\,a(x)
-\end{equation}
+and the $\ln\det$ becomes, for the hierarchical matrix $A=CD+R^2$,
+\begin{equation} \label{eq:replica.det.cont}
+  \lim_{n\to0}\frac1n\ln\det A
+  =\ln(a_d-\langle a\rangle)
+  +\frac{a(0)}{a_d-\langle a\rangle}
+  -\int_0^1\frac{dx}{x^2}\ln\left(\frac{a_d-\langle a\rangle+\int_0^xdy\,a(y)-xa(x)}{a_d-\langle a\rangle}\right)
+\end{equation}
+The saddle point equations then take the form
+\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}
 
 \subsection{Supersymmetric complexity}
 
-Using standard
-manipulations (Appendix B), one finds also a continuous version of the
+Using standard manipulations, 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\left(\hat\beta r_df'(1)+r_d^2f''(1)+\ln 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
+where $\chi(q)=\int_1^qdq'\int_0^{q'}dq''\,P(q)$ for $P(q)$ the distribution of elements in a row of $C$, as in the equilibrium case. Like in the equilibrium case, $\chi$ must be concave, monotonically decreasing, and have $\chi(1)=0$, $\chi'(1)=-1$.
+
+First, we use this solution to inspect the ground state of a full RSB system.
+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
+over all functions $\chi$. This gives
 \begin{equation}
-  \chi(q)=\frac1{\hat\beta}\left(f''(q)^{-1/2}-r_d\right)
+  \chi_0(q\mid\hat\beta,r_d)=\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.
-
+Satisfying the boundary condition $\chi(1)=0$ requires $r_d=f''(1)^{-1/2}$ and $\hat\beta=\frac12f'''(1)/f''(1)^{3/2}$.
+This in turn implies $\mu^*=\frac1{r_d}+f''(1)r_d=\sqrt{4f''(1)}=\mu_m$.
+Therefore, the FRSB ground state is always marginal, as excepted. 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.
+\emph{no} energies at which minima are most numerous; saddles always dominate.
 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)
+Away from the ground state, this expression still correctly counts a class of
+non-dominant minima. However, like in the equilibrium solution, the function
+$\chi$ which produces an extremal value is not smooth in the entire range
+$[0,1]$, but adopts a piecewise form
+\begin{equation}
+  \chi(q)=\begin{cases}
+    \chi_0(q\mid\hat\beta,r_d) & q\leq q_\textrm{max} \\
+    1-q & \text{otherwise}
+  \end{cases}
+\end{equation}
+With this ansatz, the complexity must be extremized with respect to $r_d$ and
+$\hat\beta$, while simultaneously ensuring that $q_\textrm{max}$ is such that
+$\chi(q)$ is continuous, that is, that
+$\chi_0(q_\textrm{max}\mid\hat\beta,r_d)=1-q_\textrm{max}$. The significance of
+the minima counted by this method is unclear, but they do represent a nodal
+line in the off-diagonal parts of $R$ and $D$. Since, as usual, $\chi(q)$ is
+related to $c(x)$ by $-\chi'(c(x))=x$, there is a corresponding $x_\textrm{max}$ given by
+\begin{equation}
+  x_\textrm{max}=-\chi'(q_\textrm{max})=\frac1{2\hat\beta}\frac{f'''(q_\textrm{max})}{f''(q_\textrm{max})^{3/2}}
 \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}
@@ -771,76 +839,75 @@ 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$
+expected to be small. In particular, we expect the functions $c(x)$, $r(x)$, and $d(x)$
 to approach zero at the transition, and moreover take the form
-\begin{equation}
+\begin{align}
   c(x)=\begin{cases}\bar cx&x\leq x_\mathrm{max}\\\bar cx_\mathrm{max}&\text{otherwise}\end{cases}
-\end{equation}
+                           &&
+  r(x)=\begin{cases}\bar rx&x\leq x_\mathrm{max}\\\bar rx_\mathrm{max}&\text{otherwise}\end{cases}
+                           &&
+  d(x)=\begin{cases}\bar dx&x\leq x_\mathrm{max}\\\bar dx_\mathrm{max}&\text{otherwise}\end{cases}
+\end{align}
 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.
+numerics of \S\ref{sec:examples}
 
 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
+\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 linear 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)
+    -\hat\beta-\frac1{f'(1)+f''(0)}\left(r_d(f''(0)+f''(1))-\mu^*\right)+O(x_\mathrm{max})
   \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)
+    &=-2r_d\frac{\bar r}{\bar c}-\frac1{f'(1)}(r_d^2f''(0)+d_d(f'(1)+f''(0))-1)+O(x_\mathrm{max})
   \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
+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, using \eqref{eq:replica.det.cont} to
+treat the determinant. 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) \\
+  =\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\ln(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
+The spectrum of the Hessian of $\Sigma$ evaluated at the replica symmetric
+solution with $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}
+\begin{equation} \label{eq:mu.transition}
   \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.
+that this is the line of stability for the replica symmetric solution. The
+numerics in \S\ref{sec:examples} bear this out.
 
 
 \section{General solution: examples}
@@ -858,9 +925,9 @@ state. With these parameters in hand, small steps are then made in energy $E$
 or stability $\mu$, after which known values are used as the initial condition
 for a saddle-finding problem. In this section, we use this basic numeric idea
 to map out the complexity for two representative examples: a model with a
-$2RSB$ equilibrium ground state and therefore $1RSB$ complexity in its
-vicinity, and a model with a $FRSB$ equilibrium ground state, and therefore
-$FRSB$ complexity as well.
+2RSB equilibrium ground state and therefore 1RSB complexity in its
+vicinity, and a model with a FRSB equilibrium ground state, and therefore
+FRSB complexity as well.
 
 \subsection{1RSB complexity}
 
@@ -889,7 +956,7 @@ 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
+topology of the energy function might be relevant to where this algorithmic threshold
 lies. For the $3+16$ model at hand, $E_\mathrm{alg}=-1.275\,140\,128\ldots$.
 
 In this model, the RS complexity gives an inconsistent answer for the
@@ -900,12 +967,21 @@ 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.287\,605\,527\ldots$, which is very slightly greater than $E_0$. Saddles
 become dominant over minima at a higher energy $E_\mathrm{th}=-1.287\,575\,114\ldots$.
-The 1RSB complexity transitions to a RS description for dominant stationary points at an energy
-$E_1=-1.273\,886\,852\ldots$. The highest energy for which the 1RSB description exists is $E_\mathrm{max}=-0.886\,029\,051\ldots$
-
-The complexity along interesting trajectories is plotted in Fig.~\ref{fig:2rsb.complexity}.
-
-
+The 1RSB complexity transitions to a RS description for dominant stationary
+points at an energy $E_1=-1.273\,886\,852\ldots$. The highest energy for which
+the 1RSB description exists is $E_\mathrm{max}=-0.886\,029\,051\ldots$
+
+The complexity as a function of energy difference from the ground state is
+plotted in Fig.~\ref{fig:2rsb.complexity}. In that figure, the complexity is
+plotted for dominant minima and saddles, marginal minima, and supersymmetric
+minima. A contour plot of the complexity as a function of energy $E$ and
+stability $\mu$ is shown in Fig.~\ref{fig:2rsb.contour}. That plot also shows
+the RS--1RSB transition line in the complexity. For minima, the complexity does
+not inherit a 1RSB description until the energy is with in a close vicinity of
+the ground state. On the other hand, for high-index saddles the complexity
+becomes described by 1RSB at quite high energies. This suggests that when
+sampling a landscape at high energies, high index saddles may show a sign of
+replica symmetry breaking when minima do not.
 
 \begin{figure}
   \includegraphics{figs/316_complexity.pdf}
@@ -930,15 +1006,20 @@ The complexity along interesting trajectories is plotted in Fig.~\ref{fig:2rsb.c
 
   \caption{
     Complexity of the $3+16$ model in the energy $E$ and stability $\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.
-  }
+    plane. The right shows a detail of the left.
+  } \label{fig:2rsb.contour}
 \end{figure}
 
+Fig.~\ref{fig:2rsb.phases} shows a different detail of the complexity in the
+vicinity of the ground state, now as functions of the energy difference and
+stability difference from the ground state. Several of the landmark energies
+described above are plotted, alongside the boundaries between the `phases.'
+Though $E_\mathrm{alg}$ looks quite close to the energy at which dominant
+saddles transition from 1RSB to RS, they differ by roughly $10^{-3}$, as
+evidenced by the numbers cited above. Likewise, though $\langle E\rangle_1$
+looks very close to $E_\mathrm{max}$, where the 1RSB transition line
+terminates, they too differ.
+
 \begin{figure}
   \centering
   \includegraphics{figs/316_detail.pdf}
@@ -958,7 +1039,7 @@ The complexity along interesting trajectories is plotted in Fig.~\ref{fig:2rsb.c
     E\rangle_2$ and $\langle E\rangle_1$, respectively. Though the figure is
     suggestive, $E_\mathrm{alg}$ lies slightly below the termination of the RS
     -- 1RSB transition line.
-  }
+  } \label{fig:2rsb.phases}
 \end{figure}
 
 Fig.~\ref{fig:2rsb.comparison} shows the saddle parameters for the $3+16$
@@ -975,7 +1056,8 @@ complexity meets the zero complexity line. Here, $d_1$ diverges like
   d_1=-\left(\frac1{f'(1)}-(d_d+r_d^2)\right)(1-x_1)^{-1}+O(1)
 \end{equation}
 while $x_1$ and $q_1$ both go to one. Note that this is the only place along
-the phase boundary where $q_1$ goes to one.
+the phase boundary where $q_1$ goes to one. The significance of this critical
+point in the complexity of high-index saddles in worth further study.
 
 \begin{figure}
   \centering
@@ -989,13 +1071,51 @@ the phase boundary where $q_1$ goes to one.
   \includegraphics{figs/316_comparison_b.pdf}
 
   \caption{
+    Comparison of the saddle point parameters for the $3+16$ model along
+    different trajectories in the energy and stability space, and with the
+    equilibrium values (when they exist) at the same value of $\langle
+    E\rangle$.
   } \label{fig:2rsb.comparison}
 \end{figure}
 
 
 \subsection{Full RSB complexity}
 
-$\beta_\infty=1$ $\langle E\rangle_\infty=-0.53125$
+On the other hand, if the covariance $f$ is chosen to be concave, then one develops FRSB in equilibrium. To this purpose, we choose
+\begin{equation}
+  f(q)=\frac12\left(q^2+\frac1{16}q^4\right)
+\end{equation}
+also studied before in equilibrium \cite{Crisanti_2004_Spherical, Crisanti_2006_Spherical}. Because the ground state is FRSB, for this model
+\begin{equation}
+  E_0=E_\mathrm{alg}=E_\mathrm{th}=-\int_0^1dq\sqrt{f''(q)}=-1.059\,384\,319\ldots
+\end{equation}
+In the equilibrium solution, the transition temperature from RS to FRSB is $\beta_\infty=1$, with corresponding average energy $\langle E\rangle_\infty=-0.53125\ldots$.
+
+Along the supersymmetric line, the FRSB solution can be found in full, exact
+function form.  To treat the FRSB away from this line numerically, we result to
+finite $k$RSB approximations.  Since we are not trying to find the actual
+$k$RSB solution, but approximate the FRSB one, we drop the extremal condition
+\eqref{eq:cond.x} for $x_1,\ldots,x_k$ and instead set
+\begin{equation}
+  x_i=\left(\frac i{k+1}\right)x_\textrm{max}
+\end{equation}
+and extremize over $x_\textrm{max}$ alone. This dramatically simplifies the
+equations that must be solved to find solutions. In the results that follow, a
+20RSB solution is used to trace the dominant saddles and marginal minima, while
+a 5RSB solution is used to trace the (much longer) boundaries of the
+complexity.
+
+Fig.~\ref{fig:frsb.complexity} shows the complexity for this model as a
+function of energy difference from the ground state for several notable
+trajectories in the energy and stability plane. Fig.~\ref{fig:frsb.phases}
+shows these trajectories, along with the phase boundaries of the complexity in
+this plane. Notably, the phase boundary predicted by \eqref{eq:mu.transition}
+correctly predicts where all of the finite $k$RSB approximations terminate.
+Like the 1RSB model in the previous subsection, this phase boundary is oriented
+such that very few, low energy, minima are described by a FRSB solution, while
+relatively high energy saddles of high index are also. Again, this suggests
+that studying the mutual distribution of high-index saddle points might give
+insight into lower-energy symmetry breaking in more general contexts.
 
 \begin{figure}
   \centering
@@ -1010,13 +1130,25 @@ $\beta_\infty=1$ $\langle E\rangle_\infty=-0.53125$
   } \label{fig:frsb.complexity}
 \end{figure}
 
+
 \begin{figure}
   \centering
   \includegraphics{figs/24_phases.pdf}
   \caption{
+    `Phases' of the complexity for the $2+4$ model in the energy $E$ and
+    stability $\mu^*$ plane. The region shaded gray shows where the RS solution
+    is correct, while the region shaded red shows that where the FRSB solution
+    is correct. The white region shows where the complexity is zero.
   } \label{fig:frsb.phases}
 \end{figure}
 
+Fig.~\ref{fig:24.func} shows the value of $x_\textrm{max}$ along several
+trajectories of interest. Everywhere along the transition line,
+$x_\textrm{max}$ continuously goes to zero. Examples of our 20RSB
+approximations of the continuous functions $c(x)$, $r(x)$, and $d(x)$ are also
+shown. As expected, these functions approach linear ones as $x_\textrm{max}$
+goes to zero with finite slopes.
+
 \begin{figure}
   \raggedright
   \hspace{2em}\includegraphics{figs/24_opt_q(x).pdf}
@@ -1028,31 +1160,13 @@ $\beta_\infty=1$ $\langle E\rangle_\infty=-0.53125$
   \includegraphics{figs/24_opt_d(x).pdf}
 
   \caption{
+    $x_\textrm{max}$ as a function of $E$ for several trajectories of interest,
+    along with examples of the 20RSB approximations of the functions $c(x)$,
+    $r(x)$, and $d(x)$ along the dominant saddles. Colors of the approximate
+    functions correspond to the points on the $x_\textrm{max}$ plot.
   } \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:frsb.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}
@@ -1262,146 +1376,11 @@ A first and very important application of the method here is to perform the calc
 a clear understanding of what happens in a low-temperature realistic jamming dynamics \cite{Maimbourg_2016_Solution}. 
 
 \paragraph{Acknowledgements}
-J K-D and J K would like to thank Valentina Ros for helpful discussions.
+The authors would like to thank Valentina Ros for helpful discussions.
 
 \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}