Superposition of Random Plane Waves in High Spatial Dimensions: Random Matrix Approach to Landscape Complexity

Motivated by current interest in understanding statistical properties of random landscapes in high-dimensional spaces, we consider a model of the landscape in $\mathbb{R}^N$ obtained by superimposing $M>N$ plane waves of random wavevectors and amplitudes. 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\to \infty$ at fixed ratio $\alpha=M/N>1$. The framework of this computation requires us to study spectral properties of $N\times N$ matrices $W=KTK^T$, where $T$ is diagonal with $M$ mean zero 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 and related matrices.


I. INTRODUCTION
The problem of characterising the statistical properties of high-dimensional random functions, frequently called "random landscapes", originated in the theory of disordered systems such as spin glasses, see [1] for an accessible introduction, and gradually became popular beyond the original setting, finding numerous applications in such diverse fields as cosmology [2,3], machine learning via deep neural networks [4,5], and large-size inference problems in statistics [6][7][8][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, 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][11][12][13][14], and became a focus of vigorous activity in recent years, see [15][16][17][18][19][20][21][22][23][24] and 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 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 W i j (x) = −∂ 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 appear also 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 allowsed 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 certain class of potentials V (x) [18,22], though in general they may differ.
The simplest, yet informative characteristics of the landscape is just the total annealed complexity of all possible stationary points. In particular, studying such complexity revealed that the invariant landscapes (1)-(2) undergo a "topology trivialisation transition" from the glassy phase with exponentially many stationary points for small values of the ratio µ/ Γ (0) < 1, to the phase where there is typically only a single minimum for µ/ Γ (0) > 1. Similar transitions are also operative in random potentials confined to high-dimensional spheres [7,[25][26][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 the u (a) 1,2 are M i.i.d. normal (mean zero, variance one) random variables and the M vectors k a = (k 1,a , · · · , k N,a ) are all i.i.d., with N real independent, normally distributed (mean zero, variance 1/N) random components k i,a , i = 1, . . . , 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 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 stationary points in not unrelated random superpositions of plane waves in low spatial dimension N = 2 is currently under intensive study as a universal model for high-energy eigenfunctions of the Laplace operator on generic compact Riemannian manifolds, see [28] and references therein.
Finally, for low values of M, N the model (3) describes periodic systems such as crystals in 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. [29,30].
We believe that the suggested landscape in high dimensions represents a new interesting type worth consideration. In particular 1. such landscapes still allow 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 µ.
2. In contrast to the 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 = KT K T and turns out to belong to the "Gaussian Marchenko-Pastur Ensemble" (GMPE) that 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 · · · T the average with respect to the diagonal matrix T , · · · K the average with respect to the matrix K and · · · the average with respect to both K and T . The suggested name for this random matrix ensemble [31] seems appropriate to us as this type of matrices appeared already in the seminal Marchenko-Pastur paper [32], though does 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 to 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 the resulting ensemble GMPE is interesting in its own right and we further study its mean eigenvalue density as N → ∞ as well as 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 [33]. 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 . We will provide some explicit results for the 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 [34][35][36][37][38][39][40][41][42][43][44][45].
The rest of the paper is organised as follows. In section II, we summarise the main results obtained in this paper. In section III, we consider the total annealed complexity of the landscape (1) with the random potential defined in (3). In the following sections, we characterise the Gaussian Marchenko-Pastur Ensemble associated to this counting problem by computing its mean spectral density in section IV and correlation function of its characteristic polynomial defined in (6) in section V. Finally, in section VI, we discuss and interpret our results and give some perspective on further investigation of this random landscape. Some details of the computations are relegated to Appendices A, B and C.

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 the landscape (1) with the random potential (3) can be expressed as FIG. 1. Comparison between the annealed total complexity Ξ tot (µ; α) for α = 2 (black solid line) obtained by inserting the solution of Eqs. (8)(9) into Eq.(7) and numerical simulations of (orange triangles) for N = 200, plotted as a function of µ. It is expected that these three quantities coincide in the N → ∞ limit. The agreement is good for low values of N but there is some discrepancy 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. in terms of a function m r (µ) := m r which together with its counterpart m i (µ) := m i satisfies the 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 trivialisation transition" and stays positive for any finite value 0 < µ < ∞ of the strength of the confinement.
Still, it decreases very rapidly as µ → ∞: 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 [12] and related studies. Mean asymptotic eigenvalue density for GMPE: Consider the mean spectral density ρ(λ ) = lim N→∞ 1 N ∑ N k=1 δ (λ − λ 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 ρ(λ ) = ℑ[m(λ )] π where for GMPE the resolvent m(λ ) is obtained by solving the following transcendental equation: with Ψ(u) = π 2 e u 2 2 (1 + erf(u/ √ 2)). The resulting mean spectral density can be shown to be even: ρ(λ ) = ρ(−λ ) and supported on the whole real line, with asymptotic behaviour 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 ρ(λ ) ≈ 1/ (2π) 3/2 |λ |.
Let us finally mention that for α < 1 the density ρ(λ ) displays a delta peak at the origin part tends to a finite value at the origin: α/( √ 2π(1 −α)). In Fig. 2, we plot a comparison between our analytical prediction for the average density of eigenvalue ρ(λ ) in the N → ∞ limit and results from numerical simulation for N = 500, showing excellent agreement.

Moments and correlation functions of characteristic polynomials:
For a fixed matrix T the positive moments of the characteristic polynomial Z β ,p (λ ; T ) ≡ Z β ,p (λ 1; T ) defined in Eq.(6) with setting λ = λ 1 = λ (1, · · · , 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 2p × 2p Pfaffian for β = 1 In these expressions, the functions g N,m (λ ; T ) and q N,m are defined as where e l (X 1 , · · · , X n ) = ∑ 1≤i 1 <i 2 <···<i l ≤n ∏ l j=1 X i j 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. [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 = in terms of the hypergeometric function defined as In the large N limit, we study the rate of growth of the mean characteristic polynomial and show that where 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 = dt t 2 n(t) and reads In particular in the large N limit we define and find that it only depends on the rescaled positive parameter u = |λ 1 λ 2 |/ T 2 ≥ 0: Note a somewhat surprising feature of the two-point function in the Marchenko-Pastur ensembles with random matrix T : it depends only on the product λ 1 λ 2 , in sharp contrast with earlier studied cases [34] where in the large-N limit it rather depended on the spectral difference |λ 1 − λ 2 |.

ACKNOWLEDGMENTS
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 subsequent development of the field.
We would like to thank Pierre Le Doussal for his interest in the work, and in particular for EP/L015854/1.

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 the case (1) reads where W i j (x) = −∂ x i ,x j V (x) and the expectation · · · in the case of potential (3) has been taken over both the components of the vector k a , and the coefficients u (a) 1,2 for a = 1, · · · , M. It turns out to be expedient to pass from the original random variables u (a) 1,2 , a = 1, · · · , 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 dependent, the resulting random variables are statistically x independent, with zero mean and covariance The random force 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 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 to take explicitly the expectation with respect to the 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 · · · K is the remaining expectation over the Gaussian 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 = KT K 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 the λ i 's being the random eigenvalues of W = KT K T . Note that in the model with 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 asymptics we find it expedient to introduce the normalised limiting density of T a as and replace the M-dimensional integral over the 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.
This form is now amenable to using the Laplace method justified by the large parameter N, yielding for the annealed complexity where we have denoted n den , n tot and l tot , l den the values of the limiting density and Lagrange parameters that minimise S den , S tot respectively. 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 [12], however instead of maximising with respect to a single real variable one has to optimize with respect to the limiting 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 it can be recast rigorously using the Large Deviations framework.
Let us first consider the stationarity conditions for the action functional S den . Requiring its variation with respect to n(t) and the derivative with respect to the Lagrange parameter to be zero, one obtains the equations One can check easily that the distribution n den (t) is then given by the properly normalized Gaussian for which S den [n den ; Similarly, 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 normalisation 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 density ρ(λ ) satisfies for a fixed limiting density n(t) a self-consistency equation derived by Marchenko and Pastur in [32] which reads The resolvent can be introduced in the saddle-point equation (47) using the identity, where m r (ν) is the real part of the resolvent m(ν) at the real spectral parameter ν and we have used that z m r (z) = z− ∞ −∞ dλ ρ(λ )/(z − λ ) → 1 for any real or complex spectral parameter z in the limit |z| → ∞. Using the self-consistency equation (48), 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 (47) implies that the optimal density n tot (t) satisfies the equation which is immediately solved in terms of the real and imaginary parts of the resolvent at position µ with Z(µ) = e l tot +1 / √ 2π = e l tot −l den being the µ−dependent normalisation factor given explicitly by Finally, at a real value of µ the self-consistency equation (48) 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 term 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 density ρ(λ ) of W = KT K T is also unbounded, see [32]. Hence m i (µ) > 0 for any 0 < µ < ∞ and one can divide by m i on both sides of Eq. (58). Inserting then the expression of n tot (t) one obtains the more explicit equations (8)(9). This pair of equations allows to express m i (µ) and m r (µ) at the optimal solution. One now needs to relate these values to the complexity expression (43) 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 minimises 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 m i (µ) ∝ 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. (56) and (57) 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 behaviour is in stark contrast with previously considered models with isotropic confinement, where the complexity was found to vanish for µ > µ c , signalling 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 where Ψ(u) = π 2 e u 2 2 (1 + erf(u/ √ 2)). Taking the paramter z = λ on the real axis, we introduce as in the previous section 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 equation (65) allows an efficient numerical evaluation of the density at a given real position z = λ , the properties of the mean density ρ(λ ) are most easily analysed using the set of equations for the real and imaginary part of the resolvent: For an unbounded limiting density n(t) of T , such as the Gaussian density in GMPE, the associated mean spectral density ρ(λ ) for matrices W = KT K T is also unbounded [32]. In that limit, one expects that m i (λ ) = πρ(λ ) ∝ n(λ ) m r (λ ) as λ → ∞. Using this approximation and introducing in the integrals of Eqs. (67) and (68) the change of variable t = 1/m r + m i u/m 3 r , one obtains to leading order in m i In particular, for a Gaussian density n(t), it yields the asymptotic behaviour For α > 1, the density ρ(λ = 0) = ρ α is finite and in the present case it is obtained by setting z = 0 and m(z) = iπρ α in Eq. (65), yielding As α → 1 + , one obtains that the value ρ α ≈ [ √ 2π(1 − α)] −1 diverges.
For α < 1 instead, the mean 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 Finally, in the special case α = 1, the density is continuous but diverges on approaching the origin as Let us finally consider the moments λ n = ∞ −∞ dλ λ n ρ(λ ). As ρ(λ ) is even, it trivially yields all odd moments vanishing: λ 2n+1 = 0. For any α, the even moments can be computed using that for z → ∞, Using Eq. (48), one obtains that where T k = ∞ −∞ dt t k n(t) and we have considered a symmetric distribution n(−t) = n(t). In particular for the Gaussian distribution T 2n = 2 n Γ(n + 1/2)/ √ π 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 descibe 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 a lot of 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 which 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 where Q (s) = Q + Q T is the symmetric part of Q = ∑ p j=1 ψ jψ T j , and in the second line we have used that · · · k a is the expectation with respect to the Gaussian density p(k a ) ∝ e −Nβ k † a k a /2 for a th column of K, which we evaluated in the third line, Eq.(79). Introducing the 2p × 2p matrix satisfying Tr(Q (s) ) k = − TrQ k , one can show that As the integrand in Eq. (79) only depends onQ one can use the so-called "bosonization" trick proposed in [46] and replace the integral over anticommuting variables by the invariant integration over unitary matrices U belonging to the group C(4/β )E(p): where Θ a,β (U) = det(I (2/β )p − T a N U) β /2 and 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 that satisfies 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 = λ (1, · · · , 1). In that case, both integrands in the numerator and the denominator of Eq.(82) 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 det 1≤i, j≤p (x together with the Cauchy-Binet-Andréief formula [47], Eq.(86) 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 the 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 β ,p (λ ; T ) as a function of λ given in Eq. (94) and (88) respectively for β = 1, 2 taking fixed diagonal matrix T with numerical simulations, obtaining an excellent agreement.
Let us finally mention that in the case λ = 0 one can use the following identities from the book with and where the hypergeometric function of matrix argument is defined as with κ standing for a partition of k, (a) (α) κ being the generalised Pochhammer symbol associated to the partition κ = (κ 1 , · · · , κ m ), and the function C κ (X) is the so-called C-normalisation of Jack polynomial [48]. This allows to express the function Z β ,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 β ,p (λ; T ) over realisations 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 formulae are universal in the sense that Z β ,p (λ) depends on the whole distribution of T only via the lowest moments T k = ∞ −∞ dt t k n(t), k = 1, 2 of the i.i.d. elements of T , assuming those moments are finite. In such a case we can directly use the results in Eq. (88) and (94) 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 (101) simplifies drastically: Z β ,1 (λ ) = det(λ I −W ) = λ N . If T = 0 instead and assuming | T | < N, where we remind the definition of the standard hypergeometric funtion with (a) n = ∏ n−1 i=0 (a + i) 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 = 0, changing variable from z to w = ( T /N)z in Eq.(101), one obtains which can then be conveniently analysed using Laplace's method. The saddle points of L 1 (w; x, α) in the w domain are given by w ± (x, α) ] where x ± (α) = (1 ± √ α) 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.
Let us consider the large N limit of Eq. (108). To this purpose, we re-write the ratio of Eq.(108) by (λ 1 λ 2 ) N as the finite sum 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 behaviour 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 maximises the action L 2 (η; u, α) where 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.(63), i.e.
On the other hand, one can obtain that the function Z β ,2 (λ) achieves a value exponentially large In Fig. 6, we show a comparison between the expression of Ξ 2 (u; α) in Eq. (116) and numerical simulations of N −1 ln |Z 1,2 ( √ u, √ 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→∞ N −1 ln |Z 1,2 (λ 1 , λ 2 )/(λ 1 λ 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 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. [12], this system does not display a "topological trivialisation transition", i.e. the total annealed complexity is strictly positive for any finite µ.
The random part of the Jacobian associated to such a random landscape defines a random matrix W = KT K T belonging to a new random matrix ensemble that we denote "Gaussian Marchencko 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,49].
The unforeseen behaviour of the total annealed complexity naturally invites further study of the Gibbs measure associated to 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 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, as well as on 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.
Appendix A: Laplace approximation for Z β ,1 Let us start from Eq.(104) that we reproduce here where x = λ / T . The saddle points of L 1 (w; x, α) are given by The contribution of each of them depends on the sign of 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 − √ α) −1 , (1 + √ α) −1 ). These two roots translate in terms of x as the two edges x ± (α) of the usual Marchencko Pastur distribution and one can check that for x > x + (α), one has w − ∈ ((1 − √ α) −1 , (1 + √ α) −1 ) (and w + outside of this interval) and conversely for x < x − (α), one has 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 sign for x = (α − 1) 2 /(α + 1).
For β = 1, we impose T 2 = 1 to simplify the computations. The general case is retrieved by mul- while now the matrix U is parametrized with a 2×2 block matrix. Lastly, we introduce 6 variables: h 1 , h 2 ∈ (0, 1) and φ 11 , φ 21 , φ 22 , φ 14 ∈ (0, 2π). Therefore (see [50]), and Similarly to Eq.(B2), in order to retrieve the normalized Haar measure, the Eucledean line element In details we have that: To simplify the calculation, we can re-arrange the differential of the 6 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 Working out the remaining integration, the expression above becomes and Eq(108) is retrieved with the introduction of the hypergeometric function.
Appendix C: Laplace approximation for Z β ,2 We consider the case of centred random variables, i.e. T i = 0. In order to obtain the asymptotics of Z β ,2 , it is convenient to start from Eq.(108) and to introduce x = λ 1 λ 2 / T 2 . Expanding the terms within the hypergeometric function, we notice where (b) k = Γ(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 The maximum of L 2 is located in η * ∈ [0, 1]. The latter is the solution of the equation ∂ η L 2 η=η * = 0 , ⇒ ln(x − η * x) − ln(α + η * − 1) − 2 ln(η * ) = 0 .
Hence, we need to retrieve the roots of p(η) = η 2 (α − 1 + η) − x(1 − η). Firstly, one observes 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 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 = 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 substract 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.