Big rewrite for clarity.

3 files changed, 202 insertions, 192 deletions

diff --git a/bezout.tex b/bezout.tex
index f950b8c..f050707 100644
--- a/bezout.tex
+++ b/bezout.tex
@@ -42,11 +42,11 @@
 
 \maketitle
 
-Spin-glasses have long been considered the paradigm of many variable `complex
-landscapes,' a subject that includes neural networks and optimization problems,
-most notably 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
+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}
@@ -58,11 +58,10 @@ 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 the
-gradient-descent---or more generally Langevin---dynamics staring 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.
+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
@@ -77,7 +76,7 @@ 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, sheding
+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.
 
@@ -86,15 +85,16 @@ $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, the features of landscape
+Scorzato_2016_The}.  In order to do this correctly, features of landscape
 of the action in complex space---like the relative position of its
 saddles---must be understood. Such landscapes are in general not random: here
 we propose to follow the strategy of computer science of understanding the
-generic features of random instances, expecting that this sheds light on 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 energy is
+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}
@@ -102,35 +102,35 @@ 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 would be 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 must be neglected if we are to draw conclusions about
-generic functions, and the constraint surface $z^Tz=N$ is the unique surface
+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
+must be neglected, and the constraint surface $z^Tz=N$ is the unique surface
 whose normal is parallel to $z$ and which contains the configuration space of
 the real $p$-spin model as a subspace.
 
 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
+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.
 
 A consequence of the constraint is that the model's configuration space is not
-compact, nor is its energy bounded. This is not necessarily a problem, as many
+compact, nor is its energy bounded. This is not necessarily problematic, as many
 related problems have similar properties but are concerned with subspaces on
-which the energy is bounded. (In fact, identifying the appropriate subspace on
-which to define one's model often requires the study of critical points in the
-whole space.) Where it might be a problem, we introduce the additional
-constraint $z^\dagger z\leq Nr^2$. The resulting configuration space is a complex
-manifold with boundary. We shall see that the `radius' $r$ proves an insightful
-knob in our present problem, revealing structure as it is 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 the set of equations
+which the energy is bounded. (In fact, identifying the appropriate subspace
+often requires the study of critical points in the whole space.) Where it might
+become problematic, we introduce an additional constraint that bounds the
+`radius' per spin $r^2\equiv z^\dagger z/N\leq R^2$.  The resulting
+configuration space is a complex manifold with boundary. We shall see that the
+`radius' $r$ and its upper bound $R$ prove to be 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
@@ -140,7 +140,7 @@ 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/N$, the average energy.
+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
@@ -149,7 +149,7 @@ function of $2N$ real variables. The number of critical points of $H$ is thus gi
 usual Kac--Rice formula applied to $\operatorname{Re}H$:
 \begin{equation} \label{eq:real.kac-rice}
   \begin{aligned}
-    \mathcal N&(\kappa,\epsilon,r)
+    \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 \\
@@ -157,11 +157,18 @@ usual Kac--Rice formula applied to $\operatorname{Re}H$:
         \end{bmatrix}\right|.
   \end{aligned}
 \end{equation}
-The Cauchy--Riemann equations may be used to write this in a manifestly complex
-way.  With the Wirtinger derivative $\partial=\frac12(\partial_x-i\partial_y)$,
-one can write $\partial_x\operatorname{Re}H=\operatorname{Re}\partial H$ and
+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. In most of the parameter space that we shall study here, the
+\emph{annealed approximation} $N \Sigma \sim \log \overline{ \mathcal N}$ is
+exact. 
+
+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, we have
+transformations through, one finds
 \begin{equation} \label{eq:complex.kac-rice}
   \begin{aligned}
     \mathcal N&(\kappa,\epsilon,r)
@@ -177,8 +184,9 @@ transformations through, we have
   \end{aligned}
 \end{equation}
 This gives three equivalent expressions for the determinant of the Hessian: as
-that of a $2N\times 2N$ real matrix, that of an $N\times N$ Hermitian matrix,
-or the norm squared of that of an $N\times N$ complex symmetric matrix.
+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
@@ -189,52 +197,13 @@ 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$.
 
-The expression \eqref{eq:complex.kac-rice} 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. In most of the parameter-space that we shall
-study here, the \emph{annealed approximation} $N \Sigma \sim \log \overline{
-\mathcal N}$ is exact. 
-
 A useful property of the Gaussian $J$ is that gradient and Hessian at fixed
-$\epsilon$ are statistically independent \cite{Bray_2007_Statistics,
+energy $\epsilon$ are statistically independent \cite{Bray_2007_Statistics,
 Fyodorov_2004_Complexity}, so that the $\delta$-functions and the Hessian may
-be averaged independently. The $\delta$-functions are converted to exponentials
-by the introduction of auxiliary fields $\hat z=\hat x+i\hat y$.  The average
-of those factors 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 $Na=z^\dagger z$, $N\hat a=\hat z^\dagger\hat z$, $N\hat c=\hat z^T\hat z$, $Nb=\hat z^\dagger z$, and $Nd=\hat z^Tz$ (and their
-conjugates). The result, to leading order in $N$, is
-\begin{equation} \label{eq:saddle}
-    \overline{\mathcal N}(\kappa,\epsilon,r)
-        = \int da\,d\hat a\,db\,db^*d\hat c\,d\hat c^*dd\,dd^*e^{Nf(a,\hat a,b,\hat c,d)},
-\end{equation}
-where the argument of the exponential is
-\begin{widetext}
-  \begin{equation}
-    f=2+\frac12\log\det\frac12\begin{bmatrix}
-      1 & a & d & b \\
-      a & 1 & b^* & d^* \\
-      d & b^* & \hat c & \hat a \\
-      b & d^* & \hat a & \hat c^*
-    \end{bmatrix}
-    +\int d\lambda\,\rho(\lambda)\log|\lambda|^2
-    +p\operatorname{Re}\left\{
-      \frac18\left[\hat aa^{p-1}+(p-1)|d|^2a^{p-2}+\kappa(\hat c^*+(p-1)b^2)\right]-\epsilon b
-    \right\}.
-  \end{equation}
-  The integral of the distribution $\rho$ of eigenvalues of $\partial\partial
-  H$ comes from the Hessian and is dependant on $a$ alone. This function has an
-  extremum in $\hat a$, $b$, $\hat c$, and $d$ at which its value is
-  \begin{equation} \label{eq:free.energy.a}
-      f(a)=1+\frac12\log\left(\frac4{p^2}\frac{a^2-1}{a^{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{a^p(1+p(a^2-1))\mp a^2|\kappa|}{a^{2p}\pm(p-1)a^p(a^2-1)|\kappa|-a^2|\kappa|^2}.
-\end{equation}
-This leaves a single parameter, $a$, which dictates the norm of $z$. 
+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
@@ -246,28 +215,28 @@ Hessian of the unconstrained Hamiltonian is
   =\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
-direction along the constraint is neglected. Given its variances
-$\overline{|\partial_i\partial_j H_0|^2}=p(p-1)a^{p-2}/2N$ and
+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})}{a^{p-2}+|\kappa|}\right)^2+
-  \left(\frac{\operatorname{Im}(\lambda e^{i\theta})}{a^{p-2}-|\kappa|}\right)^2
-  <\frac{p(p-1)}{2a^{p-2}}
+  \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 different from the
-spectrum $\rho(\lambda)$ of $\partial\partial H$, but rather equivalent to the
+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 this spectrum using the replica method.
+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
@@ -280,25 +249,24 @@ Green function \cite{Livan_2018_Introduction} gives
       \right]
     \right\},
   \end{equation}
-  with sums taken over repeated Latin indices.  The average is then made over
+  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)})^*\cdot\zeta^{(\beta)}$ and
-  $N\chi_{\alpha\beta}=\zeta^{(\alpha)}\cdot\zeta^{(\beta)}$ and a series of
-  replica vectors. The replica-symmetric ansatz leaves all off-diagonal
-  elements and vectors zero, and
-  $\alpha_{\alpha\beta}=\alpha_0\delta_{\alpha\beta}$,
+  $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}a^{p-2}\alpha_0^2-\frac{\alpha_0\sigma}2+\frac12\log(\alpha_0^2-|\chi_0|^2)
+      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}[b]
+\begin{figure}
   \centering
 
   \includegraphics{fig/spectra_0.0.pdf}
@@ -307,25 +275,27 @@ Green function \cite{Livan_2018_Introduction} gives
   \includegraphics{fig/spectra_1.5.pdf}
 
   \caption{
-    Eigenvalue and singular value spectra of the matrix $\partial\partial H$
-    for $p=3$, $a=\frac54$, and $\kappa=\frac34e^{-i3\pi/4}$ with (a)
-    $\epsilon=0$, (b) $\epsilon=-\frac12|\epsilon_{\mathrm{th}}|$, (c)
+    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. The solid line on each
-    plot shows the distribution of singular values, 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$.
+    $\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$ in all the cases we studied gives the correct
-solution. A detailed analysis of the saddle point integration is needed to
-understand why this is so. Given such $\alpha_0$, the density of singular
-values follows from the jump across the cut, or
+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)
@@ -335,70 +305,109 @@ values follows from the jump across the cut, or
 Examples can be seen in Fig.~\ref{fig:spectra} compared with numeric
 experiments.
 
-The transition from a one-cut to two-cut 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 then reduced to a
-geometry problem, and yields
+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)^2a^{p-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 a^{-(p-2)}$.
+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}.
 
-Given $\rho$, the integral in \eqref{eq:free.energy.a} may be preformed for
-arbitrary $a$. The resulting expression is maximized for $a=r^2$ for all
-values of $\kappa$ and $\epsilon$. Evaluating the complexity at this saddle, in
-the limit of unbounded spins, gives
+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)
+  \lim_{R\to\infty}\log\overline{\mathcal N}(\kappa,\epsilon,R)
   =N\log(p-1).
 \end{equation}
-This is, to this order,  precisely the Bézout bound, the maximum number of
-solutions that $N$ equations of degree $p-1$ may have
-\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}.
+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{fig/complexity.pdf}
   \caption{
-    The complexity of the 3-spin model at $\epsilon=0$ as a function of
-    the maximum `radius' $r$ at several values of $\kappa$. The dashed line
-    shows $\frac12\log(p-1)$, while the dotted shows $\log(p-1)$.
+    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$, giving
+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).
+  \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, given that the $y$ contribution to
-the volume shrinks to zero as that of an $N$-dimensional sphere $\sum_i y_i^2=N(r^2-1)$ with volume
-$\sim(r^2-1)^N$.  However, when the result is analytically continued to
+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$. Since the $r$-dependence gives a
-cumulative count, this implies a $\delta$-function density of critical points
-along the line $y=0$.  The number of critical points contained within is
+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)
+  \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
+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.
 
@@ -406,7 +415,7 @@ $\epsilon$ is varied.
   \centering
   \includegraphics{fig/desert.pdf}
   \caption{
-    The value of `radius' $r$ for which $\Sigma(\kappa,\epsilon,r)=0$ as a
+    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
@@ -416,25 +425,25 @@ $\epsilon$ is varied.
   } \label{fig:desert}
 \end{figure}
 
-In the thermodynamic limit \eqref{eq:complexity.zero.energy} implies that most
-critical points are concentrated at infinite radius, i.e., at complex vectors with
-very large squared norm. For finite $N$ the average radius of critical points is likewise finite. By differentiating $\overline{\mathcal N}$ with
-respect to $r$ and normalizing, one has the distribution
-of critical points as a function of $r$. The average radius this yields is
-$\propto N^{1/4}+O(N^{-3/4})$. One therefore expects typical
-critical points to have a norm that grows modestly with system size.
+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 of $r$
-close to one for which there are no solutions. 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$.  Another
-remarkable feature of this limit is that there is still a gap without solutions
-around `deep' real energies where there is no real solution. A moment's thought
-tells us that this is a necessity: 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.
+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
@@ -447,37 +456,38 @@ this should not happen.
   \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) $a=2$, (b) $a=1.325$, (c) $a=1.125$,
-    and (d) $a=1$. No shaded region is shown in (d) because no states exist at
+    $\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
+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 always come at smaller magnitude than the extremal state, or
-always come at larger magnitude, or cross as a function of complex argument.
-For sufficiently large $r$ the threshold always comes at larger magnitude than
-the extremal state. If this were to happen in the real case, it would likely
-imply our replica symmetric computation is 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 obvious contradiction in
-the complex case. The relationship between the threshold, i.e., where the gap
-appears, and the dynamics of, e.g., a minimization algorithm or physical
-dynamics, are a problem we hope to address in future work.
-
- This paper provides a first step towards the study of a complex landscape with
- complex variables. The next obvious one is to study the topology of 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, will
- play a crucial role as it does in the real case.
+ 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.
diff --git a/fig/complexity.pdf b/fig/complexity.pdf
index f9336bb..cd4864d 100644
--- a/fig/complexity.pdf
+++ b/fig/complexity.pdf
diff --git a/fig/desert.pdf b/fig/desert.pdf
index 08d8f41..c9e03df 100644
--- a/fig/desert.pdf
+++ b/fig/desert.pdf