\documentclass[12pt]{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 \usepackage{pifont} \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 among the lowest minima 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{Witten_2011_Analytic}. 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 oscillate and the results continued to the desired model \cite{Alexandru_2022_Complex}. 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 track. 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 a phase space average of some observable $\mathcal O$ of the form \begin{equation} \label{eq:observable} \langle\mathcal O\rangle=\frac1Z\int_\Omega ds\,e^{-\beta\mathcal S(s)}\mathcal O(s) \end{equation} where the partition function $Z$ normalizes the average as \begin{equation} \label{eq:partition.function} Z=\int_\Omega ds\,e^{-\beta\mathcal S(s)} \end{equation} Rather than focus on any specific observable, we will study the partition function itself, since it exhibits the essential features that readily generalize to arbitrary observable averages. We've defined $Z$ in a way that 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 on which we will focus). In this paper we will consider only analytic continuation of the parameter $\beta$, but any other parameter 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. An example action used throughout this section is shown in Fig.~\ref{fig:example.action}. \begin{figure} \hspace{5pc} \includegraphics{figs/circle.pdf}\hfill \includegraphics{figs/action.pdf}\hfill \includegraphics{figs/stationaryPoints.pdf} \caption{ An example of a simple action and its stationary points. \textbf{Left:} The configuration space of the $N=2$ spherical (or circular) model, defined for $s\in\mathbb R^N$ restricted to the circle $N=s^2$. It can be parameterized by one angle $\theta=\arctan(s_2/s_1)$. Its natural complex extension takes instead $s\in\mathbb C^N$ restricted to the hyperbola $N=s^2=(\operatorname{Re}s)^2-(\operatorname{Im}s)^2$. The (now complex) angle $\theta$ is still a good parameterization of phase space. \textbf{Center:} An action $\mathcal S$ for circular $3$-spin model, defined 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$, plotted as a function of $\theta$. \textbf{Right:} The stationary points of $\mathcal S$ in the complex-$\theta$ plane. In this example, $\Sigma=\{\mbox{\ding{117}},{\mbox{\ding{72}}},{\mbox{\ding{115}}},{\mbox{\ding{116}}},{\mbox{\ding{108}}},{\mbox{\ding{110}}}\}$ and $\Sigma_0=\{\mbox{\ding{117}},{\mbox{\ding{116}}}\}$. 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} In order to analytically continue \eref{eq:partition.function}, $\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, and replacing its real arguments with complex ones. 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 often 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$. 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 ds\,e^{-\beta\mathcal S(s)} \end{equation} For the moment this translation has only changed one 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$. This is illustrated in Fig.~\ref{fig:contour.deformation}. \begin{figure} \hspace{5pc} \includegraphics{figs/hyperbola_1.pdf}\hfill \includegraphics{figs/hyperbola_2.pdf}\hfill \includegraphics{figs/hyperbola_3.pdf} \hspace{5pc} \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. } \label{fig:contour.deformation} \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 each is defined by all points that approach the stationary point $s_\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(s_\sigma)$ at the stationary point, since all other points on the thimble must descend to reach it. And, as we will see in the following subsection, 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 its integral is to converge, since these are the places where the real part of the action is very large and the real part of 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} \hspace{5pc} \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 $s_\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$. \begin{figure} \hspace{5pc} \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 $\mbox{\ding{117}}$ in the direction implied by its orientation, and then follow the thimble from the maximum ${\mbox{\ding{116}}}$ \emph{against} the direction implied by its orientation, from top to bottom. Therefore, $\mathcal C=\mathcal J_{\mbox{\ding{117}}}-\mathcal J_{\mbox{\ding{116}}}$. \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_{\mbox{\ding{117}}}+\mathcal J_{\mbox{\ding{116}}}$. } \label{fig:thimble.orientation} \end{figure} 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}ds\,e^{-\beta\mathcal S(s)}. \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} \hspace{5pc} \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_{\mbox{\ding{117}}}$ intersects the stationary point ${\mbox{\ding{115}}}$ 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_{\mbox{\ding{115}}}$ must be included. Notice that in this figure, because of the symmetry of the pure models, the thimble $\mathcal J_{\mbox{\ding{110}}}$ 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} 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. All one needs to work out the projection operator is a coordinate system. Here we sketch the derivation for the spherical models. 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} D_\alpha z^i=\frac{\partial z^i}{\partial u^\alpha}=\delta^i_\alpha-\delta_N^{i}\frac{u^\alpha}{\sqrt{N-u^2}} \end{equation} Taking the appropriate combination of $Dz$ yields \begin{equation} P=I-\frac{zz^\dagger}{|z|^2} \end{equation} One can quickly verify that this operator indeed projects the dynamics onto the manifold: the vector perpendicular to the manifold 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$, the identity. \begin{figure} \hspace{5pc} \hfill \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{equation} \eqalign{ \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{equation} 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{equation} \eqalign{ \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{equation} Using these relationships, the hessian \eref{eq:real.hessian} can be written in the more manifestly complex way \begin{equation} \eqalign{ \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{equation} 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{equation} \eqalign{ &\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{equation} 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, whose solutions are called \emph{Takagi values and vectors} \cite{Takagi_1924_On}. 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{equation} \eqalign{ \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{equation} 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 ${\mbox{\ding{115}}}$ have exactly the same imaginary energy, and yet they do not share a thimble. \begin{figure} \hspace{5pc} \hfill \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 ${\mbox{\ding{115}}}$ are exactly the same (and, as a result of conjugation symmetry, the points ${\mbox{\ding{72}}}$ and ${\mbox{\ding{110}}}$). } \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} \hspace{5pc} \hfill \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. For thorough account of evaluating these integrals (including \emph{at} Stokes points), see Howls \cite{Howls_1997_Hyperasymptotics}. 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{equation} \eqalign{ \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{equation} 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{equation} \eqalign{ 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{equation} 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 before any Stokes points, \begin{equation} \label{eq:real.thimble.partition.function} Z_\sigma(\beta)\simeq\left(\frac{2\pi}\beta\right)^{D/2}i^{k_\sigma}|\det\operatorname{Hess}\mathcal S(s_\sigma)|^{-\frac12}e^{-\beta\mathcal S(s_\sigma)} \end{equation} \begin{equation} 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{equation} \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 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{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} \\ &=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} \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\} \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 \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')\log2|\epsilon-\epsilon'| =-\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+\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} \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} \fl \eqalign{ 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) \\ &\hspace{28em}+\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{equation} \eqalign{ 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{equation} 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 $(\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 \begin{equation} \Sigma = \log(p-1)-\frac12\log\left( \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}) \end{equation} where \begin{equation} R_\pm^2=\frac{p-1}2\frac{(r^{p-2}\pm1)\left[ r^{2(p-1)}\pm(p-1)r^{p-2}(r^2-1)-1 \right]}{ 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 \begin{equation} \eqalign{ I_p(u) &=\left(\frac12+\frac1{r^{p-2}-1}\right)^{-1}(\operatorname{Re}u)^2 -\left(\frac12-\frac1{r^{p-2}+1}\right)^{-1}(\operatorname{Im}u)^2 \\ &\qquad-\log\left( r^{p-2}\left| u+\sqrt{u^2-1} \right|^2 \right)+2\operatorname{Re} \left( u\sqrt{u^2-1} \right) } \end{equation} otherwise. The branch of the square roots are chosen such that the real part of the root has the opposite sign as the real part of $u$, e.g., if $\operatorname{Re}u<0$ then $\operatorname{Re}\sqrt{u^2-1}>0$. If the real part is zero, then the sign is taken so that the imaginary part of the root has the opposite sign of the imaginary part of $u$. \begin{figure} \hspace{4pc} \includegraphics{figs/re_complexity.pdf} \hspace{-2pc} \includegraphics{figs/im_complexity.pdf} \includegraphics{figs/leg_complexity.pdf} \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 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} \end{figure} Contours of this complexity for the pure 3-spin are plotted in Fig.~\ref{fig:p-spin.complexity} for pure and imaginary energy. The thick black 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. \begin{figure} \hspace{2pc} \includegraphics{figs/ground_complexity.pdf} \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$ 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 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. } \label{fig:ground.complexity} \end{figure} 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 complexity of rank one saddles becomes zero \cite{Auffinger_2012_Random}. There are several interesting features of the complexity. First is this inequivalence between the real complexity and the limit of the complex complexity to zero complex part. It implies, among other things, a desert of stationary points in the complex plan surrounding the lowest minima, something we shall see more explicitly in the next section. Second, there is only a small collection of stationary points that appear with positive complexity and a gapped spectrum: the small region in Fig.~\ref{fig:ground.complexity} that is both to the right of the thick line and brightly shaded. We suspect that these are the only stationary points that are somewhat protected from participation in Stokes points. \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 dynamic 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{equation} \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{equation} \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} \eqalign{ 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) \right. \\ &\qquad\qquad\left.-\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} \eqalign{ 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) \\ &\qquad\qquad-\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{equation} \eqalign{ 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) \\ &\qquad\qquad\qquad\qquad -\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{equation} 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{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}) \\ &=2(1-\sqrt{1-|y_2|^2}\cos\theta_{xx}-i|y_2|\cos\theta_{xy}) } \end{equation} 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} \hspace{5pc} \includegraphics{figs/bound.pdf} \hfill \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} A lot of information is contained in the full two-replica complexity, but we will focus on the following question: what does the population of stationary points nearby a given real stationary point look like? We think this is a relevant question for the tendency for Stokes lines, for the following reason. To determine whether two given stationary points, when tuned to have the same imaginary energy, will share a Stokes line, one needs to solve what is known as the global connection problem. As we have seen, this as a question of a kind of 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} 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 \begin{equation} \varphi=\arctan\sqrt{\frac{1+\operatorname{Re}\delta}{1-\operatorname{Re}\delta}} \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. \begin{figure} \hspace{5pc} \includegraphics{figs/neighbor_geometry.pdf} \caption{ The geometric definition of the angle $\varphi$, between the displacement between two stationary points and the real configuration space. } \end{figure} First, we examine the important 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 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 becomes extremely flat. \begin{figure} \includegraphics{figs/neighbor_thres.pdf} \caption{ The scaled two-replica complexity $\Upsilon$ 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$. \begin{figure} \includegraphics{figs/neighbor_limit_thres_above.pdf} \hspace{-1em} \includegraphics{figs/neighbor_limit_thres_at.pdf} \hspace{-1em} \includegraphics{figs/neighbor_limit_thres_below.pdf} \hfill \includegraphics{figs/neighbor_limit_thres_legend.pdf} \caption{ The scaled two-replica complexity $\Upsilon$ 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:} $\epsilon=\epsilon_\mathrm{th}-0.001$. All lines have been normalized by the complexity $\Sigma$ of the real 3-spin model at the same energy, or (where relevant) by the complexity $\Sigma_2$ of rank-two saddles of the real 3-spin model. } \end{figure} Thus, there is an abrupt geometric transition in the population of nearest neighbors as the threshold is crossed: above they are broadly distributed at all angles, while below they are highly concentrated around $90^\circ$. From this analysis it appears that the complexity of the nearest neighbors, at zero distance, behaves as that of the index-2 saddles at all angles, which would imply that the nearest neighbors vanish at the same point as the index-2 saddles. However, this is not the case: we have only shown that this is how the neighbors at \emph{identical energy} scale, which is correct above the threshold, but no longer underneath. If an energy is taken under the threshold and the two-replica complexity maximized with respect to both $\epsilon_2$ and $\varphi$, one finds that as $\Delta\to0$, $\epsilon_2\to\epsilon_1$, as must be the case the find a 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 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 angle $\varphi$ for various $\Delta$ at $\epsilon_1=\mathcal E_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$ which leads to the largest maximum value. As $\Delta$ varies this varies like $\epsilon_2-\epsilon_1\propto\Delta^2$. \textbf{Dashed lines:} The complexity evaluated at $\epsilon_2=\epsilon_1$. } \end{figure} Below $\mathcal E_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) \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 $\mathcal E_1$, $\epsilon_2-\epsilon_1\propto -|\epsilon_1-\mathcal E_1|^2$, $\varphi-90^\circ\propto-|\epsilon_1-\mathcal E_1|^{1/2}$, and $\Delta_\textrm{min}\propto|\epsilon_1-\mathcal E_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 involve upward flow from a point that already has weight, neighbors that have a lower energy won't be eligible to be involved in a Stokes line that causes a change of weight until the phase of $\beta$ has rotated almost $180^\circ$. \begin{figure} \includegraphics{figs/neighbor_closest_energy.pdf} \caption{ The energy $\epsilon_2$ of the nearest neighbor stationary points in the complex plane to a given real stationary point of energy $\epsilon_1$. The dashed line shows $\epsilon_2=\epsilon_1$. The nearest neighbor energy coincides with the dashed line until $\mathcal E_1$, the energy where rank-one saddles vanish, where it peels off. } \end{figure} \begin{figure} \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 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 all angles, but there are exponentially more found at $90^\circ$ than at any other angle. Below $\mathcal E_2$ but above $\mathcal E_1$, stationary points are found at arbitrarily close distance but only at $90^\circ$. Below $\mathcal E_1$, neighboring stationary points are separated by a minimum squared distance $\Delta_\textrm{min}$, and the angle they are found at drifts. } \label{fig:nearest.properties} \end{figure} \section{The {\it p}-spin spherical models: 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_iP_i^{(1,1)}(2t-1) \end{equation} where the $g$s are undetermined complex vectors and the $P_i^{(1,1)}(x)$ are the Jacobi polynomials, orthogonal on the interval $[-1,1]$ under the weight $(1-x)(1+x)$. The Jacobi polynomials are used because they are orthogonal with respect to integration over precisely the term they appear inside above. 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))^\dagger z'(t)]}{|\dot z(z(t))||z'(t)|} \end{equation} where $\dot z$ is the flow at $z$ given by \eqref{eq:flow}, this cost is given by \begin{equation} \label{eq:cost} \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. 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 the \emph{real} configuration space of the $p$-spin model. Then, a saddle-finding routine is run on the complex configuration space in the close vicinity of the real saddle, using random initial conditions in a slowly increasing radius of the real stationary point. When this process finds a new distinct stationary point, it is done. This method of sampling pairs heavily biases the statistics we report here in favor of seeing Stokes points. Once a pair of nearby stationary points has been found, one real and one in the complex plane, their energies are used to compute the phase $\theta$ necessary to give $\beta$ in order to set their imaginary energies to the same value, the necessary conditions for a Stokes line. A straight line (ignoring even the constraint) is thrown between them and then minimized using the cost function \eqref{eq:cost} for some initial $m=5$. Once a minimum is found, $m$ is iteratively increased several times, each time minimizing the cost in between, until $m=20$. If at some point in this process the cost blows up, indicating that the solution is running away, the pair is thrown out; this happens infrequently. At the end, there are several ways to asses whether a given minimized line is a Stokes line: the value of the cost, the integrated deviation from the constraint, the integrated deviation from the phase of the two stationary points. Among minimized lines these values fall into doubly-peaked histograms which well-separated prospective Stokes lines into `good' and `bad' values for the given level of approximation $m$. One cannot explicitly study the effect of crossing various landmark energies on the $p$-spin in the system sizes that were accessible to our study, up to around $N=64$, as the presence of, e.g., the threshold energy separating saddles from minima is not noticeable until much larger size \cite{Folena_2020_Rethinking}. However, we are able to examine the effect of its symptoms: namely, the influence of the spectrum of the stationary point in question on the likelihood that a randomly chosen neighbor will share a Stokes line. \begin{figure} \includegraphics{figs/numerics_prob_eigenvalue.pdf} \caption{ The probability $P_\mathrm{Stokes}$ that a real stationary point will share a Stokes line with its randomly chosen neighbor as a function of $|\lambda_\textrm{min}|$, the magnitude of the minimum eigenvalue of the hessian at the real stationary point. The horizontal axis has been rescaled to collapse the data at different system sizes $N$. } \label{fig:numeric.prob.eigenvalue} \end{figure} Data for the likelihood of a Stokes line as a function of the empirical gap $|\lambda_\textrm{min}|$ of the real stationary point is shown in Fig.~\ref{fig:numeric.prob.eigenvalue}. There, one sees that the probability of finding a Stokes line with a near neighbor falls off as an exponential in the magnitude of the smallest eigenvalue. As a function of system size, the tail contracts like $N^{-1/2}$, which means that in the thermodynamic limit one expects the probability of finding such a Stokes line will approach zero everywhere expect where $\lambda_\textrm{min}\ll1$. This supports the idea that gapped minima are unlikely to see Stokes lines. \begin{figure} \includegraphics{figs/numerics_angle_gap_32.pdf} \caption{ The probability density function for identified Stokes points as a function of $|\theta|$, the magnitude of the phase necessary to add to $\beta$ to reach the Stokes point, at $N=32$ and for several binned $|\lambda_\textrm{min}|$. As the empirical gap is increased, the population of discovered Stokes points becomes more concentrated around $|\theta|=\pi$. } \label{fig:numeric.angle.gap} \end{figure} We can also see that as the empirical gap is increased, Stokes points tend to occur at very large phases. This can be seen for $N=32$ in Fig.~\ref{fig:numeric.angle.gap}, which shows the probability distribution of Stokes lines discovered as a function of phase $|\theta|$ necessary to reach them. The curves are broken into sets representing different bins of the empirical gap $|\lambda_\textrm{min}|$. As the empirical gap grows, Stokes points become depleted around small phases and concentrate on very large ones. This supports the idea that around the gapped minima, Stokes points will be concentrated at phases that are nearly $180^\circ$, where the two-replica calculation shows that almost all of their nearest neighbors will lie. \subsection{Pure {\it p}-spin: is analytic continuation possible?} After this work, one is motivated to ask: can analytic continuation be done in even a simple complex model like the pure $p$-spin? Numeric and analytic evidence indicates that the project is hopeless if ungapped stationary points take a significant weight in the partition function, since for these Stokes lines proliferate at even small continuation and there is no hope of tracking them. However, for gapped stationary points we have seen compelling evidence that suggests they will not participate in Stokes points, at least until a large phase rotation of the parameter being continued. This gives some hope for continuation of the low-temperature thermodynamic phase of the $p$-spin, where weight is concentrated in precisely these points. Recalling our expression of the single-thimble contribution to the partition function for a real stationary point of a real action expanded to lowest order in large $|\beta|$ \eqref{eq:real.thimble.partition.function}, we can write for the $p$-spin after an infintesimal rotation of $\beta$ into the complex plane (before any Stokes points have been encountered) \begin{equation} \eqalign{ Z(\beta) &=\sum_{\sigma\in\Sigma_0}n_\sigma Z_\sigma(\beta) \\ &\simeq\sum_{\sigma\in\Sigma_0}\left(\frac{2\pi}\beta\right)^{D/2}i^{k_\sigma} |\det\operatorname{Hess}\mathcal S(s_\sigma)|^{-\frac12} e^{-\beta\mathcal S(s_\sigma)} \\ &\simeq\sum_{k=0}^D\int d\epsilon\,\mathcal N_\mathrm{typ}(\epsilon,k) \left(\frac{2\pi}\beta\right)^{D/2}i^k |\det\operatorname{Hess}\mathcal S(s_\sigma)|^{-\frac12} e^{-\beta N\epsilon} } \end{equation} where $\mathcal N_\mathrm{typ}(\epsilon,k)$ is the typical number of stationary points in a sample of the real $p$-spin model in the energy range $\epsilon$ to $\epsilon+d\epsilon$ and with index $k$. Following Derrida \cite{Derrida_1991_The}, this is related to the \emph{average} number of stationary points in this range at large $N$ by \begin{equation} \mathcal N_\mathrm{typ}(\epsilon,k)=\overline\mathcal N(\epsilon,k)+\eta(\epsilon,k)\overline\mathcal N(\epsilon,k)^{1/2} \end{equation} where $\eta$ is a random, sample-dependant number of order one. This gives two terms to the typical partition function \begin{equation} Z_\mathrm{typ}=Z_A+Z_B \end{equation} where \begin{eqnarray} \fl 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 \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} \\ =\int d\epsilon\,\tilde\eta(\epsilon)e^{Nf_B(\epsilon)} \end{eqnarray} for functions $f_A$ and $f_B$ defined by \begin{eqnarray} f_A &=-\beta\epsilon+\Sigma(\epsilon)-\frac12\int d\lambda\,\rho(\lambda\mid\epsilon)|\lambda|+\frac12\log\frac{2\pi}\beta +i\frac\pi2P(\lambda<0\mid\epsilon) \\ f_B &=-\beta\epsilon+\frac12\Sigma(\epsilon)-\frac12\int d\lambda\,\rho(\lambda\mid\epsilon)|\lambda|+\frac12\log\frac{2\pi}\beta +i\frac\pi2P(\lambda<0\mid\epsilon) \end{eqnarray} and $P(\lambda<0\mid\epsilon)$ is the cumulative probability distribution of the eigenvalues of the spectrum given $\epsilon$, \begin{equation} P(\lambda<0\mid\epsilon)=\int_{-\infty}^0 d\lambda'\,\rho(\lambda'\mid\epsilon) \end{equation} and produces the macroscopic index $k/N$. Each integral will be dominated by its value near the maximum of the real part of the exponential argument. Assuming that $\epsilon<\epsilon_\mathrm{th}$, this maximum occurs at \begin{equation} 0=\frac{\partial}{\partial\epsilon}\operatorname{Re}f_A=-\operatorname{Re}\beta-\frac12\frac{3p-4}{p-1}\epsilon+\frac12\frac p{p-1}\sqrt{\epsilon^2-\epsilon_\mathrm{th}^2} \end{equation} \begin{equation} 0=\frac{\partial}{\partial\epsilon}\operatorname{Re}f_B=-\operatorname{Re}\beta-\epsilon \end{equation} As with the 2-spin model, the integral over $\epsilon$ is oscillatory and can only be reliably evaluated with a saddle point when either the period of oscillation diverges \emph{or} when the maximum lies at a cusp. We therefore expect changes in behavior when $\epsilon=\mathcal E_0$, the ground state energy. The temperature at which this happens is \begin{eqnarray} \operatorname{Re}\beta_A&=-\frac12\frac{3p-4}{p-1}\mathcal E_0+\frac12\frac p{p-1}\sqrt{\mathcal E_0^2-\mathcal E_\mathrm{th}^2}\\ \operatorname{Re}\beta_B&=-\mathcal E_0 \end{eqnarray} which for all $p\geq2$ has $\operatorname{Re}\beta_A>\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(\mathcal E_0)$). It is only at $\operatorname{Re}\beta_B=-\mathcal E_0$ where both terms contributing to the partition function at large $N$ involve incoherent integrals near the maximum, and only here where the density of zeros is expected to become nonzero. In fact, in the limit of $|\beta|\to\infty$, $\operatorname{Re}\beta_B$ is precisely the transition found in \cite{Obuchi_2012_Partition-function} between phases with and without a density of zeros. This value is an underestimate for the transition for finite $|\beta|$, which likely results from the invalidity of our large-$\beta$ approximation. More of the phase diagram might be constructed by continuing the series for individual thimbles to higher powers in $\beta$, which would be equivalent to allowing non-constant terms in the Jacobian over the thimble. \begin{figure} \hspace{5pc} \hfill\includegraphics{figs/obuchi_3-spin.pdf} \caption{ Phases of the 3-spin model in the complex-$\beta$, following Obuchi \& Takahashi \cite{Obuchi_2012_Partition-function}. The phase P2 contains a nonzero density of zeros of the partition function, while the `spin-glass' phase SG does not. Analytic continuation via thimbles correctly predicts the boundary between these two phases when $|\beta|\gg1$ to be $\operatorname{Re}\beta=-\epsilon_0$, shown with a dashed line. } \label{fig:obuchi_3-spin} \end{figure} This zeroth-order analysis for the $p$-spin suggests that analytic continuation can be sometimes done despite the presence of a great many complex stationary points. In particular, when weight is concentrated in certain minima Stokes lines do not appear to interrupt the proceedings. How bad the situation in is other regimes, like for smaller $|\beta|$, remains to be seen: our analysis cannot tell between the effects of Stokes points changing the contour and the large-$|\beta|$ saddle-point used to evaluate the thimble integrals. Taking the thimbles to the next order in $\beta$ may reveal more explicitly where Stokes points become important. \section{Conclusion} \section*{References} \bibliographystyle{unsrt} \bibliography{stokes} \end{document}