One step replica symmetry breaking and extreme order statistics of logarithmic REMs

Building upon the one-step replica symmetry breaking formalism, duly understood and ramified, we show that the sequence of ordered extreme values of a general class of Euclidean-space logarithmically correlated random energy models (logREMs) behave in the thermodynamic limit as a randomly shifted decorated exponential Poisson point process. The distribution of the random shift is determined solely by the large-distance ("infra-red", IR) limit of the model, and is equal to the free energy distribution at the critical temperature up to a translation. the decoration process is determined solely by the small-distance ("ultraviolet", UV) limit, in terms of the biased minimal process. Our approach provides connections of the replica framework to results in the probability literature and sheds further light on the freezing/duality conjecture which was the source of many previous results for log-REMs. In this way we derive the general and explicit formulae for the joint probability density of depths of the first and second minima (as well its higher-order generalizations) in terms of model-specific contributions from UV as well as IR limits. In particular, we show that the second min statistics is largely independent of details of UV data, whose influence is seen only through the mean value of the gap. For a given log-correlated field this parameter can be evaluated numerically, and we provide several numerical tests of our theory using the circular model of $1/f$-noise.


I. INTRODUCTION
Extreme value statistics of logarithmically correlated random energy models (logREMs) is a subject of active recent research by both physicists and mathematicians. LogREMs (and their ancestors, uncorrelated REM and simplycorrelated GREM) originated in physics [1][2][3][4] as simplified models of mean field spin glasses, and then were singled out as a particularly interesting case of the general problem, the thermodynamics of a particle in random potential [5][6][7]. When approached with the replica trick of disordered statistical mechanics, these models are known to belong to the class that can be solved by one-step replica symmetry breaking (1RSB) scheme. Such solution is certainly much simpler than the full replica symmetry breaking scheme deemed necessary for description of the low-temperature phase of other classes of the mean field spin-glass models [8,9], see [10] for the modern rigorous approach, as well as of general manifolds in disordered landscapes [11,12]. Yet one can argue that among all models of 1RSB class, logREMs are a special marginal case precisely at the boundary with the full RSB models [6,7], and that at the mean-field level they share such features of the full RSB as presence of massless replicon modes in their fluctuation spectrum in the whole low-temperature phase, see e.g. Appendix D of [13].
Another line of research initiated in [4] heavily relies on relations between LogREMs and traveling wave equations of Fisher-Kolmogorov-Petrovsky-Piscounov (FKPP) [14,15] type. Such a relation is exact for a particular instance of LogREM, the branching Brownian motion (bBm), and the deep mathematical analysis in the FKPP framework due to Bramson [16] provided a paradigm for understanding generic logREMs. In physics literature these ideas were promoted to greater generality by a renormalization group argument in [5]. More recently, in the mathematical literature adopting and generalizing Bramson's methods and ideas to other log-correlated processes and fields, most notably to the two-dimensional Gaussian Free Field (2DGFF), helped to develop understanding of extreme and high values of those objects to a great depth , see [17][18][19] and references later in this section. Besides providing mathematic tractability, the logREM-FKPP interplay looked upon from the physical viewpoint [4,5] revealed the fundamental freezing phenomenon manifesting itself via the temperature-independence of extensive free energy in the whole low-temperature, 1RSB phase of logREMs. Basic thermodynamics implies then that the entropy vanishes, i.e., the Boltzmann-Gibbs measure is dominated by a few extremely low energies. The freezing phenomenon is the central feature of logREMs, with consequences in other areas of physics and mathematics. To that end one may mention the multi-fractal properties of electronic wave-functions in disordered gauge field [20,21], as well as understanding of high values of multifractal patterns of more general nature [22,23]. Most recently freezing was employed to conjecture the statistics of high values of characteristic polynomials of random matrices and eventually high values of the Riemann zeta-function in intervals along the critical line [24][25][26] which generated a new stream of activity [27][28][29][30]. As such, the whole subject of freezing and logREM extrema was under an intensive research in recent years and was addressed from several directions and viewpoints. Not attempting a comprehensive review, we very briefly survey this body of work below: -Building upon the picture of the freezing phenomenon developed in [5] for Euclidean-space logREMs (also known as logREMs with statistical translational invariance), the paper [7] proposed for the first time exact free-energy distribution for a particular regularized version of the 2DGFF sampled along a circle, which is also known as (log-)circular model of 1/f -noise. The treatment was based on assuming the validity of the freezing scenario conjecturing that in the thermodynamic limit the whole distribution of the Gumbel-convoluted free-energy F − G/β, with F being the free energy at temperature 1/β, and G being a standard Gumbel variable independent of F , is freezing, that is becomes temperature-independent below the transition temperature, modulo a translation. Further work [31] observed the co-existence of freezing with the duality-invariance property, and extended the above freezing scenario to the freezing-duality conjecture (FDC), which has been checked numerically for several observables in various logREMs [31][32][33][34][35]. Although a relation between FDC and 1RSB was hinted briefly in [32], up to now very little has been worked out explicitly in this direction. Intriguingly, recent mathematical work [36,37] provided a proof of the freezing of the distribution for F − G/β without any reference to duality, though physical approach points towards an intimate relation between the two, at least for all the Gaussian models that are exactly solved (these models are sometimes referred to as integrable, although quantum/classical integrable structures behind them are not clearly identified yet). In fact, a lot of general and rigorous results concerning the extreme values of 2d GFF and the associated Boltzmann-Gibbs-Liouville measure have been obtained in recent years in mathematical literature. Note that the 2d GFF is a random generalized function [38] so that studying its value distribution necessarily involves a regularization scheme. Among most popular rigorous schemes are the 'multiplicative chaos' framework [39], see applications to LogREM-related models in [40,41], and various variants of the discretization/lattice models [19,[42][43][44][45][46][47]. Although providing so far no explicit distribution of free-energy or the global minimum (see however [48]), these studies corroborated the freezing picture of FDC.
-In a closely related development, an understanding of the full process of minimal values of the branching Brownian motion (bBm) was initiated by [49] and elaborated in full rigour in [50][51][52][53]. The bBm is the prototype of hierarchical logREMs. Such models are exactly described by FKPP equations, which describes only approximately the Euclidean-space ones discussed in the precedent paragraph. It is now known that in the thermodynamic limit, the minima process tends to a randomly shifted, decorated Gumbel Poisson point process (SDPPP) [37]. With respect to the Gumbel Poisson point process which describes the minima of the uncorrelated REM, the novel ingredient of SDPPP, i.e., the decoration process, describes the internal structure of blocks of extreme values which share a near ancestor, and thus are highly correlated. Such a picture is expected to apply also in the context of non-hierarchical logREMs, in particular those related to 2d GFF and 1/f -noise, in which the blocks are those of extreme values given by nearby points. However, until recently [54], this expected picture has not been established in worked out in sufficient quantitative details which can be compared to numerical simulations of particular instances of logREMs. Also, despite the insightful discussions made in [37], we feel that the relation between SDPPP, freezing and 1RSB needs to be spelled out in a more explicit and comprehensive manner in the context of GFF-type logREMs. In particular, we feel necessary to start bridging the remaining gaps between physicists' and mathematicians' description of this common object of interest.
The goal of our present work is to provide such a detailed account at the physical level of rigour by employing the standard and powerful, albeit heuristic framework and language of the replica approach to a general class of logREMs. In doing so we will assume a few conditions, similar to [19], which allow us to define what physicists usually call the ultraviolet (UV, short-distance) and infrared (IR, long-distance) data of a given logREM model in the thermodynamic limit. Under these assumptions, using 1RSB scheme, we are able to show that the ordered sequence of minimal values forms a SDPPP. More precisely, the random shift is equal to the free energy at critical temperature (modulo a translation), which depends only on the IR data. The free energy distribution at the critical temperature is shown to determine also the distribution in the whole low temperature phase (and thus the minimum) by the freezing scenario, of which we spell out the relation with the duality invariance property in the 1RSB framework. The decorating process depends, on the other hand, solely on the UV data, and is equal to the process of biased minima. In this way we derive the general and explicit formulae for the joint probability density of depths of the first and second minima (as well its higher-order generalizations) in terms of model-specific contributions from UV as well as IR limits. In particular, we show that the second min statistics is largely independent of details of UV data, whose influence is seen only through the mean value of the gap in a universal way. The paper is organised as follows. In Sect. II, after defining the class of logREMs that we will consider and specifying their IR and UV data, we will state our results restricted to first and second minima for clarity, and then eventually describe the whole sequence of minimal values. In Sect. III we test some of our analytical predictions numerically using the circular 1/f -noise model. Sect. IV we provide basic information on the 1RSB framework. After describing in general the 1RSB approach for logREM in Sect. IV A, we apply the method to retrieve the known freezing result of the free energy/minimum of the potential (Sect. IV B), and elaborate on the 1RSB-FDC relation in Sect. IV B 2. After this warm up, we spell out the derivation for the second minima (eq. (24)) in Sect. IV C. For a reader interested only in getting flavour of the method rather than in full technical detail Sections IV B and IV C (without sub-subsections) constitute a summary of the derivation of (24).
Sect. V generalizes the above results to the whole sequence of minima values, where the structure of the decorated Poisson point process becomes evident. Sect. V A gives the derivation of (41), which generalises (24), and discusses its consequences, including the joint distribution of all order statistics, which coincides with SDPPP (this will be shown in Appendix C). Sect. D is devoted to the calculation of some marginal distributions (e.g., that of k-th minimum) of SDPPP. Sect. VI concentrates on the UV sector; we shall first study a variant of the uncorrelated REM that mimics the UV structure of logREMs, and confirms the 1RSB treatment of the UV sector by a replica-free calculation; this will be followed by a discussion on the nature of the decoration process.
A few appendices are devoted to more technical topics. Appendix A gives a full replica symmetry breaking analysis of logREM and its multi-scale generalisations. Appendix B computes by 1RSB the Edwards-Anderson order parameter (defined in [35]) for a general class of spherical logREMs. Appendix C recalls the definition of SDPPP and calculates its joint minima value distribution, of which Appendix D computes some marginal distributions. Motivated by recent work [54], Appendix E takes the minima positions into the picture. Finally, Appendix F provides a derivation of a technical result on uncorrelated REM which is used in Sect. VI A.

II. SYNOPSIS
In this section, after defining the UV and IR data of logREMs in Sect. II A, we give the main results in Sect. II B restricting ourselves initially to the value distribution of the first and second minima which will be eventually the scope of direct numerical checks (in Sect. III). That particular case captures already all the most essential structures of the problem. Finally, in Sect. II C we provide our results for the general case as well.

A. Log-REMs and their UV and IR limits
In this paper we are going to consider a general log-REM associated with a real random Gaussian potential V j,M , with j = 1, ..M, indexing M sites of a lattice (which we assume to be either d-dimensional or a hierarchical tree). The covariance V i,M V j,M c will be assumed to decay logarithmically with the distance between sites, e.g. in the case of d-dimensional lattice we may informally assume V (x)V (y) c ∼ −2d ln (|x − y|/M ), a more detailed description following later on. The property of logarithmic decay can be conveniently characterized quantitatively as follows where |A| is the size of the set A. Throughout this work, [. . . ] means averaging over disorder, i.e., over the (V j,M ) M j=1 in the above case, unless otherwise stated. In particular, we require the (log)REM variance a key assumption to be supplemented below in (12). Note that this definition applies to both hierarchical logREMs and logREMs in Euclidean space, i.e. with Statistical Translational Invariance. The partition function of the logREM is defined as usual where β is the inverse temperature. The definition (1) has the advantage of fixing the critical temperature The general criterion (1) specifies only the covariance in the scaling regime 0 < q < 1. Although this is enough to guarantee the super-universal properties of logREM predicted in [5], to fully define the model and calculate the distribution of minima values one needs to specify the model behaviour at the system-size scale and the lattice-spacing scale, in terms of IR and UV limit data, respectively, which we define below. An illustration is given in Fig. 1.

IR data
The IR data consist first of a d-dimensional manifold X discretized by a family of point meshes ξ i,M ∈ X that tend (as M → ∞) to a point density measure: The other part of the IR data is a continuous mean value (sometimes called the background charge) c(ξ) and a continuous covariance kernel C(ξ, η), ξ, η ∈ X satisfying: To ensure matching with (2), the kernel C should be logarithmically diverging when ξ → η: more precisely, if X is d-dimensional, The IR data allow to define the continuous replica partition function The reason for which we write b instead of β in (10) will be clear later on. The replicated partition function is usually given a Coulomb-gas integral, except that the (identical) charges attract each other. This causes the integral to diverge, generically for nb 2 ≥ 1. Nevertheless we shall assume that by analytical continuation, (10) can be defined for generic n, b ∈ C (this has been carried out for several concrete models starting from [7,31], see [55][56][57][58][59][60] for the relevant rigorous mathematical developments). Note that here the [. . . ] is rather a metaphor since we do not define the continuum random partition function Z at the first place; heuristically it could be defined as Then (10) can be "derived" by Wick theorem.
Here v(ξ) can be seen as a regularized Gaussian Free Field with zero mean. It does not make probabilistic sense directly, since e.g., its variance is negative, however a rigorous construction of the integrand of the RHS is possible in the multiplicative chaos framework developed originally in [39], see [61] for recent developments. We will return to the meaning of Z b later in E.
When approaching the critical temperature b → 1 − , it is now known [40] that the well-behaved object is not the partition function Z b but its derivative (∂ b Z b )| b=1 . This echoes the fact that, in all known exactly solved cases [7,[31][32][33][34][35], the analytically continued continuum integral Z n b has a zero of order n at b = β c = 1. Therefore we define The same object (modulo a constant) was called z n op. cit.. We shall assume that (11) exists for generic n ∈ C.

UV data
The UV data describe the limit of covariance between points i and j separated by a fixed distance in units of lattice spacing in the M → ∞ limit. From the IR point of view, ξ i,M and ξ j,M tend to the same point ξ ∈ X in this limit, and V i,M V j,M ≈ C(ξ, ξ) has a divergence. However, if we zoom in to the lattice spacing scale, the points ξ i,M close to ξ will form a d-dimensional hyper-cubic lattice where the grid side a ∝ (M exp(−l(ξ))) −1/d because of (5). Thus we expect from (9) that C(ξ i,M , ξ j,M ) ≈ −2d ln |ξ i,M − ξ j,M | ≈ 2(ln M − l(ξ)); note that this agrees with, and supplements, eq. (2). Now the UV data are defined as the limit of all further corrections, which we express as a symmetric matrix (f ij ): We assume that the limit exists, does not depend on ξ (homogeneity of UV data), and locally translation invariant in the sense that f ij = f (x i − x j ). This implies in particular that diagonal elements f ii are all equal. The notion of UV limit data will play a crucial rôle in the correct calculation of the higher order statistics because different minima can be at the distance of order of lattice spacing from each other, and probe the covariance matrix at that scale. Note that only the IR data are needed for calculating the minimum value distribution. Obtained by different limit procedures, the UV and IR data are independent of each other, in the sense that one can modify the IR data for a given model while keeping its UV data intact, and vice versa. Let us the latter point by reviewing some concrete 1-d models.

Example: Circular model
The circular model of 1/f -noise [7,22] is a Gaussian log-REM defined by the following mean and covariance matrix: Its IR data are the following, as one can easily check by (8): The points ξ j,M provide thus uniform meshes of the unit circle, whereas C is the Green function of the planar 2d Gaussian free field. The associated continuum Coulomb integral (10) is the Dyson integral [62], which has a closed form expression that can be analytically continued: In particular, (11) is satisfied. Its UV data can be calculated from (12) using l(ξ) = ln(2π) from (14) where Ci(x) = − ∞ x cos t/tdt is the cosine integral. As |i − j| → ∞, f ij ≈ −2 ln |i − j| grows logarithmically, which is a general feature of UV asymptotics for logREMs.
A few variants of the circular model have been studied. One example is a generalization based on 2D field in a disc with Dirichlet boundary conditions [31], which modifies the covariance of (13) by replacing λ k = 1/k with λ k = (1 − q k )/k, where q ∈ (0, 1). One can check that the UV data (16) remain unchanged, while the continuous covariance matrix of IR data is modified to C(ξ, η) = −2 ln |ξ − η|+2 ln |ξ − qη| whereas l(ξ) and c(ξ) are not changed. For these models the analytical continuation of the continuum integral (10) is known [34] only as an infinite sum, not as a closed-form formula similar to (15).
One can also change solely the UV data. A such example is the "long range" circular model [31], defined as , different from (16). However, the IR limit remains the same, mainly because lim M→∞ λ k = k −1 as in (13).

Example: Interval model
An example of less trivial l(ξ) and c(ξ) is the interval model [31,33]. Its general discrete-version definition depends on four parameters c 0 , c 1 , d 0 , d 1 1 . One starts by defining the non-uniform mesh points whose density follows a Beta distribution: in terms of which the mean and covariance of V j,M are given by Here f 0 is a constant, which can be fixed arbitrarily, as long as the covariance matrix above is semi positive definite for all M . The IR data of the model do not depend on f 0 , and are given by c(ξ) = c 0 ln ξ + c 1 ln(1 − ξ) and l(ξ) = −d 0 ln ξ − d 1 ln(1 − ξ), and its continuum limit representation is given by the Selberg integral [63]: The continuation of its exact solution (due to Selberg) to complex values of n is discussed in [31] and [55][56][57][58][59][60]. In particular, one knows that it has a Γ(1 − b 2 ) −n factor, like (15), so that (11) makes sense. On the other hand, the UV data depends only on f 0 in (18), and is given by f ii = f 0 and f ij = −2 ln |i − j| for i = j.
B. First and second minima

Observables
From now on, we will omit the system-size label, i.e., ...,M will be dropped in the subscripts to ease the notations, so V i = V i,M denotes the values of some logREM with well-defined IR and UV data. We consider the following observables related to the first and second minima: The first observable is the generating function of the partition function (i.e., the Laplace transform of its probability density) commonly used in the disordered system context; see (53) and (52). In particular, it decreases from 1 to 0 as y goes from −∞ to ∞. H β can be also interpreted in a similar way, see (75). They have the following zero-temperature limits: where here and below we denote ∂ y1,...,yp def = ∂ y1 . . . ∂ yp . Namely, 1 − G ∞ is the cumulative density function of the minimum value, and derivatives of H ∞ (y 0 , y) provide the joint probability distribution of the first and second minima values. Note that if we had not summed over j in (21), the zero temperature limit (23) would contain in addition the minimum position.

Main result
The main result is encoded in the equation which holds in the M → ∞ limit for all β : where K β (∆) is the UV correction factor definition of which will be given below. In spite of the product form, we do not have an obvious interpretation of (24) in terms of statistical independence; in particular, the minimum and the first gap V min,1 − V min are not uncorrelated, see (32) below. Nevertheless, the two terms of the RHS have contrasted nature, as we detail below: -The shape of G β (y) (i.e. modulo a translation) depends only on the IR data (10) and manifests freezing: In (26), F (M, β, f ij ) contains the M dependence of the free energy, and a M -independent, UV-data-dependent correction. The latter can be exactly calculated in the β < 1 phase (recall that f ii is independent of i), see the beginning of Sect. IV B 2, eq. (66). In the other phases this correction remains an unknown constant (a 0 , a 1 above); this is related to unsolved issues when Z in (10) is replaced by the renormalized value Z in (11). As a consequence of (25), the distribution of the minimum (modulo a UV dependent translation) is determined by the critical partition function Z −t of eq. (11).
-The UV correction factor K β (∆) depends only on the UV data f ij . In particular, distribution of the gap between the first and second minima only depends on UV data (f ij ) defined in (12). To describe how the latter determine K β (∆), we define the local logREM, which is a sequence of Gaussian variables u i with zero mean and the following covariance: Here the over-line [. . . ] means averaging over the local disorder (u i ) N i=1 , with C N being a large positive number chosen to make the covariance matrix positive-definite (whose precise value turns out to be irrelevant, see below). Then K β (∆) is given in terms of the partition functions of the local logREM: The numerical convergence as N → ∞ is expected to be much faster than in the M → ∞ limit in the originally defined logREM, eq (1); the reason is that the observable K β (∆) captures essentially microscopic information, even when N → ∞. We will discuss this in more detail and generality in Sect. VI B.

High temperature phase
In the high-temperature β < 1 phase, where m = 1, (28) and (29) become easy to evaluate: Therefore (24) simplifies to The equation is (super) universal since K β<1 , (30), depends on neither UV nor IR data. At the same time, in the low-temperature β > 1 phase K β will have in general a non-trivial temperature dependence, in contrast to G β>1 , which manifests freezing there. The non-freezing of K β (and its higher order statistics extensions, see (98)), is the essentially new feature that prevents us from extending naïvely freezing to conclude "H β>1 = H β=1 ", which as we find is false in general. Nevertheless, the naïve freezing does hold for a modified version of H β (see (82)), whose zero temperature limit gives the statistics of the second minimum whose position is far away from that of the global minimum, see (50) for precise definition.

Zero-temperature limit: distribution of the first and second minima values
In the zero temperature β → ∞ limit, (24) and (23) give the joint distribution of V min and V min,1 : where K ∞ can be obtained by taking the β → ∞ limit of K β defined in (28). In that limit, the local partition functions (29) are dominated by u min and u min,1 , the first and second minima of the local logREM, respectively. More precisely, we have z β → e −βumin for any δ ≥ 0 and can check that w β (j, δ) → e −βumin if u j = u min , and w β (j, δ) → e −β min(umin+δ,umin,1) otherwise. Now inserting these into (28) with δ = ∆ and δ = +∞ and observing that only the case where u j = u min yields a non-vanishing contribution in (28) we obtain K ∞ (∆) = exp(− min(u min + ∆, u min,1 )) − exp(−u min,1 ) exp(−u min ) .
Differentiating, we get We would like to emphasize that the factors exp(−u min ) in the numerator and denominator do not cancel each other unless assuming independence of u min and the local gap u min,1 − u min , so K ∆ is not a function of only the local gap u min,1 − u min . Nevertheless, (34) implies that −K ′ (∆) decreases from −K ′ ∞ (0) = 1 to K ′ ∞ (∆ → ∞) = 0. This allows us to define a new pair of the biased extrema v min < v min,1 for which the following relation holds: The biased extrema are defined up to an overall translation, i.e. only the gaps are meaningful observables and one may as well fix v min = 0. They are biased in the sense that events in which u min is more negative have dominating (with weight exp(−u min )) contribution to the statistics of the gap u min,1 − u min .
Combining (35) and (32) we see that the first gap g 1 = V min,1 − V min satisfies For the second minimum V min,1 we have In other words, possible distributions of V min,1 lie in a one-parameter family determined solely by the IR data, with the parameter being the mean value of the gap depending only on the UV data. The above results are only a fragment of a larger picture, i.e. that of the randomly shifted decorated Gumbel Poisson point process, with the decoration process being the biased local minima process (v min , v min,1 , . . . ). The best way to see how the latter emerges is to generalize the results to the full minima value process, which we are going to describe now.

C. Higher order statistics
Here we summarize the main results on the higher order statistics, referring to Sect. V, Appendix D and E for the detailed derivation and discussions, as well as for additional results. They are all obtained from the 1RSB calculation of the following observable that generalizes H β (y 0 , y) of (21): Here y 0 < y 1 < · · · < y k−1 < y and j 1 , . . . j k−1 are all distinct; a position-value couple (j s , y s ) will be called a marker.
Some explanations are in order. First, the equality of (41) holds in the M → ∞ limit. On the LHS, the notation * j0...j k−1 means that the sum is over distinct positions j 0 , . . . , j k−1 . On the RHS, the sum P(k) is over all the partitions of the set of k first minima {y 0 , . . . , y k−1 } = p q=1 {y q,0 < · · · < y q,kq−1 }; here 1 ≤ p ≤ k is the number of parts and k 1 , . . . , k p , with p q=1 k q = k, are the sizes of the parts (see Sect. V A for detailed introduction to partitions). The term corresponding to a partition of p parts contains a product of p UV correction factors; they are an infinite family of functions {K β (∆ 1 , . . . , ∆ l )} l=1,2,... , defined independently of k in (41); their explicit definition will be given in (96), here it suffices to mention that they depend only on the UV data f ij as well as on the temperature β, and reduce to (28) when l = 1.
In the β → ∞ limit (41) implies the following joint distribution of the ordered minima V min = V min,0 < V min,1 < V min,2 < . . . of the logREM (i.e., V min,k is the (k + 1)-st minimum): This equation is equivalent to saying that all the minima V min , V min,1 , . . . are generated by a randomly shifted decorated Gumbel Poisson point process (SDPPP), such that the random shift has the same distribution (modulo a translation) as the free energy at critical temperature β = 1, and the decoration process is given by the biased minima v min , v min,1 , . . . . Their statistics are characterized in terms of minima of the local logREM u min , u min,1 , . . . by (44), which is a straightforward generalization of (35). SDPPPs have recently been shown to describe the minima of several models related to Branching Brownian motion [52,53], and are expected to describe minima of general logREMs. For logREMs generated by discrete GFF this was very recently shown rigorously in [54]. Our analysis provides a replicaapproach derivation of this result in great generality, and provides an alternative characterization of the decoration process. Furthermore, the results concerning the first and second minima will be tested numerically to high precision in the next Section III.
In Appendix D, we generalize (39) to the marginal distribution of all higher order minima in terms of the following generating function: Similarly, a generalization of (38) is obtained for the marginal distribution of all gaps g k in terms of the biased minima process See Appendix D 1 for more discussion. Note that our main formulae (40)-(46) are expected to be valid for any SDPPP, so that any of its possible characterizations (e.g. like one used in [37]) should necessarily imply them. Nevertheless, to the best of our knowledge such expressions did not seem to appear in the mathematical literature in an explicit form.

III. NUMERICAL STUDY
In this section, we test our predictions on the circular model defined in (13). Recall that its covariance matrix is cyclic: the Fourier mode k has the eigenvalue λ can be generated efficiently by the fast Fourier transform [31]. Among all Euclidean-space logREMs this model is arguably the most amenable to both analytical and numerical analysis, and is also singled out by its relevance in the random matrix context, see [23,31]. However, despite intensive interest all previous studies focused only on its IR data, in particular on the distribution of the (shifted) global minimum derived originally in [7]: Here K n is the Bessel K-function and the unknown shift y M can be fixed by the average value: The expression (39) of the present paper allows us to predict the following explicit form of the distribution of the second deepest minimum: Let us test this prediction directly. For this we generate S independent samples of the circular model of size M = 2 8 , . . . , 2 23 , with S = 10 7 for M ≤ 2 15 and S ≥ 10 6 for M ≤ 2 20 and S ≥ 10 5 for even larger M . We first observe, in Fig. 2 (Left panel) that the mean value of the gap V min,1 − V min has a strong M dependence which should be taken into account to extract the thermodynamic limit value g 1 . A good numerical fit is provided by a quadratic function of 1/ ln M (the same Ansatz is known [35] to work in other observables of this model). The result is g 1 = 0.70(1), and is fed into the RHS of (48). On the other hand, we collect the full distributions of the shifted second minima V min,1 − V min + 2γ E where V min is the numerically measured mean value of the minimum of the same system size.
The results, with the prediction (48) subtracted, and shown in Fig. 2 (Right). We observe again a slow convergence; nevertheless, the extrapolation to M → ∞ with the same quadratic 1/ ln M Ansatz has an excellent agreement with the prediction (48). To compare, we plotted also the RHS of (48), with two other values g = 0 and g = 1; it is clear that both would give wrong predictions.  Next, we focus on the UV sector of our theory, and study the temperature dependence of the factor K β (∆). We focus on the cases ∆ = 0 for simplicity. To probe it numerically, we shall use the following equation, which we shall derive in Sect. IV C 1: Now the right hand side can be directly measured numerically. We do this for 40 different values of β = 0.05, . . . , ≅ 4.5, each in various system sizes M = 2 10 , . . . , 2 20 in order to extrapolate to thermodynamic limit using a quadratic 1/ ln M Ansatz similarly to the procedure described above. The result is shown in Fig. 3. In the β < 1 phase, we have the analytical prediction K β (0) = 1, (30). In that phase, away from β = 1, the numerics agrees perfectly with the prediction, with invisible finite-size effects. Finite-size effects become more pronounced as β 0.8, and have an intriguing non-monotone behaviour around β = 1. In the β > 1 phase the naïve freezing K β>1 (0) = 1 is clearly ruled out, and so is the "analytical continuation" K β>1 (0) = 1/β, which accounts for the explicit min(1, 1/β) nonanalyticity in (49). The β-dependence of K β (0) in the β > 1 phase is non-trivial and we do not have an educated guess to fit the data. To check consistency we look at the β → ∞ limit; there, (38) predicts that K ∞ (0) = g 1 , the latter being the mean value of the first gap, which we measure numerically in Fig. 2 (Left panel) to be g 1 ≅ 0.7, and indicated in Fig. 3 as the black horizontal line. As expected, we observe that the latter coincides with the β → ∞ asymptotic value of the extrapolated K β (0) data.
In the β > 1 phase, K β (∆) is determined by the UV data of the model (given in (16)) via the local logREM, see eqs. (27) and (28). We now test these predictions in the zero temperature β → ∞ limit, where the gap distribution is related to the biased gap defined in (35), see eq. (37). Thanks to the translation-invariant structure of the covariance, the local logREM (27) of the circular model can also be simulated efficiently using the fast Fourier transform. We do this for N = 2 1 , . . . , 2 12 , with 10 8 independent samples for each size. We measure, using (35), the cumulative distribution function of the biased gap θ(v min,1 − v min − ∆), and extrapolate to the N → ∞ limit using the quadratic 1/ ln N Ansatz. We then compare this with θ(V min,1 − V min − ∆)e ∆ measured independently in the original circular model , and extrapolated to M → ∞ using the same Ansatz in 1/ log M . The result, Fig.4 (Left), shows an excellent agreement between the two limit distributions, validating the prediction (37). We remark that the finite size corrections of the local logREM are much smaller than those of the full model: even N = 8 provides an approximation as decent as M = 256. On the other hand, the measure of the biased minima (35) is sensitive to contributions of rare events, and requires to collect better statistics to be reliably verified.
Despite being unable to deal with non-trivial UV data analytically, we can get rid of them by defining the second minimum in a different way, as   The results (blue curves) are multiplied by exp(∆) (note that an exponential distribution would give a horizontal line at height 1). The first gap of the biased process is measured by (35), in local logREMs of sizes 2 ≤ N ≤ 2 12 (red circles), and extrapolated using the same Ansatz to N → ∞ (black circles). The good agreement between the latter curve (without multiplying by exp ∆) and the M → ∞ curve confirms the prediction (37). Right: the far-gap distribution g far 1 , defined in (50), with N = √ M , compares well to the exponential, (51) .
where j min the position of the minimum, whereas both N → ∞ and M → ∞ keeping N ≪ M in the thermodynamic limit; numerically, N = √ M is used. The idea is that when looking for the second minimum located far from the global minimum the UV data trivialize and one retrieves the exponential gap distribution in the thermodynamic limit: which is the case for the random energy model with uncorrelated potential. This prediction is well verified by the numerical data, which are shown in Fig. 4 (Right). A derivation of this result within the 1RSB framework will be outlined in Sect. IV C 2.

IV. FIRST AND SECOND MINIMA
In this section we will derive the main results concerning the first and second minima value distribution in Sect. II B. First we are going to give an outline of the essential ingredients of the 1RSB approach to logREM models, with calculation of the minimum value distribution presented as a good pedagogical example. Computations of the higher order statistics involve natural, albeit non-trivial, ramifications of that procedure.
As a starting point for our analysis we consider a finite temperature generating function (20) related to the global minimum value distribution: where Z = M j=1 exp(−βV j,M ) is the standard partition function of the logREM as defined in (3). In the replica approach, it is also convenient to present (52) as a formal power series: A direct consequence of the representation (53) which will be used later on is that differentiation on the LHS corresponds to multiplication on the RHS: A. Overview of the 1RSB Ansatz Z n is a sum over positions of n replicas, which "interact" via disorder. Here, the interaction is pair-wise as a result of the Wick theorem: This is the prototype of replica averages to which replica symmetry breaking (RSB) Ansätze are applied. For logREMs (as well as uncorrelated REMs), it is known that a one-step RSB (1RSB) Ansatz is correct. For hierarchical logREMs, this fact was first realized [4] using KPP equation. For Euclidean-space logREMs, 1RSB was argued to hold in [5] (Sect. III.E.2), by building upon a functional renormalisation group analysis, from which KPP-type equations emerge again, as flow equations. In Appendix A, we provide a full replica symmetry breaking analysis of logREMs, which show 1RSB results from the general definition (1). In the main text, we will take for granted the 1RSB Ansatz as the starting point of our derivations. The Ansatz claims the following (such a description of 1RSB, which is a ramification of the standard one, first appeared in [32]): In the thermodynamic limit the replica sums such as (55) are dominated by configurations such as depicted in Fig. 5. More precisely, the n replicas form n/m groups of equal size m. Within each individual group of size m different replicas should be thought of as being within the mutual distances of the order of lattice spacing, which is the definition of the UV scale; therefore, their interaction should be calculated using (12). On the other hand, distinct groups are separated typically by IR (system-size) scale, so one can treat the sum over group positions as continuum integral, and calculate inter-group interactions using IR data. Now, the group size m is a non-analytical function of β: The unique transition at β = 1 (4) corresponds to the non-analyticity of (56). When β > 1 one has m = 1/β < 1 which is the standard, though odd, feature of the replica limit n → 0 which forces one to consider all integer parameters as promoted to variables which can take in general any real values (strictly speaking, continuation to the complex plane is necessary, e.g., when performing inverse Laplace transform). We shall see how the above applies to the replica sums Z n in the next subsection. We argue that the same method applies for partition sums used to address higher order statistics. Indeed, the object in question here is a discrete partition function of a particle in a potential (U j ), obtained from that of a (log)REM by shifting a finite number of values: Here δ is the Kronecker delta, and (j 1 ,

B. 1RSB for minimum
The general idea of the 1RSB calculation of Z n as described in Sect. IV A has been briefly outlined in [32], and we describe it in detail below for the sake of clarity and completeness. For this, we divide the system into blocks of N sites, where N is any intermediate scale such that N ≪ M , whereas considering that both N → ∞ and M → ∞ in the thermodynamic limit. Each block will be labelled by their position in C (because the ξ j,M ∈ C of a same block will converge to the same point); positions inside a block labelled 1, . . . , N ; in other words, we have a one-to-one correspondence j(ξ, i) between the global address and the hierarchical one. Now, n replicas form n/m groups of size m; when n is continued to complex values we will have m = max(1, 1/β) by eq. (56). Let us index the groups by g = 1, . . . , n/m and let g(a), a = 1, . . . , n be the group that replica a belongs to. Due to the permutation symmetry in replica space we need only consider, for instance, a particular g(a) = ⌈a/m⌉ (⌈x⌉ standing for the smallest integer larger or equal to x) and multiply by the combinatorial factor counting the number of different groupings 2 Each group is confined in a block with coordinates concentrated at UV distances around a macroscopic position ξ g . Different groups are separated by a distance of system-size order so that the sum over block-positions can be approximated by a continuous integral: Compared with (6), we have M/N factors on the RHS because the LHS sums over block positions and there are M/N blocks. The continuum limits also apply to mean values V j → c(ξ), as well as to inter-group Wick contractions . On the other hand, the intra-group Wick contractions (depicted in Fig. 5 by violet lines in the zoomed-in version) should be treated using the UV limit data (12). We will organize them into intra-group interactions, see (61) below; note however that the factors resulting from the shifts by 2l(ξ) in the LHS of (12) will be systematically combined with other instances of l(ξ) appearing in (59). Carrying out the above treatment, we reduce (55) in the M → ∞ limit to the following expression: Here E(ξ) is the intra-group interaction of the group at block ξ; we note that the bar in z m β denotes averaged moments of the local-logREM partition function z β defined in (29). Note that to alleviate formulae like (61) and below we will tacitly assume that the limit N → ∞ is taken along with that of M → ∞ in the order explained above. The expression (60) can be rearranged in terms of the continuous partition sum defined in (10) with parameters b → mβ and n → n/m Here n should be understood as continued to a complex-valued variable, so that (62) yields the Laplace transform: where Eq. (64)  and a shift S [UV] which depends on β and the UV data. In the high temperature phase, β < 1, by (56), m = 1 and b = β, so C n,m = 1 by (58) and z m β = N exp(β 2 (C N + f ii )/2) by (29) and (27). Plugging this into (64), we have which is equivalent to the β < 1 case of (25) and (26) via the Laplace transform (52). The β > 1 case is more involved and more interesting, because it is directly related to the freezing scenario; we will discuss it in Sect. IV B 1. Sect. IV B 2 compares the current approach to the freezing and duality conjecture. The goal of those sections is to give some simple relations and arguments which, to our knowledge, have not been clearly spelt out, and which clarify the respective rôle of 1RSB and the duality invariance property in the freezing phenomenon.

Freezing scenario
In this subsection, we relate the result (62) to the freezing scenario (25) in the β > 1 phase. Our arguments show quite generally how the 1RSB leads to the freezing scenario.
Throughout this subsection we shall consider the low temperature phase β > 1. In that case (56) gives m = 1/β, and the continuum integral factor in (62), Z nβ mβ is a Coulomb gas integral (10) over nβ = −t charges (see (63)), and with the renormalised temperature b = mβ = 1. Yet we recall from (11) that Z b is ill-behaved at b = 1 and we shall (11), giving In the RHS, the factor Z −t accumulates all the IR-data dependent contributions to the free energy. The factor ǫ −t = exp(t ln(1/ǫ)) corresponds to a shift of − ln ǫ → +∞ in the free energy. For finite system sizes M this divergence is regularized and as argued in [23] should give rise to the ln ln M corrections in (26).
Let us now discuss the problem of analytical continuation of the combinatorial factor C n,m defined a priori only when n, m and n/m are all natural numbers. To choose the appropriate continuation we use the reflection formula for Γ-functions and, denoting s(x) = sin πx, transform (58) as follows: In the last equality we discarded the fraction ms( n m ) s(n) . This is justifiable because for a fixed m, lim n→n0 ms( n m ) s(n) = 1, provided n 0 /m is an even integer or m is an odd integer. We claim that the last form of (68) is the correct continuation of the combinatorial factor. To see this, we plug (68) back to (62) (and substitute m = 1/β for the β > 1 phase), and separate terms of deterministic shift and fluctuation types: see also eqs (67) and (63)  it has no poles) if Z t is; had we chosen the naïve continuation of (58), this would no longer be the case since we would have instead " exp(tF ) = Γ(1−tβ −1 ) Γ(1−t) Z −t e tc ′′ . Now, it is easy to see that the prediction for the fluctuating part of the free energy F encoded in (70), which was obtained solely from arguments based on 1RSB and analytical continuations, is equivalent to the freezing scenario (25): Indeed, the RHS of (25) (for β > 1) and (70) coincide. By (52), the LHS side of (70) equals e tβ −1 ln Z0 R −∂ y G β (y)e ty dy, which is just (25), provided we make the identification of F (M, β > 1, f ij ) in (26) with the unknown and formally divergent constant in (69): Unfortunately, the present level of development of the replica formalism is insufficient for determining this object to O(1) precision in the β > 1 phase, because the relation between the formal divergence in (69) and the ln ln M corrections of (26) is not yet well understood with full quantitative control. More generally, size-dependent corrections of logREM are difficult to get access to analytically, especially using replicas; however, for simple versions of REM there is a recent progress using alternative approach [64].

Freezing-duality conjecture
This subsection discusses the freezing-duality conjecture (FDC), from a 1RSB viewpoint. We will show that duality invariance extends the optimization of the 1RSB parameter m (which is the mechanism behind the formula (56)) to the level of fluctuations, and demonstrate that predictions of FDC so far are fully compatible with 1RSB. Our discussion will be quite general; for a non-trivial concrete example, we reproduce and substantially generalize the calculation of the Edwards-Anderson's order parameter of the circular model [35] using 1RSB in Appendix B.
Let us first recall the essence of FDC. To motivate it, consider the extensive free energy of logREMs, given by (26). In the β < 1 phase, the free energy density f (β) def = lim M→∞ F (M, β, f ij )/ ln M = −(β + β −1 ), which shows that as a function of β, f (β) depends on β only through the combination β + β −1 . This is usually called the duality invariance property; the name refers to the β → β −1 duality transform, but we stress that it concerns only the high temperature phase. In the β > 1 phase, f (β) = −2 as a manifestation of freezing. Freezing can be seen as a consequence of duality, as the latter implies β = 1 is a stationary point of f : ∂ β f | β=1 = 0; then, by convexity and monotonicity, f (β > 1) can only be a constant. FDC promotes this idea further by boldly suggesting that any duality invariant thermodynamic observable must freeze beyond the self-dual point β = 1. It has been successfully applied to calculate free energy distribution [31], moments of minimum position distribution [32,33], and more involved quantities [35].
For the sake of concreteness, let us first restrict our discussion to the distribution of the free energy. It has been observed that for all known exactly solved Gaussian logREMs the continuum partition function (10) enjoys the following duality invariance 4 : where A(Q, t) is a model-dependent function of two arguments, while the factor Γ(1 + t/β) is model independent, and r 1 is a function of β only; in fact, in all known cases, we have r 1 = β −1 ln Γ(1 − β 2 ) + c where c is a numerical constant. For example, for the circular model (15), we have A = Γ(1 + tβ)Γ(1 + t/β) and c = β ln(2π) (clearly, A is invariant by β → 1/β, hence via power series expansion is only a function of β + 1/β and t. To discuss the rôle of the equation (72) in freezing, we reconsider (64) (with the continuation (68) for C n,m ) in the β > 1 phase, but regarding m as a free parameter (i.e., forgetting for the moment (56)): where In the last equality we have used (72) with β → b. It follows that, when β > 1, the solution m = β −1 ⇔ b = 1 of (56) is also a stationary point for the extensive (proportional to ln M ) part of the first moment, as well as for all higher cumulants of F : The equation (74) shows that the duality invariance property (72) ensures that the choice m = β −1 of (56) satisfies an infinite series of selection criteria 5 . The first of them concerns the optimization of the extensive free energy and is familiar in the disordered system theories, while the others extend the RSB optimization principle to fluctuations. As a consequence, the renormalised temperature b = mβ = 1 freezes in the β > 1 phase, and so does the dual invariant quantity (73): this is an instance of FDC. In hindsight, from a 1RSB viewpoint, FDC holds here because the quantity is both duality invariant in the β < 1, replica symmetric phase (where b = β), and has its temperature effectively renormalized as β → b = mβ in the β > 1, RSB phase, with the parameter m to be optimized. Duality invariance ensures that m = 1/β ⇔ b = 1 is optimal in the thermodynamic limit M → ∞, and thus insures freezing of the b-dependent, duality invariant quantities. Clearly the latter mechanism is not restricted to the free energy but applies in general. Therefore, we have shown that FDC would follow from a "1RSB-duality scenario", stating that observables O(β) will possess two properties simultaneously: (i) when calculated using 1RSB with group size m, O(β) = O ′ (b = βm) depends only on the renormalised temperature; (ii) in the β < 1 replica symmetric phase, At the moment we are unable to verify either FDC or 1RSB-duality scenario from the first principles. Nevertheless, in all known applications of FDC, the result can be recovered from a 1RSB calculation, and the above 1RSB-duality scenario holds. Indeed, the above discussion covers all (integrable) cases where the replicated partition function Z n can be calculated explicitly, and can be extended to min-max correlation as was shown in [35]. As to the other cases, such as the distribution of position of the global minimum [32,33] or Edwards-Anderson order parameter [34], one can argue they can either be reduced to the replicated partition function case (we shall see this for the minimum position in Sect. E), or to the variation thereof with respect to some parameter (the Edwards-Anderson case).
The above analysis also reveals the limits of when to expect the freezing phenomena. We see that such phenomena are in fine reducible to the Z n case, for which the UV contribution is well confined to a O(1) first moment shift (eq. (73)). This will no longer be the case for the observables related to higher minima, to which we turn our attention now.

C. Second minimum by 1RSB
From this section on we consider the higher order statistics proper, starting with the 1RSB derivation of the main result for the second minimum (24). The first step is similar to that of the global minimum value case, eq. (53): we start with writing H β (21) as the generating function (inverse Fourier transform) of replica averages: (75) Observe that W(j, ∆) is the partition function of the potential obtained from V j (the log-REM) by shifting a single value V j → V j + ∆. We shall call (j, ∆) a marker, of which j is the position and ∆ the shift-value. Therefore W is a special case of (57) in section (IV A), so as argued there, for any j fixed we can employ the 1RSB Ansatz in the replica sums over j 1 , . . . , j n . The latter form groups of m replicas, each group occupying a block. Now (75) contains one more sum over the marker position j. Yet, because of the presence of the difference [. . . ]| y−y0 ∆=+∞ , some j a must be equal to j to give a non-vanishing contribution. This implies that one (and only one) group out of n/m will be at the block containing j, and the sum over j reduces to choosing one group amongst the n/m and summing over the UV positions of the block that the chosen group occupies. By permutation symmetry between the groups, we can assume that the marker is near the group g = 1 and multiply the answer by the combinatorial factor n m C n,m . As observed above, W(j, ∆), W(j, ∞) and Z differ only by one Boltzmann factor, which affects only the intra-group interaction of the marked (g = 1) group. After collecting and factoring out all similar terms the (difference of) modified intra-group interaction isolated in D(ξ, ∆) takes the form given below: where w j (β, ∆) was defined in eq. (29). In (76), we recall that E(ξ), defined in (61), is the intra-group interaction of a group at block ξ without marker. We now compare (76) to (60). The two expressions share the same continuum integral part and most of combinatorial factors, though (76) acquires an extra factor n/m. As for the intra-group interactions, (60) has a product of n/m factors E(ξ g )'s, whereas in (76) one replaces the factor E(ξ g=1 ) with D(ξ g=1 ).
In summary, we have j W n (j, ∆) Now we apply (54) to express the above in terms of the generating functions (53) and (75): To bring (78) to the form (24) we need to write down the last ratio factor in (78) explicitly. Combining (77) and (61), we have: by recalling the definition of K β in (28). This finishes the derivation of (24). (49) We now derive the relation (49) used in our numerical study. For this we set y 0 = y in (24), and then integrate over y:

Remote second minimum
In Sect. III, we observed that the gap is exponentially distributed when measured between the global minimum V min and the remote second minimum V far min,1 defined in (50) (rather than the true second minimum V min,1 ). For this we note that the joint distribution of V min and V far min,1 can be accessed by the following modified version of (21) with y > y 0 (see also (23)): To treat this, we will follow the derivation just presented for H β , and indicate only the necessary modifications. First, Compared to (75), the potential values V k for 0 < |k − j| ≤ N are always shifted by ∆ = y − y 0 , in addition to V j which is still shifted by δ = y − y 0 and δ = +∞ (respectively in the two terms). Since 1 ≪ N ≪ M in eq. (50), when implementing 1RSB procedure, we can arrange the blocks in a way such that all the shifted sites form exactly one block, of which the site j = j(ξ, i = 0) summed over in (82) is at the centre. Note that this makes the division of the system into boxes dependent on the marker position, while the division is fixed once and for all in the previous derivations. Working through the arguments leading to (97), mutatis mutandis 6 , one arrives to a version of that equation, in which the intra-group energy of the marker's block (see (77)) is replaced by In terms of the generating functions (78) takes the following form: (86) 6 One encounters here the following technical subtlety: here, one sums over the marker's position j = 1, . . . , M , and constrains one replica group to be in the block centred at j by the subtraction argument analogous to (IV C); each of the other replicas groups chooses one block from N/M . Hence, the (M/N ) n/m factor in (76) and (60) will be replaced by M × (M/N ) n/m−1 here, giving rise to the factor N in (85) which replaces the sum N i=1 in (77). Now, the approach here can be used in the second minimum calculation above, and by doing that, we would obtain a different expression of (77), obtained by replacing N i=1 . . . (i) → N . . . (i = 0). However, we argue that such a difference becomes irrelevant in the N → ∞ limit, in which the local logREM is statistically translation invariant, see (27) and (12).
The remaining job is to calculate the last ratio. For this, observe from (85) that (i)D far (ξ, δ, ∆)/D far (ξ, ∆, ∆) is a function of δ − ∆; (ii) when δ = ∆, all potential values in the block are shifted by ∆, so D far (ξ, ∆, ∆)/E(ξ) = N e −mβ∆ . Combining the two steps, we have where F is some unknown function. Plugging this into (86), we obtain K far β (∆) = e −mβ∆ N (F (0)−F (+∞)) . Applying the analogue of (36) gives δ(g far ). Normalisation of the probability fixes the constant N (F (0) − F (+∞)) = 1, giving where the parameter m as a function of the temperature is defined in the equation (56). In the zero-temperature limit, this recovers the simple exponential gap distribution (51) for the gap with a remote second minimum. As one could have expected, for such a situation the UV correction factor becomes UV-data independent for any temperature.
From that angle such quantity reflects only uncorrelated REM-like features of the model, rather than the full logREM structure. Note that (88) provides an example of freezing that seems not to be accompanied by any obvious duality invariance, the feature shared by the generating function (20) of the uncorrelated REM [7].

V. GENERALIZATION TO HIGHER ORDERS
To generalize the calculations of the previous section to higher order statistics, we shall consider the following observable, which generalizes to the case of k ≥ 1 markers the function H β (y 0 , y) defined for k = 1 in (21): β→∞ −→ P(V j0 < y 0 , . . . , V j k−1 < y k−1 , all other V j > y) .
Here and below, we assume y 0 < y 1 < · · · < y k−1 < y. Taking derivatives with respect to y 1 , . . . , y k−1 and suppressing the argument for compactness, we obtain: where V min,j stands for the value of the (j + 1)-th ordered minimum, with V min,0 < V min,1 < V min,2 < . . . (as for Gaussian sequences the probability of two minima values to coincide is strictly zero); in particular, with such notations the global minimum is V min = V min,0 . In the above equation, the position of (k + 1)-th minimum is not fixed, while that of the first, second, . . . , and k-th minima are fixed at j 0 , . . . , j k−1 (with the order determined by y 0 , . . . , y k−1 ).
where * j0,...,j k−1 sums over distinct positions (j 0 , . . . , j k−1 ). A. 1RSB calculation of (89) In this section we demonstrate how to evaluate (89) (with marker positions summed over) within 1RSB framework. Similarly to (75), (89) generates partition sums with markers at positions j 0 , . . . , j k−1 : Note that W is identical to the partition function in (57). So, as we justified in Sect. IV A, we can treat the sum over replicas using 1RSB, i.e., considering only the dominant configurations where the replicas form groups of size m = min(1, 1/β) (eq. (56)). We then follow exactly the recipe (59) to replace the sum over replica positions by integrals over group macroscopic coordinates times sums over microscopic intra-group positions. As before, each marker position must coincide with some replica positions to give non-vanishing contribution to H β , due to the cancellations in (95).This imposes severe constraints on the sum over marker positions in (93); in particular, at the macroscopic level it is reduced to summing over all possible ways of assigning each marker to a replica-group, which is described by a mapping {0, . . . , k − 1} → {1, . . . , n/m}. Since we will eventually continue n and m to non-integer values, we need to work out the combinatorics explicitly. For this, we recall that any above mapping can be constructed by partitioning {y 0 , . . . , y k−1 } into p parts followed by the assignment of the p parts to distinct groups; the latter results in a combinatorial factor (n/m)(n/m − 1) . . . (n/m − p + 1) = (n/m) p in terms of the Pochhammer's symbol 7 . So the replica average in the RHS of (93) takes the form of a sum P(k) (n/m) p [. . . ], where P(k) denotes the sum over all partitions of {y 0 , . . . , y k−1 }, and [. . . ] denotes the contribution of each partition. The notion of a partition we use here should not be confused with the standard partitions of an integer: in the latter, the elements as well as the parts are indistinguishable objects, while in our case, only the parts are indistinguishable, but the elements are distinguishable. For example, when k = 3, p = 2 (see Fig. 6 for illustration), {{1, 2}, {3}} and {{1, 3}, {2}} are distinct partitions (although they give the same integer partition 3 = 2 + 1). In general, we describe a partition as follows: p denotes the number of parts, (k 1 , . . . , k p ) denote the numbers of elements in each part, so that k 1 + · · · + k p = k; for each part, let {y q,0 < · · · < y q,kq −1 } be the list in ascending order of its elements; since the parts are indistinguishable, we also order them by their first element: y 1,0 < y 2,0 · · · < y p,0 . In the following, the set of all partitions of k elements will be denoted P(k). Clearly, knowing a partition of {y 0 , . . . , y k−1 } is equivalent to that of {0, . . . , k − 1} = {j q,s : q = 1, . . . , p, s = 0, . . . , k q − 1}, and for any set of variables {x 0 , . . . , x k−1 } indexed similarly, we will adopt the short-hand x q,s = x jq,s .
Let us now consider the contribution of a given partition. Similarly to the second minimum case, the presence of markers does not change the Coulomb gas integral of the IR sector, but does modify the intra-group interactions.
Moreover, since more than one marker can be at one block we need to introduce the higher order generalizations of (29) to express the ratios of intra-group energies: With the above ingredients at our disposal, it is not hard to generalize the calculation of the second minimum and obtain the following: * j0,...,j k−1 Note that (98) defines an infinite series of UV correction factors depending on ℓ = 1, 2, . . . arguments and generalizing (28), to which (98) reduces for ℓ = 1. The RHS of (98) is symmetric in its arguments ∆ 0 , . . . , ∆ ℓ−1 ; it will prove helpful to assume the ordering ∆ 0 > · · · > ∆ ℓ−1 , in agreement with the ordering of the arguments appearing in (97).
In terms of the generating functions (53) and (93) Here we used the Laplace transform (54) that helps to transform (n/m) p = (nβ/b) p in (97) into (b −1 ∂ y ) p = (−1) p D β,p (the signs (−1) p are inserted for convenience later). This generalizes (24) (recovered for k = 1, p = 1) and is the main result of this work in full generality. We first consider implications of (100) in the β < 1 phase with the parameter m = 1. Plugging (96) into (98), we have Note that the formula (95) has a nice feature, namely that f (x 0 , . . . , x ℓ−1 )| bs xs=as = 0 whenever f does not depend on all of its arguments. So the term i / ∈{is} in the above equation always vanishes since it does not depend on any δ s . Each of the other terms depends on one (and only one) of the variables δ s , and thus vanishes when ℓ > 1. So we have This implies that only the partition with p = k and ∀k q = 1 survives in (100) (e.g. for k = 3 it is the rightmost one in Fig. 6) which makes the latter to simplify to * This relation generalises (31), and is also (super) universal, i.e. shows neither IR nor UV data dependence.

B. Extreme values process
Let us now discuss the implications of (100) at zero temperature. In all the β > 1 phase the differential operator D β>1,p (101) becomes To obtain minima value distribution we should combine (92) to (100), hence the need to calculate the derivatives ∂ y0...y k−1 . . . . Noting that the product p q=1 K β (y − y q,0 , . . . , y − y q,kq−1 ) contains every y k−1 and only once for each, we obtain where ∂ 0...ℓ−1 = ∂ ∆0 . . . ∂ ∆ ℓ−1 . The last mixed partial derivative can be evaluated more explicitly. For this we first calculate β → ∞ limit of K β from (98) and (96): where u min = u min,0 , u min,1 , . . . , u min,ℓ are the first (ℓ + 1) minima of the UV block and R denotes the 2 ℓ − 1 other terms generated by (95), which will be eventually annihilated by the mixed partial derivative ∂ 0...ℓ−1 (nonetheless, it will be useful to write them out explicitly for other purposes, see Appendix D 1, eq. (D14)). The situation is similar to (34): although one cannot cancel out the denominator, one can define a biased minima process (v min,s ) ∞ s=0 , v min = v min,0 < v min,1 < . . . , uniquely up to a translation, by the following property for any reasonable function of f of arbitrary number of variables. Note that (108) is generalizing (35) to the multivariable case; we will discuss it in more detail in VI B. In terms of (v min,s ) (107) is expressed as Note that (109) is only valid if ∆ 0 > ∆ 1 > · · · > ∆ ℓ−1 are ordered, and its RHS is not symmetric in ∆ 0 , . . . , ∆ ℓ−1 .
Nevertheless, it is a symmetric function of variables In particular, x 0 = ∆ 0 . This ensures that ∂ 0,...,ℓ−1 = ∂ x0,...,x ℓ−1 is well-defined independently of the differentiation order, so we can select any suitable order: Applying the above identity with (110) to (109) we obtain Here D(∆ 0 , . . . , ∆ ℓ−1 ), ∆ 0 > · · · > ∆ ℓ−1 is an observable of the process (v min,s ) and is the probability (density) of the event described by the last line. In particular, gives the pdf of the first ℓ − 1 gaps of (v min , s). Thus the knowledge of D(∆ 0 , . . . , ∆ ℓ−1 ) for all ℓ = 1, 2, . . . determines the process (v min,s ) completely. In fact the information in the set of D(∆ 0 , . . . , ∆ ℓ−1 ) is excessive, since e.g., by applying ℓ−1 s=0 ∂ ∆s to (113) and comparing with (115) (with ℓ ℓ + 1), we obtain a series of functional relations between D(. . . )'s Now plugging (112) back into (106), we obtain the final expression for the joint distributions of minimum values, quoted in (43) e yq,0−y D(y − y q,0 , . . . , y − y q,kq−1 ) . (117) As a check, let us retrieve the special case (32) from (117). A possible way to proceed starts with setting k = 1 so (112); applying then −∂ y we get directly (32). The second, perhaps more instructive, way is to set k = 2 and y = y k−1 . Then there are two partitions of {y 0 , y 1 }: {{y 0 , y 1 }} and {{y 0 }, {y 1 }}, so (117) can be calculated explicitly (we denote ∆ = y 1 − y 0 , and use repeatedly (113)): which is equivalent to (32) combined with (37) and (36). At this point an initiated reader may have recognized that (117) is precisely the order statistics of a randomly shifted decorated Gumbel Poisson point process (SDPPP), with the decoration process given by (108) (we recall the definition of SDPPP in Appendix C and why its extreme value statistics coincides with (117)). Thus, by arriving at (117) we have achieved one of the main conceptual goals of this work. Namely, we demonstrated that the 1-step replica symmetry breaking mechanism, properly understood and generalized, implies that the sequence of ordered minima of a general Euclidean-space logREM is described by SDPPP. To our knowledge, the approach is new and is a first clear demonstration of the "1RSB ⇒ SDPPP" relation. Our result, although not unexpected, does not seem to be addressed from such a viewpoint and in such generality in the literature; the characterisation of the decoration process (108) seems, to the best of our knowledge, new.
It is apparent from (117) that the order and gap statistics are primary objects of interest in our approach, and we feel they were not paid due attention in the existing SDPPP literature. To fill in this gap we shall study some implications of (117) in the Appendix D, which contains the derivation of (46) and (45).

A. A toy model
One of the most essential points we hope to have clarified in our article is the way to take UV limit data into account when performing replica calculations. Although the contribution of UV sector is captured by (98) in a somewhat sophisticated form, the predictions are amenable to direct numerical verifications, see Sect. III. Therefore, it is instructive too see how eq. (98) can arise in an ab initio calculation without any use of the replica-trick. For this purpose we consider below the following toy model. It is defined by M values of the random potential that form blocks of a fixed size N , V j = V (g,i) , g = 1, . . . , M/N, i = 1, . . . , N , such that different blocks are statistically independent and the correlations inside a block mimic those of the logREM: where f ii ′ is a symmetric positive definite N × N matrix, with N being fixed for simplicity.
This model retains the UV-limit correlations of general logREMs, but has trivial "IR" correlations. So the notion of the local logREM (u i ) N i=1 and its partition functions, defined in (27) and (96), still makes sense here: The toy model is exactly M/N independent copies of local logREMs, to each of which we added a constant centered Gaussian variable shift C g of variance 2 ln M ; shifts of different copies are independent. The independence among different blocks makes it possible to perform calculations without replicas. For this we need some technical result standard in the study of uncorrelated REM. Let C be a centred Gaussian of variance 2 ln M , and consider the following average over C: For convenience of the readers the asymptotic behaviour of (121) is analysed in Appendix F. The results can be summarized as follows: The above equations hold asymptotically as M → +∞ with x − f ∼ O(1) fixed, and m = min(1, 1/β) is the same as defined in (56). Coming back to our toy model, we start with the observable for the minimum value (20). By independence, it is given by M/N -th power of the contribution of an individual block (recall θ β (x) = exp(−e −βx )): Note that g β (y) is independent of the block index g, since all blocks are identical. The average [. . . ] here goes over the local logREM and over the local C g . Clearly, by (121), g β (y) = γ(y + β −1 ln z) (where the over-line in the RHS refers to the random variable z alone). Applying (122), we reduce (123) in the limit M → ∞ to In the β → ∞ limit this implies that the the minimum value for the toy model is distributed, as expected, according to a shifted Gumbel law, with the shift depending on the minimum value of the local REM potential z 1/β → β→+∞ exp(− min i u i ). Note that in comparison with logREMs (26) the shift for the toy model displays a different coefficient in front of ln ln M term, consistent with the fact that our toy model is not a logREM, according to the definition (1). Now to probe the higher order statistics, we calculate the observable (89), summed over marker positions. Observe first that the sum over markers positions in the LHS is factored into a sum over partitions indicating the occupation of the blocks, and the combinatorial factor M N p that chooses p blocks for the parts. The blocks without markers each yield the same contribution as in the RHS of (123), while the blocks with markers contribute a similar factor, with z replaced with the local partition function with a few value shifts in the arguments, i.e., w((∆ 0 , i 0 ), . . . , (∆ k−1 , i k−1 )). Finally, we get Note that the above equations are exact for finite M . In the M → ∞ limit, by applying (121) and (122) with x = y + β −1 ln w((i 0 , δ 0 ), . . . , (i k−1 , δ k−1 )) to (126) and using (95), we can simplify them to the following form: j0,...,j k−1 H β (y 0 , j 0 , . . . , y k−1 , j k−1 , y) = (−1) k (y − y q,0 , . . . , y − y q,kq−1 ) , To compare with (100), we shall apply the differential operator D p = −∂ y (1 − ∂ y ) . . . (p − 1 − ∂ y ) (eq. (101)) to G β (y) (eq. (124)). One can check by recursion that Combining this with (127) and (128), one shows with some rearrangement that j0,...,j k−1 H β (y 0 , j 0 , . . . , y k−1 , j k−1 , y) = (−1) p+k These equations are identical to (100), with the same UV sector contribution (98), except that the UV block size N is fixed here. This result is obtained without any use of the replica trick, by a well-controlled calculation.on a variant of uncorrelated REM that mimics the correlations of logREMs at the UV scale, so the coincidence of the results confirms the validity of the replica approach to UV sector.

B. Biased minima: local-global equivalence
In a log-correlated model the scale separation between UV and IR is only emergent in the thermodynamic limit. Note that, in contrast to the toy model considered above, for log-correlated models the size of UV block N is not fixed and one should consider it infinite in the thermodynamic limit as well. As we have however seen in Sect. II A, in general, the covariance matrix of the local logREM has a logarithmic behavior at large distance.
Therefore, when N → ∞ the local minimum v min will have the same log-REM signatures as V min . Moreover, the minima in the local REM can be separated by distances diverging in the thermodynamic limit, so the gap process of local minima (u min,k − u min ) k≥1 can be in practice hard to distinguished from that of the global minima (V min ) k≥1 of some logREM. The worst scenario seems to arise for the case of the interval model (with c 1,2 = d 1,2 = 0, see Sect. II A 4), where the local logREM looks identical to the global one with an appropriate choice of C N . These facts may cast legitimate doubts on the consistency of our strategy of separating UV contributions.
To this end, we are able to demonstrate that statistical properties of the biased minima are unchanged by a local-global replacement (u min,k ) (V min,k ). To show this, it suffices to calculate which differs from (113) by the local-global replacement (and by setting ∆ k−1 = 0, which causes no loss of generality since we cover all k's). The statistics of (V min,k ) is given by (117) (where we take now y = y k−1 ) which implies e −∆q,0 D(∆ q,0 , . . . , ∆ q,kq −1 ) (133) Here, y stands for V min,k−1 and is integrated over. To proceed further, recall by (105) that D p = −∂ y (1 − ∂ y ) . . . (p − 1 − ∂ y ). Note that for any function f (y) satisfying f (y)e −y → 0 when y → ±∞, we have R e −y (1 − ∂ y )f (y) = − R ∂ y (e −y f (y)) = 0. Applying this to f = −∂ y (2 − ∂ y ) . . . (p − 1 − ∂ y )G ∞ , and assuming fast enough decay of G ∞ (y) with y → ±∞ (we will come back to this crucial point below), we obtain So the sum (133) has in fact only one surviving term given by the trivial partition, for which we have p = 1 and ∆ 1,s = ∆ s , s = 0, . . . , k − 1, with ∆ k−1 = 0. Now by (22) we have R e −y D 1 G ∞ (y) = exp(−V min ). Incorporating these into the RHS of (133) simplifies it to This means precisely that statistics of the biased minima can be obtained from either local logREM minima following (113), or from the original logREM minima following (132), with the results being identical. Moreover, the calculation tells us that the biased minima process contains solely UV information: indeed, the only contributing partition corresponds to the configurations in which all the minima in consideration are in one UV block. Therefore, even if one supposes that further scale separation emerged inside the blocks themselves, so that the local minima (u min,k ) acquire an addition partitioning structure of the global minima, only the trivial partition would contribute to the biased minima process (v min,k ). This fact makes it clear that the assumption of an additional scale separation inside the UV blocks is de facto immaterial. As promised, we will now justify the assumption of a fast decay of G ∞ (y → ±∞). It is clear that the sufficient and necessary condition of (134) is e −y ∂ y G ∞ (y) → 0 as y → ±∞. The right tail +∞ poses no problem, but the left tail does: it is well known that the asymptotic Carpentier-Le Doussal (CLD) tail ∂ y G ∞ (y) ∼ Aye y holds for generic logREM [5], so e −y ∂ y G ∞ (y) ∝ y does not vanish. However, the CLD tail holds down to y = −∞ only in the thermodynamic M → ∞ limit; for any system of finite size M , when y/(2 ln M ) ≫ 1, the tail crosses over to a Gaussian tail ∼ B exp(−Cy 2 ) 8 , which still decays to 0 after multiplying by e −y . This argument justifies eq. (134), and provides a further insight: the exponential weight exp(−u min ) of the biased minima process probes configurations with atypically negative values beyond the logREM scaling regime y = −2 ln M + o(ln M ). Given the correlation structure of logREM, multiple occurrence of values much more negative than the typical scale −2 ln M is much more likely to happen at nearby points compared to that of ∼ −2 ln M values. This reinforces the claim that the biased minima process concerns exclusively the UV data, and explains the numerical observations made in Fig. 4 (Left): its numerical measure has faster convergence with respect to N , but the relevance of atypical events means that higher statistics are necessary.

VII. CONCLUSION AND OUTLOOK
In this paper we systematically applied the replica symmetry breaking (RSB) formalism to the extreme order statistics of a generic (see Sect. II A for definition of scope) logarithmic random energy models defined on Euclidean spaces. As is well-known, for all such models the one-step replica-symmetry broken solution (1RSB) characterized by the familiar overlap distribution P (q) = mδ(q) + (1 − m)δ(q − 1), m = min(β, 1) (see Sect. IV A and Appendix A, in particular (A5) and (56)) captures correctly the thermodynamics of the model. Our main finding is that after appropriate ramifications the same 1RSB mechanism allows one to efficiently analyse fluctuations in the ordered sequence of the deepest minima in the model. These fluctuations turn out to be crucially dependent on the IR and UV scaling limits (Sect. II A) of particular models, and in a peculiar way (see Sect. V). Namely, the process of minimal values form a randomly shifted and decorated Gumbel Poisson point process (see Appendix C). The random shift depends only on the IR data (up to a translation, it has the same law as the critical temperature free energy). The decoration process depends only on the UV data, and is identified with the biased minimum process described in Section VI B. This picture is in full agreement with that obtained in [49] and rigorously in [52,53] for hierarchical logREMs. Motivated by recent work of [46,54], we also extend the results to include minimum positions in Appendix E.
We worked out in a considerable detail the implications of this picture for the form of the joint distribution of the first and second minima in Section II B, and provided two numerical tests for the circular model in Section III. For the first gap, see Fig. 4, we tested direct numerical simulations of the original model vs. our analytical prediction 8 To understand this, it is instructive to have a brief look at a similar tail for a uncorrelated REM, which has the shape e −y . To this end we calculate θ(V min − y) for a uncorrelated REM. Letting y = −2 ln M + y 0 , we have θ(V min − y) = exp(− exp(y 0 − y 2 0 /4 ln M )) ∼ 1 − exp(y 0 − y 2 0 /4 ln M ) as y 0 → −∞. If y 0 / ln M does not vanish in the M → ∞ limit, the quadratic term cannot be neglected.
(combined with numerically evaluated UV contributions), and also analyzed a situation where the decoration process is trivial and the first gap is exponentially distributed. For the fluctuation of the second minima, see Fig. 4, our prediction (48) (in general (39)) constrains its distribution to a one-parameter family, where the parameter is given by the mean value of the first gap. We also derive a generalization of (39) to arbitrary order statistics: the distribution of the k-th minimum is determined by that of the first minimum and the mean value of k − 1 first successive gaps, in a universal way encoded in the generating function (45). We also clarified the relation of 1RSB approach to that based on freezing, in particular, the freezing-duality conjecture, see Section IV B 2. The main message is that FDC holds true for observables probing exclusively IR data. When the UV data is relevant, non-trivial temperature dependence takes place in the β > 1 phase.
What is left for future work is the remaining major challenge of performing analytical calculations of the decoration/biased minimum process. Although this process is fully characterized for BBM in [52,53] and for discrete GFF in [54], we are not aware of any quantitative prediction of the gap distributions, besides the trivial decoration process v min = 0, v min,1 = v min,2 = · · · = +∞, which is obtained when the minima are constrained to be far away from each other (see Sect. III for precise definition). At least two approaches may be anticipated towards that goal: performing perturbation expansion around the trivial case to obtain approximate solutions, or relating the decoration process to some discrete integrable model (recall that explicit predictions of the minimum value distribution have so far relied on exact evaluation and subsequent analytical continuation of the Coulomb gas integrals (10)). Implementing either of the approaches would yield new and interesting results.

Optimization of the RSB order parameter
Now, we can rewrite the replica sum Z n (52) as an integral over overlaps between pairs of replicas, which form the overlap matrix Q ab : Here DQ = a<b dQ ab and the integral is taken over all matrices satisfying Q aa = 1 and Q ab = Q ba ∈ [0, 1].
Comparing to the treatment in the main text all the UV and IR data give o(ln M ) corrections which are not explicitly written here. This is done in order to focus on the extensive quantities (note, however, that those corrections are the most essential ingredient when one studies fluctuations we will resort to the assumptions of the RSB Ansatz. The first assumption is that Q ab has the Parisi's RSB form, determined by the overlap distribution function P > (q): This is the (continuous) full-RSB Ansatz. Note that P > (q) depends on n as well as q. It is known (and can be shown by the standard replica trick) that P > (q)| n→0 gives the probability that two particles thermalized according to the Boltzmann-Gibbs measure generated by the same random potential occupy positions with overlap > q. So P ′ > (q) ≤ 0 for n → 0.
In terms of the overlap distribution P > , the energy part of the action is Note that the energy part depends on the definition of the overlap in terms of the covariance (A1), but not directly on f (q). Now we consider the entropy part S[Q]. To simplify, let us first restrict our attention to the discrete RSB (and then take the continuous limit), where P > (q) is a staircase-like function: P > (q) = p i , q ∈ (q i , q i+1 ) , 0 ≤ i ≤ s , 1 = p 0 > p 1 > · · · > p s > p s+1 = 0 , 0 = q 0 < q 1 < · · · < q s ≤ q s+1 = 1 , (A7) s being the number of steps. The discrete RSB has an equivalent description in terms of hierarchy of replica groups. To spell it out, let us choose an arbitrary replica index a as a reference. Let i = 0, . . . , s + 1. Among the n − 1 other replicas (n − 1)p i of them have overlap ≥ q i with a (for n large the inequality is saturated by (A7)). Therefore, the i-group formed by these elements and the added replica a is of the size For i = s + 1 the group consists of a single replica: m i = 1; while for i = 0 the group contains all the replicas. Now, by (A2), the members of the i-th group occupy a volume of M 1−f (qi) (number of points j such that q(j, j a ) > q i ). So seen as a whole, the group's position is determined only up to M 1−f (qi) . On the other hand, for any i > 0 the level-i group is contained in a level-(i − 1) group. The entropy contribution of a level-i group is the log-volume available for different positioning of the group as a whole (i.e. the volume of the level-(i − 1) group containing it), minus the log of its position resolution (i.e. its own volume; what happens at smaller scales will be taken care of in contributions of smaller groups): a being arbitrary, what we said holds for any i-group; in particular, there are n/m i of them in total. So we can calculate S[Q], given by the sum of (A9) over all groups, as follows: Now we can take the s → ∞ limit where P > (q) tends to any non-increasing function; in that case the RHS still makes sense. Combining with (A6), we have where m(. . . ) is defined in (A8). Considering ln M as a large parameter we can employ the saddle-point method to evaluate the integral in (A3), and thus have to minimize H[Q] with respect to P > (q). By (A8), we can equivalently optimize with respect to m = m(P > (q)) for n integer, and then continue the answer to n → 0 + . This operation involves two subtleties, both well known in general replica calculations [67]. The first concerns the range of the group size m. Naïvely, m ∈ [1, n] when n is an integer; however, when n is continued to n → 0 + , we demand m be still between 1 and n, i.e., m ∈ (0, 1]. The second subtlety is that after the continuation the minimization becomes maximization. To see this one may use (A8) and integration by parts to rewrite eq. (A11) as where we recalled f (0) = 0, f (1) = 1 (A2) and P > (1) = 0 from (A5). Then we observe that the P > (q) dependent-part is proportional to n(n − 1), which changes sign when going from n > 1 to the range 0 < n < 1.
With these technical points in mind, the actual optimization of (A11) is rather straightforward. In fact, the problem amounts to the maximization of h[m] = −n(β 2 m + f ′ (q)m −1 ) with respect to m ∈ (0, 1], for each q ∈ [0, 1]. The solution is Note that the convexity assumption f ′′ (q) ≥ 0 ensures dm/dq ≥ 0, so that P ′ > (q) ≤ 0 for n → 0 (since m(P ) = 1 + (n − 1)P ). Applying (56) to (A11) and (A3), and using the saddle point approximation, we have the extensive free energy We note that such a formula for free energy appears also for the generalized Random Energy Model [2,3,68].

LogREM
In this subsection we set f (q) = q to specify the above results to logREM. By (A13), m(P > (q)) = min 1, β −1 is a constant (with respect to q), which coincides with (56). Eq. (A14) recovers also the logREM free energy density (26). Finally, applying (A8) with n → 0 + , we obtain The above result is then the well-known zero-one law of overlap in logREM: the two particles have zero overlap with probability min(β −1 , 1), and have unity overlap with probability max(1−β −1 , 0). It was first calculated on hierarchical logREMs [4]. Later it is argued to hold for Euclidean ones in [5], Sect. III.E.2. More recently, it has been proved quite generally [44]. The RSB analysis here does not assume hierarchical structure of the logREM. Nevertheless, inherent to any RSB Ansatz is the assumption of hierarchical organisation (overlap structure) of the replicas. Let us identify the logREM result with the 1RSB Ansatz described in Sect. IV A. When β < 1, P > (q) is of the form (A7) with s = 0, q 0 = 0: this is the replica symmetric solution such that distinct replicas have overlap 0, i.e., their separations are of the order of IR scale. When β > 1, we have s = 1, q 0 = 0, q 1 = 1 and p 0 = m−1 n−1 = 1/β−1 n−1 (note that p 0 ∈ [0, 1] only when n → 0). This is the 1-step RSB solution: the replicas form groups of size m given by (56). Distinct groups have overlap 0, i.e., are separated by system size, while replicas of the same group have their mutual overlap equal to unity, so their separation is of UV (lattice spacing) scale.

Multi-scale generalization
Going beyond logREM, let us discuss the physical significance of general f (q), by relating to the works [65,66], which considered multi-scale generalisations of Euclidean logREMs. We recall that these models have the following covariance matrix (see eq. (13) of [66]) : where Φ satisfies Φ ′ (y) ≥ 0, Φ ′′ (y) ≥ 0. We demand moreover Φ(0) = −1 and Φ(1) = 0 so that V j,M c = 2 ln M (2) still holds, and that points separated by the system size are almost non-correlated . One can check that for logREM, Φ(y) = y − 1. Using (A1), it is direct to show that (A2) is satisfied with i.e., they are essentially inverse of each other. Note that the properties of Φ(y) above are equivalent to those of f (q).
(A13) implies that every value of f ′ (q) is a critical inverse temperature ,i.e., a non-analyticity of m.
It is helpful to consider first the case of discrete s-step RSB, which can be retrieved by a piece-wise linear function determined by its cusps: where we define β −2 s+1 = 0 by convention. One may check that (A19) and (A17) imply (A18). Note that s i=1 g i = β −2 1 . V j,M can be seen as the independent sum of s "logREM"s, of which the i-th is smooth below the scale M 1−fi and log-correlated above. As the temperature T = 1/β decreases to < 1/β 1 , the particle starts to be confined in a block of length M 1−f1 , inside of which the first "log-REM" becomes smooth and no longer affects the thermodynamics. The same repeats at scale M 1−f2 as one continues to cool down to T = 1/β 2 , and so on. After the s-th (last) transition, the confinement scale becomes M 1−fs = 1, i.e., the UV scale. Now, the generic continuous RSB can be seen as the limit in which s → ∞ and the transition points β i and β i+1 are separated by O(1/s). One has a replica symmetry phase for 0 ≤ β ≤ β min = f ′ (0), beyond which we have full RSB phases with a continuum of phase transitions, until β = β max = f ′ (1). Finally, if β max < +∞, for β > β max one has a frozen phase: from (A14) free energy density is The free energy density becomes independent of β and gives also the leading behaviour of the minimum. Leaving detailed analysis to future work, we expect the full minima process to be a randomly shifted decorated Ruelle cascade [69], where Ruelle cascade is the analogue of Gumbel Poisson point process for full RSB, and the random shift and the decoration behave as in the 1-step case considered in this work. However, the frozen phase can disappear if β max = +∞, which is the case for example when there is a power law f (q) ∼ 1 − (1 − q) α , q → 1 with 0 < α < 1. In this case the zero temperature becomes critical, and describing the minima process is an interesting question. In fact, even the absolute minimum is expected [66] to have different universal behaviour from logREMs, e.g., the 3 2 ln ln M correction is expected to be no longer true. To this end one may mention also a mathematical work [70] investigating related questions, though presented from a somewhat different perspective. The Edwards-Anderson's (EA) order parameter for the circular model (see Sect. II A 3 for its basics) has been calculated in [35] using the freezing-duality conjecture (FDC). Its numerical check made the most convincing case for FDC; yet, the way of applying FDC there was quite tricky. Here we reproduce the result using 1RSB. In fact, we shall consider a more general class of models that we call spherical logREMs. They are Euclidean logREMs whose geometric manifold X (see Sect. II A 1) is contained in the unit sphere of a Euclidean space where (ξ.ξ ′ ) denotes the Euclidean inner product. This is a vast class of models, since there is no further restriction on IR data 11 . For example, our analysis applies to the 2d GFF-logREM on the sphere, or variants of circular models related to Morris integral. We consider the following replica sum Here O = j O j e −Vj /Z is the thermal average. Recall that according to the 1RSB Ansatz, the sum (B2) is dominated by configurations where the replicas form groups of size m = min(β −1 , 1): each group occupying a O(1) block corresponding to a continuum position ξ g (g = 1, . . . , m) on the spherical manifold X , and the sum ja → comb. dµ(ξ 1 ) . . . dµ(ξ n/m ) is replaced by a sum over replica grouping configurations (involving a combinatorial factor) and a continuum integral over group positions (times UV contributions). For the simple replica partition function Z n (60), the former sum gives merely a factor C n,m (58), the only novelty of X β (n) with respect to Z n is that there are two kinds of grouping configurations: for m−1 n−1 C n,m configurations, the replicas 1 and 2 can be in the same group, so ξ j1 = ξ j2 ⇒ (ξ j1 .ξ j2 ) = 1 and the continuum integral is that of Z n ; for the remaining n−m n−1 C n,m configurations a new continuum integral defined in (B4) (with n/m variables and with renormalized temperature b = βm) replaces that of Z n . The combinatorial factors above can be understood by a "probabilistic" argument: if all C n,m groupings are equi-probable, the probability that replicas 1 and 2 be in the same group is m−1 n−1 . Note also that the inserted observable (ξ j1 .ξ j2 ) does not probe any UV structure (i.e., the scalar product is unity if the two replica belong to the same group), so the UV contributions to X β (n) is identical to that of Z n . Summarizing, we have the following: Here the two terms in the brackets represent contributions of two cases discussed above, and X β (k) is the ratio between the new continuum integral and the standard continuum integral of the logREM, Z k b . In the case of the circular model we recall that the latter is given by (15). Using the variable t = −nβ, and defining free energy as F = −β −1 ln Z, one can rewrite (B3) as: In particular, a pleasing general result can be obtained by setting t = 0 in (B5): In other words, we predict that for all spherical logREMs the EA order parameter is linear in temperature T = 1/β in the whole β > 1 phase. Since at T = 0 the EA order parameter is always 1, only the slope is model-dependent. Now we specify our consideration to the particular case of the circular model, to compare our result here with that of [35], Eq. (42). Its LHS is equal to that of (B5) above. In [35] an exact Coulomb gas integral is used to solve the high temperature phase. We quote the result for β = 1 (assuming continuity at that point): X 1 (n) = (2 − n) −1 . Plugging this into the β > 1 case of (B5), we obtain in agreement with eq. (42) of [35]. Given the β < 1-phase solution, obtained due to exact evaluation of the Coulomb gas integrals, a non-trivial application of the FDC was needed to obtain the correct β > 1 phase answer in [35]. Here, the first-principle 1RSB calculation provides the correct non-analytical continuation, with the knowledge of the β = 1 solution alone.
Appendix C: Joint order statistics from randomly shifted decorated Gumbel Poisson point process Here we recall the notion of a randomly shifted decorated Gumbel Poisson point process (SDPPP) and calculate its minima value distribution. A reader recognising that (117) describes a SDPPP may safely skip this subsection.
A SDPPP is by definition generated by three steps. First, one generates a Gumbel Poisson point process (PPP) with the density e x dx . Let us recall a defining property of the PPP: the joint probability of the number of points N I1 = n 1 , ..N I k = n k in a set of disjoint intervals I 1 , ..I k is given by k j=1 1 nj ! (N j ) nj e −Nj where N j = x∈Ij e x dx. From this, denoting the sequence of ordered minima values as X 1 < X 2 < . . . , it is easy to write down the following joint probability density, for x 1 < · · · < x p < x a property which completely describes PPP and can equivalently be taken as its definition. In particular, −X 1 has a Gumbel distribution, hence the name.
To define the decorated Poisson point process (dPPP) without random shift we need to specify a decoration process. It must be a point process that contains 0 and contains only non-negative points, but can be otherwise arbitrary, so long as we can order its elements 0 = Y 0 < Y 1 < Y 2 < . . . . Similarly as in (113), the decoration process (Y s ) ∞ s=0 is completely determined by the following observables: As is discussed in the main text (around (115)), this is equivalent to giving the joint pdfs of Y 0 , . . . , Y ℓ−1 for any ℓ = 1, 2, . . . . We now make an independent copy of (Y s ) ∞ s=0 for each q = 1, 2, . . . , denoted ((Y q,s ) ∞ s=0 ) ∞ q=1 , i.e., for any q, (C2) would hold if one replaces Y s Y q,s . The copies are also jointly independent from the PPP. The dPPP is then obtained by replacing each point of the PPP X q by the points X q + Y q,0 , X q + Y q,1 , X q + Y q,2 , . . . , i.e., the q-th copy of the decoration process shifted by X q .
Finally, the SDPPP is obtained by shifting all the points so far by one random variable F independent of everything else.
Let us denote X min,0 < X min,1 < . . . the minima of the dPPP, and thus X min,0 +F < X min,1 +F < . . . the minima of the SDPPP. We first set the random shift aside and calculate the joint distribution θ(X min,k − y) k−1 s=0 δ(X min,s − y s ), where y 0 < · · · < y k−1 < y, in order to compare with (117). By the above definition this is given by the sum over the partitions of y 0 , . . . , y k−1 , indicating what points are in the same decoration. So let us fix a partition {y q,0 , . . . , y q,kq−1 } p q=1 and calculate its contribution. First one has the factor arising due to PPP, which is equal to the probability (density) that X q+1 > y and that X q = y q,0 , q = 1, . . . , p, which is given by (C1) with x = y and x q = y q,0 , q =, . . . , p. Then one has a factor coming from each decoration group, which is the probability that x q + Y q,kq > y and that x q + Y q,s = y q,s with s = 0, 1, . . . , ℓ = k q − 1. This is given by (C2) with ∆ s = y − y q,s , s = 0, 1, . . . , ℓ = k q − 1.
In summary, we have: δ(X min,s − y s ) = P(k) P(y 1,0 , . . . , y p,0 , y) using that D p exp(− exp(y)) = exp(py − exp(y)) where D p is defined in (105). The above formula is identical to (117) and (113), with G ∞ = Gum (the REM-Gumbel G ∞ ) and the decoration functions D = d. Now we incorporate the random shift, that is replace y ... → y ... − F in the above calculation and average over F . Namely, denoting the probability measure δ(F − f )df = dP(f ), we have Comparing with (117), we identify Gum(y − f )dP(f ) = G ∞ (y) . Since Gum(y) = θ(−Gum > y) where Gum is a standard Gumbel distribution and G ∞ (y) = θ(V min > y) (V min here is the minimum of the log-REM giving G ∞ ), we get the convolution relation V min = F − Gum. Comparing the latter with the freezing scenario (see (25) and Sect. IV B 1), we identify the random shift F of the SDPPP with the critical-temperature free energy of the logREM. Clearly, since the random shift F has no effect on the gap distributions, the latter will be shared by the log-REM and the SDPPP with the same UV data/decoration process, and given by (D2).
where we have used that for any K(j) and a p the following holds as an identity of formal power series (i.e. term by term): In particular, taking a p = (p − 1)!, which gives F a (x) = − ln(1 − x), and recalling the equation (D5), we have: which provides the marginal distribution of all global gaps, obtained by differentiation: δ(g k − ǫ) = ∂ 2 ǫ (g k − ǫ) + (recall K(j) depends on ǫ by (D7)).
We now show that (D12) is equivalent to (46) in Sect. II, by calculating explicitly K(j). For this we shall write down explicitly (107), by expressing the 2 j terms generated by equation (95)  The last equality follows from the definition of the biased minima process (108). Combining with (D12), we have which is just (46). In particular, at order t 1 , we have (exp(−ǫ) − exp(v min − v min,1 )) + = (g 1 − ǫ) + , which recovers (38). Therefore, (D15) is a non-trivial generalisation of (35), which gives the distribution of all gaps in terms of the decoration process.

Order statistics
Now we turn our attention to the marginal distribution of (k + 1)-th minimum V min,k . For this, it is useful to observe that (D1) provides two different expressions of joint pdf (a special case of this observation for first and second minima is in Sect. V B): -The first is obtained by considering (D1) with y = y k−1 . Performing such a substitution makes the θ-term in (D1) always 1, so we are left with the k − 1 δ's. The result is: Now, to obtain the marginal distribution of the (k + 1)-th minimum, we integrate out ∆ s = y − y s , s = 0, . . . , k − 1 in the domain S k = {∆ 0 > · · · > ∆ k−1 > 0} (in this subsection we set ǫ = 0) in (D17). The result can be rearranged into δ(V min,k − y) + where the second to third line integrates out ∆ k−1 while keeping ∆ s −∆ k−1 , s = 0, . . . , k−2 fixed (by applying Newton-Leibniz). On the other hand, it is not hard to see that if we integrate (D16) over ∆ s = y k−1 − y s , s = 0, . . . , k − 2 while keeping ∆ k−1 = 0, we also obtain the last line of (D18). This implies the following recursion relation: e −∆q,0 D(∆ q,0 , . . . , ∆ q,kq −1 ) , Writing the generating function ∞ k=1 t k . . . of both sides, and then applying the generating function relation (D10) with a p = D p G ∞ (y), and after some rearrangement we obtain a power series generating all higher order statistics: Now since D p = −∂ y (1 − ∂ y ) . . . (p − 1 − ∂ y ) (see eq. (105)) we can apply the binomial formula to the RHS of above, regarding ∂ y formally as a variable: and then apply (D12) (with ǫ = 0) to arrive at ∞ k=0 θ(V min,k − y)t k = 1 which is (45) in Sect. II C. Expanding it out as a power series with respect to both t and in ∂ y , we see that the cdf of any V min,k is a linear combination of derivatives of G ∞ (y): Using computer algebra software this procedure can be easily continued to arbitrary order. Note that (D24) and (D25) (for k = 1) covers (39). Although from (D16) it is not hard to see that (D24) must hold for some set of coefficients c k,s that depend only on the biased minima process, it is not obvious that these coefficients are polynomial functions of the mean values of the gaps. As another check, let us apply R dy(−y)∂ y to both sides of (D23): The LHS gives (V min + g 1 + · · · + g k )t k .
For the RHS, recalling from (22) that −G ′ ∞ (y) is the pdf of V min , we have (V min + g 1 + · · · + g k )t k .
The two sides give the same result as expected. Although (D23) is a non-trivial general property of SDPPP, and can be used for numerical tests, it is highly unlikely to be a sufficient characterization of SDPPP, since it concerns only marginal distributions. Taking the β → ∞ limit, evaluating the mixed partial derivative ∂ y0 . . . ∂ y k−1 , and setting y to y k−1 give rise to joint position-value distribution of the first k minima: ∂ y0...y k−1 H ∞ ((η q ), (y q,s ), y) = θ(V min,k − y) k−1 s=0 δ(V min,s − y s ) p q=1 δ(Π q , η q ) (E2) In plain terms the above is the probability of the following event: the k first minimal values are given by y 0 < · · · < y k−1 , and the (k + 1)-th minimum is larger than y; moreover, the minima {y 0 , . . . , y k−1 } form p blocks, described by a partition {y q,s : q = 1, . . . , p, s = 0, . . . , k − q − 1}, such that the q-th block is located at η q . Let us emphasize that Π 1 , . . . , Π q ∈ X are not macroscopic positions of the first q minima, but the first q minima that are locally minimal.
The normalisation of the Dirac-δ's on the positions will be clarified below, see (E16 which relates the random measure ρ with I ∞ which is in turn related to continuum Coulomb gas integrals (E8) via (E13). Let us briefly discuss some consequences of this correspondence. First, we look at the total mass Z[ρ].
Integrating the above over all η q and comparing with (E15) we obtain Applying Laplace transform to this formula and comparing to (25) gives us where F (. . . ), (26) is the non-fluctuating part of the minimum, and contributes a deterministic shift to ln Z[ρ]. Up to translation, − ln Z[ρ] has the same distribution as the critical temperature free energy, as one can see from the above equation and (70). Note that up to now in the present paper the over-line in Z t (11) was a formal symbol and did not correspond to any genuine average, so (E20) gives it a probability interpretation, in terms of the random measure ρ. Thus, from now on we shall identify Z with Z[ρ] exp(F (M, ∞, f ij )). Next by applying the Laplace transform directly to (E18) and comparing with (E13) gives after cancelling out the extensive free energy factors by (E20) and (71), we obtain Eqs. (E21) and (E20) express properties of the random measure ρdµ 1 in terms of analytical continuation of Coulomb gas integrals (E5), (E8). Many general transformation properties of ρ shown in [47] can be understood in terms of change of variables in (E5). In particular, if we set −t = p, and apply (E5), (E8), (E21) andρ = ρ/Z[ρ] in (E18), we obtain p q=1 ρ(η q ) ρ = q<r exp(C(η q , η r )) .

Relation with freezing approach
[32] and [33] studied also the position of the global minimum, using the replica-freezing approach, of which the derivation of (E14) above is a generalization to higher order statistics. To see this more clearly let us restrict (E14) to the particular case where k = 1 (so p = 1 and the partition is the trivial one) and y = y k−1 . Applying inverse Fourier transform to both sides and using (E13) we have Setting t = 0 we get the minimum position distribution see (E8) and (E21). By (E5), the latter is an analytically continued Coulomb gas integral over n → −1 moving charges and one charge fixed at η, and the temperature fixed to the critical value: we obtained the same replica-trick recipe of minimum position in [33], sect. 2.3, which used the freezing-duality conjecture. Unfortunately, we know of no exactly solved integrals that give non-trivial position distribution. We note however that [33] calculated integer moments of the minimum position of the interval model (see Sect II A 4 for definition). Although the results are in principle equivalent to knowing δ(Π 0 − η), an explicit expression of the latter is still lacking.
where . . . stand for the omitted contributions from other poles (when they arise), as those can be easily written down and checked to be immaterial for our goals. Let us again consider x = f + δ where δ is of the order of unity but f < 0 is defined by the identity: Note, that such f together with x < −2β ln M immediately implies 0 < β < 1. We also easily check that the last term in (F8) is subdominant. Finally we get which is the β < 1 case of (122).