A bunch of writing in the p-spin sections.

1 files changed, 181 insertions, 283 deletions

diff --git a/stokes.tex b/stokes.tex
index 84757b9..b6f3e36 100644
--- a/stokes.tex
+++ b/stokes.tex
@@ -29,13 +29,13 @@
   In this paper we follow up the study of `complex complex landscapes'
   \cite{Kent-Dobias_2021_Complex}, rugged landscapes of many complex
   variables.  Unlike real landscapes, there is no useful classification of
-  saddles by index. Instead, the spectrum at critical points determines their
+  saddles by index. Instead, the spectrum at stationary points determines their
   tendency to trade topological numbers under analytic continuation of the
   theory. These trades, which occur at Stokes points, proliferate when the
   spectrum includes marginal directions and are exponentially suppressed
   otherwise. This gives a direct interpretation of the `threshold' energy---which
   in the real case separates saddles from minima---where the spectrum of
-  typical critical points develops a gap. This leads to different consequences
+  typical stationary points develops a gap. This leads to different consequences
   for the analytic continuation of real landscapes with different structures:
   the global minima of ``one step replica-symmetry broken'' landscapes lie
   beyond a threshold and are locally protected from Stokes points, whereas
@@ -161,7 +161,7 @@ action potentially has more stationary points. We'll call $\Sigma$ the set of
   \includegraphics{figs/stationaryPoints.pdf}
 
   \caption{
-    An example of a simple action and its critical points. \textbf{Left:} An
+    An example of a simple action and its stationary points. \textbf{Left:} An
     action $\mathcal S$ for the $N=2$ spherical (or circular) $3$-spin model,
     defined for $s_1,s_2\in\mathbb R$ on the circle $s_1^2+s_2^2=2$ by
     $\mathcal S(s_1,s_2)=-1.051s_1^3-1.180s_1^2s_2-0.823s_1s_2^2-1.045s_2^3$. In
@@ -363,7 +363,7 @@ behavior can be seen in Fig.~\ref{fig:1d.stokes}.
 \end{figure}
 
 The prevalence (or not) of Stokes points in a given continuation, and whether
-those that do appear affect the weights of critical points of interest, is a
+those that do appear affect the weights of stationary points of interest, is a
 concern for the analytic continuation of theories. If they do not occur or
 occur order-one times, one could reasonably hope to perform such a procedure.
 If they occur exponentially often in the system size, there is little hope of
@@ -484,7 +484,7 @@ multiplier $\mu$, we define the constrained action
 \end{equation}
 The condition for a stationary point is that $\partial\tilde\mathcal S=0$. This implies that, at any stationary point,
 $\partial\mathcal S=\mu\partial g$. In particular, if $\partial g^T\partial g\neq0$, we find the value for $\mu$ as
-\begin{equation}
+\begin{equation} \label{eq:multiplier}
   \mu=\frac{\partial g^T\partial\mathcal S}{\partial g^T\partial g}
 \end{equation}
 As a condition for a stationary point, this can be intuited as projecting out
@@ -687,7 +687,7 @@ After all the work of decomposing an integral into a sub over thimbles, one
 eventually wants to actually evaluate it. For large $|\beta|$ and in the
 absence of any Stokes points, one can come to a nice asymptotic expression.
 
-Suppose that $\sigma\in\Sigma$ is a critical point at $s_\sigma\in\tilde\Omega$
+Suppose that $\sigma\in\Sigma$ is a stationary point at $s_\sigma\in\tilde\Omega$
 with a thimble $\mathcal J_\sigma$ that is not involved in any upstream Stokes points.
 Define its contribution to the partition function (neglecting the integer
 weight) as
@@ -696,7 +696,7 @@ weight) as
 \end{equation}
 To evaluate this contour integral in the limit of large $|\beta|$, we will make
 use of the saddle point method, since the integral will be dominated by its
-value at and around the critical point, where the real part of the action is by
+value at and around the stationary point, where the real part of the action is by
 construction at its minimum on the thimble and the integrand is therefore
 largest.
 
@@ -706,7 +706,7 @@ We will make a change of coordinates $u(s)\in\mathbb R^D$ such that
 \end{equation}
 \emph{and} the direction of each $\partial u/\partial s$ is along the direction
 of the contour. This is possible because, in the absence of any Stokes points,
-the eigenvectors of the hessian at the critical point associated with positive
+the eigenvectors of the hessian at the stationary point associated with positive
 eigenvalues provide a basis for the thimble. The coordinates $u$ can be real
 because the imaginary part of the action is constant on the thimble, and
 therefore stays with the value it holds at the stationary point, and the real
@@ -764,14 +764,14 @@ However, because we are dealing with a path integral, the directions matter,
 and there is not an absolute value around the determinant. Therefore, we must
 determine the phase that it contributes.
 
-This is difficult in general, but for real critical points it can be reasoned
+This is difficult in general, but for real stationary points it can be reasoned
 out easily. Take the same convention we used earlier, that the direction of
 contours along the real line is in the conventional directions. Then, a
-critical point of index $k$ has $k$ real eigenvectors and $D-k$ purely
+stationary point of index $k$ has $k$ real eigenvectors and $D-k$ purely
 imaginary eigenvectors that contribute to its thimble. The matrix of
 eigenvectors can therefore be written $U=i^kO$ for an orthogonal matrix $O$,
 and with all eigenvectors canonically oriented $\det O=1$. We therefore have
-$\det U=i^k$. As the argument of $\beta$ is changed, we know how the eigenvectors change: by a factor of $e^{-i\phi/2}$. Therefore, the contribution more generally is $(e^{-i\phi/2})^Di^k$. We therefore have, for real critical points of a real action,
+$\det U=i^k$. As the argument of $\beta$ is changed, we know how the eigenvectors change: by a factor of $e^{-i\phi/2}$. Therefore, the contribution more generally is $(e^{-i\phi/2})^Di^k$. We therefore have, for real stationary points of a real action,
 \begin{equation}
   Z_\sigma(\beta)\simeq\left(\frac{2\pi}\beta\right)^{D/2}i^{k_\sigma}\prod_{\lambda_0>0}\lambda_0^{-\frac12}e^{-\beta\mathcal S(s_\sigma)}
 \end{equation}
@@ -959,50 +959,71 @@ for $\delta=C_0/A_0$.
 
 \section{The \textit{p}-spin spherical models}
 
-The $p$-spin spherical models are statistical mechanics models defined by the
-action $\mathcal S=-\beta H$ for the Hamiltonian
+The $p$-spin spherical models are defined by the action
 \begin{equation} \label{eq:p-spin.hamiltonian}
-  H(x)=\sum_{p=2}^\infty\frac{a_p}{p!}\sum_{i_1\cdots i_p}^NJ_{i_1\cdots i_p}x_{i_1}\cdots x_{i_p}
+  \mathcal S(x)=\sum_{p=2}^\infty a_p\mathcal S_p(x)
+\end{equation}
+which is a sum of the `pure' actions
+\begin{equation} \label{eq:pure.p-spin.hamiltonian}
+  \mathcal S_p(x)=\frac1{p!}\sum_{i_1\cdots i_p}J_{i_1\cdots i_p}x_{i_1}\cdots x_{i_p}
 \end{equation}
-where the $x\in\mathbb R^N$ are constrained to lie on the sphere $x^2=N$.  The tensor
-components $J$ are complex normally distributed with zero mean and variances
-$\overline{|J|^2}=p!/2N^{p-1}$ and $\overline{J^2}=\kappa\overline{|J|^2}$, and
-the numbers $a$ define the model. The pure real $p$-spin model has
-$a_i=\delta_{ip}$ and $\kappa=1$.
+The variables $x\in\mathbb R^N$ are constrained to lie on the sphere $x^2=N$,
+making the model $D=N-1$ dimensional. The couplings $J$ form a totally symmetry
+$p$-tensor whose components are normally distributed with zero mean and
+variance $\overline{J^2}=p!/2N^{p-1}$. The `pure' $p$-spin models have
+$a_i=\delta_{ip}$, while the mixed have some more complicated coefficients $a$.
 
 The phase space manifold $\Omega=\{x\mid x^2=N, x\in\mathbb R^N\}$ has a
 complex extension $\tilde\Omega=\{z\mid z^2=N, z\in\mathbb C^N\}$. The natural
-extension of the hamiltonian \eref{eq:p-spin.hamiltonian} to this complex manifold is
-holomorphic. The normal to this manifold at any point $z\in\tilde\Omega$ is
-always in the direction $z$.  The projection operator onto the tangent space of
-this manifold is given by
+extension of the hamiltonian \eref{eq:p-spin.hamiltonian} to this complex
+manifold by replacing $x$ with $z\in\mathbb C^N$ is holomorphic. The normal to
+this manifold at any point $z\in\tilde\Omega$ is always in the direction $z$.
+The projection operator onto the tangent space of this manifold is given by
 \begin{equation}
   P=I-\frac{zz^\dagger}{|z|^2},
 \end{equation}
 where indeed $Pz=z-z|z|^2/|z|^2=0$ and $Pz'=z'$ for any $z'$ orthogonal to $z$.
+When studying stationary points, the constraint can be added to the action in
+using a Lagrange multiplier $\mu$ by writing
+\begin{equation}
+  \tilde\mathcal S(z)=\mathcal S(z)-\frac\mu2(z^Tz-N)
+\end{equation}
+The gradient of the constraint is simple with $\partial g=z$, and \eqref{eq:multiplier} implies that
+\begin{eqnarray}
+  \mu
+  &=\frac1Nz^T\partial\mathcal S
+  =\frac1Nz_{i_p}\sum_{p=2}^\infty a_p\frac p{p!}\sum_{i_1\cdots i_p}J_{i_1\cdots i_p}z_{i_1}\cdots z_{i_{p-1}} \\
+  &=\sum_{p=2}^\infty a_pp\frac{\mathcal S_p(z)}N
+\end{eqnarray}
+which for the pure $p$-spin in particular implies that $\mu=p\epsilon$ for specific energy $\epsilon$.
 
 \subsection{2-spin}
 
-The Hamiltonian of the $2$-spin model is defined for $z\in\mathbb C^N$ by
+The Hamiltonian of the pure $2$-spin model is defined by
 \begin{equation}
-  H_0=\frac12z^TJz.
+  \mathcal S_2(z)=\frac12z^TJz.
 \end{equation}
-$J$ is generically diagonalizable by a complex orthogonal matrix. In a diagonal basis, $J_{ij}=\lambda_i\delta_{ij}$. Then $\partial_i H=\lambda_iz_i$. We will henceforth assume to be working in this basis. To find the critical points, we use the method of Lagrange multipliers. Introducing $\epsilon$, the constrained Hamiltonian is
+$J$ is generically diagonalizable. In a diagonal basis,
+$J_{ij}=\lambda_i\delta_{ij}$. Then $\partial_i H=\lambda_iz_i$. We will
+henceforth assume to be working in this basis. The constrained action is
 \begin{equation}
-  H=H_0+\epsilon(N-z^2)
+  \tilde\mathcal S(z)=\mathcal S_2(z)-\epsilon(z^Tz-N)
 \end{equation}
-As usual, $\epsilon$ is equivalent to the energy per spin at any critical point.
-Critical points must satisfy
+Stationary points must satisfy
 \begin{equation}
-  0=\partial_iH=(\lambda_i-2\epsilon)z_i
+  0=\partial_i\tilde\mathcal S=(\lambda_i-2\epsilon)z_i
 \end{equation}
-which is only possible for $z_i=0$ or $\epsilon=\frac12\lambda_i$. Generically the $\lambda_i$ will all differ, so this can only be satisfied for one $\lambda_i$ at a time, and to be a critical point all other $z_j$ must be zero. In the direction in question,
+which is only possible for $z_i=0$ or $\epsilon=\frac12\lambda_i$. Generically
+the $\lambda_i$ will all differ, so this can only be satisfied for one
+$\lambda_i$ at a time, and to be a stationary point all other $z_j$ must be
+zero. In the direction in question,
 \begin{equation}
   \frac1N\frac12\lambda_iz_i^2=\epsilon=\frac12\lambda_i,
 \end{equation}
-whence $z_i=\pm\sqrt{N}$. Thus there are $2N$ critical points, each corresponding to $\pm$ the cardinal directions in the diagonalized basis.
+whence $z_i=\pm\sqrt{N}$. Thus there are $2N$ stationary points, each
+corresponding to $\pm$ the cardinal directions in the diagonalized basis.
 
-Suppose that two critical points have the same imaginary energy; without loss
+Suppose that two stationary points have the same imaginary energy; without loss
 of generality, assume these are associated with the first and second
 cardinal directions. Since the gradient is proportional to $z$, any components that are
 zero at some time will be zero at all times. The dynamics for the components of
@@ -1014,14 +1035,14 @@ interest assuming all others are zero are
 \end{eqnarray}
 and the same for $z_2$ with all indices swapped.  Since $\Delta=d_1-d_2$ is
 real, if $z_1$ begins real it remains real, with the same for $z_2$. Since the
-critical points are at real $z$, we make this restriction, and find
+stationary points are at real $z$, we make this restriction, and find
 \begin{equation}
   \frac d{dt}(z_1^2+z_2^2)=0 \qquad
   \frac d{dt}\frac{z_2}{z_1}=\Delta\frac{z_2}{z_1}
 \end{equation}
 Therefore $z_2/z_1=e^{\Delta t}$, with $z_1^2+z_2^2=N$ as necessary.  Depending
-on the sign of $\Delta$, $z$ flows from one critical point to the other over
-infinite time. This is a Stokes line, and establishes that any two critical
+on the sign of $\Delta$, $z$ flows from one stationary point to the other over
+infinite time. This is a Stokes line, and establishes that any two stationary
 points in the 2-spin model with the same imaginary energy will possess one.
 These trajectories are plotted in Fig.~\ref{fig:two-spin}.
 
@@ -1034,18 +1055,18 @@ These trajectories are plotted in Fig.~\ref{fig:two-spin}.
          1 / sqrt(1 + exp(- 2 * x)) t '$z_2$'
   \end{gnuplot}
   \caption{
-    The Stokes line in the 2-spin model when the critical points associated
+    The Stokes line in the 2-spin model when the stationary points associated
     with the first and second cardinal directions are brought to the same
     imaginary energy. $\Delta$ is proportional to the difference between the
-    real energies of the first and the second critical point; when $\Delta >0$
+    real energies of the first and the second stationary point; when $\Delta >0$
     flow is from first to second, while when $\Delta < 0$ it is reversed.
   } \label{fig:two-spin}
 \end{figure}
 
-Since they sit at the corners of a simplex, the critical points of the 2-spin
-model are all adjacent: no critical point is separated from another by the
+Since they sit at the corners of a simplex, the stationary points of the 2-spin
+model are all adjacent: no stationary point is separated from another by the
 separatrix of a third. This means that when the imaginary energies of two
-critical points are brought to the same value, their surfaces of constant
+stationary points are brought to the same value, their surfaces of constant
 imaginary energy join.
 
 \begin{eqnarray}
@@ -1117,18 +1138,44 @@ maximum again becomes coherent. These conditions correspond precisely to those
 found for the density of zeros in the 2-spin model found previously
 \cite{Obuchi_2012_Partition-function, Takahashi_2013_Zeros}.
 
-
 \subsection{Pure \textit{p}-spin: where are the saddles?}
 
-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$:
+We studied the distribution of stationary points in the pure $p$-spin models in
+a previous work \cite{Kent-Dobias_2021_Complex}. Here, we will review the
+method and elaborate on some of the results relevant to their analytic
+continuation.
+
+In the real case, the $p$-spin models posses a threshold energy
+$\epsilon_{\mathrm{th}}$, below which there are exponentially many minima
+compared to saddles, and above which vice versa. This threshold persists in a
+more generic form in the complex case, where now the threshold separates mostly
+gapped from mostly ungapped saddles. Since the $p$-spin model has a Hessian
+that consists of a symmetric complex matrix with a shifted diagonal, we can use
+the results of \S\ref{sec:stationary.hessian} appropriately scaled. The variance of the $p$-spin hessian without shift is
+\begin{equation}
+  \overline{|\partial\partial\mathcal S_p|^2}
+  =\frac{p(p-1)(\frac1Nz^\dagger z)^{p-2}}{2N}
+  =\frac{p(p-1)}{2N}(1+Y)^{p-2}
+\end{equation}
+\begin{equation}
+  \overline{(\partial\partial\mathcal S_p)^2}
+  =\frac{p(p-1)(\frac1Nz^Tz)^{p-2}}{2N}
+  =\frac{p(p-1)}{2N}
+\end{equation}
+As expected for a real problem, the two variances coincide when $Y=0$.
+The diagonal shift is $-p\epsilon$. In the language of
+\S\ref{sec:stationary.hessian}, this means that $A_0=p(p-1)(1+Y)^{p-2}/2N$,
+$C_0=p(p-1)2N$, and $\lambda_0=-p\epsilon$.
+
+\begin{equation}
+  |\epsilon_{\mathrm{th}}|^2=\frac{p-1}{2p}
+\end{equation}
+
+The location of stationary points can be determined by the Kac--Rice formula. Any stationary point of the action is a stationary point of the real part of the action, and we can write
 \begin{equation} \label{eq:real.kac-rice}
-  \mathcal N(\kappa,\epsilon,R)
-    = \int dx\,dy\,\delta(\partial_x\operatorname{Re}H)\delta(\partial_y\operatorname{Re}H)
-            \left|\det\operatorname{Hess}_{x,y}\operatorname{Re}H\right|.
+  \mathcal N
+    = \int dx\,dy\,\delta(\partial_x\operatorname{Re}\tilde\mathcal S_p)\delta(\partial_y\operatorname{Re}\tilde\mathcal S_p)
+            \left|\det\operatorname{Hess}_{x,y}\operatorname{Re}\mathcal S_p\right|
 \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
@@ -1136,244 +1183,95 @@ trick. Based on the experience from similar problems \cite{Castellani_2005_Spin-
 \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}
-  \mathcal N(\kappa,\epsilon,r)
-    = \int dx\,dy\,\delta(\operatorname{Re}\partial H)\delta(\operatorname{Im}\partial H)
-    |\det\operatorname{Hess}H|^2.
-\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
-
-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{equation}
-  f=2+\frac12\log\det\frac12\left[\matrix{
-    1 & r & b & a \cr
-    r & 1 & a^* & b^* \cr
-    b & a^* & \hat c & \hat r \cr
-    a & b^* & \hat r & \hat c^*
-  }\right]
-  +\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}
-where $\theta=\frac12\arg\kappa$ and
+As in \S\ref{sec:stationary.hessian}, these can be bright into a manifestly complex form using Cauchy--Riemann relations. This gives
+\begin{eqnarray} \label{eq:real.kac-rice}
+  \mathcal N
+  &=\int dz^*dz\,d\hat z^*d\hat z\,d\eta^*d\eta\,d\gamma^*d\gamma\exp\left\{
+    -\operatorname{Re}\left(
+      \hat z^T\partial\tilde\mathcal S_p+\eta^T\partial\partial\mathcal S_p\gamma
+    \right)
+  \right\} \\
+  &=\int dz^*dz\,d\hat z^*d\hat z\,d\eta^*d\eta\,d\gamma^*\,d\gamma\exp\left\{
+    -\frac12\left(
+      \hat z^T\partial\tilde\mathcal S_p+\hat z^\dagger(\partial\tilde\mathcal S_p)^*+
+      \eta^T\partial\partial\mathcal S_p\gamma+
+      \eta^\dagger(\partial\partial\mathcal S_p)^*\gamma^*
+    \right)
+  \right\}
+\end{eqnarray}
+where $\eta$ and $\gamma$ are Grassmann variables. This can be more
+conveniently studied using the method of superfields. Introducing the
+one-component Grassman variables $\theta$ and $\bar\theta$, define the
+superfield
+\begin{eqnarray}
+  \zeta(i)&=z+\bar\theta(i)\eta^*+\gamma\theta(i)+\hat z\bar\theta(i)\theta(i)
+\end{eqnarray}
+Then the expression for the number of stationary points can be written in a very compact form, as
 \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{figs/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}
+  \mathcal N=\int d\zeta^*d\zeta\,\exp\left\{
+    -\frac12\int d1\,\left(
+      \tilde\mathcal S_p(\zeta(1))+\tilde\mathcal S_p(\zeta^*(1))
+    \right)
+  \right\}
+\end{equation}
+where $d1=d\bar\theta(1)\,d\theta(1)$ is a shorthand for integration over the auxiliary grassman variables.
+This can be related to the previous expression by expansion with respect to
+the Grassman variables, recognizing that $\theta^2=\bar\theta^2=0$ restricts
+the series to two derivatives.
 
-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
+From here the process can be treated as usual, averaging over the couplings and
+replacing bilinear combinations of the fields with their own variables via a
+Hubbard--Stratonovich transformation. Defining the supermatrix
 \begin{equation}
-  \lim_{R\to1}\lim_{\kappa\to1}\log\overline{\mathcal N}(\kappa,0,R)
-  = \frac12N\log(p-1),
+  A(i,j)=\frac1N\{\zeta(i),\zeta^*(i)\}\otimes\{\zeta(j),\zeta^*(j)\}
 \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{figs/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{figs/threshold_2.000.pdf}
-  \includegraphics{figs/threshold_1.325.pdf} \\
-  \includegraphics{figs/threshold_1.125.pdf}
-  \includegraphics{figs/threshold_1.000.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.
-
+the result can be written
 \begin{equation}
-  H_p=\frac1{p!}\sum_{i_1\cdots i_p}J_{i_1\cdots i_p}z_{i_1}\cdots z_{i_p}
+  \Sigma=
+  \int d1\,d2\,\operatorname{Tr}\left(
+      \frac1{16}\left[
+        \matrix{1&1\cr1&1}
+      \right]A^{(p)}(1,2)+\frac p4\left[
+        \matrix{\epsilon&0\cr0&\epsilon^*}
+      \right](I-A(1,1))\delta(1,2)
+    \right)
+    +\frac12\log\det A
 \end{equation}
-
-\begin{figure}
-  \begin{gnuplot}[terminal=epslatex, terminaloptions={size 8.65cm,5.35cm}]
-    set parametric
-    set hidden3d
-    set isosamples 100,25
-    set samples 100,100
-    unset key
-    set dummy u,r
-    set urange [-pi:pi]
-    set vrange [1:1.5]
-    set cbrange [0:2]
-    set xyplane 0
-
-    set xlabel '$\operatorname{Re}\epsilon$'
-    set ylabel '$\operatorname{Im}\epsilon$'
-    set zlabel '$r$'
-    set cblabel '$\frac\epsilon{\epsilon_{\mathrm{th}}}$'
-
-    p = 4
-    set palette defined (0 "blue", 0.99 "blue", 1.0 "white", 1.01 "red", 2 "red")
-    set pm3d depthorder border linewidth 0.5
-
-    s(r) = sqrt(0.75 * log(9 * r**4 / (1 + r**2 + r**4)) / (8 * r**4 - r**2 - 1))
-    x(u, r) = cos(u) * s(r) * sqrt(1 + 5 * r**2 + 5 * r**4 + r**6)
-    y(u, r) = sin(u) * s(r) * sqrt((r**2 - 1)**3)
-    thres(u, r) = ((x(u,r) / (r**(p - 2) + 1))**2 + (y(u,r) / (r**(p - 2) - 1))**2) / ((p - 1) / (2 * p * r**(p - 2)))
-
-    splot "++" using (x(u, r)):(y(u, r)):2:(thres(u, r)) with pm3d lc palette
-  \end{gnuplot}
-  \caption{
-    The surface of extant states for the 4-spin model, that is, those for which
-    the complexity is zero.
-  }
-\end{figure}
+where the exponent in parentheses denotes element-wise exponentiation, and
+\begin{equation}
+  \delta(i,j)=(\theta(i)-\theta(j))(\bar\theta(i)-\bar\theta(j))
+\end{equation}
+is the super Dirac-$\delta$, and the determinant is a superdeterminant. This leads to the condition for a saddle point of
+\begin{equation}
+  0
+  =\frac{\partial\Sigma}{\partial A(1,2)}
+  =\frac p{16}A^{(p-1)}(1,2)-\frac p4\left[
+    \matrix{\epsilon&0\cr0&\epsilon^*}
+  \right]\delta(1,2)
+  +\frac12A^{-1}(1,2)
+\end{equation}
+where $\odot$ denotes element-wise multiplication and the inverse superfield is defined by
+\begin{equation}
+  I\delta(1,2)=\int d3\,A^{-1}(1,3)A(3,2)
+\end{equation}
+Making such a transformation, we arrive at the saddle point equations
+\begin{eqnarray}
+  0
+  &=\int d3\,\frac{\partial\Sigma}{\partial A(1,3)}A(3,2) \\
+  &=\frac p{16}\int d3A^{(p-1)}(1,3)A(3,2)-\frac p4\left[
+    \matrix{\epsilon&0\cr0&\epsilon^*}
+  \right]A(1,2)+\frac12I\delta(1,2)
+\end{eqnarray}
+When expanded, the supermatrix $A$ contains nine independent bilinear
+combinations of the original variables: $z^\dagger z$, $\hat z^T z$, $\hat
+z^\dagger z$, $\hat z^T\hat z$, $\hat z^\dagger\hat z$, $\eta^\dagger\eta$,
+$\gamma^\dagger\gamma$, $\eta^\dagger\gamma$, and $\eta^T\gamma$. The saddle
+point equations can be used to eliminate all but one of these, the `radius'
+like term $z^\dagger z$. When combined with the constraint, this term can be
+related directly to the magnitude of the imaginary part of $z$, since
+$z^\dagger z=x^Tx+y^Ty=N+2y^Ty$. We therefore define $Y=\frac1N(z^\dagger
+z-N)$, the specific measure of the distance into the complex plane from the
+real sphere. The complexity can then be written
 
 \subsection{Pure \textit{p}-spin: where are my neighbors?}
 
@@ -1392,7 +1290,7 @@ gives a sense of whether many Stokes lines should be expected, and when.
 
 To determine this, we perform the same Kac--Rice produce as in the previous
 section, but now with two probe points, or replicas of the system. The number of
-critical points with given energies $\epsilon_1$ and $\epsilon_2$ are
+stationary points with given energies $\epsilon_1$ and $\epsilon_2$ are
 \begin{equation}
   \mathcal N(\epsilon_1,\epsilon_2)
   =\int dx\,dz\,dz^*\,\delta(\partial H(x))\,\delta(\operatorname{Re}\partial H(z))\delta(\operatorname{Im}\partial H(z))|\det\operatorname{Hess}H(z)|^2|\det\operatorname{Hess}H(x)|