Complex complex landscapes

Jaron Kent-Dobias and Jorge Kurchan

Laboratoire de Physique de l'Ecole Normale Supérieure, 75005 Paris, France

<figure><img src="image_1.png" /></figure>

(Received 21 December 2020; accepted 1 April 2021; published 22 April 2021)

We study the saddle points of the p-spin model—the best understood example of a “complex” (rugged) landscape—when its N variables are complex. These points are the solutions to a system of N random equations of degree p − 1. We solve for $\overline{N}$, the number of solutions averaged over randomness in the $N \to \infty$ limit. We find that it saturates the Bézout bound $\ln \overline{N} \sim N \ln(p-1)$. The Hessian of each saddle is given by a random matrix of the form $C^†C$, where C is a complex symmetric Gaussian matrix with a shift to its diagonal. Its spectrum has a transition where a gap develops that generalizes the notion of “threshold level” well known in the real problem. The results from the real problem are recovered in the limit of real parameters. In this case, only the square root of the total number of solutions are real. In terms of the complex energy, the solutions are divided into sectors where the saddles have different topological properties.

DOI: 10.1103/PhysRevResearch.3.023064

I. INTRODUCTION

Spin glasses are the paradigm of many-variable “complex landscapes,” a category that also includes neural networks and optimization problems like constraint satisfaction [1]. The most tractable family of these are the mean-field spherical p-spin models [2] (for a review see [3]) defined by the energy

$$
H_0 = \frac{1}{p!} \sum_{i_1 \cdots i_p}^{N} J_{i_1 \cdots i_p} z_{i_1} \cdots z_{i_p}, \quad (1)
$$

where J is a symmetric tensor whose elements are real Gaus-
sian variables and z ∈ R N is constrained to the sphere z T z =
N. This problem has been studied in the algebra [4] and proba-
bility literature [5,6]. It has been attacked from several angles:
the replica trick to compute the Boltzmann–Gibbs distribution
[2], a Kac–Rice [7–9] procedure (similar to the Fadeev–Popov
integral) to compute the number of saddle points of the en-
ergy function [10], and gradient-descent (or more generally
Langevin) dynamics starting from a high-energy configura-
tion [11]. Thanks to the simplicity of the energy, all these
approaches yield analytic results in the large- N limit.

In this paper we extend the study to complex variables: we shall take $z \in C^N$ and $J$ to be a symmetric tensor whose elements are complex normal, with $|J|^2 = p!/(2N^{p-1})$ and $\bar{J}^2 = \kappa|J|^2$ for complex parameter $|\kappa| < 1$. The constraint remains $z^T z = N$.

The motivations for this paper are of three types. On the
practical side, there are indeed situations in which complex
variables appear naturally in disordered problems: such is

the case in which the variables are *phases*, as in random
laser problems [12]. Quiver Hamiltonians—used to model
black hole horizons in the zero-temperature limit—also have
a Hamiltonian very close to ours [13]. A second reason is that,
as we know from experience, extending a real problem to the
complex plane often uncovers underlying simplicity that is
otherwise hidden, shedding light on the original real problem,
e.g., as in the radius of convergence of a series.

Finally, deforming an integral in N real variables to a surface of dimension N in 2N-dimensional complex space has turned out to be necessary for correctly defining and analyzing path integrals with complex action (see [14,15]), and as a useful palliative for the sign problem [16–18]. In order to do this correctly, features of the action’s landscape in complex space—such as the relative position of saddles and the existence of Stokes lines joining them—must be understood. This is typically done for simple actions with few saddles, or for a target phenomenology with symmetries that restrict the set of saddles to few candidates. Given the recent proliferation of “glassiness” in condensed matter and high energy physics, it is inevitable that someone will want to apply these methods to a system with a complex landscape, and will find they cannot use approaches that rely on such assumptions. Their landscape may not be random: here we follow the standard strategy of computer science by understanding the generic features of random instances of a simple case, expecting that this sheds light on practical, nonrandom problems. While in this paper we do not yet address analytic continuation of integrals, understanding the distribution and spectra of critical points is an essential first step.

Returning to our problem, the spherical constraint is en-
forced using the method of Lagrange multipliers: introducing
$\epsilon \in C$, our constrained energy is

$$
H = H_0 + \frac{p}{2} \epsilon \left( N - \sum_{i}^{N} z_i^2 \right). \qquad (2)
$$

<footnote>Published by the American Physical Society under the terms of the
Creative Commons Attribution 4.0 International license. Further
distribution of this work must maintain attribution to the author(s)
and the published article's title, journal citation, and DOI.</footnote>
One might balk at the constraint $z^T z = N$—which could appropriately be called a *hyperbolic* constraint—by comparison with $z^\dagger z = N$. The reasoning behind the choice is twofold.

First, we seek to draw conclusions from our model that are applicable to generic holomorphic functions without any symmetry. Samples of $H_0$ nearly provide this, save for a single anomaly: the value of the energy and its gradient at any point $z$ correlate along the $z$ direction, with $\overline{H_0 \partial H_0} \propto \overline{H_0 (\partial H_0)^*} \propto z$. This anomalous direction should thus be forbidden, and the constraint surface $z^T z = N$ accomplishes this.

Second, taking the constraint to be the level set of a holomorphic function means the resulting configuration space is a *bona fide* complex manifold, and therefore permits easy generalization of the integration techniques referenced above. The same cannot be said for the space defined by $z^\dagger z = N$, which is topologically the $(2N - 1)$-sphere and cannot admit a complex structure.

Imposing the constraint with a holomorphic function makes the resulting configuration space a *bona fide* complex manifold, which is, as we mentioned, the situation we wish to model. The same cannot be said for the space defined by $z^\dagger z = N$, which is topologically the $(2N - 1)$-sphere, does not admit a complex structure, and thus yields a trivial structure of saddles. However, we will introduce the bound $r^2 \equiv z^\dagger z/N \le R^2$ on the “radius” per spin as a device to classify saddles. We shall see that this “radius” $r$ and its upper bound $R$ are insightful parameters in our present problem, revealing structure as they are varied. Note that taking $R = 1$ reduces the problem to that of the ordinary $p$ spin.

The critical points are of $H$ given by the solutions to

$$ \frac{p}{p!} \sum_{j_1 \cdots j_{p-1}}^{N} J_{i_{j_1} \cdots j_{p-1}} z_{j_1} \cdots z_{j_{p-1}} = p \epsilon z_i \quad (3) $$

for all $i = \{1, \dots, N\}$, which for fixed $\epsilon$ is a set of $N$ equations of degree $p-1$, to which one must add the constraint. In this sense this study also provides a complement to the work on the distribution of zeros of random polynomials [19], which are for $N=1$ and $p \to \infty$. We see from (3) that at any critical point $\epsilon = H_0/N$, the average energy.

Since $H$ is holomorphic, any critical point of $\mathrm{Re} H$ is also one of $\mathrm{Im} H$, and therefore of $H$ itself. Writing $z = x + iy$ for $x, y \in \mathbb{R}^N$, $\mathrm{Re} H$ can be considered a real-valued function of $2N$ real variables. The number of critical points of $H$ is thus given by the usual Kac–Rice formula applied to $\mathrm{Re} H$:

$$ \mathcal{N}(\kappa, \epsilon, R) = \int dx dy \delta(\partial_x \mathrm{Re} H) \delta(\partial_y \mathrm{Re} H) \times \left| \det \begin{bmatrix} \partial_x \partial_x \mathrm{Re} H & \partial_x \partial_y \mathrm{Re} H \\ \partial_y \partial_x \mathrm{Re} H & \partial_y \partial_y \mathrm{Re} H \end{bmatrix} \right|. \quad (4) $$

This expression is to be averaged over $J$ to give the complexity $\Sigma$ as $N\Sigma = \ln N$, a calculation that involves the replica trick. Based on the experience from similar problems [3], the *annealed approximation* $N\Sigma \sim \ln N$ is expected to be exact wherever the complexity is positive.

The Cauchy–Riemann equations may be used to write (4) in a manifestly complex way. With the Wirtinger derivative $\partial = \frac{1}{2}(\partial_x - i\partial_y)$, one can write $\partial_x \mathrm{Re} H = \mathrm{Re} \partial H$ and $\partial_y \mathrm{Re} H = -\mathrm{Im} \partial H$. Carrying these transformations through,

one finds

$$ \begin{aligned} \mathcal{N}(\kappa, \epsilon, R) &= \int dx dy \delta(\mathrm{Re} \partial H) \delta(\mathrm{Im} \partial H) \\ & \quad \times \left| \det \begin{bmatrix} \mathrm{Re} \partial \partial H & -\mathrm{Im} \partial \partial H \\ -\mathrm{Im} \partial \partial H & -\mathrm{Re} \partial \partial H \end{bmatrix} \right| \\ &= \int dx dy \delta(\mathrm{Re} \partial H) \delta(\mathrm{Im} \partial H) \left| \det[(\partial \partial H)^\dagger \partial \partial H] \right| \\ &= \int dx dy \delta(\mathrm{Re} \partial H) \delta(\mathrm{Im} \partial H) \left| \det \partial \partial H \right|^2. \end{aligned} \quad (5) $$

This gives three equivalent expressions for the determinant of the Hessian: as that of a $2N \times 2N$ real symmetric matrix, that of the $N \times N$ Hermitian matrix $(\partial \partial H)^\dagger \partial \partial H$, or the norm squared of that of the $N \times N$ complex symmetric matrix $\partial \partial H$.

These equivalences belie a deeper connection between the spectra of the corresponding matrices. Each positive eigenvalue of the real matrix has a negative partner. For each such pair $\pm\lambda$, $\lambda^2$ is an eigenvalue of the Hermitian matrix and $|\lambda|$ is a *singular value* of the complex symmetric matrix. The distribution of positive eigenvalues of the Hessian is therefore the same as the distribution of singular values of $\partial \partial H$, or the distribution of square-rooted eigenvalues of $(\partial \partial H)^\dagger \partial \partial H$.

## II. CALCULATION

A useful property of the Gaussian $J$ is that gradient and Hessian at fixed energy $\epsilon$ are statistically independent [9,20], so that the $\delta$ functions and the Hessian may be averaged independently. First we shall compute the spectrum of the Hessian, which can in turn be used to compute the determinant. Then we will treat the $\delta$ functions and the resulting saddle point equations.

The Hessian $\partial \partial H = \partial \partial H_0 - p\epsilon I$ is equal to the unconstrained Hessian with a constant added to its diagonal. The eigenvalue distribution $\rho$ is therefore related to the unconstrained distribution $\rho_0$ by a similar shift: $\rho(\lambda) = \rho_0(\lambda + p\epsilon)$. The Hessian of the unconstrained Hamiltonian is

$$ \partial_i \partial_j H_0 = \frac{p(p-1)}{p!} \sum_{k_1 \cdots k_{p-2}}^{N} J_{i j k_1 \cdots k_{p-2}} z_{k_1} \cdots z_{k_{p-2}}, \quad (6) $$

which makes its ensemble that of Gaussian complex symmetric matrices, when the anomalous direction normal to the constraint surface is neglected. Given its variances $|\partial_i \partial_j H_0|^2 = p(p-1)r^{p-2}/2N$ and $(\partial_i \partial_j H_0)^2 = p(p-1)\kappa/2N$, $\rho_0(\lambda)$ is constant inside the ellipse

$$ \left( \frac{\mathrm{Re}(\lambda e^{i\theta})}{r^{p-2} + |\kappa|} \right)^2 + \left( \frac{\mathrm{Im}(\lambda e^{i\theta})}{r^{p-2} - |\kappa|} \right)^2 < \frac{p(p-1)}{2r^{p-2}}, \quad (7) $$

where $\theta = \frac{1}{2} \arg \kappa$ [21]. The eigenvalue spectrum of $\partial \partial H$ is therefore constant inside the same ellipse translated so that its center lies at $-p\epsilon$. Examples of these distributions are shown in the insets of Fig. 1.

The eigenvalue spectrum of the Hessian of the real part is not the spectrum $\rho(\lambda)$ of $\partial \partial H$, but instead the square-root eigenvalue spectrum of $(\partial \partial H)^\dagger \partial \partial H$; in other words, the singular value spectrum $\rho(\sigma)$ of $\partial \partial H$. When $\kappa = 0$ and the elements of $J$ are standard complex normal, this is a complex
<figure><img src="image_2.png" /><figcaption>FIG. 1. Eigenvalue and singular value spectra of the Hessian $\partial\partial H$ of the 3-spin model with $\kappa = \frac{3}{4}e^{-i3\pi/4}$. Pictured distributions are for critical points at 'radius' $r = \sqrt{5}/4$ and with energy per spin (a) $\epsilon = 0$, (b) $\epsilon = -\frac{1}{2}|\epsilon_{th}|$, (c) $\epsilon = -|\epsilon_{th}|$, and (d) $\epsilon = -\frac{3}{2}|\epsilon_{th}|$. The shaded region of each inset shows the support of the eigenvalue distribution (7). The solid line on each plot shows the distribution of singular values (10), while the overlaid histogram shows the empirical distribution from $2^{10} \times 2^{10}$ complex normal matrices with the same covariance and diagonal shift as $\partial\partial H$.</figcaption></figure>

Wishart distribution. For $\kappa \neq 0$ 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.

Introducing replicas to bring the partition function into the numerator of the Green function [22] gives

$$G(\sigma) = \lim_{n \to 0} \int d\zeta d\zeta^* (\zeta_i^{(0)})^* \zeta_i^{(0)} \\ \times \exp \left\{ \frac{1}{2} \sum_{\alpha}^{n} \left[ (\zeta_i^{(\alpha)})^* \zeta_i^{(\alpha)} \sigma - \text{Re} (\zeta_i^{(\alpha)} \zeta_j^{(\alpha)} \partial_i \partial_j H) \right] \right\}, \quad (8)$$

with sums taken over repeated latin indices. The average is then made over $J$ and Hubbard–Stratonovich is used to change variables to the replica matrices $N\alpha_{\alpha\beta} = (\zeta^{(\alpha)})^{\dagger}\zeta^{(\beta)}$ and $N\chi_{\alpha\beta} = (\zeta^{(\alpha)})^{\dagger}\zeta^{(\beta)}$, and a series of replica vectors. The replica-symmetric ansatz leaves all replica vectors zero, and $\alpha_{\alpha\beta} = \alpha_0\delta_{\alpha\beta}$, $\chi_{\alpha\beta} = \chi_0\delta_{\alpha\beta}$. The result is

$$\bar{G}(\sigma) = N \lim_{n \to 0} \int d\alpha_0 d\chi_0 d\chi_0^* \alpha_0 \\ \times \exp \left\{ nN \left[ 1 + \frac{p(p-1)}{16} r^{p-2} \alpha_0^2 - \frac{\alpha_0 \sigma}{2} \right. \right. \\ \left. + \frac{1}{2} \ln (\alpha_0^2 - |\chi_0|^2) \right. \\ \left. + \frac{p}{4} \text{Re} \left( \frac{(p-1)}{8} \kappa^* \chi_0^2 - \epsilon^* \chi_0 \right) \right] \right\}. \quad (9)$$

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

$$\rho(\sigma) = \frac{1}{i\pi N} \left( \lim_{\text{Im}\sigma \to 0^+} \bar{G}(\sigma) - \lim_{\text{Im}\sigma \to 0^-} \bar{G}(\sigma) \right). \quad (10)$$

Examples can be seen in Fig. 1 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 [23]. 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 (7) intersects the origin, and yields

$$|\epsilon_{\text{th}}|^2 = \frac{p-1}{2p} \frac{(1 - |\delta|^2)^2 r^{p-2}}{1 + |\delta|^2 - 2|\delta| \cos(\arg \kappa + 2 \arg \epsilon)} \quad (11)$$

for $\delta = \kappa r^{-(p-2)}$. Notice that the threshold depends on both the energy per spin $\epsilon$ on the “radius” $r$ of the saddle.

We will now address the $\delta$ functions of (5). These are converted to exponentials by the introduction of auxiliary fields $\hat{z} = \hat{x} + i\hat{y}$. The average over $J$ can then be performed. A generalized Hubbard–Stratonovich allows a change of variables from the $4N$ original and auxiliary fields to eight bilinears defined by $Nr = z^\dagger z$, $N\hat{r} = \hat{z}^\dagger \hat{z}$, $Na = \hat{z}^T \hat{z}$, $Nb = \hat{z}^T z$, and $N\hat{c} = \hat{z}^T \hat{z}$ (and their conjugates). The result, to leading order in $N$, is

$$\bar{N}(\kappa, \epsilon, R) = \int dr d\hat{r} da da^* db db^* d\hat{c} d\hat{c}^* e^{Nf(r, \hat{r}, a, b, \hat{c})}, \quad (12)$$

where the argument of the exponential is

$$f = 2 + \frac{1}{2} \ln \det \frac{1}{2} \begin{bmatrix} 1 & r & b & a \\ r & 1 & a^* & b^* \\ b & a^* & \hat{c} & \hat{r} \\ a & b^* & \hat{r} & \hat{c}^* \end{bmatrix} \\ + \int d\lambda d\lambda^* \rho(\lambda) \ln |\lambda|^2 + p \text{Re} \left\{ \frac{1}{8} [\hat{r} r^{p-1} + (p-1)|b|^2 r^{p-2} + \kappa(\hat{c}^* + (p-1)a^2)] - \epsilon a \right\}. \quad (13)$$

The spectrum $\rho$ is given in (7) and is dependant on $r$ alone. This function has an extremum in $\hat{r}$, $a$, $b$, and $\hat{c}$ at which its value is

$$f = 1 + \frac{1}{2} \ln \left( \frac{4}{p^2} \frac{r^2 - 1}{r^{2(p-1)} - |\kappa|^2} \right) \\ + \int d\lambda \rho(\lambda) \ln |\lambda|^2 - 2C_+[\text{Re}(\epsilon e^{-i\theta})]^2 \\ - 2C_-[\text{Im}(\epsilon e^{-i\theta})]^2, \quad (14)$$
<figure><img src="image_3.png" /><figcaption>FIG. 2. The complexity of the 3-spin model as a function of the maximum “radius” R at zero energy and several values of κ. The dashed line shows 1/2 ln(p-1), while the dotted shows ln(p-1).</figcaption></figure>

where $\theta = \frac{1}{2} \arg \kappa$ and

$$C_{\pm} = \frac{r^p(1 + p(r^2 - 1)) \mp r^2 |\kappa|}{r^{2p} \pm (p-1)r^p(r^2-1)|\kappa| - r^2|\kappa|^2} \quad (15)$$

Notice that level sets of f in energy $\epsilon$ also give ellipses, but of different form from the ellipse in (7).

### III. RESULTS

This expression is maximized for $r = R$, its value at the boundary, for all values of $\kappa$ and $\epsilon$. Evaluating the complexity at this saddle, in the limit of unbounded spins, gives

$$\lim_{R \to \infty} \ln \bar{N}(\kappa, \epsilon, R) = N \ln(p-1). \quad (16)$$

This is, to leading order, precisely the Bézout bound, the maximum number of solutions to N equations of degree $p-1$ [24]. That we saturate this bound is perhaps not surprising, since the coefficients of our polynomial equations (3) are complex and have no symmetries. Reaching Bézout in (16) is not our main result, but it provides a good check. Analogous asymptotic scaling has been found for the number of pure Higgs states in supersymmetric quiver theories [25].

For finite R, everything is analytically tractable at $\epsilon = 0$:

$$\Sigma(\kappa, 0, R) = \ln(p-1) - \frac{1}{2} \ln \left( \frac{1 - |\kappa|^2 R^{-4(p-1)}}{1 - R^{-4}} \right). \quad (17)$$

This is plotted as a function of R for several values of $\kappa$ in Fig. 2. For any $|\kappa| < 1$, the complexity goes to negative infinity as $R \to 1$, i.e., as the spins are restricted to the reals. This is natural, since volume of configuration space vanishes in this limit like $(R^2 - 1)^N$. However, when the result is analytically continued to $\kappa = 1$ (which corresponds to real J) something novel occurs: the complexity has a finite value at $R = 1$. This implies a $\delta$-function density of critical points on the $r = 1$ (or $y = 0$) boundary. The number of critical points contained there is

$$\lim_{R \to 1} \lim_{\kappa \to 1} \ln \bar{N}(\kappa, 0, R) = \frac{1}{2} N \ln(p-1), \quad (18)$$

half of (16) and corresponding precisely to the number of critical points of the real p-spin model. (Note the role of conjugation symmetry, already underlined in [19].) The full

<figure><img src="image_4.png" /><figcaption>FIG. 3. The value of bounding “radius” R for which $\Sigma(\kappa, \epsilon, R) = 0$ as a function of (real) energy per spin $\epsilon$ for the 3-spin model at several values of $\kappa$. Above each line the complexity is positive and critical points proliferate, while below it the complexity is negative and critical points are exponentially suppressed. The dotted black lines show the location of the ground and highest exited state energies for the real 3-spin model.</figcaption></figure>

$\epsilon$ dependence of the real p-spin is recovered by this limit as $\epsilon$ is varied.

In the thermodynamic limit, (17) implies that most critical points are concentrated at infinite radius $r$. For finite N the average radius of critical points is likewise finite. By differentiating $\bar{N}$ with respect to R and normalizing, one obtains the distribution of critical points as a function of $r$. This yields an average radius proportional to $N^{1/4}$. One therefore expects typical critical points to have a norm that grows modestly with system size.

These qualitative features carry over to nonzero $\epsilon$. In Fig. 3 we show that for $\kappa < 1$ there is always a gap in $r$ close to one in which solutions are exponentially suppressed. When $\kappa = 1$—the analytic continuation to the real computation—the situation is more interesting. In the range of energies where there are real solutions this gap closes, which is only possible if the density of solutions diverges at $r = 1$. Outside this range, around “deep” real energies where real solutions are exponentially suppressed, the gap remains. A moment’s thought tells us that this is necessary: otherwise a small perturbation of the J's could produce an unusually deep solution to the real problem, in a region where this should not happen.

The relationship between the threshold and ground, or extremal, state energies is richer than in the real case. In Fig. 4 these are shown in the complex-$\epsilon$ plane for several examples. Depending on the parameters, the threshold might have a smaller or larger magnitude than the extremal state, or cross as a function of complex argument. For sufficiently large $r$ the threshold is always at a larger magnitude. If this were to happen in the real case, it would likely imply our replica symmetric computation were unstable, since having a ground state above the threshold implies a ground state Hessian with many negative eigenvalues, a contradiction. However, this is not an contradiction in the complex case, where the energy is not bounded from below. The relationship between the threshold, i.e., where the gap appears, and the dynamics of, e.g., a minimization algorithm, deformed integration cycle, or physical dynamics, are a problem we hope to address in future work.
<figure><img src="image_5.png" /><figcaption>FIG. 4. Energies at which states exist (green shaded region) and threshold energies (black solid line) for the 3-spin model with $\kappa = \frac{3}{4}e^{-i3\pi/4}$ and (a) $r = \sqrt{2}$, (b) $r = \sqrt{1.325}$, (c) $r = \sqrt{1.125}$, and (d) $r = 1$. No shaded region is shown in (d) because no states exist at any energy.</figcaption></figure>

This paper provides a first step towards the study of complex landscapes with complex variables. The next obvious step is to study the topology of the critical points, the sets reached following gradient descent (the Lefschetz thimbles), and ascent (the antithimbles) [14–16,18,26], which act as constant-phase integrating “contours.” Locating and counting the saddles that are joined by gradient lines—the Stokes points, which play an important role in the theory—is also well within reach of the present-day spin-glass literature techniques. We anticipate that the threshold level, where the system develops a mid-spectrum gap, plays a crucial role in determining whether these Stokes points proliferate under some continuous change of parameters.

ACKNOWLEDGMENTS

We wish to thank Alexander Altland, Satya Majumdar, and Gregory Schehr for a useful suggestions. J.K.-D. and J.K. are supported by the Simons Foundation Grant No. 454943.

[1] M. Mézard and A. Montanari, *Information, Physics, and Computation*, Oxford Graduate Texts (Oxford University Press, Oxford, 2009).

[2] A. Crisanti and H.-J. Sommers, The spherical p-spin interaction spin glass model: the statics, *Zeitschrift für Physik B Condensed Matter* **87**, 341 (1992).

[3] T. Castellani and A. Cavagna, Spin-glass theory for pedestrians, *J. Stat. Mech.* (2005) P05012.

[4] D. Cartwright and B. Sturmfels, The number of eigenvalues of a tensor, *Linear Algebra Appl.* **438**, 942 (2013).

[5] A. Auffinger, G. Ben Arous, and J. Černý, Random matrices and complexity of spin glasses, *Commun. Pure Appl. Math.* **66**, 165 (2012).

[6] A. Auffinger and G. Ben Arous, Complexity of random smooth functions on the high-dimensional sphere, *Ann. Probab.* **41**, 4214 (2013).

[7] M. Kac, On the average number of real roots of a random algebraic equation, *Bull. Am. Math. Soc.* **49**, 314 (1943).

[8] S. O. Rice, The distribution of the maxima of a random curve, *Am. J. Math.* **61**, 409 (1939).

[9] Y. V. Fyodorov, Complexity of Random Energy Landscapes, Glass Transition, and Absolute Value of the Spectral Determinant of Random Matrices, *Phys. Rev. Lett.* **92**, 240601 (2004).

[10] A. Crisanti and H.-J. Sommers, Thouless-Anderson-Palmer approach to the spherical p-spin spin glass model, *J. Phys. I France* **5**, 805 (1995).

[11] L. F. Cugliandolo and J. Kurchan, Analytical Solution of the Off-Equilibrium Dynamics of a Long-Range Spin-Glass Model, *Phys. Rev. Lett.* **71**, 173 (1993).

[12] F. Antenucci, A. Crisanti, and L. Leuzzi, Complex spherical 2 + 4 spin glass: A model for nonlinear optics in random media, *Phys. Rev. A* **91**, 053816 (2015).

[13] D. Anninos, T. Anous, and F. Denef, Disordered quivers and cold horizons, *J. High Energy Phys.* **12** (2016) 071.

[14] E. Witten, A new look at the path integral of quantum mechanics, *Surv. Differ. Geom.* **15**, 345 (2010).

[15] E. Witten, Analytic continuation of Chern-Simons theory, in *Chern-Simons Gauge Theory: 20 Years After*, AMS/IP Studies in Advanced Mathematics, edited by J. E. Andersen, H. U. Boden, A. Hahn, and B. Himpel (American Mathematical Society, Providence, 2011), Vol. 50, pp. 347–446.

[16] M. Cristoforetti, F. Di Renzo, and L. Scorzato, New approach to the sign problem in quantum field theories: High density QCD on a Lefschetz thimble, *Phys. Rev. D* **86**, 074506 (2012).

[17] Y. Tanizaki, H. Nishimura, and J. J. M. Verbaarschot, Gradient flows without blow-up for Lefschetz thimbles, *J. High Energy Phys.* **10** (2017) 100.

[18] L. Scorzato, The Lefschetz thimble and the sign problem, *PoS (LATTICE 2015)*, 016 (2016).

[19] E. Bogomolny, O. Bohigas, and P. Leboeuf, Distribution of Roots of Random Polynomials, *Phys. Rev. Lett.* **68**, 2726 (1992).

[20] A. J. Bray and D. S. Dean, Statistics of Critical Points of Gaussian Fields on Large-Dimensional Spaces, *Phys. Rev. Lett.* **98**, 150201 (2007).

[21] H. H. Nguyen and S. O’Rourke, The elliptic law, *Int. Math. Res. Not.* **2015**, 7620 (2015).

[22] G. Livan, M. Novaes, and P. Vivo, *Introduction to Random Matrices*, Springer Briefs in Mathematical Physics, Vol. 26 (Springer International, New York, 2018).

[23] H. Weyl, Das asymptotische verteilungsgesetz der eigenwerte linearer partieller differentialgleichungen (mit einer anwendung auf die theorie der hohlraumstrahlung), *Math. Ann.* **71**, 441 (1912).

[24] E. Bézout, *Théorie générale des équations algébriques* (L’imprimerie de Ph.-D. Pierres, 1779).

[25] J. Manschot, B. Pioline, and A. Sen, From black holes to quivers, *J. High Energy Phys.* **11** (2012) 023.

[26] A. Behtash, G. V. Dunne, T. Schäfer, T. Sulejmanpasic, and M. Ünsal, Toward Picard–Lefschetz theory of path integrals, complex saddles and resurgence, *Ann. Math. Sci. Appl.* **2**, 95 (2017).