\documentclass[aps,prl,reprint,longbibliography,floatfix]{revtex4-2}

\usepackage[utf8]{inputenc} % why not type "Bézout" with unicode?
\usepackage[T1]{fontenc} % vector fonts plz
\usepackage{newtxtext,newtxmath} % Times for PR
\usepackage[
  colorlinks=true,
  urlcolor=purple,
  citecolor=purple,
  filecolor=purple,
  linkcolor=purple
]{hyperref} % ref and cite links with pretty colors
\usepackage{amsmath, graphicx, xcolor} % standard packages

\begin{document}

\title{Complex complex landscapes}

\author{Jaron Kent-Dobias}
\author{Jorge Kurchan}

\affiliation{Laboratoire de Physique de l'Ecole Normale Supérieure, Paris, France}

\date\today

\begin{abstract}
  We study the saddle-points of the $p$-spin model -- the best understood
  example of a `complex' (rugged) landscape -- when its $N$ variables are
  complex. These points are the solutions to a system of $N$ random equations
  of degree $p-1$.  We solve for $\overline{\mathcal{N}}$, the number of
  solutions averaged over randomness in the $N\to\infty$ limit.  We find that
  it saturates the Bézout bound $\log\overline{\mathcal{N}}\sim N \log(p-1)$.
  The Hessian of each saddle is given by a random matrix of the form $C^\dagger
  C$, where $C$ is a complex symmetric Gaussian matrix with a shift to its diagonal.  Its
  spectrum has a transition where a gap develops that generalizes the notion of
  `threshold level' well-known in the real problem.  The results from the real
  problem are recovered in the limit of real parameters. In this case, only the
  square-root of the total number of solutions are real.  In terms of the
  complex energy, the solutions are divided into sectors where the saddles have
  different topological properties.
\end{abstract}

\maketitle

Spin-glasses are the paradigm of many-variable `complex landscapes,' a category
that also includes neural networks and optimization problems like constraint
satisfaction \cite{Mezard_2009_Information}.  The most tractable family of
these are the mean-field spherical $p$-spin models \cite{Crisanti_1992_The}
(for a review see \cite{Castellani_2005_Spin-glass}) defined by the energy
\begin{equation} \label{eq:bare.hamiltonian}
  H_0 = \frac1{p!}\sum_{i_1\cdots i_p}^NJ_{i_1\cdots i_p}z_{i_1}\cdots z_{i_p},
\end{equation}
where $J$ is a symmetric tensor whose elements are real Gaussian variables and
$z\in\mathbb R^N$ is constrained to the sphere $z^Tz=N$. This problem has been
studied in the algebra \cite{Cartwright_2013_The} and probability literature
\cite{Auffinger_2012_Random, Auffinger_2013_Complexity}.  It has been attacked
from several angles: the replica trick to compute the Boltzmann--Gibbs
distribution \cite{Crisanti_1992_The}, a Kac--Rice \cite{Kac_1943_On,
Rice_1939_The, Fyodorov_2004_Complexity} procedure (similar to the
Fadeev--Popov integral) to compute the number of saddle-points of the energy
function \cite{Crisanti_1995_Thouless-Anderson-Palmer}, and gradient-descent
(or more generally Langevin) dynamics starting from a high-energy configuration
\cite{Cugliandolo_1993_Analytical}. Thanks to the simplicity of the energy, all
these approaches yield analytic results in the large-$N$ limit.

In this paper we extend the study to complex variables: we shall take
$z\in\mathbb C^N$ and $J$ to be a symmetric tensor whose elements are
\emph{complex} normal, with $\overline{|J|^2}=p!/2N^{p-1}$ and
$\overline{J^2}=\kappa\overline{|J|^2}$ for complex parameter $|\kappa|<1$. The
constraint remains $z^Tz=N$.

The motivations for this paper are of three types. On the practical side, there
are indeed situations in which complex variables appear naturally in disordered
problems: such is the case in which the variables are \emph{phases}, as in
random laser problems \cite{Antenucci_2015_Complex}. Quiver Hamiltonians---used
to model black hole horizons in the zero-temperature limit---also have a
Hamiltonian very close to ours \cite{Anninos_2016_Disordered}.  A second reason
is that, as we know from experience, extending a real problem to the complex
plane often uncovers underlying simplicity that is otherwise hidden, shedding
light on the original real problem, e.g., as in the radius of convergence of a
series.

Finally, deforming an integral in $N$ real variables to a surface of dimension $N$ in
$2N$-dimensional complex space has turned out to be necessary for correctly
defining and analyzing path integrals with complex action (see
\cite{Witten_2010_A, Witten_2011_Analytic}), and as a useful palliative for the
sign problem \cite{Cristoforetti_2012_New, Tanizaki_2017_Gradient,
Scorzato_2016_The}.  In order to do this correctly, features of landscape
of the action in complex space---such as the relative position of saddles and the existence of Stokes lines joining them---must be understood. Such landscapes are in general not random: here
we follow the standard strategy of computer science by understanding the
generic features of random instances, expecting that this sheds light on
practical, nonrandom problems.

Returning to our problem, the spherical constraint is enforced using the method
of Lagrange multipliers: introducing $\epsilon\in\mathbb C$, our constrained
energy is
\begin{equation} \label{eq:constrained.hamiltonian}
  H = H_0+\frac p2\epsilon\left(N-\sum_i^Nz_i^2\right).
\end{equation}
One might balk at the constraint $z^Tz=N$---which could appropriately be called
a \emph{hyperbolic} constraint---by comparison with $z^\dagger z=N$. The reasoning
behind the choice is twofold.

First, we seek draw conclusions from our model that are applicable to generic
holomorphic functions without any symmetry. Samples of $H_0$ nearly provide
this, save for a single anomaly: the value of the energy and its gradient at
any point $z$ correlate along the $z$ direction, with $\overline{H_0\partial
H_0}\propto \overline{H_0(\partial H_0)^*}\propto z$. This anomalous direction
should thus be forbidden, and the constraint surface $z^Tz=N$ accomplishes this.

Second, taking the constraint to be the level set of a holomorphic function
means the resulting configuration space is a \emph{bone fide} complex manifold,
and therefore permits easy generalization of the integration techniques
referenced above. The same cannot be said for the space defined by $z^\dagger
z=N$, which is topologically the $(2N-1)$-sphere and cannot admit a complex
structure.

Imposing the constraint with  a holomorphic function
makes the resulting configuration space a \emph{bone fide} complex manifold, which is, as we mentioned, the
situation we wish to model. The same cannot be said for the space defined by $z^\dagger
z=N$, which is topologically the $(2N-1)$-sphere, does not admit a complex
structure, and thus yields a trivial structure of saddles.
However, we will introduce the bound $r^2\equiv z^\dagger z/N\leq R^2$
on the `radius' per spin as a device to classify saddles.   We shall see that this
`radius' $r$ and its upper bound $R$ are insightful knobs in our present
problem, revealing structure as they are varied. Note that taking $R=1$ reduces
the problem to that of the ordinary $p$-spin.

The critical points are of $H$ given by the solutions to
\begin{equation} \label{eq:polynomial}
  \frac{p}{p!}\sum_{j_1\cdots j_{p-1}}^NJ_{ij_1\cdots j_{p-1}}z_{j_1}\cdots z_{j_{p-1}}
  = p\epsilon z_i
\end{equation}
for all $i=\{1,\ldots,N\}$, which for fixed $\epsilon$ is a set of $N$
equations of degree $p-1$, to which one must add the constraint.  In this sense
this study also provides a complement to the work on the distribution of zeroes
of random polynomials \cite{Bogomolny_1992_Distribution}, which are for $N=1$
and $p\to\infty$.  We see from \eqref{eq:polynomial} that at any critical
point $\epsilon=H_0/N$, the average energy.

Since $H$ is holomorphic, any critical point of $\operatorname{Re}H$ is also
one of $\operatorname{Im}H$, and therefore of $H$ itself. Writing $z=x+iy$ for
$x,y\in\mathbb R^N$, $\operatorname{Re}H$ can be considered a real-valued
function of $2N$ real variables. The number of critical points of $H$ is thus given by the
usual Kac--Rice formula applied to $\operatorname{Re}H$:
\begin{equation} \label{eq:real.kac-rice}
  \begin{aligned}
    \mathcal N&(\kappa,\epsilon,R)
      = \int dx\,dy\,\delta(\partial_x\operatorname{Re}H)\delta(\partial_y\operatorname{Re}H) \\
      &\hspace{6pc}\times\left|\det\begin{bmatrix}
          \partial_x\partial_x\operatorname{Re}H & \partial_x\partial_y\operatorname{Re}H \\
          \partial_y\partial_x\operatorname{Re}H & \partial_y\partial_y\operatorname{Re}H
        \end{bmatrix}\right|.
  \end{aligned}
\end{equation}
This expression is to be averaged over $J$ to give the complexity $\Sigma$ as
$N \Sigma= \overline{\log\mathcal N}$, a calculation that involves the replica
trick. Based on the experience from similar problems \cite{Castellani_2005_Spin-glass},  the
\emph{annealed approximation} $N \Sigma \sim \log \overline{ \mathcal N}$ is
expected to be exact wherever the complexity is positive.

The Cauchy--Riemann equations may be used to write \eqref{eq:real.kac-rice} in
a manifestly complex way.  With the Wirtinger derivative
$\partial\equiv\frac12(\partial_x-i\partial_y)$, one can write
$\partial_x\operatorname{Re}H=\operatorname{Re}\partial H$ and
$\partial_y\operatorname{Re}H=-\operatorname{Im}\partial H$. Carrying these
transformations through, one finds
\begin{equation} \label{eq:complex.kac-rice}
  \begin{aligned}
    \mathcal N&(\kappa,\epsilon,r)
      = \int dx\,dy\,\delta(\operatorname{Re}\partial H)\delta(\operatorname{Im}\partial H) \\
                &\hspace{6pc}\times\left|\det\begin{bmatrix}
            \operatorname{Re}\partial\partial H & -\operatorname{Im}\partial\partial H \\
            -\operatorname{Im}\partial\partial H & -\operatorname{Re}\partial\partial H
          \end{bmatrix}\right| \\
      &= \int dx\,dy\,\delta(\operatorname{Re}\partial H)\delta(\operatorname{Im}\partial H)
        \left|\det[(\partial\partial H)^\dagger\partial\partial H]\right| \\
      &= \int dx\,dy\,\delta(\operatorname{Re}\partial H)\delta(\operatorname{Im}\partial H)
        |\det\partial\partial H|^2.
  \end{aligned}
\end{equation}
This gives three equivalent expressions for the determinant of the Hessian: as
that of a $2N\times 2N$ real symmetric matrix, that of the $N\times N$ Hermitian
matrix $(\partial\partial H)^\dagger\partial\partial H$, or the norm squared of
that of the $N\times N$ complex symmetric matrix $\partial\partial H$.

These equivalences belie a deeper connection between the spectra of the
corresponding matrices. Each positive eigenvalue of the real matrix has a
negative partner. For each such pair $\pm\lambda$, $\lambda^2$ is an eigenvalue
of the Hermitian matrix and $|\lambda|$ is a \emph{singular value} of the
complex symmetric matrix. The distribution of positive eigenvalues of the
Hessian is therefore the same as the distribution of singular values of
$\partial\partial H$, or the distribution of square-rooted eigenvalues of
$(\partial\partial H)^\dagger\partial\partial H$.

A useful property of the Gaussian $J$ is that gradient and Hessian at fixed
energy $\epsilon$ are statistically independent \cite{Bray_2007_Statistics,
Fyodorov_2004_Complexity}, so that the $\delta$-functions and the Hessian may
be averaged independently. First we shall compute the spectrum of the Hessian,
which can in turn be used to compute the determinant. Then we will treat the
$\delta$-functions and the resulting saddle point equations. The results of
these calculations begin around \eqref{eq:bezout}.

The Hessian $\partial\partial H=\partial\partial H_0-p\epsilon I$ is equal to
the unconstrained Hessian with a constant added to its diagonal. The eigenvalue
distribution $\rho$ is therefore related to the unconstrained distribution
$\rho_0$ by a similar shift: $\rho(\lambda)=\rho_0(\lambda+p\epsilon)$. The
Hessian of the unconstrained Hamiltonian is
\begin{equation} \label{eq:bare.hessian}
  \partial_i\partial_jH_0
  =\frac{p(p-1)}{p!}\sum_{k_1\cdots k_{p-2}}^NJ_{ijk_1\cdots k_{p-2}}z_{k_1}\cdots z_{k_{p-2}},
\end{equation}
which makes its ensemble that of Gaussian complex symmetric matrices, when the
anomalous direction normal to the constraint surface is neglected. Given its variances
$\overline{|\partial_i\partial_j H_0|^2}=p(p-1)r^{p-2}/2N$ and
$\overline{(\partial_i\partial_j H_0)^2}=p(p-1)\kappa/2N$, $\rho_0(\lambda)$ is
constant inside the ellipse
\begin{equation} \label{eq:ellipse}
  \left(\frac{\operatorname{Re}(\lambda e^{i\theta})}{r^{p-2}+|\kappa|}\right)^2+
  \left(\frac{\operatorname{Im}(\lambda e^{i\theta})}{r^{p-2}-|\kappa|}\right)^2
  <\frac{p(p-1)}{2r^{p-2}}
\end{equation}
where $\theta=\frac12\arg\kappa$ \cite{Nguyen_2014_The}. The eigenvalue
spectrum of $\partial\partial H$ is therefore constant inside the same ellipse
translated so that its center lies at $-p\epsilon$.  Examples of these
distributions are shown in the insets of Fig.~\ref{fig:spectra}.

The eigenvalue spectrum of the Hessian of the real part is not the
spectrum $\rho(\lambda)$ of $\partial\partial H$, but instead the
square-root eigenvalue spectrum of $(\partial\partial H)^\dagger\partial\partial H$;
in other words, the singular value spectrum $\rho(\sigma)$ of $\partial\partial
H$. When $\kappa=0$ and the elements of $J$ are standard complex normal, this
is a complex Wishart distribution. For $\kappa\neq0$ the problem changes, and
to our knowledge a closed form is not in the literature.  We have worked out an
implicit form for the singular value spectrum using the replica method.

Introducing replicas to bring the partition function into the numerator of the
Green function \cite{Livan_2018_Introduction} gives
\begin{widetext}
  \begin{equation} \label{eq:green.replicas}
    G(\sigma)=\lim_{n\to0}\int d\zeta\,d\zeta^*\,(\zeta_i^{(0)})^*\zeta_i^{(0)}
      \exp\left\{
      \frac12\sum_\alpha^n\left[(\zeta_i^{(\alpha)})^*\zeta_i^{(\alpha)}\sigma
        -\operatorname{Re}\left(\zeta_i^{(\alpha)}\zeta_j^{(\alpha)}\partial_i\partial_jH\right)
      \right]
    \right\},
  \end{equation}
  with sums taken over repeated Latin indices. The average is then made over
  $J$ and Hubbard--Stratonovich is used to change variables to the replica matrices
  $N\alpha_{\alpha\beta}=(\zeta^{(\alpha)})^\dagger\zeta^{(\beta)}$ and
  $N\chi_{\alpha\beta}=(\zeta^{(\alpha)})^T\zeta^{(\beta)}$, and a series of
  replica vectors. The replica-symmetric ansatz leaves all replica vectors
  zero, and $\alpha_{\alpha\beta}=\alpha_0\delta_{\alpha\beta}$,
  $\chi_{\alpha\beta}=\chi_0\delta_{\alpha\beta}$. The result is
  \begin{equation}\label{eq:green.saddle}
    \overline G(\sigma)=N\lim_{n\to0}\int d\alpha_0\,d\chi_0\,d\chi_0^*\,\alpha_0
    \exp\left\{nN\left[
      1+\frac{p(p-1)}{16}r^{p-2}\alpha_0^2-\frac{\alpha_0\sigma}2+\frac12\log(\alpha_0^2-|\chi_0|^2)
      +\frac p4\operatorname{Re}\left(\frac{(p-1)}8\kappa^*\chi_0^2-\epsilon^*\chi_0\right)
    \right]\right\}.
    \nonumber % He's too long, and we don't cite him (now)!
  \end{equation}
\end{widetext}

\begin{figure}
  \centering

  \includegraphics{spectra_00.pdf}
  \includegraphics{spectra_05.pdf}\\
  \includegraphics{spectra_10.pdf}
  \includegraphics{spectra_15.pdf}

  \caption{
    Eigenvalue and singular value spectra of the Hessian $\partial\partial H$
    of the $3$-spin model with $\kappa=\frac34e^{-i3\pi/4}$. Pictured
    distributions are for critical points at `radius' $r=\sqrt{5/4}$ and with
    energy per spin (a) $\epsilon=0$, (b)
    $\epsilon=-\frac12|\epsilon_{\mathrm{th}}|$, (c)
    $\epsilon=-|\epsilon_{\mathrm{th}}|$, and (d)
    $\epsilon=-\frac32|\epsilon_{\mathrm{th}}|$. The shaded region of each
    inset shows the support of the eigenvalue distribution \eqref{eq:ellipse}.
    The solid line on each plot shows the distribution of singular values
    \eqref{eq:spectral.density}, while the overlaid histogram shows the
    empirical distribution from $2^{10}\times2^{10}$ complex normal matrices
    with the same covariance and diagonal shift as $\partial\partial H$.
  } \label{fig:spectra}
\end{figure}

The argument of the exponential has several saddles. The solutions $\alpha_0$
are the roots of a sixth-order polynomial, and the root with the smallest value
of $\operatorname{Re}\alpha_0$ gives the correct solution in all the cases we
studied. A detailed analysis of the saddle point integration is needed to
understand why this is so. Evaluated at such a solution, the density of
singular values follows from the jump across the cut, or
\begin{equation} \label{eq:spectral.density}
  \rho(\sigma)=\frac1{i\pi N}\left(
    \lim_{\operatorname{Im}\sigma\to0^+}\overline G(\sigma)
    -\lim_{\operatorname{Im}\sigma\to0^-}\overline G(\sigma)
  \right)
\end{equation}
Examples can be seen in Fig.~\ref{fig:spectra} compared with numeric
experiments.

The formation of a gap in the singular value spectrum naturally corresponds to
the origin leaving the support of the eigenvalue spectrum.  Weyl's theorem
requires that the product over the norm of all eigenvalues must not be greater
than the product over all singular values \cite{Weyl_1912_Das}.  Therefore, the
absence of zero eigenvalues implies the absence of zero singular values. The
determination of the threshold energy---the energy at which the distribution
of singular values becomes gapped---is reduced to the geometry problem of
determining when the boundary of the ellipse defined in \eqref{eq:ellipse}
intersects the origin, and yields
\begin{equation} \label{eq:threshold.energy}
  |\epsilon_{\mathrm{th}}|^2
  =\frac{p-1}{2p}\frac{(1-|\delta|^2)^2r^{p-2}}
  {1+|\delta|^2-2|\delta|\cos(\arg\kappa+2\arg\epsilon)}
\end{equation}
for $\delta=\kappa r^{-(p-2)}$. Notice that the threshold depends on both the
energy per spin $\epsilon$ on the `radius' $r$ of the saddle.

We will now address the $\delta$-functions of \eqref{eq:complex.kac-rice}.
These are converted to exponentials by the introduction of auxiliary fields
$\hat z=\hat x+i\hat y$.  The average over $J$ can then be performed. A
generalized Hubbard--Stratonovich allows a change of variables from the $4N$
original and auxiliary fields to eight bilinears defined by $Nr=z^\dagger z$,
$N\hat r=\hat z^\dagger\hat z$, $Na=\hat z^\dagger z$, $Nb=\hat z^Tz$, and
$N\hat c=\hat z^T\hat z$ (and their conjugates). The result, to leading order
in $N$, is
\begin{equation} \label{eq:saddle}
    \overline{\mathcal N}(\kappa,\epsilon,R)
        = \int dr\,d\hat r\,da\,da^*db\,db^*d\hat c\,d\hat c^*e^{Nf(r,\hat r,a,b,\hat c)},
\end{equation}
where the argument of the exponential is
\begin{widetext}
  \begin{equation}
    f=2+\frac12\log\det\frac12\begin{bmatrix}
      1 & r & b & a \\
      r & 1 & a^* & b^* \\
      b & a^* & \hat c & \hat r \\
      a & b^* & \hat r & \hat c^*
    \end{bmatrix}
    +\int d\lambda\,d\lambda^*\rho(\lambda)\log|\lambda|^2
    +p\operatorname{Re}\left\{
      \frac18\left[\hat rr^{p-1}+(p-1)|b|^2r^{p-2}+\kappa(\hat c^*+(p-1)a^2)\right]-\epsilon a
    \right\}.
  \end{equation}
  The spectrum $\rho$ is given in \eqref{eq:ellipse} and is dependant on $r$ alone. This function has an
  extremum in $\hat r$, $a$, $b$, and $\hat c$ at which its value is
  \begin{equation} \label{eq:free.energy.a}
      f=1+\frac12\log\left(\frac4{p^2}\frac{r^2-1}{r^{2(p-1)}-|\kappa|^2}\right)+\int d\lambda\,\rho(\lambda)\log|\lambda|^2
      -2C_+[\operatorname{Re}(\epsilon e^{-i\theta})]^2-2C_-[\operatorname{Im}(\epsilon e^{-i\theta})]^2,
  \end{equation}
\end{widetext}
where $\theta=\frac12\arg\kappa$ and
\begin{equation}
  C_{\pm}=\frac{r^p(1+p(r^2-1))\mp r^2|\kappa|}{r^{2p}\pm(p-1)r^p(r^2-1)|\kappa|-r^2|\kappa|^2}.
\end{equation}
Notice that level sets of $f$ in energy $\epsilon$ also give ellipses, but of
different form from the ellipse in \eqref{eq:ellipse}.

This expression is maximized for $r=R$, its value at the boundary, for
all values of $\kappa$ and $\epsilon$. Evaluating the complexity at this
saddle, in the limit of unbounded spins, gives
\begin{equation} \label{eq:bezout}
  \lim_{R\to\infty}\log\overline{\mathcal N}(\kappa,\epsilon,R)
  =N\log(p-1).
\end{equation}
This is, to leading order, precisely the Bézout bound, the maximum number of
solutions to $N$ equations of degree $p-1$ \cite{Bezout_1779_Theorie}. That we
saturate this bound is perhaps not surprising, since the coefficients of our
polynomial equations \eqref{eq:polynomial} are complex and have no symmetries.
Reaching Bézout in \eqref{eq:bezout} is not our main result, but it provides a
good check.  Analogous asymptotic scaling has been found for the number of pure
Higgs states in supersymmetric quiver theories \cite{Manschot_2012_From}.

\begin{figure}[htpb]
  \centering
  \includegraphics{complexity.pdf}
  \caption{
    The complexity of the 3-spin model as a function of the maximum `radius'
    $R$ at zero energy and several values of $\kappa$. The dashed line shows
    $\frac12\log(p-1)$, while the dotted shows $\log(p-1)$.
  } \label{fig:complexity}
\end{figure}

For finite $R$, everything is analytically tractable at $\epsilon=0$:
\begin{equation} \label{eq:complexity.zero.energy}
  \Sigma(\kappa,0,R)
  =\log(p-1)-\frac12\log\left(\frac{1-|\kappa|^2R^{-4(p-1)}}{1-R^{-4}}\right).
\end{equation}
This is plotted as a function of $R$ for several values of $\kappa$ in
Fig.~\ref{fig:complexity}. For any $|\kappa|<1$, the complexity goes to
negative infinity as $R\to1$, i.e., as the spins are restricted to the reals.
This is natural, since volume of configuration space vanishes in this limit
like $(R^2-1)^N$. However, when the result is analytically continued to
$\kappa=1$ (which corresponds to real $J$) something novel occurs: the
complexity has a finite value at $R=1$. This implies a $\delta$-function
density of critical points on the $r=1$ (or $y=0$) boundary.  The number of
critical points contained there is
\begin{equation}
  \lim_{R\to1}\lim_{\kappa\to1}\log\overline{\mathcal N}(\kappa,0,R)
  = \frac12N\log(p-1),
\end{equation}
half of \eqref{eq:bezout} and corresponding precisely to the number of critical
points of the real $p$-spin model. (Note the role of conjugation symmetry,
already underlined in \cite{Bogomolny_1992_Distribution}.) The full
$\epsilon$-dependence of the real $p$-spin is recovered by this limit as
$\epsilon$ is varied.

\begin{figure}[b]
  \centering
  \includegraphics{desert.pdf}
  \caption{
    The value of bounding `radius' $R$ for which $\Sigma(\kappa,\epsilon,R)=0$ as a
    function of (real) energy per spin $\epsilon$  for the 3-spin model at
    several values of $\kappa$. Above each line the complexity is positive and
    critical points proliferate, while below it the complexity is negative and
    critical points are exponentially suppressed. The dotted black lines show
    the location of the ground and highest exited state energies for the real
    3-spin model.
  } \label{fig:desert}
\end{figure}

In the thermodynamic limit, \eqref{eq:complexity.zero.energy} implies that most
critical points are concentrated at infinite radius $r$. For finite $N$ the
average radius of critical points is likewise finite. By differentiating
$\overline{\mathcal N}$ with respect to $R$ and normalizing, one obtains the
distribution of critical points as a function of $r$. This yields an average
radius proportional to $N^{1/4}$.  One therefore expects typical critical
points to have a norm that grows modestly with system size.

These qualitative features carry over to nonzero $\epsilon$. In
Fig.~\ref{fig:desert} we show that for $\kappa<1$ there is always a gap in $r$
close to one in which solutions are exponentially suppressed. When
$\kappa=1$---the analytic continuation to the real computation---the situation
is more interesting. In the range of energies where there are real solutions
this gap closes, which is only possible if the density of solutions diverges at
$r=1$. Outside this range, around `deep' real energies where real solutions are
exponentially suppressed, the gap remains. A moment's thought tells us that
this is necessary: otherwise a small perturbation of the $J$s could produce
an unusually deep solution to the real problem, in a region where this should
not happen.

\begin{figure}[t]
  \centering

  \includegraphics{threshold_2000.pdf}
  \includegraphics{threshold_1325.pdf} \\
  \includegraphics{threshold_1125.pdf}
  \includegraphics{threshold_1000.pdf}

  \caption{
    Energies at which states exist (green shaded region) and threshold energies
    (black solid line) for the 3-spin model with
    $\kappa=\frac34e^{-i3\pi/4}$ and (a) $r=\sqrt2$, (b) $r=\sqrt{1.325}$, (c) $r=\sqrt{1.125}$,
    and (d) $r=1$. No shaded region is shown in (d) because no states exist at
    any energy.
  } \label{fig:eggs}
\end{figure}

The relationship between the threshold and ground, or extremal, state energies
is richer than in the real case.  In Fig.~\ref{fig:eggs} these are shown in the
complex-$\epsilon$ plane for several examples. Depending on the parameters, the
threshold might have a smaller or larger magnitude than the extremal state, or
cross as a function of complex argument.  For sufficiently large $r$ the
threshold is always at a larger magnitude. If this were to happen in the real
case, it would likely imply our replica symmetric computation were unstable,
since having a ground state above the threshold implies a ground state Hessian
with many negative eigenvalues, a contradiction. However, this is not an
contradiction in the complex case, where the energy is not bounded from below.
The relationship between the threshold, i.e., where the gap appears, and the
dynamics of, e.g., a minimization algorithm, deformed integration cycle, or
physical dynamics, are a problem we hope to address in future work.

 This paper provides a first step towards the study of complex landscapes with
 complex variables. The next obvious step is to study the topology of the
 critical points, the sets reached following  gradient descent (the
 Lefschetz thimbles), and ascent (the anti-thimbles) \cite{Witten_2010_A,
 Witten_2011_Analytic, Cristoforetti_2012_New, Behtash_2017_Toward,
 Scorzato_2016_The}, which act as constant-phase integrating `contours.'
 Locating and counting the saddles that are joined by gradient lines---the
 Stokes points, which play an important role in the theory---is also well within
 reach of the present-day spin-glass literature techniques. We anticipate
 that the threshold level, where the system develops a mid-spectrum gap, plays
 a crucial role in determining whether these Stokes points proliferate under
 some continuous change of parameters.

\begin{acknowledgments}
  We wish to thank Alexander Altland, Satya Majumdar and Gregory Schehr for a useful suggestions.
  JK-D and JK are supported by the Simons Foundation Grant No.~454943.
\end{acknowledgments}

\bibliography{bezout}

\end{document}