Motivated by current interest in understanding statistical properties of random landscapes in high-dimensional spaces, we consider a model of the landscape in $RN$ 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

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,

In such a case, the associated Hessian $Wij(x)=\u2212\u2202xi,xjV(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 $\mu /\Gamma \u2032\u2032(0)<1$ to the phase where there is typically only a single minimum for $\mu /\Gamma \u2032\u2032(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,

where $u1,2(a)$ are 2*M* 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*orrandom vectors on the unit

*N*-dimensional sphere. Such a vector**k**_{a}can be built as $kia=(Oa)ii$ 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^{28} 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**^{2}. 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,

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

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 $\cdots T$ as the average with respect to the diagonal matrix *T*, $\cdots K$ as the average with respect to the matrix *K*, and $\cdots $ as the average with respect to both *K* and *T*. The suggested name for this random matrix ensemble^{32} seems appropriate to us as these types of matrices appeared already in the seminal Marchenko–Pastur paper,^{33} 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

where we take the vector ** λ** = (

*λ*

_{1}, …,

*λ*

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

*λ*

_{1}= ⋯ =

*λ*

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

*λ*

_{1}≠

*λ*

_{2}, 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

in terms of a function *m*_{r}(*μ*) ≔ *m*_{r}, which together with its counterpart *m*_{i}(*μ*) ≔ *m*_{i} satisfies the following pair of equations:

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 *μ* → *∞*,

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.

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 $\rho (\lambda )=limN\u2192\u221e1N\u2211k=1N\delta (\lambda \u2212\lambda 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

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 $\rho (\lambda )=I[m(\lambda )]\pi $ where for GMPE, the resolvent *m*(*λ*) is obtained by solving the following transcendental equation:

with $\Psi (u)=\pi 2eu22(1+erf(u/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

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 $\rho (\lambda )\u22481/(2\pi )3/2|\lambda |$.

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: $\alpha /(2\pi (1\u2212\alpha ))$. 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.

#### 2. Moments and correlation functions of characteristic polynomials

For a fixed matrix *T*, the positive moments of the characteristic polynomial $Z\beta ,p(\lambda ;T)\u2261Z\beta ,p(\lambda 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,

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

In these expressions, the functions *g*_{N,m}(*λ*; *T*) and *q*_{N,m} are defined as

where $O$ is the matrix whose elements are all zero and $el(X1,\u2026,Xn)=\u22111\u2264i1<i2<\cdots <il\u2264n\u220fj=1lXij$ is the *l*th 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=\u222bdttn(t)$, where $n(t)\u2254limM\u2192\u221e1M\u2211a=1M\delta (t\u2212Ta)$ is the associated density function. For non-zero $T\u22600$,

in terms of the hypergeometric function defined as $Fqp(a;b;x)=\u2211n=0\u221e\u220fi=1p(ai)n\u220fj=1q(bj)nn!xn$. For $T=0$, we instead find $Z\beta ,1(\lambda )=\lambda N$. Note that in any case, $Z\beta ,1(\lambda )$ is independent of *β*. In the large *N* limit, we study the rate of growth of the mean characteristic polynomial and show that

where we denoted $x=\lambda /T$ and $x\xb1(\alpha )=(1\xb1\alpha )2$ and introduced

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 $T2=\u222bdtt2n(t)$ and reads

In particular, in the large *N* limit, we define

and find that it only depends on the rescaled positive parameter $u=|\lambda 1\lambda 2|/T2\u22650$,

Note a somewhat surprising feature of the two-point function in Marchenko–Pastur ensembles with random matrix *T*: it only depends on the product *λ*_{1}*λ*_{2}, in sharp contrast to earlier studied cases^{34} where in the large-*N* limit, it rather depended on the spectral difference |*λ*_{1} − *λ*_{2}|.

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

where $Wij(x)=\u2212\u2202xi,xjV(x)$ and the expectation $\cdots $ in the case of potential (3) has been taken over both the components of the vector **k**_{a} and the coefficients $u1,2(a)$ for *a* = 1, …, *M*. It turns out to be expedient to pass from the original random variables $u1,2(a),a=1,\u2026,M$, to the new set of random variables obtained by rotation as

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

The random force $fi=\u2212\u2202xiV(x)$ and the Hessian $Wij=\u2212\u2202xi,xjV(x)$ associated with the potential *V*(**x**) are conveniently expressed in terms of these random variables as

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,

which allows us to take explicitly the expectation with respect to random variables *G*_{a}’s, yielding

Further integration over **x** yields

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

where $\cdots 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

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

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

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

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

To this end, we replace the *M*-dimensional integral over components *T*_{a}’s by a functional integral over this density as

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,

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

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

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

where *m*_{r}(*ν*) is the real part of the resolvent *m*(*ν*) at the real spectral parameter *ν*, and we have used that *zm*_{r}(*z*) = *z**dλ**ρ*(*λ*)/(*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,

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

one obtains a simple identity

This identity can, then, be used to obtain

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:

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

with $Z(\mu )=eltot+1/2\pi $ being the *μ*-dependent normalization factor given explicitly by

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}(*μ*),

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

where we have used that *n*_{tot} is the density that minimizes *S*_{tot}. Finally, using again that as *μ* → *∞*, *m*_{r}(*μ*) = *dλρ*(*λ*)/(*μ* − *λ*) → *μ*^{−1} such that *m*_{r}(*μ*) − *μ*^{−1} → 0, one arrives at the final expression for the annealed complexity,

In the limit *μ* → *∞*, we expect that *m*_{r}(*μ*) ∼ *μ*^{−1}, while $mi(\mu )\u221de\u2212\mu 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

In particular, from the second equation, one obtains that

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

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/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(\u2212t2/2)/2\pi $. Computing explicitly the integral over *t*, one obtains the following transcendental equation for *m*(*z*):

where $\Psi (u)=\pi 2eu22(1+erf(u/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

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,

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.^{33} 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/mr+miu/mr3$, one obtains to leading order in *m*_{i},

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

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

As *α* → 1_{+}, one obtains that the value $\rho \alpha \u2248[2\pi (1\u2212\alpha )]\u22121$ 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

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

Let us finally consider the moments $\lambda n=\u222b\u2212\u221e\u221ed\lambda \lambda n\rho (\lambda )$. As *ρ*(*λ*) is even, it trivially yields all odd moments vanishing: $\lambda 2n+1=0$. For any *α*, the even moments can be computed using that for *z* → *∞*,

Using Eq. (44), one obtains that

where $Tk=\u222b\u2212\u221e\u221edttkn(t)$, and we have considered a symmetric density profile *n*(−*t*) = *n*(*t*). In particular, for the Gaussian distribution, $T2n=2n\Gamma (n+1/2)/\pi $ and the lowest moments can be obtained explicitly as

## 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\beta ,p(\lambda ;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\beta ,p(\lambda ;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 ${\psi \u0303j,\psi 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)=\u222bD\psi D\psi \u0303exp\psi \u0303TA\psi =\u222bD\psi D\psi \u0303exp\u2212Tr(A\psi \u0303T\psi )$, yielding

where *Q*^{(s)} = *Q* + *Q*^{T} is the symmetric part of $Q=\u2211j=1p\psi j\psi \u0303jT$, and in the second line, we have used that $\cdots ka$ is the expectation with respect to the Gaussian density $p(ka)\u221de\u2212N\beta ka\u2020ka/2$ for the *a*th column of *K*, which we evaluated in the second line of Eq. (75). Introducing the 2*p* × 2*p* matrix

satisfying $Tr(Q(s))k=\u2212TrQ\u0302k$, one can show that

As the integrand in Eq. (75) only depends on $Q\u0302$, 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*),

where $\Theta a,\beta (U)=det(I(2/\beta )p\u2212TaNU)\beta /2$ and

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

that satisfy

On the other hand, for *β* = 2, the group *C*2*E*(*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:

where we introduced a function

For *β* = 2, using that $det1\u2264i,j\u2264p(xij\u22121)=\u220fi<jp(xi\u2212xj)$ together with the Cauchy–Binet–Andréief formula,^{48} Eq. (82) can be re-written as the ratio of determinants,

where

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

where *e*_{l}(*X*) is the *l*th elementary symmetric polynomial defined as

Similarly, for *β* = 1, using the identity

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

In Fig. 3, we have compared our exact analytical results for $Z\beta ,p(\lambda ;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.

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

with

where the hypergeometric function of matrix argument is defined as

with *κ* standing for a partition of *k* and $(a)\kappa (\alpha )$ being the generalized Pochhammer symbol associated with the partition *κ* = (*κ*_{1}, …, *κ*_{m}),

and the function $C\kappa (\alpha )(X)$ is the so-called *C*-normalization of the Jack polynomial.^{49} This allows us to express the function $Z\beta ,p(0;T)$ via a ratio of the hypergeometric functions as

### B. Random diagonal matrix *T*

Let us now consider the case of random diagonal matrix *T* and average the *p*-point characteristic polynomial $Z\beta ,p(\lambda ;T)$ over realizations of *T*, denoting these expectations as $Z\beta ,p(\lambda )=Z\beta ,p(\lambda ;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\beta ,p(\lambda )$ depends on the whole distribution of *T* only via the lowest moments $Tk=\u222b\u2212\u221e\u221edttkn(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

In the particular case of a symmetric distribution *n*(*t*) = *n*(−*t*), we have $T=0$ and (97) simplifies drastically: $Z\beta ,1(\lambda )=det(\lambda I\u2212W)=\lambda N$. If $T\u22600$ instead and assuming $|T|<N$, $Z\beta ,1(\lambda )$ is expressed for finite *N* as

where we remind the definition of the standard hypergeometric function,

with $(a)n=\u220fi=0n\u22121(a+i)$ being the Pochhammer symbol.

Note that the ratio $Z\beta ,1(\lambda )/\lambda N$ only depends on the rescaled variable $x=\lambda /T$. In the limit *N* → *∞* and for $T\u22600$, changing the variable from *z* to $w=(T/N)z$ in Eq. (97), one obtains

which can, then, be conveniently analyzed using Laplace’s method. The saddle points of $L1(w;x,\alpha )$ in the *w* domain are given by $w\xb1(x,\alpha )=(1\u2212\alpha +x\xb1(\alpha \u22121\u2212x)2\u22124x)/(2x)$ and are complex conjugate for (*α* − 1 − *x*)^{2} < 4*x*, i.e., *x* ∈ [*x*_{−}(*α*), *x*_{+}(*α*)], where $x\xb1(\alpha )=(1\xb1\alpha )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\xb1(\alpha )=(1\xb1\alpha )2$. In the large *N* limit, a careful analysis of Eq. (100) (see Appendix A for details) allows us to obtain

The function Ξ_{1}(*x*; *α*) changes the sign over the real interval and behaves in the limit *x* → ±*∞* as

which matches smoothly the expression $ln|Z\beta ,1(\lambda )/\lambda 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 Ξ_{1}(*x*, *α*) for *α* = 2 and the finite *N* = 20 expression of $N\u22121\u2061ln|Z1,1(x)/xN|$ with $T=1$, showing excellent agreement.

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

In order to obtain an exact finite *N* expression of the two-point correlation function $Z\beta ,2(\lambda )=Z\beta ,2(\lambda ;T)T$ for ** λ** = (

*λ*

_{1},

*λ*

_{2}), 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:

Note that $Z\beta ,2(\lambda )$ only depends on ** λ** = (

*λ*

_{1},

*λ*

_{2}) through the product

*λ*

_{1}

*λ*

_{2}, and the ratio $Z\beta ,2(\lambda =(\lambda 1,\lambda 2))/(\lambda 1\lambda 2)N$ only depends on the rescaled variable $y=\lambda 1\lambda 2/T2$. This is highly unusual as similar objects for Gaussian or Wishart matrices become in the large-

*N*limit a function of the difference |

*λ*

_{1}−

*λ*

_{2}|; 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.

Let us consider the large *N* limit of Eq. (104). To this purpose, we re-write the ratio of Eq. (104) by $(\lambda 1\lambda 2)N$ as the finite sum,

where $y=\lambda 1\lambda 2/T2$. 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

where (see Appendix C for details)

The saddle-point equation yields a third order algebraic equation for the value *η* that maximizes the action $L2(\eta ;u,\alpha )$,

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

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

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

On the other hand, one can obtain that the function $Z\beta ,2(\lambda )$ achieves a value exponentially large in *N* for *λ*_{1}*λ*_{2} = 0, with

In Fig. 6, we show a comparison between the expression of Ξ_{2}(*u*; *α*) in Eq. (112) and numerical simulations of $N\u22121\u2061ln|Z1,2(u,u)/uN|$ and $N\u22121\u2061ln|Z1,2(1,u)/uN|$. In the large *N* limit, it is expected that $limN\u2192\u221eN\u22121\u2061ln|Z1,2(\lambda 1,\lambda 2)/(\lambda 1\lambda 2)N|$ only depends on the rescaled parameter *u* = |*λ*_{1}*λ*_{2}| 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*.

## 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,^{14} 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\beta ,1$

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

where $x=\lambda /T$. The saddle points of $L1(w;x,\alpha )$ are given by

The contribution of each of them depends on the sign of Δ_{1} = (*α* − 1 − *x*)^{2} − 4*x*.

For Δ_{1} > 0, corresponding either to *x* < *x*_{−}(*α*) or *x* > *x*_{+}(*α*) with $x\xb1(\alpha )=(1\xb1\alpha )2$, the saddle points *w*_{±} both lie on the real line. In order to obtain the stability of the fixed points, let us compute

which is real and positive only for *w*_{±} real and inside the interval between its two roots $((1\u2212\alpha )\u22121,(1+\alpha )\u22121)$. 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\u2212\u2208((1\u2212\alpha )\u22121,(1+\alpha )\u22121)$ (and *w*_{+} outside of this interval), and conversely for *x* < *x*_{−}(*α*), one has $w+\u2208((1\u2212\alpha )\u22121,(1+\alpha )\u22121)$ (and *w*_{−} outside).

Conversely, for Δ_{1} < 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)^{2}/(*α* + 1).

Summarizing, one obtains for *N* ≫ 1,

Finally, for *x* = *x*_{±}(*α*), the saddle points are unique and of order 2:

- For
*x*=*x*_{+}(*α*), one has that $w+=(1+\alpha )\u22121$ such thatTherefore, the steepest descent directions at the point(A7)$\u2202w2L1(w;x,\alpha )w=w+=0,\u2202w3L1(w;x,\alpha )w=w+=\u22122(1+\alpha )4\alpha .$*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π/3}*t*+*w*_{+}and*w*(*t*) =*e*^{i2π/3}*t*+*w*_{+}, respectively, yielding(A8)$Z\beta ,1(\lambda )\lambda N=eNL1(w+;x,\alpha )2i\pi w+ei2\pi 3\u222b0+\u221edte\u2212N2(1+\alpha )43!\alpha t3+ei4\pi 3\u222b+\u221e0dte\u2212N2(1+\alpha )43!\alpha t3=\Gamma (1/3)2\pi 31/6N(1+\alpha )\alpha 1/3eNL1(w+;x,\alpha ).$ - For
*x*=*x*_{−}(*α*), $w\u2212=(1\u2212\alpha )\u22121$ such thatThe steepest descent directions at(A9)$\u2202w2L1(w;x,\alpha )w=w\u2212=0,\u2202w3L1(w;x,\alpha )w=w\u2212=2(1\u2212\alpha )4\alpha .$*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, yieldingIn each case, one obtains that the term(A10)$Z\beta ,1(\lambda )\lambda N=eNL1(w\u2212;x,\alpha )2i\pi w\u2212ei\pi 3\u222b0+\u221edte\u2212N2(1\u2212\alpha )43!\alpha t3+ei5\pi 3\u222b+\u221e0dte\u2212N2(1\u2212\alpha )43!\alpha t3=\Gamma (1/3)(\u22121)N2\pi 31/6N(\alpha \u22121)\alpha 1/3eNL1(w\u2212;x,\alpha ).$where $x=\lambda /T$.(A11)$\Xi 1(x=x\xb1(\alpha );\alpha )=limN\u2192+\u221e1NlnZ\beta ,1(\lambda )\lambda N=L1(w\xb1;x\xb1(\alpha ),\alpha ),$

### APPENDIX B: FINITE *N* EXPRESSION OF $Z\beta ,2$

Let us first consider *β* = 2. In this case, Λ = diag(*λ*_{1}, *λ*_{2}), 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:^{51}

The Euclidean line element (d*s*)^{2} = *∑*_{i,j}|d*U*_{ij}|^{2} reads in this parametrization

The associated Jacobian is given by $2detM$, where

The normalized Haar measure is obtained calculating $\u222bdetMd\psi 11d\psi 12d\psi 22d\varphi $ and is given by $dU=14\pi 3cos\varphi sin\varphi d\psi 11d\psi 12d\psi 22d\varphi $. With the further assumption that $T=0$, from what stated above follows

The denominator of Eq. (78) reads

while the numerator reads

For *β* = 1, we impose $T2=1$ to simplify the computations. The general case is retrieved by multiplying $Z2,\beta =1(\lambda )$ by $T2N$ and rescaling $\lambda 1,2\u2192\lambda 1,2/T2$. We rename Λ = diag(*λ*_{1}, *λ*_{2}, *λ*_{1}, *λ*_{2}), while now the matrix *U* is parameterized with a 2 × 2 block matrix. Finally, we introduce six variables: *h*_{1}, *h*_{2} ∈ (0, 1) and *ϕ*_{11}, *ϕ*_{21}, *ϕ*_{22}, *ϕ*_{14} ∈ (0, 2*π*). Therefore (see Ref. 51),

with

and

Similarly to Eq. (B2), in order to retrieve the normalized Haar measure, the Euclidean line element is (d*s*)^{2} = *∑*_{i,j}|d*U*_{ij}|^{2} = 2*∑*_{i,j}|d*U*_{11;i,j}|^{2} + *∑*_{i,j}|d*U*_{12;i,j}|^{2} + *∑*_{i,j}|d*U*_{21;i,j}|^{2}. In detail, we have that

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*_{1} and *A*_{2}. The latter are

and

We obtain the Haar measure by calculating

Integrating out the six linear and angular variables, one obtains

The denominator appearing in Eq. (78) is

whereas the numerator in Eq. (78) reads

The angular variables *ϕ*_{22} and *ϕ*_{21} in the integral above can be integrated out with the auxiliary variable *ϕ* = *ϕ*_{22} + *ϕ*_{21}. The integral appearing in Eq. (B21) is simplified to

Working out the remaining integration, the expression above becomes

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

### APPENDIX C: LAPLACE APPROXIMATION FOR $Z\beta ,2$

We consider the case of centered random variables, i.e., $Ti=0$. In order to obtain the asymptotics of $Z\beta ,2$, it is convenient to start from Eq. (104) and to introduce $x=\lambda 1\lambda 2/T2$. Expanding the terms within the hypergeometric function, we note

where $(b)k=\Gamma (b+k)\Gamma (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

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

Hence, we need to retrieve the roots of *p*(*η*) = *η*^{2}(*α* − 1 + *η*) − *x*(1 − *η*). First, one observes that lim_{η→±∞}*p*(*η*) = ±*∞* and *p*(0) = −*x*. For *η* > 0, *p*(*η*) is strictly increasing since *η*^{2} + 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

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

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

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

We may now in each term replace the sum over *k* with an integral over, respectively, *η*_{e} = 2*k*/*N* and *η*_{o} = (2*k* + 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

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.