Superposition of random plane waves in high spatial dimensions: Random matrix approach to landscape complexity

Lacroix-A-Chez-Toine, Bertrand; Fyodorov, Yan V.; Fedeli, Sirio Belga

doi:10.1063/5.0086919

Motivated by current interest in understanding statistical properties of random landscapes in high-dimensional spaces, we consider a model of the landscape in $R^{N}$ obtained by superimposing M > N plane waves of random wavevectors and amplitudes and further restricted by a uniform parabolic confinement in all directions. For this landscape, we show how to compute the “annealed complexity,” controlling the asymptotic growth rate of the mean number of stationary points as N → ∞ at fixed ratio α = M/N > 1. The framework of this computation requires us to study spectral properties of N × N matrices W = KTK^T, where T is a diagonal matrix with M mean zero independent and identically distributed (i.i.d.) real normally distributed entries, and all MN entries of K are also i.i.d. real normal random variables. We suggest to call the latter Gaussian Marchenko–Pastur ensemble as such matrices appeared in the seminal 1967 paper by those authors. We compute the associated mean spectral density and evaluate some moments and correlation functions involving products of characteristic polynomials for such matrices.

I. INTRODUCTION

The problem of characterizing statistical properties of high-dimensional random functions, frequently called ”random landscapes,” originated in the theory of disordered systems, such as spin glasses (see Ref. 1 for an accessible introduction), and gradually became popular beyond the original setting, finding numerous applications in diverse fields, such as cosmology,^2,3 machine learning via deep neural networks,^4,5 and large-size inference problems in statistics.^6–9 One of the simplest yet nontrivial landscape characteristics are the so-called landscape complexities given by the rates of exponential growth of the number of stationary points (minima, maxima, and saddles of a particular index) with the number N of dimensions of the associated landscape. Extracting those rates for various models of random spin-glass landscapes initially attracted attention in theoretical physics literature^10–14 and became a focus of vigorous activity in recent years; see Refs. 15–24 and the references therein.

The fully controlled (and eventually rigorous) evaluation of complexities has been so far demonstrated only for a few types of landscapes whose random parts possess a high level of statistical symmetry in the underlying N-dimensional space. One of the most paradigmatic examples is the landscape

V (x) = \frac{μ}{2} | x |^{2} + V (x), x \in R^{N},

(1)

where μ is the strength of a isotropic harmonic confinement and V(x) is a Gaussian-distributed random function, which is translationally and rotationally invariant and characterized by the covariance structure,

⟨V (x_{1}) V (x_{2})⟩ = N Γ (\frac{{(x_{1} - x_{2})}^{2}}{2 N}) .

(2)

In such a case, the associated Hessian $W_{i j} (x) = - \partial_{x_{i}, x_{j}} V (x)$ turns out to be given by a matrix simply related^12,25 to the classical Gaussian orthogonal ensemble of random matrices. Similar Hessians also appear in studies of rotationally invariant random functions confined to the surface of a high-dimensional sphere.^15,16 In both situations, the exploitation of the so-called Kac–Rice formula allowed to perform explicit computations of the so-called “annealed” (i.e., related to the mean number of stationary points) complexities. Moreover, one can further show that annealed and quenched (i.e., characterizing typical realizations of the landscape) complexities coincide for a certain class of potentials V(x)^18,22 though, in general, they may differ.

Studying this total annealed complexity revealed that the invariant landscapes (1) and (2) undergo a “topology trivialization transition” from the glassy phase with exponentially many stationary points for small values of the ratio $μ / \sqrt{Γ^{''} (0)} < 1$ to the phase where there is typically only a single minimum for $μ / \sqrt{Γ^{''} (0)} > 1$ ⁠. Similar transitions are also operative in random potentials confined to high-dimensional spheres.^7,25–27

The goal of this article is to suggest another class of random high-dimensional landscapes amenable to well-controlled treatment of the annealed complexity. Our construction starts with the same representation (1) but proceeds with choosing the random part in the form of a superposition of random plane waves,

V (x) = \sum_{a = 1}^{M} [u_{1}^{(a)} \cos (k_{a} \cdot x) + u_{2}^{(a)} \sin (k_{a} \cdot x)],

(3)

where $u_{1, 2}^{(a)}$ are 2M independent and identically distributed (i.i.d.) normal (mean zero, variance one) random variables.

In Eq. (3), M i.i.d. random vectors k_a are drawn as either

random vectors of real Gaussian random variables with zero average and variance 1/N or
random vectors on the unit N-dimensional sphere. Such a vector k_a can be built as $k_{i a} = {(O_{a})}_{i i}$ for i = 1, …, N, where the N × N random matrix O_a is drawn from the orthogonal group.

We will show that both cases lead to the same asymptotic result in the limit N → ∞. We, henceforth, denote α = M/N. Note that the random functions defined via (3) are not Gaussian-distributed due to additional randomness in the choice of wavevectors k_a. Note also that for α < 1 at every realization, the function V(x) is constant in the M − N subspace orthogonal to the one spanned by M random wavevectors k_a. To avoid this situation and to have V(x) varying in the whole space, we, therefore, restrict our computation of the landscape complexity to the case α ≥ 1.

It is also worth mentioning that superpositions of random plane waves naturally emerge in the context of quantum chaos through Berry’s conjecture²⁸ and are currently under intensive study as a universal model for high-energy eigenfunctions of the Laplace operator on generic compact Riemannian manifolds with fixed eigenvalue E = k². Stationary points in that context have intensively been studied in low spatial dimension N = 2; see Ref. 29 and the references therein. In contrast, we are interested in the case of asymptotically high dimension N ≫ 1.

Finally, for low values of M, N, model (3) describes periodic systems, such as crystals in the presence of quenched impurities, and was much studied mostly for its thermodynamics and dynamics, especially in the context of vortex lattices in superconductors; for reviews, see, e.g., Refs. 30 and 31.

We believe that the suggested landscape in high dimensions represents a new interesting type worth consideration. In particular, we have the following:

Such landscapes still allow for an explicit analytical computation of the total annealed complexity beyond the class of universality of previously studied potentials satisfying (2). We will see that the arising μ-dependence will be quite different. In particular, the landscape will not show the sharp topology trivialization transition at any finite value of μ.
In contrast to translationally invariant random fields, which are extremely difficult to simulate even for moderately big N ∼ 10, a superposition of random plane waves is easily constructed numerically.

As will be demonstrated below, the mean total number of stationary points in the landscape can be expressed in terms of the random matrix average,

⟨N_{tot} (μ; α)⟩ = \frac{⟨| \det (μ I - W) |⟩}{μ^{N}}, μ > 0,

(4)

where the Hessian matrix W can be represented as W = KTK^T. For M i.i.d. Gaussian vectors k_a’s, this matrix turns out to belong to the “Gaussian Marchenko–Pastur ensemble (GMPE), which we define as

W_{i j} = \sum_{a = 1}^{M} T_{a} k_{i a} k_{j a}, 1 \leq i, j \leq N,

(5)

where the variables k_ia for 1 ≤ i ≤ N and 1 ≤ a ≤ M (considered as the entries of matrix K) are independent real normal random variables (mean zero, variance 1/N), whereas the M × M matrix T is diagonal with N(0, 1) normal random elements T_a, a = 1, …, M. Here and in the following, we denote ${⟨\dots⟩}_{T}$ as the average with respect to the diagonal matrix T, ${⟨\dots⟩}_{K}$ as the average with respect to the matrix K, and $⟨\dots⟩$ as the average with respect to both K and T. The suggested name for this random matrix ensemble³² seems appropriate to us as these types of matrices appeared already in the seminal Marchenko–Pastur paper,³³ although they do not seem to be under much study ever since.

Note that had we replaced the diagonal matrix T with the identity matrix, one would arrive at the standard Wishart ensemble whose eigenvalue density at N → ∞ is given by the famous Marchenko–Pastur distribution. Taking any positive-definite T would provide a natural generalization of the Wishart ensemble. However, the sign-indefinite nature of T in our case makes the properties of the resulting matrices W very different from the Wishart case. We believe that the resulting ensemble GMPE is interesting in its own right, and we further study its mean eigenvalue density as N → ∞ and correlation properties of its characteristic polynomial for both fixed and random diagonal matrix T. While the main focus of this article is on the landscape complexity, hence on the real symmetric matrices with the Dyson index β = 1 (corresponding to K being a matrix with real independent Gaussian elements), from the point of view of Random Matrix Theory (RMT), it is natural to extend the definition of the Gaussian Marchenko–Pastur ensemble to the Dyson index β = 2 by considering K with normal complex independent elements. (One can also consider taking K with quaternion independent normally distributed elements to generate the ensemble with Dyson index β = 4.) The p-point correlation function of the associated characteristic polynomial for a fixed diagonal matrix T is, then, defined as

Z_{β, p} (λ; T) = {⟨\prod_{j = 1}^{p} \det (λ_{j} I - W)⟩}_{K},

(6)

where we take the vector λ = (λ₁, …, λ_p) of spectral parameters to have real elements. For random T, the correlation function in (6) is defined with additionally taking the expectation with respect to the distribution of T_a’s. We will provide some explicit results for correlation functions for fixed T at coinciding points λ₁ = ⋯ = λ_p for any integer p and N, whereas for random T, we will restrict ourselves to p = 1, 2. Treating in the latter case for non-coinciding points λ₁ ≠ λ₂, we demonstrate that the correlation structure emerging after averaging over random T_a is quite different from the one known for both the standard invariant ensembles of random matrices and for ensembles with i.i.d. entries. The correlations of characteristic polynomials in the latter two types of RMT ensembles have been under intensive studies in the last two decades; see Refs. 34–45.

The rest of this paper is organized as follows. In Sec. II, we summarize the main results obtained in this paper. In Sec. III, we consider the total annealed complexity of landscape (1) with the random potential defined in (3). In Secs. IV–V, we characterize the Gaussian Marchenko–Pastur ensemble associated with this counting problem by computing its mean spectral density in Sec. IV and the correlation function of its characteristic polynomial defined in (6) in Sec. V. Finally, in Sec. VI, we discuss and interpret our results and give some perspective on the further investigation of this random landscape. Some details of the computations are relegated to Appendixes A–C.

II. MAIN RESULTS

A. Results on the counting problem

Our main result is that in the limit M, N → ∞ with fixed ratio α = M/N > 1, the total annealed complexity for landscape (1) with the random potential (3) can be expressed as

Ξ_{tot} (μ; α) = \lim_{N \to \infty} \frac{1}{N} \ln ⟨N_{tot} (μ; α)⟩ = \int_{μ}^{\infty} (\frac{1}{ν} - m_{r} (ν)) d ν

(7)

in terms of a function m_r(μ) ≔ m_r, which together with its counterpart m_i(μ) ≔ m_i satisfies the following pair of equations:

μ \int_{- \infty}^{\infty} d t e^{- \frac{t^{2}}{2}} \sqrt{t^{2} m_{i}^{2} + {(1 - t m_{r})}^{2}} = α \int_{- \infty}^{\infty} d t \frac{t e^{- \frac{t^{2}}{2}}}{\sqrt{t^{2} m_{i}^{2} + {(1 - t m_{r})}^{2}}},

(8)

\int_{- \infty}^{\infty} d t e^{- \frac{t^{2}}{2}} \sqrt{t^{2} m_{i}^{2} + {(1 - t m_{r})}^{2}} = α (m_{i}^{2} + m_{r}^{2}) \int_{- \infty}^{\infty} d t \frac{t^{2} e^{- \frac{t^{2}}{2}}}{\sqrt{t^{2} m_{i}^{2} + {(1 - t m_{r})}^{2}}} .

(9)

Using this expression, we show that in sharp contrast to previously studied models with isotropic harmonic confinement, the total annealed complexity does not undergo a “topology trivialization transition” and stays positive for any finite value 0 < μ < ∞ of the strength of the confinement. Still, it decreases very rapidly as μ → ∞,

Ξ_{tot} (μ; α) \approx \frac{α}{μ^{3}} \sqrt{\frac{2}{π}} e^{- \frac{μ^{2}}{2}}, μ \to \infty .

(10)

In Fig. 1, we show a plot of our analytical prediction for the annealed total complexity in the limit N → ∞ together with numerical simulation for finite N, showing good agreement.

FIG. 1.

View large Download slide

Comparison between the annealed total complexity Ξ_tot(μ; α) for α = 2 (black solid line) obtained by inserting the solution of Eqs. (8) and (9) into Eq. (7) and numerical simulations of $N^{- 1} \ln ⟨| \det (I - μ^{- 1} W) |⟩$ (blue circle) plotted as a function of μ. It is expected that this quantity converges to the annealed complexity in the N → ∞ limit. These quantities have good agreement for low values of μ, but the numerical data seem to diverge away from its theoretical prediction for larger values of μ. This discrepancy results both from finite size effects and from the fact that the average is dominated by rare events, preventing an efficient numerical evaluation.

It is worth noting that methodologically the main difference with previously studied cases of Eq. (2) is that the annealed complexity is obtained using the functional variational principle involving the mean spectral density n(t) of the eigenvalues of the random matrix T instead of a variation principle over a single real random variable, as in Ref. 12 and related studies.

B. Results on the general Marchenko–Pastur ensemble

In this section, we give the main outcomes of our studies of spectral characteristics of general matrices of the form W = KTK^T related to the Marchenko–Pastur ensemble, not only those directly involved in the computation of annealed complexity.

1. Mean asymptotic eigenvalue density for GMPE

Consider the mean spectral density $ρ (λ) = \lim_{N \to \infty} \frac{1}{N} \sum_{k = 1}^{N} ⟨δ (λ - λ_{k})⟩$ for the Gaussian MPE matrices W defined in (5). It turns out to be convenient to define the associated resolvent/Stiltjes transform, defined here for a real spectral parameter λ as

(11)

where m_r ≡ m_r(λ) and m_i ≡ m_i(λ) are, respectively, its real and imaginary parts. The mean spectral density can be expressed in terms of the resolvent m(λ) as $ρ (λ) = \frac{I [m (λ)]}{π}$ where for GMPE, the resolvent m(λ) is obtained by solving the following transcendental equation:

λ = \frac{1 - α}{m (λ)} + \frac{i α}{m^{2} (λ)} Ψ (\frac{1}{i m (λ)}),

(12)

with $Ψ (u) = \sqrt{\frac{π}{2}} e^{\frac{u^{2}}{2}} (1 + erf (u / \sqrt{2}))$ ⁠. The resulting mean spectral density can be shown to be even: ρ(λ) = ρ(−λ), and supported on the whole real line, with the asymptotic behavior as λ → ±∞ given by

ρ (λ) \approx \frac{α}{\sqrt{2 π}} e^{- \frac{λ^{2}}{2}}, λ \to \pm \infty .

(13)

For α > 1, the density ρ(λ) takes its finite maximum value at the origin λ = 0, where ρ_α ≔ ρ(λ = 0) is the solution of the transcendental equation obtained by inserting m(0) = iπρ_α in Eq. (12). For α = 1, the density diverges close to the origin as $ρ (λ) \approx 1 / \sqrt{{(2 π)}^{3 / 2} | λ |}$ ⁠.

Let us finally mention that for α < 1, the density ρ(λ) displays a delta peak at the origin (1 − α)δ(λ), reflecting the presence of M − N zero eigenvalues. At the same time, its continuous part tends to a finite value at the origin: $α / (\sqrt{2 π} (1 - α))$ ⁠. In Fig. 2, we plot a comparison between our analytical prediction for the average density of eigenvalues ρ(λ) in the N → ∞ limit and results from numerical simulation for N = 500, showing excellent agreement.

FIG. 2.

View large Download slide

Comparison between the average density ρ(λ) for α = 1, 5 obtained analytically in the limit N → ∞ from Eq. (12) by using that $ρ (λ) = \frac{I [m (λ)]}{π}$ (respectively, blue solid and orange dashed lines) and numerically for a matrix of finite size N = 500 (respectively, blue circles and orange triangles).

2. Moments and correlation functions of characteristic polynomials

For a fixed matrix T, the positive moments of the characteristic polynomial $Z_{β, p} (λ; T) \equiv Z_{β, p} (λ 1; T)$ defined in Eq. (6) with setting λ = λ1 = (λ, …, λ) can be expressed for any integer p > 0 and any size N in terms of a p × p determinant for β = 2,

Z_{β = 2, p} (λ; T) = \frac{\det {[g_{N, k - j} (λ; T)]}_{j, k = 1}^{p}}{\det {[q_{N, k - j}]}_{j, k = 1}^{p}},

(14)

and in terms of a 2p × 2p Pfaffian for β = 1,

Z_{β = 1, p} (λ; T) = \frac{Pf {[(j - k) g_{N, 2 p + 1 - (k + j)} (λ; T)]}_{j, k = 1}^{2 p}}{Pf {[(j - k) q_{N, 2 p + 1 - (k + j)}]}_{j, k = 1}^{2 p}} .

(15)

In these expressions, the functions g_N,m(λ; T) and q_N,m are defined as

g_{N, m} (λ; T) = λ^{N + m} \sum_{l = 0}^{N + m} \frac{{(- 1)}^{l}}{(N + m - l)!} e_{l} (\frac{T}{N λ}), q_{N, m} = g_{N, m} (1; O) = \frac{1}{(N + m)!},

(16)

where $O$ is the matrix whose elements are all zero and $e_{l} (X_{1}, \dots, X_{n}) = \sum_{1 \leq i_{1} < i_{2} < \dots < i_{l} \leq n} \prod_{j = 1}^{l} X_{i_{j}}$ is the lth elementary symmetric polynomial. Note that taking $T = I$ ⁠, one recovers the known results for the Wishart–Laguerre ensemble; see, e.g., Ref. 43.

For a random diagonal matrix T (not necessarily normally distributed), we only managed to consider the mean of the characteristic polynomials and its two-point correlation function, p = 2. We find that the mean characteristic polynomial only depends on the first moment $⟨T⟩ = \int d t t n (t)$ ⁠, where $n (t) ≔ \lim_{M \to \infty} \frac{1}{M} \sum_{a = 1}^{M} δ (t - T_{a})$ is the associated density function. For non-zero $⟨T⟩ \neq 0$ ⁠,

Z_{β, 1} (λ) = {⟨Z_{β, 1} (λ; T)⟩}_{T} = {(\frac{⟨T⟩}{N})}^{N} \frac{M!}{(M - N)!} {}_{1}F_{1} (- N; 1 + M - N; \frac{λ N}{⟨T⟩})

(17)

in terms of the hypergeometric function defined as ${}_{p}F_{q} (a; b; x) = \sum_{n = 0}^{\infty} \frac{\prod_{i = 1}^{p} {(a_{i})}_{n}}{\prod_{j = 1}^{q} {(b_{j})}_{n} n!} x^{n}$ ⁠. For $⟨T⟩ = 0$ ⁠, we instead find $Z_{β, 1} (λ) = λ^{N}$ ⁠. Note that in any case, $Z_{β, 1} (λ)$ is independent of β. In the large N limit, we study the rate of growth of the mean characteristic polynomial and show that

Ξ_{1} (x; α) = \lim_{N \to \infty} \frac{1}{N} \ln |\frac{Z_{β, 1} (λ)}{λ^{N}}| = \{\begin{cases} L_{1} (w_{+}; x, α), & x < x_{-} (α), \\ R [L_{1} (w_{\pm}; x, α)], & x_{-} (α) \leq x \leq x_{+} (α), \\ L_{1} (w_{-}; x, α), & x > x_{+} (α), \end{cases})

(18)

where we denoted $x = λ / ⟨T⟩$ and $x_{\pm} (α) = {(1 \pm \sqrt{α})}^{2}$ and introduced

L_{1} (w; x, α) = x w - \ln x w + α \ln (1 - w) - 1,

(19)

w_{\pm} (x, α) = \frac{1 - α + x \pm \sqrt{{(α - 1 - x)}^{2} - 4 x}}{2 x} .

(20)

The two-point spectral correlation function of the characteristic polynomials has been computed for any size N in the case of vanishing mean $⟨T⟩ = 0$ ⁠. It, then, turns out to depend on T only through its variance/second moment $⟨T^{2}⟩ = \int d t t^{2} n (t)$ and reads

\begin{align} Z_{β, 2} (λ = (λ_{1}, λ_{2})) = {(\frac{⟨T^{2}⟩}{N^{2}})}^{N} \frac{β}{2} Γ (N + \frac{2}{β} + 1) \frac{M!}{(M - N)!} {}_{1}F_{2} (- N; \frac{2}{β} + 1, 1 + M - N; - \frac{λ_{1} λ_{2} N^{2}}{⟨T^{2}⟩}) . \end{align}

(21)

In particular, in the large N limit, we define

Ξ_{2} (u; α) = \lim_{N \to \infty} \frac{1}{N} \ln |\frac{Z_{β, 2} (λ)}{{(λ_{1} λ_{2})}^{N}}| = L_{2} (η_{*}; u, α)

(22)

and find that it only depends on the rescaled positive parameter $u = | λ_{1} λ_{2} | / ⟨T^{2}⟩ \geq 0$ ⁠,

L_{2} (η; u, α) = α \ln (α) - (α + η - 1) \ln (α + η - 1)

(23)

+ η \ln (u (1 - η)) + 2 η - 2 η \ln (η) - \ln (1 - η) - 2 - \ln u,

\begin{align} η_{*} = \frac{1 - α}{3} & - \frac{\sqrt[3]{2} Θ}{3 \sqrt[3]{Δ + \sqrt{Δ^{2} + 4 Θ^{3}}}} + \frac{\sqrt[3]{Δ + \sqrt{Δ^{2} + 4 Θ^{3}}}}{3 \sqrt[3]{2}}, \\ w h e r e Δ & = 9 u (2 + α) - 2 {(α - 1)}^{3}, Θ = 3 u - {(α - 1)}^{2} . \end{align}

(24)

Note a somewhat surprising feature of the two-point function in Marchenko–Pastur ensembles with random matrix T: it only depends on the product λ₁λ₂, in sharp contrast to earlier studied cases³⁴ where in the large-N limit, it rather depended on the spectral difference |λ₁ − λ₂|.

III. SOLUTION OF THE ANNEALED COUNTING PROBLEM

In this section, we show how to solve the annealed counting problem discussed in the Introduction. The starting point of our consideration is the standard Kac–Rice formula for the total number of equilibria, which in case (1) reads

⟨N_{tot} (μ; α)⟩ = \int_{R^{N}} d^{N} x ⟨\prod_{i = 1}^{N} δ (μ x_{i} + \partial_{x_{i}} V (x)) | \det (μ I - W (x)) |⟩,

(25)

where $W_{i j} (x) = - \partial_{x_{i}, x_{j}} V (x)$ and the expectation $⟨\dots⟩$ in the case of potential (3) has been taken over both the components of the vector k_a and the coefficients $u_{1, 2}^{(a)}$ for a = 1, …, M. It turns out to be expedient to pass from the original random variables $u_{1, 2}^{(a)}, a = 1, \dots, M$ ⁠, to the new set of random variables obtained by rotation as

T_{a} = u_{1}^{(a)} \cos (k_{a}^{T} \cdot x) + u_{2}^{(a)} \sin (k_{a}^{T} \cdot x),

(26)

G_{a} = u_{1}^{(a)} \sin (k_{a}^{T} \cdot x) - u_{2}^{(a)} \cos (k_{a}^{T} \cdot x) .

(27)

Despite the rotation used to define T_a and G_a being explicitly x and k_a dependent, the resulting random variables are statistically x and k_a independent with zero mean and covariance

⟨T_{a} T_{b}⟩ = ⟨G_{a} G_{b}⟩ = δ_{a, b}, ⟨T_{a} G_{b}⟩ = 0, a, b = 1, \dots, M .

(28)

The random force $f_{i} = - \partial_{x_{i}} V (x)$ and the Hessian $W_{i j} = - \partial_{x_{i}, x_{j}} V (x)$ associated with the potential V(x) are conveniently expressed in terms of these random variables as

f_{i} \equiv - \partial_{x_{i}} V (x) = \sum_{a = 1}^{M} k_{a, i} G_{a},

(29)

W_{i j} \equiv - \partial_{x_{i}, x_{j}} V (x) = \sum_{a = 1}^{M} k_{a, i} k_{a, j} T_{a},

(30)

and using the properties of G_a, T_a, and K, it is easy to see that statistical properties of these random variables are also independent of the position x. Now, we may use in Eq. (25) the Fourier integral representation for the Dirac delta function,

δ (μ x_{j} - f_{j}) = \int_{- \infty}^{\infty} \frac{d q_{j}}{2 π} e^{i q_{j} (μ x_{j} - f_{j})} = \int_{- \infty}^{\infty} \frac{d q_{j}}{2 π} e^{i q_{j} (μ x_{j} - \sum_{a = 1}^{M} k_{a, j} G_{a})},

(31)

which allows us to take explicitly the expectation with respect to random variables G_a’s, yielding

⟨N_{tot} (μ; α)⟩ = \int_{R^{N}} d^{N} x \int_{R^{N}} \frac{d^{N} q}{{(2 π)}^{N}} e^{i μ \sum_{j = 1}^{N} q_{j} x_{j}} {⟨e^{- \sum_{a = 1}^{M} \frac{{(k_{a}^{T} \cdot q)}^{2}}{2}} | \det (μ I - K T K^{T}) |⟩}_{K, T} .

(32)

Further integration over x yields

\begin{align} ⟨N_{tot} (μ; α)⟩ & = \int_{R^{N}} d^{N} q \prod_{j = 1}^{N} δ (μ q_{j}) {⟨e^{- \sum_{a = 1}^{M} \frac{{(k_{a}^{T} \cdot q)}^{2}}{2}} | \det (μ I - K T K^{T}) |⟩}_{K, T} \\ = \frac{{⟨| \det (μ I - K T K^{T}) |⟩}_{K, T}}{μ^{N}}, \end{align}

(33)

and after writing explicitly the expectation over the diagonal matrix T, the total number of equilibria finally reads

⟨N_{tot} (μ; α)⟩ = \frac{1}{μ^{N}} \int_{R^{M}} \prod_{a = 1}^{M} \frac{d T_{a}}{\sqrt{2 π}} e^{- \sum_{a = 1}^{M} \frac{T_{a}^{2}}{2}} {⟨| \det (μ I - K T K^{T}) |⟩}_{K},

(34)

where ${⟨\dots⟩}_{K}$ is the remaining expectation over the random matrix K.

We will now focus on the limit N, M → ∞ with 1 ≤ M/N = α < ∞. In this particular limit, we define the annealed total complexity as

Ξ_{tot} (μ; α) = \lim_{N \to \infty} \frac{1}{N} \ln ⟨N_{tot} (μ; α)⟩ = \lim_{N \to \infty} \frac{1}{N} \ln \frac{⟨| \det (μ I - W) |⟩}{μ^{N}},

(35)

where we remind that W = KTK^T. In order to make progress in this computation, we assume that for a fixed diagonal matrix T holds the strong self-averaging property, i.e.,

\lim_{N \to \infty} \frac{1}{N} \ln {⟨| \det (μ I - K T K^{T}) |⟩}_{K} = \lim_{N \to \infty} \frac{1}{N} Tr {⟨\ln | μ I - K T K^{T} |⟩}_{K} = \int d λ ρ (λ) \ln | λ - μ |,

(36)

where the mean limiting spectral density ρ(λ) is defined via

ρ (λ) = \lim_{N \to \infty} \frac{1}{N} \sum_{i = 1}^{N} {⟨δ (λ - λ_{i})⟩}_{K},

(37)

with λ_i’s being the random eigenvalues of W = KTK^T. This density obviously depends on the diagonal matrix T. Note that in the model with the Gaussian random potential satisfying Eq. (2), a similar strong self-averaging property has been proven rigorously^15,16 in considerable generality. With the assumption of its validity for the present model, we may express the average total number of equilibria/stationary points for finite but large N as

⟨N_{tot} (μ; α)⟩ \approx \frac{1}{μ^{N}} \int_{R^{M}} \prod_{a = 1}^{M} \frac{d T_{a}}{\sqrt{2 π}} e^{N (- \frac{α}{M} \sum_{a = 1}^{M} \frac{T_{a}^{2}}{2} + \int_{- \infty}^{\infty} d λ ρ (λ) \ln | λ - μ |)} .

(38)

To facilitate extracting the large-N asymptotics, we find it expedient to introduce a variational problem over the normalized limiting density of T_a as

n (t) = \lim_{M \to \infty} \frac{1}{M} \sum_{a = 1}^{M} δ (T_{a} - t) .

(39)

To this end, we replace the M-dimensional integral over components T_a’s by a functional integral over this density as

⟨N_{tot} (μ; α)⟩ \approx \frac{1}{μ^{N} {(2 π)}^{M / 2}} \int \frac{d l}{2 π} \int D n (t) e^{- N S_{tot} [n; μ, α, l]},

(40)

\begin{align} S_{tot} [n; μ, α, l] & = α [\frac{1}{2} \int_{- \infty}^{\infty} d t n (t) t^{2} + \int_{- \infty}^{\infty} d t n (t) \ln n (t) + l (\int_{- \infty}^{\infty} d t n (t) - 1)] \\ - \int_{- \infty}^{\infty} d λ ρ (λ) \ln | μ - λ |, \end{align}

(41)

where l is the Lagrange multiplier necessary to ensure the correct normalization of the density. In the expression of S_tot, the second term corresponds to an entropic contribution resulting from the change of variables from the M variables T_a’s to their density n(t) (see Refs. 14 and 46 for a detailed derivation). This form is now amenable to using the Laplace method justified by the large parameter N, yielding for the annealed complexity,

Ξ_{tot} (μ; α) = \lim_{N \to \infty} \frac{1}{N} \ln ⟨N_{tot} (μ; α)⟩ = - \min_{l, n (t)} S_{tot} [n; μ, α, l] - \ln μ - α \ln (2 π) .

(42)

Let us denote n_tot and l_tot as the values of the limiting density and Lagrange parameters that minimize S_tot. Note that the whole procedure is not unlike to what has been done when evaluating the complexity for Gaussian translationally and rotationally invariant model in Ref. 12; however, instead of maximizing with respect to a single real variable, one has to optimize with respect to the limiting spectral density n(t) in the present case. We also note that although the optimization is presented here in an informal style of theoretical physics, we believe that it can be recast rigorously using the large deviations framework.

Varying the action functional S_tot with respect to n(t) yields the stationarity condition

{(\frac{δ S_{tot}}{δ n (t)}|}_{n = n_{tot}} = α [\frac{t^{2}}{2} + l_{tot} + 1 + \ln n_{tot} (t)] - \int_{- \infty}^{\infty} d λ {(\frac{δ ρ (λ)}{δ n (t)}|}_{n = n_{tot}} \ln | μ - λ | = 0,

(43)

while the differentiation with respect to the Lagrange parameter l ensures again the normalization for n_tot(t). To find n_tot(t) and the associated action value, one needs to compute explicitly the functional derivative of ρ(λ) with respect to n(t). We will now show that for the specific type of matrix considered here, this problem is exactly solvable. The resolvent/Stiltjes transform function m(z) associated with the spectral density ρ(λ) satisfies for a fixed limiting density n(t) a self-consistency equation derived by Marchenko and Pastur in Ref. 33, which reads

m (z) = \int_{- \infty}^{\infty} d λ \frac{ρ (λ)}{z - λ} = {(z - α \int_{- \infty}^{\infty} d t \frac{t n (t)}{1 - t m (z)})}^{- 1} .

(44)

This equation holds in the asymptotic limit N, M → ∞ with fixed α = M/N = O(1) for a matrix K whose columns k_a for a = 1, …, M are either i.i.d. Gaussian vectors or i.i.d. vectors uniformly drawn on the unit N-sphere. The resolvent can be introduced in the saddle-point Eq. (43) using the identity

\begin{align} \int_{- \infty}^{\infty} d λ ρ (λ) (\ln | μ - λ | - \ln | μ |) & = \int_{\infty}^{μ} d ν [R (\int_{- \infty}^{\infty} d λ \frac{ρ (λ)}{ν - λ}) - \frac{1}{ν}] \\ = \int_{μ}^{\infty} d ν (\frac{1}{ν} - m_{r} (ν)), \end{align}

(45)

where m_r(ν) is the real part of the resolvent m(ν) at the real spectral parameter ν, and we have used that zm_r(z) = zdλρ(λ)/(z − λ) → 1 for any real or complex spectral parameter z in the limit |z| → ∞. Using the self-consistency equation [Eq. (44)], one can compute the functional derivative of m(z) at any complex spectral parameter z with respect to n(t) explicitly,

\frac{δ m}{δ n (t)} = {[1 - α m^{2} \int_{- \infty}^{\infty} d t^{'} \frac{{t^{'}}^{2} n (t^{'})}{{(1 - t^{'} m)}^{2}}]}^{- 1} \frac{α t m^{2}}{1 - t m} .

(46)

Comparing this expression with the derivative of the resolvent m(z) with respect to z,

\partial_{z} m = - {[1 - α m^{2} \int_{- \infty}^{\infty} d t^{'} \frac{{t^{'}}^{2} n (t^{'})}{{(1 - t^{'} m)}^{2}}]}^{- 1} m^{2},

(47)

one obtains a simple identity

\frac{δ m}{δ n (t)} = - \frac{α t}{1 - t m} \partial_{z} m = α \partial_{z} \ln [1 - t m] .

(48)

This identity can, then, be used to obtain

\begin{align} \frac{δ}{δ n (t)} \int_{- \infty}^{\infty} d λ ρ (λ) (\ln | μ - λ | - \ln μ) & = \int_{\infty}^{μ} d ν \frac{δ}{δ n (t)} [R (\int_{- \infty}^{\infty} d λ \frac{ρ (λ)}{ν - λ}) - \frac{1}{ν}] \\ = α R (\int_{\infty}^{μ} d ν \partial_{ν} \ln [1 - t m (z)]) \\ = \frac{α}{2} \ln [{(1 - t m_{r} (μ))}^{2} + t^{2} m_{i} {(μ)}^{2}], \end{align}

(49)

where we introduced m_r(μ) and m_i(μ) as the real and imaginary parts of the resolvent m(μ) at the real spectral parameter μ, respectively. Now, the stationarity condition (43) implies that the optimal density n_tot(t) satisfies the following equation:

\frac{t^{2}}{2} + l_{tot} + 1 + \ln n_{tot} = \frac{1}{2} \ln [{(1 - t m_{r} (μ))}^{2} + t^{2} m_{i} {(μ)}^{2}],

(50)

which is immediately solved in terms of the real and imaginary parts of the resolvent at position μ, yielding

n_{tot} (t) = \frac{e^{- \frac{t^{2}}{2}}}{\sqrt{2 π} Z (μ)} \sqrt{{(1 - t m_{r} (μ))}^{2} + t^{2} m_{i} {(μ)}^{2}},

(51)

with $Z (μ) = e^{l_{tot} + 1} / \sqrt{2 π}$ being the μ-dependent normalization factor given explicitly by

Z (μ) = \int_{- \infty}^{\infty} d t \frac{e^{- \frac{t^{2}}{2}}}{\sqrt{2 π}} \sqrt{{(1 - t m_{r} (μ))}^{2} + t^{2} m_{i} {(μ)}^{2}} .

(52)

Finally, at a real value of μ, the self-consistency equation Eq. (44) for the full resolvent m(μ) can be rewritten as a pair of equations for the real and imaginary parts m_r ≡ m_r(μ) and m_i ≡ m_i(μ),

\frac{m_{r}}{m_{i}^{2} + m_{r}^{2}} = μ - α \int_{- \infty}^{\infty} d t \frac{t n_{tot} (t) (1 - t m_{r})}{m_{i}^{2} t^{2} + {(1 - t m_{r})}^{2}},

(53)

\frac{m_{i}}{m_{i}^{2} + m_{r}^{2}} = α \int_{- \infty}^{\infty} d t \frac{t^{2} n_{tot} (t) m_{i}}{m_{i}^{2} t^{2} + {(1 - t m_{r})}^{2}} .

(54)

Recall that the imaginary part of the resolvent is expressed in terms of the limiting spectral density as m_i(μ) = πρ(μ). For the optimal density n_tot(t) of T being unbounded, it is expected that the limiting spectral density ρ(λ) of W = KTK^T is also unbounded; see Ref. 33. Hence, m_i(μ) > 0 for any 0 < μ < ∞, and one can divide by m_i on both sides of Eq. (54). Inserting, then, the expression of n_tot(t), one obtains the more explicit equations [Eqs. (8) and (9)]. This pair of equation allows us to express m_i(μ) and m_r(μ) at the optimal solution. One now needs to relate these values to the complexity expression (42) for Ξ_tot(μ; α). In order to obtain a convenient representation for the total complexity, we first compute its derivative with respect to μ as

(55)

where we have used that n_tot is the density that minimizes S_tot. Finally, using again that as μ → ∞, m_r(μ) = dλρ(λ)/(μ − λ) → μ⁻¹ such that m_r(μ) − μ⁻¹ → 0, one arrives at the final expression for the annealed complexity,

Ξ_{tot} (μ; α) = \int_{μ}^{\infty} d ν (\frac{1}{ν} - m_{r} (ν)) .

(56)

In the limit μ → ∞, we expect that m_r(μ) ∼ μ⁻¹, while $m_{i} (μ) \propto e^{- μ^{2} / 2}$ ⁠, which is inherited from the Gaussian decay of n_tot(t) (see also a related discussion on the asymptotic of spectral density below). In that limit, one can safely set m_i(μ) → 0 in Eqs. (52) and (53), which yields

\begin{align} Z (μ) & \approx \int_{- \infty}^{\infty} d t \frac{e^{- \frac{t^{2}}{2}}}{\sqrt{2 π}} | 1 - t m_{r} (μ) | = \sqrt{\frac{2}{π}} m_{r} e^{- \frac{1}{2 m_{r}^{2}}} + erf (\frac{1}{\sqrt{2} m_{r}}) \\ \approx 1 + \sqrt{\frac{2}{π}} m_{r}^{3} e^{- 1 / m_{r}^{2}}, μ \to \infty, \end{align}

(57)

\begin{align} \frac{1}{m_{r} (μ)} & \approx μ - α \int_{- \infty}^{\infty} d t \frac{t e^{- \frac{t^{2}}{2}}}{\sqrt{2 π} Z (μ)} sign (1 - t m_{r} (μ)), \\ = μ + α \sqrt{\frac{2}{π}} e^{- 1 / (2 m_{r}^{2})}, μ \to \infty . \end{align}

(58)

In particular, from the second equation, one obtains that

Ξ_{tot} (μ; α) = \int_{μ}^{\infty} d ν (\frac{1}{ν} - m_{r} (ν)) \approx α \sqrt{\frac{2}{π}} \int_{μ}^{\infty} \frac{d ν}{ν^{2}} e^{- ν^{2} / 2} \approx \frac{α}{μ^{3}} \sqrt{\frac{2}{π}} e^{- \frac{μ^{2}}{2}}, μ \to \infty,

(59)

showing that the total complexity is never zero for any finite strength of the isotropic confinement μ, although it vanishes very rapidly at large μ. Such behavior is in stark contrast with previously considered models with isotropic confinement, where the complexity was found to vanish for μ > μ_c, signaling of the total “landscape topology trivialization” beyond a finite value of the confining parameter.

IV. SPECTRAL DENSITY OF THE GAUSSIAN MARCHENKO–PASTUR ENSEMBLE

In this section, we study the simplest spectral characteristic of GMPE: its mean spectral density in the limit N, M → ∞ with α = M/N = O(1) defined as

ρ (λ) = \lim_{N \to \infty} \frac{1}{N} \sum_{i = 1}^{N} ⟨δ (λ - λ_{i})⟩,

(60)

where λ_i’s are the random eigenvalues of W = KTK^T. In this expression, the expectation is taken over both the random N × M matrix K with independent $N (0, 1 / \sqrt{N})$ entries and the random M × M diagonal matrix T with real i.i.d. N(0, 1) random elements. The resolvent m(z) = ∫dλρ(λ)/(z − λ), defined as the Stieltjes transform of the spectral density ρ(λ), is solution of (44) with the Gaussian limiting distribution $n (t) = \exp (- t^{2} / 2) / \sqrt{2 π}$ ⁠. Computing explicitly the integral over t, one obtains the following transcendental equation for m(z):

z = \frac{1 - α}{m (z)} + \frac{i α}{m^{2} (z)} Ψ (\frac{1}{i m (z)}),

(61)

where $Ψ (u) = \sqrt{\frac{π}{2}} e^{\frac{u^{2}}{2}} (1 + erf (u / \sqrt{2}))$ ⁠. Taking the parameter z = λ on the real axis, we introduce as in Sec. III the real and imaginary part of the resolvent m(λ) = m_r(λ) + im_i(λ). The spectral density is, then, extracted from the resolvent m(λ) by using the identity

ρ (λ) = \frac{d}{d λ} \lim_{y \to 0^{+}} \int_{λ_{0}}^{λ} \frac{d u}{π} m_{i} (u + i y) = \frac{m_{i} (λ)}{π} .

(62)

While Eq. (61) allows for an efficient numerical evaluation of the density at a given real position z = λ, the properties of the mean density ρ(λ) are most easily analyzed using the set of equations for the real and imaginary part of the resolvent,

\frac{m_{r}}{m_{i}^{2} + m_{r}^{2}} = λ - α \int_{- \infty}^{\infty} d t \frac{t n (t) (1 - t m_{r})}{m_{i}^{2} t^{2} + {(1 - t m_{r})}^{2}},

(63)

\frac{m_{i}}{m_{i}^{2} + m_{r}^{2}} = α \int_{- \infty}^{\infty} d t \frac{t^{2} n (t) m_{i}}{m_{i}^{2} t^{2} + {(1 - t m_{r})}^{2}} .

(64)

In particular, for an even limiting density n(−t) = n(t), one can show from (63) and (64) that the mean spectral density ρ(−λ) = ρ(λ) is even as well. Indeed, the pair of Eqs. (63) and (64) is invariant under the symmetry (λ, m_r, m_i) → (−λ, −m_r, m_i) together with the change of variable t → −t in the integrals. As ρ(λ) = m_i(λ)/π, the mean spectral density is also even.

For an unbounded limiting density n(t) of T, such as the Gaussian density in GMPE, the associated mean spectral density ρ(λ) for matrices W = KTK^T is also unbounded.³³ In that limit, one expects that m_i(λ) = πρ(λ) ∝ n(λ) ≪ m_r(λ) as λ → ∞. Using this approximation and introducing in the integrals of Eqs. (63) and (64) the change of variable $t = 1 / m_{r} + m_{i} u / m_{r}^{3}$ ⁠, one obtains to leading order in m_i,

\frac{1}{m_{r} (λ)} \approx λ + \frac{α}{m_{r}^{2}} n (\frac{1}{m_{r}}) \int_{- \infty}^{\infty} \frac{u d u}{u^{2} + m_{r}^{2}} = λ,

(65)

m_{i} (λ) \approx α n (\frac{1}{m_{r}}) \int_{- \infty}^{\infty} \frac{m_{r} d u}{m_{r}^{2} + u^{2}} \approx π α n (λ) .

(66)

In particular, for a Gaussian-shaped density n(t), it yields the asymptotic behavior

ρ (λ) = \frac{m_{i} (λ)}{π} \approx \frac{α}{\sqrt{2 π}} e^{- \frac{λ^{2}}{2}} for | λ | ≫ 1 .

(67)

For α > 1, the spectral density ρ(λ = 0) = ρ_α is finite, and in the present case, it is obtained by setting z = 0 and m(z) = iπρ_α in Eq. (61), yielding

ρ (λ = 0) = ρ_{α}, α > 1, (1 - \frac{1}{α}) \sqrt{2 π} ρ_{α} = e^{\frac{1}{2 π^{2} ρ_{α}^{2}}} erfc (\frac{1}{\sqrt{2} π ρ_{α}}) .

(68)

As α → 1₊, one obtains that the value $ρ_{α} \approx {[\sqrt{2 π} (1 - α)]}^{- 1}$ diverges.

For α < 1, instead, the mean spectral density displays both a continuous and a discontinuous part. The fraction of zero eigenvalues is given by 1 − M/N = 1 − α and is independent of the limiting density n(t) (as long as T does not have a macroscopic number of zero eigenvalues). On the other hand, the continuous part of the spectral density reaches a finite value ρ_α for λ → 0. The value of this constant can be obtained by inserting m(λ) = (1 − α)/λ + πρ_α in Eq. (61) and expanding for λ → 0. This yields

ρ_{α} = \lim_{λ \to 0} [ρ (λ) - (1 - α) δ (λ)] = \frac{α}{(1 - α) \sqrt{2 π}}, α < 1 .

(69)

Finally, in the special case α = 1, the spectral density is continuous but diverges on approaching the origin as

ρ (λ) \approx \frac{1}{{(2 π)}^{3 / 4} \sqrt{| λ |}}, λ \to 0, α = 1 .

(70)

Let us finally consider the moments $⟨λ^{n}⟩ = \int_{- \infty}^{\infty} d λ λ^{n} ρ (λ)$ ⁠. As ρ(λ) is even, it trivially yields all odd moments vanishing: $⟨λ^{2 n + 1}⟩ = 0$ ⁠. For any α, the even moments can be computed using that for z → ∞,

m (z) = \int_{- \infty}^{\infty} d λ \frac{ρ (λ)}{z - λ} \approx \frac{1}{z} \sum_{k = 0}^{\infty} \int_{- \infty}^{\infty} d λ {(\frac{λ}{z})}^{k} ρ (λ) = \sum_{k = 0}^{\infty} \frac{⟨λ^{k}⟩}{z^{k + 1}} = \sum_{n = 0}^{\infty} \frac{⟨λ^{2 n}⟩}{z^{2 n + 1}} .

(71)

Using Eq. (44), one obtains that

z - \frac{1}{m (z)} = α \int_{- \infty}^{\infty} d t \frac{t n (t)}{1 - t m (z)} \approx α \int_{- \infty}^{\infty} d t t n (t) \sum_{k = 0}^{\infty} {(t m (z))}^{k} \approx α \sum_{n = 1}^{\infty} ⟨T^{2 n}⟩ m {(z)}^{2 n - 1},

(72)

where $⟨T^{k}⟩ = \int_{- \infty}^{\infty} d t t^{k} n (t)$ ⁠, and we have considered a symmetric density profile n(−t) = n(t). In particular, for the Gaussian distribution, $⟨T^{2 n}⟩ = 2^{n} Γ (n + 1 / 2) / \sqrt{π}$ and the lowest moments can be obtained explicitly as

⟨λ^{2}⟩ = \int_{- \infty}^{\infty} d λ λ^{2} ρ (λ) = α,

(73)

⟨λ^{4}⟩ = \int_{- \infty}^{\infty} d λ λ^{4} ρ (λ) = α (3 + 2 α) .

(74)

V. MOMENTS AND CORRELATION FUNCTIONS OF CHARACTERISTIC POLYNOMIALS OF MPE AND GMPE

In this section, we describe a computation of the p-point correlation function of the characteristic polynomial defined in Eq. (6) for Marchenko–Pastur ensembles. Such objects are known to provide insights into spectral correlations and, as such, are interesting on their own and attract much attention. Note also that the mean total number of equilibria is expressed in Eq. (4) in terms of the object not dissimilar to $Z_{β, p} (λ; T)$ defined in Eq. (6). However, the presence of the absolute value forced us to resort only to asymptotic analysis for N → ∞ relying on the strong self-averaging hypothesis. In contrast, we will show that the moments that do not contain absolute values can be computed exactly for any finite N without any approximation.

Let us now obtain an explicit expression for $Z_{β, p} (λ; T)$ ⁠. We first consider the case of a fixed diagonal matrix T with elements that are not necessarily positive. Introduce p pairs of N-dimensional vectors ${{\tilde{ψ}}_{j}, ψ_{j}}$ with anticommuting/Grassmann components such that ψ_iψ_j = −ψ_jψ_i, ∫dψ = 0, and ∫ψdψ = 1. Then, each of the p determinant factors in the right-hand side of Eq. (6) can be written as the standard Berezin integral $\det (A) = \int D ψ D \tilde{ψ} \exp ({\tilde{ψ}}^{T} A ψ) = \int D ψ D \tilde{ψ} \exp (- Tr (A {\tilde{ψ}}^{T} ψ))$ ⁠, yielding

\begin{align} Z_{β, p} (λ; T) & = {⟨\prod_{j = 1}^{p} \det (λ_{j} I_{N} - W)⟩}_{K} = \int D ψ D \tilde{ψ} e^{- \sum_{j = 1}^{p} λ_{j} Tr ψ_{j} {\tilde{ψ}}_{j}^{T}} \prod_{a = 1}^{M} {⟨e^{\frac{T_{a}}{2} Tr Q^{(s)} k_{a} k_{a}^{T}}⟩}_{k_{a}} \\ = \int D ψ D \tilde{ψ} e^{- \sum_{j = 1}^{p} λ_{j} Tr ψ_{j} {\tilde{ψ}}_{j}^{T}} \prod_{a = 1}^{M} \det {(I_{N} - \frac{T_{a}}{N} Q^{(s)})}^{- β / 2}, \end{align}

(75)

where Q^(s) = Q + Q^T is the symmetric part of $Q = \sum_{j = 1}^{p} ψ_{j} {\tilde{ψ}}_{j}^{T}$ ⁠, and in the second line, we have used that ${⟨\dots⟩}_{k_{a}}$ is the expectation with respect to the Gaussian density $p (k_{a}) \propto e^{- N β k_{a}^{†} k_{a} / 2}$ for the ath column of K, which we evaluated in the second line of Eq. (75). Introducing the 2p × 2p matrix

\hat{Q} = [\begin{matrix} {\tilde{ψ}}_{i}^{T} ψ_{j} & {\tilde{ψ}}_{i}^{T} {\tilde{ψ}}_{j} \\ - ψ_{i}^{T} ψ_{j} & - ψ_{i}^{T} {\tilde{ψ}}_{j} \end{matrix}],

(76)

satisfying $Tr {(Q^{(s)})}^{k} = - Tr {\hat{Q}}^{k}$ ⁠, one can show that

\det {(I_{N} - \frac{T_{a}}{N} Q^{(s)})}^{- β / 2} = \det {(I_{(2 / β) p} - \frac{T_{a}}{N} \hat{Q})}^{β / 2} .

(77)

As the integrand in Eq. (75) only depends on $\hat{Q}$ ⁠, one can use the so-called “bosonization” trick proposed in Ref. 47 and replace the integral over anticommuting variables by the invariant integration over unitary matrices U belonging to the group C(4/β)E(p),

Z_{β, p} (λ; T) = \frac{\int_{C (4 / β) E (p)} (U^{†} d U) \det U^{- \frac{β N}{2}} e^{\frac{β}{2} Tr Λ U} \prod_{a = 1}^{M} Θ_{a} (U)}{\int_{C (4 / β) E (p)} (U^{†} d U) \det U^{- \frac{β N}{2}} e^{\frac{β}{2} Tr U}},

(78)

where $Θ_{a, β} (U) = \det {(I_{(2 / β) p} - \frac{T_{a}}{N} U)}^{β / 2}$ and

Λ = I_{2 / β} \otimes [\begin{matrix} λ_{1} & 0 & \dots & 0 \\ 0 & λ_{2} & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & \dots & λ_{p} \end{matrix}] .

(79)

Note that for β = 1, the integration over the group C4E(p) = CSE(p) defines the circular symplectic ensemble, i.e., the ensemble of unitary matrices U with skew symmetric sub-blocks,

U = [\begin{matrix} U_{11} & U_{12} \\ U_{21} & U_{11}^{T} \end{matrix}],

(80)

that satisfy

P^{- 1} U P = U^{T}, P = [\begin{matrix} O_{p} & - I_{p} \\ I_{p} & O_{p} \end{matrix}] .

(81)

On the other hand, for β = 2, the group C2E(p) = CUE(p) gives the standard Circular Unitary Ensemble (CUE).

A. Fixed diagonal matrix T

For a fixed diagonal matrix T, we focus our attention on the special case where all λ_j’s are identical, i.e., λ = λ1 = (λ, …, λ). In that case, both integrands in the numerator and the denominator of Eq. (78) are unitary invariant, and hence, their ratio can be expressed solely in terms of the integrals over eigenvalues of the unitary matrix U. This yields the following expression:

Z_{β, p} (λ; T) \equiv Z_{β, p} (λ 1_{p}; T) = \frac{\int_{{[0, 2 π]}^{p}} \prod_{j = 1}^{p} d θ_{j} \prod_{j = 1}^{p} e^{- i N θ} h (λ; θ_{j}) \prod_{k < j} {|e^{i θ_{j}} - e^{i θ_{k}}|}^{\frac{4}{β}}}{\int_{{[0, 2 π]}^{p}} \prod_{j = 1}^{p} d θ_{j} \prod_{j = 1}^{p} e^{- i N θ + e^{i θ_{j}}} \prod_{k < j} {|e^{i θ_{j}} - e^{i θ_{k}}|}^{\frac{4}{β}}},

(82)

where we introduced a function

h (λ; θ) = \exp (λ e^{i θ}) \prod_{a = 1}^{M} (1 - \frac{T_{a}}{N} e^{i θ}) .

(83)

For β = 2, using that $\det_{1 \leq i, j \leq p} (x_{i}^{j - 1}) = \prod_{i < j}^{p} (x_{i} - x_{j})$ together with the Cauchy–Binet–Andréief formula,⁴⁸ Eq. (82) can be re-written as the ratio of determinants,

Z_{β = 2, p} (λ; T) = {⟨\det {(λ I - K T K^{T})}^{p}⟩}_{K} = \frac{\det {[g_{N, k - j} (λ; T)]}_{j, k = 1}^{p}}{\det {[q_{N, k - j}]}_{j, k = 1}^{p}},

(84)

where

g_{N, m} (λ; T) = \int \frac{d θ}{2 π} e^{- i N θ} h (λ; θ) e^{i m θ} {= \oint}_{| z | = 1} \frac{d z}{2 i π z^{N + m + 1}} e^{λ z} \prod_{a = 1}^{M} (1 - \frac{T_{a}}{N} z),

(85)

q_{N, m} = g_{N, m} (1; O) = \int \frac{d θ}{2 π} e^{- i N θ} e^{e^{i θ}} e^{i m θ} {= \oint}_{| z | = 1} \frac{d z}{2 i π z^{N + m + 1}} e^{z} .

(86)

For each integral here, there is a unique pole of order N + m + 1 at z = 0, and using the residue theorem yields

g_{N, m} (λ; T) = λ^{N + m} \sum_{l = 0}^{N + m} \frac{{(- 1)}^{l}}{(N + m - l)!} e_{l} (\frac{T}{N λ}), q_{N, m} = \frac{1}{(N + m)!},

(87)

where e_l(X) is the lth elementary symmetric polynomial defined as

e_{l} (X_{1}, \dots, X_{n}) = \sum_{1 \leq i_{1} < i_{2} < \dots < i_{l} \leq n} \prod_{j = 1}^{l} X_{i_{j}} .

(88)

Similarly, for β = 1, using the identity

\prod_{j < k} | e^{i θ_{j}} - e^{i θ_{k}} |^{4} = \prod_{j = 1}^{p} e^{- 2 i (p - 1) θ_{j}} \det_{1 \leq k \leq 2 p, 1 \leq j \leq p} [e^{i (k - 1) θ_{j}}; (k - 1) e^{i (k - 2) θ_{j}}]

(89)

and applying further De Bruijn’s identity, the pth moment of the characteristic polynomial can be expressed as a ratio of Pfaffians,

Z_{β = 1, p} (λ; T) = {⟨\det {(λ I - K T K^{T})}^{p}⟩}_{K} = \frac{Pf {[(j - k) g_{N, 2 p + 1 - (k + j)} (λ; T)]}_{j, k = 1}^{2 p}}{Pf {[(j - k) q_{N, 2 p + 1 - (k + j)}]}_{j, k = 1}^{2 p}} .

(90)

In Fig. 3, we have compared our exact analytical results for $Z_{β, p} (λ; T)$ as a function of λ given in Eqs. (84) and (90), respectively, for β = 1, 2 taking fixed diagonal matrix T with numerical simulations, obtaining excellent agreement.

FIG. 3.

View large Download slide

Left: comparison between the analytical expression for $Z_{2, 2} (λ, T)$ in Eq. (84) (dark line) for N = 5, α = M/N = 2 and fixed T⁽¹⁾ = (0.4339, 1.0562, −0.3710, 0.3090, −0.8655, −1.0482, −1.6036, −0.7595, 1.1370, −1.0643)^T and numerical simulations (markers), plotted as a function of λ. Right: comparison between the analytical expression for $Z_{1, 2} (λ, T)$ in Eq. (90) (dark line) for N = 4, α = M/N = 2 and fixed T⁽²⁾ = (0.4339, 1.0562, −0.3710, 0.3090, −0.8655, −1.0482, −1.6036, −0.7595)^T and numerical simulations (markers), plotted as a function of λ. For each panel, the number of samples is 10⁵.

Let us finally mention that in the case λ = 0, one can use the following identities from Ref. 49:

\begin{align} \int_{{[0, 2 π]}^{p}} \prod_{j = 1}^{p} (\frac{d θ_{j} e^{- i N θ_{j}}}{2 π} e^{e^{i θ_{j}}}) \prod_{j < k} | e^{i θ_{k}} - e^{i θ_{j}} |^{\frac{4}{β}} \\ = M_{p} (- \frac{N}{2}, \frac{N}{2}, \frac{2}{β}) {}_{1}F_{1}^{(\frac{β}{2})} (- \frac{N}{2}, - \frac{N}{2} + 1 + \frac{2}{β} (p - 1); - 1), \end{align}

(91)

\begin{align} \int_{{[0, 2 π]}^{p}} \prod_{j = 1}^{p} (\frac{d θ_{j} e^{- i N θ_{j}}}{2 π} \prod_{a = 1}^{M} (1 - \frac{T_{a}}{N} e^{i θ_{j}})) \prod_{j < k} | e^{i θ_{k}} - e^{i θ_{j}} |^{\frac{4}{β}} \\ = M_{p} (- \frac{N}{2}, \frac{N}{2}, \frac{2}{β}) {}_{2}F_{1}^{(\frac{2}{β})} (- p, \frac{β N}{4}; (1 - p) + \frac{β}{2} (\frac{N}{2} - 1); - \frac{T}{N}), \end{align}

(92)

with

M_{p} (a, b, c) = \prod_{j = 0}^{p - 1} \frac{Γ (c j + a + b + 1) Γ (c (j + 1) + 1)}{Γ (c j + a + 1) Γ (c j + b + 1) Γ (1 + c)},

(93)

where the hypergeometric function of matrix argument is defined as

{}_{p}F_{q}^{(α)} (a; b; X) = \sum_{k = 0} \sum_{κ ⊢ k} \frac{\prod_{j = 1}^{p} {(a_{j})}_{κ}^{(α)}}{\prod_{j = 1}^{q} {(b_{j})}_{κ}^{(α)}} \frac{C_{κ}^{(α)} {(X)}^{k}}{k!},

(94)

with κ standing for a partition of k and ${(a)}_{κ}^{(α)}$ being the generalized Pochhammer symbol associated with the partition κ = (κ₁, …, κ_m),

{(a)}_{κ}^{(α)} = \prod_{i = 1}^{m} \prod_{j = 1}^{κ_{i}} (a - \frac{i - 1}{α} + j - 1),

(95)

and the function $C_{κ}^{(α)} (X)$ is the so-called C-normalization of the Jack polynomial.⁴⁹ This allows us to express the function $Z_{β, p} (0; T)$ via a ratio of the hypergeometric functions as

Z_{β, p} (0; T) = {⟨\det {(- W)}^{p}⟩}_{K} = \frac{{}_{2}F_{1}^{(\frac{2}{β})} (- p, \frac{β N}{4}; (1 - p) + \frac{β}{2} (\frac{N}{2} - 1); - \frac{T}{N})}{{}_{1}F_{1}^{(\frac{β}{2})} (- \frac{N}{2}, - \frac{N}{2} + 1 + \frac{2}{β} (p - 1); - 1)} .

(96)

B. Random diagonal matrix T

Let us now consider the case of random diagonal matrix T and average the p-point characteristic polynomial $Z_{β, p} (λ; T)$ over realizations of T, denoting these expectations as $Z_{β, p} (λ) = {⟨Z_{β, p} (λ; T)⟩}_{T}$ ⁠. In contrast to the case of fixed T, we now consider λ_j’s that are not necessarily identical but focus only on the low values p = 1 and p = 2 for β = 1, 2. Interestingly, we find that the ensuing formulas are universal in the sense that $Z_{β, p} (λ)$ depends on the whole distribution of T only via the lowest moments $⟨T^{k}⟩ = \int_{- \infty}^{\infty} d t t^{k} n (t), k = 1, 2$ ⁠, of the i.i.d. elements of T, assuming that those moments are finite.

1. The case p = 1

In such a case, we can directly use the results in Eqs. (84) and (90) for fixed T and take the average value with respect to the i.i.d. diagonal elements T_a of T. This yields

Z_{β, 1} (λ) = {⟨Z_{β, 1} (λ; T)⟩}_{T} = N {! \oint}_{| z | = 1} \frac{d z}{2 π i z^{N + 1}} e^{λ z} {(1 - \frac{⟨T⟩}{N} z)}^{M} .

(97)

In the particular case of a symmetric distribution n(t) = n(−t), we have $⟨T⟩ = 0$ and (97) simplifies drastically: $Z_{β, 1} (λ) = ⟨\det (λ I - W)⟩ = λ^{N}$ ⁠. If $⟨T⟩ \neq 0$ instead and assuming $| ⟨T⟩ | < N$ ⁠, $Z_{β, 1} (λ)$ is expressed for finite N as

Z_{β, 1} (λ) = {(\frac{⟨T⟩}{N})}^{N} \frac{M!}{(M - N)!} {}_{1}F_{1} (- N; 1 + M - N; \frac{λ N}{⟨T⟩}),

(98)

where we remind the definition of the standard hypergeometric function,

{}_{p}F_{q} (a; b; x) = \sum_{k = 0} \frac{\prod_{j = 1}^{p} {(a_{j})}_{k}}{\prod_{j = 1}^{q} {(b_{j})}_{k}} \frac{x^{k}}{k!},

(99)

with ${(a)}_{n} = \prod_{i = 0}^{n - 1} (a + i)$ being the Pochhammer symbol.

Note that the ratio $Z_{β, 1} (λ) / λ^{N}$ only depends on the rescaled variable $x = λ / ⟨T⟩$ ⁠. In the limit N → ∞ and for $⟨T⟩ \neq 0$ ⁠, changing the variable from z to $w = (⟨T⟩ / N) z$ in Eq. (97), one obtains

\frac{Z_{β, 1} (λ)}{λ^{N}} = \frac{N!}{N^{N} e^{- N}} \oint \frac{d w}{2 π i w} e^{N L_{1} (w; x, α)}, x = \frac{λ}{⟨T⟩},

(100)

L_{1} (w; x, α) = x w - \ln x w + α \ln (1 - w) - 1,

(101)

which can, then, be conveniently analyzed using Laplace’s method. The saddle points of $L_{1} (w; x, α)$ in the w domain are given by $w_{\pm} (x, α) = (1 - α + x \pm \sqrt{{(α - 1 - x)}^{2} - 4 x}) / (2 x)$ and are complex conjugate for (α − 1 − x)² < 4x, i.e., x ∈ [x₋(α), x₊(α)], where $x_{\pm} (α) = {(1 \pm \sqrt{α})}^{2}$ or both real in the converse case. Interestingly, the boundaries of the interval where w_± are complex conjugate correspond to the boundaries of the support of the classic Marchenko–Pastur distribution, i.e., $x_{\pm} (α) = {(1 \pm \sqrt{α})}^{2}$ ⁠. In the large N limit, a careful analysis of Eq. (100) (see Appendix A for details) allows us to obtain

Ξ_{1} (x; α) = \lim_{N \to \infty} \frac{1}{N} \ln |\frac{Z_{β, 1} (λ)}{λ^{N}}| = \{\begin{cases} L_{1} (w_{+}; x, α), & x < x_{-} (α), \\ R [L_{1} (w_{\pm}; x, α)], & x_{-} (α) \leq x \leq x_{+} (α), \\ L_{1} (w_{-}; x, α), & x > x_{+} (α) . \end{cases})

(102)

The function Ξ₁(x; α) changes the sign over the real interval and behaves in the limit x → ±∞ as

Ξ_{1} (x; α) \approx - \frac{α}{x}, x = \frac{λ}{⟨T⟩} \to \pm \infty,

(103)

which matches smoothly the expression $\ln | Z_{β, 1} (λ) / λ^{N} | = 0$ for $⟨T⟩ = 0$ ⁠. This result is in stark contrast with the behavior of Ξ_tot(μ; α), which is positive for any value of μ > 0 and vanishes as μ → ∞ as in Eq. (59), i.e., much faster than the algebraic decay, emphasizing the crucial role played by the absolute value in definition (35). In Fig. 4, we show a comparison between the analytical expression of Ξ₁(x, α) for α = 2 and the finite N = 20 expression of $N^{- 1} \ln | Z_{1, 1} (x) / x^{N} |$ with $⟨T⟩ = 1$ ⁠, showing excellent agreement.

FIG. 4.

View large Download slide

Comparison between the analytical expression of Ξ₁(x, α) for α = 2 (black solid line) obtained from Eq. (102) and the expression of $\frac{1}{N} \ln | Z_{1, 1} (x) / x^{N} |$ (blue circles) for $⟨T⟩ = 1$ ⁠, N = 20, and α = 2, plotted as a function of x. The number of samples is 10⁷.

2. Calculation for p = 2 and β = 1, 2

In order to obtain an exact finite N expression of the two-point correlation function $Z_{β, 2} (λ) = {⟨Z_{β, 2} (λ; T)⟩}_{T}$ for λ = (λ₁, λ₂), we start from Eq. (78) and introduce an explicit representation of matrices U ∈ C(4/β)E(2), separately for β = 1 and β = 2. In the special case of vanishing moment $⟨T⟩ = 0$ ⁠, one is able to perform the integration over matrices U (see Appendix B for details) and obtain the following expression:

\begin{align} Z_{β, 2} (λ = (λ_{1}, λ_{2})) = {(\frac{⟨T^{2}⟩}{N^{2}})}^{N} \frac{β}{2} Γ (N + \frac{2}{β} + 1) \frac{M!}{(M - N)!} {}_{1}F_{2} (- N; \frac{2}{β} + 1, 1 + M - N; - \frac{λ_{1} λ_{2} N^{2}}{⟨T^{2}⟩}) . \end{align}

(104)

Note that $Z_{β, 2} (λ)$ only depends on λ = (λ₁, λ₂) through the product λ₁λ₂, and the ratio $Z_{β, 2} (λ = (λ_{1}, λ_{2})) / {(λ_{1} λ_{2})}^{N}$ only depends on the rescaled variable $y = λ_{1} λ_{2} / ⟨T^{2}⟩$ ⁠. This is highly unusual as similar objects for Gaussian or Wishart matrices become in the large-N limit a function of the difference |λ₁ − λ₂|; see Refs. 34 and 50. In Fig. 5, we show a comparison between our analytical result in Eq. (104) and numerical simulations, showing excellent agreement.

FIG. 5.

View large Download slide

Comparison between the analytical expression of ${⟨Z_{1, 2}⟩}_{T}$ in Eq. (104) (solid line) and numerical simulations (markers), plotted as a function of the product of spectral parameters x = λ₁λ₂. The left panel corresponds to α = 1, while the right one corresponds to α = 2. For each panel, we have fixed $⟨T^{2}⟩ = 0.5, ⟨T⟩ = 0, N = 6$ and the number of samples is 10⁶.

Let us consider the large N limit of Eq. (104). To this purpose, we re-write the ratio of Eq. (104) by ${(λ_{1} λ_{2})}^{N}$ as the finite sum,

\frac{Z_{β, 2} (λ)}{{(λ_{1} λ_{2})}^{N}} = \sum_{k = 0}^{N} \frac{β}{2 {(N^{2} y)}^{N}} \frac{{(N - k + 1)}_{k} Γ (M + 1) Γ (N + \frac{2}{β} + 1)}{{(\frac{2}{β} + 1)}_{k} Γ (M - N + 1) {(M - N + 1)}_{k}} \frac{{(N^{2} y)}^{k}}{k!},

(105)

where $y = λ_{1} λ_{2} / ⟨T^{2}⟩$ ⁠. Replacing M = αN with α ≥ 1 and transforming the sum over k into a Riemann integral over η = k/N, one can, then, evaluate the large N behavior of this expression using the Laplace method such that

Ξ_{2} (u; α) = \lim_{N \to \infty} \frac{1}{N} \ln |\frac{Z_{β, 2} (λ)}{{(λ_{1} λ_{2})}^{N}}| = \max_{η \in [0, 1]} L_{2} (η; u, α), u = | y | = \frac{| λ_{1} λ_{2} |}{⟨T^{2}⟩},

(106)

where (see Appendix C for details)

\begin{align} L_{2} (η; u, α) & = α \ln (α) - (α + η - 1) \ln (α + η - 1) \\ + η \ln (u (1 - η)) + 2 η - 2 η \ln (η) - \ln (1 - η) - 2 - \ln u . \end{align}

(107)

The saddle-point equation yields a third order algebraic equation for the value η that maximizes the action $L_{2} (η; u, α)$ ⁠,

{(\partial_{η} L_{2}|}_{η = η_{*}} = - \ln (α + η_{*} - 1) + \ln (1 - η_{*}) + \ln (u (1 - η_{*})) - 2 \ln η_{*} = 0,

(108)

\Rightarrow η_{*} = \frac{1 - α}{3} + r_{*} w i t h - \frac{Δ}{27} + \frac{Θ}{3} r_{*} + r_{*}^{3} = 0,

(109)

w h e r e Δ = 9 u (2 + α) - 2 {(α - 1)}^{3}, Θ = 3 u - {(α - 1)}^{2}

(110)

and where the correct solution is obtained using Cardano’s formula as (see Appendix C)

η_{*} = \frac{1 - α}{3} - \frac{\sqrt[3]{2} Θ}{3 \sqrt[3]{Δ + \sqrt{Δ^{2} + 4 Θ^{3}}}} + \frac{\sqrt[3]{Δ + \sqrt{Δ^{2} + 4 Θ^{3}}}}{3 \sqrt[3]{2}} .

(111)

Thus, for any value of u and α, one can show that

Ξ_{2} (u; α) = \lim_{N \to \infty} \frac{1}{N} \ln |\frac{Z_{β, 2} (λ)}{{(λ_{1} λ_{2})}^{N}}| = L_{2} (η_{*}; u, α), u = \frac{| λ_{1} λ_{2} |}{⟨T^{2}⟩} .

(112)

Similarly as for the total complexity Ξ_tot(μ; α), the function Ξ₂(u; α) is positive for any finite u and vanishes only as u → ∞ although in a much slower fashion compared to Eq. (59), i.e.,

Ξ_{2} (u; α) \approx \frac{α}{u}, u \to \infty .

(113)

On the other hand, one can obtain that the function $Z_{β, 2} (λ)$ achieves a value exponentially large in N for λ₁λ₂ = 0, with

\lim_{N \to \infty} \frac{1}{N} \ln Z_{β, 2} (λ) = α \ln α + (α - 1) \ln (α - 1) - 2 .

(114)

In Fig. 6, we show a comparison between the expression of Ξ₂(u; α) in Eq. (112) and numerical simulations of $N^{- 1} \ln | Z_{1, 2} (\sqrt{u}, \sqrt{u}) / u^{N} |$ and $N^{- 1} \ln | Z_{1, 2} (1, u) / u^{N} |$ ⁠. In the large N limit, it is expected that $\lim_{N \to \infty} N^{- 1} \ln | Z_{1, 2} (λ_{1}, λ_{2}) / {(λ_{1} λ_{2})}^{N} |$ only depends on the rescaled parameter u = |λ₁λ₂| and these quantities should, thus, coincide. The agreement is reasonably good for low values of u, but there is some discrepancy (probably resulting from finite N corrections) for larger values of u.

FIG. 6.

View large Download slide

Comparison between the analytical expression of Ξ₂(x, α) for α = 2 (black solid line) obtained from Eq. (112) and the expression of $N^{- 1} \ln | Z_{1, 2} (\sqrt{u}, \sqrt{u}) / u^{N} |$ (blue circles) and $N^{- 1} \ln | Z_{1, 2} (1, u) / u^{N} |$ (orange triangles) for $⟨T⟩ = 0$ ⁠, $⟨T^{2}⟩ = 0$ ⁠, N = 20, and α = 2, plotted as a function of x. The number of samples is 10⁷. It is expected that all three curves coincide in the limit N → ∞. The agreement is rather good for low values of u, but there is some discrepancy for larger values of u, most probably due to finite N effects.

VI. CONCLUSION

In this article, we have introduced a model of the complex landscape resulting from a superposition of a large number M of random plane waves (larger than the dimension N of the Euclidean space). We have computed exactly the total annealed complexity of this random landscape confined by an isotropic harmonic potential of strength μ in the limit N, M → ∞ with α = M/N > 1. In contrast to previously studied cases, see, e.g., Ref. 12, this system does not display a “topological trivialization transition,” i.e., the total annealed complexity is strictly positive for any finite μ.

The random part of the Jacobian associated with such a random landscape defines a random matrix W = KTK^T belonging to a new random matrix ensemble that we denote “Gaussian Marchenko–Pastur Ensemble” (GMPE). We have studied the properties of this ensemble by computing both the average density of eigenvalues in the large N limit and some specific correlation functions of the characteristic polynomial. Rather surprisingly, we have shown that in the large N limit, the two-point correlation function depends on the product of the two spectral parameters instead of their distance, as is the case for standard random matrix ensembles.^34,50

The unforeseen behavior of the total annealed complexity naturally invites a further study of the Gibbs measure associated with this random landscape at fixed inverse temperature β = 1/T. Indeed, for translationally and rotationally invariant Gaussian models,¹⁴ it was shown that the trivial phase naturally corresponds to the replica symmetric phase in the Parisi framework. The absence of such a phase in the complexity might indicate that the replica symmetry breaking is operative for any μ at T = 0 and possibly also for positive values of the inverse temperature in this system. This will be the subject of a future publication.

Finally, it would be worth performing some numerical studies on the number of stationary points, their type, and other landscape statistics for the type of landscapes studied in the present paper. We hope that for moderately big values of N, such simulations should be more feasible than for landscapes based on more standard rotationally and translationally invariant Gaussian models.

DEDICATION

Submitting this article to the collection celebrating Freeman Dyson, one of the founding fathers of the modern random matrix theory, we would also like to dedicate it to Vladimir Marchenko who turns 100 years old in July 2022 and whose pioneering 1967 paper with Leonid Pastur deeply influenced the subsequent development of the field.

ACKNOWLEDGMENTS

We would like to thank Pierre Le Doussal for his interest in the work and, in particular, for pointing out a few relevant references. The research by Bertrand Lacroix-A-Chez-Toine and Yan V. Fyodorov was supported by the EPSRC under Grant No. EP/V002473/1 (Random Hessians and Jacobians: theory and applications). Sirio Belga Fedeli acknowledges the support from the EPSRC Centre for Doctoral Training in Cross-Disciplinary Approaches to Non-Equilibrium Systems through Grant No. EP/L015854/1.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Bertrand Lacroix-A-Chez-Toine: Formal analysis (equal); Writing – original draft (equal); Writing – review & editing (equal). Yan V. Fyodorov: Conceptualization (equal); Supervision (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Sirio Belga Fedeli: Formal analysis (equal); Writing – original draft (equal); Writing – review & editing (equal).

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

APPENDIX A: LAPLACE APPROXIMATION FOR $Z_{β, 1}$

Let us start from Eq. (100) that we reproduce here,

\frac{Z_{β, 1} (λ)}{λ^{N}} = \frac{N!}{N^{N} e^{- N}} \oint \frac{d w}{2 π i w} e^{N L_{1} (w; x, α)},

(A1)

L_{1} (w; x, α) = x w - \ln x w + α \ln (1 - w) - 1,

(A2)

where $x = λ / ⟨T⟩$ ⁠. The saddle points of $L_{1} (w; x, α)$ are given by

w_{\pm} = \frac{1 - α + x \pm \sqrt{{(α - 1 - x)}^{2} - 4 x}}{2 x} .

(A3)

The contribution of each of them depends on the sign of Δ₁ = (α − 1 − x)² − 4x.

For Δ₁ > 0, corresponding either to x < x₋(α) or x > x₊(α) with $x_{\pm} (α) = {(1 \pm \sqrt{α})}^{2}$ ⁠, the saddle points w_± both lie on the real line. In order to obtain the stability of the fixed points, let us compute

{(\partial_{w}^{2} L_{1} (w; x, α)|}_{w = w_{\pm}} = \frac{1}{w_{\pm}^{2}} - \frac{α}{{(1 - w_{\pm})}^{2}},

(A4)

which is real and positive only for w_± real and inside the interval between its two roots $({(1 - \sqrt{α})}^{- 1}, {(1 + \sqrt{α})}^{- 1})$ ⁠. These two roots translate in terms of x as the two edges x_±(α) of the usual Marchenko–Pastur distribution, and one can check that for x > x₊(α), one has $w_{-} \in ({(1 - \sqrt{α})}^{- 1}, {(1 + \sqrt{α})}^{- 1})$ (and w₊ outside of this interval), and conversely for x < x₋(α), one has $w_{+} \in ({(1 - \sqrt{α})}^{- 1}, {(1 + \sqrt{α})}^{- 1})$ (and w₋ outside).

Conversely, for Δ₁ < 0, i.e., x ∈ (x₋(α), x₊(α)), w_± are complex conjugate complex numbers and both contribute in encircling the origin. In this region, the real part of the term in Eq. (A4) is positive and the same for both w_± such that both saddle-points have equal contribution. The imaginary part of the term in Eq. (A4) changes its sign for x = (α − 1)²/(α + 1).

Summarizing, one obtains for N ≫ 1,

\frac{Z_{β, 1} (λ)}{λ^{N}} = \{\begin{cases} I (w_{+}), & x < x_{-} (α), \\ - I (w_{+}) - I (w_{-}), & x_{-} (α) < x < \frac{{(α - 1)}^{2}}{α + 1}, \\ I (w_{+}) + I (w_{-}), & \frac{{(α - 1)}^{2}}{α + 1} < x < x_{+} (α), \\ I (w_{-}), & x > x_{+} (α), \end{cases})

(A5)

w h e r e I (w_{\pm}) = \frac{e^{N L_{1} (w_{\pm}; x, α)}}{\sqrt{2 π N w_{\pm}^{2} \partial_{w}^{2} L_{1} (w_{\pm}; x, α)}} .

(A6)

Finally, for x = x_±(α), the saddle points are unique and of order 2:

For x = x₊(α), one has that $w_{+} = {(1 + \sqrt{α})}^{- 1}$ such that
${(\partial_{w}^{2} L_{1} (w; x, α)|}_{w = w_{+}} = 0, {(\partial_{w}^{3} L_{1} (w; x, α)|}_{w = w_{+}} = - 2 \frac{{(1 + \sqrt{α})}^{4}}{\sqrt{α}} .$
(A7)
Therefore, the steepest descent directions at the point w₊ are {2π/3, 4π/3, 0}. One can locally deform the contour, introducing the scalar t ∈ (−∞, +∞) and approaching and leaving w₊ with w(t) = e^i4π/3t + w₊ and w(t) = e^i2π/3t + w₊, respectively, yielding
$\begin{align} \frac{Z_{β, 1} (λ)}{λ^{N}} & = \frac{e^{N L_{1} (w_{+}; x, α)}}{2 i π w_{+}} (e^{i \frac{2 π}{3}} \int_{0}^{+ \infty} d t e^{- N \frac{2 {(1 + \sqrt{α})}^{4}}{3! \sqrt{α}} t^{3}} + e^{i \frac{4 π}{3}} \int_{+ \infty}^{0} d t e^{- N \frac{2 {(1 + \sqrt{α})}^{4}}{3! \sqrt{α}} t^{3}}) \\ = \frac{Γ (1 / 3)}{2 π 3^{1 / 6} {(N \frac{(1 + \sqrt{α})}{\sqrt{α}})}^{1 / 3}} e^{N L_{1} (w_{+}; x, α)} . \end{align}$
(A8)
For x = x₋(α), $w_{-} = {(1 - \sqrt{α})}^{- 1}$ such that
${(\partial_{w}^{2} L_{1} (w; x, α)|}_{w = w_{-}} = 0, {(\partial_{w}^{3} L_{1} (w; x, α)|}_{w = w_{-}} = 2 \frac{{(1 - \sqrt{α})}^{4}}{\sqrt{α}} .$
(A9)
The steepest descent directions at w₋ are {π/3, π, 5π/3}. The contour is deformed by parametrizing the path and approaching and leaving w₋ with w(t) = e^i1/3πt + w₋ and w(t) = e^i5/3πt + w₋, respectively, yielding
$\begin{align} \frac{Z_{β, 1} (λ)}{λ^{N}} & = \frac{e^{N L_{1} (w_{-}; x, α)}}{2 i π w_{-}} (e^{i \frac{π}{3}} \int_{0}^{+ \infty} d t e^{- N \frac{2 {(1 - \sqrt{α})}^{4}}{3! \sqrt{α}} t^{3}} + e^{i \frac{5 π}{3}} \int_{+ \infty}^{0} d t e^{- N \frac{2 {(1 - \sqrt{α})}^{4}}{3! \sqrt{α}} t^{3}}) \\ = \frac{Γ (1 / 3) {(- 1)}^{N}}{2 π 3^{1 / 6} {(N \frac{(\sqrt{α} - 1)}{\sqrt{α}})}^{1 / 3}} e^{N L_{1} (w_{-}; x, α)} . \end{align}$
(A10)
In each case, one obtains that the term
$Ξ_{1} (x = x_{\pm} (α); α) = \lim_{N \to + \infty} \frac{1}{N} \ln |\frac{Z_{β, 1} (λ)}{λ^{N}}| = L_{1} (w_{\pm}; x_{\pm} (α), α),$
(A11)
where $x = λ / ⟨T⟩$ ⁠.

APPENDIX B: FINITE N EXPRESSION OF $Z_{β, 2}$

Let us first consider β = 2. In this case, Λ = diag(λ₁, λ₂), and we introduce four variables ψ_ij with i ≤ j and i, j = 1, 2, and ϕ such that 0 ≤ ψ_ij ≤ 2π and 0 ≤ ϕ ≤ π/2. Therefore, a possible parametrization of U is given as follows:⁵¹

U = [\begin{matrix} e^{i (ψ_{11} + ψ_{12})} \cos ϕ & e^{i ψ_{12}} \sin ϕ \\ - e^{i (ψ_{11} + ψ_{22})} \sin ϕ & e^{i ψ_{22}} \cos ϕ \end{matrix}] .

(B1)

The Euclidean line element (ds)² = ∑_i,j|dU_ij|² reads in this parametrization

\begin{align} {(d s)}^{2} & = {(d ψ_{11})}^{2} + {(d ψ_{12})}^{2} + {(d ψ_{22})}^{2} + 2 {(d ϕ)}^{2} \\ + 2 {(\cos ϕ)}^{2} (d ψ_{11}) (d ψ_{12}) + 2 {(\sin ϕ)}^{2} (d ψ_{11}) (d ψ_{22}) . \end{align}

(B2)

The associated Jacobian is given by $\sqrt{2 \det M}$ ⁠, where

M = [\begin{matrix} 1 & {(\cos ϕ)}^{2} & {(\sin ϕ)}^{2} \\ {(\cos ϕ)}^{2} & 1 & 0 \\ {(\sin ϕ)}^{2} & 0 & 1 \end{matrix}] .

(B3)

The normalized Haar measure is obtained calculating $\int \sqrt{\det M} d ψ_{11} d ψ_{12} d ψ_{22} d ϕ$ and is given by $d U = \frac{1}{4 π^{3}} \cos ϕ \sin ϕ d ψ_{11} d ψ_{12} d ψ_{22} d ϕ$ ⁠. With the further assumption that $⟨T⟩ = 0$ ⁠, from what stated above follows

\det U = e^{i (ψ_{11} + ψ_{12} + ψ_{22})},

(B4)

Tr Λ U = \cos ϕ (e^{i (ψ_{11} + ψ_{12})} λ_{1} + e^{i (ψ_{22})} λ_{2}),

(B5)

⟨\det (I_{2} - \frac{T_{a}}{N} U)⟩ = 1 + \frac{⟨T^{2}⟩}{N^{2}} e^{i (ψ_{11} + ψ_{12} + ψ_{22})} .

(B6)

The denominator of Eq. (78) reads

\begin{align} \frac{1}{4 π^{3}} \int_{0}^{2 π} d ψ_{11} \int_{0}^{2 π} d ψ_{12} \int_{0}^{2 π} d ψ_{22} d ϕ \cos ϕ \sin ϕ e^{\cos ϕ (e^{i (ψ_{11} + ψ_{12})} + e^{i ψ_{22}})} e^{- i N (ψ_{11} + ψ_{12} + ψ_{22})} \\ = \frac{2}{{(N!)}^{2}} \int_{0}^{\frac{π}{2}} d ϕ {(\cos ϕ)}^{2 N + 1} \sin ϕ = \frac{1}{{(N!)}^{2}} \frac{1}{N + 1}, \end{align}

(B7)

while the numerator reads

\begin{align} \frac{1}{4 π^{3}} \int_{0}^{2 π} d ψ_{11} \int_{0}^{2 π} d ψ_{12} \int_{0}^{2 π} d ψ_{22} d ϕ \cos ϕ \sin ϕ \\ \times e^{\cos ϕ (e^{i (ψ_{11} + ψ_{12})} λ_{1} + e^{i (ψ_{22})} λ_{2})} e^{- i N (ψ_{11} + ψ_{12} + ψ_{22})} {(1 + \frac{⟨T^{2}⟩}{N^{2}} e^{i (ψ_{11} + ψ_{12} + ψ_{22})})}^{M} \\ = \sum_{k = 0}^{N} \frac{2}{k! (N - k!)!} \frac{{(λ_{1} λ_{2})}^{k}}{(N - k)!} \frac{α N!}{(α N - k)!} {(\frac{⟨T^{2}⟩}{N^{2}})}^{k} \frac{1}{2 - 2 k + 2 N} . \end{align}

(B8)

Therefore, for β = 2, Eq. (104) is given by the ratio of Eqs. (B7) and (B8).

For β = 1, we impose $⟨T^{2}⟩ = 1$ to simplify the computations. The general case is retrieved by multiplying $Z_{2, β = 1} (λ)$ by ${⟨T^{2}⟩}^{N}$ and rescaling $λ_{1, 2} \to λ_{1, 2} / \sqrt{⟨T^{2}⟩}$ ⁠. We rename Λ = diag(λ₁, λ₂, λ₁, λ₂), while now the matrix U is parameterized with a 2 × 2 block matrix. Finally, we introduce six variables: h₁, h₂ ∈ (0, 1) and ϕ₁₁, ϕ₂₁, ϕ₂₂, ϕ₁₄ ∈ (0, 2π). Therefore (see Ref. 51),

U = [\begin{matrix} U_{11} & U_{12} \\ - U_{12} e^{i (ϕ_{12} + ϕ_{21} - 2 ϕ_{14})} & U_{11}^{t} \end{matrix}]

(B9)

with

U_{11} = [\begin{matrix} e^{i ϕ_{11}} h_{1} & e^{i ϕ_{22}} h_{2} \sqrt{1 - h_{1}^{2}} \\ e^{i ϕ_{21}} h_{2} \sqrt{1 - h_{1}^{2}} & - e^{- i (ϕ_{11} - ϕ_{21} - ϕ_{22})} h_{1} \end{matrix}]

(B10)

and

U_{12} = e^{i ϕ_{14}} \sqrt{1 - h_{1}^{2}} \sqrt{1 - h_{2}^{2}} [\begin{matrix} 0 & - 1 \\ 1 & 0 \end{matrix}] .

(B11)

\begin{align} {(d s)}^{2} & = 2 ({(d h_{1})}^{2} + h_{1}^{2} {(d ϕ_{11})}^{2} + 2 (1 - h_{1}^{2}) {(d h_{2})}^{2} - 4 d h (d h_{1}) (d h_{2})) \\ + \frac{2 {(d h_{1})}^{2}}{(1 - h_{1}^{2})} {(d h_{1})}^{2} + h_{2}^{2} (1 - h_{1}^{2}) {(d ϕ_{22})}^{2} + h_{2}^{2} (1 - h_{1}^{2}) {(d ϕ_{21})}^{2} + {(d h_{1})}^{2} \\ + h_{1}^{2} {(d ϕ_{22} + d ϕ_{21} - d ϕ_{11})}^{2} + \frac{2 {(h_{2} (1 - h_{1}^{2}) (d h_{2}) + h_{1} (1 - h_{2}^{2}) (d h_{1}))}^{2}}{(1 - h_{2}^{2}) (1 - h_{1}^{2})} \\ (+ (1 - h_{2}^{2}) (1 - h_{1}^{2}) {(d ϕ_{14})}^{2} + (1 - h_{2}^{2}) (1 - h_{1}^{2}) {(d ϕ_{22} + d ϕ_{21} - d ϕ_{14})}^{2}) . \end{align}

(B12)

To simplify the calculation, we can re-arrange the differential of the six variables such that the Jacobian is given by the product of the square roots of the matrices A₁ and A₂. The latter are

A_{1} = 4 [\begin{matrix} 1 + \frac{h_{1}^{2} (1 - h_{2}^{2}) + h_{2}^{2} h^{2}}{(1 - h_{1}^{2})} & 0 \\ 0 & (1 - h_{1}^{2}) + \frac{h_{2}^{2} (1 - h_{1}^{2})}{(1 - h_{2}^{2})} \end{matrix}]

(B13)

and

A_{2} = [\begin{matrix} 4 h^{2} & - 2 h^{2} & - 2 h^{2} & 0 \\ - 2 h^{2} & 2 & 2 - 2 h_{2}^{2} + 2 d^{2} h^{2} & - 2 (1 - h_{2}^{2}) (1 - h_{1}^{2}) \\ - 2 h^{2} & 2 - 2 h_{2}^{2} + 2 h_{2}^{2} h^{2} & 2 & - 2 (1 - h_{2}^{2}) (1 - h_{1}^{2}) \\ 0 & - 2 (1 - h_{2}^{2}) (1 - h_{1}^{2}) & - 2 (1 - h_{2}^{2}) (1 - h_{1}^{2}) & 4 (1 - h_{2}^{2} - h_{1}^{2} + h_{2}^{2} h^{2}) \end{matrix}] .

(B14)

We obtain the Haar measure by calculating

\int \sqrt{\det (A_{1}) \det (A_{2})} d h_{1} d h_{2} d ϕ_{11} d ϕ_{21} d ϕ_{22} d ϕ_{14},

(B15)

\det (A_{1}) = \frac{16}{1 - h_{2}^{2}},

(B16)

\det (A_{2}) = 64 h_{2}^{2} (1 - h_{2}^{2}) h_{1}^{2} {(1 - h_{1}^{2})}^{2} .

(B17)

Integrating out the six linear and angular variables, one obtains

d U = \frac{1}{2 π^{4}} h_{1} h_{2} (1 - h_{1}^{2}) (d h_{1} d h_{2} d ϕ_{11} d ϕ_{21} d ϕ_{22} d ϕ_{14}) .

(B18)

The denominator appearing in Eq. (78) is

\frac{1}{2 π^{3}} \int_{0}^{1} d h_{1} \int_{0}^{2 π} d ϕ_{11} \int_{0}^{2 π} d ϕ_{21} \int_{0}^{2 π} d ϕ_{22} e^{h e^{- i ϕ_{11}} (e^{2 i ϕ_{11}} - e^{i (ϕ_{21} + ϕ_{22})})} e^{- i N (ϕ_{21} + ϕ_{22})}

(B19)

= \frac{4}{{(N!)}^{2}} \int_{0}^{1} d h_{1} h_{1}^{2 n + 1} (1 - h_{1}^{2}) = \frac{2}{{(N!)}^{2}} \frac{1}{(2 + 3 N + N^{2})},

(B20)

whereas the numerator in Eq. (78) reads

\begin{align} \frac{1}{2 π^{3}} \int_{0}^{1} d h_{1} \int_{0}^{2 π} d ϕ_{11} \int_{0}^{2 π} d ϕ_{21} \int_{0}^{2 π} d ϕ_{22} h_{1} (1 - h_{1}^{2}) \\ \times e^{- i N (ϕ_{22} + ϕ_{21})} e^{h_{1} (λ_{1} e^{i ϕ_{11}} - λ_{2} e^{i (ϕ_{22} + ϕ_{21} - ϕ_{11})})} {(1 - \frac{e^{i (ϕ_{22} + ϕ_{21})}}{N^{2}})}^{M} . \end{align}

(B21)

The angular variables ϕ₂₂ and ϕ₂₁ in the integral above can be integrated out with the auxiliary variable ϕ = ϕ₂₂ + ϕ₂₁. The integral appearing in Eq. (B21) is simplified to

\sum_{k = 0}^{N} \frac{2 (α N!) N^{2 (α N - k)}}{π N^{2 α N} k! (N - k)! (α N - k)!} \int_{0}^{1} d h_{1} h_{1} (1 - h_{1}^{2}) {(h_{1} λ_{2})}^{N - k} \int_{0}^{2 π} d ϕ_{11} e^{h_{1} λ_{1} e^{i ϕ_{11}} - i (N - k) ϕ_{11}} .

(B22)

Working out the remaining integration, the expression above becomes

\sum_{k = 0}^{N} \frac{2 (α N!) {(λ_{1} λ_{2})}^{N}}{(k! {(λ_{1} λ_{2} N^{2})}^{k}) {((N - k)!)}^{2} (α N - k) (2 + 3 (N - k) + {(N - k)}^{2})!},

(B23)

and Eq. (104) is retrieved with the introduction of the hypergeometric function.

APPENDIX C: LAPLACE APPROXIMATION FOR $Z_{β, 2}$

We consider the case of centered random variables, i.e., ${⟨T⟩}_{i} = 0$ ⁠. In order to obtain the asymptotics of $Z_{β, 2}$ ⁠, it is convenient to start from Eq. (104) and to introduce $x = λ_{1} λ_{2} / ⟨T^{2}⟩$ ⁠. Expanding the terms within the hypergeometric function, we note

\begin{align} \frac{Z_{β, 2} (λ)}{{(λ_{1} λ_{2})}^{N}} & = \sum_{k = 0}^{N} \frac{β}{2} \frac{{(N - k + 1)}_{k} Γ (α N + 1) Γ (N + \frac{2}{β} + 1)}{{(\frac{2}{β} + 1)}_{k} Γ ((α - 1) N + 1) {((α - 1) N + 1)}_{k}} \frac{{(N^{2} x)}^{k - N}}{k!} \\ = \sum_{k = 0}^{N} a_{N} (x; k), \end{align}

(C1)

where ${(b)}_{k} = \frac{Γ (b + k)}{Γ (b)}$ is the Pochhammer symbol.

Let us first consider the case where x > 0. In that case, we can safely replace in the large N limit the summation above with an integral by rescaling k = Nη with η ∈ [0, 1]. Therefore, for N sufficiently large, one observes

a_{N} (x; N η) \approx \frac{β \sqrt{\frac{2 α}{π}} Γ (\frac{2}{β} + 1) η^{- \frac{2}{β} - 1}}{4 \sqrt{N} \sqrt{1 - η} \sqrt{(α + η - 1)}} e^{N L_{2} (η; x, α)},

(C2)

\begin{align} L_{2} (η; x, α) & = α \ln (α) - (α + η - 1) \ln (α + η - 1) \\ + η \ln (x (1 - η)) + 2 η - 2 η \ln (η) - \ln (1 - η) - 2 - \ln u . \end{align}

(C3)

The maximum of $L_{2}$ is located in η_* ∈ [0, 1]. The latter is the solution of the following equation:

{(\partial_{η} L_{2}|}_{η = η_{*}} = 0, \Rightarrow \ln (x - η_{*} x) - \ln (α + η_{*} - 1) - 2 \ln (η_{*}) = 0 .

(C4)

Hence, we need to retrieve the roots of p(η) = η²(α − 1 + η) − x(1 − η). First, one observes that lim_η→±∞p(η) = ±∞ and p(0) = −x. For η > 0, p(η) is strictly increasing since η² + 2η(α + η − 1) + x > 0. For x → 0⁺, the zero of p tends to 0. For x → +∞, it tends to 1. With the help of Cardano’s formula, the root is

\begin{align} η_{*} = \frac{1 - α}{3} - \frac{\sqrt[3]{2} Θ}{3 \sqrt[3]{Δ + \sqrt{Δ^{2} + 4 Θ^{3}}}} + \frac{\sqrt[3]{Δ + \sqrt{Δ^{2} + 4 Θ^{3}}}}{3 \sqrt[3]{2}}, \\ Δ = 9 x (2 + α) - 2 {(α - 1)}^{3}, Θ = 3 x - {(α - 1)}^{2} . \end{align}

(C5)

One can check that this solution η_* ∈ (0, 1), and in this interval, it corresponds to a maximum as

{(\frac{d^{2} L}{d η^{2}} (η)|}_{η = η_{*}} = - \frac{2}{η_{*}} + \frac{α}{(η_{*} - 1) (η_{*} + α - 1)} < 0 .

(C6)

Finally, one obtains that in the regime N ≫ 1,

\frac{Z_{β, 2} (λ)}{{(λ_{1} λ_{2})}^{N}} \approx \frac{β \sqrt{\frac{2}{π}} Γ (\frac{2}{β} + 1) {η^{*}}^{- \frac{2}{β} - 1} \sqrt{α N}}{4 \sqrt{1 - η_{*}} \sqrt{(α + η_{*} - 1)}} e^{N L (η_{*})} \int_{- \infty}^{+ \infty} d η e^{\frac{N}{2} \frac{d^{2} L}{d η^{2}} (η_{*}) {(η - η_{*})}^{2}} .

(C7)

Let us now consider the more involved case where x < 0. In that case, one needs to consider separately the terms for even and odd values of k by re-writing the expression as

\begin{align} \frac{Z_{β, 2} (λ)}{{(λ_{1} λ_{2})}^{N}} & = \sum_{k = 0}^{N / 2} \frac{{(N - 2 k + 1)}_{2 k + 1}}{{(1 + \frac{2}{β})}_{2 k} {((α - 1) N + 1)}_{2 k}} \frac{{(N^{2} | x |)}^{2 k - N}}{(2 k)!} \\ - \sum_{k = 0}^{N / 2 - 1} \frac{{(N - 2 k)}_{2 k + 1}}{{(1 + \frac{2}{β})}_{2 k + 1} {((α - 1) N + 1)}_{2 k + 1}} \frac{{(N^{2} | x |)}^{2 k + 1 - N}}{(2 k + 1)!} . \end{align}

(C8)

We may now in each term replace the sum over k with an integral over, respectively, η_e = 2k/N and η_o = (2k + 1)/N with η_o, η_e ∈ [0, 1]. The saddle point is the same as for x > 0, but one now needs to subtract the contribution from even and odd values of k. This is equivalent to doing a Taylor expansion in the parameter η_* of the expression in Eq. (C7). The result at exponential order in N is the same, but the finite N corrections are, thus, different.

REFERENCES

1.

T.

Castellani

and

A.

Cavagna

, “

Spin-glass theory for pedestrians

,”

J. Stat. Mech.

2005

(

05

),

P05012

.

https://doi.org/10.1088/1742-5468/2005/05/p05012

Crossref

2.

A.

Masoumi

,

A.

Vilenkin

, and

M.

Yamada

, “

Inflation in multi-field random Gaussian landscapes

,”

J. Cosmol. Astropart. Phys.

2017

(

12

),

035

.

https://doi.org/10.1088/1475-7516/2017/12/035

Crossref

3.

L. L.

Feng

,

S.

Hotchkiss

, and

R.

Easther

, “

The distribution of vacua in random landscape potentials

,”

J. Cosmol. Astropart. Phys.

2021

(

01

),

029

.

https://doi.org/10.1088/1475-7516/2021/01/029

Crossref

4.

A.

Choromanska

,

Y.

LeCun

, and

G.

Ben Arous

, “

Open problem: The landscape of the loss surfaces of multilayer networks

,” in

Proceedings of the 28th Conference on Learning Theory, PMLR

(

PMLR

,

2015

), Vol. 40, pp.

1756

–

1760

.

Google Scholar

5.

N. P.

Baskerville

,

J. P.

Keating

,

F.

Mezzadri

, and

J.

Najnudel

, “

The loss surfaces of neural networks with general activation functions

,”

J. Stat. Mech.: Theory Exp.

2021

,

064001

.

https://doi.org/10.1088/1742-5468/abfa1e

Crossref

6.

G. B.

Arous

,

S.

Mei

,

A.

Montanari

, and

M.

Nica

, “

The landscape of the spiked tensor model

,”

Commun. Pure Appl. Math.

72

(

11

),

2282

–

2330

(

2019

).

https://doi.org/10.1002/cpa.21861

Google Scholar

Crossref

7.

V.

Ros

,

G. B.

Arous

,

G.

Biroli

, and

C.

Cammarota

, “

Complex energy landscapes in spiked-tensor and simple glassy models: Ruggedness, arrangements of local minima, and phase transitions

,”

Phys. Rev. X

9

(

1

),

011003

(

2020

).

https://doi.org/10.1103/physrevx.9.011003

Google Scholar

Crossref

8.

A.

Maillard

,

G.

Ben Arous

, and

G.

Biroli

, “

Landscape complexity for the empirical risk of generalized linear models

,” in

Proceedings of Machine Learning Research

(

PMLR

,

2020

), Vol. 107, pp.

287

–

327

.

Google Scholar

9.

Y. V.

Fyodorov

and

R.

Tublin

, “

Optimization landscape in the simplest constrained random least-square problem

,”

J. Phys. A: Math. Theor.

55

(

24

),

244008

(

2022

).

https://doi.org/10.1088/1751-8121/ac6d8e

Google Scholar

Crossref

10.

A.

Cavagna

,

I.

Giardina

, and

G.

Parisi

, “

Stationary points of the Thouless-Anderson-Palmer free energy

,”

Phys. Rev. B

57

(

18

),

11251

–

11257

(

1998

).

https://doi.org/10.1103/physrevb.57.11251

Google Scholar

Crossref

11.

A.

Cavagna

,

J. P.

Garrahan

, and

I.

Giardina

, “

Quenched complexity of the mean-field p-spin spherical model with external magnetic field

,”

J. Phys. A: Math. Gen.

32

,

711

–

723

(

1999

).

https://doi.org/10.1088/0305-4470/32/5/004

Google Scholar

Crossref

12.

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

);

https://doi.org/10.1103/physrevlett.92.240601

Crossref

[PubMed]

Google Scholar

Erratum

93

,

149901(E)

(

2004

).

https://doi.org/10.1103/PhysRevLett.93.149901

13.

A. J.

Bray

and

D. S.

Dean

, “

Statistics of critical points of Gaussian fields on large-dimensional spaces

,”

Phys. Rev. Lett.

98

,

150201

(

2007

).

https://doi.org/10.1103/physrevlett.98.150201

Google Scholar

Crossref

PubMed

14.

Y. V.

Fyodorov

and

I.

Williams

, “

Replica symmetry breaking condition exposed by random matrix calculation of landscape complexity

,”

J. Stat. Phys.

129

,

1081

–

1116

(

2007

).

https://doi.org/10.1007/s10955-007-9386-x

Google Scholar

Crossref

15.

A.

Auffinger

,

G. B.

Arous

, and

J.

Černý

, “

Random matrices and complexity of spin glasses

,”

Commun. Pure Appl. Math.

66

(

2

),

165

–

201

(

2013

).

https://doi.org/10.1002/cpa.21422

Google Scholar

Crossref

16.

A.

Auffinger

and

G.

Ben Arous

, “

Complexity of random smooth functions on the high-dimensional sphere

,”

Ann. Probab.

41

(

6

),

4214

–

4247

(

2013

).

https://doi.org/10.1214/13-aop862

Google Scholar

Crossref

17.

Y. V.

Fyodorov

and

C.

Nadal

, “

Critical behavior of the number of minima of a random landscape at the glass transition point and the Tracy-Widom distribution

,”

Phys. Rev. Lett.

109

,

167203

(

2012

).

https://doi.org/10.1103/physrevlett.109.167203

Google Scholar

Crossref

PubMed

18.

E.

Subag

, “

The complexity of spherical p-spin model—A second moment approach

,”

Ann. Probab.

45

(

5

),

3385

–

3450

(

2017

).

https://doi.org/10.1214/16-aop1139

Google Scholar

Crossref

19.

Y. V.

Fyodorov

,

P.

Le Doussal

,

A.

Rosso

, and

C.

Texier

, “

Exponential number of equilibria and depinning threshold for a directed polymer in a random potential

,”

Ann. Phys.

397

,

1

–

64

(

2018

).

https://doi.org/10.1016/j.aop.2018.07.029

Google Scholar

Crossref

20.

V.

Ros

, “

Distribution of rare saddles in the p-spin energy landscape

,”

J. Phys. A: Math. Theor.

53

(

12

),

125002

(

2020

).

https://doi.org/10.1088/1751-8121/ab73ac

Google Scholar

Crossref

21.

Y. V.

Fyodorov

and

R.

Tublin

, “

Counting stationary points of the loss function in the simplest constrained least-square optimization

,”

Acta Phys. Pol., B

51

,

1663

–

1672

(

2020

).

https://doi.org/10.5506/aphyspolb.51.1663

Google Scholar

Crossref

22.

E.

Subag

and

O.

Zeitouni

, “

Concentration of the complexity of spherical pure p-spin models at arbitrary energies

,”

J. Math. Phys.

62

(

12

),

123301

(

2021

).

https://doi.org/10.1063/5.0070582

Google Scholar

Crossref

23.

G.

Ben Arous

,

P.

Bourgade

, and

B.

McKenna

, “

Landscape complexity beyond invariance and the elastic manifold

,” arXiv:2105.05051.

24.

J.

Grela

and

B. A.

Khoruzhenko

, “

Glass-like transition described by toppling of stability hierarchy

,”

J. Phys. A: Math. Theor.

55

(

15

),

154001

(

2022

).

https://doi.org/10.1088/1751-8121/ac56aa

Google Scholar

Crossref

25.

Y. V.

Fyodorov

, “

High-dimensional random fields and random matrix theory

,”

Markov Processes Relat. Fields

21

(

3

),

483

–

518

(

2015

), http://math-mprf.org/journal/articles/id1385/.

Google Scholar

26.

Y. V.

Fyodorov

and

P.

Le Doussal

, “

Topology trivialization and large deviations for the minimum in the simplest random optimization

,”

J. Stat. Phys.

154

(

1

),

466

–

490

(

2014

).

https://doi.org/10.1007/s10955-013-0838-1

Google Scholar

Crossref

27.

D.

Belius

,

J.

Černý

,

S.

Nakajima

, and

M.

Schmidt

, “

Triviality of the geometry of mixed p-spin spherical Hamiltonians with external field

,”

J. Stat. Phys.

186

(

1

),

12

(

2022

).

https://doi.org/10.1007/s10955-021-02855-6

Google Scholar

Crossref

28.

M. V.

Berry

, “

Regular and irregular semiclassical wavefunctions

,”

J. Phys. A: Math. Gen.

10

(

10

),

2083

–

2091

(

1977

).

https://doi.org/10.1088/0305-4470/10/12/016

Google Scholar

Crossref

29.

D.

Beliaev

,

V.

Cammarota

, and

I.

Wigman

, “

Two point function for critical points of a random plane wave

,”

Int. Math. Res. Not.

2019

(

9

),

2661

–

2689

.

https://doi.org/10.1093/imrn/rnx197

Crossref

30.

G.

Blatter

,

M. V.

Feigel’man

,

V. B.

Geshkenbein

,

A. I.

Larkin

, and

V. M.

Vinokur

, “

Vortices in high-temperature superconductors

,”

Rev. Mod. Phys.

66

,

1125

–

1388

(

1994

).

https://doi.org/10.1103/revmodphys.66.1125

Google Scholar

Crossref

31.

P.

Le Doussal

, “

Novel phases of vortices in superconductors

,”

Int. J. Mod. Phys. B

24

,

3855

–

3914

(

2010

).

https://doi.org/10.1142/s0217979210056384

Google Scholar

Crossref

32.

Note that in the Ph.D. thesis of one of the authors: S. Belga Fedeli. Random Matrix Methods in Complex Systems Analysis, (2021), the same ensemble was called “the Generalized Wishart.” We, however, feel that the name used in this paper is more appropriate.

33.

V. A.

Marčenko

and

L. A.

Pastur

, “

Distribution of eigenvalues for some sets of random matrices

,”

Math. USSR-Sb.

1

,

457

–

483

(

1967

).

https://doi.org/10.1070/sm1967v001n04abeh001994

Google Scholar

Crossref

34.

E.

Brézin

and

S.

Hikami

, “

Characteristic polynomials of random matrices

,”

Commun. Math. Phys.

214

,

111

–

135

(

2000

).

https://doi.org/10.1007/s002200000256

Google Scholar

Crossref

35.

E.

Brézin

and

S.

Hikami

, “

Characteristic polynomials of real symmetric random matrices

,”

Commun. Math. Phys.

223

,

363

–

382

(

2001

).

https://doi.org/10.1007/s002200100547

Google Scholar

Crossref

36.

Y. V.

Fyodorov

and

E.

Strahov

, “

An exact formula for general spectral correlation function of random Hermitian matrices

,”

J. Phys. A: Math. Gen.

36

(

12

),

3203

–

3214

(

2003

).

https://doi.org/10.1088/0305-4470/36/12/320

Google Scholar

Crossref

37.

E.

Strahov

and

Y. V.

Fyodorov

, “

Universal results for correlations of characteristic polynomials: Riemann-Hilbert approach

,”

Commun. Math. Phys.

241

(

2–3

),

343

–

382

(

2003

).

https://doi.org/10.1007/s00220-003-0938-x

Google Scholar

Crossref

38.

J.

Baik

,

P.

Deift

, and

E.

Strahov

, “

Products and ratios of characteristic polynomials of random Hermitian matrices

,”

J. Math. Phys.

44

,

3657

–

3670

(

2003

).

https://doi.org/10.1063/1.1587875

Google Scholar

Crossref

39.

A.

Borodin

and

E.

Strahov

, “

Averages of characteristic polynomials in random matrix theory

,”

Commun. Pure Appl. Math.

59

(

2

),

161

–

253

(

2006

).

https://doi.org/10.1002/cpa.20092

Google Scholar

Crossref

40.

H.

Köesters

, “

On the second-order correlation function of the characteristic polynomial of a real symmetric Wigner matrix

,”

Electron. Commun. Probab.

13

,

435

–

447

(

2008

).

https://doi.org/10.1214/ecp.v13-1400

Google Scholar

Crossref

41.

H.

Köesters

, “

Characteristic polynomials of sample covariance matrices

,”

J. Theor. Probab.

24

,

545

–

576

(

2011

).

https://doi.org/10.1007/s10959-010-0290-y

Google Scholar

Crossref

42.

T.

Shcherbina

, “

On the correlation function of the characteristic polynomials of the Hermitian Wigner ensemble

,”

Commun. Math. Phys.

308

,

1

–

21

(

2011

).

https://doi.org/10.1007/s00220-011-1316-8

Google Scholar

Crossref

43.

T.

Shcherbina

, “

On the correlation functions of the characteristic polynomials of the Hermitian sample covariance matrices

,”

Probab. Theory Relat. Fields

156

,

449

–

482

(

2013

).

https://doi.org/10.1007/s00440-012-0433-4

Google Scholar

Crossref

44.

I.

Afanasiev

, “

On the correlation functions of the characteristic polynomials of the sparse Hermitian random matrices

,”

J. Stat. Phys.

163

(

2

),

324

–

356

(

2016

).

https://doi.org/10.1007/s10955-016-1486-z

Google Scholar

Crossref

45.

I.

Afanasiev

, “

On the correlation functions of the characteristic polynomials of real random matrices with independent entries

,”

Z. Mat. Fiz., Anal., Geom.

16

(

2

),

91

–

118

(

2020

).

https://doi.org/10.15407/mag16.02.091

Google Scholar

Crossref

46.

D. S.

Dean

and

S. N.

Majumdar

, “

Large deviations of extreme eigenvalues of random matrices

,”

Phys. Rev. Lett.

97

,

160201

(

2006

).

https://doi.org/10.1103/physrevlett.97.160201

Google Scholar

Crossref

PubMed

47.

P.

Littelmann

,

H.-J.

Sommers

, and

M. R.

Zirnbauer

, “

Superbosonization of invariant random matrix ensembles

,”

Commun. Math. Phys.

283

(

2

),

343

–

395

(

2008

).

https://doi.org/10.1007/s00220-008-0535-0

Google Scholar

Crossref

48.

C.

Andréief

, “

Note sur une relation les intégrales définies des produits des fonctions

,”

Mém. Soc. Sci. Bordeaux

2

(

1

),

1

–

14

(

1883

).

Google Scholar

49.

P. J.

Forrester

,

Log-Gases and Random Matrices

(

Princeton University Press

,

2010

).

Google Scholar

Crossref

50.

The Oxford Handbook of Random Matrix Theory

, edited by

G.

Akemann

,

J.

Baik

, and

P.

Di Francesco

(

Oxford University Press

,

2018

), Vol. 1.

Google Scholar

51.

P.

Dita

, “

Parametrisation of unitary matrices

,”

J. Phys. A: Math. Gen.

15

,

3465

–

3473

(

1982

).

https://doi.org/10.1088/0305-4470/15/11/023

Google Scholar

Crossref

2022

Author(s)

Superposition of random plane waves in high spatial dimensions: Random matrix approach to landscape complexity

I. INTRODUCTION

II. MAIN RESULTS

A. Results on the counting problem

B. Results on the general Marchenko–Pastur ensemble

1. Mean asymptotic eigenvalue density for GMPE

2. Moments and correlation functions of characteristic polynomials

III. SOLUTION OF THE ANNEALED COUNTING PROBLEM

IV. SPECTRAL DENSITY OF THE GAUSSIAN MARCHENKO–PASTUR ENSEMBLE

V. MOMENTS AND CORRELATION FUNCTIONS OF CHARACTERISTIC POLYNOMIALS OF MPE AND GMPE

A. Fixed diagonal matrix T

B. Random diagonal matrix T

1. The case p = 1

2. Calculation for p = 2 and β = 1, 2

VI. CONCLUSION

DEDICATION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: LAPLACE APPROXIMATION FOR $Z_{β, 1}$

APPENDIX B: FINITE N EXPRESSION OF $Z_{β, 2}$

APPENDIX C: LAPLACE APPROXIMATION FOR $Z_{β, 2}$

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Superposition of random plane waves in high spatial dimensions: Random matrix approach to landscape complexity

I. INTRODUCTION

II. MAIN RESULTS

A. Results on the counting problem

B. Results on the general Marchenko–Pastur ensemble

1. Mean asymptotic eigenvalue density for GMPE

2. Moments and correlation functions of characteristic polynomials

III. SOLUTION OF THE ANNEALED COUNTING PROBLEM

IV. SPECTRAL DENSITY OF THE GAUSSIAN MARCHENKO–PASTUR ENSEMBLE

V. MOMENTS AND CORRELATION FUNCTIONS OF CHARACTERISTIC POLYNOMIALS OF MPE AND GMPE

A. Fixed diagonal matrix T

B. Random diagonal matrix T

1. The case p = 1

2. Calculation for p = 2 and β = 1, 2

VI. CONCLUSION

DEDICATION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: LAPLACE APPROXIMATION FOR Zβ,1

APPENDIX B: FINITE N EXPRESSION OF Zβ,2

APPENDIX C: LAPLACE APPROXIMATION FOR Zβ,2

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

APPENDIX A: LAPLACE APPROXIMATION FOR $Z_{β, 1}$

APPENDIX B: FINITE N EXPRESSION OF $Z_{β, 2}$

APPENDIX C: LAPLACE APPROXIMATION FOR $Z_{β, 2}$