Manifolds pinned by a high-dimensional random landscape: Hessian at the global energy minimum

We consider an elastic manifold of internal dimension $d$ and length $L$ pinned in a $N$ dimensional random potential and confined by an additional parabolic potential of curvature $\mu$. We are interested in the mean spectral density $\rho(\lambda)$ of the Hessian matrix $K$ at the absolute minimum of the total energy. We use the replica approach to derive the system of equations for $\rho(\lambda)$ for a fixed $L^d$ in the $N \to \infty$ limit extending $d=0$ results of our previous work. A particular attention is devoted to analyzing the limit of extended lattice systems by letting $L\to \infty$. In all cases we show that for a confinement curvature $\mu$ exceeding a critical value $\mu_c$, the so-called"Larkin mass", the system is replica-symmetric and the Hessian spectrum is always gapped (from zero). The gap vanishes quadratically at $\mu\to \mu_c$. For $\mu<\mu_c$ the replica symmetry breaking (RSB) occurs and the Hessian spectrum is either gapped or extends down to zero, depending on whether RSB is 1-step or full. In the 1-RSB case the gap vanishes in all $d$ as $(\mu_c-\mu)^4$ near the transition. In the full RSB case the gap is identically zero. A set of specific landscapes realize the so-called"marginal cases"in $d=1,2$ which share both feature of the 1-step and the full RSB solution, and exhibit some scale invariance. We also obtain the average Green function associated to the Hessian and find that at the edge of the spectrum it decays exponentially in the distance within the internal space of the manifold with a length scale equal in all cases to the Larkin length introduced in the theory of pinning.


The random manifold model and some known results
Numerous physical systems can be modeled by a collection of points or particles coupled by an elastic energy, usually called an elastic manifold, submitted to a random potential (see [1][2][3] for reviews). They are often called "disordered elastic systems" and generically exhibit pinning in their statics and depinning transitions and avalanches in their driven dynamics [4][5][6][7][8][9]. Their energy landscape is complex leading to glassy behavior.
The manifold is usually parameterized by a N -component real displacement field u(x) ∈ R N , where x belongs to an internal space x ∈ Ω. Ω can be either a finite collection of points, such as a subset L d of an internal space of dimension d, Ω ⊂ Z d , for discrete models, or Ω ⊂ R d in a continous setting. The case d = 1 corresponds to a line in N dimensions and for N = 1 was studied in the present context in [10]. The case d = 0 usually refers below to Ω being a single point, previously studied in [11] in the large N limit, and the present study can be seen as its generalization to a manifold. There are two terms in the total energy. First the points in Ω are coupled via an elastic energy, which is a quadratic form in the fields u(x). We also include in this quadratic term a parabolic confining potential of curvature µ > 0. The absolute minimum of this first term is thus the flat, undisturbed, configuration u(x) = 0. The second term is the quenched disorder, modeled by a random potential energy which couples directly to u(x). We thus consider the following model of an elastic manifold in a random potential given by its energy functional where here x ∈ L d ⊂ Z d , 1 is an appropriate identity operator, and the matrix −t∆ xy is required to be positive definite. Here ∆ can be chosen as the discrete Laplacian in the hypercube L d with periodic boundary conditions. In that case its eigenmodes are plane waves ∼ e ikx and we denote ∆(k) its eigenvalues, i.e. in d = 1, ∆(k) = 2(cos k − 1) with k = 2πn/L, n = 0, ..L − 1. For general d similar formula holds and t must be positive, t > 0. All formula below extend immediately to more general functions t∆(k), e.g. to more general elasticity (such as long range elasticity). They also extend to cases where t∆ xy is a quadratic form defined on any graph Ω. Finally, they also extend to the limit of the continuum manifold model, e.g. with the standard Laplacian ∆ = d i=1 ∂ 2 ∂x 2 i whose spectrum is given by ∆(k) = −k 2 . We thus use the notation k = 1 (2π) d so our main formula are valid both for discrete and continuum (in the continuum x ≡ d d x). We see from (1) that µ acts as a "mass" which, for the continuum model, leads to reducing the fluctuations beyond the scale L µ = t/µ.
Here we will consider V (u, x) to be a mean-zero Gaussian-distributed random potential in R N × Z d with a rotational and translational invariant covariance (also called the correlator in the physical context) such that potential values are uncorrelated for different points in the internal space, but correlated for different displacements: ‡ In Eq.(2) and henceforth the notation · · · stands for the quantities averaged over the random potential.
The equilibrium statics of this model has been much studied. From the competition between the elastic and the disorder energy, the minimal energy configuration u 0 (x) (ground state) is non trivial and exhibits interesting statistically self-affine properties characterized by a roughness exponent: u 0 (x) − u 0 (0) ∼ |x| ζ . The sample to sample fluctuations of the ground state energy (and, at finite temperature, of the free energy) grow with the scale as ∼ L θ , with θ = d − 2 + 2ζ as a consequence of the symmetries of the model (1)- (2). In addition the manifold is pinned, i.e. its macroscopic response to an external force is non-linear. The early (and partly phenomenological) theory of pinning is due to Larkin and Ovchinnikov (see [1] for a review). Below the so-called Larkin length scale L c , with L c ∼ (B (0)) −1/(4−d) for weak disorder (small B (0)) the deformations are elastic, the response is linear, deformations can be calculated from perturbation theory, leading to the roughness exponent ζ = (4−d)/2. Above L c metastability sets in, the response to perturbations involves jumps (shocks), with non-trivial roughness ζ of minimal energy configurations. Describing that regime has been a challenge, and progress was later achieved using the bag of tools of the statistical mechanics of disordered systems, most notably replica methods. Exact results have been obtained, but in only a few analytically tractable cases. The first of such cases are mean field type models, notably the model (1) in the limit N → ∞. Saddle point equations in replica space [14][15][16][17] lead to solutions exhibiting replica symmetry breaking (RSB) for µ < µ c , which describe the glass phase where the manifold is pinned. The critical mass µ c corresponds to the Larkin scale L c = t/µ c and the glass phase appears at scales exceeding L c . A second set of results were obtained using the functional renormalization group [2,[18][19][20] and are valid in an expansion around d = 4 − (for any N ). While the resulting physical picture is somewhat different, these could be reconciled [14,15]. Note also that the Larkin picture was fully confirmed by these studies. Finally, for d = 1 the problem can be mapped to stirred Burgers and Kardar-Parisi-Zhang growth (see Ref. [21] for review of earlier works). For N = 1, a number of exact results were obtained recently from an emerging integrability structure of the theory, both in physics and mathematics. Besides proving the exact roughness ζ = 2/3 and free energy fluctuation exponent, θ = 1/3, it was shown, e.g., that the probability density of the free energy for a long polymer converges to the famous Tracy-Widom distribution both at zero temperature [22], and at finite temperature in the continuum [23][24][25][26][27]. Finally, note that the model (1)- (2) also arises in the study of the decaying Burgers equation with a random initial conditions in dimension N , which exhibits interesting transitions and regimes, see e.g. for N = 1 [28] and for large N [29].

Motivation and goals of the paper
While these results predict large scale properties of the low energy configurations, little is known about the detailed statistical structure of the complex energy landscape of pinned manifolds. This relates to the broad effort of understanding the statistical structure of stationary points (minima, maxima and saddles) of random landscapes which is of steady interest in theoretical physics [30][31][32][33][34][35][36][37][38][39], with recent applications to statistical physics [10,[34][35][36][38][39][40][41], neural networks and complex dynamics [42][43][44][45][46], string theory [47,48] and cosmology [49,50]. It is also of active current interest in pure and applied mathematics [51][52][53][54][55][56][57][58][59][60], For the model (1)- (2) in the simplest case d = 0 (x is a single point), the mean number of stationary points and of minima of the energy function was investigated in the limit of large N 1 in [35,38,39], see also [37,50,52]. It was found that a sharp transition occurs from a 'simple' landscape for µ > µ c (the same µ c as given by the onset of RSB, see above), with typically only a single stationary point (the minimum) to a complex ('glassy') landscapes for µ < µ c with exponentially many stationary points. Similar transitions were found in related systems upon applying various external perturbations [40,41,44] in particular in the mean number of stationary points which was also studied recently for the case of an elastic string d = 1 in dimension N = 1 [10]. Relations with Anderson localization was discussed there in this context.
An important quantity which characterizes the stability of local equilibria, and is crucial both for equilibrium and slowly driven dynamics, is the Hessian matrix. In particular, the question of whether the spectrum of the Hessian at low lying local minima is gapped (away from zero) or not, the behavior of its mean density of eigenvalues near zero, and the nature of the associated low lying modes, has been identified as a crucial feature to describe classical [61][62][63] and quantum glasses [64][65][66][67][68]. Clearly, a 'gapless' spectrum reflects the existence of very 'flat' directions in configuration space along which moving away from the local minimum incurs very little 'cost'. This flatness, also known as a 'marginal stability', is ubiquitous in various types of glasses [62,63] and appears naturally in models exhibiting a hierarchical structure of the energy landscapes [69,70]. The Hessian matrix was studied recently numerically in the context of the depinning of an elastic line d = 1 in a one dimension random potential, N = 1, in an effort to identify the "soft modes" which trigger the avalanches. It was found that in the stationary state reached upon quasi-static driving, the low-lying modes of the Hessian are localized, with a localization length directly related to the Larkin pinning length [71]. Although studying the Hessian at equilibrium, and specifically at the global minimum would be also very interesting, it is analytically challenging for d or N small.
Recently, by combining methods of random matrix theory with methods of statistical mechanics of disordered systems, we were able to study the Hessian at the absolute minimum for the particle model (d = 0) in the limit of large N → ∞ [11]. The main goal of the present paper is to extend this study to the pinned elastic manifold. Hence we will study the N L d × N L d Hessian matrix in particular its density of eigenvalues ρ(λ) normalized as ρ(λ) dλ = 1. An important feature of such matrix is its (block-)band structure visualized below: Structure of the (N L × N L) Hessian matrix K in 1d discrete lattice model with L internal sites and periodic boundary conditions.
Only non-zero N × N blocks are indicated.
where for r = 1, . . . , L we have introduced N ×N random matrices W (r) with entries W (r) Our main focus here is the problem where the Hessian K ix,jy [u 0 ] is chosen at the global minimal energy configuration u 0 ≡ u 0 (x). At the same time it is worth noting another interesting problem, where the Hessian is not conditioned by the global energy minimum, but instead chosen at a generic point in configuration space, i.e. at an arbitrary fixed u(x). It is easy to see from (3) and from the statistical translational invariance of the correlator in (2) that the Hessian is then statistically independent of the choice of u(x), i.e. we may as well chose it at u(x) = 0. The covariance structure of the random potential (2) implies, after a simple differentiation that entries of the matrices W (r) are mean-zero Gaussiandistributed, independent for different r and have the following covariance structure: The matrices of such block-band type, with W (r) in diagonal blocks replaced with GOE matrices with i.i.d. entries, were introduced by Wegner [72] in his famous studies of the Anderson localization, and are now known by the general name of Wigner orbital models. Various instances of the models kept attracting attention in Theoretical and Mathematical Physics literature over the years, see e.g. recent paper [73] and references therein. In particular, the mean eigenvalue density for such type of models as N → ∞ is known to be determined by the deformed semicircle equation rigorously derived in [74]. That equation naturally generalizes the so-called Pastur equation of random matrix theory [75]. We will see below that the difference between GOE covariance and our choice (4) is immaterial for the calculation of the mean eigenvalue density which will be found to satisfy exactly the same deformed semicircle equation. Moreover, when we condition the Hessian by being at the global energy minimum the equation retains its validity, albeit with the renormalized curvature parameter µ → µ ef f , which should be determined by a separate minimization procedure. The replacement µ → µ ef f is crucial in determining the global position of the support of the density of states, i.e. the position of the edge(s) and the value of the gap, for the Hessian at the global energy minimum, but the general form of the density can be already determined without that knowledge by studying the above mentioned equation. An analysis of the density profile following from that equation hence comes naturally in our problem as by-product of its solution. Surprisingly, we were not able to trace it in the literature for the most interesting case of infinite system of size L → ∞. As it may have a separate interest and is quite instructive, we are going to fill in that gap in the present paper and provide such an analysis for d = 1.
In the case of a continuous manifold the Hessian matrix K becomes a matrixvalued differential operator K acting in the space of N −component vectors f (x) := (f 1 (x), . . . , f N (x)) T where, e.g. x ∈ [0, L] d , by the following rule: with appropriate boundary conditions (e.g. periodic, or Dirichlet). Without conditioning by global minimum the covariance structure ofŴ is a natural analogue of (4): In particular, for d = 1 the operator can be visualized in the following form of an N × N matrix §: Models of such type are sometimes called the matrix Anderson models, and are essentially continuous versions of Wegner orbital models. Again, we will show below that the associated 'deformed semicircle' equation for the mean eigenvalue density of such problem can be solved as long as L → ∞ and yields an explicit form of the density profile.

Summary of the main results
In this paper our main object of interest is the disorder-averaged resolvent (Green's function) of the Hessian, calculated at the absolute minimum u 0 of the total energy: as well as its limit at coinciding points G(x, x; λ, u 0 ), which relates to the mean spectral density of the Hessian as Employing the replica trick, we first show that for N → ∞ (the limit being taken for a fixed value of L d ) the average Green's function is given by where the value of the parameter p is determined by the following self-consistent 'deformed semicircle' equation for the diagonal part § Note that in such a case the spectral density in the infinite-volume limit can not be which is essentially of the same form as one for the orbital model with lattice Laplacian [74]. The only quantity which contains all the information about the optimization leading to the ground state u 0 is the parameter µ eff . Below µ eff will be calculated in the various cases (replica-symmetric, 1RSB and FRSB) in the framework of the replica theory. We recall that the notation k applies both to the discrete models k = 1 L d k and the continuum limit k = d d k (2π) d . These equations are quite general and apply to basically arbitrary graph Laplacian matrices t∆ ix,jy (even not translationally invariant ones, provided the formula are are generalized by replacing k 1 A → trA −1 , i.e. the trace in internal space of the inverse matrix).

Spectral Density of the Hessian at a generic point
As has been already mentioned, with setting µ ef f = µ the above expressions (10,11) provides the mean resolvent and the mean spectral density ρ(λ) for the manifold Hessian around a generic point of the disordered landscape. Such object is interesting by itself and we study the shapes of the spectral density in detail for several examples. Its generic feature is the square-root singularity at the spectral edges, which is thus a universal characteristics of the mean-field type spectral densities for disordered elastic systems of any dimension d. The shape as a whole is not universal and essentially depends on the dimension and the type of the Laplacian matrix (discrete or continuous).
As relatively few explicit formulas are available in the literature for eigenvalue densities of disordered matrices and operators beyond Wigner semicircular, Marchenko-Pastur and 1D chains ( see the book [76] for those and further examples) we want to emphasize that in our model it turns out to be possible to find explicitly the spectral density for the 1D matrix Anderson model (7), of infinite length L → ∞ and the Laplacian spectrum −∆(k) = k 2 , −∞ < k < ∞: We have plotted in the Fig. 1 the parameter free scaling function r c (Λ). The spectral edge Λ e is given in this case by Λ e = −1. The function r c (Λ) reaches its maximum at Λ = 0 and then decays at Λ 1 as r c (Λ 1) ∼ 1 √ 3Λ . The latter regime corresponds to the spectral density ρ(λ) = 1 of the disorder-free operator µ − d 2 dx 2 with the spectrum λ = µ + tk 2 .
normalized. E.g. recall the density for the disorder-free model with N = 1 given by ρ(λ) = In the case of 1D disordered elastic discrete chain with −∆(k) = 2(1−cos k), 0 ≤ k ≤ 2π the shape of the spectral density for the associated banded Hessian (and hence for the related Wegner orbital model) can be shown to be of the form but the function r(Λ, y) does not have a simple form for y ∼ 1. However, in the limiting case of weak disorder y 1 a very explicit characterisation is again possible. In this case the graph r(Λ, y) has two spectral edges at Λ (−) e = − 3 2 y −2/3 and Λ (+) e = 2 + 3 2 y −2/3 and the density profile in the vicinity of the edges is simply related to the density profile r c (Λ) of 1D continuous system. Namely, in the vicinity of the left edge r(Λ, y) ≈ y −2/3 r c 2 3 y 2/3 Λ , |Λ| ∼ y −2/3 (15) and essentially the same profile in the vicinity of the upper edge |Λ − 2| ∼ y −2/3 . In between the edges, for any finite 0 < Λ < 2 the profile for y >> 1 is given by the "disorder-free" shape r(Λ, y) 1/(y Λ(2 − Λ)). The numerically calculated spectral density for y = 10 is presented in Fig. 2 and shows all those features.
After this digression about the Hessian spectral densities at a generic point of the disordered landscape, we return to our main task of analysing the Hessian spectra conditioned by the requirement of sampling at the global minimum in the landscape, which requires the determination of µ ef f . Before briefly summarizing our main results, we need to be more specific about the correlations of the landscape, i.e. the choice of B(q) in (2). The corresponding discussion is given below. . Blue: scaling function for the Hessian spectral density, r(Λ, y) versus Λ = λ−µ 2t for the infinite discrete 1D chain given by Eq. (14), for y = t 2 B (0) = 10 (weak disorder). In the weak disorder limit, the central part converges to the spectral density without disorder (indicated here in orange), while the two parts around the edges converge, upon rescaling, to the density for the continuum model plotted in Fig. 1, according to Eq. (15).

Correlations of the random landscape and main features of the phase diagram
For a general classification of the functions B(q) corresponding to allowed covariances of isotropic stationary Gaussian fields we refer to [11] and references therein. Here, for applications to elastic manifolds we mainly consider the powerlaw class when the derivative B (q) can be written as in [15] : As special limiting cases this class also includes the (i) exponential B(z) = B 0 e −z/r 2 f as the limit γ → +∞, and (ii) the log-correlated case for γ → 1. Here r f is the correlation length of the random potential which enters in Larkin's theory, and B 0 has dimension of energy square. For notational simplicity we will consider Let us recall the main features of the replica solution [12,15] for N → ∞ (restricting for simplicity to d 4). Let us first define the "Flory" roughness and the free energy fluctuation exponents We follow the definitions and notations of [15] (section II A and B) which are consistent with the original paper [12]. The parameter γ is thus identical to the one defined in our previous work [11].
Then it was found that for µ < µ c (T ), full replica symmetry breaking, FRSB, occurs whenever θ F (γ) > 0, and 1-step replica symmetry breaking, 1RSB, occurs when θ F (γ) 0. The first case, FRSB, thus always occurs for manifold of dimensions 2 < d < 4, whereas for 0 ≤ d < 2 it is possible whenever 0 < γ < γ c (d) = 2 2−d . In that case the exponents ζ, θ (which are defined in the limit µ → 0) are given by their Flory values. In the limit µ → 0 the system was shown to remain in the glass FRSB phase at any temperature T (no transition). The second case, 1RSB, occurs for d < 2 and γ > γ c (d). In that case there is a phase transition at T c (µ) which survives for µ = 0. It is worth mentioning that in the marginal case γ = γ c (d) this transition is of a continuous nature.
The exponents are θ = 0 and ζ = 2−d 2 in both the high-T phase and the low-T 1RSB phase, with however different amplitudes ¶. The special case γ = γ c (d) is called marginal and exhibits features of both 1RSB and FRSB. Note that it also includes as a special limit the case of d = 2 and the disorder with exponential covariance.
In [11], for the case of a single particle d = 0, we have distinguished longrange correlated (Full RSB) 0 < γ < 1, and short-range correlated (1-RSB) γ > 1 landscapes. For the manifold such as distinction thus also holds, however the critical value of γ is not unity anymore, but equal to γ c (d) = 2 2−d . In particular for d > 2 one is always in the LRC case + . This is because the total energy now also includes the elastic energy, which increases the correlations of the effective random landscape seen by the manifold.

Hessian spectrum at the point of global energy minimum
Our results here extend the ones of [11], which are recovered in the special case of d = 0. There are many similarities with that case. The most important parameter in the theory is the "Larkin mass" µ c > 0 which controls the value of the parabolic confinement µ below which the replica symmetry breaking (RSB) occurs. Its value turns out to be given by the positive solution of which is controlled both by disorder strength and the elasticity matrix. For example, for 1D continuous system a simple calculation gives . Our analysis shows that in the replica symmetric phase the lower spectral edge λ (−) e of the Hessian (which we associate with the spectral gap) as a function of µ is given by ¶ This is an artefact of N = ∞ for 1RSB and may not survive for N = 1, except maybe in the boundary case γ = γ c (d) (certainly it survives for d = 0, γ c = 1 the log-correlated case). + Note however that for d > 4 there is again a RS phase for weak disorder.
This formula immediately shows that for µ > µ c the Hessian spectrum is always gapped (from zero). Upon expanding for µ → µ c and using (19) one immediately finds the gap vanishing quadratically at µ c . For µ < µ c the Hessian spectrum is either gapped or extends down to zero, depending if 1-step RSB or full-RSB occurs. In the first case, the gap vanishes as (µ c − µ) 4 near the transition from below, with the super-universal exponent. For example, for the continuum model in dimension In the second case of full RSB the gap is identically zero everywhere for µ ≤ µ c .
We also obtain the average two point Green function (8) and we find that at the edge of the spectrum it decays exponentially as ∼ e −|x−y|/Lc , with the characteristic length precisely equal in all cases to the Larkin length L c introduced in the theory of pinning. For the continuum model with short-range elasticity and weak disorder, c . This is thus reminiscent of the results of [71] although obtained there in a slightly different context (depinning). Remarkably, this property holds also for µ > µ c , i.e. in the RS phase.
As a by product of these studies we arrived to a very precise criterion which allows to determine which types of covariance functions B(q) in a given manifold dimension d will lead to the full-RSB solution. It reads which generalizes the criterion given in [77] for d = 0. Inserting B(q) for the power law models (16) gives a criterion in agreement with one given in [12], namely that the full RSB solution holds (i) for any value of γ if d 2 and (ii) for γ γ c (d) Finally, for d = 1 the above criterion specifies the covariance B(q) = A c+q as the special marginal case which shares simultaneously the features of 1RSB and FRSB, and for d = 2 the exponential B(q) ∼ e −aq plays a similar role. In particular, the Hessian spectrum is gapless in those potentials. We study both cases in much detail and show that for them the Parisi equations can be solved exactly and explicitly. Note that in N → ∞ class of models these cases play the same role for d = 1 and d = 2 as the logarithmically correlated case identified as marginal in d = 0 [77]. It is worth mentioning here that due to marginality many special properties of logarithmically correlated potential in d = 0 survive for finite N , as was originally suggested in [78] and much studied in the last decade, see e.g. [79][80][81]. It would be interesting to investigate whether some universality holds for the finite-N elastic disordered systems in the above marginally correlated cases for d = 1, 2 as well.
Let us mention here some works on related models, although they are more similar to the case d = 0, and not the manifold. In [41,50] the Hessian statistics is sampled over all saddle-points or minima at a given value of the potential H(u) = E = const, a priori quite different from imposing the absolute minimum. The spectrum of the soft modes was also calculated in a mean-field model of the jamming transition, the 'soft spherical perceptron'. The Hessian matrix in that model has the shape of a (uniformly shifted) Wishart matrix, whose spectrum is given by the (shifted) Marchenko-Pastur law, while in [11] the Hessian spectrum is given by a shifted Wigner semicircle. The model has two phases: ' RS simple' and 'FRSB complex' and the Marchenko-Pastur spectrum in that model was demonstrated to undergo a transition from gapped to gapless, similar to what we find here for Gaussian landscapes. Finally, it is worth mentioning a quite detailed recent characterization of the energy landscape of spherical p−spinglass in full-RSB phase close to the global minimum, see [82] and references therein.
The outline of this paper is as follows. In Section 3 we provide a derivation of the average Green function, resolvent and the spectral density of the Hessian using two sets of replica. The second set is necessary to specify that the Hessian is considered at the absolute energy minimum. We obtain the general saddle point equations which determine these quantities. In Section 4 we analyze the results. In the first subsection 4.1 we obtain the spectral density and the Green function keeping µ eff as a free parameter. The general results only weakly depend on this parameter, which simply globally shifts the support of the spectral density. In the second part 4.2 we complete the study by calculating µ eff from the explicit solution of the replica saddle point equations. This leads to the determination of the spectral edges and of the gap, in the three main distinct cases: replica symmetric, FRSB and 1RSB. The case of marginal 1RSB is given a special attention. Finally Section 5 contains the conclusion.

Derivation of the average Green function using replica
Below we use the following notational conventions. The sums over the internal points of the manifold x, y, .. are denoted as x ≡ L d x=1 , the sum over the first set of replica indices α, γ, .. are denoted α ≡ m α=1 , the sums over the second set of replica indices a, b, c.. are denoted a ≡ m a=1 , and similarly for the products. The indices i = 1, ..N and the dot product is used in R N .
The notation Tr is the trace over all indices x, i and a or α, i.e. over R m × L d or R n × L d , e.g. Tr A = xa A xa,xa . The notation tr is reserved for the traces over a or α only, i.e. over R m or R n , i.e. tr A = a A aa .

Green's function and the first set of replica
As the starting point of our approach, we introduce the resolvent of the Hessian K(u) defined in (3), for a given generic configuration u(x) (not necessarily the minimum of the total energy) and in a given realization of the random potential V (u(x), x).

The associated Green's function is then defined via
Such Green's function admits then the following representation in terms of m replicated Gaussian integrals over N −component real-valued vector fields φ α (x), with α = 1, . . . , m: where we assumed that Im λ < 0 and set the factor ( i π ) m/2 → 1 for m = 0. From this we calculate the mean spectral density of the Hessian eigenvalues "at a temperature T ", defined as where the thermal averaged value of any functional g(u) of a configuration u( being the Boltzmann-Gibbs weights associated with the configurations via the energy functional (1). Our final aim is then to obtain the mean spectral density of Hessian eigenvalues at the absolute minimum by setting temperature to zero: The problem therefore amounts to first calculating the disorder and thermal average where e i 2 x,y,α φα(x)·K(u)·φα(y) and then, by performing the zero-temperature limit, to capture the contribution from the global minimum configuration only.

Average Green function and second set of replica
In the framework of the replica trick we represent the normalization factor Z −1 β in Eq.(28) formally as 1/Z β = lim n→0 Z n−1 β and treat the parameter n before the limit as a positive integer. After this is done, averaging the product of n integrals over the Gaussian potential V (u) is an easy task. The calculation is very similar to the one for d = 0 in [11], (apart from an additional factor of 2 for each derivative of B arising due to a slightly different normalization of the covariance used in [11] ) and we simply quote the result referring to [11] for more details. We obtain e i 2 x,y,α φα(x)·K(u)·φα(y) with where the u-independent part of the action is given by We can thus rewrite Now we introduce auxiliary fields and their conjugate fields. We define for each value of the argument x (which we omit for brevity) the differentials and then define DQ(x) = x dQ(x), etc. This allows us to use the formal identity where the contours of integration are duly chosen. This identity can be then inserted inside (33) effectively allowing to replace all scalar products and u a (x) · φ α (x)) by Q ab , P αβ and R a,α respectively inside L n,m [u, φ], leaving a simple quadratic form in the fields u a (x) and φ α (x), which can be integrated out.
Restricting for now, for simplicity, to the diagonal element of the resolvent we obtain where we have defined the action The last piece δL[Q, σ, P, τ, R, η] is given in Appendix A. Since it vanishes at the saddle point we do not need to give it here.
We can now write the saddle point equations. It is easy to check that the equations admit an x-invariant solution in all variables: i.e. Q ab (x) = Q ab , σ ab (x) = σ ab , P αβ (x) = P αβ , τ ab (x) = τ ab at the saddle point. We will consider only this solution on physical grounds. Let us define again the quantity Taking first the functional derivatives w.r.t. τ αβ (x) and P αβ (x) we arrive at the saddle-point equations Moreover, similarly to d = 0 case treated in detail in [11], it is easy to check that the invariance of the action under rotating matrices P and τ in the replica space implies that the corresponding saddle point solutions must be actually proportional to the identity matrix: P αβ = pδ αβ and τ αβ = τ δ αβ . In the limit m → 0 one then finds that τ satisfies the equation which when substituted to the corresponding equation for p yields the closed selfconsistency equation for the latter: This condition has exactly the form of the 'deformed semicircle' equation derived in [74] for the block-banded Wegner orbital model, assuming the random matrices W (r) , r = 1, . . . , L on the main diagonal to be of standard GOE type. To that end it is worth noting that in the action eq.(39) eventually responsible for fixing the shape of the self-consistency equation the difference between our choice for W (r) , see (4), and the GOE appears only via the term (tr P ) 2 absent for GOE case. However for m → 0 that term gives contribution of the order m 2 , hence is negligible in comparison with the dominant contributions of order m. Hence, from the point of view of calculating the profile of the mean eigenvalue density the difference between our Hessians and the Wegner orbital model is immaterial. In particular, in the case of d = 0 one recovers the self-consistency condition found in [11] ip = 1 λ − µ eff − 4ipB (0) (46) whose solution is the genuine semicircular density as typical for d = 0 random matrix problems. The analysis of the solution to the self-consistency equation (45) for higher d, and especially for d = 1 will be provided in detail below. Along similar lines one can derive the average Green's function at two different points x = y. Starting from its definition and using a source term, it is not difficult to see that in the limit of large N employing the same saddle point solutions one arrives at the following representation: where p is an apriori complex number determined by the (self-consistent) equation for the diagonal part. Finally, using this saddle point, we see that the term which couples P and Q is proportional to m and hence as m → 0 can be neglected in the saddle point equation for Q. The resulting equations are identical to those obtained in [12,14,15] and we now briefly recall them here. Taking the functional derivatives w.r.t. σ ab (x) and τ αβ (x) yields where we define, as in [12,14,15] We will recall briefly the analysis of these equations below, as needed.

Analysis of the results
We now analyze the results from these saddle-point equations in two stages. First we analyse the general form of the average Green function and the spectral density by simply assuming that µ eff takes some value at T = 0. That gives the shape of the spectral density ρ(λ), up to a global shift of λ. In a second part, we recall the analysis leading to the various phases (RS, FRSB and 1RSB) and obtain from it the corresponding possible values of µ eff as a function of µ, which allows to determine the location of the edge of the spectrum.
There are usually multiple solutions for p and we must choose the branch such that for λ → ±∞ one has ip ∼ 1 λ . For a discrete model the spectrum of the perturbed Laplacian is bounded and large |λ| necessarily correspond to being outside of the spectrum. In the continuum model the same holds for large negative λ. In the range of λ outside of the spectrum, p is necessarily pure imaginary. When λ reaches the edges of the spectrum and goes inside the spectral support, p develops a real part proportional to the mean spectral density. Hence we can write for real p 1 , p 2 which converts (52) after separating the real and imaginary parts, into two coupled equations which determine p 1 and p 2 as functions of λ: The edge of the spectrum is for λ = λ e such that p 1 acquire a non zero value, hence it is determined by eliminating p 2 = p e 2 in the system since at the edge one can set p 1 = p e 1 = 0. Note that there can be more than one edge, i.e. more than one solution to this system. Note also that assuming the right-hand-sides in (54) and (55) are analytic functions F, G of all arguments λ, p 2 and u = p 2 1 , a straightforward expansion in powers of λ − λ e shows that just above the lower edge implying a square-root singularity of the density of eigenvalues at the thresholdin the generic case (where neither the numerator or denominator vanishes in (58)).
To further analyze these equations we introduce the Larkin mass µ c > 0 defined as the positive solution of Anticipating a little on the subsequent analysis, (59) precisely determines the range of curvatures when the replica-symmetric solution becomes unstable. Namely, it becomes unstable in the interval 0 ≤ µ < µ c , with µ c determined by (59). The Larkin mass exists whenever B (0) > 1/(4 k 1 −t∆(k) ) and when this is the case, it is unique. In the opposite case, the RS solution is stable for all values of µ. Note that µ c depends only on B (0) and on the graph Laplacian elasticity matrix.
We now assume that we are in the first case and there exists a finite Larkin mass µ c > 0. It is then easy to find a solution for the spectral edge λ e . One sees that p e 2 is now determined in terms of µ c and that (56) is equivalent to which determines λ e as a function of µ c . It turns out (see below) that this is always the lower edge, hence we denoted it λ − e . We discuss below how to obtain the other edge(s) when they exist. (i) First recall that for a single-site (equivalently, zero dimensional d = 0) system with L d = 1 the equation (52) gives Hence the Hessian spectral density from (53) is given by the semicircular law and zero otherwise. This is precisely the result obtained in [11]. On the other hand we can determine the edge using (56)- (57). First let us examine the equation (59). In that case it reads The positive root is µ c = 2 B (0) and from (60) we find and recover the lower threshold (i). If we now use the negative root of (63) and diverges at the edges of this interval. Hence one expects two roots, one for µ c = µ + c > 0, and, by symmetry, one for µ c = µ − c = −4t − µ + c . In the following we will always associate the positive root µ + c = µ c with the Larkin mass. If the set −D consists of several intervals, or several points, where the r.h.s. of (59) is infinite, there can be several additional solution to (59) besides one associated with the Larkin mass. That one we know must be the largest one since −t∆ is required to be positive definite. Hence it corresponds to the lower edge of the Hessian.
(ii) As soon as we have L 2 the spectral density for the Hessian is not a semicircle as we now discuss. For a line d = 1 with L points and periodic boundary conditions the eigenvalues of −t∆ are 2t(1 − cos 2πj L ), j = 1, . . . , L − 1. The equation (59) becomes It is easy to see that for very weak disorder, B (0) 1, there are 2L roots to this equation, which we denote µ c = µ j,± c , j = 0, . . . , L − 1 with which can be found by considering successively all the quadratic divergences in each term in the sum in the left-hand side of (65) and approximating the sum accordingly. The formula (60) for the corresponding edge becomes Substituting here the values of µ j,± c found above we arrive, up to subdominant terms at weak disorder, to the corresponding edge values For L = 1 the above approximation is exact and one recovers the formula (i) valid for any disorder. For L 2 there are 2L edges and L bands at weak disorder. It is easy to see why. When disorder is zero the Hessian is simply the Hessian of the elastic matrix and its spectrum is the set of delta peaks at 2t(1 − cos 2πj L ) + µ (in that case µ eff = µ). As disorder increases, each of these delta peaks broadens, leading to a band, as described by (68). One can expect that these bands will remain well separated as long as their width 8 B (0) L is much smaller than their separations ≈ 4t L . This gives the criterion to have separated bands. It is reasonable to expect that in that situation each band will have a semi-circle form, since each basically solves independently the d = 0 equation.
To study the merging of such bands, let us consider the case L = 2 in more details (the eigenmodes are then k = 0, π). The equation becomes Hence there are two cases. Either disorder is weak 2t 2 B (0) > 2, and there are 4 real roots with |z ±,+ | > 1 and |z ±,− | < 1 always. Or disorder is strong and only the two roots z ±,+ exist. These roots correspond to edges of the spectrum of the Hessian There are thus 4 edges (weak disorder) and 2 edges (strong disorder). The lowest edge is located at We can now study the spectral density for the case L = 2. It is given by ρ(λ) = 1 π Im(ip) where ip satisfies the cubic equation The resulting density of states is plotted in Fig. 3 and 4. The evolution described above from disjoint supports (weak disorder) to a single support (strong disorder), as well as the transition at B (0)/t 2 = 1 is clearly visible.
(iii) Our next example is the continuum 1D line of infinite length L → +∞. Since the Laplacian spectrum for such a system is given by −t∆(k) = tk 2 with k ∈ [0, +∞[ for such a case there is only one spectral edge in the system with disorder, the lower edge λ − e , determined from the Larkin mass µ c , i.e. the unique positive solution of Hessian at the minimum for manifolds in a high-dimensional random landscape 23   leading to The spectral density in the interval λ > λ − e is then given by ρ(λ) = 1 π Im(ip) where the complex p is obtained by solving the following equation: Introducing new scaled variables y,p,λ via and rescaling the integration variable as k → k y −1/3 the equation (78) takes the form Hessian at the minimum for manifolds in a high-dimensional random landscape 24 Taking the square and further introducing the variable w and parameter δ bỹ the equation for w attains especially simple form: The general theory of cubic equations then dictates that for δ < −1 the equation (82) has only real solutions, hence p ∼ −iw −1 will be purely imaginary implying zero density of eigenvalues. This parameter range fully agrees with the position of lower spectral threshold λ − e found in (77), which in the new variables reads λ − e = −3y −1/6 . As long as δ > −1 there is one real root given by the Cardano formula in the form which is positive and decreases from 2 to 0 as δ increases from −1 to +∞. There are also two complex conjugated solutions w and w. To find their imaginary part we use the Vieta's formulas: which gives Re(w) = − wr 2 and Im(w) = ± 2 wr − w 2 r 4 . It is then easy to see that the spectral density can be found explicitly and is given by: Im(w) Re 2 (w) + Im 2 (w) where we have to choose the sign of Im(w) which ensures positivity of the mean density. Recalling that in terms of the original variables We have plotted in the Fig. 1 the parameter free scaling function r c (Λ) = √ U − U 4 , U = wr 2 with w r given by (83) and δ = Λ 3 . It has the following asymptotics for large Λ and for Λ near the edge Λ e = −1: In particular, equation (84) then implies that showing the expected square-root singularity close to the spectral edge λ = λ − e . Moreover, it is easy to see that drc dU = 0 for U = 2 −2/3 , hence w r = 2 1/3 corresponding according to (83) to δ = 0. We conclude that Λ = 0 is exactly the point of the maximum for the scaled density profile r c (Λ), which is readily seen from Fig. 1. (iv) Our last example is an infinite discrete lattice d = 1 with the number of sites L → ∞. The Laplacian spectrum for such a system is given by To determine the spectral edges we use the integrals 2π 0 dk 2π for real |x| > 1. This reduces finding the roots of (59) to solving the equation for µ c > 0 or µ c < −4t. Denoting r = 1 + µc 2t and y = t 2 B (0) we rewrite the above equation as our equation has only a single real root w = w c given by the Cardano formula: which is obviously positive as needed for our goals. For the parameter µ c we then have two solutions. The positive one corresponds to the Larkin length: and the second solution µ − c = −2t 1 + w c y −2/3 + 1 ≡ −4t − µ c as expected by symmetry. In the case 0 ≤ y ≤ 2 3 √ 3 the cubic equations has all three real roots. Introducing the angle θ ∈ [0, π 2 ] such that cos θ = 3 √ 3 2 y the roots can be conveniently written in the so-called trigonometric form: c can be used for the above procedure and yields µ ± c . Finally, this gives us the two spectral edges as Let us give a simple example: t = 2 1/2 , B (0) = 3 3/2 which gives y = 2 3 3/2 , hence ∆ = 1/2 and w c = 2 2/3 which eventually gives for the Larkin length µ + c = 2 To calculate the spectral density profile one needs the following generalization of (88): valid for any real t > 0 and complex a such that |a + √ a 2 − 4t 2 | − 2t = 0, excluding a real in the interval [−2t, 2t]. The spectral density in the interval λ − e < λ < λ + e is then given by ρ(λ) = 1 π Im(ip) where the complex p is obtained by solving the following equation, cf. : where we denoted a(p) = λ−2t−µ eff −4ipB (0). We define the scaled variables which brings Eq. (93) in the dimensionless form . The roots of (96) must satisfy the following equation forP = iP The spectral density is then given in terms of the function r(Λ, y) = Im(P ) as The parameter y reflects the strength of the disorder relative to the elasticity, so that the larger is y the weaker is the disorder. For strong disorder (or vanishing elasticity), y → 0, the system decouples in non-interacting zerodimensional units and the spectral density is given by the semicircular law, as can be seen e.g. setting t = 0 in (93). For a moderate disorder y ∼ 1 the shape is not a semicircle any longer, but is qualitatively similar, as can be seen in Fig. 5 where the scaling function r(Λ, y) is plotted for y = 1. However with decreasing disorder/increasing elasticity the shape of the spectral density changes qualitatively and develops a characteristic form with two maxima and a minimum in between, see plot for y = 10 in Fig.2. Some hints towards the origin of such shape can be obtained by considering the limit of vanishing disorder B (0) → 0, i.e. y → +∞. In this limit one expects that the spectral density should in a certain sense converge to the one of the purely elastic 1d system: which implies that for y 1, r(Λ, y) 1/(y Λ(2 − Λ)). The correspondence with the disorder-free result is visible on the Fig. 2. in the central part around the minimum. To understand the two-maxima shape we investigate analytically the case of large but finite y 1 more accurately. RescalingP = q/y the equation (97) takes the form: Now it is obvious that letting y → ∞ for a fixed 0 < Λ < 2 the above is reduced to the quadratic equation with purely imaginary roots . This solution yields precisely the density for the pure elastic case (99). However, it is also evident that in the vicinity of the points Λ = 0 or Λ = 2 such naive limit breaks down and requires a separate treatment. We illustrate it by providing analysis in the vicinity of Λ = 0, one for the region around Λ = 2 being fully analogous. A simple scaling argument demonstrates that the relevant vicinity of Λ = 0 is of the width |Λ| ∼ y −2/3 so that it makes sense to introduce a new parameter δ = 2 3 y 2/3 Λ 3 and also introduce the scaled variable w via q = y 1/3 w . Substituting this into (100) and taking the limit y → ∞ one finds the equation for w exactly given by the equation (82) studied by us in much detail above in our analysis of the 1d disordered continuum problem. We therefore conclude that for y → ∞ and around Λ = 0 the scaled spectral density profile r(Λ, y) for the 1D discrete model is simply given in terms of the scaled density profile of the continuum model r c (Λ) obtained in (84) as r(Λ, y) = y −2/3 r c 2 3 y 2/3 Λ . In particular, recalling that r c (Λ = −1) = 0 in the continuum case, we find that the position of the left spectral threshold in the discrete case for y >> 1 is given by Λ − e = − 3 2 y −2/3 . Close to this threshold the density increases as the square root Λ − Λ − e , eventually reaching its maximal value r(0, y) = y −2/3 2 −4/3 √ 3 exactly at Λ = 0 and then decaying for larger Λ >> y −2/3 in agreement with the asymptotics (85) as r(δ 1) ∼ y −2/3 1 which precisely matches the Λ << 1 behaviour from the 'central part' r(Λ, y) 1/(y Λ(2 − Λ)). This demonstrates that for weak disorder, y 1, the shape of the spectral density for the 1D infinite elastic lattice is given by (i) a central part which converges to the pure, "disorder-free", density of states (ii) two edge regions, |Λ| ∼ y −2/3 and |2 − Λ| ∼ y −2/3 , where the divergent density of states of the pure system is converted into a finite profile, identical upon rescaling to the one of the 1D disordered continuous elastic line.

Phases from replica, determination of µ eff and of the gap
In the previous section we have obtained the spectral density and its support, in particular the lower edge, for various cases. The formulas however contained a single as yet unknown parameter µ eff which corresponds to a global shift of the support of the Hessian spectral density.
In this section our aim is to determine µ eff , the missing information about the global position of the Hessian spectrum. As µ is varied the system can be in different phases (RS, 1RSB, FRSB) and the formula leading to µ eff must be detemined accordingly. In each case we first recall briefly the known replica saddle point solutions for the random manifold problem [12,15], i.e. the solutions to the equations (48)-(50).
4.2.1. Replica-symmetric phase Let us start with the replica symmetric (RS) phase, which occurs for µ > µ c . Let us look for a replica symmetric solution of the saddle point equations (48)- (50).
Note that the condition b σ ab = 0 in this parametrization reads σ c + nσ = 0 which in the replica limit n → 0 implies that we can choose σ c = 0. We also have n) and the inversion of the RS matrix gives: or equivalently in the replica limit n → 0 with χ aa = 0 by definition. The saddle point equation (50) leads then to the explicit formula for σ, which determines completely the solution and Q ab = T k G ab (k). As is well known [12] the RS solution is valid for µ > µ c (T ) with (see Eq. (17) in [15]) which gives, in the T = 0 limit, µ c (T = 0) = µ c > 0, i.e. the Larkin mass determined by as anticipated in the previous Section.
We can now determine µ eff and the edges of the Hessian in the RS phase. From (41) we obtain (for n = 0) Substituting this value of µ eff in (60) we thus obtain the final formula for the lower spectral edge λ (−) e of the Hessian (which we associate with the spectral gap) as a function of µ in the RS phase This formula immediately shows that the gap vanishes quadratically at µ c , i.e. upon expanding for µ > µ c where the linear term cancels as a consequence of (106). In d = 0 we recover the similar formula obtained in [11]. For the continuum model −t∆(k) = tk 2 there is only one edge, the lower edge, which we just determined. However, as extensively discussed in the previous Section, for other models (e.g. discrete models) there may be several edges (and bands). As discussed there extensively all the edges λ α e are obtained by considering all the real roots µ α c of (106) and inserting them in the formula (108). We refer to that Section for details. Let us simply note the case of the discrete d = 1 model, where by symmetry the two roots of (106) are µ c and −4t − µ c . That gives the upper edge for that model, and one can then write a single formula for both edges λ (±),1d,discrete (110) which gives simple formula for the midpoint and the band width.

Full RSB phase
Let us now discuss the full RSB solution. We choose not to give the finite-n hierarchical structure here as it is relatively cumbersome, but rather simply follow the n → 0 analysis * in [12][13][14][15]. The off-diagonal part σ a =b is represented by a Parisi function where from RSB replica matrix inversion one has 2 (112) * in particular the Section VIII B of [14] (note the small misprint in the definition of σ at the very Taking in (111) the derivatives w.r.t. v and exploiting (112), one finds that for any interval of v either (i) σ(v) is constant or (ii) it satisfies the marginality condition In particular at the breakpoint v = v c one has which by comparison with RS stability condition (105) implies that and at T = 0, as a function of µ c the Larkin mass determined by (106), one has For use here and in the next Section let us define the notations for l 1 and somewhat abusively One can calculate the solution for [σ](v) for arbitrary covariance B (see e.g. formula (8.16) in [14]) as we now show. We assume that B is a monotonous decreasing function (B < 0). Inverting the marginality condition (114) and inserting in (111) leads to Taking a derivative of the above w.r.t. v with the help of the identities where f −1 is the functional inverse of f , we obtain after rearranging and assuming beginning of the section there, paragraph above (8.4): it should be T G −1 ab (k)−(k 2 +m 2 )δ ab = −σ ab ). Note that f =f = −B in [12] and we follow [12] for the definition of G(k) i.e. we do not absorb which determines by inversion [σ](v) as a function of v, which must be an increasing function. It is now convenient to introduce and define the function F (b) via the relation As the result, the relation (121) can be rewritten as Taking yet another derivative w.r.t. v and noticing that dA dv < 0, this leads to the following condition for FRSB solution to exist d dq where we used that dq/dv < 0. More precisely the condition for FRSB to hold in an interval [v µ , v c ] is that (123) holds for q ∈ [q(v c ), q(v µ )].
One may now notice that for a d−dimensional continuum model with Laplacian spectrum −t∆(k) = t k 2 and dk ∼ |k| d−1 d|k| the behaviour of the integrals I l>d/2 (x) in (118) for x → 0 is dominated by the infrared (|k| → 0) limit and is given by . Taking d < 4 we then see that (122) implies in the limit of small µ and small b the behaviour The same behavior also holds for discretized models on an infinite d dimensional lattice. Replacing F (B (q)) in (123) with the small-argument asymptotics (124) leads then to the full-RSB condition, which we gave in the Introduction, see (22).
Let us now calculate µ eff for the FRSB solution. For this we first set v = v c in (112) and (113) and get the relations Now inserting the FRSB form into the definition (41) we get which can be further rewritten using (125) and the definition of σ(v) in (111) as: In the limit T → 0, and recalling that Σ c = µ c − µ we find This has the same form as the RS formula (107) where one replaces µ by the Larkin mass µ c , i.e. it can be interpreted as the mass µ freezing at µ c , that is retaining for µ < µ c its critical value.
Let us now determine the lower edge of the Hessian. From (60) we obtain upon inserting (130) Hence the lower edge of the Hessian remains frozen at zero within the FRSB phase for all values of µ. For models with more than one edge, their positions can be found from the other roots of the equation (59) as discussed in the previous Section. One should then insert them in (60) and use (130) without inserting them in (130), the latter being defined in terms of the Larkin length µ c .

1-step replica symmetry breaking phase
We now study SRC potentials which exhibit the 1RSB solution. For the continuum models, the 1RSB solution holds for d 2 and γ γ c (d) = 2/(2 − d).
Let us give a brief account of the 1RSB parametrization and the ensuing procedure. We start with introducing two parameters σ 1 and Σ c in terms of which we construct a v c × v c matrix σ d with entries (σ d ) ab = −Σ c δ ab + σ 1 . The full n × n matrix σ has n/v c identical diagonal blocks σ d , all entries being equal to the value σ 0 outside those blocks. The constraint Σ b σ ab = 0 then yields the relation in the n → 0 limit: Inversion of the matrix G −1 = µ − t∆(k) + σ produces n × n matrix G with the diagonal v c × v c blocks G d having entries (G d ) ab = (G − G 1 δ ab + G 1 and outside those blocks G has identical entries G 0 . The entries in the limit n → 0 remembering (132) are given by relations: which according to (51) leads to To determine the equilibrium values of the parameters involved we rely upon the expression for the free energy Φ(T ) associated with the model given Taking a derivative of the free energy w.r.t. Σ c leads to Let us consider the T = 0 limit. Denoting v c = vT , Q : in terms of the integrals defined in (118) and noticing that in this limit χ 0 → Q and B(χ 1 )/2T → B (0)I 1 (µ + Σ c ) one gets where we have defined with F µ (0) = F µ (0) = 0 and F µ (0) = I 2 (µ). Upon derivation of the zero-temperature free energy w.r.t. Σ c (cancelling the the T inside it as done in [14,15]. see Eq. (170) in Section III.D p17 or the arXiv version. in the paper [15] (up to an irrelevant constant). common factor I 2 ) and v one obtains the following system of equations: which should be augmented with the definition of Q in (138). For small Q > 0 we have Σ c 2vB (0)Q and substituting in (138) we then find that the transition to the phase with nonzero value of Q occurs at µ = µ c determined by which identifies with the definition of the Larkin mass, cf. (59).
We can now give the formula for µ eff in the 1RSB phase. From (41) we have, inserting the one-step RSB ansatz which in the limit T → 0 yields Recalling from (60) that we finally obtain, within the 1RSB phase providing the expression for the position of the lower spectral edge in the 1RSB phase. We now expand below and near the transition: we insert Σ c from the first equation in the second and third, which gives two coupled equations for Q and v. In these equations we insert, with δ > 0, and solve order by order. It is convenient in the calculation to use that I l (x) = (−1) l−1 I (l) 0 (x)/(l − 1)! for l 1. We give only the lowest order recalling that B (0) < 0. To this order one finds that the edge λ e vanishes to order O(δ 2 ). To find the first non-vanishing order, O(δ 4 ) one needs to calculate v 1 , Q 2 , v 2 , Q 3 iteratively. Performing the calculation using the Mathematica software we finally find, after some rearrangments using (143), up to O (δ 5 ) terms (150) This result is very general, for any discrete or continuum model. For the d = 0 'single particle' model, I l (µ c ) = µ l c for l 1 and (150) reduces exactly to the formula (76) obtained in our previous work [11].
We can now specify to the continuum model by setting −t∆(k) = k 2 (we set t = 1 for simplicity) and recall that k denotes Restricting our consideration to d < 4 one can see that F and G are defined by a UV convergent integral, so that we can set the UV cutoff k max to infinity. We can check that the ratio entering in (150) are hence we obtain for the continuum model in dimension d Let us study d = 1, 2 in more details. In d = 2 for the continuum model we have and the transition occurs at From the last equation in (141) we find that Σ c = µ(e 2πQv − 1), and substituting into the other two equations we obtain the system The second equation allows to obtain v = v(Q) and reporting in the first one it leads to an equation for Q.
Let us now specify to the exponential case B(q) = e −cq which we expect to be marginal at the boundary with full RSB. One finds that the solution is remarkably Inserting into (147) we find, for the exponential case, that the edge of the spectrum of the Hessian is exactly at zero λ e = 0 (160) for all µ µ c . This is the confirmation of the case being marginal for d = 2, i.e. it can be obtained as a limiting case from the FRSB side. It is interesting to note that its exact solution is also very simple.
Let us now consider the marginality for d = 1 when Then We must now solve the equations It is convenient to introduce the following variables and parameters: in terms of which the above system takes the form Substituting the last of those equations to the second one and remembering that for Σ c > 0 we have y > √ µ we see that the second equation takes the form y = Ω µ x 2 implying further that z = 1 √ µ − µ Ωx 2 . Substituting these relations into the first equation we see that it can be brought to the form The first solution x 2 = µ 3/2 Ω is however not admissible since it corresponds to y = √ µ and therefore to Σ c = 0. The only nontrivial solution as µ is decreased below µ c is provided then by the remaining root x = µ 1/2 Ω 1/3 and in the original variables finally yields the relations Substituting in (147) we again find λ e = 0, confirming marginality for this case.

Spatial structure of the Green function, pinning and localization
One of the interest in the manifold problem compared to the point (d = 0) is the rich internal space structure. The hierarchical construction of the Gibbs measure encoded in the RSB solution was discussed in the context of the manifold in the Appendix of [12] (see also discussion in [2,3]). In that picture the Gibbs measure is a superposition of Gaussians, with power law distribution of weights, each centered around distinct seed configurations u α (x), with fluctuations controled by an "effective mass" (each Fourier mode has its own decomposition into states).
The picture is either one-step (1RSB) or is hierarchically repeated (FRSB). It was shown that the closeby states (at v = v c ) correspond typically to the scale of the Larkin length (with effective mass µ + Σ c = µ c ), while the large scale statistics (e.g. of u(x) − u(0) at large x) is controled, throughout the glass phase, by the small [σ](v) bevahior (with effective mass µ + [σ](v)). Hence in the FRSB phase it is the small v behavior, corresponding to distant states, which to the non-trivial roughness exponent. Here we study the Hessian at the global minimum and its "soft modes" contain information the structure of the states (we saw in particular that the gap is zero from the marginality condition). We can thus now ask about the spatial structure contained in the averaged Green function. Let us examine again the formula (47). If we choose λ = λ − e , i.e. at the lower edge we can write where we used (60). Hence we see that the averaged Green function decays exponentially ∼ e −|x−y|/Lc , with the characteristic length given by the Larkin length L c . For the continuum model with short-range elasticity and weak disorder, c . This is very reminiscent of the result found numerically in [71] in the context of depinning. Remarkably however, here this property holds also for µ > µ c , i.e. in the RS phase.
Note that in the standard interpretation of the localization theory the decay rate of the disorder-averaged Green's function in the bulk of spectrum defines the so-called mean-free path and generically has little to do with the true localization length. The situation at the spectral edge may however be different as, in contrast to the bulk, in that region of spectrum the Green function is not expected to show fast oscillations with random phase in every disorder realization, whose averaging gives rise to the decaying mean. We therefore expect that the decay rate at the edge may have relation to the localization properties of the lowest eigenmode of the Hessian.

Conclusion
In this paper we have extended our previous work on the spectrum of the Hessian matrix at the global minimum of a high dimensional random potential, to the case of many points coupled by an elastic matrix. This is of interest in several contexts, in particular for disordered elastic systems pinned in a random environment. We have calculated the averaged Green function and its imaginary part, the spectral density, of the Hessian matrix. Technically this was achieved using a saddle point method and two sets of replica, one to express the Green function, the second to impose the constraint of global minimum. The latter requires a replica symmetry breaking solution for the saddle point equations, either of the 1 step kind (1RSB) or with full replica symmetry breaking (FRSB). We have derived the criterion according to which one has the former or the latter, which generalizes the concept of short range (leading to 1RSB) or long range (FRSB) disorder to the case of the elastic manifold.
The main difference with the case of the particle d = 0 in a random potential is that the spectral density of the Hessian is not a semi-circle anymore. We have calculated its form in a number of examples and obtain the values of the edges. We have shown how it can evolve from a many band to a single band as the disorder is increased. In all generic cases however it retains a semi-circle shape near its edges. Especially complete and explicit characterisation of the arising spectral density has been achieved in the 1D continuous system of infinite length.
Concerning the position of the lower edge, we have shown that qualitatively the scenario found for the particle remains valid for the manifolds. For short range disorder cases and µ > µ c the Hessian spectrum is gapped away from zero. At µ = µ c the gap vanishes, i.e. the lowest eigenvalue is zero. For µ < µ c the saddle-point solution is 1RSB and we find that the gap is non zero and vanishes as ∼ (µ c − µ) 4 near the transition. For long range disorder cases we find that the gap vanishes identically for µ ≤ µ c , reflecting the marginality of the FRSB solution. We also identified and studied the cases of marginally correlated disorder in d = 1 and d = 2 which can be of separate interest.
A new feature which emerges in the study of the manifold is the information about the internal spatial dependence of the averaged Green function. We found that near the edge it decays over a length scale identical to the so-called Larkin length, related to µ c , which plays a central role in the theory of pinning. Below the Larkin scale the system responds elastically, while above the Larkin scale, metastability sets in leading to glassy non linear response. Our result in the high embedding dimension limit, are reminiscent to what was found in numerical simulations for elastic strings at the depinning transition, where the localization length of the low lying modes of the Hessian was found to be equal to the Larkin length.
Many questions remain. One is to understand the statistics of the lowest eigenvalues. Clearly it cannot be of the Tracy Widom type since it is bounded by zero. The question of its universality remains open. One possible way to tackle this difficult problem is to study the large deviations for the minimal eigenvalue. Another interesting problem is to generalize counting analysis of the minima and saddle-points from the particle case d = 0 [35,38,39] to the present manifold model with d ≥ 1. Work is presently in progress in those directions.
Finally, the most interesting but very challenging problem is to study the present model by taking limits N → ∞ and L → ∞ in a coordinated way, and scaling the coupling t accordingly to enter the regime when Anderson localization effects in Hessian spectrum should be dominant. It remains to be seen if fieldtheoretical/supersymmetric methods which proved to be instrumental in getting insights into spectra and eigenvectors of matrices of banded type [83,84] could be used successfully in the present problem.