\documentclass[]{iopart}

\usepackage[utf8]{inputenc} % why not type "Stokes" with unicode?
\usepackage[T1]{fontenc} % vector fonts
\usepackage[
  colorlinks=true,
  urlcolor=purple,
  citecolor=purple,
  filecolor=purple,
  linkcolor=purple
]{hyperref} % ref and cite links with pretty colors
\usepackage{amsopn, amssymb, graphicx, xcolor} % standard packages
\usepackage[subfolder]{gnuplottex} % need to compile separately for APS

\begin{document}

\newcommand\eqref[1]{\eref{#1}}
\makeatletter
\renewcommand\tableofcontents{\@starttoc{toc}}
\makeatother

\title{Analytic continuation over complex landscapes}

\author{Jaron Kent-Dobias and Jorge Kurchan}

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

\begin{abstract}
  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 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 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
  those of ``many step replica-symmetry broken'' lie at the threshold and
  Stokes points immediately proliferate.
  A new matrix ensemble is found, playing the role that  GUE plays for real landscapes in determining
  the topological nature of saddles.
\end{abstract}

\maketitle

\tableofcontents

\section{Preamble}

Complex landscapes are basically functions of many variables having many minima and, inevitably,
many saddles of all `indices' (their number of unstable directions). Optimization problems 
require us to find the deepest minima, often a difficult task.
For example, particles with a repulsive mutual potential enclosed in a box will have many stable
configurations, and we are asked to find the one with lowest energy.

An aim of complexity theory is to be able to classify these landscapes in families having
common properties. Two simplifications make the task potentially tractable. The first is to consider
the limit of many variables. In the example of the particles, the limit of many particles (i.e. the thermodynamic limit)
may be expected to bring about simplifications.
The second simplification is of more technical nature:  we consider  functions that contain some random element to them, and we
study the average of an ensemble. The paradigm of this is the spin-glass, where the interactions are random, and we are asked to find the
ground state energy {\em on average over randomness}.

Spin glass theory gave a surprise: random   landscapes come in two kinds:  those that
have a `threshold level' of energy, below which there are many minima but almost no saddles, 
separated by high barriers, and those that have all sorts of saddles all the way down to the lowest 
energy levels, and local minima are separated by relatively small barriers.
The latter are still complex, but good solutions are easier to find.
This classification is closely related to the structure of their Replica Trick solutions.
Armed with this solvable (random) example, it was easy to find non-random examples 
that behave, at least approximately,  in these two ways (e.g. sphere packings and the Travelling salesman Problem,
belong to first and second classes, respectively).

What about systems whose variables are not real, but rather complex? 
Recalling the Cauchy-Riemann conditions, we immediately see a difficulty: if our cost is, say, the real part of a function of
$N$ complex variables, in terms of the corresponding $2N$ real variables it has only saddles of index $N$!
Even worse: often not all saddles are equally interesting, so simply finding the lowest is not usually what we
need to do (more about this below).
As it turns out, there is a set of interesting questions to ask, as we describe below. For each saddle, there is a `thimble' 
spanned by the
lines along which the cost function decreases. The way in which these thimbles fill the complex space is crucial for many
problems of analytic continuation, and is thus what we need to study. The central role played by saddles
in a real landscape, the `barriers', is now played by the Stokes lines, where thimbles exchange their properties. 
Perhaps not surprisingly, the two classes of real landscapes described above retain their role in the complex case, but now  
the distinction is that while in the first class the Stokes lines are rare, in the second class they proliferate.

In  this paper we shall start from a many-variable integral of a real function, and deform it in the many variable complex space.
The landscape one faces is the full one in this space, and we shall see that this is an example where the proliferation -- or lack of it -- 
of Stokes lines is the interesting quantity  in this context.


\section{Introduction}

Analytic continuation of physical theories is sometimes useful. Some theories
have a well-motivated hamiltonian or action that nevertheless results in a
divergent partition function, and can only be properly defined by continuation
from a parameter regime where everything is well-defined \cite{}. Others result
in oscillatory phase space measures that spoil the use of Monte Carlo or saddle
point techniques, but can be treated in a regime where the measure does not
oscillated and the results continued to the desired model \cite{}.

In any case, the nicest modern technique (which we will describe in some detail
later) consists of deforming the phase space integral into a complex phase
space and then breaking it into pieces associated with stationary points of the
action. Each of these pieces, known as \emph{thimbles}, has wonderful
properties that guarantee convergence and prevent oscillations. Once such a
decomposition is made, analytic continuation is mostly easy, save for instances
where the thimbles interact, which must be accounted for.

When your action has a manageable set of stationary points, this process is
usually tractable. However, many actions of interest are complex, having many
stationary points with no simple symmetry relating them, far too many to
individually monitor. Besides appearing in classical descriptions of structural
and spin glasses, complex landscapes are recently become important objects of
study in the computer science of machine learning, the condensed matter theory
of strange metals, and the high energy physics of black holes. What becomes of
analytic continuation under these conditions?

\section{Thimble integration and analytic continuation}

\subsection{Decomposition of the partition function into thimbles}

Consider an action $\mathcal S$ defined on the (real) phase space $\Omega$. A
typical calculation stems from the partition function
\begin{equation} \label{eq:partition.function}
  Z(\beta)=\int_\Omega ds\,e^{-\beta\mathcal S(s)}
\end{equation}
We've defined $Z$ in a way that strongly suggests application in statistical
mechanics, but everything here is general: the action can be complex- or even
imaginary-valued, and $\Omega$ could be infinite-dimensional. In typical
contexts, $\Omega$ will be the euclidean real space $\mathbb R^N$ or some
subspace of this like the sphere $S^{N-1}$ (as in the $p$-spin spherical models
we will study later). We will consider only the analytic continuation of the
parameter $\beta$, but any other would work equally well, e.g., of some
parameter inside the action. The action will have some stationary points, e.g.,
minima, maxima, saddles, and the set of those points in $\Omega$ we will call
$\Sigma_0$, the set of real stationary points.

In order to analytically continue \eref{eq:partition.function} by the method we
will describe, $\mathcal S$ must have an extension to a holomorphic function on
a larger complex phase space $\tilde\Omega$ containing $\Omega$. In many cases
this is accomplished by simply noticing that the action is some sum or product
of holomorphic functions, e.g., polynomials. For $\mathbb R^N$ the complex
phase space $\tilde\Omega$ will be $\mathbb C^N$, while for the sphere
$S^{N-1}$ it takes a little more effort. $S^{N-1}$ can be defined by all points
$x\in\mathbb R^N$ such that $x^Tx=1$. A complex extension of the sphere is made
by extending this constraint: all points $z\in \mathbb C^N$ such that $z^Tz=1$.
Both cases are complex manifolds and moreover Kähler manifolds, since they
are defined by holomorphic constraints, and therefore admit a hermitian
metric and a symplectic structure. In the extended complex phase space, the
action potentially has more stationary points. We'll call $\Sigma$ the set of
\emph{all} stationary points of the action, which naturally contains the set of
\emph{real} stationary points $\Sigma_0$.

\begin{figure}
  \includegraphics{figs/action.pdf}\hfill
  \includegraphics{figs/stationaryPoints.pdf}

  \caption{
    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
    the example figures in this section, we will mostly use the single angular
    variable $\theta$ defined by $s_1=\sqrt2\cos\theta$,
    $s_2=\sqrt2\sin\theta$, which parameterizes the unit circle and its complex
    extension, as $\cos^2\theta+\sin^2\theta=1$ is true even for complex
    $\theta$. \textbf{Right:} The stationary points of $\mathcal S$ in the
    complex-$\theta$ plane. In this example,
    $\Sigma=\{\blacklozenge,\bigstar,\blacktriangle,\blacktriangledown,\bullet,\blacksquare\}$
    and $\Sigma_0=\{\blacklozenge,\blacktriangledown\}$. Symmetries exist
    between the stationary points both as a result of the conjugation symmetry
    of $\mathcal S$, which produces the vertical reflection, and because in the
    pure 3-spin models $\mathcal S(-s)=-\mathcal S(s)$, which produces the
    simultaneous translation and inversion symmetry.
  } \label{fig:example.action}
\end{figure}

Assuming $\mathcal S$ is holomorphic (and that the phase space $\Omega$ is
orientable, which is usually true) the integral in \eref{eq:partition.function}
can be considered an integral over a contour in the complex phase space $\tilde\Omega$,
or
\begin{equation} \label{eq:contour.partition.function}
  Z(\beta)=\oint_\Omega dz\,e^{-\beta\mathcal S(z)}
\end{equation}
For the moment this translation has only changed some of our symbols from
\eref{eq:partition.function}, but conceptually it is very important: contour
integrals can have their contour freely deformed (under some constraints)
without changing their value. This means that we are free to choose a nicer
contour than our initial phase space $\Omega$.

\begin{figure}
  \includegraphics{figs/hyperbola_1.pdf}\hfill
  \includegraphics{figs/hyperbola_2.pdf}\hfill
  \includegraphics{figs/hyperbola_3.pdf}\\
  \includegraphics{figs/anglepath_1.pdf}\hfill
  \includegraphics{figs/anglepath_2.pdf}\hfill
  \includegraphics{figs/anglepath_3.pdf}

  \caption{
    A schematic picture of the complex phase space for the circular $p$-spin
    model and its standard integration contour. \textbf{Top:} For real variables,
    the model is a circle, and its analytic continuation is a kind of complex
    hyperbola, here shown schematically in three dimensions. \textbf{Bottom:}
    Since the real manifold (the circle) is one-dimensional, the complex
    manifold has one complex dimension, here parameterized by the angle
    $\theta$ on the circle. \textbf{Left:} The integration contour over the real phase
    space of the circular model. \textbf{Center:} Complex analysis implies that the
    contour can be freely deformed without changing the value of the integral.
    \textbf{Right:} A funny deformation of the contour in which pieces have been
    pinched off to infinity. So long as no poles have been crossed, even this
    is legal.
  }
\end{figure}

What contour properties are desirable? Consider the two main motivations cited
in the introduction for performing analytic continuation in the first place: we
want our partition function to be well-defined, e.g., for the phase space
integral to converge, and we want to avoid oscillations in the phase of the
integrand. The first condition, convergence, necessitates that the real part of
the action $\operatorname{Re}\beta\mathcal S$ be bounded from below, and that it
approach infinity in any limiting direction along the contour. The second,
constant phase, necessitates that the imaginary part of the action
$\operatorname{Im}\beta\mathcal S$ be constant.

Remarkably, there is an elegant recipe for accomplishing both these criteria at
once, courtesy of Morse theory. For a more thorough review, see
\cite{Witten_2011_Analytic}. We are going to construct our deformed contour out
of a collection of pieces called \emph{Lefschetz thimbles}, or just thimbles.
There is one thimble $\mathcal J_\sigma$ associated with each of the stationary
points $\sigma\in\Sigma$ of the action, and it is defined by all points that
approach the stationary point $z_\sigma$ under gradient descent on
$\operatorname{Re}\beta\mathcal S$.

Thimbles guarantee convergent integrals by construction: the value of
$\operatorname{Re}\beta\mathcal S$ is bounded from below on the thimble $\mathcal J_\sigma$ by its value
$\operatorname{Re}\beta\mathcal S(z_\sigma)$ at the stationary point,
since all other points on the thimble must descend to reach it. And, as we will
see in a moment, thimbles guarantee constant phase for the integrand as well, a
result of the underlying complex geometry of the problem.

What thimbles are necessary to reproduce our original contour, $\Omega$? The
answer is, we need the minimal set which produces a contour between the same
places. Simply stated, if $\Omega=\mathbb R$ produced a phase space integral
running along the real line from left to right, then our contour must likewise
take one continuously from left to right, perhaps with detours to well-behaved
places at infinity (see Fig.~\ref{fig:thimble.homology}). The less simply stated versions follows.

Let $\tilde\Omega_T$ be the set of all points $z\in\tilde\Omega$ such that
$\operatorname{Re}\beta\mathcal S(z)\geq T$, where we will take $T$ to be a very,
very large number. $\tilde\Omega_T$ is then the parts of the manifold where it
is safe for any contour to end up if it wants its integral to converge, since
these are the places where the real part of the action is very large and the
integrand vanishes exponentially. The relative homology group
$H_N(\tilde\Omega,\tilde\Omega_T)$ describes the homology of cycles which begin
and end in $\Omega_T$, i.e., are well-behaved. Therefore, any well-behaved
cycle must represent an element of $H_N(\tilde\Omega,\tilde\Omega_T)$. In order
for our collection of thimbles to produce the correct contour, the composition
of the thimbles must represent the same element of this relative homology
group.

\begin{figure}
  \includegraphics{figs/thimble_homology.pdf}
  \hfill
  \includegraphics{figs/antithimble_homology.pdf}

  \caption{
    A demonstration of the rules of thimble homology. Both figures depict the
    complex-$\theta$ plane of action $\mathcal S$ featured in
    Fig.~\ref{fig:example.action} with $\arg\beta=0.4$. The black symbols lie
    on the stationary points of the action, and the grey regions depict the
    sets $\tilde\Omega_T$ of well-behaved regions at infinity (here $T=5$).
    \textbf{Left:} Lines show the thimbles of each stationary point. The
    thimbles necessary to recreate the cyclic path from left to right are
    darkly shaded, while those unnecessary for the task are lightly shaded.
    Notice that all thimbles come and go from the well-behaved regions.
    \textbf{Right:} Lines show the antithimbles of each stationary point.
    Notice that those of the stationary points involved in the contour (shaded
    darkly) all intersect the desired contour (the real axis), while those not
    involved do not intersect it.
  } \label{fig:thimble.homology}
\end{figure}

Each thimble represents an element of the relative homology, since each thimble
is a contour on which the real part of the action diverges in any direction.
And, thankfully for us, Morse theory on our complex manifold $\tilde\Omega$
implies that the set of all thimbles produces a basis for this relative
homology group, and therefore any contour can be represented by some
composition of thimbles! There is even a systematic way to determine the
contribution from each thimble: for the stationary point $\sigma\in\Sigma$, let
$\mathcal K_\sigma$ be its \emph{antithimble}, defined by all points brought to
$z_\sigma$ by gradient \emph{ascent} (and representing an element of the
relative homology group $H_N(\tilde\Omega,\tilde\Omega_{-T})$). Then each
thimble $\mathcal J_\sigma$ contributes to the contour with a weight given by
its intersection pairing $n_\sigma=\langle\mathcal C,\mathcal K_\sigma\rangle$.

With these tools in hands, we can finally write the partition function as a sum
over contributions from each thimble, or
\begin{equation} \label{eq:thimble.integral}
  Z(\beta)=\sum_{\sigma\in\Sigma}n_\sigma\oint_{\mathcal J_\sigma}dz\,e^{-\beta\mathcal S(z)}.
\end{equation}
Under analytic continuation, the form of \eref{eq:thimble.integral}
generically persists. When the relative homology of the thimbles is unchanged
by the continuation, the integer weights are likewise unchanged, and one can
therefore use the knowledge of these weights in one regime to compute the
partition function in the other. However, their relative homology can change,
and when this happens the integer weights can be traded between stationary
points. These trades occur when two thimbles intersect, or alternatively when
one stationary point lies in the gradient descent of another. These places are
called \emph{Stokes points}, and the gradient descent trajectories that join
two stationary points are called \emph{Stokes lines}. An example of this
behavior can be seen in Fig.~\ref{fig:1d.stokes}.

\begin{figure}
  \includegraphics{figs/thimble_stokes_1.pdf}\hfill
  \includegraphics{figs/thimble_stokes_2.pdf}\hfill
  \includegraphics{figs/thimble_stokes_3.pdf}

  \caption{
    An example of a Stokes point in the continuation of the phase space
    integral involving the action $\mathcal S$ featured in
    Fig.~\ref{fig:example.action}. \textbf{Left:} $\arg\beta=1.176$. The collection of
    thimbles necessary to progress around from left to right, highlighted in a
    darker color, is the same as it was in Fig.~\ref{fig:thimble.homology}.
    \textbf{Center:} $\arg\beta=1.336$. The thimble $\mathcal J_\blacklozenge$
    intersects the stationary point $\blacktriangle$ and its thimble, leading
    to a situation where the contour is not easily defined using thimbles. This
    is a Stokes point. \textbf{Right:} $\arg\beta=1.496$. The Stokes point has passed,
    and the collection of thimbles necessary to produce the path has changed:
    now $\mathcal J_\blacktriangle$ must be included. Notice that in this
    figure, because of the symmetry of the pure models, the thimble $\mathcal
    J_\blacksquare$ also experiences a Stokes point, but this does not result
    in a change to the path involving that thimble.
  } \label{fig:1d.stokes}
\end{figure}

\begin{figure}
  \includegraphics{figs/thimble_orientation_1.pdf}\hfill
  \includegraphics{figs/thimble_orientation_2.pdf}\hfill
  \includegraphics{figs/thimble_orientation_3.pdf}

  \caption{
    The behavior of thimble contours near $\arg\beta=0$ for real actions. In all
    pictures, green arrows depict a canonical orientation of the thimbles
    relative to the real axis, while purple arrows show the direction of
    integration implied by the orientation. \textbf{Left:} $\arg\beta=-0.1$. To
    progress from left to right, one must follow the thimble from the minimum
    $\blacklozenge$ in the direction implied by its orientation, and then
    follow the thimble from the maximum $\blacktriangledown$ \emph{against} the
    direction implied by its orientation, from top to bottom. Therefore,
    $\mathcal C=\mathcal J_\blacklozenge-\mathcal J_\blacktriangledown$.
    \textbf{Center:} $\arg\beta=0$. Here the thimble of the minimum covers
    almost all of the real axis, reducing the problem to the real phase space
    integral. This is also a Stokes point. \textbf{Right:} $\arg\beta=0.1$. Here, one follows the thimble of
    the minimum from left to right again, but now follows that of the maximum
    in the direction implied by its orientation, from bottom to top. Therefore,
    $\mathcal C=\mathcal J_\blacklozenge+\mathcal J_\blacktriangledown$.
  } \label{fig:thimble.orientation}
\end{figure}

The prevalence (or not) of Stokes points in a given continuation, and whether
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
keeping track of the resulting weights, and analytic continuation is intractable.

\subsection{Gradient flow}

The `dynamics' describing thimbles is defined by gradient descent on the real
part of the action, with a given thimble incorporating all trajectories which
asymptotically flow to its associated stationary point. Since our phase space
is not necessary flat (as for the \emph{spherical} $p$-spin models), we will
have to do a bit of differential geometry to work out their form. Gradient
descent on a complex (Kähler) manifold is given by
\begin{equation} \label{eq:flow.coordinate.free}
  \dot s
  =-\operatorname{grad}\operatorname{Re}\beta\mathcal S
  =-\left(\frac\partial{\partial s^*}\operatorname{Re}\beta\mathcal S\right)^\sharp
  =-\frac{\beta^*}2\frac{\partial\mathcal S^*}{\partial s^*}g^{-1}\frac\partial{\partial s}
\end{equation}
where $g$ is the metric and the holomorphicity of the action was used to set
$\partial\mathcal S/\partial s^*=0$. If the complex phase space is $\mathbb C^N$ and the
metric is diagonal, this means that the flow is proportional to the conjugate
of the gradient, or $\dot s\propto-\beta^*(\partial S/\partial s)^*$.

In the cases we will consider here (namely, that of the spherical models), it
will be more convenient to work in terms of coordinates in a flat embedding
space than in terms of local coordinates in the curved space, e.g., in terms of
$z\in\mathbb C^N$ instead of $s\in S^{N-1}$. Let $z:\tilde\Omega\to\mathbb C^N$
be an embedding of complex phase space into complex euclidean space. The
dynamics in the embedding space is given by
\begin{equation}\label{eq:flow.raw}
  \dot z
  =-\frac{\beta^*}2\frac{\partial\mathcal S^*}{\partial z^*}(Dz)^* g^{-1}(Dz)^T\frac\partial{\partial z}
\end{equation}
where $Dz=\partial z/\partial s$ is the Jacobian of the embedding.
The embedding induces a metric on $\tilde\Omega$ by $g=(Dz)^\dagger Dz$.
Writing $\partial=\partial/\partial z$, this gives
\begin{equation} \label{eq:flow}
  \dot z=-\frac{\beta^*}2(\partial\mathcal S)^\dagger(Dz)^*[(Dz)^\dagger(Dz)]^{-1}(Dz)^T
  =-\frac12(\partial \mathcal S)^\dagger P
\end{equation}
which is nothing but the projection of $(\partial\mathcal S)^*$ into the
tangent space of the manifold, with the projection operator
$P=(Dz)^*[(Dz)^\dagger(Dz)]^{-1}(Dz)^T$. Note that $P$ is hermitian. For the spherical models, where $\tilde\Omega$ is the complex phase spaced defined by all points $z\in\mathbb C^N$ such that $z^Tz=1$, the projection operator is given by
\begin{equation}
  P=I-\frac{zz^\dagger}{|z|^2}
\end{equation}
something that we be worked out in detail in a following section. One can
quickly verify that this operator indeed projects the dynamics onto the
manifold: its tangent at any point $z$ is given by $\partial(z^Tz)=z$, and
$Pz=z-z|z|^2/|z|^2=0$. For any vector $u$ perpendicular to $z$, i.e.,
$z^\dagger u=0$, $Pu=u$.

\begin{figure}
  \includegraphics{figs/thimble_flow.pdf}

  \caption{Example of gradient descent flow on the action $\mathcal S$ featured
    in Fig.~\ref{fig:example.action} in the complex-$\theta$ plane, with
    $\arg\beta=0.4$. Symbols denote the stationary points, while thick blue and
    red lines depict the thimbles and antithimbles, respectively. Streamlines
    of the flow equations are plotted in a color set by their value of
    $\operatorname{Im}\beta\mathcal S$; notice that the color is constant along
    each streamline.
  } \label{fig:flow.example}
\end{figure}

Gradient descent on $\operatorname{Re}\beta\mathcal S$ is equivalent to
Hamiltonian dynamics with the Hamiltonian $\operatorname{Im}\beta\mathcal S$
and conjugate coordinates given by the real and imaginary parts of each complex
coordinate. This is because $(\tilde\Omega, g)$ is Kähler and therefore admits
a symplectic structure, but that the flow conserves
$\operatorname{Im}\beta\mathcal S$ can be shown using \eref{eq:flow} and the
holomorphic property of $\mathcal S$:
\begin{eqnarray}
  \frac d{dt}\operatorname{Im}\beta\mathcal S
  &=\dot z\partial\operatorname{Im}\beta\mathcal S+\dot z^*\partial^*\operatorname{Im}\beta\mathcal S \\
  &=\frac i4\left(
    (\beta\partial \mathcal S)^\dagger P\beta\partial\mathcal S-(\beta\partial\mathcal S)^TP^*(\beta\partial\mathcal S)^*
  \right) \\
  &=\frac{i|\beta|^2}4\left(
    (\partial\mathcal S)^\dagger P\partial\mathcal S-[(\partial\mathcal S)^\dagger P\partial\mathcal S]^*
  \right) \\
  &=\frac{i|\beta|^2}4\left(
    \|\partial\mathcal S\|^2-(\|\partial\mathcal S\|^*)^2
  \right)=0.
\end{eqnarray}
A consequence of this conservation is that the flow in the action takes a
simple form:
\begin{equation}
  \dot{\mathcal S}
  =\dot z\partial\mathcal S
  =-\frac{\beta^*}2(\partial\mathcal S)^\dagger P\partial\mathcal S
  =-\frac{\beta^*}2\|\partial\mathcal S\|^2.
\end{equation}
In the complex-$\mathcal S$ plane, dynamics is occurs along straight lines in
a direction set by the argument of $\beta$.

\subsection{The structure of stationary points}
\label{sec:stationary.hessian}

The shape of each thimble in the vicinity of a stationary point can be
described using an analysis of the hessian of the real part of the action at
the stationary point. Here we'll review some general properties of this
hessian, which because the action is holomorphic has rich structure.

Writing down the hessian using the complex geometry of the previous section
would be quite arduous. Luckily, we are only interested in the hessian at
stationary points, and our manifolds of interest are all constraint surfaces.
These two facts allow us to find the hessian at stationary points using a
simpler technique, that of Lagrange multipliers.

Suppose that our complex manifold $\tilde\Omega$ is defined by all points
$z\in\mathbb C^N$ such that $g(z)=0$ for some holomorphic function $z$. In the
case of the spherical models, $g(z)=\frac12(z^Tz-N)$. Introducing the Lagrange
multiplier $\mu$, we define the constrained action
\begin{equation}
  \tilde\mathcal S(z)=\mathcal S(z)-\mu g(z)
\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} \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
the normal to the constraint surface $\partial g$ from the gradient of the
unconstrained action. It implies that the hessian with respect to the complex
embedding coordinates $z$ at any stationary point is
\begin{equation} \label{eq:complex.hessian}
  \operatorname{Hess}\mathcal S
  =\partial\partial\tilde\mathcal S
  =\partial\partial\mathcal S-\frac{\partial g^T\partial\mathcal S}{\partial g^T\partial g}\partial\partial g
\end{equation}
In practice one must neglect the directions normal to the constraint surface by
projecting them out using $P$ from the previous section, i.e.,
$P\operatorname{Hess}\mathcal SP^T$. For notational simplicity we will not
include this here.

In order to describe the structure of thimbles, one must study the Hessian of
$\operatorname{Re}\beta\mathcal S$.  We first pose the problem problem as one
of $2N$ real variables $x,y\in\mathbb R^N$ with $z=x+iy$, the hessian of the
real part of the action with respect to these real variables is
\begin{equation} \label{eq:real.hessian}
  \operatorname{Hess}_{x,y}\operatorname{Re}\beta\mathcal S
  =\left[\matrix{
    \partial_x\partial_x\operatorname{Re}\beta\tilde\mathcal S &
    \partial_y\partial_x\operatorname{Re}\beta\tilde\mathcal S \cr
    \partial_x\partial_y\operatorname{Re}\beta\tilde\mathcal S &
    \partial_y\partial_y\operatorname{Re}\beta\tilde\mathcal S
  }\right]
\end{equation}
This can be simplified using the fact that the action is holomorphic, which
means that it obeys the Cauchy--Riemann equations
\begin{equation}
  \partial_x\operatorname{Re}\tilde\mathcal S=\partial_y\operatorname{Im}\tilde\mathcal S
  \qquad
  \partial_y\operatorname{Re}\tilde\mathcal S=-\partial_x\operatorname{Im}\tilde\mathcal S
\end{equation}
Using these relationships alongside the Wirtinger derivative
$\partial\equiv\frac12(\partial_x-i\partial_y)$ allows the order of the
derivatives and the real or imaginary parts to be commuted, with
\begin{eqnarray}
  \partial_x\operatorname{Re}\tilde\mathcal S=\operatorname{Re}\partial\tilde\mathcal S
  \qquad
  \partial_y\operatorname{Re}\tilde\mathcal S=-\operatorname{Im}\partial\tilde\mathcal S \\
  \partial_x\operatorname{Im}\tilde\mathcal S=\operatorname{Im}\partial\tilde\mathcal S
  \qquad
  \partial_y\operatorname{Im}\tilde\mathcal S=\operatorname{Re}\partial\tilde\mathcal S
\end{eqnarray}
Using these relationships, the hessian \eref{eq:real.hessian} can be written in
the more manifestly complex way
\begin{eqnarray}
  \operatorname{Hess}_{x,y}\operatorname{Re}\beta\mathcal S
  &=\left[\matrix{
    \hphantom{-}\operatorname{Re}\beta\partial\partial\tilde\mathcal S &
    -\operatorname{Im}\beta\partial\partial\tilde\mathcal S \cr
    -\operatorname{Im}\beta\partial\partial\tilde\mathcal S &
    -\operatorname{Re}\beta\partial\partial\tilde\mathcal S
  }\right] \\
  &=\left[\matrix{
      \hphantom{-}\operatorname{Re}\beta\operatorname{Hess}\mathcal S &
    -\operatorname{Im}\beta\operatorname{Hess}\mathcal S \cr
    -\operatorname{Im}\beta\operatorname{Hess}\mathcal S &
    -\operatorname{Re}\beta\operatorname{Hess}\mathcal S
  }\right]
\end{eqnarray}
where $\operatorname{Hess}\mathcal S$ is the hessian with respect to $z$ given
in \eqref{eq:complex.hessian}.

The eigenvalues and eigenvectors of the Hessian are important for evaluating
thimble integrals, because those associated with upward directions provide a
local basis for the surface of the thimble. Suppose that $v_x,v_y\in\mathbb
R^N$ are such that
\begin{equation}
  (\operatorname{Hess}_{x,y}\operatorname{Re}\beta\mathcal S)\left[\matrix{v_x \cr v_y}\right]
  =\lambda\left[\matrix{v_x \cr v_y}\right]
\end{equation}
where the eigenvalue $\lambda$ must be real because the hessian is real symmetric. The problem can be put into a more obviously complex form by a change of basis. Writing $v=v_x+iv_y$, we find
\begin{eqnarray}
  &\left[\matrix{0&(i\beta\operatorname{Hess}\mathcal S)^*\cr i\beta\operatorname{Hess}\mathcal S&0}\right]
  \left[\matrix{v \cr iv^*}\right]\\
  &\qquad=\left[\matrix{1&i\cr i&1}\right]
  (\operatorname{Hess}_{x,y}\operatorname{Re}\beta\mathcal S)
  \left[\matrix{1&i\cr i&1}\right]^{-1}
  \left[\matrix{1&i\cr i&1}\right]
  \left[\matrix{v_x \cr v_y}\right] \\
  &\qquad=\lambda\left[\matrix{1&i\cr i&1}\right]\left[\matrix{v_x \cr v_y}\right]
  =\lambda\left[\matrix{v \cr iv^*}\right]
\end{eqnarray}
It therefore follows that the eigenvalues and vectors of the real hessian satisfy the equation
\begin{equation} \label{eq:generalized.eigenproblem}
  \beta\operatorname{Hess}\mathcal S v=\lambda v^*
\end{equation}
a sort of generalized
eigenvalue problem. If we did not know the eigenvalues were real, we could
still see it from the second implied equation,
$(\beta\operatorname{Hess}\mathcal S)^*v^*=\lambda v$, which is the conjugate
of the first if $\lambda^*=\lambda$.

Something somewhat hidden in the structure of the real hessian but more clear
in its complex form is that each eigenvalue comes in a pair, since
\begin{equation}
  \beta\operatorname{Hess}\mathcal S(iv)=i\lambda v^*=-\lambda(iv)
\end{equation}
Therefore, if $\lambda$ is an eigenvalue of the hessian with eigenvector $v$,
than so is $-\lambda$, with associated eigenvector $iv$, rotated in the complex
plane. It follows that each stationary point has an equal number of descending
and ascending directions, e.g. the index of each stationary point is $N$. For a
stationary point in a real problem this might seem strange, because there are
clear differences between minima, maxima, and saddles of different index.
However, we can quickly see here that for a such a stationary point, its $N$
real eigenvectors which determine its index in the real problem are accompanied
by $N$ purely imaginary eigenvectors, pointing into the complex plane and each
with the negative eigenvalue of its partner.

The effect of changing the phase of $\beta$ is revealed by
\eqref{eq:generalized.eigenproblem}. Writing $\beta=|\beta|e^{i\phi}$ and
dividing both sides by $|\beta|e^{i\phi/2}$, one finds
\begin{equation}
  \operatorname{Hess}\mathcal S(e^{i\phi/2}v)
  =\frac{\lambda}{|\beta|}e^{-i\phi/2}v^*
  =\frac{\lambda}{|\beta|}(e^{i\phi/2}v)^*
\end{equation}
Therefore, one only needs to consider solutions to the eigenproblem for the
action alone, $\operatorname{Hess}\mathcal Sv_0=\lambda_0 v_0^*$, and then rotate the
resulting vectors by a constant phase corresponding to half the phase of
$\beta$ or $v(\phi)=v_0e^{-i\phi/2}$. One can see this in the examples of Figs. \ref{fig:1d.stokes} and
\ref{fig:thimble.orientation}, where increasing the argument of $\beta$ from
left to right produces a clockwise rotation in the thimbles in the
complex-$\theta$ plane.

These eigenvalues and vectors can be further related to properties of the
complex symmetric matrix $\beta\operatorname{Hess}\mathcal S$. Suppose that
$u\in\mathbb R^N$ satisfies the eigenvalue equation
\begin{equation}
  (\beta\operatorname{Hess} S)^\dagger(\beta\operatorname{Hess} S)u
  =\sigma u
\end{equation}
for some positive real $\sigma$ (real because $(\beta\operatorname{Hess}
S)^\dagger(\beta\operatorname{Hess} S)$ is self-adjoint). The square root of these
numbers, $\sqrt{\sigma}$, are the definition of the \emph{singular values} of
$\beta\operatorname{Hess}\mathcal S$. A direct relationship between these singular
values and the eigenvalues of the hessian immediately follows by taking an
eigenvector $v\in\mathbb C$ that satisfies \eref{eq:generalized.eigenproblem},
and writing
\begin{eqnarray}
  \sigma v^\dagger u
  &=v^\dagger(\beta\operatorname{Hess} S)^\dagger(\beta\operatorname{Hess} S)u
  =(\beta\operatorname{Hess} Sv)^\dagger(\beta\operatorname{Hess} S)u\\
  &=(\lambda v^*)^\dagger(\beta\operatorname{Hess} S)u
  =\lambda v^T(\beta\operatorname{Hess} S)u
  =\lambda^2 v^\dagger u
\end{eqnarray}
Thus if $v^\dagger u\neq0$, $\lambda^2=\sigma$. It follows that the eigenvalues
of the real hessian are the singular values of the complex matrix
$\beta\operatorname{Hess}\mathcal S$, and their eigenvectors coincide up to a
constant complex factor.

\subsection{The conditions for Stokes points}

As we have seen, gradient descent on the real part of the action results in a
flow which preserves the imaginary part of the action. Therefore, a necessary
condition for the existence of a Stokes line between two stationary points is
for those points to have the same imaginary action. However, this is not a
sufficient condition. This can be seen in Fig.~\ref{fig:4.spin}, which shows
the thimbles of the circular 6-spin model. The argument of $\beta$ has been
chosen such that the stationary points marked by $\clubsuit$ and
$\blacktriangle$ have exactly the same imaginary energy, and yet they do not
share a thimble.

\begin{figure}
  \includegraphics{figs/6_spin.pdf}
  \caption{
  Some thimbles of the circular 6-spin model, where the argument of $\beta$ has
  been chosen such that the imaginary parts of the action at the stationary
  points $\clubsuit$ and $\blacktriangle$ are exactly the same (and, as a
  result of conjugation symmetry, the points $\bigstar$ and $\blacksquare$).
  } \label{fig:4.spin}
\end{figure}

This is because these stationary points are not adjacent: they are separated
from each other by the thimbles of other stationary points. This is a
consistent story in one complex dimension, since the codimension of the
thimbles is the same as the codimension of the constant imaginary energy
surface is one, and such a surface can divide space into regions. However, in higher dimensions thimbles do not have codimension high enough to divide space into regions.

\begin{figure}
  \includegraphics{figs/2_spin_thimbles.pdf}
  \caption{
    Thimbles of the $N=3$ spherical 2-spin model projected into the
    $\operatorname{Re}\theta$, $\operatorname{Re}\phi$,
    $\operatorname{Im}\theta$ space. The blue and green lines trace gradient
    descent of the two minima, while the red and orange lines trace those of
    the two saddles. The location of the maxima are marked as points, but their
    thimbles are not shown.
  }
\end{figure}

\subsection{Evaluating thimble integrals}

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 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
\begin{equation}
  Z_\sigma(\beta)=\oint_{\mathcal J_\sigma}ds\,e^{-\beta\mathcal S(s)}
\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 stationary point, where the real part of the action is by
construction at its minimum on the thimble and the integrand is therefore
largest.

We will make a change of coordinates $u(s)\in\mathbb R^D$ such that
\begin{equation} \label{eq:thimble.integration.def}
  \beta\mathcal S(s)=\beta\mathcal S(s_\sigma)+\frac{|\beta|}2 u(s)^Tu(s)
\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 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
part is at its minimum.

The coordinates $u$ can be constructed implicitly in the close vicinity of the stationary point, with
\begin{equation}
  s(u)=s_\sigma+\sum_{i=1}^{D}\sqrt{\frac{|\beta|}{\lambda^{(i)}}}v^{(i)}u_i+O(u^2)
\end{equation}
where the sum is over pairs $(\lambda, v)$ which satisfy
\eqref{eq:generalized.eigenproblem} and have $\lambda>0$. It is straightforward
to confirm that these coordinates satisfy \eqref{eq:thimble.integration.def} asymptotically close to the stationary point, as
\begin{eqnarray}
  \beta\mathcal S(s(u))
  &=\beta\mathcal S(s_\sigma)
  +\frac12(s(u)-s_\sigma)^T(\beta\operatorname{Hess}\mathcal S)(s(u)-s_\sigma)+\cdots \\
  &=\beta\mathcal S(s_\sigma)
  +\frac{|\beta|}2\sum_{ij}\frac{v^{(i)}_k}{\sqrt{\lambda^{(i)}}}(\beta[\operatorname{Hess}\mathcal S]_{k\ell})\frac{v^{(j)}_\ell}{\sqrt{\lambda^{(j)}}}u_iu_j+\cdots \\
  &=\beta\mathcal S(s_\sigma)
  +\frac{|\beta|}2\sum_{ij}\frac{v^{(i)}_k}{\sqrt{\lambda^{(i)}}}\frac{\lambda^{(j)}(v^{(j)}_k)^*}{\sqrt{\lambda^{(j)}}}u_iu_j+\cdots \\
  &=\beta\mathcal S(s_\sigma)
  +\frac{|\beta|}2\sum_{ij}\frac{\sqrt{\lambda^{(j)}}}{\sqrt{\lambda^{(i)}}}\delta_{ij}u_iu_j+\cdots \\
  &=\beta\mathcal S(s_\sigma)
  +\frac{|\beta|}2\sum_iu_i^2+\cdots
\end{eqnarray}
The Jacobian of this transformation is
\begin{equation}
  \frac{\partial s_i}{\partial u_j}=\sqrt{\frac{|\beta|}{\lambda^{(j)}}}v^{(j)}_i
  =\sqrt{\frac1{\lambda_0^{(j)}}}v^{(j)}_i
\end{equation}
where $\lambda_0=\lambda/|\beta|$ is the eigenvalue of the hessian for
$\beta=1$. This is a $N\times D$ matrix. Since the eigenvectors $v$ are
mutually complex orthogonal, $v_i^{(j)}$ is nearly a unitary matrix, and it can
be made one by including $v^{(N)}=\widehat{\partial g}$, the unit normal to the
constraint manifold. This lets us write $U_{ij}=v_i^{(j)}$ an $N\times N$
unitary matrix, whose determinant will give the correct phase for the measure.

We therefore have
\begin{equation}
Z_\sigma(\beta)=e^{-\beta\mathcal S(s_\sigma)}\int du\,\det\frac{ds}{du}e^{-\frac{|\beta|}2u^Tu}
\end{equation}
Now we take the saddle point approximation, assuming the integral is dominated
by its value at the stationary point such that the determinant can be
approximated by its value at the stationary point. This gives
\begin{eqnarray}
  Z_\sigma(\beta)
    &\simeq e^{-\beta\mathcal S(s_\sigma)}\left.\det\frac{ds}{du}\right|_{s=s_\sigma}\int du\,e^{-\frac{|\beta|}2u^Tu} \\
    &=e^{-\beta\mathcal S(s_\sigma)}\left(\prod_i^D\sqrt{\frac1{\lambda_0^{(i)}}}\right)\det U\left(\frac{2\pi}{|\beta|}\right)^{D/2}
\end{eqnarray}
We are left with evaluating the determinant of the unitary part of the coordinate transformation.
In circumstances you may be used to, only the absolute value of the determinant
from the coordinate transformation is relevant, and since the determinant of a
unitary matrix is always magnitude one, it doesn't enter the computation.
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 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
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 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}

\begin{eqnarray}
  Z(\beta)^*
  =\sum_{\sigma\in\Sigma_0}n_\sigma Z_\sigma(\beta)^*
  =\sum_{\sigma\in\Sigma_0}n_\sigma(-1)^{k_\sigma}Z_\sigma(\beta^*)
  =Z(\beta^*)
\end{eqnarray}

\section{The ensemble of symmetric complex-normal matrices}

We will now begin dealing with the implications of actions defined in very many
dimensions. We saw in \S\ref{sec:stationary.hessian} that the singular values
of the complex hessian of the action at any stationary point are important in
the study of thimbles. Hessians are symmetric matrices by construction. For
real actions of real variables, the study random symmetric matrices with
Gaussian entries provides insight into a wide variety of problems. In our case,
we will find the relevant ensemble is that of random symmetric matrices with
\emph{complex-normal} entries. In this section, we will introduce this
distribution, review its known properties, and derive its singular value
distribution in the large-matrix limit.

The complex normal distribution with zero mean is the unique Gaussian
distribution in one complex variable $Z$ whose variances are
$\overline{Z^*Z}=\overline{|Z|^2}=\Gamma$ and $\overline{Z^2}=C$. $\Gamma$ is
positive, and $|C|\leq\Gamma$. The special case of $C=\Gamma$, where the
variance of the complex variable and its covariance with its conjugate are the
same, reduces to the ordinary normal distribution. The case where $C=0$ results
in the real and imaginary parts of $Z$ being uncorrelated, in what is known as
the standard complex normal distribution. Its probability density function is
defined by
\begin{equation}
  p(z\mid\Gamma,C)=
  \frac1{\pi\sqrt{\Gamma^2-|C|^2}}\exp\left\{
    \frac12\left[\matrix{z^*&z}\right]\left[\matrix{
        \Gamma & C \cr C^* & \Gamma
    }\right]^{-1}\left[\matrix{z\cr z^*}\right]
  \right\}
\end{equation}

We will consider an ensemble of random $N\times N$ matrices $B=A+\lambda_0I$, where the
entries of $A$ are complex-normal distributed with variances
$\overline{|A|^2}=A_0/N$ and $\overline{A^2}=C_0/N$, and $\lambda_0$ is a
constant shift to its diagonal. The eigenvalue distribution of these matrices
is already known to take the form of an elliptical ensemble, with constant
support inside the ellipse defined by
\begin{equation} \label{eq:ellipse}
  \left(\frac{\operatorname{Re}(\lambda e^{i\theta})}{1+|C_0|/A_0}\right)^2+
  \left(\frac{\operatorname{Im}(\lambda e^{i\theta})}{1-|C_0|/A_0}\right)^2
  <A_0
\end{equation}
where $\theta=\frac12\arg C_0$ \cite{Nguyen_2014_The}. The eigenvalue
spectrum of $B$ is therefore constant inside the same ellipse
translated so that its center lies at $\lambda_0$.  Examples of these
distributions are shown in the insets of Fig.~\ref{fig:spectra}.

When $C=0$ and the elements of $A$ are standard complex normal, the singular
value distribution of $B$ is a complex Wishart distribution. For $C\neq0$ the
problem changes, and to our knowledge a closed form is not in the literature.
We have worked out an implicit form for the singular value spectrum using the
replica method.

The singular values of $B$ correspond with the square-root of the eigenvalues
of $B^\dagger B$, but also they correspond to the absolute value of the
eigenvalues of the real $2N\times2N$ block matrix
\begin{equation}
  \left[\matrix{\operatorname{Re}B&-\operatorname{Im}B\cr-\operatorname{Im}B&-\operatorname{Re}B}\right]
\end{equation}
as we saw in \S\ref{sec:stationary.hessian}. The eigenvalue spectrum of this
block matrix can be studied by ordinary means.  Defining the `partition function'
\begin{equation}
  Z(\sigma)=\int dx\,dy\,\exp\left\{
    -\frac12\left[\matrix{x&y}\right]
    \left(\sigma I-
    \left[\matrix{\hphantom{-}\operatorname{Re}B&-\operatorname{Im}B\cr-\operatorname{Im}B&-\operatorname{Re}B}\right]
  \right)
    \left[\matrix{x\cr y}\right]
  \right\}
\end{equation}
implies a Green function
\begin{equation}
  G(\sigma)=\frac\partial{\partial\sigma}\log Z(\sigma)
\end{equation}
This can be put into a manifestly complex form in the same way it was done in \S\ref{sec:stationary.hessian}, using the same linear transformation of $x$ and $y$ into $z$ and $z^*$. This gives
\begin{eqnarray}
  Z(\sigma)
  &=\int dz\,dz^*\,\exp\left\{
    -\frac12\left[\matrix{z^*&-iz}\right]
    \left(\sigma I-
    \left[\matrix{0&(iB)^*\cr iB&0}\right]
  \right)
    \left[\matrix{z\cr iz^*}\right]
  \right\} \\
  &=\int dz\,dz^*\,\exp\left\{
    -\frac12\left(
      2z^\dagger z\sigma-z^\dagger B^*z^*-z^TBz
    \right)
  \right\} \\
  &=\int dz\,dz^*\,\exp\left\{
    -z^\dagger z\sigma+\operatorname{Re}(z^TBz)
  \right\}
\end{eqnarray}
which is a general expression for the singular values $\sigma$ of a symmetric
complex matrix $B$.


Introducing replicas to bring the partition function into the numerator of the
Green function \cite{Livan_2018_Introduction} gives
\begin{equation} \label{eq:green.replicas}
  G(\sigma)=\lim_{n\to0}\int dz\,dz^*\,(z^{(0)})^\dagger z^{(0)}
    \exp\left\{
    -\sum_\alpha^n\left[(z^{(\alpha)})^\dagger z^{(\alpha)}\sigma
      +\operatorname{Re}\left((z^{(\alpha)})^TBz^{(\alpha)}\right)
    \right]
  \right\},
\end{equation}
The average is then made over
$J$ and Hubbard--Stratonovich is used to change variables to the replica matrices
$N\alpha_{\alpha\beta}=(z^{(\alpha)})^\dagger z^{(\beta)}$ and
$N\chi_{\alpha\beta}=(z^{(\alpha)})^Tz^{(\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+\frac18A_0\alpha_0^2-\frac{\alpha_0\sigma}2+\frac12\log(\alpha_0^2-|\chi_0|^2)
    +\frac14\operatorname{Re}\left(\frac14C_0^*\chi_0^2+\lambda_0^*\chi_0\right)
  \right]\right\}.
\end{equation}

\begin{figure}
  \centering

  \includegraphics{figs/spectra_0.0.pdf}
  \includegraphics{figs/spectra_0.5.pdf}\\
  \includegraphics{figs/spectra_1.0.pdf}
  \includegraphics{figs/spectra_1.5.pdf}

  \caption{
    Eigenvalue and singular value spectra of a random matrix $A+\lambda_0I$, where the entries of $A$ are complex-normal distributed with $\overline{|A|^2}=A_0=5/4$ and $\overline{A^2}=C_0=\frac34e^{-i3\pi/4}$.
    The diaginal shifts differ in each plot, with (a) $\lambda_0=0$, (b)
    $\lambda_0=\frac12|\lambda_{\mathrm{th}}|$, (c)
    $\lambda_0=|\lambda_{\mathrm{th}}|$, and (d)
    $\lambda_0=\frac32|\lambda_{\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.
  } \label{fig:spectra}
\end{figure}

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

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

\section{The \textit{p}-spin spherical models}

The $p$-spin spherical models are defined by the action
\begin{equation} \label{eq:p-spin.hamiltonian}
  \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}
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 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 pure $2$-spin model is defined by
\begin{equation}
  \mathcal S_2(z)=\frac12z^TJz.
\end{equation}
$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}
  \tilde\mathcal S(z)=\mathcal S_2(z)-\epsilon(z^Tz-N)
\end{equation}
Stationary points must satisfy
\begin{equation}
  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 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$ stationary points, each
corresponding to $\pm$ the cardinal directions in the diagonalized basis.

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
interest assuming all others are zero are
\begin{eqnarray}
  \dot z_1
  &=-z_1^*\left(d_1^*-\frac{d_1^*z_1^*z_1+d_2^*z_2^*z_2}{|z_1|^2+|z_2|^2}\right) \\
  &=-(d_1-d_2)^*z_1^*\frac{|z_2|^2}{|z_1|^2+|z_2|^2}
\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
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 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}.

\begin{figure}
  \begin{gnuplot}[terminal=epslatex, terminaloptions={size 8.65cm,5.35cm}]
    set xlabel '$\Delta t$'
    set ylabel '$z(t) / \sqrt{N}$'

    plot 1 / sqrt(1 + exp(2 * x)) t '$z_1$', \
         1 / sqrt(1 + exp(- 2 * x)) t '$z_2$'
  \end{gnuplot}
  \caption{
    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 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 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
stationary points are brought to the same value, their surfaces of constant
imaginary energy join.

\begin{eqnarray}
  Z(\beta)
  &=\int_{S^{N-1}}ds\,e^{-\beta H(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} \\
  &=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}\log\frac12|\epsilon_k-\epsilon_\ell|
  \right\}
\end{eqnarray}
\begin{eqnarray}
  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')\log\frac12|\epsilon-\epsilon'|
  \right\}
\end{eqnarray}

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
\begin{equation}
  \rho(\epsilon\mid\epsilon_{\mathrm{th}})=\frac2{\pi\epsilon_{\mathrm{th}}^2}\sqrt{\epsilon_{\mathrm{th}}^2-\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}}
  \right)
\end{equation}
Finally, the product over the singular values corresponding to descending directions gives
\begin{equation}
  \frac12\int d\epsilon'\,\rho(\epsilon')\log\frac12|\epsilon-\epsilon'|
  =\log\frac{\epsilon_{\mathrm{th}}}2-\frac14+\frac12\left(\frac{\epsilon}{\epsilon_{\mathrm{th}}}\right)^2
\end{equation}
for $\epsilon<\epsilon_{\mathrm{th}}$. Then
\begin{equation}
  \operatorname{Re}f=-\epsilon\operatorname{Re}\beta-\log\frac12+\frac14-\frac12\epsilon^2
\end{equation}
\begin{equation}
  \operatorname{Im}f=-\epsilon\operatorname{Im}\beta+\frac12\left(
    \epsilon\sqrt{1-\epsilon^2}+2\tan^{-1}\frac{1+\epsilon}{\sqrt{1-\epsilon^2}}
  \right)
\end{equation}
The value of the integral will be dominated by the contribution near the
maximum of the real part of $f$, which is
\begin{equation}
  \epsilon_{\mathrm{max}}=\left\{\matrix{-\operatorname{Re}\beta & \operatorname{Re}\beta\leq1\cr -1 & \mathrm{otherwise}}\right.
\end{equation}
For $\operatorname{Re}\beta>1$, the maximum is concentrated in the ground state
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}}$.

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
to zero. This happens for
\begin{equation}
  0
  =\frac{\partial}{\partial\epsilon}\operatorname{Im}f\Big|_{\epsilon=\epsilon_{\mathrm{max}}}
  =-\operatorname{Im}\beta+\sqrt{1-\epsilon_{\mathrm{max}}^2}
  =-\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
\cite{Obuchi_2012_Partition-function, Takahashi_2013_Zeros}.

\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.

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
    = \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
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
\begin{equation} \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\tilde\mathcal S_p\gamma
    \right)
  \right\}
\end{equation}
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{equation}
  \phi(1)=z+\bar\theta(1)\eta+\gamma\theta(1)+\hat z\bar\theta(1)\theta(1)
\end{equation}
and its measure $d\phi=dz\,d\hat z\,d\eta\,d\gamma$.
Then the expression for the number of stationary points can be written in a very compact form, as
\begin{equation}
  \mathcal N=\int d\phi^*d\phi\,\exp\left\{
    \int d1\,\operatorname{Re}
      \tilde\mathcal S_p(\phi(1))
  \right\}
\end{equation}
where $d1=d\bar\theta(1)\,d\theta(1)$ denotes the integration over the grasssman 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.

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}
  Q(1,2)=\frac1N\left[\matrix{
      \phi(1)^T\phi(2)&\phi(1)^T\phi(2)^*\cr
      \phi(1)^\dagger\phi(2)&\phi(1)^\dagger\phi(2)^*
  }\right]
\end{equation}
the result can be written, neglecting constant factors,
\begin{equation}
  \overline\mathcal N\simeq\int dQ\,e^{NS_\mathrm{eff}(Q)}
\end{equation}
for an effective action functional of the supermatrix $Q$
\begin{equation}
  S_{\mathrm{eff}}=
  \int d1\,d2\,\operatorname{Tr}\left(
      \frac14\left[
        \matrix{\frac14&\frac14\cr\frac14&\frac14}
      \right]Q^{(p)}(1,2)-\frac p2\left[
        \matrix{\frac\epsilon2&0\cr0&\frac{\epsilon^*}2}
      \right](Q(1,1)-I)\delta(1,2)
    \right)
    +\frac12\log\det Q
\end{equation}
where 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.
This leads to the condition for a saddle point of
\begin{equation}
  0
  =\frac{\partial S_\mathrm{eff}}{\partial Q(1,2)}
  =\frac p{16}Q^{(p-1)}(1,2)-\frac p2\left[
    \matrix{\frac\epsilon2&0\cr0&\frac{\epsilon^*}2}
  \right]\delta(1,2)
  +\frac12Q^{-1}(1,2)
\end{equation}
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
\begin{eqnarray}
  0
  &=\int d3\,\frac{\partial S_\mathrm{eff}}{\partial Q(1,3)}Q(3,2) \\
  &=\frac p{16}\int d3\,Q^{(p-1)}(1,3)Q(3,2)-\frac p2\left[
  \matrix{\frac\epsilon2&0\cr0&\frac{\epsilon^*}2}
  \right]Q(1,2)+\frac12I\delta(1,2)
\end{eqnarray}
When expanded, the supermatrix $Q$ 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?}

The problem of counting the density of Stokes points in an analytic
continuation of the spherical models is quite challenging, as the problem of
finding dyramic trajectories with endpoints at stationary points is already
difficult, and once made complex the problem has twice the number of fields
squared.

In this section, we begin to address the problem heuristically by instead
asking: if you are at a stationary point, where are your neighbors? The
stationary points geometrically nearest to a given stationary point should make
up the bulk of its adjacent points in the sense of being susceptible to Stokes
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
\begin{eqnarray}
  \mathcal N
  &=\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{eqnarray}
\begin{equation}
  Q(1,2)
  =\left[
    \matrix{
      \phi_1(1)^T\phi_1(2) & \phi_1(1)^T\phi_2(2) & \phi_1(1)^T\phi_2(2)^* \cr
      \phi_2(1)^T\phi_1(2) & \phi_2(1)^T\phi_2(2) & \phi_2(1)^T\phi_2(2)^* \cr
      \phi_2(1)^\dagger\phi_1(2) & \phi_2(1)^\dagger\phi_2(2) & \phi_2(1)^\dagger\phi_2(2)^*
    }
  \right]
\end{equation}
\begin{equation}
  S_\mathrm{eff}
  =\int d1\,d2\,\operatorname{Tr}\left\{
    \frac14\left[\matrix{1&\frac12&\frac12\cr\frac12&\frac14&\frac14\cr\frac12&\frac14&\frac14}\right]Q^{(p)}(1,2)
    -\frac p2\left[
    \matrix{\epsilon_1&0&0\cr0&\frac12\epsilon_2&0\cr0&0&\frac12\epsilon_2^*}
      \right]\left(
      Q(1,1)-I
    \right)\delta(1,2)
  \right\}+\frac12\det Q
\end{equation}
\begin{equation}
  0=\frac{\partial S_\mathrm{eff}}{\partial Q(1,2)}
  =
    \frac p4\left[\matrix{1&\frac12&\frac12\cr\frac12&\frac14&\frac14\cr\frac12&\frac14&\frac14}\right]\odot Q^{(p-1)}(1,2)
    -\frac p2\left[
  \matrix{\epsilon_1&0&0\cr0&\frac12\epsilon_2&0\cr0&0&\frac12\epsilon_2^*}\right]\delta(1,2)
  +\frac12Q^{-1}(1,2)
\end{equation}
where $\odot$ denotes element-wise multiplication.
\begin{eqnarray}
  0
  &=\int d3\,\frac{\partial S_\mathrm{eff}}{\partial Q(1,3)}Q(3,2) \\
  &=\frac p4\int d3\,
    \left\{\left[\matrix{1&\frac12&\frac12\cr\frac12&\frac14&\frac14\cr\frac12&\frac14&\frac14}\right]\odot Q^{(p-1)}(1,3)\right\}Q(3,2)
    -\frac p2 \left[
  \matrix{\epsilon_1&0&0\cr0&\frac12\epsilon_2&0\cr0&0&\frac12\epsilon_2^*}\right]Q(1,2)
  +\frac12I\delta(1,2)
\end{eqnarray}
Despite being able to pose the saddle point problem in a compact way, a great
deal of complexity lies within. The supermatrix $Q$ depends on 35 independent
bilinear products, and when the superfields are expanded produces 48 (not entirely independent) equations.
These equations can be split into 30 involving bilinear products of the
fermionic fields and 18 without them. The 18 equations without fermionic
bilinear products can be solved with a computer algebra package to eliminate 17
of the 20 non-fermionic bilinear products. The fermionic equations are
unfortunately more complicated.

They can be simplified somewhat by examination of the real two-replica problem.
There, all bilinear products involving fermionic fields from different
replicas, like $\eta_1^T\eta_2$, vanish. This is related to the influence of
the relative position of the two replicas to their spectra, with the vanishing
being equivalent to having no influence, e.g., the value of the determinant at
each stationary point is exactly what it would be in the one-replica problem
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
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}
\begin{equation}
  \Delta=2(1+|y_2|^2-\sqrt{1-|y_2|^2}\cos\theta_{xx})
\end{equation}
\begin{eqnarray}
  \delta\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{eqnarray}
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$.

\begin{figure}
  \includegraphics{figs/bound.pdf}
  \includegraphics{figs/example_bound.pdf}

  \caption{
    The line bounding $\delta$ 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
    $\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.
  }
\end{figure}

\section{The $p$-spin spherical models: numerics}

To confirm the presence of Stokes lines under certain processes in the $p$-spin, we studied the problem numerically.


\bibliographystyle{unsrt}
\bibliography{stokes}

\appendix

\section{Geometry}

The surface $M\subset\mathbb C^N$ defined by $N=f(z)=z^2$ is an $N-1$ dimensional
\emph{Stein manifold}, a type of complex manifold defined by the level set of a
holomorphic function \cite{Forstneric_2017_Stein}. One can define a Hermitian
metric $h$ on $M$ by taking the restriction of the standard metric of $\mathbb
C^N$ to vectors tangent along $M$. For any smooth function $\phi:M\to\mathbb
R$, its gradient is a holomorphic vector field given by
\begin{equation}
  \operatorname{grad}\phi=\bar\partial^\sharp\phi
\end{equation}
Suppose $u:M\to\mathbb C^{N-1}$ is a coordinate system. Then
\begin{equation}
  \operatorname{grad}\phi=h^{\bar\beta\alpha}\bar\partial_{\bar\beta}\phi\frac{\partial}{\partial u^\alpha}
\end{equation}
Let $z=u^{-1}$.
\begin{equation}
  \frac\partial{\partial u^\alpha}=\frac{\partial z^i}{\partial u^\alpha}\frac\partial{\partial z^i}
\end{equation}
\begin{equation}
  \bar\partial_{\bar\beta}\phi=\frac{\partial\bar z^i}{\partial\bar u^{\bar\beta}}\frac{\partial\phi}{\partial\bar z^i}
\end{equation}
\begin{equation}
  \operatorname{grad}\phi
  =\frac{\partial\bar z^{\bar\jmath}}{\partial\bar u^{\bar\beta}}h^{\bar\beta\alpha}\frac{\partial z^i}{\partial u^\alpha}\frac{\partial\phi}{\partial\bar z^{\bar\jmath}}\frac\partial{\partial z^i}
\end{equation}
At any point $z_0\in M$, $z_0$ is parallel to the normal to $M$. Without loss of generality, take $z_0=e^N$. In a neighborhood of $z_0$, the map $u^\alpha=z^\alpha$ is a coordinate system.
its inverse is $z^i=u^i$ for $1\leq i\leq N-1$ and
\begin{equation}
  z^N=\sqrt{N-u^2}.
\end{equation}
The Jacobian is
\begin{equation}
  \frac{\partial z^i}{\partial u^\alpha}=\delta^i_\alpha-\delta_N^{i}\frac{u^\alpha}{\sqrt{N-u^2}}
\end{equation}
and therefore the Hermitian metric induced by the map is
\begin{equation}
  h_{\alpha\bar\beta}=\frac{\partial\bar z^i}{\partial\bar u^\alpha}\frac{\partial z^{\bar\jmath}}{\partial u^{\bar\beta}}\delta_{i\bar\jmath}
  =\delta_{\alpha\bar\beta}+\frac{\bar u^{\alpha}u^{\bar\beta}}{|N-u^2|}
\end{equation}
The metric can be inverted explicitly:
\begin{equation}
  h^{\bar\beta\alpha}
  =\delta^{\bar\beta\alpha}-\frac{\bar u^{\bar\beta}u^\alpha}{|N-u^2|+|u|^2}.
\end{equation}
Putting these pieces together, we find
\begin{equation}
  \frac{\partial\bar z^{\bar\jmath}}{\partial\bar u^{\bar\beta}}h^{\bar\beta\alpha}\frac{\partial z^i}{\partial u^\alpha}
  =\delta^{\bar\jmath i}-\frac{z^{\bar\jmath}\bar z^i}{|z|^2}
\end{equation}
\begin{equation}
  \operatorname{grad}\phi
  =\left(\delta^{\bar\jmath i}-\frac{z^{\bar\jmath}\bar z^i}{|z|^2}\right)
  \frac{\partial\phi}{\partial\bar z^{\bar\jmath}}\frac\partial{\partial z^i}
\end{equation}

\section{Numerics}

To study Stokes lines numerically, we approximated them by parametric curves.
If $z_0$ and $z_1$ are two stationary points of the action with
$\operatorname{Re}\mathcal S(z_0)>\operatorname{Re}\mathcal S(z_1)$, then we
take the curve
\begin{equation}
  z(t)
  =(1-t)z_0+tz_1+(1-t)t\sum_{i=0}^mg_it^i
\end{equation}
where the $g$s are undetermined complex vectors. These are fixed by minimizing
a cost function, which has a global minimum only for Stokes lines. Defining
\begin{equation}
  \mathcal L(t)
  = 1-\frac{\operatorname{Re}[\dot z^*(z(t))\cdot z'(t)]}{|\dot z(z(t))||z'(t)|}
\end{equation}
this cost is given by
\begin{equation}
  \mathcal C=\int_0^1 dt\,\mathcal L(t)
\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
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$,
and approximating the cost integral by a finite sum.

\end{document}