path: root/topology.tex

\documentclass{SciPost}

% Prevent all line breaks in inline equations.
\binoppenalty=10000
\relpenalty=10000

\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath,latexsym,graphicx}
\usepackage[bitstream-charter]{mathdesign}
\usepackage[dvipsnames]{xcolor}
\usepackage{anyfontsize,authblk}
\usepackage{xfrac} % sfrac

\urlstyle{sf}

\fancypagestyle{SPstyle}{
\fancyhf{}
\lhead{\colorbox{scipostblue}{\bf \color{white} ~SciPost Physics }}
\rhead{{\bf \color{scipostdeepblue} ~Submission }}
\renewcommand{\headrulewidth}{1pt}
\fancyfoot[C]{\textbf{\thepage}}
}

% Fix \cal and \mathcal characters look (so it's not the same as \mathscr)
\DeclareSymbolFont{usualmathcal}{OMS}{cmsy}{m}{n}
\DeclareSymbolFontAlphabet{\mathcal}{usualmathcal}

\hypersetup{
    colorlinks,
    linkcolor={red!50!black},
    citecolor={blue!50!black},
    urlcolor={blue!80!black}
}

\author{\footnote{\url{jaron.kent-dobias@roma1.infn.it}}}
\affil{Istituto Nazionale di Fisica Nucleare, Sezione di Roma I, Italy}

\begin{document}

\pagestyle{SPstyle}

\begin{center}{\Large \textbf{\color{scipostdeepblue}{
  On the topology of solutions to random continuous constraint satisfaction problems\\
}}}\end{center}

\begin{center}
  \textbf{Jaron Kent-Dobias$^\star$}
\end{center}

\begin{center}
Istituto Nazionale di Fisica Nucleare, Sezione di Roma I, Italy
\\[\baselineskip]
$\star$ \href{mailto:jaron.kent-dobias@roma1.infn.it}{\small jaron.kent-dobias@roma1.infn.it}
\end{center}

\section*{\color{scipostdeepblue}{Abstract}}
\textbf{\boldmath{%
We consider the set of solutions to $M$ random polynomial equations of $N$
variables. Each equation has independent Gaussian coefficients and a target
value $V_0$, and their variables are restricted to the $(N-1)$-sphere. When
solutions exist, they form a manifold. We compute the average Euler
characteristic of this manifold in the limit of large $N$, and find different
behavior depending on the target value $V_0$, the ratio $\alpha=M/N$, and the
variances of the polynomial coefficients. We divide the behavior in four
phases: a connected phase, an onset phase, a shattered phase, and an
\textsc{unsat} phase. In the connected phase, the average characteristic is 2
and there is a single extensive connected component, while in the onset phase
the Euler characteristic is exponentially large in $N$. In the shattered phase
the characteristic remains exponentially large but subextensive components
appear, while in the \textsc{unsat} phase the manifold vanishes. When $M=1$
there is a correspondence between this problem and the topology of
energy level sets in the spherical spin glasses. We conjecture that
the transition from the onset to shattered phase corresponds to the asymptotic
limit of gradient descent from a random initial condition.
}
}

\vspace{\baselineskip}

%%%%%%%%%% BLOCK: Copyright information
% This block will be filled during the proof stage, and finilized just before publication.
% It exists here only as a placeholder, and should not be modified by authors.
\noindent\textcolor{white!90!black}{%
\fbox{\parbox{0.975\linewidth}{%
\textcolor{white!40!black}{\begin{tabular}{lr}%
  \begin{minipage}{0.6\textwidth}%
    {\small Copyright attribution to authors. \newline
    This work is a submission to SciPost Physics. \newline
    License information to appear upon publication. \newline
    Publication information to appear upon publication.}
  \end{minipage} & \begin{minipage}{0.4\textwidth}
    {\small Received Date \newline Accepted Date \newline Published Date}%
  \end{minipage}
\end{tabular}}
}}
}
%%%%%%%%%% BLOCK: Copyright information


%%%%%%%%%% TODO: LINENO
% For convenience during refereeing we turn on line numbers:
%\linenumbers
% You should run LaTeX twice in order for the line numbers to appear.
%%%%%%%%%% END TODO: LINENO

%%%%%%%%%% TODO: TOC 
% Guideline: if your paper is longer that 6 pages, include a TOC
% To remove the TOC, simply cut the following block
\vspace{10pt}
\noindent\rule{\textwidth}{1pt}
\tableofcontents
\noindent\rule{\textwidth}{1pt}
\vspace{10pt}
%%%%%%%%%% END TODO: TOC

\section{Introduction}

Constraint satisfaction problems seek configurations that simultaneously
satisfy a set of equations, and form a basis for thinking about problems as
diverse as neural networks \cite{Mezard_2009_Constraint}, granular materials
\cite{Franz_2017_Universality}, ecosystems \cite{Altieri_2019_Constraint}, and
confluent tissues \cite{Urbani_2023_A}. All but the last of these examples deal
with sets of inequalities, while the last considers a set of equality
constraints. Inequality constraints are familiar in situations like zero-cost
solutions in neural networks with ReLu activations and stable equilibrium in the
forces between physical objects. Equality constraints naturally appear in the
zero-gradient solutions to overparameterized smooth neural networks and in vertex models of tissues.

In such problems, there is great interest in characterizing structure in the
set of solutions, which can be influential in how algorithms behave when trying
to solve them \cite{Baldassi_2016_Unreasonable, Baldassi_2019_Properties,
Beneventano_2023_On}. Here, we show how topological information about
the set of solutions can be calculated in a simple model of satisfying random
nonlinear equalities. This allows us to reason about the connectivity of this
solution set. The topological properties revealed by this calculation yield
surprising results for the well-studied spherical spin glasses, where a
topological transition thought to occur at a threshold energy $E_\text{th}$
where marginal minima are dominant is shown to occur at a different energy
$E_\text{sh}$. We conjecture that this difference resolves an outstanding
problem in gradient descent dynamics in these systems.

We consider the problem of finding configurations $\mathbf x\in\mathbb R^N$
lying on the $(N-1)$-sphere $\|\mathbf x\|^2=N$ that simultaneously satisfy $M$
nonlinear constraints $V_k(\mathbf x)=V_0$ for $1\leq k\leq M$ and some
constant $V_0\in\mathbb R$. The nonlinear constraints are taken to be centered
Gaussian random functions with covariance
\begin{equation} \label{eq:covariance}
  \overline{V_i(\mathbf x)V_j(\mathbf x')}
  =\delta_{ij}f\left(\frac{\mathbf x\cdot\mathbf x'}N\right)
\end{equation}
for some choice of function $f$. When the covariance function $f$ is polynomial, the
$V_k$ are also polynomial, with a term of degree $p$ in $f$ corresponding to
all possible terms of degree $p$ in $V_k$. In particular, taking
\begin{equation}
  V_k(\mathbf x)
  =\sum_{p=0}^\infty\frac1{p!}\sqrt{\frac{f^{(p)}(0)}{N^p}}
  \sum_{i_1\cdots i_p}^NJ^{(k,p)}_{i_1\cdots i_p}x_{i_1}\cdots x_{i_p}
\end{equation}
with the elements of the tensors $J^{(k,p)}$ as independently distributed
unit normal random variables satisfies \eqref{eq:covariance}. The size of the
series coefficients of $f$ therefore control the variances in the coefficients
of random polynomial constraints. When $M=1$, this problem corresponds to the
level set of a spherical spin glass with energy density $E=V_0/\sqrt{N}$.


This problem or small variations thereof have attracted attention recently for
their resemblance to encryption, optimization, and vertex models of confluent
tissues \cite{Fyodorov_2019_A, Fyodorov_2020_Counting,
  Fyodorov_2022_Optimization, Tublin_2022_A, Vivo_2024_Random, Urbani_2023_A, Kamali_2023_Dynamical,
  Kamali_2023_Stochastic, Urbani_2024_Statistical, Montanari_2023_Solving,
Montanari_2024_On, Kent-Dobias_2024_Conditioning, Kent-Dobias_2024_Algorithm-independent}. In each of these cases, the authors studied properties of
the cost function
\begin{equation} \label{eq:cost}
  \mathscr C(\mathbf x)=\frac12\sum_{k=1}^M\big[V_k(\mathbf x)-V_0\big]^2
\end{equation}
which achieves zero only for configurations that satisfy all the constraints.
Here we dispense with the cost function and study the set of solutions
directly.
This set
can be written as
\begin{equation}
  \Omega=\big\{\mathbf x\in\mathbb R^N\mid \|\mathbf x\|^2=N,V_k(\mathbf x)=V_0
  \;\forall\;k=1,\ldots,M\big\}
\end{equation}
Because the constraints are all smooth functions, $\Omega$ is almost always a manifold without singular points.\footnote{The conditions for a singular point are that
  $0=\frac\partial{\partial\mathbf x}V_k(\mathbf x)$ for all $k$. This is
  equivalent to asking that the constraints $V_k$ all have a stationary point at
  the same place. When the $V_k$ are independent and random, this is vanishingly
  unlikely, requiring $NM+1$ independent equations to be simultaneously satisfied.
  This means that different connected components of the set of solutions do not
intersect, nor are there self-intersections, without extraordinary fine-tuning.}
We study the topology of the manifold $\Omega$ by two related means: its
average Euler characteristic, and the average number of stationary points of a
linear height function restricted to the manifold. These measures tell us
complementary pieces of information. We find that for the varied cases
we study, these two always coincide at the largest exponential order in $N$,
putting strong constraints on the resulting topology and geometry.

\section{The average Euler characteristic}

The Euler characteristic $\chi$ of a manifold is a topological invariant \cite{Hatcher_2002_Algebraic}. It is
perhaps most familiar in the context of connected compact orientable surfaces, where it
characterizes the number of handles in the surface: $\chi=2(1-\#)$ for $\#$
handles. For general $d$, the Euler characteristic of the $d$-sphere is $2$ if $d$ is even and 0 if $d$ is odd. The canonical method for computing the Euler characteristic is done by
defining a complex on the manifold in question, essentially a
higher-dimensional generalization of a polygonal tiling. Then $\chi$ is given
by an alternating sum over the number of cells of increasing dimension, which
for 2-manifolds corresponds to the number of vertices, minus the number of
edges, plus the number of faces.

Morse theory offers another way to compute the Euler characteristic of a manifold $\Omega$ using the
statistics of stationary points in a function $H:\Omega\to\mathbb R$ \cite{Audin_2014_Morse}. For
functions $H$ without any symmetries with respect to the manifold, the surfaces
of gradient flow between adjacent stationary points form a complex. The
alternating sum over cells to compute $\chi$ becomes an alternating sum over
the count of stationary points of $H$ with increasing index, or
\begin{equation}
  \chi(\Omega)=\sum_{i=0}^N(-1)^i\mathcal N_H(\text{index}=i)
\end{equation}
Conveniently, we can express this sum as an integral over the manifold
using a small variation on the Kac--Rice formula for counting stationary
points \cite{Kac_1943_On, Rice_1939_The}. Since the sign of the determinant of the Hessian matrix of $H$ at a
stationary point is equal to its index, if we count stationary points including
the sign of the determinant, we arrive at the Euler characteristic, or
\begin{equation} \label{eq:kac-rice}
  \chi(\Omega)=\int_\Omega d\mathbf x\,\delta\big(\nabla H(\mathbf x)\big)\det\operatorname{Hess}H(\mathbf x)
\end{equation}
When the Kac--Rice formula is used to \emph{count} stationary points, the sign
of the determinant is a nuisance that one must take pains to preserve
\cite{Fyodorov_2004_Complexity}. Here we are correct to exclude it.

We need to choose a function $H$ for our calculation. Because $\chi$ is
a topological invariant, any choice will work so long as it does not share some
symmetry with the underlying manifold, i.e., that $H$ satisfies the Smale condition. Because our manifold of random
constraints has no symmetries, we can take a simple height function $H(\mathbf
x)=\mathbf x_0\cdot\mathbf x$ for some $\mathbf x_0\in\mathbb R^N$ with
$\|\mathbf x_0\|^2=N$. $H$ is a height function because when $\mathbf x_0$ is
used as the polar axis, $H$ gives the height on the sphere.

We treat the integral over the implicitly defined manifold $\Omega$ using the
method of Lagrange multipliers. We introduce one multiplier $\omega_0$ to
enforce the spherical constraint and $M$ multipliers $\omega_k$ to enforce the vanishing of
each of the $V_k$, resulting in the Lagrangian
\begin{equation} \label{eq:lagrangian}
  L(\mathbf x,\pmb\omega)
  =H(\mathbf x)+\frac12\omega_0\big(\|\mathbf x\|^2-N\big)
  +\sum_{k=1}^M\omega_k\big(V_k(\mathbf x)-V_0\big)
\end{equation}
The integral over $\Omega$ in \eqref{eq:kac-rice} then becomes
\begin{equation} \label{eq:kac-rice.lagrange}
  \chi(\Omega)=\int_{\mathbb R^N} d\mathbf x\int_{\mathbb R^{M+1}}d\pmb\omega
  \,\delta\big(\partial L(\mathbf x,\pmb\omega)\big)
  \det\partial\partial L(\mathbf x,\pmb\omega)
\end{equation}
where $\partial=[\frac\partial{\partial\mathbf x},\frac\partial{\partial\pmb\omega}]$
is the vector of partial derivatives with respect to all $N+M+1$ variables.
This integral is now in a form where standard techniques from mean-field theory
can be applied to calculate it.

To evaluate the average of $\chi$ over the constraints, we first translate the $\delta$ functions and determinant to integral form, with
\begin{align}
  \label{eq:delta.exp}
  \delta\big(\partial L(\mathbf x,\pmb\omega)\big)
  &=\int\frac{d\hat{\mathbf x}}{(2\pi)^N}\frac{d\hat{\pmb\omega}}{(2\pi)^{M+1}}
  e^{i[\hat{\mathbf x},\hat{\pmb\omega}]\cdot\partial L(\mathbf x,\pmb\omega)}
  \\
  \label{eq:det.exp}
  \det\partial\partial L(\mathbf x,\pmb\omega)
  &=\int d\bar{\pmb\eta}\,d\pmb\eta\,d\bar{\pmb\gamma}\,d\pmb\gamma\,
  e^{-[\bar{\pmb\eta},\bar{\pmb\gamma}]^T\partial\partial L(\mathbf x,\pmb\omega)[\pmb\eta,\pmb\gamma]}
\end{align}
where $\hat{\mathbf x}$ and $\hat{\pmb\omega}$ are ordinary vectors and
$\bar{\pmb\eta}$, $\pmb\eta$, $\bar{\pmb\gamma}$, and $\pmb\gamma$ are
Grassmann vectors. With these expressions substituted into
\eqref{eq:kac-rice.lagrange}, the result is a integral over an exponential
whose argument is linear in the random functions $V_k$. These functions can
therefore be averaged over, and the resulting expression treated with standard
methods. Details of this calculation can be found in Appendix~\ref{sec:euler}.
The result is the reduction of the average Euler characteristic to an expression of the
form
\begin{equation} \label{eq:pre-saddle.characteristic}
  \overline{\chi(\Omega)}
  =\left(\frac N{2\pi}\right)^2\int dR\,dD\,dm\,d\hat m\,g(R,D,m,\hat m)\,e^{N\mathcal S_\chi(R,D,m,\hat m)}
\end{equation}
where $g$ is a prefactor of $o(N^0)$, and $\mathcal S_\chi$ is an effective action defined by
\begin{equation} \label{eq:euler.action}
  \begin{aligned}
    \mathcal S_\chi(R,D,m,\hat m\mid\alpha,V_0)
    &=\hat m-\frac\alpha2\left[
      \log\left(1+\frac{f(1)D}{f'(1)R^2}\right)
      +\frac{V_0^2}{f(1)}\left(1+\frac{f'(1)R^2}{f(1)D}\right)^{-1}
    \right] \\
    &\hspace{7em}+\frac12\log\left(
      1+\frac{(1-m^2)D+\hat m^2-2Rm\hat m}{R^2}
    \right)
  \end{aligned}
\end{equation}
The remaining order parameters defined by the scalar products
\begin{align}
  R=-i\frac1N\mathbf x\cdot\hat{\mathbf x}
  &&
  D=\frac1N\hat{\mathbf x}\cdot\hat{\mathbf x}
  &&
  m=\frac1N\mathbf x\cdot\mathbf x_0
  &&
  \hat m=-i\frac1N\hat {\mathbf x}\cdot\mathbf x_0
\end{align}

This integral can be evaluated to leading order by a saddle point method. For reasons we will
see, it is best to extremize with respect to $R$, $D$, and $\hat m$, leaving a
new effective action of $m$ alone. This can be solved to give
\begin{equation}
  D=-\frac{m+R^*}{1-m^2}R^* \qquad \hat m=0
\end{equation}
\begin{equation}
  \begin{aligned}
    R^*
    =\frac{-m(1-m^2)}{2[f(1)-(1-m^2)f'(1)]^2}
    \Bigg[
      \alpha V_0^2f'(1)
      +(2-\alpha)f(1)\left(\frac{f(1)}{1-m^2}-f'(1)\right) \quad \\
    \quad+\alpha\sqrt{
      \tfrac{4V_0^2}\alpha f(1)f'(1)\left[\tfrac{f(1)}{1-m^2}-f'(1)\right]
      +\left[\tfrac{f(1)^2}{1-m^2}-\big(V_0^2+f(1)\big)f'(1)\right]^2
    }
  \Bigg]
  \end{aligned}
\end{equation}
with the resulting effective action as a function of $m$ alone given by
\begin{equation}
  \mathcal S_\chi(m)
  =-\frac\alpha2\bigg[
    \log\left(
      1-\frac{f(1)}{f'(1)}\frac{1+\frac m{R^*}}{1-m^2}
    \right)
    +\frac{V_0^2}{f(1)}\left(
      1-\frac{f'(1)}{f(1)}\frac{1-m^2}{1+\frac m{R^*}}
    \right)^{-1}
  \bigg]
  +\frac12\log\left(-\frac m{R^*}\right)
\end{equation}
This function is plotted as a function of $m$ in Fig.~\ref{fig:action} for a variety of $V_0$ and $f$.
To finish evaluating the integral, this expression should be maximized with
respect to $m$. The order parameter $m$ is both physical and interpretable, as
it gives the overlap of the configuration $\mathbf x$ with the height axis
$\mathbf x_0$. Therefore, the value $m^*$ that maximizes this action can be
understood as the latitude on the sphere where most of the contribution to the
Euler characteristic is made.

\begin{figure}
  \includegraphics{figs/action_1.pdf}
  \hspace{-3.5em}
  \includegraphics{figs/action_3.pdf}

  \caption{
    \textbf{Effective action for the Euler characteristic.}
    The effective action governing the average Euler characteristic as a function of the overlap
    $m=\frac1N\mathbf x\cdot\mathbf x_0$ with the height direction for two
    different homogeneous polynomial functions and a variety of target values $V_0$. In both
    plots $\alpha=\frac12$. \textbf{Left:} With linear functions there are two
    regimes. For small $V_0$, there are maxima at $m=\pm m^*$ where the action
    is zero, while after the satisfiability transition at $V_0=V_\text{\textsc{sat}}=1$, $m^*$
    goes to zero and the action becomes negative. \textbf{Right:} With nonlinear
    functions there are four regimes. For small $V_0$ the behavior is the same
    as in the linear case, with zero action. After an onset transition at
    $V_0=V_\text{on}\simeq1.099$ the maxima are at the edge of validity of the
    action and the action is positive. At a shattering transition at
    $V_0=V_\text{sh}\simeq1.394$, the maximum goes to zero and the action is positive.
    Finally, at the satisfiability transition
    $V_0=V_\text{\textsc{sat}}\simeq1.440$ the action becomes negative.
  } \label{fig:action}
\end{figure}

The action $\mathcal S_\chi$ is extremized with respect to $m$ at $m=0$ or $m=\pm m^*=\mp R^*(m^*)$.
In the latter case, $m^*$ takes the value
\begin{equation}
  m^*=\sqrt{1-\frac{\alpha}{f'(1)}\big(V_0^2+f(1)\big)}
\end{equation}
and $\mathcal S_\chi(m^*)=0$. If this solution were always well-defined, it would vanish when the argument of the square root vanishes for
\begin{equation}
  V_0^2>V_{\text{\textsc{sat}}\ast}^2\equiv\frac{f'(1)}\alpha-f(1)
\end{equation}
This corresponds precisely to the satisfiability transition found in previous work by a replica symmetric analysis of the cost function \eqref{eq:cost} \cite{Fyodorov_2019_A, Fyodorov_2020_Counting, Fyodorov_2022_Optimization, Tublin_2022_A, Vivo_2024_Random}.
However, when the magnitude of $V_0$ is sufficiently large, with
\begin{equation}
  V_0^2>V_\text{on}^2\equiv\frac{1-\alpha+\sqrt{1-\alpha}}\alpha f(1)
\end{equation}
$R^*(m^*)$ becomes complex and this solution is no longer valid. Since
\begin{equation}
  V_\text{on}^2-V_{\text{\textsc{sat}}\ast}^2
  =\frac{f'(1)-f(1)}\alpha-\frac{\sqrt{1-\alpha}}\alpha f(1)
\end{equation}
If $f(q)$ is purely linear, then $f'(1)=f(1)$ and
$V_\text{on}^2>V_{\text{\textsc{sat}}\ast}^2$, so the naive satisfiability
transition happens first. On the other hand, when $f(q)$ contains powers of $q$
strictly greater than 1, then $f'(1)\geq 2f(1)$ and $V_\text{on}^2\leq
V_{\text{\textsc{sat}}\ast}^2$, so the onset happens first. In situations with
mixed constant, linear, and nonlinear terms in $f$, the order of the
transitions depends on the precise form of $f$.

Likewise, the solution at $m=0$ is not always valid. For $V_0$ of magnitude sufficiently small, with
\begin{equation}
  V_0^2<V_\text{sh}^2\equiv\frac{2(1+\sqrt{1-\alpha})-\alpha}{\alpha}f(1)\left(1-\frac{f(1)}{f'(1)}\right)
\end{equation}
the maximum at $m=0$ becomes complex and that solution is invalid. For purely
linear $f(q)$, $V_\text{sh}=0$ and the solution at $m=0$ is always real, though
for $V_0^2<V_{\text{\textsc{sat}}\ast}^2$ it is a minimum rather than a
maximum.
These transition values of the target $V_0$ correspond with transition values in $\alpha$ of
\begin{align}
  \alpha_\text{on}=1-\left(\frac{V_0^2}{V_0^2+f(1)}\right)^2
  &&
  \alpha_\text{sh}=4V_0^2f(1)f'(1)\frac{f'(1)-f(1)}{\big((V_0^2+f(1))f'(1)-f(1)^2\big)^2}
\end{align}
Finally, there is another satisfiability transition at
$V=V_\text{\textsc{sat}}$ corresponding to the vanishing of the effective
action at the $m=0$ solution. For generic covariance function $f$ it is not
possible to write an explicit formula for $V_\text{\textsc{sat}}$, and we
usually calculate it through a numeric root-finding algorithm. However, for
linear $f(q)$ one can see that
$V_\text{\textsc{sat}}=V_{\text{\textsc{sat}}\ast}$, and indeed this is the
case whenever $V_{\text{\textsc{sat}}\ast}^2<V_\text{on}^2$.

\begin{figure}
  \includegraphics{figs/phases_1.pdf}
  \hspace{-3em}
  \includegraphics{figs/phases_2.pdf}
  \hspace{-3em}
  \includegraphics{figs/phases_3.pdf}

  \caption{
    \textbf{Topological phase diagram.}
    Topological phases of the model for three different homogeneous covariance
    functions. The onset transition $V_\text{on}$, shattering transition
    $V_\text{sh}$, and satisfiability transition $V_\text{\textsc{sat}}$
    are indicated when they exist. In the limit of $\alpha\to0$, the behavior
    of level sets of the spherical spin glasses are recovered: the final plot
    shows how in the pure cubic model the ground state energy $E_\text{gs}$ and threshold energy
    $E_\text{th}$ correspond with the limits of the satisfiability and
    shattering transitions, respectively. Note that for mixed models with
    inhomogeneous covariance functions, $E_\text{th}$ is not the lower limit of
    $V_\text{sh}$.
  } \label{fig:phases}
\end{figure}

The phase diagram implied by these transitions is shown in
Fig.~\ref{fig:phases} for three different homogeneous $f(q)$. One interesting
feature occurs in the limit of $\alpha$ to zero. If $V_0$ is likewise rescaled
in the correct way, the limit of these phase boundaries approaches known
landmark energy values in the pure spherical spin glasses. In particular, the
limit to zero $\alpha$ of the scaled satisfiability transition
$V_\text{\textsc{sat}}\sqrt\alpha$ approaches the ground state energy
$E_\text{gs}$, while the limit to zero $\alpha$ of the scaled shattering
transition $V_\text{sh}\sqrt\alpha$ approaches the threshold energy
$E_\text{th}$. The correspondence between ground state and satisfiability is
expected: when the energy of a level set is greater in magnitude than the
ground state, the level set will usually be empty. The correspondence between
the threshold and shattering energies is also sensible, since the threshold
energy is typically understood as the point where the landscape fractures into
pieces. However, this second correspondence is only true for the pure spherical
models with homogeneous $f(q)$. For any other model with an inhomogeneous
$f(q)$, $E_\text{sh}^2<E_\text{th}^2$. This may have implications for dynamics
in such mixed models, see the discussion in Section~\ref{sec:ssg}.

Interpreting the specific topology implied by the Euler characteristic is difficult when the prediction is positive, but
a larger problem lies with the many negative values it can take. A manifold with many
handles will have a very negative Euler characteristic, and our calculation
relying on the saddle point of an exponential integral would be invalid. In the
regime of $\mathcal S_\chi(m)$ where the action is complex-valued, is this
breakdown the result of a negative average characteristic, or is it the result
of another reason? It is difficult to answer this with the previous calculation
alone. However, in the next section we will make a complementary calculation
that rules out the negative Euler characteristic picture. Instead, we will see
that in the complex region, the large-deviation function in $m$ that $\mathcal
S_\chi(m)$ represents breaks down due to the vanishingly small probability of
finding any stationary points.



\section{Complexity of a linear test function}

One way to definitely narrow possible interpretations of the average Euler
characteristic is to compute a complementary average. The Euler characteristic
is the alternating sum of numbers of critical points of different index. If
instead we make the direct sum
\begin{equation} \label{eq:abs.kac-rice}
  \mathcal N_H(\Omega)=\sum_{i=0}^N\mathcal N_H(\text{index}=i)
  =\int_\Omega d\mathbf x\,\delta\big(\nabla H(\mathbf x)\big)
  \,\big|\det\operatorname{Hess}H(\mathbf x)\big|
\end{equation}
we find the total number of stationary points. The formula is exactly the same
as that for the average Euler characteristic except for an absolute value sign
around the determinant of the Hessian.

Understanding the number of stationary points as a function of latitude $m$
will clarify the meaning of our effective action for the average Euler
characteristic in the range of overlaps $m$ where it takes a complex value. This is because the average number of stationary points is a
nonnegative number. If the region of complex $\mathcal S_\chi$ has a
well-defined number of stationary points, it indicates that we are looking at a
situation with a negative average Euler characteristic. On the other hand, if
the average number of stationary points yields a complex value at some latitude
$m$, it must be because it is either too large or small in $N$ to be captured by
the calculation, e.g., that it behaves like $e^{N^2\Sigma}$ or
$e^{-N^2\Sigma}$. The following calculation indicates this second situation:
the region of complex action is due to a lack of stationary points to
contribute to the Euler characteristic at those latitudes.

To compute the complexity, we follow a similar procedure to the Euler
characteristic. The main difference lies in how we treat the absolute value
function around the determinant. Following \cite{Fyodorov_2004_Complexity}, we
make use of the identity
\begin{equation} \label{eq:fyodorov.nightmare}
  \begin{aligned}
    |\det A|
  &=\lim_{\epsilon\to0}\frac{(\det A)^2}{\sqrt{\det(A+i\epsilon I)}\sqrt{\det(A-i\epsilon I)}}
  \\
  &=\frac1{(2\pi)^N}\int d\bar{\pmb\eta}_1\,d\pmb\eta_1\,d\bar{\pmb\eta}_2\,d\pmb\eta_2\,d\mathbf a_+\,d\mathbf a_2\,
  e^{
    -\bar{\pmb\eta}_1^TA\pmb\eta_1-\bar{\pmb\eta}_2^TA\pmb\eta_2
    -\frac12\mathbf a_+^T(A+i\epsilon I)\mathbf a_+-\frac12\mathbf a_2^T(A-i\epsilon I)\mathbf a_2
  }
\end{aligned}
\end{equation}
for an $N\times N$ matrix $A$. Here $\bar{\pmb\eta}_1$, $\pmb\eta_1$,
$\bar{\pmb\eta}_2$, and $\pmb\eta_2$ are Grassmann vectors and $\mathbf a$ and
$\mathbf b$ are regular vectors. This introduces many new order parameters into
the problem, but this is a difficulty of scale rather than principle. With this
identity substituted for the usual determinant one, the problem can be solved
much as before. The details of this solution are relegated to
Appendix~\ref{sec:complexity.details}. The result is that, to largest order in
$N$, the logarithm of the average number of stationary points is the same as
the logarithm of the average Euler characteristic.

\section{Implications for the topology of solutions}

It is not straightforward to directly use the average Euler characteristic to
infer something about the number of connected components in the set of
solutions. To understand why, a simple example is helpful. Consider the set of
solutions on the sphere $\|\mathbf x\|^2=N$ that satisfy the single quadratic
constraint
\begin{equation}
  0=\sum_{i=1}^N\sigma_ix_i^2
\end{equation}
where each $\sigma_i$ is taken to be $\pm1$ with equal probability. If we take $\mathbf x$ to be ordered such that all terms with $\sigma_i=+1$ come first, this gives
\begin{equation}
  0=\sum_{i=1}^{N_+}x_i^2-\sum_{i=N_++1}^Nx_i^2
\end{equation}
where $N_+$ is the number of terms with $\sigma_i=+1$. The topology of the resulting manifold can be found by adding and subtracting this constraint from the spherical one, which gives
\begin{align}
  \frac12=\sum_{i=1}^{N_+}x_i^2
  \qquad
  \frac12=\sum_{i=N_++1}^{N}x_i^2
\end{align}
These are two independent equations for spheres of radius $1/\sqrt2$, one of
dimension $N_+$ and the other of dimension $N-N_+$. Therefore, the topology of
the configuration space is that of $S^{N_+-1}\times S^{N-N_+-1}$. The Euler
characteristic of a product space is the product of the Euler characteristics,
and so we have $\chi(\Omega)=\chi(S^{N_+-1})\chi(S^{N-N_+-1})$.

What is the average value of the Euler characteristic over values of
$\sigma_i$? First, recall that the Euler characteristic of a sphere $S^d$ is 2
when $d$ is even and 0 when $d$ is odd. When $N$ is odd, any value
of $N_+$ will result in one of the two spheres in the product to be
odd-dimensional, and therefore $\chi(\Omega)=0$, as is always true for
odd-dimensional manifolds. When $N$ is even, there are two possibilities: when $N_+$ is even then both spheres are odd-dimensional, while when $N_+$ is odd then both spheres are even-dimensional.
The number of terms $N_+$ with $\sigma_i=+1$ is distributed with the binomial distribution
\begin{equation}
  P(N_+)=\frac1{2^N}\binom{N}{N_+}
\end{equation}
Therefore, the average Euler characteristic for even $N$ is
\begin{equation}
  \overline{\chi(\Omega)}
  =\sum_{N_+=0}^NP(N_+)\chi(S^{N_+-1})\chi(S^{N-N_+-1})
  =\frac4{2^N}\sum_{n=0}^{N/2}\binom{N}{2n}
  =2
\end{equation}
Thus we find the average Euler characteristic in this simple example is 2
despite the fact that the possible manifolds resulting from the constraints
have characteristics of either 0 or 4.

\begin{figure}
  \includegraphics{figs/regime_1.pdf}
  \hspace{-3.5em}
  \includegraphics{figs/regime_2.pdf}
  \hspace{-3.5em}
  \includegraphics{figs/regime_3.pdf}
  \hspace{-3.5em}
  \includegraphics{figs/regime_4.pdf}

  \caption{
    \textbf{Behavior of the action in four regimes.} The effective action $\mathcal S_\chi$ as a function of overlap $m$ with
    the height axis for our model with $f(q)=\frac12q^3$ and $\alpha=\frac12$
    at three different target values $V_0$. \textbf{Left: the connected
    regime.} The action is maximized with $\mathcal S_\chi(m^*)=0$, and no
    stationary points are found with overlap less than $m_\text{min}$ with a
    random point on the sphere. \textbf{Center: the onset regime.} The action
    is maximized with $\mathcal S_\chi(m_\text{min})>0$, and is positive up
    to $m_\text{max}$. No stationary points are found with overlap less than $m_\text{min}$. \textbf{Right: the shattered regime.} The action is
    maximized with $\mathcal S_\chi(0)>0$, and is positive up to
    $m_\text{max}$.
  }
\end{figure}

With these caveats in mind, we can combine the information from the previous
calculations to reason about what the topology and geometry of $\Omega$ should
be. We know what the average Euler characteristic is, and we know that it
corresponds to the average number of stationary points. We also know these two
averages correspond with each other when restricted to arbitrary latitudes $m$
associated with an arbitrary height axis $\mathbf x_0$, and we know their value
at each latitude. From this information, we draw the following inferences about
the form of the solution manifold in our four possible regimes.

\paragraph{The connected regime: \boldmath{$V_0^2<V_\text{on}^2$}.}

In our calculation above, $\overline{\chi(\Omega)}=2$ could mean a fine-tuned
average like this, or it could indicate the presence of manifold homeomorphic
to $S^{N-M-1}$ for even $N-M-1$. In either case it strongly indicates a single
connected component, whose few minima and maxima are usually found at the
latitude $m^*$. Randomly chosen points on the sphere have a typical nearest
overlap $m^*$ with the solution manifold, but can never have a smaller overlap than
$m_\text{min}$, indicating that the manifold is extensive.

\paragraph{The onset regime: \boldmath{$V_\text{on}^2<V_0^2<V_\text{sh}^2$}.}

In this regime $\log\overline{\chi(\Omega)}=O(N)$ but a minimum overlap
$m_\text{min}>0$ still exists. The minimum overlap indicates that the solution
manifold is exclusively made up of extensive components, because the existence
of small components would lead to stationary points near the equator with
respect to a randomly chosen axis $\mathbf x_0$. The solution manifold is
homeomorphic to a topological space with large Euler characteristic like the
product or disjoint union of many spheres. In the former case we would have one
topologically nontrivial connected component, while in the latter we would have
many simple disconnected components; the reality could be a combination of the
two. In the framework of this calculation, it is not possible to distinguish between
these scenarios. In any case, the minima and maxima of the height on the
solution manifold are typically found at latitude $m_\text{min}$ but are found
in exponential number up to the latitude $m_\text{max}$.

\paragraph{The shattered regime: \boldmath{$V_\text{sh}^2<V_0^2<V_\text{\textsc{sat}}^2$}.}

Here $\log\overline{\chi(\Omega)}=O(N)$ and the primary contribution to the
number of stationary points and to the Euler characteristic comes from the
equator. This indicates the presence of a large number of small disconnected
components to the solution manifold. While most minima and maxima of the height
are located at the equator $m=0$, they are found in exponential number up to
the latitude $m_\text{max}$.

\paragraph{The \textsc{unsat} regime: \boldmath{$V_\text{\textsc{sat}}^2<V_0^2$}.}

In this regime $\log\overline{\chi(\Omega)}<0$, indicating that in most
realizations of the functions $V_k$ the set $\Omega$ is empty.

\begin{figure}
  \includegraphics[width=0.245\textwidth]{figs/connected.pdf}
  \includegraphics[width=0.245\textwidth]{figs/coexist.pdf}
  \includegraphics[width=0.245\textwidth]{figs/shattered.pdf}
  \includegraphics[width=0.245\textwidth]{figs/gone.pdf}

  \includegraphics{figs/bar.pdf}

  \caption{
    \textbf{Cartoon of the solution manifold.}
    A low-dimensional sketch of the solution manifold in black, with stationary
    points of the height function as red points.
    The arrow shows the vector $\mathbf x_0$ defining the height
    function. For $V_0<V_\text{on}$, the manifold has a single connected
    component. Above the onset with $V_\text{on}<V_0<V_\text{sh}$, the manifold
    has large connected components with nontrivial topology. Above the shattering transition, or
    $V_\text{sh}<V_0<V_\text{\textsc{sat}}$, small disconnected pieces appear. Finally, above the satisfiability transition
    $V_\text{\textsc{sat}}$ the manifold vanishes.
  } \label{fig:cartoons}
\end{figure}

\paragraph{}
Fig.~\ref{fig:cartoons} shows a low-dimensional facsimile of what these
difference regimes look like. Surprisingly, this approach sees no signal of a
replica symmetry breaking (\textsc{rsb}) transition previously found in
\cite{Urbani_2023_A}. The
instability is predicted to occur when
\begin{equation} \label{eq:vrsb}
  V_0^2>V_\text{\textsc{rsb}}^2
  \equiv\frac{[f(1)-f(0)]^2}{\alpha f''(0)}
  -f(0)-\frac{f'(0)}{f''(0)}
\end{equation}
Some preliminary attempts to find an instability towards finite-\textsc{rsb} in
this region in the calculation of the Euler characteristic did not turn up anything. We conjecture that the \textsc{rsb}
instability found in \cite{Urbani_2023_A} is a trait of the cost function
\eqref{eq:cost}, and is not inherent to the structure of the solution manifold.
Perhaps the best evidence for this is to consider the limit of $M=1$, or
$\alpha\to0$ with $E=V_0\sqrt\alpha$ held fixed, where this problem reduces to
the level sets of the spherical spin glasses. The instability \eqref{eq:vrsb}
implies for the pure spherical 2-spin model with $f(q)=\frac12q^2$ that
$E_\textsc{rsb}=\frac12$, though nothing of note is known to occur in the level
sets of 2-spin model at such an energy.

\section{Implications for the dynamics of spherical spin glasses}
\label{sec:ssg}


As indicated earlier, for $M=1$ the solution manifold corresponds to the energy
level set of a spherical spin glass with energy density $E=\sqrt NV_0$. All the
results from the previous sections follow, and can be translated to the spin
glasses by taking the limit $\alpha\to0$ while keeping $E=V_0\alpha^{-1/2}$ fixed. With a little algebra this procedure yields
\begin{align}
  E_\text{on}=\pm\sqrt{2f(1)}
  &&
  E_\text{sh}=\pm\sqrt{4f(1)\left(1-\frac{f(1)}{f'(1)}\right)}
  \label{eq:ssg.energies}
\end{align}
for the energies at which level sets of the spherical spin glasses have
disconnected pieces appear, and that at which a large connected component
vanishes. For the pure $p$-spin spherical spin glasses with $f(q)=\frac12q^p$,
$E_\text{sh}=\sqrt{2(p-1)/p}$, precisely the threshold energy in these
models. This is expected, since threshold energy, defined as the place where
marginal minima are dominant in the landscape, is widely understood as the
place where level sets are broken into pieces.

However, for general mixed models the threshold energy is
\begin{equation}
  E_\mathrm{th}=\pm\frac{f(1)[f''(1)-f'(1)]+f'(1)^2}{f'(1)\sqrt{f''(1)}}
\end{equation}
which satisfies $|E_\text{sh}|\leq|E_\text{th}|$. Therefore, as one
descends in energy in generic models one will meet the shattering energy before
the threshold energy. This is perhaps unexpected, since one wight imagine that
where level sets of the energy break into many pieces would coincide with the
largest concentration of shallow minima in the landscape. We see here that this isn't the case.

This fact mirrors another another that was made clear
recently: when gradient decent dynamics are run on these models, they will
asymptotically reach an energy above the threshold energy
\cite{Folena_2020_Rethinking, Folena_2021_Gradient, Folena_2023_On}. The old
belief that the threshold energy qualitatively coincides with a kind of
shattering is the origin of the expectation that the dynamic limit should be
the threshold energy. Given that our work shows the actual shattering energy is
different, here we compare it with data on the asymptotic limits of dynamics to
see if they coincide.

\begin{figure}
  \includegraphics{figs/dynamics_2.pdf}
  \hspace{-0.5em}
  \includegraphics{figs/dynamics_3.pdf}

  \caption{
    \textbf{Is the shattering energy a dynamic threshold?} Comparison of the shattering energy $E_\text{sh}$ with the asymptotic
    performance of gradient descent from a random initial condition in $p+s$
    models with $p=2$ and $p=3$ and varying $s$. The values of $\lambda$ depend on $p$ and $s$ and are taken from \cite{Folena_2023_On}. The points show the asymptotic performance
    extrapolated using two different methods and have unknown uncertainty, from \cite{Folena_2023_On}. Also
    shown is the annealed threshold energy $E_\text{th}$, where marginal minima
    are the most common type of stationary point. The section of $E_\text{sh}$
    that is dashed on the left plot indicates the continuation of the annealed
    result, whereas the solid portion gives the value calculated with a
    {\oldstylenums 1}\textsc{frsb} ansatz.
  } \label{fig:ssg}
\end{figure}

Measurements of the asymptotic limits of dynamics were recently taken in
\cite{Folena_2023_On} for two different classes of models with inhomogeneous
$f(q)$, with
\begin{equation}
  f(q)=\frac12\big[\lambda q^p+(1-\lambda)q^s\big]
\end{equation}
They studied models with this covariance for $p=2$ and $p=3$ while varying $s$.
In both cases, the relative weight $\lambda$ varies with $s$ and was chosen to
maximize a heuristic to increase the chances of seeing nontrivial behavior. The
authors numerically integrated the dynamic mean field theory (\textsc{dmft})
equations for gradient descent on these models from a random initial condition to large but finite time, then attempted to
extrapolate the infinite-time behavior by two different methods. The black
symbols in Fig.~\ref{fig:ssg} show the measurements taken from
\cite{Folena_2023_On}. The difference between the two extrapolations is not
critical here, see the original paper for details. We simply note that the
authors of \cite{Folena_2023_On} did not associate an uncertainty with them,
nor were they confident that they are unbiased estimates of the asymptotic value.

Fig.~\ref{fig:ssg} also shows the shattering and {\oldstylenums1}\textsc{rsb}
threshold energies as a function of $s$. The solid lines come from using
\textsl{Mathematica}'s \texttt{Interpolation} function to create a smooth
function $\lambda(s)$ through the values used in \cite{Folena_2023_On}. For the
$2+s$ models for sufficiently large $s$, the ground state is described by a
{\oldstylenums1}\textsc{frsb} order and both the threshold energy and
shattering energy calculated using an annealed average are likely inaccurate
\cite{Auffinger_2022_The}. In Appendix~\ref{sec:1frsb} we calculate the
quenched ground state and shattering energies for these models consistent with
the {\oldstylenums1}\textsc{frsb} equilibrium order. In the left panel of Fig.~\ref{fig:ssg}, the solid line shows the quenched calculation, while the dashed line shows the annealed formula \eqref{eq:ssg.energies}.

Is the shattering energy consistent with the dynamic threshold for gradient
descent from a random initial condition? The evidence in Fig.~\ref{fig:ssg} is
compelling but inconclusive. The difference between the shattering energy and
the extrapolated \textsc{dmft} data is about the same as the difference between
the values predicted by the two extrapolation methods. If both extrapolation
methods suffer from similar systematic biases, it is plausible they true value
corresponds with the shattering energy. However, better estimates of the
asymptotic values are needed to support or refute this conjecture. This
motivates working to integrate the \textsc{dmft} equations to longer times, or
else look for analytic asymptotic solutions that approach $E_\text{sh}$.

\section{Conclusion}

We have shown how to calculate the average Euler characteristic of the solution
manifold in a simple model of random constraint satisfaction. The results
constrain the topology and geometry of this manifold, revealing when it is
connected and trivial, when it is extensive but topologically nontrivial, and
when it is shattered into disconnected pieces. These inferences are supported
by a auxiliary calculation of the complexity of stationary points for a test
function on the same solution manifold.

This calculation has novel implications for the geometry of the energy
landscape in the spherical spin glasses, where it reveals a previously unknown
landmark energy $E_\text{sh}$. This shattering energy is where the topological
calculation implies that the level set of the energy breaks into disconnected
pieces, and differs from the threshold energy $E_\text{th}$ in mixed models.
It's possible that $E_\text{sh}$ is the asymptotic limit of gradient descent
from a random initial condition in such models, but the quality of the
currently available data makes this conjecture inconclusive.

This paper has focused on equality constraints, while most existing studies of
constraint satisfaction study inequality constraints \cite{Franz_2016_The,
Franz_2017_Universality, Franz_2019_Critical, Annesi_2023_Star-shaped,
Baldassi_2023_Typical}. To generalize the technique developed in this paper to
such cases is not a trivial extension. The set of solutions to such problems
are manifolds with boundary, and these boundaries are often not smooth. To
study such cases with these techniques will require using extensions of the
Morse theory for manifolds with boundary, and will be the subject of future work.

\paragraph{Acknowledgements}
The authors thank Pierfrancesco Urbani for helpful conversations on these
topics, and Giampaolo Folena for supplying his \textsc{dmft} data for the
spherical spin glasses.

\paragraph{Funding information}
JK-D is supported by a \textsc{DynSysMath} Specific Initiative of the INFN.

\appendix

\section{Details of the calculation of the average Euler characteristic}
\label{sec:euler}

Our starting point is the expression \eqref{eq:kac-rice.lagrange} with the
substitutions of the $\delta$-function and determinant \eqref{eq:delta.exp} and
\eqref{eq:det.exp} made. To make the calculation compact, we introduce
superspace coordinates \cite{DeWitt_1992_Supermanifolds}. Introducing the Grassmann indices $\bar\theta_1$
and $\theta_1$, we define the supervectors
\begin{align}
  \pmb\phi(1)=\mathbf x+\bar\theta_1\pmb\eta+\bar{\pmb\eta}\theta_1+\bar\theta_1\theta_1i\hat{\mathbf x}
  &&
  \sigma_k(1)=\omega_k+\bar\theta_1\gamma_k+\bar\gamma_k\theta_1+\bar\theta_1\theta_1i\hat\omega_k
\end{align}
with associated measures
\begin{align}
  d\pmb\phi=d\mathbf x\,\frac{d\hat{\mathbf x}}{(2\pi)^N}\,d\bar{\pmb\eta}\,d\pmb\eta
  &&
  d\pmb\sigma=d\pmb\omega\,\frac{d\hat{\pmb\omega}}{(2\pi)^{M+1}}\,d\bar{\pmb\gamma}\,d\pmb\gamma
\end{align}
The Euler characteristic can be expressed using these supervectors as
\begin{align}
  &\chi(\Omega)
  =\int d\pmb\phi\,d\pmb\sigma\,e^{\int d1\,L(\pmb\phi(1),\pmb\sigma(1))} \\
  &=\int d\pmb\phi\,d\pmb\sigma\,\exp\left\{
    \int d1\left[
      H\big(\pmb\phi(1)\big)
      +\frac12\sigma_0(1)\left(\|\pmb\phi(1)\|^2-N\right)
      +\sum_{k=1}^M\sigma_k(1)\Big(V_k\big(\pmb\phi(1)\big)-V_0\Big)
    \right]
  \right\} \notag
\end{align}
where $d1=d\bar\theta_1\,d\theta_1$ is the integral over both Grassmann
indices. Since this is an exponential integrand linear in the Gaussian
functions $V_k$, we can take their average to find
\begin{equation} \label{eq:χ.post-average}
  \begin{aligned}
    \overline{\chi(\Omega)}
    =\int d\pmb\phi\,d\pmb\sigma\,\exp\Bigg\{
      \int d1\left[
        H(\pmb\phi(1))
        +\frac12\sigma_0(1)\big(\|\pmb\phi(1)\|^2-N\big)
        -V_0\sum_{k=1}^M\sigma_k(1)
      \right] \\
      +\frac12\int d1\,d2\,\sum_{k=1}^M\sigma_k(1)\sigma_k(2)f\left(\frac{\pmb\phi(1)\cdot\pmb\phi(2)}N\right)
    \Bigg\}
  \end{aligned}
\end{equation}
This is a super-Gaussian integral in the super-Lagrange multipliers $\sigma_k$ with $1\leq k\leq M$.
Performing that integral yields
\begin{equation}
  \begin{aligned}
    \overline{\chi(\Omega)}
    &=\int d\pmb\phi\,d\sigma_0\,\exp\Bigg\{
      \int d1\left[
        H(\pmb\phi(1))
        +\frac12\sigma_0(1)\big(\|\pmb\phi(1)\|^2-N\big)
      \right] \\
    &\hspace{5em}-\frac M2V_0^2\int d1\,d2\,f\left(\frac{\pmb\phi(1)\cdot\pmb\phi(2)}N\right)^{-1}
    -\frac M2\log\operatorname{sdet}f\left(\frac{\pmb\phi(1)\cdot\pmb\phi(2)}N\right)
  \Bigg\}
  \end{aligned}
\end{equation}
The supervector $\pmb\phi$ enters this expression as a function only of the
scalar product with itself and with the vector $\mathbf x_0$ inside the
height function $H(\mathbf x)=\mathbf x_0\cdot\mathbf x$. We therefore make a change
of variables to the superoperator $\mathbb Q$ and the supervector $\mathbb M$
defined by
\begin{equation}
  \mathbb Q(1,2)=\frac{\pmb\phi(1)\cdot\pmb\phi(2)}N
  \qquad
  \mathbb M(1)=\frac{\pmb\phi(1)\cdot\mathbf x_0}N
\end{equation}
These new variables can replace $\pmb\phi$ in the integral using a generalized Hubbard--Stratonovich transformation, which yields
\begin{align}
  &\overline{\chi(\Omega)}
    =\frac12\int d\mathbb Q\,d\mathbb M\,d\sigma_0\,
    \left([\operatorname{sdet}(\mathbb Q-\mathbb M\mathbb M^T)]^\frac12+O(N^{-1})\right)
    \,\exp\Bigg\{
    \frac N2\log\operatorname{sdet}(\mathbb Q-\mathbb M\mathbb M^T)
  \label{eq:post.hubbard-strat}
    \\
  &+N\int d1\left[
        \mathbb M(1)
        +\frac12\sigma_0(1)\big(\mathbb Q(1,1)-1\big)
      \right]
    -\frac M2V_0^2\int d1\,d2\,f(\mathbb Q)^{-1}(1,2)
    -\frac M2\log\operatorname{sdet}f(\mathbb Q)
  \Bigg\}
  \notag
\end{align}
where we show the asymptotic value of the prefactor in Appendix~\ref{sec:prefactor}.
To move on from this expression,
we need to expand the superspace notation. We can write
\begin{equation} \label{eq:ops.q}
  \begin{aligned}
    \mathbb Q(1,2)
    &=C-R(\bar\theta_1\theta_1+\bar\theta_2\theta_2)
    -G(\bar\theta_1\theta_2+\bar\theta_2\theta_1)
    -D\bar\theta_1\theta_1\bar\theta_2\theta_2 \\
    &\qquad
    +(\bar\theta_1+\bar\theta_2)H
    +\bar H(\theta_1+\theta_2)
    -(\bar\theta_1\theta_1\bar\theta_2+\bar\theta_2\theta_2\bar\theta_1)i\hat H
    -\bar{\hat H}(\theta_1\bar\theta_2\theta_2+\theta_1\bar\theta_1\theta_1)
  \end{aligned}
\end{equation}
and
\begin{equation}
  \mathbb M(1)
  =m+\bar\theta_1H_0+\bar H_0\theta_1+i\hat m\bar\theta_1\theta_1
  \label{eq:ops.m}
\end{equation}
with associated measures
\begin{align}
  d\mathbb Q
  =dC\,dR\,dG\,\frac{dD}{(2\pi)^2}\,d\bar H\,dH\,d\bar{\hat H}\,d\hat H
  &&
  d\mathbb M=dm\,\frac{d\hat m}{2\pi}\,d\bar H_0\,dH_0
  \label{eq:op.measures}
\end{align}
The order parameters $C$, $R$, $G$, $D$, $m$, and $\hat m$ are ordinary numbers defined by
\begin{align} \label{eq:ops.bos}
  C=\frac{\mathbf x\cdot\mathbf x}N
  &&
  R=-i\frac{\mathbf x\cdot\hat{\mathbf x}}N
  &&
  G=\frac{\bar{\pmb\eta}\cdot\pmb\eta}N
  &&
  D=\frac{\hat{\mathbf x}\cdot\hat{\mathbf x}}N
  &&
  m=\frac{\mathbf x_0\cdot\mathbf x}N
  &&
  \hat m=-i\frac{\mathbf x_0\cdot\hat{\mathbf x}}N
\end{align}
while $\bar H$, $H$, $\bar{\hat H}$, $\hat H$, $\bar H_0$ and $H_0$ are Grassmann numbers defined by
\begin{align} \label{eq:ops.ferm}
  \bar H=\frac{\bar{\pmb\eta}\cdot\mathbf x}N
  &&
  H=\frac{\pmb\eta\cdot\mathbf x}N
  &&
  \bar{\hat H}=-i\frac{\bar{\pmb\eta}\cdot\hat{\mathbf x}}N
  &&
  \hat H=-i\frac{\pmb\eta\cdot\hat{\mathbf x}}N
  &&
  \bar H_0=\frac{\bar{\pmb\eta}\cdot\mathbf x_0}N
  &&
  H_0=\frac{\pmb\eta\cdot\mathbf x_0}N
\end{align}
We can treat the integral over $\sigma_0$ immediately. It gives
\begin{equation} \label{eq:sigma0.integral}
  \int d\sigma_0\,e^{N\int d1\,\frac12\sigma_0(1)(\mathbb Q(1,1)-1)}
  =2\times2\pi\delta(C-1)\,\delta(G+R)\,\bar HH
\end{equation}
This therefore sets $C=1$ and $G=-R$ in the remainder of the integrand, as well
as removing all dependence on $\bar H$ and $H$. With these solutions inserted, the remaining terms in the exponential give
\begin{align} \label{eq:sdet.q}
  &\operatorname{sdet}(\mathbb Q-\mathbb M\mathbb M^T)
  =1+\frac{(1-m^2)D+\hat m^2-2Rm\hat m}{R^2}
  -\frac6{R^4}\bar H_0H_0\bar{\hat H}\hat H
  \\
  &\hspace{5pc}+\frac2{R^3}\left[
    (mR-\hat m)(\bar{\hat H}H_0+\bar H_0\hat H)
    -(D+R^2)\bar H_0H_0
    +(1-m^2)\bar{\hat H}\hat H
  \right] \notag
  \\ \label{eq:sdet.fq}
  &\operatorname{sdet}f(\mathbb Q)
  =1+\frac{Df(1)}{R^2f'(1)}
  +\frac{2f(1)}{R^3f'(1)}\bar{\hat H}\hat H
  \\ \label{eq:inv.q}
  &\int d1\,d2\,f(\mathbb Q)^{-1}(1,2)
  =\left(f(1)+\frac{R^2f'(1)}{D}\right)^{-1}
  +2\frac{Rf'(1)}{(Df(1)+R^2f'(1))^2}\bar{\hat H}\hat H
\end{align}
The Grassmann terms in these expressions do not contribute to the effective
action, but will be important in our derivation of the prefactor for the
calculation in the connected regime. The substitution of these expressions into \eqref{eq:post.hubbard-strat} without the Grassmann terms yields \eqref{eq:euler.action} from the main text.

\section{Calculation of the prefactor of the average Euler characteristic}
\label{sec:prefactor}


Because of our convention of including the appropriate factors of $2\pi$ in the
superspace measure, super-Gaussian integrals do not produce such factors in our
derivation. Prefactors to our calculation come from three sources: the
introduction of $\delta$-functions to define the order parameters, integrals
over Grassmann order parameters, and from the saddle point approximation to the
large-$N$ integral. In addition, there are important contributions of a sign of the magnetization at the solution that arise from our super-Gaussian integrations.

\subsection{Contribution from the Hubbard--Stratonovich transformations}
\label{sec:prefactor.hs}

First, we examine the factors arising from the definition of order parameters. This begins by introducing to the integral the factor of one
\begin{equation}
  1=(2\pi)^3\int d\mathbb Q\,d\mathbb M\,
  \delta\big(N\mathbb Q(1,2)-\pmb\phi(1)\cdot\pmb\phi(2)\big)\,
  \delta\big(N\mathbb M(1)-\mathbf x_0\cdot\pmb\phi(1)\big)
\end{equation}
where three factors of $2\pi$ come from the measures as defined in
\eqref{eq:op.measures}. Converting the $\delta$-function into an exponential
integral yields
\begin{equation}
  \begin{aligned}
    1=\frac12\int d\mathbb Q\,d\mathbb M\,d\tilde{\mathbb Q}\,d\tilde{\mathbb M}\,
    \exp\Bigg\{
      \frac12\int d1\,d2\,\tilde{\mathbb Q}(1,2)\big(N\mathbb Q(1,2)-\pmb\phi(1)\cdot\pmb\phi(2)\big) \qquad \\
      +\int d1\,i\tilde{\mathbb M}(1)\big(N\mathbb M(1)-\mathbf x_0\cdot\pmb\phi(1)\big)
    \Bigg\}
  \end{aligned}
\end{equation}
where the supervectors and measures for $\tilde{\mathbb Q}$ and $\tilde{\mathbb M}$ are defined
analogously to those of $\mathbb Q$ and $\mathbb M$. This is now a super-Gaussian integral in $\pmb\phi$, which can be performed to yield
\begin{equation}
  \begin{aligned}
  \int d\pmb\phi\,1=\frac12\int d\mathbb Q\,d\mathbb M\,d\tilde{\mathbb Q}\,d\tilde{\mathbb M}\,
  \exp\Bigg\{
    \frac N2\int d1\,d2\,\tilde{\mathbb Q}(1,2)\mathbb Q(1,2)
    +N\int d1\,i\tilde{\mathbb M}(1)\mathbb M(1) \quad \\
    -\frac N2\log\operatorname{sdet}\tilde{\mathbb Q}
    -\frac N2\int d1\,d2\,\tilde{\mathbb M}(1)\tilde{\mathbb Q}^{-1}(1,2)\tilde{\mathbb M}(2)
  \Bigg\}
  \end{aligned}
\end{equation}
We can perform the remaining Gaussian integral in $\tilde{\mathbb M}$ to find
\begin{equation}
  \begin{aligned}
    \int d\pmb\phi\,1=
    \frac12\int d\mathbb Q\,d\mathbb M\,d\tilde{\mathbb Q}\,
    (\operatorname{sdet}\tilde{\mathbb Q}^{-1})^{-\frac12}
    \exp\Bigg\{
      -\frac N2\log\operatorname{sdet}\tilde{\mathbb Q}
      \hspace{5em} \\
      \frac N2\int d1\,d2\,\tilde{\mathbb Q}(1,2)\big[
        \mathbb Q(1,2)-\mathbb M(1)\mathbb M(2)
      \big]
    \Bigg\}
  \end{aligned}
\end{equation}
The integral over $\tilde{\mathbb Q}$ can be evaluated to leading order using the saddle point
method. The integrand is stationary at $\tilde{\mathbb Q}=(\mathbb Q-\mathbb
M\mathbb M^T)^{-1}$, and substituting this into the above expression results in the term
$\frac12\log\det(\mathbb Q-\mathbb M\mathbb M^T)$ in the effective action from
\eqref{eq:post.hubbard-strat}. The saddle point also yields a prefactor of the form
\begin{equation}
  \left(\operatorname{sdet}_{\{1,2\},\{3,4\}}\frac{\partial^2\frac12\log\operatorname{sdet}\tilde{\mathbb Q}}{\partial\tilde{\mathbb Q}(1,2)\partial\tilde{\mathbb Q}(3,4)}\right)^{-\frac12}
  =\left(\operatorname{sdet}_{\{1,2\},\{3,4\}}\tilde{\mathbb Q}^{-1}(3,1)\tilde{\mathbb Q}^{-1}(2,4)\right)^{-\frac12}
  =1
\end{equation}
where the final superdeterminant is identically 1 for any superoperator $\tilde{\mathbb Q}$, not just its saddle-point value.
The Hubbard--Stratonovich transformation therefore contributes a factor of
\begin{equation}
  \frac12\operatorname{sdet}(\mathbb Q-\mathbb M\mathbb M^T)^{\frac12}
  =\frac12[(C-m^2)(D+\hat m^2)+(R-m\hat m)^2]^\frac12G^{-1}
\end{equation}
to the prefactor at the largest order in $N$.

\subsection{Sign of the prefactor}

The superspace notation papers over some analytic differences between branches
of the logarithm that are not important for determining the saddle point but
are important to getting correctly the sign of the prefactor. For instance, consider the superdeterminant of $\mathbb Q$ from \eqref{eq:ops.q} (dropping the fermionic order parameters for a moment for brevity),
\begin{equation}
  \operatorname{sdet}\mathbb Q=\frac{CD+R^2}{G^2}
\end{equation}
The numerator and denominator arise from the determinant in the sector of number and Grassmann number basis elements for the superoperator, respectively. In our calculation, such superdeterminants appear after Gaussian integrals, which produce
\begin{equation}
  \int d\pmb\phi\,\exp\left\{-\frac12\int d1\,d2\,\pmb\phi(1)\mathbb Q(1,2)\pmb\phi(2)\right\}
  =(\operatorname{sdet}\mathbb Q)^{-\frac12}
  =(CD+R^2)^{-\frac12}G
\end{equation}
Here we emphasize that in the expanded result of the integral, the term from
the denominator of the square root enters not as $(G^2)^{\frac12}=|G|$ but as
$G$, including its sign. Therefore, when we write in the effective action
$\frac12\log\operatorname{sdet}\mathbb Q$, we should really be writing
\begin{equation}
  \int d\pmb\phi\,\exp\left\{-\frac12\int d1\,d2\,\pmb\phi(1)\mathbb Q(1,2)\pmb\phi(2)\right\}
  =\operatorname{sign}(G)e^{-\frac12\log\operatorname{sdet}\mathbb Q}
\end{equation}
In our calculation in Appendix \ref{sec:euler} we elide this several times,
and accumulate $M$ factors of $\operatorname{sign}(-Gf'(C))=\operatorname{sign}(-G)$ from the Gaussian
integral over Lagrange multipliers and $N$ factors of $\operatorname{sign}(-G)$
from the Hubbard--Stratonovich transformation. Since at all saddle points $G=-R$, we have
\begin{equation}
  \operatorname{sign}(R)^{N+M}e^{N\mathcal S_\chi(\tilde{\mathbb Q},\mathbb Q,\mathbb M)}
\end{equation}

\subsection{Contribution from integrating the Grassmann order parameters}

After integrating out the Lagrange multiplier enforcing the spherical
constraint in \eqref{eq:sigma0.integral}, the Grassmann variables $\bar H$ and
$H$ are eliminated from the integrand. This leaves dependence on $\bar{\hat
H}$, $\hat H$, $\bar H_0$, and $H_0$. Expanding the contributions from
\eqref{eq:sdet.q}, \eqref{eq:sdet.fq}, and \eqref{eq:inv.q}, the contribution to the action is given by
\begin{equation}
  \int d\bar{\hat H}\,d\hat H\,d\bar H_0\,dH_0\,\exp\left\{
    N\begin{bmatrix}
      \bar{\hat H} & \bar H_0
    \end{bmatrix}
    \begin{bmatrix}
      h_1 & h_2 \\
      h_2 & h_3
    \end{bmatrix}
    \begin{bmatrix}
      \hat H \\
      H_0
    \end{bmatrix}
    +Nh_4\bar{\hat H}\hat H\bar H_0H_0
  \right\}
  =N^2(h_1h_3-h_2^2)+Nh_4
\end{equation}
where
\begin{align}
  &h_1=\frac1R\left(
    \frac{1-m^2}{D(1-m^2)+R^2-2Rm\hat m+\hat m^2}
    -\alpha\frac{Df(1)^2+R^2f'(1)[V_0^2+f(1)]}{[Df(1)+R^2f'(1)]^2}
  \right)
  \\
  &h_2=\frac1R\frac{Rm-\hat m}{D(1-m^2)+R^2-2Rm\hat m+\hat m^2}
  \qquad
  h_3=-\frac1R\frac{D+R^2}{D(1-m^2)+R^2-2Rm\hat m+\hat m^2}
  \\
  &h_4=-\frac1{R^2}\frac1{D(1-m^2)+R^2-2Rm\hat m+\hat m^2}
\end{align}
The contribution at leading order in $N$ is therefore
\begin{equation}
  \frac{N^2}{R^2[D(1-m^2)+R^2-2Rm\hat m+\hat m^2]}\left(
    \alpha\frac{(D+R^2)[Df(1)^2+R^2f'(1)[V_0^2+f(1)]]}{
      [Df(1)+R^2f'(1)]^2
    }
    -1
  \right)
\end{equation}

\subsection{Contribution from the saddle point approximation}

We now want to evaluate the prefactor for the asymptotic value of
$\overline{\chi(\Omega)}$. From the previous sections, the definition of the
measures $d\mathbb Q$ and $d\mathbb M$ in \eqref{eq:op.measures}, and the
integral over $\sigma_0$ of \eqref{eq:sigma0.integral}, we can now see that the
function $g(R,D,m,\hat m)$ of \eqref{eq:pre-saddle.characteristic} is given by
\begin{equation}
  \begin{aligned}
    g(R,D,m,\hat m)
    =-
    \frac{\operatorname{sign}(R)^{N+M}}{R^3[D(1-m^2)+R^2-2Rm\hat m+\hat m^2]^{\frac12}}
    \hspace{3pc} \\
    \times\left(
      \alpha\frac{(D+R^2)[Df(1)^2+R^2f'(1)[V_0^2+f(1)]]}{
        [Df(1)+R^2f'(1)]^2
      }
      -1
    \right)
  \end{aligned}
\end{equation}
In the connected regime, there are two saddle points of the integrand that
contribute to the asymptotic value of the integral, at $m=\pm m^*$ with $R=-m^*$, $D=0$, and $\hat m=0$. At this saddle point $\mathcal S_\chi=0$. We can therefore write
\begin{equation}
  \overline{\chi(\Omega)}
  =\sum_{m=\pm m^*}
    g(-m,0,m,0)\big[\det\partial\partial\mathcal S_\chi(-m,0,m,0)\big]^{-\frac12}
\end{equation}
where here $\partial=[\frac{\partial}{\partial R},\frac{\partial}{\partial D},\frac\partial{\partial m},\frac\partial{\partial \hat m}]$ is the
vector of derivatives with respect to the remaining order parameters. For both of the two saddle points, the determinant of the Hessian of the effective action evaluates to
\begin{equation}
  \det\partial\partial\mathcal S_\chi
  =\left[\frac1{(m^*)^4}\left(1-\frac{\alpha[V_0^2+f(1)]}{f'(1)}\right)\right]^2
\end{equation}
whereas
\begin{equation}
  g(\mp m^*,0,\pm m^*,0)
  =\frac{(\mp1)^{N+M+1}}{|m^*|^4}\left(1-\frac{\alpha[V_0^2+f(1)]}{f'(1)}\right)
\end{equation}
The saddle point at $m=-m^*$, characterized by minima of the height function, always contributes with a positive term. On the other hand, the saddle point with $m=+m^*$, characterized by maxima of the height function, contributes with a sign depending on if $N+M+1$ is even or odd. This follows from the fact that minima, with an index of 0, have a positive contribution to the sum over stationary points, while maxima, with an index of $N-M-1$, have a contribution that depends on the dimension of the manifold.

We have finally that, in the connected regime,
\begin{equation}
  \overline{\chi(\Omega)}=1+(-1)^{N+M+1}+O(N^{-1})
\end{equation}
When $N+M+1$ is odd, this evaluates to zero. In fact it must be zero to all
orders in $N$, since for odd-dimensional manifolds the Euler characteristic is
always zero. When $N+M+1$ is even, we have $\overline{\chi(\Omega)}=2$ to
leading order in $N$, as specified in the main text.

\section{Details for the average number of stationary points}
\label{sec:complexity.details}

Starting from \eqref{eq:abs.kac-rice}, we make the substitution
\eqref{eq:delta.exp} to treat the Dirac $\delta$-function and use
\eqref{eq:fyodorov.nightmare} to write the absolute value of the determinant as
\begin{align}
  &\big|\det\partial\partial L(\mathbf x,\pmb\omega)\big|
  =\lim_{\epsilon\to0}\frac1{(2\pi)^{N+M+1}}\int d\bar{\pmb\eta}_1\,d\pmb\eta_1\,d\bar{\pmb\eta}_2\,d\pmb\eta_2\,
  d\bar{\pmb\gamma}_1\,d\pmb\gamma_1\,d\bar{\pmb\gamma}_2\,d\pmb\gamma_2\,
  d\mathbf a_1\,d\mathbf a_2\,d\mathbf b_1\,d\mathbf b_2\, \notag \\
  &\qquad\times\exp\Bigg\{
    -\frac12\begin{bmatrix}\mathbf a_1^T&\mathbf b_1^T\end{bmatrix}\big(\partial\partial L(\mathbf x,\pmb\omega)+i\epsilon I\big)\begin{bmatrix}\mathbf a_1\\\mathbf b_1\end{bmatrix}
    -\frac12\begin{bmatrix}\mathbf a_2^T&\mathbf b_2^T\end{bmatrix}\big(\partial\partial L(\mathbf x,\pmb\omega)-i\epsilon I\big)\begin{bmatrix}\mathbf a_2\\\mathbf b_2\end{bmatrix}
    \notag \\
  &\hspace{6em}-\begin{bmatrix}\bar{\pmb\eta}_1^T&\bar{\pmb\gamma}_1^T\end{bmatrix}\partial\partial L(\mathbf x,\pmb\omega)\begin{bmatrix}\pmb\eta_1\\\pmb\gamma_1\end{bmatrix}
    -\begin{bmatrix}\bar{\pmb\eta}_2^T&\bar{\pmb\gamma}_2^T\end{bmatrix}\partial\partial L(\mathbf x,\pmb\omega)\begin{bmatrix}\pmb\eta_2\\\pmb\gamma_2\end{bmatrix}
  \Bigg\}
  \label{eq:abs.det.exp}
\end{align}
where the $\mathbf a_{\sfrac12}$ are in $\mathbb R^N$, the $\mathbf b_{\sfrac12}$ are in
$\mathbb R^M$, and the $\pmb\eta_{\sfrac12}$ and $\pmb\gamma_{\sfrac12}$ are Grassmann vectors just as in
\eqref{eq:det.exp}, except with an extra copy of each. This zoo of vectors
quickly becomes tiring. Thankfully, there is a way to compactly represent this
calculation again using superspace vectors.

Consider the vectors
\begin{align}
  \pmb\phi(1,2)
  &=\mathbf x
  +\bar\theta_1\pmb\eta_1+\bar{\pmb\eta}_1\theta_1\bar\theta_2\theta_2
  +\bar\theta_2\pmb\eta_2+\bar{\pmb\eta}_2\theta_2\bar\theta_1\theta_1 \\
  &\qquad+\frac1{\sqrt2}(\bar\theta_1\theta_2+\bar\theta_2\theta_1)\mathbf a_1
  +\frac i{\sqrt2}(\bar\theta_1\theta_1+\bar\theta_2\theta_2)\mathbf a_2
  +\bar\theta_1\theta_1\bar\theta_2\theta_2i\hat{\mathbf x}
  \notag \\
  \sigma_k(1,2)
  &=\omega_k
  +\bar\theta_1\gamma_{1k}+\bar{\gamma}_{1k}\theta_1\bar\theta_2\theta_2
  +\bar\theta_2\gamma_{2k}+\bar{\gamma}_{2k}\theta_2\bar\theta_1\theta_1 \\
  &\qquad+\frac1{\sqrt2}(\bar\theta_1\theta_2+\bar\theta_2\theta_1)b_{1k}
  +\frac i{\sqrt2}(\bar\theta_1\theta_1+\bar\theta_2\theta_2)b_{2k}
  +\bar\theta_1\theta_1\bar\theta_2\theta_2\hat\omega_k
  \notag
\end{align}
The entire expression for the complexity of stationary points combining
\eqref{eq:delta.exp} and \eqref{eq:abs.det.exp} can be expressed in the compact
form
\begin{equation}
  \mathcal N_H(\Omega)
  =\lim_{\epsilon\to0}\int d\pmb\phi\,d\pmb\sigma\,e^{
   \int d1\,d2\,L(\pmb\phi(1,2),\pmb\sigma(1,2))
   -\frac{i\epsilon}2
   (\|\mathbf a_1\|^2-\|\mathbf a_2\|^2+\|\mathbf b_1\|^2-\|\mathbf b_2\|^2)
 }
\end{equation}
Note that unlike in the derivation of the average Euler characteristic,
$\pmb\phi(1,2)$ does not span the entire superspace. Therefore, we will not be
able to write expressions like the superdeterminant of
$f(\pmb\phi^T\pmb\phi/N)$, which are formally zero.

Following similar steps as for the Euler characteristic, we first take the average over the random constraint functions $V$. This yields
\begin{equation}
  \begin{aligned}
    \overline{\mathcal N_H(\Omega)}
    =\int d\pmb\phi\,d\pmb\sigma\,\exp\Bigg\{
      \frac12\int d1\,d2\,d3\,d4\,\sum_{k=1}^M\sigma_k(1,2)\sigma_k(3,4)f\left(\frac{\pmb\phi(1,2)\cdot\pmb\phi(3,4)}N\right)
      \\
      +\int d1\,d2\left[
        H(\pmb\phi(1,2))
        +\frac12\sigma_0(1,2)\big(\|\pmb\phi(1,2)\|^2-N\big)
        -V_0\sum_{k=1}^M\sigma_k(1,2)
      \right]
    \Bigg\}
  \end{aligned}
\end{equation}
We once again have a Gaussian integral in the parameters inside the Lagrange
multipliers $\sigma_k$ for $k=1,\ldots,M$. However, here we must expand the
superfields at this stage. Before that, we introduce order parameters analogous
to before. Alongside those in \eqref{eq:ops.bos}, we have
\begin{align}
  A_{ij}=\frac{\mathbf a_i\cdot\mathbf a_j}N
  &&
  X_i=\frac{\mathbf x\cdot\mathbf a_i}N
  &&
  \hat X_i=-i\frac{\hat{\mathbf x}\cdot\mathbf a_i}N
  &&
  H_{ij}=\frac{\bar{\pmb\eta}_i\cdot\pmb\eta_j}N
  &&
  J=\frac{\pmb\eta_1\cdot\pmb\eta_2}N
  &&
  m_i = \frac{\mathbf x_0\cdot\mathbf a_i}N
\end{align}
where the $i$ and $j$ run over $1$ and $2$.
These come with a volume term in the action
\begin{equation}
  \begin{aligned}
    \mathcal V=\frac12\log\det\left(
      \begin{bmatrix}
        C   & iR        & X_1       & X_2       \\
        iR  & D         & i\hat X_1 & i\hat X_2 \\
        X_1 & i\hat X_1 & A_{11}    & A_{12}    \\
        X_2 & i\hat X_2 & A_{12}    & A_{22}
      \end{bmatrix}
      -\begin{bmatrix}
        m \\
        i\hat m \\
        m_1 \\
        m_2
      \end{bmatrix}
      \begin{bmatrix}
        m &
        i\hat m &
        m_1 &
        m_2
      \end{bmatrix}
    \right) \qquad \\
    -\frac12\log\det\begin{bmatrix}
      0 & G_{11} & \bar J & G_{12} \\
      -G_{11} & 0 & -G_{21} & J \\
      -\bar J & G_{21} & 0 & G_{22} \\
      -G_{12} & -J & -G_{22} & 0
    \end{bmatrix}
  \end{aligned}
\end{equation}
As before, we can immediately integrate out $\sigma_0$, which fixes certain of the order parameters. In particular, we find
\begin{align}
  C=1 &&
  X_1=0 &&
  X_2=0 &&
  G_+=-R-\frac12A_+
\end{align}
where we have defined the symmetric and antisymmetric combinations of $G_{11}$ and $G_{22}$
\begin{align}
  G_+=G_{11}+G_{22}
  &&
  G_-=G_{11}-G_{22}
  &&
  A_+=A_{11}+A_{22}
  &&
  A_-=A_{11}-A_{22}
\end{align}
Once the first three substitutions have been made, the result for the Gaussian
integral in the ordinary variables $\omega_k$, $\hat\omega_k$, $b_{1k}$, and
$b_{2k}$ is a kernel with
\begin{equation}
  K=\begin{bmatrix}
    K_{11} & iRf'(1) & -\hat X_1 f'(1) & -\hat X_2 f'(1) \\
    iRf'(1) & f(1) & 0 & 0 \\
    -\hat X_1 f'(1) & 0 & i\epsilon-\tfrac12f'(1)(A_++A_-) & -A_{12}f'(1) \\
    -\hat X_2 f'(1) & 0 & -A_{12}f'(1) & -i\epsilon-\tfrac12f'(1)(A_+-A_-)
  \end{bmatrix}
\end{equation}
\begin{equation}
  K_{11}=Df'(1)-\frac12f''(1)\left[
      (R-\tfrac12A_+)^2+\tfrac12A_-^2+2A_{12}^2-(G_-^2+4G_{12}G_{21}+4\bar JJ)
    \right]
\end{equation}
and likewise a Gaussian integral in the Grassmann variables with kernel
\begin{equation}
  K'=f'(1)\begin{bmatrix}
    0 & G_{11} & J & G_{21} \\
    -G_{11} & 0 & -G_{12} & \bar J \\
    -J & G_{12} & 0 & G_{22} \\
    -G_{21} & -\bar J & -G_{22} & 0
  \end{bmatrix}
\end{equation}
The effective action then has the form
\begin{equation}
  \mathcal S
  =\hat m+\frac12i\epsilon A_--\frac\alpha2\left(
    \log\det K-\log\det K'+V_0^2(K^{-1})_{22}
  \right)+\mathcal V
\end{equation}
It is not helpful to write out this entire expression, which is quite large.
However, we only find saddle points of this action with $\bar
J=J=G_{12}=G_{21}=G_-=\hat X_1=\hat X_2=m_1=m_2=A_{12}=0$. Setting these variables to zero, we find
\begin{align}
  \notag
  \mathcal S
  =\hat m
  +\frac12i\epsilon A_-
  -\alpha\frac12\log\frac{
    f'(1)(Df(1)+R^2f'(1))-\frac12((R-\frac12A_+)^2+\frac12 A_-^2)f(1)f''(1)
  }{
    \frac14(R+\frac12A_+)^2f'(1)^2
  }
  \\
  \notag
  -\alpha\frac12\log\frac{
    \epsilon^2+i\epsilon f'(1)A_-+\frac14(A_+^2-A_-^2)f'(1)^2
  }{
    \frac14(R+\frac12A_+)^2f'(1)^2
  }
  +\frac12\log\frac{D(1-m^2)+R^2-2Rm\hat m+\hat m^2}{\frac14(R+\frac12A_+)^2}
  \\
  +\frac12\log\frac{A_+^2-A_-^2}{(R+\frac12A_+)^2}
  -\frac12\alpha V_0^2\left(
    f(1)+\frac{R^2f'(1)^2}{Df'(1)-\frac12((R-\frac12A_+)^2+\frac12A_-^2)f''(1)}
  \right)^{-1}
\end{align}
One solution to these equations at $\epsilon=0$ is $A_-=0$ and $A_+=2R^*$ with
$D=\hat m=0$ and $R=R^*$, exactly as for the Euler characteristic. The resulting effective
action as a function of $m$ is also exactly the same. We find two other
solutions that differ from this one, but with formulae too complex to share in
the text. One alternate solution has $A_-=0$ and one has $A_-\neq0$.

\begin{figure}
  \includegraphics{figs/alt_sols_1.pdf}
  \hfill
  \includegraphics{figs/alt_sols_2.pdf}

  \caption{
    \textbf{Choosing the correct saddle-point for the complexity.} In the
    calculation of the complexity of stationary points in the height function,
    there are three plausible solutions in some regimes. \textbf{Left:} Plot of
    the three solutions for $M=1$, $f(q)=\frac12q^3$, and
    $E=E_\text{on}=3^{-1/2}$. The inset shows detail of the region where these
    solutions are nonnegative. The solid and dashed lines show the values of
    $m$ where the $A_-=0$ and $A_-\neq0$ alternative solutions are maximized,
    respectively. The value $m^*$ maximizing the action corresponding also to
    $\mathcal S_\chi$ is also marked. \textbf{Right:} Numeric measurement of the
    average value of $m$ for stationary points found using Newton's method as a
    function of model size $N$. The thick solid line is a power-law fit to
    $m-m^*$, while the thin solid and dashed lines show the same values of $m$
    as in the lefthand plot.
  } \label{fig:alt.sols}
\end{figure}

For both alternative solutions, there is a regime of small $V_0$ or $E$ for
which they are unphysical. In this case the action is everywhere negative
despite our knowledge from the calculation of the Euler characteristic that the
manifold exists. However, there is a second regime where they cross the action
$\mathcal S_\chi$ and become positive, both reaching values of $m$ that are
less than $m_\text{min}$. An example of this regime for the case with
$f(q)=\frac12q^3$ and $M=1$ (the spherical spin glass) is shown in the left
panel of Fig.~\ref{fig:alt.sols}. In order to investigate which solution is
valid in this regime, we turn to comparison with numeric experiments. We create
samples of this model at the specified parameters, choose random axis $\mathbf
x_0$, and use Newton's method with mild damping to find stationary points of
the Lagrangian starting from random initial conditions. When we find a
stationary point, we record the value of the overlap $m=\frac1N\mathbf
x\cdot\mathbf x_0$ between it and the height axis. The right panel of
Fig.~\ref{fig:alt.sols} shows that the average value of $m$ attained by this
method as a function of $N$ asymptotically approaches $m^*$, the value implied
by the Euler characteristic action $\mathcal S_\chi$, and not either of the
values implied by the alternative solutions. This indicates that $\mathcal
S_\mathcal N(m)=\mathcal S_\chi(m)$ in all regimes.

\section{The quenched shattering energy}
\label{sec:1frsb}

Here we share how the quenched shattering energy is calculated under a
{\oldstylenums1}\textsc{frsb} ansatz. To best make contact with prior work on
the spherical spin glasses, we start with \eqref{eq:χ.post-average}. The
formula in a quenched calculation is almost the same as that for the annealed,
but the order parameters $C$, $R$, $D$, and $G$ must be understood as $n\times
n$ matrices rather than scalars. In principle $m$, $\hat m$, $\omega_0$, $\hat\omega_0$, $\omega_1$, and $\hat\omega_1$ should be considered $n$-dimensional vectors, but since in our ansatz replica vectors are constant we can take them to be constant from the start. We have
\begin{align}
  &\overline{\log\chi(\Omega)}
  =\lim_{n\to0}\frac\partial{\partial n}\int dC\,dR\,dD\,dG\,dm\,d\hat m\,d\omega_0\,d\hat\omega_0\,d\omega_1\,d\hat\omega_1\,
  \exp N\Bigg\{
    n\hat m
    +\frac i2\hat\omega_0\operatorname{Tr}(C-I)
    \notag \\
  &\hspace{2em}-\omega_0\operatorname{Tr}(G+R)
    -in\hat\omega_1E_0
    +\frac12\log\det\begin{bmatrix}
      C-m^2 & i(R-m\hat m) \\ i(R-m\hat m) & D-\hat m^2
    \end{bmatrix}
    -\frac12\log G^2
    \notag \\
    &\hspace{2em}-\frac12\sum_{ab}^n\left[
      \hat\omega_1^2f(C_{ab})
      +(2i\omega_1\hat\omega_1R_{ab}+\omega_1^2D_{ab})f'(C_{ab})
      +\omega_1^2(G_{ab}^2-R_{ab}^2)f''(C_{ab})
    \right]
  \Bigg\}
\end{align}
which is completely general for the spherical spin glasses with $M=1$. We now
make a series of simplifications. Ward identities associated with the BRST
symmetry possessed by the original action \cite{Annibale_2003_The, Annibale_2003_Supersymmetric, Annibale_2004_Coexistence} indicate that
\begin{align}
  \omega_1D=-i\hat\omega_1R
  &&
  G=-R
  &&
  \hat m=0
\end{align}
Moreover, this problem with $m=0$ has a close resemblance to the complexity of
the spherical spin glasses. In both, at the supersymmetric saddle point the
matrix $R$ is diagonal with $R=r_dI$ \cite{Kent-Dobias_2023_How}.
To investigate the shattering energy, we can restrict to solutions with $m=0$
and look for the place where such solutions vanish. Inserting these simplifications, we have up to highest order in $N$
\begin{equation}
  \begin{aligned}
    \overline{\log\chi(\Omega)}
    =\lim_{n\to0}\frac\partial{\partial n}\int dC\,dR\,d\hat\omega_0\,d\omega_1\,d\hat\omega_1\,
    \exp N\Bigg\{
      \frac i2\hat\omega_0\operatorname{Tr}(C-I)
      -in\hat\omega_1E
      \qquad\\
      -i\frac12n\omega_1\hat\omega_1r_df'(1)
      -\frac12\sum_{ab}^n
        \hat\omega_1^2f(C_{ab})
      +\frac12\log\det
      \left(\frac{-i\hat\omega_1}{\omega_1r_d}C+I\right)
    \Bigg\}
  \end{aligned}
\end{equation}
If we redefine $\hat\beta=-i\hat\omega_1$ and $\tilde r_d=\omega_1 r_d$, we find
\begin{equation}
  \begin{aligned}
    \overline{\log\chi(\Omega)}
    =\lim_{n\to0}\frac\partial{\partial n}\int dC\,d\hat\beta\,d\tilde r_d\,\hat\omega_0\,
    \exp N\Bigg\{
      \frac i2\hat\omega_0\operatorname{Tr}(C-I)
      +n\hat\beta E
      \qquad\\
      +n\frac12\hat\beta\tilde r_df'(1)
      +\frac12\sum_{ab}^n
        \hat\beta^2f(C_{ab})
      +\frac12\log\det
      \left(\frac{\hat\beta}{\tilde r_d}C+I\right)
    \Bigg\}
  \end{aligned}
\end{equation}
which is exactly the effective action for the supersymmetric complexity in the
spherical spin glasses when in the regime where minima dominate
\cite{Kent-Dobias_2023_How}. As the effective action for the Euler characteristic, this expression is valid whether minima dominate or not. Following the same steps as in
\cite{Kent-Dobias_2023_How}, we can write the continuum version of this action
for arbitrary \textsc{rsb} structure in the matrix $C$ as
\begin{equation} \label{eq:cont.action}
  \frac1N\overline{\log\chi(\Omega)}=\hat\beta E+\frac12\hat\beta\tilde r_df'(1)
  +\frac12\int_0^1dq\,\left[
    \hat\beta^2f''(q)\chi(q)+\frac1{\chi(q)+\tilde r_d\hat\beta^{-1}}
  \right]
\end{equation}
where $\chi(q)=\int_1^qdq'\int_0^{q'}dq''P(q'')$ and $P(q)$ is the distribution of
off-diagonal elements of the matrix $C$. This action must be extremized over
the function $\chi$ and the variables $\hat\beta$ and $\tilde r_d$, under the
constraint that $\chi(q)$ is continuous, that it has $\chi'(1)=-1$, and
$\chi(1)=0$, necessary for $P$ to be a well-defined probability distribution.

Now the specific form of replica symmetry breaking we expect to see is
important. We want to study the mixed $2+s$ models in the regime where they may
have 1-full \textsc{rsb} in equilibrium \cite{Auffinger_2022_The}. For the Euler characteristic like the
complexity, this will correspond to full \textsc{rsb}, in an analogous way to
{\oldstylenums1}\textsc{rsb} equilibria give a \textsc{rs} complexity. Such order is characterized by a piecewise smooth $\chi$ of the form
\begin{equation}
  \chi(q)=\begin{cases}
    \chi_0(q) & q < q_0 \\
    1-q & q \geq q_0
  \end{cases}
\end{equation}
where $\chi_0$ is
\begin{equation}
  \chi_0(q)=\frac1{\hat\beta}[f''(q)^{-1/2}-\tilde r_d]
\end{equation}
the function implied by extremizing \eqref{eq:cont.action} over $\chi$. The
variable $q_0$ must be chosen so that $\chi$ is continuous. The key difference
between \textsc{frsb} and {\oldstylenums1}\textsc{frsb} in this setting is that
in the former case the ground state has $q_0=1$, while in the latter the ground
state has $q_0<1$.

We use this action to find the shattering energy in the following way. First,
we know that the ground state energy is the place where the manifold and therefore the average Euler characteristic vanishes. Therefore, setting $\overline{\log\chi(\Omega)}=0$ and solving for $E$ yields a formula
for the ground state energy
\begin{equation}
  E_\text{gs}=-\frac1{\hat\beta}\left\{
    \frac12\hat\beta\tilde r_df'(1)
    +\frac12\int_0^1dq\,\left[
      \hat\beta^2f''(q)\chi(q)+\frac1{\chi(q)+\tilde r_d\hat\beta^{-1}}
    \right]
  \right\}
\end{equation}
This expression can be maximized over $\hat\beta$ and $\tilde r_d$ to find the
correct parameters at the ground state for a particular model. Then, the
shattering energy is found by slowly lowering $q_0$ and solving the combined
extremal and continuity problem for $\hat\beta$, $\tilde r_d$, and $E$ until
$E$ reaches a maximum value and starts to decrease. This maximum is the
shattering energy, since it is the point where the $m=0$ solution vanishes.
Starting from this point, we take small steps in $s$ and $\lambda_s$, again
simultaneously extremizing, ensuring continuity, and maximizing $E$. This draws
out the shattering energy across the entire range of $s$ plotted in
Fig.~\ref{fig:ssg}. The transition to the \textsc{rs} solution occurs when the value $q_0$ that maximizes $E$ hits zero.

\bibliography{topology}

\end{document}