# Manifolds Pinned by a High-Dimensional Random Landscape: Hessian at the Global Energy Minimum

Yan Fyodorov, Pierre Le Doussal

► To cite this version:

Yan Fyodorov, Pierre Le Doussal. Manifolds Pinned by a High-Dimensional Random Landscape: Hessian at the Global Energy Minimum. Journal of Statistical Physics, 2020, 179 (1), pp.176-215. 10.1007/s10955-020-02522-2. hal-02552325

HAL Id: hal-02552325

https://hal.sorbonne-universite.fr/hal-02552325

Submitted on 23 Apr 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Manifolds Pinned by a High-Dimensional Random Landscape: Hessian at the Global Energy Minimum

Yan V. Fyodorov<fnref n="1" /> Ⓟ · Pierre Le Doussal<fnref n="2" />

Received: 19 October 2019 / Accepted: 6 March 2020 / Published online: 19 March 2020
© The Author(s) 2020

## Abstract

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 (Fyodorov et al. in Ann Phys 397:1–64, 2018). 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.

**Keywords** First keyword · Second keyword · More

Communicated by Mehran Kardar.

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

<footnote>Yan V. Fyodorov
yan.fyodorov@kcl.ac.uk</footnote>

<footnote n="1">¹ Department of Mathematics, King’s College London, London WC2R 2LS, UK</footnote>

<footnote n="2">² Laboratoire de Physique de l’Ecole Normale Supérieure, PSL University, CNRS, Sorbonne Universités, 24 rue Lhomond, 75231 Paris, France</footnote>
# 1 Introduction

## 1.1 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–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–9]. Their energy landscape is complex leading to glassy behavior.

The manifold is usually parameterized by a $N$-component real displacement field $\mathbf{u}(x) \in \mathbb{R}^N$, where $x$ belongs to an internal space $x \in \Omega$. $\Omega$ can be either a finite collection of points, such as a subset $L^d$ of an internal space of dimension $d$, $\Omega \subset \mathbb{Z}^d$, for discrete models, or $\Omega \subset \mathbb{R}^d$ in a continuous 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 $\Omega$ 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 $\Omega$ are coupled via an elastic energy, which is a quadratic form in the fields $\mathbf{u}(x)$. We also include in this quadratic term a parabolic confining potential of curvature $\mu > 0$. The absolute minimum of this first term is thus the flat, undisturbed, configuration $\mathbf{u}(x) = 0$. The second term is the quenched disorder, modeled by a random potential energy which couples directly to $\mathbf{u}(x)$. We thus consider the following model of an elastic manifold in a random potential given by its energy functional

$$ \mathcal{H}[\mathbf{u}] = \sum_{x,y} \mathbf{u}(x) \cdot (\mu \mathbf{1} - t\Delta)_{xy} \cdot \mathbf{u}(y) + \sum_x V(\mathbf{u}(x), x) \quad (1) $$

where here $x \in L^d \subset \mathbb{Z}^d$, $\mathbf{1}$ is an appropriate identity operator, and the matrix $-t\Delta_{xy}$ is required to be positive definite. Here $\Delta$ can be chosen as the discrete Laplacian in the hypercube $L^d$ with periodic boundary conditions. In that case its eigenmodes are plane waves $\sim e^{ikx}$ and we denote $\Delta(k)$ its eigenvalues, i.e. in $d=1$, $\Delta(k) = 2(\cos k - 1)$ with $k = 2\pi n/L$, $n = 0, \dots, 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\Delta(k)$, e.g. to more general elasticity (such as long range elasticity). They also extend to cases where $t\Delta_{xy}$ is a quadratic form defined on any graph $\Omega$. Finally, they also extend to the limit of the continuum manifold model, e.g. with the standard Laplacian $\Delta = \sum_{i=1}^d \frac{\partial^2}{\partial x_i^2}$ whose spectrum is given by $\Delta(\mathbf{k}) = -\mathbf{k}^2$. We thus use the notation $\int_k = \frac{1}{L^d} \sum_k \equiv \int \frac{d^d k}{(2\pi)^d}$ so our main formula are valid both for discrete and continuum (in the continuum $\sum_x \equiv \int d^d x$). We see from (1) that $\mu$ acts as a “mass” which, for the continuum model, leads to reducing the fluctuations beyond the scale $L_\mu = \sqrt{t/\mu}$.

Here we will consider $V(\mathbf{u}, x)$ to be a mean-zero Gaussian-distributed random potential in $\mathbb{R}^N \times \mathbb{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<fnref n="1" />:

$$ \overline{V(\mathbf{u}_1, x_1)V(\mathbf{u}_2, x_2)} = N B \left( \frac{(\mathbf{u}_1 - \mathbf{u}_2)^2}{N} \right) \delta^d(x_1 - x_2), \quad (2) $$

<footnote n="1">¹ We follow here the same notations as in [12–15]. Note the factor of 2 difference with the definition of $B(q)$ in [11]. Noting $\hat{B}$ the similar function there, we have $B(q) = \hat{B}(q/2)$. $Q$ defined below is the same object in both papers.</footnote>
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 $\mathbf{u}_0(x)$ (ground state) is non trivial and exhibits interesting statistically self-affine properties characterized by a roughness exponent: $\mathbf{u}_0(x) - \mathbf{u}_0(0) \sim |x|^{\zeta}$. The sample to sample fluctuations of the ground state energy (and, at finite temperature, of the free energy) grow with the scale as $\sim L^{\theta}$, with $\theta = d - 2 + 2\zeta$ 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 \sim (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 $\zeta = (4-d)/2$. Above $L_c$ metastability sets in, the response to perturbations involves jumps (shocks), with non-trivial roughness $\zeta$ 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 \to \infty$. Saddle point equations in replica space [14–17] lead to solutions exhibiting replica symmetry breaking (RSB) for $\mu < \mu_c$, which describe the glass phase where the manifold is pinned. The critical mass $\mu_c$ corresponds to the Larkin scale $L_c = \sqrt{\Gamma/\mu_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–20] and are valid in an expansion around $d = 4 - \epsilon$ (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 [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 $\zeta = 2/3$ and free energy fluctuation exponent, $\theta = 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–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].

## 1.2 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–39], with recent applications to statistical physics [10,34–36,38–41], neural networks and complex dynamics [42–46], string theory [47,48] and cosmology [49,50]. It is also of active current interest in pure and applied mathematics [51–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 \gg 1$ in [35,38,39], see also [37,50,52]. It was found that a sharp transition occurs from a 'simple'
landscape for $\mu > \mu_c$ (the same $\mu_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 $\mu < \mu_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–63] and quantum glasses [64–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 \to \infty$ [11]. The main goal of the present paper is to extend this study to the pinned elastic manifold. Hence we will study the $NL^d \times NL^d$ Hessian matrix

$$ K_{ix,jy}[\mathbf{u}] = \frac{\partial^2}{\partial u_i(x) \partial u_j(y)} \mathcal{H}[\mathbf{u}] = \delta_{ij}(\mu \mathbf{1} - t\Delta)_{xy} + \delta_{xy} \frac{\partial^2}{\partial u_i \partial u_j} V(\mathbf{u}(x), x) \quad (3) $$

in particular its density of eigenvalues $\rho(\lambda)$ normalized as $\int \rho(\lambda) d\lambda = 1$. An important feature of such matrix is its (block-)band structure visualized below:
<figure><img src="image_2.png" /><figcaption>Structure of the (NL × NL) Hessian matrix K in 1d discrete lattice model with L internal sites and periodic boundary conditions.

Only non-zero N × N blocks are indicated.</figcaption></figure>

where for $r = 1, \dots, L$ we have introduced $N \times N$ random matrices $\mathbf{W}^{(r)}$ with entries
$$ \mathbf{W}_{ij}^{(r)} = \left. \frac{\partial^2}{\partial u_i \partial u_j} V(\mathbf{u}(x), x) \right|_{x=x_r} $$

Our main focus here is the problem where the Hessian $K_{ix,jy}[\mathbf{u}_0]$ is chosen at the global minimal energy configuration $\mathbf{u}_0 \equiv \mathbf{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 $\mathbf{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 $\mathbf{u}(x)$, i.e. we may as well chose it at $\mathbf{u}(x) = \mathbf{0}$. The covariance structure of the random potential (2) implies, after a simple differentiation that entries of the matrices $\mathbf{W}^{(r)}$ are mean-zero Gaussian-distributed, independent for different $r$ and have the following covariance structure:

$$ \overline{\mathbf{W}}_{ij}^{(r)} \overline{\mathbf{W}}_{kl}^{(s)} = \delta_{rs} \frac{4}{N} B''(0) (\delta_{ij}\delta_{lk} + \delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}) \quad (4) $$

The matrices of such block-band type, with $\mathbf{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 Wegner 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 par-
ticular, the mean eigenvalue density for such type of models as $N \to \infty$ 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 $\mu \to \mu_{eff}$, which should be determined by a separate minimization procedure. The replacement $\mu \to \mu_{eff}$ 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.

Surprisingly, we were not able to trace the analysis of the density profile arising from that equation in the literature for the most interesting case of infinite system of size $L \to \infty$. 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 matrix-valued differential operator $K$ acting in the space of $N$-component vectors $\mathbf{f}(x) := (f_1(x), \dots, f_N(x))^T$ where, e.g. $x \in [0, L]^d$, by the following rule:

$$K\mathbf{f} = (\mu\mathbf{1} - t\Delta)\mathbf{f} + \hat{W}\mathbf{f}, \quad W_{i,j}(x) = \frac{\partial^2}{\partial u_i \partial u_j} V(\mathbf{u}(x), x) \quad (5)$$

with appropriate boundary conditions (e.g. periodic, or Dirichlet).

Without conditioning by global minimum the covariance structure of $\hat{W}$ is a natural analogue of (4):

$$\overline{W_{i,j}(x_1)W_{k,l}(x_2)} = \delta(x_1 - x_2) \frac{4}{N} B''(0) (\delta_{ij}\delta_{lk} + \delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}) \quad (6)$$

In particular, for $d = 1$ the operator can be visualized in the following form of an $N \times N$ matrix:

$$K = \begin{pmatrix} -t \frac{d^2}{dx^2} + \mu + W_{1,1}(x) & W_{1,2}(x) & \dots & W_{1,N}(x) \\ W_{1,2}(x) & -t \frac{d^2}{dx^2} + \mu + W_{2,2}(x) & \dots & W_{2,N}(x) \\ \dots & \dots & \dots & \dots \\ W_{N,N}(x) & W_{N,N-1}(x) & \dots & -t \frac{d^2}{dx^2} + \mu + W_{N,N}(x) \end{pmatrix} \quad (7)$$

Models of Such type are sometimes called the matrix Anderson models, and are essentially continuous versions of Wegner orbital models. Note that in such a case the spectral density in the infinite-volume limit can not be normalized.<fnref n="2" /> 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 \to \infty$ and yields an explicit form of the density profile.

<footnote n="2">² Recall the disorder-free model with $N = 1$ where $\rho(\lambda) = \int_{\mathbb{R}^d} \frac{d\mathbf{k}}{(2\pi)^d} \delta(\lambda - \mu - t\mathbf{k}^2) \propto t^{-1}((\lambda - \mu)/t)^{\frac{d-2}{2}} \theta(\lambda - \mu)$.</footnote>
## 2 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 $\mathbf{u}_0$ of the total energy:

$$ \overline{\mathcal{G}(x, y; \lambda, \mathbf{u}_0)} = \frac{1}{N} \sum_{i=1}^{N} \left( \frac{1}{\lambda - \mathcal{K}(\mathbf{u}_0)} \right)_{xi,yi} \qquad (8) $$

as well as its limit at coinciding points $\overline{\mathcal{G}(x, x; \lambda, \mathbf{u}_0)}$, which relates to the mean spectral density of the Hessian as

$$ \rho(\lambda) = \frac{1}{\pi} \lim_{\lambda \to 0^{-}} \mathrm{Im} \overline{\mathcal{G}(\lambda, \mathbf{u}_0)} , \quad \mathcal{G}(\lambda, \mathbf{u}_0) = \frac{1}{NL^d} \sum_x \overline{\mathcal{G}(x, x; \lambda, \mathbf{u}_0)} \qquad (9) $$

Employing the replica trick, we first show that for $N \to \infty$ (the limit being taken for a fixed value of $L^d$) the average Green's function is given by

$$ \overline{\mathcal{G}(x, y; \lambda, \mathbf{u}_0)} = \int_k \frac{e^{ik(x-y)}}{\lambda - \mu_{\mathrm{eff}} + t\Delta(k) - 4ipB''(0)} \qquad (10) $$

where the value of the parameter $p$ is determined by the following self-consistent “deformed semicircle” equation for the diagonal part

$$ \overline{\mathcal{G}(x, x; \lambda, \mathbf{u}_0)} = ip = \int_k \frac{1}{\lambda - \mu_{\mathrm{eff}} + t\Delta(k) - 4ipB''(0)} \qquad (11) $$

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 $\mathbf{u}_0$ is the parameter $\mu_{\mathrm{eff}}$. Below $\mu_{\mathrm{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 $\int_k$ applies both to the discrete models $\int_k = \frac{1}{L^d} \sum_k$ and the continuum limit $\int_k = \int \frac{d^d k}{(2\pi)^d}$. These equations are quite general and apply to basically arbitrary graph Laplacian matrices $t \Delta_{ix,jy}$ (even not translationally invariant ones, provided the formula are generalized by replacing $\int_k \frac{1}{A} \to \mathrm{tr} A^{-1}$, i.e. the trace in internal space of the inverse matrix).

### 2.1 Spectral Density of the Hessian at a Generic Point

As has been already mentioned, with setting $\mu_{\mathrm{eff}} = \mu$ the above expressions (10, 11) provides the mean resolvent and the mean spectral density $\rho(\lambda)$ 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
<figure><img src="image_3.png" /><figcaption>Fig. 1 Scaling function $r_c(\Lambda)$ for the Hessian spectral density for the $d=1$ continuum model plotted versus $\Lambda = t^{1/3} \frac{\lambda - \mu_{\text{eff}}}{3B''(0)^{2/3}}$, as defined in (12)</figcaption></figure>

matrix Anderson model (7), of infinite length $L \to \infty$ and the Laplacian spectrum $-\Delta(k) = k^2, -\infty < k < \infty$:

$$ \rho(\lambda) = \frac{1}{2\pi(t B''(0))^{1/3}} r_c \left( \Lambda = t^{1/3} \frac{\lambda - \mu}{3B''(0)^{2/3}} \right), \\ r_c(\Lambda) = \frac{w_r(\Lambda)^2}{4} \sqrt{\left( \frac{2}{w_r(\Lambda)} \right)^3 - 1} \qquad (12) $$

where

$$ w_r(\Lambda) = [1 + \sqrt{1 + \Lambda^3}]^{1/3} + [1 - \sqrt{1 + \Lambda^3}]^{1/3}, \qquad (13) $$

We have plotted in the Fig. 1 the parameter free scaling function $r_c(\Lambda)$. The spectral edge $\Lambda_e$ is given in this case by $\Lambda_e = -1$. The function $r_c(\Lambda)$ reaches its maximum at $\Lambda = 0$ and then decays at $\Lambda \gg 1$ as $r_c(\Lambda \gg 1) \sim \frac{1}{\sqrt{3\Lambda}}$. The latter regime corresponds to the spectral density $\rho(\lambda) = \frac{1}{2\pi} \frac{1}{\sqrt{t(\lambda - \mu)}}$ of the disorder-free operator $\mu - \frac{d^2}{dx^2}$ with the spectrum $\lambda = \mu + tk^2$.

In the case of 1D disordered elastic discrete chain with $-\Delta(k) = 2(1 - \cos k)$, $0 \le k \le 2\pi$ 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

$$ \rho(\lambda) = \frac{t}{2\pi B''(0)} r \left( \Lambda = \frac{\lambda - \mu}{2t}, y \right), y = \frac{t^2}{B''(0)} \qquad (14) $$

but the function $r(\Lambda, y)$ does not have a simple form for $y \sim 1$. However, in the limiting case of weak disorder $y \gg 1$ a very explicit characterisation is again possible. In this case the graph $r(\Lambda, y)$ has two spectral edges at $\Lambda_e^{(-)} = -\frac{3}{2}y^{-2/3}$ and $\Lambda_e^{(+)} = 2 + \frac{3}{2}y^{-2/3}$ and the density profile in the vicinity of the edges is simply related to the density profile $r_c(\Lambda)$ of 1D continuous system. Namely, in the vicinity of the left edge
<figure><img src="image_4.png" /><figcaption>Fig. 2 Blue: scaling function for the Hessian spectral density, $r(\Lambda, y)$ versus $\Lambda = \frac{\lambda - \mu}{2t}$ for the infinite discrete 1D chain given by Eq. (14), for $y = \frac{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)</figcaption></figure>

$$r(\Lambda, y) \approx y^{-2/3} r_c \left( \frac{2}{3} y^{2/3} \Lambda \right), \quad |\Lambda| \sim y^{-2/3} \qquad (15)$$

and essentially the same profile in the vicinity of the upper edge $|\Lambda - 2| \sim y^{-2/3}$. In between the edges, for any finite $0 < \Lambda < 2$ the profile for $y \gg 1$ is given by the "disorder-free" shape $r(\Lambda, y) \simeq 1/(y\sqrt{\Lambda(2-\Lambda)})$. 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 $\mu_{eff}$. 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.

## 2.2 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 power-law class when the derivative $B'(q)$ can be written as in [15]$^3$:

$$B'(q) = - \frac{B_0}{r_f^2 \left(1 + \frac{q}{\gamma r_f^2}\right)^\gamma}, \quad \gamma > 0 \qquad (16)$$

<footnote>$^3$ We follow the definitions and notations of [15] (Section II A and B) which are consistent with the original paper [12]. The parameter $\gamma$ is thus identical to the one defined in our previous work [11].</footnote>
As special limiting cases this class also includes the (i) exponential $B(q) = B_0e^{-q/r_f^2}$ as the limit $\gamma \to +\infty$, and (ii) the log-correlated case for $\gamma \to 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

$$B(q) = A(c+q)^{-\gamma+1} \qquad (17)$$

hence choosing $c = \gamma r_f^2$ and $A = B_0 c^\gamma / ((\gamma - 1)r_f^2)$.

Let us recall the main features of the replica solution [12,15] for $N \to \infty$ (restricting for simplicity to $d \le 4$). Let us first define the “Flory” roughness and the free energy fluctuation exponents

$$\zeta = \zeta_F(\gamma) = \frac{4-d}{2+2\gamma}, \quad \theta = \theta_F(\gamma) = d-2+2\zeta_F(\gamma) \qquad (18)$$

Then it was found that for $\mu < \mu_c(T)$, full replica symmetry breaking, FRSB, occurs whenever $\theta_F(\gamma) > 0$, and 1-step replica symmetry breaking, 1RSB, occurs when $\theta_F(\gamma) \le 0$. The first case, FRSB, thus always occurs for manifold of dimensions $2 < d < 4$, whereas for $0 \le d < 2$ it is possible whenever $0 < \gamma < \gamma_c(d) = \frac{2}{2-d}$. In that case the exponents $\zeta, \theta$ (which are defined in the limit $\mu \to 0$) are given by their Flory values. In the limit $\mu \to 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 $\gamma > \gamma_c(d)$. In that case there is a phase transition at $T_c(\mu)$ which survives for $\mu = 0$. It is worth mentioning that in the marginal case $\gamma = \gamma_c(d)$ this transition is of a continuous nature.

The exponents are $\theta = 0$ and $\zeta = \frac{2-d}{2}$ in both the high-T phase and the low-T 1RSB phase, with however different amplitudes.<fnref n="4" /> The special case $\gamma = \gamma_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 long-range correlated (Full RSB) $0 < \gamma < 1$, and short-range correlated (1-RSB) $\gamma > 1$ landscapes. For the manifold such as distinction thus also holds, however the critical value of $\gamma$ is not unity anymore, but equal to $\gamma_c(d) = \frac{2}{2-d}$. In particular for $d > 2$ one is always in the LRC case.<fnref n="5" /> 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.

## 2.3 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” $\mu_c > 0$ which controls the value of the parabolic confinement $\mu$ below which the replica symmetry breaking (RSB) occurs at zero temperature. Its value turns out to be given by the *positive* solution of

$$1 = 4B''(0) \int_k \frac{1}{(-t\Delta(k) + \mu_c)^2} \qquad (19)$$

which is controlled both by disorder strength and the elasticity matrix. For example, for 1D continuous system a simple calculation gives $\mu_c = (\frac{B''(0)}{\sqrt{f}})^{2/3}$. Our analysis shows that in

<footnote n="4">⁴ This is an artefact of $N = \infty$ for 1RSB and may not survive for $N = 1$, except maybe in the boundary case $\gamma = \gamma_c(d)$ (certainly it survives for $d = 0$, $\gamma_c = 1$ the log-correlated case).</footnote>

<footnote n="5">⁵ Note however that for $d > 4$ there is again a RS phase for weak disorder.</footnote>
the replica symmetric phase the lower spectral edge $\lambda_e^{(-)}$ of the Hessian (which we associate with the spectral gap) as a function of $\mu$ is given by

$$ \lambda_e^{(-)} = \mu - \mu_c + 4B''(0) \int_k \left[ \frac{1}{\mu - t\Delta(k)} - \int_k \frac{1}{\mu_c - t\Delta(k)} \right] \quad (20) $$

This formula immediately shows that for $\mu > \mu_c$ the Hessian spectrum is always gapped (from zero). Upon expanding for $\mu \to \mu_c$ and using (19) one immediately finds the gap vanishing quadratically at $\mu_c$. For $\mu < \mu_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 $(\mu_c - \mu)^4$ near the transition from below, with the super-universal exponent. For example, for the continuum model in dimension $d$ we get for $\mu = \mu_c(1 - \delta)$

$$ \lambda_e^{(-)} = \frac{\mu_c}{36B^{(3)}(0)^4} \left(\frac{4-d}{4}\right)^3 \left(B^{(4)}(0)B''(0) - \frac{2(3-d)}{4-d}B^{(3)}(0)^2\right)^2 \delta^4 + O(\delta^5) \quad (21) $$

In the second case of full RSB the gap is identically zero everywhere for $\mu \le \mu_c$.

We also obtain the average Green function (8) and find that at the edge of the spectrum it decays exponentially as $\sim e^{-|x-y|/L_c}$, 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, $L_c \sim 1/\mu_c^{1/2}$. This is thus reminiscent of the results of [71] although obtained there in a slightly different context (depinning). Remarkably, this property holds also for $\mu > \mu_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

$$ A(q) = \frac{2(3-d)}{4-d} (B'''(q))^2 - B''(q)B''''(q) < 0 \quad \Leftrightarrow \quad \text{Full RSB} \quad (22) $$

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 $\gamma$ if $d \ge 2$ and (ii) for $\gamma \le \gamma_c(d) = 2/(2-d)$ if $d \le 2$.

Finally, for $d = 1$ the above criterion specifies the covariance $B(q) = \frac{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) \sim 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 \to \infty$ 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–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(\mathbf{u}) = E = \text{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 Sect. 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 Sect. 4 we analyze the results. In the first Sect. 4.1 we obtain the spectral density and the Green function keeping $\mu_{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 $\mu_{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 Sect. 5 contains the conclusion.

## 3 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, \dots$ are denoted as $\sum_x \equiv \sum_{x=1}^{L^d}$, the sum over the first set of replica indices $\alpha, \gamma, \dots$ are denoted $\sum_\alpha \equiv \sum_{\alpha=1}^m$, the sums over the second set of replica indices $a, b, c, \dots$ are denoted $\sum_a \equiv \sum_{a=1}^m$, and similarly for the products. The indices $i = 1, \dots, N$ and the dot product is used in $\mathbb{R}^N$.

The notation Tr is the trace over all indices $x, i$ and $a$ or $\alpha$, i.e. over $\mathbb{R}^m \times L^d$ or $\mathbb{R}^n \times L^d$, e.g. $\text{Tr}A = \sum_{xa} A_{xa,xa}$. The notation tr is reserved for the traces over $a$ or $\alpha$ only, i.e. over $\mathbb{R}^m$ or $\mathbb{R}^n$, i.e. $\text{tr}A = \sum_a A_{aa}$.

### 3.1 Green's Function and the First Set of Replica

As the starting point of our approach, we introduce the resolvent of the Hessian $K(\mathbf{u})$ defined in (3), for a given generic configuration $\mathbf{u}(x)$ (not necessarily the minimum of the total energy) and in a given realization of the random potential $V(\mathbf{u}(x), x)$. The associated Green's function is then defined via

$$ \mathcal{G}(x, y; \lambda, \mathbf{u}) = \frac{1}{N} \sum_{i=1}^{N} \left( \frac{1}{\lambda - K(\mathbf{u})} \right)_{x_i, y_i}, \quad \mathcal{G}(\lambda, \mathbf{u}) = \frac{1}{NL^d} \sum_x \mathcal{G}(x, x; \lambda, \mathbf{u}) \quad (23) $$

Such Green's function admits then the following representation in terms of $m$ replicated Gaussian integrals over $N$-component real-valued vector fields $\phi_\alpha(x)$, with $\alpha = 1, \dots, m$:

$$ \mathcal{G}(x, y; \lambda, \mathbf{u}) = \lim_{m \to 0} \int_{\mathbb{R}^{Nm}} e^{-\frac{i}{2}\lambda \sum_{x,\alpha} \phi_{\alpha}^2(x)} e^{\frac{i}{2} \sum_{x,y,\alpha} \phi_{\alpha}(x) \cdot K(\mathbf{u}) \cdot \phi_{\alpha}(y)} \\ \times \left[ \frac{i}{mN} \sum_{\gamma} \phi_{\gamma}(x) \cdot \phi_{\gamma}(y) \right] \prod_{x,\alpha} \mathcal{D}\phi_{\alpha}(x) \quad (24) $$
where we assumed that $Im \lambda < 0$ and set the factor $(i/\pi)^{m/2} \to 1$ for $m = 0$. From this we calculate the mean spectral density of the Hessian eigenvalues “at a temperature $T$”, defined as

$$ \rho_T(\lambda) = \frac{1}{\pi} \lim_{\text{Im}\lambda \to 0^-} \text{Im} \overline{\langle \mathcal{G}(\lambda, \mathbf{u}) \rangle}_T \quad (25) $$

where the thermal averaged value of any functional $g(\mathbf{u})$ of a configuration $\mathbf{u}(x)$ at a temperature $T = \beta^{-1}$ is defined as $\langle g(\mathbf{u}) \rangle_T := \int g(\mathbf{u}) \pi_{\beta}(\mathbf{u}) \mathcal{D}\mathbf{u}(x)$, with $\pi_{\beta}(\mathbf{u}) = Z_{\beta}^{-1} e^{-\beta \mathcal{H}[\mathbf{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:

$$ \rho(\lambda) = \lim_{T \to 0} \rho_T(\lambda) = \frac{1}{\pi} \lim_{\text{Im}\lambda \to 0^-} \text{Im} \overline{\mathcal{G}(\lambda, \mathbf{u}_m)} \quad (26) $$

The problem therefore amounts to first calculating the disorder and thermal average

$$ \overline{\langle \mathcal{G}(x, y; \lambda, \mathbf{u}) \rangle}_T = \lim_{m \to 0} \int_{\mathbb{R}^{N_m}} e^{-\frac{i}{2}\lambda \sum_{x, \alpha} \phi_{\alpha}^2(x)} \overline{\langle e^{\frac{i}{2} \sum_{x, y, \alpha} \phi_{\alpha}(x) \cdot K(\mathbf{u}) \cdot \phi_{\alpha}(y)} \rangle}_T \\ \times \left[ \frac{i}{mN} \sum_{\gamma} \phi_{\gamma}(x) \cdot \phi_{\gamma}(y) \right] \prod_{x, \alpha} \mathcal{D}\phi_{\alpha}(x) \quad (27) $$

where

$$ \langle e^{\frac{i}{2} \sum_{x,y,\alpha} \phi_{\alpha}(x) \cdot K(\mathbf{u}) \cdot \phi_{\alpha}(y)} \rangle_T = Z_{\beta}^{-1} \int_{\mathbb{R}^N} \mathcal{D}\mathbf{u}(x) e^{\frac{i}{2} \sum_{x,y,\alpha} \phi_{\alpha}(x) \cdot K(\mathbf{u}) \cdot \phi_{\alpha}(y) - \beta \mathcal{H}[\mathbf{u}]} \quad (28) $$

and then, by performing the zero-temperature limit, to capture the contribution from the global minimum configuration only.

## 3.2 Average Green Function and Second Set of Replica

In the framework of the replica trick we represent the normalization factor $Z_{\beta}^{-1}$ in Eq. (28) formally as $1/Z_{\beta} = \lim_{n \to 0} Z_{\beta}^{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(\mathbf{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

$$ \overline{\langle e^{\frac{i}{2} \sum_{x,y,\alpha} \phi_{\alpha}(x) \cdot K(\mathbf{u}) \cdot \phi_{\alpha}(y)} \rangle}_T = \lim_{n \to 0} \int \prod_{x,a} \mathcal{D}\mathbf{u}_a(x) e^{-L_{n,m}[\mathbf{u}, \mathcal{E}]}} \quad (29) $$

with

$$ L_{n,m}[\mathbf{u}, \phi] = \mathcal{L}_m[\phi] + \mathcal{L}_{n,m}[\mathbf{u}, \phi] \quad (30) $$
where the **u**-independent part of the action is given by

$$
\mathcal{L}_m[\phi] = \sum_x \frac{B''(0)}{2N} \left[ \left( \sum_\alpha \phi_\alpha^2(x) \right)^2 + 2 \sum_{\alpha, \beta} (\phi_\alpha(x) \cdot \phi_\beta(x))^2 \right] \\
- \frac{i}{2} \sum_{x,y,\alpha} \phi_\alpha(x) \cdot (\mu \mathbf{1} - t \Delta)_{xy} \phi_\alpha(y) \tag{31}
$$

whereas both **u**- and **φ**-dependent part is

$$
\mathcal{L}_{n,m}[\mathbf{u}, \phi] = \frac{\beta}{2} \sum_{x,y,a} \mathbf{u}_a(x) \cdot (\mu \mathbf{1} - t \Delta)_{xy} \mathbf{u}_b(y) - N \frac{\beta^2}{2} \sum_{x,a,b} B \left( \frac{(\mathbf{u}_a(x) - \mathbf{u}_b(x))^2}{N} \right) \\
+ i \beta \sum_x \left( \sum_\alpha \phi_\alpha(x)^2 \right) \sum_a B' \left( \frac{(\mathbf{u}_a(x) - \mathbf{u}_1(x))^2}{N} \right) \\
+ \frac{2i\beta}{N} \sum_x \sum_a B'' \left( \frac{(\mathbf{u}_a(x) - \mathbf{u}_1(x))^2}{N} \right) \sum_\alpha ((\mathbf{u}_a(x) - \mathbf{u}_1(x)) \cdot \phi_\alpha(x))^2 \tag{32}
$$

We can thus rewrite

$$
\overline{\langle \mathcal{G}(x, y; \lambda, \mathbf{u}) \rangle}_T = \lim_{m,n \to 0} \int_{\mathbb{R}^{Nm}} \prod_{x,\alpha} \mathcal{D}\phi_\alpha(x) \int \prod_{x,a} \mathcal{D}\mathbf{u}_a(x) e^{-\frac{i}{2}\lambda \sum_{x,a} \phi_a^2(x) - L_{n,m}[\mathbf{u}, \mathbf{E}]}
$$
$$
\times \frac{i}{mN} \left[ \sum_\gamma \phi_\gamma(x) \cdot \phi_\gamma(y) \right] \tag{33}
$$

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

$$
dQ = \prod_{1 \le a \le b \le n} dQ_{ab}, \quad dP = \prod_{1 \le \alpha \le \beta \le m} dP_{\alpha\beta}, \quad dR = \prod_{a=1}^{n} \prod_{\alpha=1}^{m} dR_{\alpha\alpha} \quad (34)
$$

$$
d\sigma = \prod_{1 \le a \le b \le n} d\sigma_{ab}, \quad d\tau = \prod_{1 \le \alpha \le \beta \le m} d\tau_{\alpha\beta}, \quad d\eta = \prod_{a=1}^{n} \prod_{\alpha=1}^{m} d\eta_{\alpha\alpha} \quad (35)
$$

and then define $\mathcal{D}Q(x) = \prod_x dQ(x)$, etc. This allows us to use the formal identity

$$
1 = \int \prod_x \mathcal{D}Q(x) \mathcal{D}\sigma(x) \mathcal{D}P(x) \mathcal{D}\tau(x) \mathcal{D}R(x) \mathcal{D}\eta(x) e^{-\frac{\beta}{2} \sum_{x,a,b} \sigma_{ab}(x) (N Q_{ab}(x) - \mathbf{u}_a(x) \cdot \mathbf{u}_b(x))} \\
\times e^{\frac{i}{2} \sum_{x,a,\beta} \tau_{\alpha\beta}(x) (N P_{\alpha\beta}(x) - \phi_{\alpha}(x) \cdot \phi_{\beta}(x)) + \frac{1}{2} \sum_{x,a,\alpha} \eta_{\alpha\beta}(x) (N R_{\alpha\beta}(x) - \mathbf{u}_a(x) \cdot \phi_{\alpha}(x))} \quad (36)
$$

where the contours of integration are duly chosen. This identity can be then inserted inside
(33) effectively allowing to replace all scalar products **u**_a(x) · **u**_b(x), φ_α(x) · φ_β(x) 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

$$ \overline{\langle \mathcal{G}(x, x; \lambda, \mathbf{u}) \rangle_T} = \lim_{m,n \to 0} \int \prod_x \mathcal{D}Q(x) \mathcal{D}\sigma(x) \mathcal{D}P(x) \mathcal{D}\tau(x) \mathcal{D}R(x) \mathcal{D}\eta(x) \left[ \frac{i \mathrm{tr} P(x)}{mN} \right] e^{-NL[Q, \sigma, P, \tau, R, \eta]} \quad (37) $$

where we have defined the action

$$ L[Q, \sigma, P, \tau, R, \eta] = L_m[P, \tau] + L_{n,m}[Q, \sigma, P] + \delta L[Q, \sigma, P, \tau, R, \eta] \quad (38) $$

$$ \begin{aligned} L_m[P, \tau] = & \sum_x \frac{B''(0)}{2} \left[ (\mathrm{tr} P(x))^2 + 2\mathrm{tr}(P(x)^2) \right] + \frac{i\lambda}{2} \mathrm{tr} P(x) \\ & + \frac{1}{2} \mathrm{Tr} \ln((\mu \mathbf{1} - t\Delta)\mathbf{1}_m - \tau \mathbf{1}) - \frac{i}{2} \sum_x \mathrm{tr}(\tau(x)P(x)) \end{aligned} \quad (39) $$

$$ \begin{aligned} L_{n,m}[Q, \sigma, P] = & \sum_x \frac{\beta\mu}{2} \mathrm{tr} Q(x) - \frac{\beta^2}{2} \sum_{a,b=1}^n B(Q_{aa}(x) + Q_{bb}(x) - 2Q_{ab}(x)) \\ & + \frac{1}{2} \mathrm{Tr} \ln((\mu \mathbf{1} - t\Delta)\mathbf{1}_n - \sigma \mathbf{1}) \\ & + \frac{\beta}{2} \sum_x \mathrm{tr}(\sigma(x)Q(x)) + i\beta \sum_x \mathrm{tr} P(x) \sum_{a=1}^n B'(Q_{aa}(x) \\ & + Q_{a1}(x) - 2Q_{11}(x)) \end{aligned} \quad (40) $$

The last piece $\delta L[Q, \sigma, P, \tau, R, \eta]$ 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}$, $\sigma_{ab}(x) = \sigma_{ab}$, $P_{\alpha\beta}(x) = P_{\alpha\beta}$, $\tau_{ab}(x) = \tau_{ab}$ at the saddle point. We will consider only this solution on physical grounds. Let us define again the quantity

$$ \mu_{\mathrm{eff}} = \mu - 2\beta \sum_{a=1}^{n} B'(Q_{aa} + Q_{11} - 2Q_{a1}) \quad (41) $$

Taking first the functional derivatives w.r.t. $\tau_{\alpha\beta}(x)$ and $P_{\alpha\beta}(x)$ we arrive at the saddle-point equations

$$ -i P_{\alpha\beta} = ((\mu \mathbf{1} - t\Delta)\mathbf{1}_m - \tau \mathbf{1})_{\alpha x, \beta x}^{-1} \quad (42) $$

$$ \tau_{\alpha\beta} = -4i B''(0) \left( P_{\alpha\beta} + \frac{1}{2} \delta_{\alpha\beta} \mathrm{tr} P \right) + (\lambda + \mu - \mu_{\mathrm{eff}}) \delta_{\alpha\beta} \quad (43) $$

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 $\tau$ in the replica space implies that the corresponding saddle point solutions must be actually proportional to the identity matrix: $P_{\alpha\beta} = p\delta_{\alpha\beta}$ and $\tau_{\alpha\beta} = \tau\delta_{\alpha\beta}$. In the limit $m \to 0$ one then finds that $\tau$ satisfies the equation

$$ \tau = -4i B''(0) p + \lambda + \mu - \mu_{\mathrm{eff}} \quad (44) $$
which when substituted to the corresponding equation for p yields the closed self-consistency equation for the latter:

$$ \begin{aligned} ip &= -((\mu \mathbf{1} - t \Delta)\mathbf{1}_m - \tau \mathbf{1})_{xx}^{-1} \\ &= -\int_k \frac{1}{\mu - t \Delta(k) - \tau} = \int_k \frac{1}{\lambda - \mu_{\text{eff}} + t \Delta(k) - 4ipB''(0)} \end{aligned} \quad (45) $$

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, \dots, 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 $(\text{tr} P)^2$ absent for GOE case. However for $m \to 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 = \frac{1}{\lambda - \mu_{\text{eff}} - 4ipB''(0)} \quad (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 Eq. (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 \neq 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:

$$ \overline{\mathcal{G}(x, y; \lambda, \mathbf{u}_0)} = \int_k \frac{e^{ik(x-y)}}{\lambda - \mu_{\text{eff}} + t\Delta(k) - 4ipB''(0)} \quad (47) $$

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 \to 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. $\sigma_{ab}(x)$ and $\tau_{\alpha\beta}(x)$ yields

$$ \beta Q_{ab} = (((\mu \mathbf{1} - t \Delta)\mathbf{1}_n - \sigma \mathbf{1})_{ax,bx}^{-1} = \int_k G_{ab}(k) \quad (48) $$

$$ \sigma_{ab} = 2\beta \left( \sum_c B'(\chi_{ac}) - B'(\chi_{ab}) \right) \Leftrightarrow \sigma_{ab} = -2\beta B'(\chi_{ab}) \quad \text{and} \quad \sum_b \sigma_{ab} = 0 \quad (49) $$

where we define, as in [12,14,15]

$$ G_{ab}(k) = ((\mu - t \Delta(k))\mathbf{1}_n - \sigma)_{ab}^{-1} \quad (50) $$

$$ \chi_{ab} = T \int_k (G_{aa}(k) + G_{bb}(k) - 2G_{ab}(k)) \quad (51) $$

We will recall briefly the analysis of these equations below, as needed.
# 4 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 $\mu_{\text{eff}}$ takes some value at $T = 0$. That gives the shape of the spectral density $\rho(\lambda)$, up to a global shift of $\lambda$. 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 $\mu_{\text{eff}}$ as a function of $\mu$, which allows to determine the location of the edge of the spectrum.

## 4.1 The Spectral Density and Green's Function

### 4.1.1 General Formula, Larkin Mass and Lower Edge

We start with recalling the self-consistency equation for the diagonal part of the Green's function, the parameter $p$:

$$ \overline{G(x, x; \lambda, \mathbf{u}_0)} = ip = \int_k \frac{1}{\lambda - \mu_{\text{eff}} + t\Delta(k) - 4ipB''(0)}. \quad (52) $$

There are usually multiple solutions for $p$ and we must choose the branch such that for $\lambda \to \pm\infty$ one has $ip \sim \frac{1}{\lambda}$. For a discrete model the spectrum of the perturbed Laplacian is bounded and large $|\lambda|$ necessarily correspond to being outside of the spectrum. In the continuum model the same holds for large negative $\lambda$. In the range of $\lambda$ outside of the spectrum, $p$ is necessarily pure imaginary. When $\lambda$ 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$

$$ p = p_1 - ip_2, \quad \rho(\lambda) = \frac{1}{\pi} \text{Im}(ip) = \frac{1}{\pi} p_1 \quad (53) $$

which converts (52) after separating the real and imaginary parts, into two coupled equations which determine $p_1$ and $p_2$ as functions of $\lambda$:

$$ p_2 = \int_k \frac{\lambda - \mu_{\text{eff}} + t\Delta(k) - 4p_2 B''(0)}{(\lambda - \mu_{\text{eff}} + t\Delta(k) - 4p_2 B''(0))^2 + [4p_1 B''(0)]^2} := F(\lambda, p_2, p_1^2) \quad (54) $$

$$ p_1 = 4p_1 B''(0) \int_k \frac{1}{(\lambda - \mu_{\text{eff}} + t\Delta(k) - 4p_2 B''(0))^2 + [4p_1 B''(0)]^2} := p_1 G(\lambda, p_2, p_1^2) \quad (55) $$

The edge of the spectrum is for $\lambda = \lambda_e$ such that $p_1$ acquire a non zero value, hence it is determined by eliminating $p_2 = p_2^e$ in the system

$$ p_2^e = \int_k \frac{1}{\lambda_e - \mu_{\text{eff}} + t\Delta(k) - 4p_2^e B''(0)} \quad (56) $$

$$ 1 = 4B''(0) \int_k \frac{1}{(\lambda_e - \mu_{\text{eff}} + t\Delta(k) - 4p_2^e B''(0))^2} \quad (57) $$

since at the edge one can set $p_1 = p_1^e = 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 $\lambda, p_2$ and $u = p_1^2$, a straightforward expansion in powers of $\lambda - \lambda_e$ shows that just above the lower edge
$$p_1^2 \approx (\lambda - \lambda_e) \frac{\frac{\partial G}{\partial p_2} |_e \frac{\partial F}{\partial \lambda} |_e - \left( \frac{\partial F}{\partial p_2} |_e - 1 \right) \frac{\partial G}{\partial \lambda} |_e}{\left( \frac{\partial F}{\partial p_2} |_e - 1 \right) \frac{\partial G}{\partial u} |_e - \frac{\partial F}{\partial u} |_e \frac{\partial G}{\partial p_2} |_e} \quad (58)$$

implying a square-root singularity of the density of eigenvalues at the threshold in the generic case (where neither the numerator or denominator vanishes in (58)).

To further analyze these equations we introduce the Larkin mass $\mu_c > 0$ defined as the positive solution of

$$1 = 4B''(0) \int_k \frac{1}{(-t\Delta(k) + \mu_c)^2} \quad (59)$$

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 \le \mu < \mu_c$, with $\mu_c$ determined by (59). The Larkin mass exists whenever $B''(0) > 1/(4 \int_k \frac{1}{-t\Delta(k)})$ and when this is the case, it is unique. In the opposite case, the RS solution is stable for all values of $\mu$. Note that $\mu_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 $\mu_c > 0$. It is then easy to find a solution for the spectral edge $\lambda_e$. One sees that $p_2^e$ is now determined in terms of $\mu_c$ and that (56) is equivalent to

$$\lambda_e := \lambda_e^- = \mu_{\text{eff}} + 4p_2^e B''(0) - \mu_c, \quad p_2^e = - \int_k \frac{1}{-t\Delta(k) + \mu_c} \quad (60)$$

which determines $\lambda_e$ as a function of $\mu_c$. It turns out (see below) that this is always the lower edge, hence we denoted it $\lambda_e^-$. We discuss below how to obtain the other edge(s) when they exist.

### 4.1.2 Some Examples: Edges and the Spectral Density Shape

Let us study some examples, remembering that we denote $\int_k = \frac{1}{L^d} \sum_k \equiv \int \frac{d^d k}{(2\pi)^d}$ for either discrete or continuum models.

1. First recall that for a single-site (equivalently, zero dimensional $d = 0$) system with $L^d = 1$ the Eq. (52) gives

$$ip = \frac{1}{\lambda - \mu_{\text{eff}} - 4ipB''(0)} \quad \Leftrightarrow \quad ip = \frac{\lambda - \mu_{\text{eff}} + i\sqrt{16B''(0) - (\lambda - \mu_{\text{eff}})^2}}{8B''(0)} \quad (61)$$

Hence the Hessian spectral density from (53) is given by the semicircular law

$$\rho(\lambda) = \frac{1}{\pi} \text{Im}(ip) = \frac{1}{8\pi B''(0)} \sqrt{16B''(0) - (\lambda - \mu_{\text{eff}})^2} \theta_{\lambda}([\lambda_e^-, \lambda_e^+]) \\ \lambda_e^{\pm} = \mu_{\text{eff}} \pm 4\sqrt{B''(0)} \quad (62)$$

where $\theta_{\lambda}([a, b]) = 1$ if $\lambda \in [a, b]$ and zero otherwise. This is precisely the result obtained in [11]. On the other hand we can determine the edge using (56) and (57). First let us examine the Eq. (59). In that case it reads

$$1 = \frac{4B''(0)}{\mu_c^2} \quad (63)$$
The positive root is $\mu_c = 2\sqrt{B''(0)}$ and from (60) we find

$$ \lambda_e^- = \mu_{\text{eff}} - 4 \frac{B''(0)}{\mu_c} - \mu_c = \mu_{\text{eff}} - 4\sqrt{B''(0)} \quad (64) $$

and recover the lower threshold (1). If we now use the negative root of (63), $\mu_c = -2\sqrt{B''(0)}$ we obtain instead $\lambda_e^+ = \mu_{\text{eff}} + 4\sqrt{B''(0)}$, i.e. the upper edge !

It is easy to see that this is a general property. In other words the Eq. (59) may have several roots. Let us call $\mathcal{D}$ the set on the real axis supporting the spectrum of $-t\Delta$. It is easy to see that for the continuum model, which has $\mathcal{D} = [0, +\infty[$, Eq. (59) may have only a single root. In contrast, consider e.g. the infinite discrete lattice in $d = 1$ with $-t\Delta(k) = 2t(1 - \cos k)$ so that its spectrum is in $\mathcal{D} = [0, 4t]$. Clearly the r.h.s of (59) is infinite for $\mu_c \in -\mathcal{D} = [-4t, 0]$ and diverges at the edges of this interval. Hence one expects two roots, one for $\mu_c = \mu_c^+ > 0$, and, by symmetry, one for $\mu_c = \mu_c^- = -4t - \mu_c^+$. In the following we will always associate the positive root $\mu_c^+ = \mu_c$ with the Larkin mass. If the set $-\mathcal{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\Delta$ is required to be positive definite. Hence it corresponds to the lower edge of the Hessian.

2. As soon as we have $L \ge 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\Delta$ are $2t(1 - \cos \frac{2\pi j}{L})$, $j = 1, \dots, L - 1$. The Eq. (59) becomes

$$ \sum_{j=0}^{L-1} \frac{1}{(\mu_c + 2t \left(1 - \cos \frac{2\pi j}{L}\right))^2} = \frac{L}{4B''(0)} \quad (65) $$

It is easy to see that for very weak disorder, $\frac{B''(0)}{t^2L} \ll 1$, there are $2L$ roots to this equation, which we denote $\mu_c = \mu_c^{j, \pm}$, $j = 0, \dots, L - 1$ with

$$ \frac{1}{t} \mu_c^{j, \pm} = -2 \left(1 - \cos \frac{2\pi j}{L}\right) \mp 2 \sqrt{\frac{B''(0)}{t^2 L}} \left(1 + O\left(\frac{B''(0)}{t^2 L}\right)\right) \quad (66) $$

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

$$ \lambda_e = \mu_{\text{eff}} - \mu_c - \frac{4B''(0)}{L} \sum_{j=0}^{L-1} \frac{1}{\mu_c + 2t \left(1 - \cos \frac{2\pi j}{L}\right)} \quad (67) $$

Substituting here the values of $\mu_c^{j, \pm}$ found above we arrive, up to subdominant terms at weak disorder, to the corresponding edge values

$$ \lambda_e^{j, \pm} = \mu_{\text{eff}} + 2t \left(1 - \cos \frac{2\pi j}{L}\right) \pm 4\sqrt{\frac{B''(0)}{L}} + \dots \quad (68) $$

For $L = 1$ the above approximation is exact and one recovers the formula (1) valid for any disorder. For $L \ge 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 \frac{2\pi j}{L}) + \mu$ (in that case $\mu_{\text{eff}} = \mu$). 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\sqrt{\frac{B''(0)}{L}}$ is much smaller than their separations $\approx \frac{4t}{L}$. This gives the criterion

$$ \frac{4LB''(0)}{t^2} \ll 1 \qquad (69) $$

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, \pi$). The equation becomes

$$ \mu_c = -2t + 2tz, \quad y = \frac{2t^2}{B''(0)}, \quad \frac{1}{(1-z)^2} + \frac{1}{(1+z)^2} = y \qquad (70) $$

Hence there are two cases. Either disorder is weak $\frac{2t^2}{B''(0)} > 2$, and there are 4 real roots

$$ z_{\pm,+} = \pm\sqrt{\frac{y + \sqrt{4y + 1} + 1}{y}}, \quad z_{\pm,-} = \pm\sqrt{\frac{y - \sqrt{4y + 1} + 1}{y}} \qquad (71) $$

with $|z_{\pm,+}| > 1$ and $|z_{\pm,-}| < 1$ always. Or disorder is strong and only the two roots $z_{\pm,+}$ exist. These roots correspond to edges of the spectrum of the Hessian

$$ \lambda_e = \mu_{\text{eff}} - \mu_c - 4B''(0) \int_k \frac{1}{-t\Delta(k) + \mu_c} = \mu_{\text{eff}} + \frac{B''(0)}{t} \left[ y(1-z) - \frac{1}{z-1} - \frac{1}{z+1} \right] \qquad (72) $$

$$ = \mu_{\text{eff}} + \frac{B''(0)}{t} \left[ y \mp \frac{1}{2} z_{+,-} (2y - 1 + \epsilon\sqrt{1+4y}) \right] \qquad (73) $$

There are thus 4 edges (weak disorder) and 2 edges (strong disorder). The lowest edge is located at

$$ \begin{aligned} \lambda_e^- &= \mu_{\text{eff}} + \frac{B''(0)}{t} \left[ y - \frac{1}{2} z_{++} (2y - 1 + \sqrt{1+4y}) \right] \\ &= \mu_{\text{eff}} + \frac{B''(0)}{4t} \left( w^2 - \sqrt{w-1}(w+3)^{3/2} - 1 \right), \quad w = \sqrt{1+4y} \\ &= \sqrt{1 + \frac{8t^2}{B''(0)}} \end{aligned} \qquad (74) $$

We can now study the spectral density for the case $L=2$. It is given by $\rho(\lambda) = \frac{1}{\pi}\text{Im}(ip)$ where $ip$ satisfies the cubic equation

$$ 2ip = \frac{1}{\lambda - \mu_{\text{eff}} - 4ipB''(0)} + \frac{1}{\lambda - \mu_{\text{eff}} - 4t - 4ipB''(0)} \qquad (75) $$

The resulting density of states is plotted in Figs. 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.
3. Our next example is the continuum 1D line of infinite length $L \to +\infty$. Since the Laplacian spectrum for such a system is given by $-t\Delta(k) = tk^2$ with $k \in [0, +\infty[$ for such a case there is only one spectral edge in the system with disorder, the lower edge $\lambda_e^-$, determined from the Larkin mass $\mu_c$, i.e. the unique positive solution of

$$ \frac{1}{4B''(0)} = \int_{-\infty}^{\infty} \frac{dk}{2\pi} \frac{1}{(tk^2 + \mu_c)^2} = \frac{1}{4t^{1/2} \mu_c^{3/2}} \quad (76) $$

leading to

$$ \mu_c = \left( \frac{B''(0)}{\sqrt{t}} \right)^{2/3}, \quad \lambda_e^- = \mu_{\text{eff}} - \mu_c - 2 \frac{B''(0)}{\sqrt{t}} \mu_c^{-1/2} = \mu_{\text{eff}} - 3 \left( \frac{B''(0)}{\sqrt{t}} \right)^{2/3} \quad (77) $$

The spectral density in the interval $\lambda > \lambda_e^-$ is then given by $\rho(\lambda) = \frac{1}{\pi}\text{Im}(ip)$ where the complex $p$ is obtained by solving the following equation:

$$ ip = \int_{-\infty}^{\infty} \frac{dk}{2\pi} \frac{1}{\lambda - \mu_{\text{eff}} - tk^2 - 4ipB''(0)}. \quad (78) $$

Introducing new scaled variables $y, \tilde{p}, \tilde{\lambda}$ via

$$ y = \frac{t^2}{B''(0)}, \quad p = \frac{\tilde{p}}{2y^{1/6}\sqrt{B''(0)}}, \quad \lambda - \mu_{\text{eff}} = \frac{\tilde{\lambda}\sqrt{B''(0)}}{y^{1/6}} \quad (79) $$

and rescaling the integration variable as $k \to k y^{-1/3}$ the Eq. (78) takes the form

$$ i\tilde{p} = \frac{1}{\pi} \int_{-\infty}^{\infty} \frac{dk}{\tilde{\lambda} - 2i\tilde{p} - k^2} = \frac{-i}{\sqrt{\tilde{\lambda} - 2i\tilde{p}}} \operatorname{sgn}\left(\operatorname{Im}\sqrt{\tilde{\lambda} - 2i\tilde{p}}\right), \quad \operatorname{Im}\sqrt{\tilde{\lambda} - 2i\tilde{p}} \neq 0 \quad (80) $$

Taking the square and further introducing the variable $w$ and parameter $\delta$ by

$$ \tilde{p} = -\frac{i}{w}, \quad \tilde{\lambda} = 3\delta^{1/3} \quad (81) $$

the equation for $w$ attains especially simple form:

$$ w^3 + 3\delta^{1/3}w - 2 = 0 \quad (82) $$

The general theory of cubic equations then dictates that for $\delta < -1$ the Eq. (82) has only real solutions, hence $p \sim -iw^{-1}$ will be purely imaginary implying zero density of eigenvalues. This parameter range fully agrees with the position of lower spectral threshold $\lambda_e^-$ found in (77), which in the new variables reads $\tilde{\lambda}_e^- = -3y^{-1/6}$. As long as $\delta > -1$ there is one real root given by the Cardano formula in the form

$$ w_r = \left[ 1 + \sqrt{1 + \delta} \right]^{1/3} + \left[ 1 - \sqrt{1 + \delta} \right]^{1/3} \quad (83) $$

which is positive and decreases from 2 to 0 as $\delta$ increases from $-1$ to $+\infty$. There are also two complex conjugated solutions $w$ and $\bar{w}$. To find their imaginary part we use the Vieta's formulas:

$$ w_r + w + \bar{w} = 0, \quad w_r w \bar{w} = 2 $$
which gives $Re(w) = -\frac{w_r}{2}$ and $Im(w) = \pm\sqrt{\frac{2}{w_r} - \frac{w_r^2}{4}}$. It is then easy to see that the spectral density can be found explicitly and is given by:

$$ \rho(\lambda) = \frac{1}{2\pi\sqrt{B''(0)}y^{1/6}} \mathrm{Im}(w^{-1}) = -\frac{1}{2\pi\sqrt{B''(0)}y^{1/6}} \frac{\mathrm{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

$$ \delta = \frac{t}{[B''(0)]^2} \frac{(\lambda - \mu_{eff})^3}{27} = \left(-1 + \frac{\lambda - \lambda_e^-}{\mu_{eff} - \lambda_e^-}\right)^3, \quad y = \frac{t^2}{B''(0)} $$

this finally implies

$$ \rho(\lambda) = \frac{1}{2\pi(t B''(0))^{1/3}} r_c \left( \Lambda = t^{1/3} \frac{\lambda - \mu_{eff}}{3B''(0)^{2/3}} \right), \quad r_c(\Lambda) = \frac{w_r^2}{4} \sqrt{\left(\frac{2}{w_r}\right)^3 - 1} \quad (84) $$

We have plotted in the Fig. 1 the parameter free scaling function $r_c(\Lambda) = \sqrt{U - U^4}$, $U = \frac{w_r}{2}$ with $w_r$ given by (83) and $\delta = \Lambda^3$. It has the following asymptotics for large $\Lambda$ and for $\Lambda$ near the edge $\Lambda_e = -1$:

$$ r_c(\Lambda \gg 1) = \frac{1}{\sqrt{3\Lambda}} - \frac{5}{54\sqrt{3}\Lambda^{7/2}} + O(\Lambda^{-13/2}) \quad (85) $$

$$ r_c(\Lambda + 1 \ll 1) = \sqrt{\Lambda + 1} - \frac{5}{18}(\Lambda + 1)^{3/2} + O((\Lambda + 1)^{5/2}) \quad (86) $$

In particular, Eq. (84) then implies that

$$ \rho(\lambda) \approx \frac{1}{2\pi(t B''(0))^{1/3}} \sqrt{\frac{\lambda - \lambda_e^-}{\mu_{eff} - \lambda_e^-}} \quad (87) $$

showing the expected square-root singularity close to the spectral edge $\lambda = \lambda_e^-$. Moreover, it is easy to see that $\frac{dr_c}{dU} = 0$ for $U = 2^{-2/3}$, hence $w_r = 2^{1/3}$ corresponding according to (83) to $\delta = 0$. We conclude that $\Lambda = 0$ is exactly the point of the maximum for the scaled density profile $r_c(\Lambda)$, which is readily seen from Fig. 1.

4. Our last example is an infinite discrete lattice $d = 1$ with the number of sites $L \to \infty$. The Laplacian spectrum for such a system is given by $-t\Delta(k) = 2t(1 - \cos k)$ with $k \in [0, 2\pi]$. To determine the spectral edges we use the integrals

$$ \int_0^{2\pi} \frac{dk}{2\pi} \frac{1}{\cos k + x} = \frac{\operatorname{sgn}(x)}{\sqrt{x^2 - 1}}, \quad \int_0^{2\pi} \frac{dk}{2\pi} \frac{1}{(\cos k + x)^2} = \frac{|x|}{(x^2 - 1)^{3/2}} \quad (88) $$

for $|x| > 1$. This reduces finding the roots of (59) to solving the equation

$$ \frac{1}{4B''(0)} = \int_0^{2\pi} \frac{dk}{2\pi} \frac{1}{(2t(1 - \cos k) + \mu_c)^2} = \frac{|\mu_c + 2t|}{(\mu_c(\mu_c + 4t))^{3/2}} \quad (89) $$

for $\mu_c > 0$ or $\mu_c < -4t$. Denoting $r = 1 + \frac{\mu_c}{2t}$ and $y = \frac{t^2}{B''(0)}$ we rewrite the above equation as $\frac{|r|}{(r^2-1)^{3/2}} = y$ which implies a simple cubic equation $w^3 - y^{-2/3}w - 1 = 0$ for $w = (r^2-1)y^{2/3}$. Note that we have $|r| > 1$ for all allowed choices of $\mu_c$, hence need to look for a real positive solution $w > 0$ of this equation. According to general properties
of cubic equations, for $y > \frac{2}{3\sqrt{3}}$ our equation has only a single real root $w = w_c$ given by the Cardano formula:

$$w_c = \Delta^{1/3} + \frac{1}{3y^{2/3}} \Delta^{-1/3}, \quad \Delta = \frac{1}{2} + \sqrt{\frac{1}{4} - \frac{1}{27y^2}}$$

which is obviously positive as needed for our goals. For the parameter $\mu_c$ we then have two solutions. The positive one corresponds to the Larkin length:

$$\mu_c = \mu_c^+ = 2t \left( \sqrt{1 + w_c y^{-2/3}} - 1 \right) > 0 \quad (90)$$

and the second solution $\mu_c^- = -2t \left( \sqrt{1 + w_c y^{-2/3}} + 1 \right) \equiv -4t - \mu_c$ as expected by symmetry.

In the case $0 \le y \le \frac{2}{3\sqrt{3}}$ the cubic equations has all three real roots. Introducing the angle $\theta \in [0, \frac{\pi}{2}]$ such that $\cos \theta = \frac{3\sqrt{3}}{2} y$ the roots can be conveniently written in the so-called trigonometric form:

$$w_c^{(1)} = \frac{2}{\sqrt{3}y^{1/3}} \cos \left(\frac{\theta}{3}\right) > 0, \quad w_c^{(2)} = \frac{2}{\sqrt{3}y^{1/3}} \cos \left(\frac{\theta + 2\pi}{3}\right) < 0$$

$$w_c^{(2)} = \frac{2}{\sqrt{3}y^{1/3}} \cos \left( \frac{\theta + 4\pi}{3} \right) < 0,$$

so only $w_c^{(1)}$ can be used for the above procedure and yields $\mu_c^{\pm}$.

Finally, this gives us the two spectral edges as

$$\lambda_e^{\pm} = \mu_{\text{eff}} - \mu_c^{\mp} - 4B''(0) \int_k \frac{1}{2t(1 - \cos k) + \mu_c^{\mp}} = \mu_{\text{eff}} - \mu_c^{\mp} \\ -4B''(0) \frac{\operatorname{sgn}(\mu_c^{\mp} + 2t)}{\sqrt{\mu_c^{\mp}(4t + \mu_c^{\mp})}} \quad (91)$$

Let us give a simple example: $t = 2^{1/2}$, $B''(0) = 3^{3/2}$ which gives $y = \frac{2}{3^{3/2}}$, hence $\Delta = 1/2$ and $w_c = 2^{2/3}$ which eventually gives for the Larkin length $\mu_c^+ = 2\sqrt{2}$ and $\mu_c^- = -6\sqrt{2}$.

To calculate the spectral density profile one needs the following generalization of (88):

$$\frac{1}{2\pi} \int_{0}^{2\pi} \frac{dk}{2t \cos k + a} = \frac{\operatorname{sgn}(|a + \sqrt{a^2 - 4t^2}| - 2t)}{\sqrt{a^2 - 4t^2}} \quad (92)$$

valid for any real $t > 0$ and complex $a$ such that $|a + \sqrt{a^2 - 4t^2}| - 2t \neq 0$, excluding $a$ real in the interval $[-2t, 2t]$.

The spectral density in the interval $\lambda_e^- < \lambda < \lambda_e^+$ is then given by $\rho(\lambda) = \frac{1}{\pi} \text{Im}(ip)$ where the complex $p$ is obtained by solving the following equation, cf. (52):

$$ip = \int_{0}^{2\pi} \frac{dk}{2\pi} \frac{1}{\lambda - \mu_{\text{eff}} + 2t(\cos k - 1) - 4ipB''(0)} \quad (93)$$

$$= \frac{\operatorname{sgn}(|a(p) + \sqrt{a^2(p) - 4t^2}| - 2t)}{\sqrt{a^2(p) - 4t^2}} \quad (94)$$
where we denoted $a(p) = \lambda - 2t - \mu_{\text{eff}} - 4ipB''(0)$. We define the scaled variables

$$y = \frac{t^2}{B''(0)}, \quad \lambda - \mu_{\text{eff}} = 2t\Lambda, \quad p = \frac{\sqrt{y}}{2\sqrt{B''(0)}} P \quad (95)$$

which brings Eq. (93) in the dimensionless form

$$iP y = \frac{s}{\sqrt{(\Lambda - 1 - iP)^2 - 1}} \quad (96)$$

where $s = \pm 1$ with $s = \text{sgn}(|\Lambda - 1 - iP + \sqrt{(\Lambda - 1 - iP)^2 - 1}| - 1)$. The roots of (96) must satisfy the following equation for $\tilde{P} = iP$

$$\tilde{P}^4 + 2(1 - \Lambda)\tilde{P}^3 - (2 - \Lambda)\Lambda\tilde{P}^2 - \frac{1}{y^2} = 0 \quad (97)$$

The spectral density is then given in terms of the function $r(\Lambda, y) = \text{Im}(\tilde{P})$ as

$$\rho(\lambda) = \frac{t}{2\pi B''(0)} r\left(\Lambda = \frac{\lambda - \mu_{\text{eff}}}{2t}, y\right), \quad (98)$$

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 \to 0$, the system decouples in non-interacting zero-dimensional 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 \sim 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(\Lambda, 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) \to 0$, i.e. $y \to +\infty$. In this limit one expects that the spectral density should in a certain sense converge to the one of the purely elastic 1d system:

$$\begin{aligned} \rho(\lambda) &= \int_0^{2\pi} \frac{dk}{2\pi} \delta(\lambda - \mu_{\text{eff}} + 2t(\cos k - 1)) \\ &= \frac{1}{2t} \int_0^{2\pi} \frac{dk}{2\pi} \delta(\Lambda - (1 - \cos k)) = \frac{1}{2\pi t} \frac{1}{\sqrt{\Lambda(2-\Lambda)}}, \end{aligned} \quad (99)$$

which implies that for $y \gg 1$, $r(\Lambda, y) \simeq 1/(y\sqrt{\Lambda(2-\Lambda)})$. 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 \gg 1$ more accurately. Rescaling $\tilde{P} = q/y$ the Eq. (97) takes the form:

$$\frac{1}{y^2} q^4 + \frac{2}{y}(1 - \Lambda)q^3 - (2 - \Lambda)\Lambda q^2 - 1 = 0 \quad (100)$$

Now it is obvious that letting $y \to \infty$ for a fixed $0 < \Lambda < 2$ the above is reduced to the quadratic equation with purely imaginary roots $q = \pm \frac{i}{\sqrt{(2-\Lambda)\Lambda}}$. This solution yields precisely the density for the pure elastic case (99). However, it is also evident that in the vicinity of the points $\Lambda = 0$ or $\Lambda = 2$ such naive limit breaks down and requires a separate treatment. We illustrate it by providing analysis in the vicinity of $\Lambda = 0$, one for the region around $\Lambda = 2$ being fully analogous. A simple scaling argument demonstrates that the relevant vicinity of $\Lambda = 0$ is of the width $|\Lambda| \sim y^{-2/3}$ so that it
<figure><img src="image_5.png" /><figcaption>Fig. 3 Density of states $\rho(\lambda)$ for two coupled points, $L = 2$. We set $t = 1$. Left: weak disorder $B''(0) = 0.3$, the support splits in two disjoint intervals where the density of states is close to a semi-circle as for $L = 1$. Right: critical disorder $B''(0) = 1$ at which the two interval touch</figcaption></figure>

makes sense to introduce a new parameter $\delta = (\frac{2}{3}y^{2/3}\Lambda)^3$ and also introduce the scaled variable $w$ via $q = \frac{y^{1/3}}{w}$. Substituting this into (100) and taking the limit $y \to \infty$ one finds the equation for $w$ exactly given by the Eq. (82) studied by us in much detail above in our analysis of the 1d disordered continuum problem. We therefore conclude that for $y \to \infty$ and around $\Lambda = 0$ the scaled spectral density profile $r(\Lambda, y)$ for the 1D discrete model is simply given in terms of the scaled density profile of the continuum model $r_c(\Lambda)$ obtained in (84) as $r(\Lambda, y) = y^{-2/3}r_c(\frac{2}{3}y^{2/3}\Lambda)$. In particular, recalling that $r_c(\Lambda = -1) = 0$ in the continuum case, we find that the position of the left spectral threshold in the discrete case for $y \gg 1$ is given by $\Lambda_e^- = -\frac{3}{2}y^{-2/3}$. Close to this threshold the density increases as the square root $\sqrt{\Lambda - \Lambda_e^-}$, eventually reaching its maximal value $r(0, y) = y^{-2/3}2^{-4/3}\sqrt{3}$ exactly at $\Lambda = 0$ and then decaying for larger $\Lambda \gg y^{-2/3}$ in agreement with the asymptotics (85) as

$$r(\delta \gg 1) \sim y^{-2/3} \frac{1}{\sqrt{3(\frac{2}{3}y^{2/3}\Lambda)}} = \frac{1}{y} \frac{1}{\sqrt{2\Lambda}}$$

which precisely matches the $\Lambda \ll 1$ behaviour from the 'central part' $r(\Lambda, y) \simeq 1/(y\sqrt{\Lambda(2-\Lambda)})$. This demonstrates that for weak disorder, $y \gg 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, $|\Lambda| \sim y^{-2/3}$ and $|2-\Lambda| \sim 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.

## 4.2 Phases from Replica, Determination of $\mu_{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 $\mu_{eff}$ which corresponds to a global shift of the support of the Hessian spectral density.

In this section our aim is to determine $\mu_{eff}$, the missing information about the global position of the Hessian spectrum. As $\mu$ is varied the system can be in different phases (RS, 1RSB, FRSB) and the formula leading to $\mu_{eff}$ must be determined accordingly. In each case we
<figure><img src="image_6.png" /><figcaption>Fig. 4 Density of states $\rho(\lambda)$ for two coupled points, $L = 2$. We set $t = 1$. Left: $B''(0) = 1.2$. Right: larger disorder $B''(0) = 4$. The vertical bar arises from a switch in the solution of the cubic equation</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5 Scaling function for the Hessian spectral density, $r(\Lambda, y)$ versus $\Lambda = \frac{\lambda-\mu_{\text{eff}}}{2t}$ given by Eq. (98), for intermediate disorder $y = \frac{t^2}{B''(0)} = 1$. The two colors correspond to two roots of the quartic Eq. (97)</figcaption></figure>

first recall briefly the known replica saddle point solutions for the random manifold problem [12,15], i.e. the solutions to the Eqs. (48)–(50).

### 4.2.1 Replica-Symmetric Phase

Let us start with the replica symmetric (RS) phase, which occurs for $\mu > \mu_c$. Let us look for a replica symmetric solution of the saddle point equations (48)–(50).

$$ \sigma_{ab} = \sigma_c \delta_{ab} + \sigma, \quad G_{aa}(k) = \tilde{G}(k), \quad G_{a \neq b}(k) = G(k) \quad (101) $$

Note that the condition $\sum_b \sigma_{ab} = 0$ in this parametrization reads $\sigma_c + n\sigma = 0$ which in the replica limit $n \to 0$ implies that we can choose $\sigma_c = 0$. We also have $[G^{-1}]_{a \neq b} = -\sigma, [G^{-1}]_{aa} = \mu - t\Delta(k) - \sigma(1-n)$ and the inversion of the RS matrix gives:
$$ 
\begin{aligned}
G_{aa}(k) &= \frac{\mu - t\Delta(k) - \sigma - \sigma(n-2)}{(\mu - t\Delta(k))(\mu - t\Delta(k) - \sigma n)}, \\
G_{a \neq b}(k) &= -\frac{-\sigma}{(\mu - t\Delta(k))(\mu - t\Delta(k) - \sigma n)},
\end{aligned}
\quad (102) $$

or equivalently in the replica limit $n \to 0$

$$ G_{ab}(k) = \frac{1}{\mu - t\Delta(k)} \delta_{ab} + \frac{\sigma}{(\mu - t\Delta(k))^2}, \quad \chi_{a \neq b} = 2T \int_k \frac{1}{\mu - t\Delta(k)} \quad (103) $$

with $\chi_{aa} = 0$ by definition. The saddle point equation (50) leads then to the explicit formula for $\sigma$, which determines completely the solution

$$ \sigma = -2\beta B' \left( 2T \int_k \frac{1}{\mu - t\Delta(k)} \right) \quad (104) $$

and $Q_{ab} = T \int_k G_{ab}(k)$. As is well known [12] the RS solution is valid for $\mu > \mu_c(T)$ with (see Eq. (17) in [15])

$$ 1 = 4 \int_k \frac{1}{(\mu_c(T) - t\Delta(k))^2} B'' \left( 2T \int_k \frac{1}{\mu_c - t\Delta(k)} \right) \quad (105) $$

which gives, in the $T = 0$ limit, $\mu_c(T = 0) = \mu_c > 0$, i.e. the Larkin mass determined by

$$ 1 = 4B''(0) \int_k \frac{1}{(\mu_c - t\Delta(k))^2} \quad (106) $$

as anticipated in the previous Section.

We can now determine $\mu_{\text{eff}}$ and the edges of the Hessian in the RS phase. From (41) we obtain (for $n = 0$)

$$ 
\begin{aligned}
\mu_{\text{eff}} = \mu &- \frac{2}{T} [B'(\chi_{11}) + (n-1)B'(\chi_{a \neq b})] \to \mu - \frac{2}{T} [B'(0) \\
&- B' \left( 2T \int_k \frac{1}{\mu - t\Delta(k)} \right)] \\
\approx & T \to 0 \mu + 4B''(0) \int_k \frac{1}{\mu - t\Delta(k)}
\end{aligned}
\quad (107) $$

Substituting this value of $\mu_{\text{eff}}$ in (60) we thus obtain the final formula for the lower spectral edge $\lambda_e^{(-)}$ of the Hessian (which we associate with the spectral gap) as a function of $\mu$ in the RS phase

$$ \lambda_e^{(-)} = \mu - \mu_c + 4B''(0) \int_k \left[ \frac{1}{\mu - t\Delta(k)} - \int_k \frac{1}{\mu_c - t\Delta(k)} \right] \quad (108) $$

This formula immediately shows that the gap vanishes quadratically at $\mu_c$, i.e. upon expanding for $\mu > \mu_c$

$$ \lambda_e^{(-)} = 4B''(0) \int_k \frac{1}{(\mu_c - t\Delta(k))^3} (\mu - \mu_c)^2 + O((\mu - \mu_c)^3) \quad (109) $$

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\Delta(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 $\lambda_e^\alpha$ are obtained by considering all the real roots $\mu_c^\alpha$ 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 $\mu_c$ and $-4t - \mu_c$. That gives the upper edge for that model, and one can then write a single formula for both edges

$$ \lambda_e^{(\pm), \text{ld,discrete}} = \mu + 2t + 4B''(0) \int_k \frac{1}{\mu - t\Delta(k)} \pm \left( 2t + \mu_c + 4B''(0) \int_k \frac{1}{\mu_c - t\Delta(k)} \right) \quad (110) $$

which gives simple formula for the midpoint and the band width.

### 4.2.2 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 \to 0$ analysis<fnref n="6" /> in [12–15]. The off-diagonal part $\sigma_{a \neq b}$ is represented by a Parisi function $\sigma(v)$, $v \in [0, 1]$, which usually is $\sigma(v) = \sigma_c$ for $v \in [1, v_c]$, varies continuously with $v$ for $v \in [v_\mu, v]$ and is constant for $0 < v < v_\mu$. It is determined by

$$ \sigma(v) = -\frac{2}{T} B' \left( 2T \int_k (\tilde{G}(k) - G(k, v)) \right) \quad (111) $$

where from RSB replica matrix inversion one has

$$ \tilde{G}(k) - G(k, v) = \int_k \frac{1}{\mu - t\Delta(k) + \Sigma_c} + \int_v^{v_c} dw \frac{\sigma'(w)}{(\mu - t\Delta(k) + [\sigma](w))^2} \quad (112) $$

where $[\sigma](v)$ is defined by

$$ [\sigma](v) = v\sigma(v) - \int_0^v dw\sigma(w), \quad [\sigma](v_c) = \Sigma_c \quad (113) $$

Taking in (111) the derivatives w.r.t. $v$ and exploiting (112), one finds that for any interval of $v$ either (i) $\sigma(v)$ is constant or (ii) it satisfies the *marginality condition*

$$ 1 = 4B'' \left( 2T \int_k (\tilde{G}(k) - G(k, v)) \right) \int_k \frac{1}{(\mu - t\Delta(k) + [\sigma](v))^2} \quad (114) $$

In particular at the breakpoint $v = v_c$ one has

$$ 1 = 4B'' \left( 2T \int_k \frac{1}{\mu - t\Delta(k) + \Sigma_c} \right) \int_k \frac{1}{(\mu - t\Delta(k) + \Sigma_c)^2} \quad (115) $$

which by comparison with RS stability condition (105) implies that

$$ \Sigma_c(T) = \mu_c(T) - \mu \quad (116) $$

and at $T = 0$, as a function of $\mu_c$ the Larkin mass determined by (106), one has

$$ \Sigma_c = \mu_c - \mu. \quad (117) $$

<footnote n="6">⁶ in particular the Section VIII B of [14] (note the small misprint in the definition of $\sigma$ at the very beginning of the section there, paragraph above (8.4): it should be $TG_{ab}^{-1}(k) - (k^2 + m^2)\delta_{ab} = -\sigma_{ab}$). Note that $f = \hat{f} = -B$ in [12] and we follow [12] for the definition of $G(k)$ i.e. we do not absorb the $T$ inside it as done in [14,15].</footnote>
For use here and in the next Section let us define the notations for $l \ge 1$

$$I_l(x) = \int_k \frac{1}{(-t\Delta(k) + x)^l}, \quad I'_l(x) = -lI_{l+1}(x) \quad (118)$$

and somewhat abusively

$$I_0(x) = \int_k \log(-t\Delta(k) + x), \quad I'_0(x) = I_1(x) \quad (119)$$

One can calculate the solution for $[\sigma](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

$$\sigma(v) = -\frac{2}{T} B' \left( (B'')^{-1} \left( \frac{1}{4I_2(\mu + [\sigma](v))} \right) \right) \quad (120)$$

Taking a derivative of the above w.r.t. $v$ with the help of the identities

$$\frac{d[\sigma](v)}{dv} = v\sigma'(v), \quad \frac{d}{dv}f^{-1}(\phi(v)) = \frac{\phi'(v)}{f'(f^{-1}(\phi(v)))}$$

where $f^{-1}$ is the functional inverse of $f$, we obtain after rearranging and assuming $\sigma'(v) \ne 0$ the formula

$$v = -4T \frac{I_2^3(\mu + [\sigma](v))}{I_3(\mu + [\sigma](v))} B''' \left( (B'')^{-1} \left( \frac{1}{4I_2(\mu + [\sigma](v))} \right) \right) \quad (121)$$

which determines by inversion $[\sigma](v)$ as a function of $v$, which must be an increasing function. It is now convenient to introduce

$$A(v) = (B'')^{-1} \left( \frac{1}{4I_2(\mu + [\sigma](v))} \right), \quad \Rightarrow \quad I_2(\mu + [\sigma](v)) = \frac{1}{4B''(A(v))}$$

and define the function $F(b)$ via the relation

$$\frac{1}{I_3(x)} = F\left(\frac{1}{4I_2(x)}\right) \quad (122)$$

As the result, the relation (121) can be rewritten as

$$v = -\frac{T B'''(A(v)) F(B''(A(v)))}{16 [B''(A(v))]^3}$$

Taking yet another derivative w.r.t. $v$ and noticing that $\frac{dA}{dv} < 0$, this leads to the following condition for FRSB solution to exist

$$\frac{d}{dq} \frac{B'''(q) F(B''(q))}{(B''(q))^3} > 0 \quad (123)$$

where we used that $dq/dv < 0$. More precisely the condition for FRSB to hold in an interval $[v_\mu, v_c]$ is that (123) holds for $q \in [q(v_c), q(v_\mu)]$.

One may now notice that for a $d$-dimensional continuum model with Laplacian spectrum $-t\Delta(k) = t\mathbf{k}^2$ and $dk \sim |\mathbf{k}|^{d-1}d|\mathbf{k}|$ the behaviour of the integrals $I_{l>d/2}(x)$ in (118) for $x \to 0$ is dominated by the infrared ($|\mathbf{k}| \to 0$) limit and is given by $I_{l>d/2}(x) \sim x^{-(l-d/2)}$. Taking $d < 4$ we then see that (122) implies in the limit of small $\mu$ and small $b$ the behaviour

$$F(b) \sim b^{\frac{6-d}{4-d}} \quad (124)$$
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).

Note that the free energy fluctuation exponent $\theta = \theta_F = d - 2 + \frac{4-d}{1+\gamma} > 0$ and the FRSB self-energy behaves as $[\sigma](v) \sim v^{2/\theta}$ at small $v$ (for $\mu = 0$).

Let us now calculate $\mu_{\text{eff}}$ for the FRSB solution. For this we first set $v = v_c$ in (112) and (113) and get the relations

$$ \tilde{G}(k) - G(k, v_c) = \int_k \frac{1}{\mu - t\Delta(k) + \Sigma_c}, \quad \Sigma_c = v_c\sigma(v_c) - \int_0^{v_c} dw\sigma(w) \quad (125) $$

Now inserting the FRSB form into the definition (41) we get

$$ \mu_{\text{eff}} = \mu - 2\beta \left( B'(0) - \int_0^{v_c} dv B' \left( \int_k \tilde{G}(k) - G(k, v) \right) \right) \quad (126) $$

$$ -(1 - v_c) B' \left( \int_k \tilde{G}(k) - G(k, v_c) \right) \quad (127) $$

which can be further rewritten using (125) and the definition of $\sigma(v)$ in (111) as:

$$ \mu_{\text{eff}} = \mu - 2\beta B'(0) - \int_0^{v_c} dv \sigma(v) + v_c \sigma(v_c) + 2\beta B' \left( 2T \int_k \frac{1}{\mu - t\Delta(k) + \Sigma_c} \right) \quad (128) $$

$$ = \mu + \Sigma_c - \frac{2}{T} \left( B'(0) - B' \left( 2T \int_k \frac{1}{\mu - t\Delta(k) + \Sigma_c} \right) \right) \quad (129) $$

In the limit $T \to 0$, and recalling that $\Sigma_c = \mu_c - \mu$ we find

$$ \mu_{\text{eff}} = \mu_c + 4B''(0) \int_k \frac{1}{\mu_c - t\Delta(k)} \quad (130) $$

This has the same form as the RS formula (107) where one replaces $\mu$ by the Larkin mass $\mu_c$, i.e. it can be interpreted as the mass $\mu$ freezing at $\mu_c$, that is retaining for $\mu < \mu_c$ its critical value.

Let us now determine the lower edge of the Hessian. From (60) we obtain upon inserting (130)

$$ \lambda_e^- = \mu_{\text{eff}} - \mu_c - 4B''(0) \int_k \frac{1}{\mu_c - t\Delta(k)} = 0 \quad (131) $$

Hence the lower edge of the Hessian remains frozen at zero within the FRSB phase for all values of $\mu$. For models with more than one edge, their positions can be found from the other roots of the Eq. (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 $\mu_c$.

### 4.2.3 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 \le 2$ and $\gamma \ge \gamma_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 $\sigma_1$ and $\Sigma_c$ in terms of which we construct a $v_c \times v_c$ matrix $\sigma_d$ with entries $(\sigma_d)_{ab} = -\Sigma_c \delta_{ab} + \sigma_1$. The full $n \times n$ matrix $\sigma$ has $n/v_c$ identical
diagonal blocks $\sigma_d$, all entries being equal to the value $\sigma_0$ outside those blocks. The constraint $\Sigma_b\sigma_{ab} = 0$ then yields the relation in the $n \to 0$ limit:

$$ - \Sigma_c + v_c\sigma_1 + (n - v_c)\sigma_0 = 0 \Rightarrow v_c(\sigma_1 - \sigma_0) = \Sigma_c \quad (132) $$

Inversion of the matrix $G^{-1} = \mu - t\Delta(k) + \sigma$ produces $n \times n$ matrix $G$ with the diagonal $v_c \times v_c$ blocks $G_d$ having entries $(G_d)_{ab} = (\tilde{G} - G_1\delta_{ab} + G_1$ and outside those blocks $G$ has identical entries $G_0$. The entries in the limit $n \to 0$ remembering (132) are given by relations:

$$ G_0 = \frac{\sigma_0}{(\mu - t\Delta(k))^2}, \quad \tilde{G} - G_0 = \left(1 - \frac{1}{v_c}\right) \frac{1}{\mu - t\Delta(k) + \Sigma_c} + \frac{1}{v_c} \frac{1}{\mu - t\Delta(k)} \quad (133) $$

and

$$ \tilde{G} - G_1 = \frac{1}{\mu - t\Delta(k) + \Sigma_c}, \quad (134) $$

which according to (51) leads to

$$ \chi_0 = \frac{2T}{v_c} \left( \int_k \frac{1}{\mu - t\Delta(k)} - \int_k \frac{1 - v_c}{\mu - t\Delta(k) + \Sigma_c} \right), \quad \chi_1 = 2T \int_k \frac{1}{\mu - t\Delta(k) + \Sigma_c} \quad (135) $$

To determine the equilibrium values of the parameters involved we rely upon the expression for the free energy $\Phi(T)$ associated with the model given$^7$

$$ \Phi(T) = \frac{1}{2T} (v_c B(\chi_0) + (1 - v_c) B(\chi_1)) + \frac{T}{2} \frac{1 - v_c}{v_c} \int_k \left( \frac{\Sigma_c}{\mu - t\Delta(k) + \Sigma_c} - \log \left( \frac{\mu - t\Delta(k) + \Sigma_c}{\mu - t\Delta(k)} \right) \right) \quad (136) $$

Taking a derivative of the free energy w.r.t. $\Sigma_c$ leads to

$$ \Sigma_c = -2 \frac{v_c}{T} (B'(\chi_1) - B'(\chi_0)) \quad (137) $$

Let us consider the $T = 0$ limit. Denoting

$$ v_c = vT, \quad Q := \frac{2}{v}(I_1(\mu) - I_1(\mu + \Sigma_c)) \quad (138) $$

in terms of the integrals defined in (118) and noticing that in this limit $\chi_0 \to Q$ and $B(\chi_1)/2T \to B'(0)I_1(\mu + \Sigma_c)$ one gets

$$ \Phi(0) = B'(0)I_1(\mu) + \frac{v}{2}(B(Q) - B(0) - QB'(0)) - \frac{1}{2v}F_{\mu}(\Sigma_c) \quad (139) $$

where we have defined

$$ F_{\mu}(x) = I_{0}(\mu + x) - I_{0}(\mu) - xI_{1}(\mu + x) \quad (140) $$

with $F_\mu(0) = F'_\mu(0) = 0$ and $F''_\mu(0) = I_2(\mu)$.

Upon derivation of the zero-temperature free energy w.r.t. $\Sigma_c$ (cancelling the common factor $I_2$) and $v$ one obtains the following system of equations:

$$ \Sigma_c = 2v(B'(Q) - B'(0)) \quad (141) $$

$$ \frac{1}{v^2} F_{\mu}(\Sigma_c) = B(0) + QB'(Q) - B(Q) \quad (142) $$

<footnote>$^7$ see Eq. (170) in Section III.D p 17 or the arXiv version, in the paper [15] (up to an irrelevant constant).</footnote>
which should be augmented with the definition of $Q$ in (138). For small $Q > 0$ we have $\Sigma_c \simeq 2vB''(0)Q$ and substituting in (138) we then find that the transition to the phase with nonzero value of $Q$ occurs at $\mu = \mu_c$ determined by

$$1 = 4B''(0)I_2(\mu_c) \qquad (143)$$

which identifies with the definition of the Larkin mass, cf. (59).

We can now give the formula for $\mu_{\text{eff}}$ in the 1RSB phase. From (41) we have, inserting the one-step RSB ansatz

$$\mu_{\text{eff}} = \mu - \frac{2}{T}(B'(0) - v_c B'(\chi_0) - (1 - v_c)B'(\chi_1)) \qquad (144)$$

which in the limit $T \to 0$ yields

$$\mu_{\text{eff}} = \mu + 4B''(0)I_1(\mu) + 2v(B'(Q) - B'(0) - QB''(0)) \qquad (145)$$

Recalling from (60) that

$$\lambda_e = \mu_{\text{eff}} - 4B''(0)I_1(\mu_c) - \mu_c \qquad (146)$$

we finally obtain, within the 1RSB phase

$$\lambda_e = \mu - \mu_c + 4B''(0)(I_1(\mu) - I_1(\mu_c)) + 2v(B'(Q) - B'(0) - QB''(0)) \qquad (147)$$

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 $\Sigma_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 $\delta > 0$,

$$\mu = \mu_c(1 - \delta), \quad Q = \sum_{n \ge 1} Q_n \delta^n, \quad v = \sum_{n \ge 0} v_n \delta^n \qquad (148)$$

and solve order by order. It is convenient in the calculation to use that $I_l(x) = (-1)^{l-1}I_0^{(l)}(x)/(l-1)!$ for $l \ge 1$. We give only the lowest order

$$v_0 = -\frac{B'''(0)}{16B''(0)^3 I_3(\mu_c)}, \quad Q_1 = -\frac{8\mu_c B''(0)^2 I_3(\mu_c)}{B'''(0)} \qquad (149)$$

recalling that $B'''(0) < 0$. To this order one finds that the edge $\lambda_e$ vanishes to order $O(\delta^2)$. To find the first non-vanishing order, $O(\delta^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 rearrangements using (143), up to $O(\delta^5)$ terms

$$\lambda_e = \frac{\mu_c}{36B^{(3)}(0)^4} \left( \frac{\mu_c I_3(\mu_c)}{I_2(\mu_c)} \right)^3 \left( B^{(4)}(0)B''(0) - 3\left(1 - \frac{I_2(\mu_c)I_4(\mu_c)}{2I_3(\mu_c)^2}\right)B^{(3)}(0)^2 \right)^2 \delta^4 \qquad (150)$$

This result is very general, for any discrete or continuum model. For the $d=0$ 'single particle' model, $I_l(\mu_c) = \mu_c^l$ for $l \ge 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\Delta(k) = k^2$ (we set $t = 1$ for simplicity) and recall that $\int_k$ denotes $\int \frac{d^d k}{(2\pi)^d}$. The Eq. (141) hold with

$$F_{\mu}(\Sigma_c) = \mu^{d/2} \mathcal{F}\left(\frac{\Sigma_c}{\mu}\right), \quad \mathcal{F}(x) = \int \frac{d^d k}{(2\pi)^d} \left(\ln\left(1 + \frac{x}{1+k^2}\right) - \frac{x}{1+k^2+x}\right) \quad (151)$$

$$Q = \frac{2}{v} \mu^{d/2-1} \mathcal{G}\left(\frac{\Sigma_c}{\mu}\right), \quad \mathcal{G}(x) = \int \frac{d^d k}{(2\pi)^d} \left(\frac{1}{k^2+1} - \frac{1}{k^2+1+x}\right) \quad (152)$$

Restricting our consideration to $d < 4$ one can see that $\mathcal{F}$ and $\mathcal{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

$$\frac{\mu_c I_3(\mu_c)}{I_2(\mu_c)} = \frac{4-d}{4}, \quad \frac{I_2(\mu_c) I_4(\mu_c)}{2 I_3(\mu_c)^2} = \frac{6-d}{3(4-d)} \quad (153)$$

hence we obtain for the continuum model in dimension $d$

$$\lambda_e = \frac{\mu_c}{36 B^{(3)}(0)^4} \left(\frac{4-d}{4}\right)^3 \left(B^{(4)}(0)B''(0) - \frac{2(3-d)}{4-d}B^{(3)}(0)^2\right)^2 \delta^4 + O(\delta^5) \quad (154)$$

Let us study $d = 1, 2$ in more details. In $d = 2$ for the continuum model we have

$$\mathcal{F}(x) = \frac{1}{4\pi}(x - \log(1 + x)), \quad \mathcal{G}(x) = \frac{1}{4\pi}\log(1 + x), \\ I_1(\mu) - I_1(\mu_c) = \frac{1}{4\pi}\log\frac{\mu_c}{\mu} \quad (155)$$

and the transition occurs at

$$\mu_c = \frac{B''(0)}{\pi} \quad (156)$$

From the last equation in (141) we find that $\Sigma_c = \mu(e^{2\pi Qv} - 1)$, and substituting into the other two equations we obtain the system

$$\mu(e^{2\pi Qv} - 1) = 2v(B'(Q) - B'(0)) \quad (157)$$

$$\frac{1}{2\pi v}(B'(Q) - B'(0) - \mu\pi Q) = B(0) + QB'(Q) - B(Q) \quad (158)$$

The second equation allows to obtain $v = v(Q)$ and inserting 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 simple. For $\mu < \mu_c = c^2/\pi$ it reads

$$v = \frac{c}{2\pi}, \quad Q = \frac{1}{c} \log \frac{c^2}{\mu\pi} \quad (159)$$

Inserting into (147) we find, for the exponential case, that the edge of the spectrum of the Hessian is exactly at zero

$$\lambda_e = 0 \quad (160)$$
for all $\mu \le \mu_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

$$ 
\begin{aligned}
\mathcal{F}(x) &= \frac{2+x}{2\sqrt{1+x}} - 1, \quad \mathcal{G}(x) = \frac{x}{2(1+x+\sqrt{1+x})}, \\
I_1(\mu) - I_1(\mu_c) &= \frac{1}{2\sqrt{\mu}} - \frac{1}{2\sqrt{\mu_c}}
\end{aligned}
\qquad (161) $$

and $I_2(\mu_c) = \frac{1}{4}\mu_c^{-3/2}$, leading to $\mu_c = (B''(0))^{2/3}$. Let us choose

$$ B(q) = \frac{A}{c+q} \qquad (162) $$

Then

$$ \mu_c = \frac{(2A)^{2/3}}{c^2}, \quad v_0 = (2A)^{-1/3}, \quad Q_1 = \frac{c}{2} \qquad (163) $$

We must now solve the equations

$$ \Sigma_c = \frac{2vA}{c^2} \left( 1 - \frac{c^2}{(c+Q)^2} \right) \qquad (164) $$

$$ \frac{1}{v^2} \left( \frac{2\mu + \Sigma_c}{2\sqrt{\mu + \Sigma_c}} - \mu^{1/2} \right) = \frac{AQ^2}{c(c+Q)^2} \qquad (165) $$

$$ Qv = \mu^{-1/2} - (\mu + \Sigma_c)^{-1/2} \qquad (166) $$

It is convenient to introduce the following variables and parameters:

$$ \frac{c}{c+Q} = x, \quad \sqrt{\mu + \Sigma_c} = y, \quad Qv = z, \quad \Omega = \frac{2A}{c^3} $$

in terms of which the above system takes the form

$$ y^2 = \mu + \Omega zx(1+x) \qquad (167) $$

$$ \frac{(y - \sqrt{\mu})^2}{y} = \Omega z^2 x^2 \qquad (168) $$

$$ z = \frac{y - \sqrt{\mu}}{\sqrt{\mu}y} \qquad (169) $$

Substituting the last of those equations to the second one and remembering that for $\Sigma_c > 0$ we have $y > \sqrt{\mu}$ we see that the second equation takes the form $y = \frac{\Omega}{\mu}x^2$ implying further that $z = \frac{1}{\sqrt{\mu}} - \frac{\mu}{\Omega x^2}$. Substituting these relations into the first equation we see that it can be brought to the form

$$ \left(\frac{\Omega}{\mu}x^2\right)^2 - \mu = \sqrt{\mu}\frac{1+x}{x}\left(\frac{\Omega}{\mu}x^2 - \sqrt{\mu}\right) \Leftrightarrow \left(\frac{\Omega}{\mu}x^2 - \sqrt{\mu}\right)\left(\frac{\Omega}{\mu}x^2 - \frac{\sqrt{\mu}}{x}\right) = 0 $$

The first solution $x^2 = \frac{\mu^{3/2}}{\Omega}$ is however not admissible since it corresponds to $y = \sqrt{\mu}$ and therefore to $\Sigma_c = 0$. The only nontrivial solution as $\mu$ is decreased below $\mu_c$ is provided then by the remaining root $x = \frac{\mu^{1/2}}{\Omega^{1/3}}$ and in the original variables finally yields the relations
$$v = (2A)^{-1/3}, \quad Q + c = \frac{(2A)^{1/3}}{\sqrt{\mu}}, \quad \mu + \Sigma = \frac{(2A)^{2/3}}{c^2} \quad (170)$$

Substituting in (147) we again find $\lambda_e = 0$, confirming marginality for this case.

## 4.3 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 $\mathbf{u}^\alpha(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 $\mu + \Sigma_c = \mu_c$), while the large scale statistics (e.g. of $\mathbf{u}(x) - \mathbf{u}(0)$ at large $x$) is controled, throughout the glass phase, by the small $[\sigma](v)$ behavior (with effective mass $\mu + [\sigma](v)$). Hence in the FRSB phase it is the small $v$ behavior, corresponding to distant states, which leads 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 $\lambda = \lambda_e^-$, i.e. at the lower edge we can write

$$\overline{G(x, y; \lambda = \lambda_e^-, \mathbf{u}_0)} = \int_k \frac{e^{ik(x-y)}}{\lambda_e^- - \mu_{\text{eff}} + t\Delta(k) - 4ipB''(0)} = - \int_k \frac{e^{ik(x-y)}}{-t\Delta(k) + \mu_c} \quad (171)$$

where we used (60). Hence we see that the averaged Green function decays exponentially $\sim e^{-|x-y|/L_c}$, with the characteristic length given by the Larkin length $L_c$. For the continuum model with short-range elasticity and weak disorder, $L_c \sim 1/\mu_c^{1/2}$. This is very reminiscent of the result found numerically in [71] in the context of depinning. Remarkably however, here this property holds also for $\mu > \mu_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.

## 5 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 (see e.g. the recent works on characterizing
the energy landscape complexity for such systems [10,83]). 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
$\mu > \mu_c$ the Hessian spectrum is gapped away from zero. At $\mu = \mu_c$ the gap vanishes, i.e.
the lowest eigenvalue is zero. For $\mu < \mu_c$ the saddle-point solution is 1RSB and we find that
the gap is non zero and vanishes as $\sim (\mu_c - \mu)^4$ near the transition. For long range disorder
cases we find that the gap vanishes identically for $\mu \le \mu_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 $\mu_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 \ge 1$. Progress in those directions is reported elsewhere
[83].

Finally, the most interesting but very challenging problem is to study the present model
by taking limits $N \to \infty$ and $L \to \infty$ 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 field-theoretical/supersymmetric methods which
proved to be instrumental in getting insights into spectra and eigenvectors of matrices of
banded type [84,85] could be used successfully in the present problem.

**Acknowledgements** YVF is grateful to Mira Shamis for a helpful discussion about the “deformed semicircle” in the context of Wegner orbital model. This research was initiated during YVF visit to Paris supported by the Philippe Meyer Institute for Theoretical Physics, which is gratefully acknowledged. The research at King’s
(YVF) was supported by EPSRC grant EP/N009436/1 “The many faces of random characteristic polynomials”. PLD acknowledges support from ANR grant ANR-17-CE30-0027-01 RaMa-TraF.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## A Analysis of $\delta L[Q, \sigma, P, \tau, R, \eta]$

We give here the last piece of the replicated action, omitted in the text.

$$ \delta L[Q, \sigma, P, \tau, R, \eta] = \frac{1}{2} \mathrm{Tr} \ln \left( \frac{(\mu \mathbf{1} - t \Delta) \mathbf{1}_n - \sigma \mathbf{1}}{\eta^T (\mu \mathbf{1} - t \Delta) \mathbf{1}_m - \tau \mathbf{1}} \right) \quad (172) $$

$$ -\frac{1}{2} \mathrm{Tr} \ln((\mu \mathbf{1} - t \Delta)\mathbf{1}_m - \tau \mathbf{1}) - \frac{1}{2} \mathrm{Tr} \ln((\mu \mathbf{1} - t \Delta)\mathbf{1}_n - \sigma \mathbf{1}) \quad (173) $$

$$ +2i\beta \sum_{a=1}^{n} B'' \left( \frac{Q_{aa} + Q_{a1} - 2Q_{11}}{2} \right) \left( (RR^T)_{aa} + (RR^T)_{11} - 2(RR^T)_{a1} \right) \quad (174) $$

$$ -\frac{1}{2} \sum_{x} \mathrm{tr}(\eta(x) R(x)) \quad (175) $$

At the saddle point $R = 0$ hence it vanishes. The main argument for that is very similar to the discussion in [11].

## References

1. Blatter, G., Feigel’man, M.V., Geshkenbein, V.B., Larkin, A.I., Vinokur, V.M.: Vortices in high-temperature superconductors. Rev. Mod. Phys. **66**, 1125–1388 (1994)

2. Le Doussal, P.: Novel phases of vortices in superconductors. Int. J. Mod. Phys. B **24**, 3855–3914 (2010). In: L. N. Cooper and D. Feldman (eds.) BCS: 50 years. World Scientific (2011)

3. Giamarchi, T., Le Doussal, P.: For review see, Statics and dynamics of disordered elastic systems. In: A.P. Young (ed.) “Spin glasses and Random fields” in Series on Directions in condensed matter physics vol 12. World Scientific (1998) [cond-mat/9705096]

4. Fisher, D.S.: Sliding charge-density waves as a dynamic critical phenomenon. Phys. Rev. B **31**, 1396–1427 (1985)

5. Le Doussal, P., Wiese, K.J., Chauve, P.: Two-loop functional renormalization group theory of the depinning transition. Phys. Rev. B **66**, 174201 (2002)

6. Rosso, A., Krauth, W.: Roughness at the depinning threshold for a long-range elastic string. Phys. Rev. E **65**, 025101 (2002)

7. Rosso, A., Le Doussal, P., Wiese, K.J.: Avalanche-size distribution at the depinning transition: a numerical test of the theory. Phys. Rev. B **80**, 144204 (2009)

8. Le Doussal, P., Wiese, K.J.: Avalanche dynamics of elastic interfaces. Phys. Rev. E **88**, 022106 (2013)

9. Nattermann, T., Stepanow, S., Tang, L.-H., Leschhorn, H.: Dynamics of interface depinning in a disordered medium. J. Phys. II (France) **2**, 1483–8 (1992)

10. Fyodorov, Y.V., Le Doussal, P., Rosso, A., Texier, C.: Exponential number of equilibria and depinning threshold for a directed polymer in a random potential. Ann. Phys. **397**, 1–64 (2018)
11. Fyodorov, Y.V., Le Doussal, P.: Hessian spectrum at the global minimum of high-dimensional random landscapes. J. Phys. A **51**(47), 474002 (2018)

12. Mezard, M., Parisi, G.: Manifolds in random media: two extreme cases. J. Phys. I France **2**, 2231–2242 (1992)

13. Giamarchi, T., Le Doussal, P.: Elastic theory of flux lattices in presence of weak disorder. Phys. Rev. B **52**, 1242–1270 (1995)

14. Le Doussal, Pierre, Wiese, Kay Joerg: Functional renormalization group at large N for disordered elastic systems, and relation to replica symmetry breaking. Phys. Rev. B **68**, 17402 (2003)

15. Le Doussal, P., Mueller, M., Wiese, K.J.: Cusps and shocks in the renormalized potential of glassy random manifolds: how functional renormalization group and replica symmetry breaking fit together. Phys. Rev. B **77**, 064203 (2008). (39 pages)

16. Balents, L., Bouchaud, J.-P., Mezard, M.: The large scale energy landscape of randomly pinned objects. J. Phys. I (France) **6**, 1007 (1996)

17. Mézard, M., Parisi, G.: Replica field theory for random manifolds. J. Phys. I (France) **1**, 809 (1991)

18. Fisher, D.S.: Interface fluctuations in disordered systems: $5 - \epsilon$ expansion and failure of dimensional reduction. Phys. Rev. Lett. **56**, 1964–1967 (1986)

19. Le Doussal, P.: Exact results and open questions in first principle functional RG. Ann. Phys. **325**, 49 (2009)

20. Le Doussal, P., Wiese, K.J., Chauve, P.: Functional renormalization group and the field theory of disordered elastic systems. Phys. Rev. E **69**, 026112 (2004). arXiv:cond-mat/0304614

21. Halpin-Healy, T., Zhang, Y.-C.: Kinetic roughening phenomena, stochastic growth, directed polymers and all that. Aspects of multidisciplinary statistical mechanics. Phys. Rep. **254**(4–6), 215–414 (1995)

22. Johansson, K.: Shape fluctuations and random matrices. Commun. Math. Phys. **209**(2), 437–476 (2000)

23. Calabrese, P., Doussal, P.L., Rosso, A.: Free-energy distribution of the directed polymer at high temperature. Europhys. Lett. **90**, 20002 (2010)

24. Dotsenko, V.: Bethe ansatz derivation of the Tracy-Widom distribution for one-dimensional directed polymers. Europhys. Lett. **90**, 20003 (2010)

25. Dotsenko, V.: Replica Bethe ansatz derivation of the Tracy-Widom distribution of the free energy fluctuations in one-dimensional directed polymers. J. Stat. Mech. **2010**, P07010 (2010)

26. Sasamoto, T., Spohn, H.: One-dimensional Kardar-Parisi-Zhang equation: an exact solution and its universality. Phys. Rev. Lett. **104**, 230602 (2010)

27. Amir, G., Corwin, I., Quastel, J.: Probability distribution of the free energy of the continuum directed random polymer in 1 + 1 dimensions. Commun. Pure Appl. Math. **64**, 466–537 (2011)

28. Fyodorov, Y.V., Le Doussal, P., Rosso, A.: Freezing transition in decaying burgers turbulence and random matrix dualities. EPL **90**, 60004 (2010)

29. Le Doussal, P., Mueller, M., Wiese, K.: Avalanches in mean-field models and the Barkhausen noise in spin-glasses. EPL **91**, 57004 (2010)

30. Longuet-Higgins, M.S.: Reflection and refraction at a random moving surface. II. Number of specular points in a Gaussian surface. J. Opt. Soc. Am. **50**, 845 (1960)

31. Halperin, B.I., Lax, M.: Impurity-band tails in the high-density limit. I. Minimum counting methods. Phys. Rev. **148**, 722–740 (1966)

32. Weinrib, A., Halperin, B.I.: Distribution of maxima, minima, and saddle points of the intensity of laser speckle patterns. Phys. Rev. B **26**, 1362–1368 (1982)

33. Freund, I.: Saddles, singularities, and extrema in random phase fields. Phys. Rev. E **52**, 2348–2360 (1995)

34. Annibale, A., Cavagna, A., Giardina, I., Parisi, G.: Supersymmetric complexity in the Sherrington-Kirkpatrick model. Phys. Rev. E **68**, 061103 (2003)

35. Fyodorov, Y.V.: Complexity of Random Energy Landscapes, Glass Transition, and Absolute Value of the Spectral Determinant of Random Matrices. Phys. Rev. Lett. **92**(24), 240601 (2004); Erratum *ibid* **93**, Issue 14, 149901(E)(2004)

36. Parisi, G.: Computing the number of metastable states in infinite-range models. In: Bovier, A., et al. (eds.) Les Houches Summer School, Session LXXXIII, vol. 83, pp. 295–329. Elsevier, Amsterdam (2005)

37. Bray, A.J., Dean, D.S.: Statistics of critical points of gaussian fields on large-dimensional spaces. Phys. Rev. Lett. **98**, 150201 (2007)

38. Fyodorov, Y.V., Williams, I.: Replica symmetry breaking condition exposed by random matrix calculation of landscape complexity. J. Stat. Phys. **129**(5–6), 1081–1116 (2007)

39. Fyodorov, Y.V., Nadal, C.: Critical behavior of the number of minima of a random landscape at the glass transition point and the Tracy-Widom distribution. Phys. Rev. Lett. **109**, 167203 (2012)

40. Fyodorov, Y.V., Le Doussal, P.: Topology trivialization and large deviations for the minimum in the simplest random optimization. J. Stat. Phys. **154**(1–2), 466–490 (2014)
41. Ros, V., Arous, G.Ben, Biroli, G., Cammarota, C.: Complex energy landscapes in spiked-tensor and simple glassy models: ruggedness, arrangements of local minima and phase transitions. Phys. Rev. X **9**, 011003 (2019)

42. Wainrib, G., Touboul, J.: Topological and Dynamical complexity of random neural networks. Phys. Rev. Lett. **110**, 118101 (2013)

43. Fyodorov, Y.V., Khoruzhenko, B.A.: Nonlinear analogue of the May-Wigner instability transition. Proc. Natl. Acad. Sci. (USA) **113**(25), 6827–6832 (2016)

44. Fyodorov, Y.V.: Topology trivialization transition in random non-gradient autonomous ODEs on a sphere. J. Stat. Mech. **2016**(12), 124003 (2016)

45. Ipsen, J.R.: May-Wigner transition in large random dynamical systems. J. Stat. Mech. **2017**, P093209 (2017)

46. Ipsen, J.R., Forrester, P.J.: Kac-Rice fixed point analysis for single- and multi-layered complex systems. J. Phys. A **51**(47), 474003 (2018)

47. Douglas, M.R., Shiffman, B., Zelditch, S.: Critical points and supersymmetric vacua, I. Commun. Math. Phys. **252**(1–3), 325–358 (2004)

48. Douglas, M.R., Shiffman, B., Zelditch, S.: Critical Points and supersymmetric vacua, III: string/M models. Commun. Math. Phys. **265**(3), 617–671 (2006)

49. Easther, R., Guth, A.H., Masoumi, A.: Counting Vacua in Random Landscapes. preprint hep-th arXiv:1612.05224 (2016)

50. Yamada, M., Vilenkin, A.: Hessian Eigenvalue distribution for a random Gaussian field. J. High Energy Phys. **2018**, 29 (2018). https://doi.org/10.1007/JHEP03(2018)029

51. Azais, J.-M., Wschebor, M.: Level Sets and Extrema of Random Processes and Fields. Wiley, Hoboken (2009)

52. Fyodorov, Y.V.: High-dimensional random fields and random matrix theory. Markov Process. Relat. Fields **21**(3), 483–518 (2015)

53. Auffinger, A., Arous, G Ben, Cerny, J.: Random matrices and complexity of spin glasses Commun. Pure. Appl. Math. **66**(2), 165–201 (2013)

54. Auffinger, A., Arous, G Ben: Complexity of random smooth functions on the high-dimensional sphere. Ann. Probab. **41**(6), 4214–4247 (2013)

55. Nicolaescu, L.I.: Complexity of random smooth functions on compact manifolds. Indiana Univ. Math. J. **63**, 1037–1065 (2014)

56. Subag, E., Zeitouni, O.: The extremal process of critical points of the pure p-spin spherical spin glass model. Probab. Theory Relat. Fields **168**(3–4), 773–820 (2017)

57. Subag, E.: The complexity of spherical p-spin models in second moment approach. Ann. Probab. **45**(5), 3385–3450 (2017)

58. Subag, E.: Free Energy Landscapes in Spherical Spin Glasses. arXiv:1804.10576 (2018)

59. Cammarota, V., Wigman, I.: Fluctuations of the total number of critical points of random spherical harmonics. Stoch. Process. Appl. **127**(12), 3825–3869 (2017)

60. Cheng, D., Schwartzman, A.: Expected number and height distribution of critical points of smooth isotropic Gaussian random fields. Bernoulli **24**(4B), 3422–3446 (2018)

61. Franz, S., Parisi, G., Urbani, P.: Universal spectrum of normal modes in low-temperature glasses: an exact solution. Proc. Nat. Acad. Sci. (USA) **112**(47), 14539–14544 (2015)

62. Mueller, M., Wyart, M.: Marginal stability in structural, spin and electron glasses. Annu. Rev. Condens. Matter Phys. **6**, 9 (2015)

63. Yan, L., Baity-Jesi, M., Mueller, M., Wyart, M.: Dynamics and correlations among soft excitations in marginally stable glasses. Phys. Rev. Lett. **114**, 247208 (2015)

64. Andreanov, A., Mueller, M.: Long-range quantum ising spin glasses at $T = 0$: gapless collective excitations and universality. Phys. Rev. Lett. **109**, 177201 (2012)

65. Giamarchi, T., Le Doussal, P.: Variational theory of elastic manifolds with correlated disorder and localization of interacting quantum particles. Phys. Rev. B **53**(22), 15206 (1996)

66. Cugliandolo, L.F., Giamarchi, T., Le Doussal, P.: Dynamic compressibility and aging in Wigner crystals and quantum glasses. Phys. Rev. Lett. **96**(21), 217203 (2006)

67. Schehr, G.: Low temperature specific heat of some quantum mean field glassy phases. Phys. Rev. B **71**(18), 184204 (2005)

68. Gurarie, V., Chalker, J.T.: Bosonic excitations in random media. Phys. Rev. B **68**, 134207 (2003)

69. Mezard, M., Parisi, G., Virasoro, M.: Spin Glass Theory and Beyond: An Introduction to the Replica Method and Its Applications. (World Scientific Lecture Notes in Physics) World Scientific Publishing Company (1986)

70. Mueller, M., Leuzzi, L., Crisanti, A.: Marginal states in mean-field glasses. Phys. Rev. B **74**(13), 134431 (2006)
71. Cao, X., Bouzat, S., Kolton, A.B., Rosso, A.: Localization of soft modes at the depinning transition. Phys. Rev. E **97**, 022118 (2018)

72. Wegner, F.: Disordered system with $n$ orbitals per site: $n = \infty$ limit. Phys. Rev. B **19**, 783-792 (1979)

73. Peled, R., Schenker, J., Shamis, M., Sodin, S.: On the Wegner orbital model. Int. Math. Res. Not. **2019**(4), 1030–1058 (2019)

74. Khorunzhy, A.M., Pastur, L.A.: Limits of infinite interaction radius, dimensionality and the number of components for random operators with off-diagonal randomness. Commun. Math. Phys. **153**(3), 605–646 (1993)

75. Pastur, L.A.: Spectra of random self adjoint operators. Russ. Math. Surv. **28**, 1–67 (1973)

76. Lifshits, I.M., Gredeskul, S.A., Pastur, L.A.: Introduction to the theory of disordered systems. Wiley, Hoboken (1988)

77. Fyodorov, Y.V., Sommers, H.-J.: Classical particle in a box with random potential: exploiting rotational symmetry of replicated Hamiltonian. Nucl. Phys. B **764**(3), 128–167 (2007)

78. Carpentier, D., Le Doussal, P.: Glass transition of a particle in a random potential, front selection in nonlinear renormalization group, and entropic phenomena in Liouville and sinh-Gordon models. Phys. Rev. E **63**, 026110 (2001)

79. Fyodorov, Y.V., Bouchaud, J.-P.: Freezing and extreme-value statistics in a random energy model with logarithmically correlated potential. J. Phys. A **41**(37), 372001 (2008)

80. Fyodorov, Y.V., Le Doussal, P., Rosso, A.: Statistical mechanics of logarithmic REM: duality, freezing and extreme value statistics of $1/f$ noises generated by Gaussian free fields. J. Stat. Mech. **2009**(10), P10005 (2009)

81. Fyodorov, Y.V., Keating, J.P.: Freezing transitions and extreme values: random-matrix theory, $\zeta(\frac{1}{2} + it)$ and disordered landscapes. Phil. Trans. R. Soc. A **372**(2007), 20120503 (2014)

82. Subag, E.: Following the ground state of full-RSB spherical spin glasses. arXiv:1812.04588

83. Fyodorov, Y.V., Le Doussal, P.: Manifolds in high dimensional random landscape: complexity of stationary points and depinning. Phys. Rev. E **101**, 020101(R) (2020)

84. Fyodorov, Y.V., Mirlin, A.D.: Scaling properties of localization in random band matrices: a $\sigma$-model approach. Phys. Rev. Lett. **67**, 2405–2409 (1991)

85. Shcherbina, M., Shcherbina, T.: Universality for 1d random band matrices: sigma-model approximation. J. Stat. Phys. **172**(2), 627–664 (2018)

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.