Lots of writing.

10 files changed, 217 insertions, 112 deletions

diff --git a/figs/ground_complexity.pdf b/figs/ground_complexity.pdf
index fe3a3b0..5bb45ce 100644
--- a/figs/ground_complexity.pdf
+++ b/figs/ground_complexity.pdf
diff --git a/figs/im_complexity.pdf b/figs/im_complexity.pdf
index 8c5bf5b..7266b8c 100644
--- a/figs/im_complexity.pdf
+++ b/figs/im_complexity.pdf
diff --git a/figs/leg_complexity.pdf b/figs/leg_complexity.pdf
index 5b093eb..6fe37f8 100644
--- a/figs/leg_complexity.pdf
+++ b/figs/leg_complexity.pdf
diff --git a/figs/neighbor_energy_limit.pdf b/figs/neighbor_energy_limit.pdf
index ee059ea..34ef5d8 100644
--- a/figs/neighbor_energy_limit.pdf
+++ b/figs/neighbor_energy_limit.pdf
diff --git a/figs/neighbor_geometry.pdf b/figs/neighbor_geometry.pdf
index c17e7b6..5e459d4 100644
--- a/figs/neighbor_geometry.pdf
+++ b/figs/neighbor_geometry.pdf
diff --git a/figs/neighbor_plot.pdf b/figs/neighbor_plot.pdf
index 9cf0a95..f7ed1be 100644
--- a/figs/neighbor_plot.pdf
+++ b/figs/neighbor_plot.pdf
diff --git a/figs/obuchi_3-spin.pdf b/figs/obuchi_3-spin.pdf
index 1ca8ac3..150b748 100644
--- a/figs/obuchi_3-spin.pdf
+++ b/figs/obuchi_3-spin.pdf
diff --git a/figs/re_complexity.pdf b/figs/re_complexity.pdf
index db4ed5a..4efac49 100644
--- a/figs/re_complexity.pdf
+++ b/figs/re_complexity.pdf
diff --git a/stokes.bib b/stokes.bib
index 5664941..02ad0bc 100644
--- a/stokes.bib
+++ b/stokes.bib
@@ -355,6 +355,20 @@
  doi = {10.1103/physrevresearch.3.023064}
 }
 
+@article{Kurchan_1992_Supersymmetry,
+ author = {Kurchan, J.},
+ title = {Supersymmetry in spin glass dynamics},
+ journal = {Journal de Physique I},
+ publisher = {EDP Sciences},
+ year = {1992},
+ month = {7},
+ number = {7},
+ volume = {2},
+ pages = {1333--1352},
+ url = {https://doi.org/10.1051%2Fjp1%3A1992214},
+ doi = {10.1051/jp1:1992214}
+}
+
 @book{Livan_2018_Introduction,
  author = {Livan, Giacomo and Novaes, Marcel and Vivo, Pierpaolo},
  title = {Introduction to Random Matrices},
diff --git a/stokes.tex b/stokes.tex
index 6081426..d16569e 100644
--- a/stokes.tex
+++ b/stokes.tex
@@ -1060,7 +1060,7 @@ determination of the constant shift $\lambda_0$ 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}
+\begin{equation} \label{eq:gap.eigenvalue}
   |\lambda_{\mathrm{gap}}|^2
   =\Gamma_0\frac{(1-|\delta|^2)^2}
   {1+|\delta|^2-2|\delta|\cos(\arg\delta+2\arg\lambda_0)}
@@ -1195,44 +1195,57 @@ by the symmetry $z\to-z$, as seen in Fig.~\ref{fig:3d.thimbles}.
 
 Since the 2-spin model with real couplings does not have any stationary points
 in the complex plane, analytic continuation can be made without any fear of
-running into Stokes points. Starting from real, large $\beta$, making an infinitesimal phase rotation into the complex plane results in a decomposition into thimbles where that of each stationary point is necessary, because all stationary points are real. The curvature of the action at the stationary points lying at $z_i=\delta_{ik}$ in the $j$th direction is given by $\lambda_k-\lambda_j=2(\epsilon_k-\epsilon_k)$. Therefore the generic case of $N$ distinct eigenvalues of the coupling matrix leads to $2N$ stationary points with $N$ distinct energies, two at each index from $0$ to $N-1$. Starting with the expression \eqref{eq:real.thimble.partition.function} valid for the partition function contribution from the thimble of a real stationary point, we have
+running into Stokes points. Starting from real, large $\beta$, making an
+infinitesimal phase rotation into the complex plane results in a decomposition
+into thimbles where that of each stationary point is necessary, because all
+stationary points are real. The curvature of the action at the stationary
+points lying at $z_i=\delta_{ik}$ in the $j$th direction is given by
+$\lambda_k-\lambda_j=2(\epsilon_k-\epsilon_k)$. Therefore the generic case of
+$N$ distinct eigenvalues of the coupling matrix leads to $2N$ stationary points
+with $N$ distinct energies, two at each index from $0$ to $D=N-1$. Starting with
+the expression \eqref{eq:real.thimble.partition.function} valid for the
+partition function contribution from the thimble of a real stationary point, we
+have
 \begin{equation}
   \eqalign{
     Z(\beta)
-    &=\int_{S^{N-1}}ds\,e^{-\beta H_2(s)}
-    =\sum_{\sigma\in\Sigma_0}n_\sigma\int_{\mathcal J_\sigma}ds\,e^{-\beta H_2(s)} \\
-    &\simeq\sum_{k=0}^D2i^k\left(\frac{2\pi}\beta\right)^{D/2}e^{-\beta H_2(s_\sigma)}\prod_{\lambda_0>0}\lambda_0^{-1/2} \\
+    &=\int_{S^{N-1}}ds\,e^{-\beta\mathcal S_2(s)}
+    =\sum_{\sigma\in\Sigma_0}n_\sigma\int_{\mathcal J_\sigma}ds\,e^{-\beta\mathcal S_2(s)} \\
+    &\simeq\sum_{\sigma\in\Sigma_0}i^{k_\sigma}\left(\frac{2\pi}\beta\right)^{D/2}e^{-\beta\mathcal S_2(s_\sigma)}|\det\operatorname{Hess}\mathcal S_2(s_\sigma)|^{-\frac12} \\
     &=2\sum_{k=0}^D\exp\left\{
       i\frac\pi2k+\frac D2\log\frac{2\pi}\beta-N\beta\epsilon_k-\frac12\sum_{\ell\neq k}\log2|\epsilon_k-\epsilon_\ell|
     \right\}
   }
 \end{equation}
+where $\epsilon_k$ is the energy of the twin stationary points of index $k$. In the large $N$ limit, we take advantage of the limiting distribution $\rho$ of these energies to write
 \begin{equation} \fl
-  Z(\beta)
-  =2\int d\epsilon\,\rho(\epsilon)\exp\left\{
-    i\frac\pi2k_\epsilon+\frac D2\log\frac{2\pi}\beta-N\beta\epsilon-\frac D2\int d\epsilon'\,\rho(\epsilon')\log2|\epsilon-\epsilon'|
-  \right\}
+  \eqalign{
+    \overline Z(\beta)
+    &=2\int d\epsilon\,\rho(\epsilon)\exp\left\{
+      i\frac\pi2k_\epsilon+\frac D2\log\frac{2\pi}\beta-N\beta\epsilon-\frac D2\int d\epsilon'\,\rho(\epsilon')\log2|\epsilon-\epsilon'|
+    \right\} \\
+    &=2\int d\epsilon\,\rho(\epsilon)e^{Nf(\epsilon)}
+  }
 \end{equation}
-
 Since the $J$ of the 2-spin model is a symmetric real matrix with variance
 $1/N$, its eigenvalues are distributed by a semicircle distribution of radius 2,
 and therefore the energies $\epsilon$ are distributed by a semicircle
-distribution of radius $\epsilon_{\mathrm{th}}=1$, with
+distribution of radius one, with
 \begin{equation}
-  \rho(\epsilon\mid\epsilon_{\mathrm{th}})=\frac2{\pi\epsilon_{\mathrm{th}}^2}\sqrt{\epsilon_{\mathrm{th}}^2-\epsilon^2}
+  \rho(\epsilon)=\frac2{\pi}\sqrt{1-\epsilon^2}
 \end{equation}
 The index as a function of energy level is given by the cumulative density function
 \begin{equation}
-  k_\epsilon=N\int_{-\infty}^\epsilon d\epsilon'\,\rho(\epsilon')=\frac N\pi\left(
-    \frac{\epsilon}{\epsilon_{\mathrm{th}}^2}\sqrt{\epsilon_{\mathrm{th}}^2-\epsilon^2}+2\tan^{-1}\frac{\epsilon_{\mathrm{th}}+\epsilon}{\sqrt{\epsilon_{\mathrm{th}}^2-\epsilon^2}}
+  k_\epsilon=D\int_{-\infty}^\epsilon d\epsilon'\,\rho(\epsilon')=\frac D\pi\left(
+    \epsilon\sqrt{1^2-\epsilon^2}+2\tan^{-1}\frac{1+\epsilon}{\sqrt{1-\epsilon^2}}
   \right)
 \end{equation}
 Finally, the product over the singular values corresponding to descending directions gives
 \begin{equation}
   \frac12\int d\epsilon'\,\rho(\epsilon')\log2|\epsilon-\epsilon'|
-  =-\frac14+\frac12\left(\frac{\epsilon}{\epsilon_{\mathrm{th}}}\right)^2
+  =-\frac14+\frac12\epsilon^2
 \end{equation}
-for $\epsilon<\epsilon_{\mathrm{th}}$. Then
+for $\epsilon^2<1$. This gives the function $f$ in the exponential as
 \begin{equation}
   \operatorname{Re}f=-\epsilon\operatorname{Re}\beta+\frac14-\frac12\epsilon^2
 \end{equation}
@@ -1251,7 +1264,8 @@ and the real part of $f$ comes to a cusp, meaning that the oscillations do not
 interfere in taking the saddle point. Once this line is crossed and the maximum
 enters the bulk of the spectrum, one expects to find cancellations caused by
 the incoherent contributions of thimbles with nearby energies to
-$\epsilon_{\mathrm{max}}$.
+$\epsilon_{\mathrm{max}}$. Therefore, one expects that $\overline Z$ enters a
+phase with no coherent average when $\operatorname{Re}\beta=1$.
 
 On the other hand, there is another point where the thimble sum becomes
 coherent. This is when the oscillation frequency near the maximum energy goes
@@ -1263,44 +1277,85 @@ to zero. This happens for
   =-\operatorname{Im}\beta+\sqrt{1-(\operatorname{Re}\beta)^2}
 \end{equation}
 or for $|\beta|=1$. Here the sum of contributions from thimbles near the
-maximum again becomes coherent. These conditions correspond precisely to those
-found for the density of zeros in the 2-spin model found previously
+maximum again becomes coherent, because the period of oscillations in
+$\epsilon$ diverges at the maximum. These conditions correspond precisely to
+the phase boundaries found for the density of zeros in the 2-spin model found previously using other methods
 \cite{Obuchi_2012_Partition-function, Takahashi_2013_Zeros}.
 
+We've seen that even in the 2-spin model, which is not complex, making a
+thimble decomposition in a theory with many saddles does not necessarily fix
+the sign problem. Instead, it takes a potentially high-dimensional sign problem
+and produces a one-dimensional one, represented by the oscillatory integral
+over $e^{Nf(\epsilon)}$. In some regimes, it can be argued that integral has a
+maximum with a coherent neighborhood, allowing computation to be made. In
+others, oscillations in the phase remain, from the sum over the many thimbles.
+We will find a similar story for the pure $p$-spin models for $p>2$ in the next
+sections, complicated by the additional presence of Stokes points in the
+continuation.
+
 \subsection{Pure \textit{p}-spin: where are the saddles?}
 
 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.
+method and elaborate on some of the results relevant to analytic continuation.
+
+The complexity of the real $p$-spin models has been studied extensively, and is even known rigorously \cite{Auffinger_2012_Random}. If $\mathcal N(\epsilon)$ is the number of stationary points with specific energy $\epsilon$, then the complexity is defined by
+\begin{equation}
+  \Sigma(\epsilon)=\lim_{N\to\infty}\frac1N\log\overline\mathcal N(\epsilon)
+\end{equation}
+a natural measure of how superextensive the average number $\overline\mathcal
+N\sim e^{N\Sigma}$ is. The complexity is also known for saddles of particular
+index, with, e.g., $\Sigma_{k=1}$ measuring the complexity of rank-one saddles
+and $\Sigma_{k=0}$ measuring that of minima. The minimum energy for which
+$\Sigma_{k=0}$ is positive corresponds to the ground state energy of the model,
+because at large $N$ below this the number of minima is expected to be
+exponentially small with $N$. We'll write the ground state energy as $\epsilon_{k=0}$, and the lowest energies at which rank $j$ saddles are found as $\epsilon_{k=j}$, so that, e.g.,
+\begin{equation}
+  0=\Sigma(\epsilon_{k=0})=\Sigma_{k=0}(\epsilon_{k=0}) \qquad
+  0=\Sigma_{k=1}(\epsilon_{k=1})
+\end{equation}
 
 In the real case, the $p$-spin models posses a threshold energy
-$\epsilon_{\mathrm{th}}$, below which there are exponentially many minima
+\begin{equation}
+  |\epsilon_{\mathrm{th}}|^2=\frac{2(p-1)}{p}
+\end{equation}
+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
+more generic form in the complex case, where now the threshold separates
+stationary points that have mostly gapped from mostly ungapped spectra. 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}. 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}
+  =\frac{p(p-1)}{2N}(1+2Y)^{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 $\Gamma_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
+where $Y=\frac1N\|\operatorname{Im}z\|^2$ is a measure of how far the
+stationary point is into the complex configuration space.  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 $\Gamma_0=p(p-1)(1+2Y)^{p-2}/2$, $C_0=p(p-1)/2$, and
+$\lambda_0=-p\epsilon$. This means that the energy at which the gap appears is,
+using \eqref{eq:gap.eigenvalue},
+\begin{equation}
+  |\epsilon_\mathrm{gap}|^2
+  =\frac{p-1}{2p}
+  \frac{[1-(1+2Y)^{2(p-2)}]^2(1+2Y)^{p-2}}
+  {1+(1+2Y)^{2(p-2)}-2(1+2Y)^{p-2}\cos(2\arg\epsilon)}
+\end{equation}
+When $\epsilon$ is real,
+$\lim_{Y\to0}|\epsilon_{\mathrm{gap}}|=|\epsilon_\mathrm{th}|$.
+
+The complexity of stationary points by their energy and location $Y$ 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
     = \int dx\,dy\,\delta(\partial_x\operatorname{Re}\tilde\mathcal S_p)\delta(\partial_y\operatorname{Re}\tilde\mathcal S_p)
@@ -1308,11 +1363,13 @@ The location of stationary points can be determined by the Kac--Rice formula. An
 \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.
+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.
 
-As in \S\ref{sec:stationary.hessian}, these can be bright into a manifestly complex form using Cauchy--Riemann relations. This gives
+As in \S\ref{sec:stationary.hessian}, this expression can be bright into a
+manifestly complex form using Cauchy--Riemann relations. This gives
 \begin{equation}
   \mathcal N
    =\int dz^*dz\,d\hat z^*d\hat z\,d\eta^*d\eta\,d\gamma^*d\gamma\exp\left\{
@@ -1321,8 +1378,13 @@ As in \S\ref{sec:stationary.hessian}, these can be bright into a manifestly comp
     \right)
   \right\}
 \end{equation}
-where $\eta$ and $\gamma$ are Grassmann variables. This can be more
-conveniently studied using the method of superfields. Introducing the
+where $\eta$ and $\gamma$ are $N$-dimensional Grassmann fields. This can be more
+conveniently studied using the method of superfields. For an overview of
+superfields applied to the $p$-spin spherical models, see
+\cite{Kurchan_1992_Supersymmetry}. Our previous work deriving the complexity
+does not use superfields \cite{Kent-Dobias_2021_Complex}, but they will be
+essential for compactly writing the \emph{two} replica complexity in the next
+section, and so we briefly introduce the technique here. Introducing the
 one-component Grassman variables $\theta$ and $\bar\theta$, define the
 superfield
 \begin{equation}
@@ -1350,11 +1412,11 @@ Hubbard--Stratonovich transformation. Defining the supermatrix
       \phi(1)^\dagger\phi(2)&\phi(1)^\dagger\phi(2)^*
   }\right]
 \end{equation}
-the result can be written, neglecting constant factors,
+the result can be written, neglecting constant factors, as an integral over $Q$ like
 \begin{equation}
   \overline\mathcal N\simeq\int dQ\,e^{NS_\mathrm{eff}(Q)}
 \end{equation}
-for an effective action functional of the supermatrix $Q$
+where the effective action functional $S_\mathrm{eff}$ of the supermatrix $Q$ is
 \begin{equation} \fl
   \eqalign{
     S_{\mathrm{eff}}&=
@@ -1368,11 +1430,11 @@ for an effective action functional of the supermatrix $Q$
                                 &\hspace{28em}+\frac12\log\det Q
     }
 \end{equation}
-where the exponent in parentheses denotes element-wise exponentiation, and
+The exponent in parentheses denotes element-wise exponentiation, and
 \begin{equation}
   \delta(1,2)=(\bar\theta(1)-\bar\theta(2))(\theta(1)-\theta(2))
 \end{equation}
-is the superspace $\delta$-function, and the determinant is a superdeterminant.
+is the superspace $\delta$-function, and the determinant and trace are a superdeterminant and supertrace, respectively. Algebraically and under calculus they behave nearly like their non-super counterparts.
 This leads to the condition for a saddle point of
 \begin{equation}
   0
@@ -1386,7 +1448,7 @@ where the inverse supermatrix is defined by
 \begin{equation}
   I\delta(1,2)=\int d3\,Q^{-1}(1,3)Q(3,2)
 \end{equation}
-Making such a transformation, we arrive at the saddle point equations
+Convolving both sides by another supermatrix to remove the inverse, we arrive at the saddle point equations
 \begin{equation}
   \eqalign{
     0
@@ -1403,9 +1465,7 @@ $\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 $(\operatorname{Im}s)^2=\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 in terms of $r=z^\dagger z/N$ as
+$z^\dagger z=x^Tx+y^Ty=N+2y^Ty=N(1+2Y)$ for $Y=\|\operatorname{Im}z\|^2/N=y^Ty/N$. The complexity can then be written in terms of $r=z^\dagger z/N=1+2Y$ as
 \begin{equation}
   \Sigma
   =
@@ -1413,7 +1473,7 @@ real sphere. The complexity can then be written in terms of $r=z^\dagger z/N$ as
     \frac{1-r^{-2(p-1)}}{1-r^{-2}}
   \right)
   -\frac{(\operatorname{Re}\epsilon)^2}{R_+^2}-\frac{(\operatorname{Im}\epsilon)^2}{R_-^2}
-  +I_p(\epsilon/\epsilon_\mathrm{th})
+  +I_p(\epsilon/|\epsilon_\mathrm{th}|)
 \end{equation}
 where
 \begin{equation}
@@ -1423,7 +1483,7 @@ where
     1+r^{2(p-2)}\left[p(p-2)(r^2-1)-1\right]
   }
 \end{equation}
-and the function $I_p(u)=0$ if $|\epsilon|^2<|\epsilon_\mathrm{th}|^2$ and
+and the function $I_p(u)=0$ if $|\epsilon|^2<|\epsilon_\mathrm{gap}|^2$ and
 \begin{equation}
   \eqalign{
     I_p(u)
@@ -1455,7 +1515,7 @@ opposite sign of the imaginary part of $u$.
   \caption{
     The complexity of the 3-spin spherical model in the complex plane, as a
     function of pure real and imaginary energy (left and right) and the
-    magnitude $(\operatorname{Im}s)^2/N$ of the distance into the complex
+    magnitude $Y=\|\operatorname{Im}z\|^2/N$ of the distance into the complex
     configuration space. The thick black contour shows the line of zero
     complexity, where stationary points become exponentially rare in $N$.
   } \label{fig:p-spin.complexity}
@@ -1467,7 +1527,10 @@ line shows the contour of zero complexity, where stationary points are no
 longer found at large $N$. As the magnitude of the imaginary part of the spin
 taken greater, more stationary points are found, and at a wider array of
 energies. This is also true in other directions into the complex energy plane,
-where the story is qualitatively the same.
+where the story is qualitatively the same. At any energy, the limit
+$Y\to\infty$ or $r\to\infty$ results in $\Sigma=\log(p-1)$, which saturates the
+Bézout bound on the number of stationary points a polynomial of order $p$ can
+have \cite{Bezout_1779_Theorie}.
 
 \begin{figure}
   \hspace{2pc}
@@ -1475,11 +1538,11 @@ where the story is qualitatively the same.
 
   \caption{
     The complexity of the 3-spin spherical model in the complex plane, as a
-    function of pure real energy and the magnitude $(\operatorname{Im}s)^2/N$
+    function of pure real energy and the magnitude $Y=\|\operatorname{Im}z\|^2/N$
     of the distance into the complex configuration space. The thick black
     contour shows the line of zero complexity, where stationary points become
     exponentially rare in $N$. The shaded region shows where stationary points
-    have a gapped spectrum. The complexity of the 3-spin model on the real
+    have an ungapped spectrum. The complexity of the 3-spin model on the real
     sphere is shown below the horizontal axis; notice that it does not
     correspond with the limiting complexity in the complex configuration space
     below the threshold energy.
@@ -1489,14 +1552,14 @@ where the story is qualitatively the same.
 Something more interesting is revealed if we zoom in on the complexity around
 the ground state, shown in Fig.~\ref{fig:ground.complexity}. Here, the region
 where most stationary points have a gapped hessian is shaded. The line
-separating gapped from ungapped distribution corresponds to the threshold
-energy $\epsilon_\mathrm{th}$ in the limit of $(\operatorname{Im}s)^2\to0$.
-Above the threshold, the limit of the complexity to zero imaginary component
-(or equivalently $r\to1$) also approaches the real complexity, plotted under
-the horizontal axis. However, below the threshold this is no longer the case:
-here the limit of $(\operatorname{Im}s)^2\to0$ of the complexity of complex
-stationary points corresponds to the complexity of \emph{rank one saddles} in
-the real problem, and their complexity becomes zero at $\epsilon_1$, where the
+$\epsilon_\mathrm{gap}$ separating gapped from ungapped distribution corresponds
+to the threshold energy $\epsilon_\mathrm{th}$ in the limit of $Y\to0$.  Above
+the threshold, the limit of the complexity as $Y\to0$ (or
+equivalently $r\to1$) also approaches the real complexity, plotted under the
+horizontal axis. However, below the threshold this is no longer the case: here
+the limit of $Y\to0$ of the complexity of complex
+stationary points corresponds to the complexity $\Sigma_{k=1}$ of \emph{rank one saddles} in
+the real problem, and their complexity becomes zero at $\epsilon_{k=1}$, where the
 complexity of rank one saddles becomes zero \cite{Auffinger_2012_Random}.
 
 There are several interesting features of the complexity. First is this
@@ -1526,16 +1589,23 @@ points. The distribution of these near neighbors in the complex plane therefore
 gives a sense of whether many Stokes lines should be expected, and when.
 
 To determine this, we perform the same Kac--Rice procedure as in the previous
-section, but now with two probe points, or replicas of the system. The number of
-stationary points with given energies $\epsilon_1$ and $\epsilon_2$ are
+section, but now with two probe points, or replicas of the system. The simplify
+things somewhat, we will examine only the case where the second probe is
+complex; the first probe will be on the real sphere. The number of stationary
+points with given energies $\epsilon_1\in\mathbb R$ and $\epsilon_2\in\mathbb C$ are, in the superfield formulation,
 \begin{equation}
-  \mathcal N
+  \mathcal N^{(2)}
   =\int d\phi_1\,d\phi_2^*\,d\phi_2\,\exp\left\{
     \int d1 \left[
       \tilde\mathcal S_p(\phi_1(1))+\operatorname{Re}\tilde\mathcal S_p(\phi_2(1))
     \right]
   \right\}
 \end{equation}
+and we expect to find a two-spin complexity counting pairs of the form
+\begin{equation}
+  \Sigma^{(2)}=\lim_{N\to\infty}\frac1N\log\mathcal N^{(2)}
+\end{equation}
+which depends on the two energies and on mutual geometric invariants of the two probe points. The calculation follows exactly as before, but with an additional field. The average over $J$ is taken, and the supermatrix
 \begin{equation}
   Q(1,2)
   =\left[
@@ -1546,6 +1616,11 @@ stationary points with given energies $\epsilon_1$ and $\epsilon_2$ are
     }
   \right]
 \end{equation}
+is inserted with a Hubbard--Stratonovich transformation. The average number of pairs can then be written in the form
+\begin{equation}
+  \overline{\mathcal N^{(2)}}\propto\int dQ\,e^{NS_\mathrm{eff}[Q]}
+\end{equation}
+for the effective action
 \begin{equation}
   \eqalign{
     S_\mathrm{eff}
@@ -1560,6 +1635,7 @@ stationary points with given energies $\epsilon_1$ and $\epsilon_2$ are
     \right\}+\frac12\det Q
   }
 \end{equation}
+Differentiating this with respect to $Q$, one finds the saddle point equations
 \begin{equation}
   \eqalign{
     0=\frac{\partial S_\mathrm{eff}}{\partial Q(1,2)}
@@ -1570,7 +1646,7 @@ stationary points with given energies $\epsilon_1$ and $\epsilon_2$ are
     +\frac12Q^{-1}(1,2)
   }
 \end{equation}
-where $\odot$ denotes element-wise multiplication.
+where $\odot$ denotes element-wise multiplication. These are simplified by convolution to remove the superinverse, finally giving
 \begin{equation}
   \eqalign{
     0
@@ -1601,38 +1677,40 @@ with the same invariants, e.g., energy and radius. Making this ansatz, the
 equations can be solved for the remaining 5 bilinear products, eliminating all
 the fermionic fields.
 
-This leaves two bilinear products: $z_2^\dagger z_2$ and $z_2^\dagger x_1$, or one real and one complex number. The first is the radius
-of the complex saddle, while the other is a generalization of the overlap in
-the real case. For us, it will be more convenient to work in terms of the
-difference $\Delta z=z-x$ and the constants which characterize it, which are
-$\Delta=\Delta z^\dagger\Delta z=\|\Delta z\|^2$ and $\delta=\frac{\Delta z^T\Delta z}{|\Delta z^\dagger\Delta z|}$. Once again
-we have one real (and strictly positive) variable $\Delta$ and one complex
-variable $\delta$.
-
-Though the value of $\delta$ is bounded  by $|\delta|\leq1$ as a result of the
-inequality $\Delta z^T\Delta z\leq\|\Delta z\|^2$, in reality this bound is not
-the relevant one, because we are confined on the manifold $N=z^2$. The relevant
+This leaves two bilinear products: $z_2^\dagger z_2$ and $z_2^\dagger z_1$, or one real and one complex number. The first is the radius
+of the complex saddle, while the other is a complex generalization of the overlap. For us, it will be more convenient to work in terms of the
+difference $\Delta z=z_2-z_1$ and the constants which characterize it, which
+are $\Delta=\Delta z^\dagger\Delta z/N=\|\Delta z\|^2/N$ and
+$\gamma=\frac{\Delta z^T\Delta z}{\|\Delta z\|}$. Once again we have one real
+(and strictly positive) variable $\Delta$ and one complex variable $\gamma$.
+
+Though the value of $\gamma$ is bounded  by $|\gamma|\leq1$ as a result of the
+inequality $|\Delta z^T\Delta z|\leq\|\Delta z\|^2$, in reality this bound is not
+the relevant one, because we are confined on the manifold $N=z^Tz$. The relevant
 bound is most easily established by returning to a $2N$-dimensional real
 problem, with $x=x_1$ and $z=x_2+iy_2$. The constraint gives $x_2^Ty_2=0$,
 $x_1^Tx_1=1$, and $x_2^Tx_2=1+y_2^Ty_2$. Then, by their definitions,
 \begin{equation}
   \Delta=1+x_2^Tx_2+y_2^Ty_2-2x_1^Tx_2=2(1+y_2^Ty_2-x_1^Tx_2)
 \end{equation}
+Define $\theta_{xx}$ as the angle between $x_1$ and $x_2$. Then $x_1^Tx_2=\|x_1\|\|x_2\|\cos\theta_{xx}=\sqrt{1-\|y_2\|}\cos\theta_{xx}$, and
 \begin{equation}
-  \Delta=2(1+|y_2|^2-\sqrt{1-|y_2|^2}\cos\theta_{xx})
+  \Delta=2(1+\|y_2\|^2-\sqrt{1-\|y_2\|^2}\cos\theta_{xx})
 \end{equation}
+The definite of $\gamma$ likewise gives
 \begin{equation}
   \eqalign{
-    \delta\Delta&=2-2x_1^Tx_2-2ix_1^Ty_2=2(1-|x_2|\cos\theta_{xx}-i|y_2|\cos\theta_{xy}) \\
+    \gamma\Delta&=2-2x_1^Tx_2-2ix_1^Ty_2=2(1-|x_2|\cos\theta_{xx}-i|y_2|\cos\theta_{xy}) \\
           &=2(1-\sqrt{1-|y_2|^2}\cos\theta_{xx}-i|y_2|\cos\theta_{xy})
   }
 \end{equation}
+where $\theta_{xy}$ is the angle between $x_1$ and $y_2$.
 There is also an inequality between the angles $\theta_{xx}$ and $\theta_{xy}$
 between $x_1$ and $x_2$ and $y_2$, respectively, which takes that form
 $\cos^2\theta_{xy}+\cos^2\theta_{xx}\leq1$. This results from the fact that
 $x_2$ and $y_2$ are orthogonal, a result of the constraint.  These equations
-along with the inequality produce the required bound on $|\delta|$ as a
-function of $\Delta$ and $\arg\delta$.
+along with the inequality produce the required bound on $|\gamma|$ as a
+function of $\Delta$ and $\arg\gamma$, which is plotted in Fig.~\ref{fig:bound}.
 
 \begin{figure}
   \hspace{5pc}
@@ -1641,13 +1719,14 @@ function of $\Delta$ and $\arg\delta$.
   \includegraphics{figs/example_bound.pdf}
 
   \caption{
-    The line bounding $\delta$ in the complex plane as a function of
+    \textbf{Left:} The line bounding $\gamma$ in the complex plane as a function of
     $\Delta=0,1,2,\ldots,6$ (outer to inner). Notice that for $\Delta\leq4$,
-    $|\delta|=1$ is saturated for positive real $\delta$, but is not for
+    $|\gamma|=1$ is saturated for positive real $\gamma$, but is not for
     $\Delta>4$, and $\Delta=4$ has a cusp in the boundary. This is due to
     $\Delta=4$ corresponding to the maximum distance between any two points on
-    the real sphere.
-  }
+    the real sphere. \textbf{Right:} The two-spin complexity for $\Delta=4$ and
+    some energy $\epsilon_1=\epsilon_2$. It approaches $-\infty$ at the boundary.
+  } \label{fig:bound}
 \end{figure}
 
 A lot of information is contained in the full two-replica complexity, but we
@@ -1661,18 +1740,25 @@ adjacency: two points will \emph{not} share a Stokes line if a third intervenes
 with its thimble between them. We reason that the number of `adjacent'
 stationary points of a given stationary point for a generic function in $D$
 complex dimensions scales linearly with $D$. Therefore, if the collection of
-nearest neighbors has a finite complexity, e.g., scales \emph{exponentially}
+nearest neighbors has a nonzero complexity, e.g., scales \emph{exponentially}
 with $D$, crowding around the stationary point in question, then these might be
 expected to overwhelm the possible adjacencies, and so doing simplify the
 problem of determining the properties of the true adjacencies. Until the
 nonlinear flow equations are solved with dynamical mean field theory as has
 been done for instantons \cite{Ros_2021_Dynamical}, this is the best heuristic.
 
-First, we find that for all displacements $\Delta$ and real energies $\epsilon_1$, the maximum complexity is found for some real values of $\epsilon_2$ and $\delta$. Therefore we can restrict our study of the most common neighbors to this. Note that the real part of $\delta$ has a very geometric interpretation in terms of the properties of the neighbors: if a stationary point sits in the complex configuration space near another, $\operatorname{Re}\delta$ can be related to the angle $\varphi$ made between the vector separating these two points and the real configuration space as
+For all displacements $\Delta$ and real energies
+$\epsilon_1$, the maximum complexity is found for some real values of
+$\epsilon_2$ and $\gamma$. Therefore we can restrict our study of the most
+common neighbors to this. Note that the real part of $\gamma$ has a geometric
+interpretation in terms of the properties of the neighbors: if a stationary
+point sits in the complex configuration space near another,
+$\operatorname{Re}\gamma$ can be related to the angle $\varphi$ made between
+the vector separating these two points and the real configuration space as
 \begin{equation}
-  \varphi=\arctan\sqrt{\frac{1+\operatorname{Re}\delta}{1-\operatorname{Re}\delta}}
+  \varphi=\arctan\sqrt{\frac{1+\operatorname{Re}\gamma}{1-\operatorname{Re}\gamma}}
 \end{equation}
-Having concluded that the most populous neighbors are confined to real $\delta$, we will make use of this angle instead of $\delta$, which has a more direct geometric interpretation.
+Having concluded that the most populous neighbors are confined to real $\gamma$, we will make use of this angle instead of $\gamma$, which has a more direct geometric interpretation.
 
 \begin{figure}
   \hspace{5pc}
@@ -1684,29 +1770,36 @@ Having concluded that the most populous neighbors are confined to real $\delta$,
   }
 \end{figure}
 
-First, we examine the important of the threshold.
+First, we examine the importance of the threshold.
 Fig.~\ref{fig:neighbor.complexity.passing.threshold} shows the two-replica
 complexity evaluated at $\Delta=2^{-4}$ and equal energy
 $\epsilon_2=\epsilon_1$ as a function of $\varphi$ for several $\epsilon_1$ as
 the threshold is passed. The curves are rescaled by the complexity
-$\Sigma_2(\epsilon_1)$ of index 2 saddles in the real problem, which is what is
+$\Sigma_{k\geq2}(\epsilon_1)$ of index 2 and greater saddles in the real problem, which is what is
 approached in the limit as $\Delta$ to zero. Below the threshold, the
 distribution of nearby saddles with the same energy by angle is broad and
 peaked around $\varphi=45^\circ$, while above the threshold it is peaked
-strongly near the minimum allowed $\varphi$. At the threshold, the function
+strongly near the maximum $\varphi$ allowed by the bound. At the threshold, the function
 becomes extremely flat.
 
 \begin{figure}
   \includegraphics{figs/neighbor_thres.pdf}
   \caption{
-    The scaled two-replica complexity $\Upsilon$ as a function of angle
+    The scaled two-replica complexity $\Sigma^{(2)}$ as a function of angle
     $\varphi$ with $\epsilon_2=\epsilon_1$, $\Delta=2^{-7}$, and various
     $\epsilon_1$. At the threshold, the function undergoes a geometric
     transition and becomes sharper with decreasing $\Delta$.
   } \label{fig:neighbor.complexity.passing.threshold}
 \end{figure}
 
-One can examine the scaling of these curves as $\Delta$ goes to zero. Both above and below the threshold, one finds a quickly-converging limit of $(\Sigma(\epsilon_1,\epsilon_1,\varphi,\Delta)/\Sigma_2(\epsilon_1)-1)/\Delta$. Above the threshold, these curves converge to a function whose peak is always precisely at $45^\circ$, while below they converge to a function with a peak that grows linearly with $\Delta^{-1}$. At the threshold, the scaling is different, and the function approaches a flat function extremely rapidly, as $\Delta^3$.
+One can examine the scaling of these curves as $\Delta$ goes to zero. Both
+above and below the threshold, one finds a quickly-converging limit of
+$(\Sigma^{(2)}/\Sigma_{k\geq2}-1)/\Delta$.
+Above the threshold, these curves converge to a function whose peak is always
+precisely at $45^\circ$, while below they converge to a function with a peak
+that grows linearly with $\Delta^{-1}$. At the threshold, the scaling is
+different, and the function approaches a flat function extremely rapidly, as
+$\Delta^3$.
 
 \begin{figure}
   \includegraphics{figs/neighbor_limit_thres_above.pdf}
@@ -1717,7 +1810,7 @@ One can examine the scaling of these curves as $\Delta$ goes to zero. Both above
   \hfill
   \includegraphics{figs/neighbor_limit_thres_legend.pdf}
   \caption{
-    The scaled two-replica complexity $\Upsilon$ as a function of angle
+    The scaled two-replica complexity $\Sigma^{(2)}$ as a function of angle
     $\varphi$ for various $\Delta$, $\epsilon_2=\epsilon_1$, and \textbf{Left:}
     $\epsilon_1=\epsilon_\mathrm{th}+0.001$ \textbf{Center:}
     $\epsilon_1=\epsilon_\mathrm{th}$ \textbf{Right:}
@@ -1745,13 +1838,13 @@ positive complexity at zero distance, but the maximum is never at
 $\epsilon_2=\epsilon_1$, but rather at a small distance $\Delta\epsilon$ that
 decreases with decreasing $\Delta$ like $\Delta^2$. When the complexity is
 maximized in both parameters, one finds that, in the limit as $\Delta\to0$, the
-peak is at $90^\circ$ but has a height equal to $\Sigma_1$, the complexity of
+peak is at $90^\circ$ but has a height equal to $\Sigma_{k=1}$, the complexity of
 rank-1 saddles.
 
 \begin{figure}
   \includegraphics{figs/neighbor_energy_limit.pdf}
   \caption{
-    The two-replica complexity $\Upsilon$ scaled by $\Sigma_1$ as a function of
+    The two-replica complexity $\Sigma^{(2)}$ scaled by $\Sigma_{k=1}$ as a function of
     angle $\varphi$ for various $\Delta$ at $\epsilon_1=\epsilon_{k=2}$, the
     point of zero complexity for rank-two saddles in the real problem.
     \textbf{Solid lines:} The complexity evaluated at the value of $\epsilon_2$
@@ -1763,13 +1856,12 @@ rank-1 saddles.
 
 Below $\epsilon_{k=1}$, where the rank-1 saddle complexity vanishes, the complexity of stationary points of any type at zero distance is negative. To find what the nearest population looks like, one must find the minimum $\Delta$ at which the complexity is nonnegative, or
 \begin{equation}
-  \Delta_\textrm{min}=\operatorname{argmin}_\Delta\left(0\leq\max_{\epsilon_2, \varphi}\Upsilon(\epsilon_1,\epsilon_2,\Delta,\varphi)\right)
+  \Delta_\textrm{min}=\operatorname{argmin}_\Delta\left(0\leq\max_{\epsilon_2, \varphi}\Sigma^{(2)}(\epsilon_1,\epsilon_2,\Delta,\varphi)\right)
 \end{equation}
 The result in $\Delta_\textrm{min}$ and the corresponding $\varphi$ that
 produces it is plotted in Fig.~\ref{fig:nearest.properties}. As the energy is
 brought below $\epsilon_{k=1}$, $\epsilon_2-\epsilon_1\propto
--|\epsilon_1-\epsilon_{k=1}|^2$, $\varphi-90^\circ\propto-|\epsilon_1-\mathcal
-E_1|^{1/2}$, and $\Delta_\textrm{min}\propto|\epsilon_1-\epsilon_{k=1}|$. The
+-|\epsilon_1-\epsilon_{k=1}|^2$, $\varphi-90^\circ\propto-|\epsilon_1-\epsilon_{k=1}|^{1/2}$, and $\Delta_\textrm{min}\propto|\epsilon_1-\epsilon_{k=1}|$. The
 fact that the population of nearest neighbors has a energy lower than the
 stationary point gives some hope for the success of continuation involving
 these points: since Stokes points only lead to a change in weight when they
@@ -1793,11 +1885,10 @@ change of weight until the phase of $\beta$ has rotated almost $180^\circ$.
   \includegraphics{figs/neighbor_plot.pdf}
 
   \caption{
-    The properties of the nearest neighbor saddles as a function of energy
-    $\epsilon$. Above the threshold energy $\mathcal E_\mathrm{th}$, stationary
+    The properties of the nearest neighbor saddles in the 3-spin model as a function of energy
+    $\epsilon$. Above the threshold energy $\epsilon_\mathrm{th}$, stationary
     points are found at arbitrarily close distance and at all angles $\varphi$
-    in the complex plane. Below $\mathcal E_\mathrm{th}$ but above $\mathcal
-    E_2$, stationary points are still found at arbitrarily close distance and
+    in the complex plane. Below $\epsilon_\mathrm{th}$ but above $\epsilon_{k=2}$, stationary points are still found at arbitrarily close distance and
     all angles, but there are exponentially more found at $90^\circ$ than at
     any other angle. Below $\epsilon_{k=2}$ but above $\epsilon_{k=1}$, stationary
     points are found at arbitrarily close distance but only at $90^\circ$.
@@ -1833,12 +1924,12 @@ where $\dot z$ is the flow at $z$ given by \eqref{eq:flow}, this cost is given b
 \end{equation}
 $\mathcal C$ has minimum of zero, which is reached only by functions $z(t)$
 whose tangent is everywhere parallel to the direction $\dot z$ of the dynamics.
-Therefore, functions that minimize $\mathcal C$ are time-reparameterized Stokes
+Therefore, functions that satisfy $\mathcal C=0$ are time-reparameterized Stokes
 lines.
 
 We explicitly computed the gradient and Hessian of $\mathcal C$ with respect to
 the parameter vectors $g$. Stokes lines are found or not between points by
-using the Levenberg--Marquardt algorithm starting from $g^{(i)}=0$ for all $i$,
+using the Levenberg--Marquardt algorithm starting from $g_i=0$ for all $i$,
 and approximating the cost integral by a finite sum. To sample nearby
 stationary points and assess their propensity for Stokes points, we do the
 following. First, a saddle-finding routine based on Newton's method is run on
@@ -1965,13 +2056,13 @@ terms to the typical partition function
 \end{equation}
 where
 \begin{eqnarray} \fl
-  Z_A
+  \overline Z_A
   \simeq\sum_{k=0}^D\int d\epsilon\,\overline\mathcal N(\epsilon,k)
     \left(\frac{2\pi}\beta\right)^{D/2}i^k
     |\det\operatorname{Hess}\mathcal S|^{-\frac12} e^{-\beta N\epsilon}
     =\int d\epsilon\,e^{Nf_A(\epsilon)}
     \\ \fl
-  Z_B
+  \overline Z_B
   \simeq\sum_{k=0}^D\int d\epsilon\,\eta(\epsilon,k)\overline\mathcal N(\epsilon,k)^{1/2}
     \left(\frac{2\pi}\beta\right)^{D/2}i^k
     |\det\operatorname{Hess}\mathcal S(s_\sigma)|^{-\frac12}e^{-\beta N\epsilon} \\
@@ -2007,7 +2098,7 @@ The temperature at which this happens is
   \operatorname{Re}\beta_A&=-\frac12\frac{3p-4}{p-1}\epsilon_{k=0}+\frac12\frac p{p-1}\sqrt{\epsilon_{k=0}^2-\epsilon_\mathrm{th}^2}\\
   \operatorname{Re}\beta_B&=-\epsilon_{k=0}
 \end{eqnarray}
-which for all $p\geq2$ has $\operatorname{Re}\beta_A>\operatorname{Re}\beta_B$.
+which for all $p\geq2$ has $\operatorname{Re}\beta_A\geq\operatorname{Re}\beta_B$.
 Therefore, the emergence of zeros in $Z_A$ does not lead to the emergence of
 zeros in the partition function as a whole, because $Z_B$ still produces a
 coherent result (despite the unknown constant factor $\eta(\epsilon_{k=0})$). It is