Ada Altieri

Félix Roy

Chiara Cammarota

Giulio Biroli

Properties of equilibria and glassy phases of the random Lotka-Volterra model with demographic noise

In this letter we study a reference model in theoretical ecology, the disordered Lotka-Volterra model for ecological communities, in the presence of finite demographic noise. Our theoretical analysis, which takes advantage of a mapping to an equilibrium disordered system, proves that for sufficiently heterogeneous interactions and low demographic noise the system displays a multiple equilibria phase, which we fully characterize. In particular, we show that in this phase the number of stable equilibria is exponential in the number of species. Upon further decreasing the demographic noise, we unveil a Gardner transition to a marginally stable phase, similar to that observed in jamming of amorphous materials. We confirm and complement our analytical results by numerical simulations. Furthermore, we extend their relevance by showing that they hold for others interacting random dynamical systems, such as the Random Replicant Model. Finally, we discuss their extension to the case of asymmetric couplings.

Introduction -Lotka-Volterra equations describing the dynamics of interacting species are a key ingredient for theoretical studies in ecology, genetics, evolution and economy [START_REF] May | Theoretical ecology: principles and applications[END_REF][START_REF] Faust | [END_REF][3][START_REF] Goodwin | Chaotic economic dynamics[END_REF][START_REF] Kessler | [END_REF][6]. Cases in which the number of species is very large are becoming of general interest in disparate fields, such as in ecology and biology, e.g. for bacteria communities [7,8], and economy where many agents trade and interact simultaneously both in financial markets and in complex economic systems [9,[START_REF] Thurner | IFP/FGS Working Paper[END_REF]. The theoretical framework used in the past for a small number of species is mainly based on the theory of dynamical systems [START_REF] Barreira | Dynamical Systems[END_REF][START_REF] Macarthur | [END_REF][START_REF] Tilman | Resource competition and community structure[END_REF][START_REF] Ruan | Delay differential equations and applications[END_REF][START_REF] Vano | [END_REF][16]. When the number of ordinary differential equations associated with the Lolta-Volterra (LV) model becomes very large, i.e. for many species, methods based on statistical physics become ideally suited. Indeed, several authors have recently investigated different aspects of community ecology, such as properties of equilibria, endogeneous dynamical fluctuations, biodiversity, using ideas and concepts rooted in statistical physics of disordered systems [START_REF] Kessler | [END_REF][17][18][19][20][21][22][START_REF] Pearce | Proceedings of the National Academy of Sciences[END_REF][START_REF] Roy | [END_REF][25][26][27]. Similar investigations have been also performed for economic systems [28]. The complexity of dealing with a large number of interacting species can actually become a welcome new ingredient both conceptually and methodologically. In fact, different collective behaviours can emerge. As it happens in physics, such phases are not tied to the specific model they come from, hence they can be characterized and characterize systems in a generic way [START_REF]Along a similar vein, models of liquids and crystals that are used in physics are disparate, approximate and often inaccurate with respect to real system[END_REF]. From this perspective, it is natural to ask which kind of different collective behaviors arise from LV models in the limit of many interacting species and what are their main properties [19,20]. These questions, which have started to attract a lot of attention recently, tie in with the analysis of the properties of equilibria [START_REF] Fyodorov | [END_REF][31][32].

Here we focus on the disordered Lotka-Volterra model of many interacting species, which is a representative model of well-mixed community ecology [33], and can be mapped or related to models used in evolutionary game theory and for economic systems [34][35][36][37][38]. We consider the case of symmetric interactions and small immigration and work out the phase diagram as a function of the degree of heterogeneity in the interactions and of the strength of the demographic noise. Compared to previous works [START_REF] Kessler | [END_REF]19,20,39] adding demographic noise not only allows us to obtain a more general picture, but also to fully characterize the phases and connect their properties to the ones of equilibria. In particular, we shall show that the number of stable equilibria in the LV model is exponential in the system size and their organization in configuration space follows general principles found for models of mean-field spin-glasses. Our findings, which are obtained for symmetric interactions, provide a useful starting point to analyze the non-symmetric case, as we shall demonstrate by drawing general conclusions on properties of equilibria in the case of small asymmetry.

Henceforth we focus on the disordered Lotka-Volterra model for ecological communities [START_REF] Kessler | [END_REF]19] defined by the equations:

dN i dt = N i   1 -N i - j,(j =i) α ij N j   + η i (t) (1)
where N i (t) is the relative abundance of species i at time t (i = 1, . . . S), and η i (t) is a Gaussian noise with zero mean and covariance η i (t)η j (t ) = 2T N i (t)δ ij δ(t -t ) (we follow Ito's convention). This noise term allows us to include the effect of demographic noise in a continuous setting [40][41][42]; the larger is the global population the smaller is the strength, T , of the demographic noise. Immigration from the mainland is modeled by a reflecting wall for the dynamics at N i = λ, since this is more practical for simulations than the usual way of adding a λ in the RHS of Eq. (1) (see the Appendix for more details).

The elements of the interaction matrix α ij are independent and identically distributed variables such that:

mean[α ij ] = µ/S var[α ij ] = σ 2 /S , (2)
that we consider in the symmetric case with α ij = α ji . As shown in [20], the stochastic process induced by eq. ( 1) admits an equilibrium-like stationary Boltzmann distribution:

P ({N i }) = exp - H({N i } T (3)
where

H =- i N i - N 2 i 2 + i<j α ij N i N j + + i [T ln N i + ln θ(N i -λ)] (4)
The before-last term is due to the demographic noise and the last one to the reflecting wall, which leads to a lower-immigration cut-off, at N i = λ (θ(x) is the Heaviside function). By taking advantage of this mapping to an equilibrium statistical mechanics problem and by using theoretical methods developed for disordered systems, we obtain the properties of the stationary states and the equilibria of the LV-model from the analysis of the equilibrium states and the local minima of the energy function H. Our theoretical framework is standard and based on the replica method [START_REF] Mézard | Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications[END_REF]; the computation is described in full details in the Appendix. Here, we present directly the results. Among the most important ones is the existence of three distinct phases for the LV-model in presence of demographic noise and small but non-zero immigration, as shown in Fig. 1 (we focus on N c = 10 -2 , similar results are obtained for smaller values of N c ). We find no sensitive dependence on the average interaction parameter, so the phase diagram has been obtained at fixed value µ = 10. More details will follow in the Appendix.

For large enough demographic noise (corresponding to high-temperature) we find that there is a single equilibrium phase, i.e. the noise is so strong that the interactions within species do not play an important role: for any initial condition the system relaxes toward a unique dynamically fluctuating stationary state. When the strength of the demographic noise decreases, multiple states emerge. We can study this transition by analyzing the stability of the thermodynamic high-temperature phase. This is performed by analyzing its free-energy Hessian matrix H. The point at which the lowest eigenvalue of H reaches zero signals the limit of stability of the high-temperature phase and the emergence of multiple equilibria.

Within the replica method that we used here, this corresponds to the breaking of replica symmetry and to the requirement of having a zero replicon eigenvalue. As explained in the Appendix, this leads to an equation for the Upon decreasing the noise three different phases can be detected: i) a single equilibrium phase; ii) a multiple equilibria regime between the light blue and the orange lines; iii) a Gardner phase, which turns out to be characterized by a hierarchical organization of the equilibria in the free-energy landscape.

transition line corresponding to the blue curve in Fig. 1:

λ R = (βσ) 2 1 -(βσ) 2 ( N 2 i -N i 2 ) 2 = 0 , (5)
where β = 1/T . The average • is the thermodynamics average taken over the effective Hamiltonian (4), while • denotes the average over the quenched disorder associated to the random interactions (i is a dummy index since statistically all species are equivalent after average over the interactions). Physically, the condition above can be shown to correspond to a diverging response function [19,20], and is a signature of the system being at the edge of stability, namely at a critical point in the parameter space.

Below the blue curve there exist multiple stateswhich one is reached dynamically depends on the initial condition. Such states correspond to dynamically fluctuating equilibria that are stable to perturbations and that have typically an overlap in configuration space given by

q 0 = 1 S i N i α N i β , (6)
where α and β denote the average within two generic states α and β. One can similarly define the intra-state overlap

q 1 = 1 S i N i 2
α . See Fig. 2 for a pictorial representation of these two quantities and the organization of equilibria in phase space. This is (in the replica jargon) the so-called one-step replica symmetry breaking phase (1RSB). In order to characterize the properties of this phase of the LV-model, we have computed the number of states, and hence of equilibria, using methods developed for structural glasses [START_REF] Monasson | [END_REF]. More specifically, we have computed the complexity Σ (see the Appendix for details), which is defined as the logarithm of the number of equilibria with a given free-energy density f normalized by the number of species S. This allow us to show that the number of equilibria below the blue line in Fig. 1 is exponential in S, i.e. there is a finite complexity Σ.

When decreasing further the demographic noise, the heterogeneity in the interactions becomes even more important and a second phase transition takes place. In order to locate it, we repeat exactly the same procedure as for the single equilibrium phase but now within one of the typical states with a given free-energy f [45]. The computation is more involved (it corresponds to analyze the stability of the 1RSB Ansatz) and leads to the condition:

λ 1rsb R = (βσ) 2 1 -(βσ) 2 ( N 2 1r -N 2 1r ) 2
m-r = 0 (7) where the two different averages correspond to i) the intra-state average • 1r and the inter-state average,

• m-r . All technical details of the calculation will follow in the Appendix. The critical temperature that results from the equation above leads to the orange line in Fig. 1. Crossing this line results in a fragmentation of each state into a fractal structure of sub-basins [46] (see the landscape on the bottom in Fig. 1): each state becomes a meta-basin that contains many equilibria, all of them marginally stable, i.e. poised at the edge of stability [20], and organized in configuration space in a hierarchical way, as it was discovered for mean-field spin glasses [START_REF] Mézard | Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications[END_REF]. This phase, which is called Gardner, plays an important role in the physics of jamming and amorphous materials [47]. Our results unveil its relevance in theoretical ecology by showing that it describes the organization of equilibria in the symmetric disordered LV-model at low enough demographic noise and for highly heterogeneous couplings.

We now present numerical simulation results that confirm and complement our analytical study. We numerically integrate the stochastic equation Eq. (1) using a specifically designed method (see the Appendix for details). The initial abundance for each species is drawn independently in [0, 1]. There are three sources of randomness for a given sample: the interactions, the initial conditions and the demographic noise. In the following, we obtain numerically the average correlation function defined by

E[N (t)N (t )] = 1 SN sample S i=1 N sample r=1 N r i (t)N r i (t ) (8)
where E[X] stands for the average over all those sources of randomness. If the system size is sufficiently large (S 1) as well as the sampling set, with N sample 1, it can be shown that the stochastic process converges in law [48]. We generically choose S ∼ 500 and N sample ∼ 50 and eventually verify there is no (S, N sample )-dependency at this scale. We find that in the high-temperature phase a time-translationally invariant (TTI) state is reached after a finite time-scale t wait :

∀t ≥ t > t wait E[N (t)N (t )] = C(t, t ) C(t -t ) . (9)
This convergence to a TTI regime is shown in Fig. 3. The long-time limit of C(t-t ) is the overlap between two generic configurations belonging to the single equilibrium state: the dashed line in Fig. 3 is the analytical prediction for lim t-t →∞ C(t -t ) which is in perfect agreement with the numerics. We have also checked that this agreement holds upon varying T and for other observables, the results are reported in the Appendix (see Fig. 11). From the time-dependence of C(t -t ) one can estimate the typical time-scale characterizing dynamical fluctuations within the single equilibrium phase. Formally, we define τ decorrel by the identity:

C(τ decorrel ) -C(∞) (C(0) -C(∞)) = 0.3 (10)
In Fig. 4, we plot τ decorrel as a function of (T -T 1RSB ), where T 1RSB is the critical value of T at which the single equilibrium phase becomes unstable (blue line in Fig. 1). We find that the thermodynamic instability is accompanied by a dynamical transition at which τ decorrel diverges as a power law with an exponent close to 0.5, see Fig. 4. For small demographic noise, i.e when T is below the blue line of Fig. 1, previous results on the dynamics of mean-field spin glasses [49][50][51] suggest that the LV-model should never reach an equilibrium stationary state, and instead it should display aging [START_REF] Cugliandolo | Slow Relaxations and nonequilibrium dynamics in condensed matter[END_REF][START_REF] Biroli | [END_REF]. In fact, one expects that among the very many equilibria the dynamics starting from high-temperature-like initial conditions falls in the basin of attraction of the most numerous and marginally stable equilibria, and display aging behaviour. This is indeed what we report in Fig. 5 which shows that the longer is the age of the system, t , the longer it takes to decorrelate. The landscape interpretation of this phenomenon is that the system approaches at long times a part of configuration space with many marginally stable equilibria. This leads to aging because the longer is the time, the smaller is the fraction of unstable directions to move, hence the slowing down of the dynamics, but the exploration never stops and eventually the system never settles down in any equilibrium [54][55][START_REF] Parisi | Proceedings of the National Academy of Sciences[END_REF]. The two dashed lines in Fig. 5 correspond to our analytical prediction for the intra-state and the inter-state overlaps of the marginally stable equilibria. The agreement is satisfactory but larger times would be needed to fully confirm it. Our characterization of the phases and the dynamics of the LV-model has important consequences on related systems, in particular on the so-called random replicant models (RRMs) that consist of an ensemble of replicants evolving according to random interactions. Given their numerous applications in biology, optimization problems [34,[START_REF] Mende | The predator-prey model: do we live in a Volterra world?[END_REF] as well as evolutionary game theory [START_REF] Smith | Evolution and the Theory of Games[END_REF][START_REF] Sigmund | Evolutionary Game Dynamics: American Mathematical Society Short Course[END_REF], RRMs still attract great theoretical interest. The RRM, which was introduced in [34] and further studied in [39], is remarkably similar to the disordered LV-model we studied. In the case of symmetric interactions, one can similarly map the problem onto an equilibrium statistical physics one with the following Hamiltonian:

H R = - S i<j=1 J ij x i x j -a S i=1 x 2 i ( 11
)
where x i /S is the concentration of the ith family in the species pool subject to the global constraint i x i = S for all x i ≥ 0. The couplings J ij are i.i.d. Gaussian variable with variance J 2 /S. Provided an appropriate rescaling of the interaction matrix of the two models, we can show that the average interaction term µ for LVstanding for a purely competitive environment -plays the same role as the Lagrange multiplier that is introduced in RRM to enforce the sum of all concentrations to be fixed. The main differences with respect to Eq. ( 4) is the absence of the logarithmic term. Our analysis can be fully extended to the RRM, as we show in the Appendix. The main result is that the three phases we found for the LV-model are present also for the RRM, and organized in a phase diagram (see Fig. 8) that is remarkably similar to the one in Fig. 1. This strengthens the generality of our results, and clarifies the nature of the glassy phase of the RRM that was first investigated in [39]. Let us finally discuss how we expect our results to change if the interactions contain a small random asymmetric component. The multiple basins structure associated with the 1RSB phase should not be affected because its basins correspond to stable stationary states, and a small non-conservative random force should not destabilize them [START_REF] Berthier | [END_REF]. On the contrary, the fractal structure and the decomposition into sub-basins are expected to be wiped out because of the marginal stability of the equilibria associated with it [32,[START_REF] Hertz | AIP Conference Proceedings[END_REF][START_REF] Hertz | Heidelberg colloquium on glassy dynamics[END_REF]. In absence of demographic noise, one therefore expects a single equilibrium at small σ, which is replaced by an exponential number of chaotic attractors at large σ. The demographic noise adds additional dynamical fluctuations to these multiple equilbria and eventually makes them merge in a single equilibrium, thus leading to a phase diagram similar to Fig. 1 but only with the blue line and two phases (single and multiple equilibria).

In conclusion, we have unveiled a complex and rich structure for the organization of equilibria in a central model for ecological communities. Our results, supported by dynamic simulations, highlight the relevance of multiple equilibria phases for the dynamics of many strongly interacting species. Moreover, our findings clarify the glassy nature of the equilibria previously studied in [START_REF] Kessler | [END_REF]19,20,22,[START_REF] Altieri | Jamming and Glass Transitions[END_REF]. As we have shown, our results carry out to more general contexts, in particular to models originating from evolutionary game theory. We expect that the collective dynamical behaviours -the phases -found in this work go beyond the LV-model itself and may play an important role in a variety of contexts from biology to economy, which can be modeled by high-dimensional dynamical systems with random couplings.

The evolution of the species abundances N i in the ecosystem (with i = 1, ..., S) is regulated by the following dynamical equation:

dN i dt = -N i   ∇ Ni V i (N i ) + j,(j =i) α ij N j   + N i η i (t) + λ (A1) where η i (t) is a white noise with covariance η i (t)η j (t ) = 2T δ ij δ(t -t )
and T is the temperature. In the presence of a (demographic) noise -an intrinsic population randomness due to birth, death and unpredictable interaction events -Eq. (A1) represents a generalized Langevin equation with a one-species quadratic potential V i (N i ):

V i (N i ) = -ρ i K i N i - N 2 i 2 , ( A2
)
which allows us to precisely recover the well-known Lotka-Volterra equations [START_REF] Lotka | [END_REF][START_REF] Volterra | Variazioni e fluttuazioni del numero d'individui in specie animali conviventi[END_REF]. The adimensional parameter ρ i = r i /K i denotes the ratio between the growth rate, r i , and the carrying capacity, K i . For simplicity, we will assume no species dependence on these parameters and set r i = 1, K i = 1 in the following. For a first-order analysis, the interaction matrix α ij is assumed to be symmetric: its elements are i.i.d. Gaussian variables with mean and variance respectively:

mean[α ij ] = µ/S var[α ij ] = σ 2 /S . (A3)
We can also suppose to add a small degree of asymmetry and perform a perturbative expansion. We should expect, even for arbitrary small asymmetry, that the critical (multiple equilibria) phase described in the main text will be suppressed, as proven long ago in spin-glass and neural network contexts [START_REF] Hertz | AIP Conference Proceedings[END_REF][START_REF] Hertz | Heidelberg colloquium on glassy dynamics[END_REF].

In the case of symmetric interactions, the dynamical equation (A1) admits an invariant probability distribution in terms of an Hamiltonian operator H, as can be proven by writing the Fokker-Planck equation of the corresponding stochastic process. More precisely, we can safely define an Hamiltonian H and solve the problem exactly within the replica formalism [START_REF] Mézard | Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications[END_REF], as we will show in detail in the next Section. We thus define

H = i V i (N i ) + i<j α ij N i N j + (T -λ) i ln N i (A4)
where λ denotes an infinitesimally small speciesindependent immigration rate. It essentially guarantees the existence of all invadable species. Furthermore, it has to satisfy the lower bound λ with λ > T such that the probability distribution P ({N i }) is correctly regularized at small N i .

The connection with the statics is now clear: the original dynamical process in Eq. (A1) describes the time evolution of a large interacting ecosystem whose thermodynamics is determined by the Hamiltonian (A4). It is also worth noticing that in the limit T → 0, λ → 0, Eq.(A4) reproduces precisely a spin-glass model [START_REF] Mézard | Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications[END_REF] where the continuous variables representing the abundances are then replaced by spin variables coupled by the α ij s. We will come back to this point in Appendix (C).

Replica symmetric ansatz

To deal with the thermodynamics of disordered systems and to give a precise characterization of possible phase transitions associated with the emergence of complex collective behaviours we can resort to the replica method [START_REF] Mézard | Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications[END_REF][START_REF] Parisi | [END_REF]. It was originally introduced to study spin-glass models but has now become a cornerstone for a vast class of complex systems. The introduction of replicas allows us to handle quantities in which the disorder plays a key role and eventually obtain the free energy by means of the identity

-βF = lim n→0 ln Z n n . (A5)
In principle, the free energy should depend on the specific realization of the disorder. However, in the thermodynamic limit we can safely ignore this issue since the free energy will converge to a unique value, thanks to its self-averaging property. Operatively, one should compute the quantity on the r.h.s for integer values of the replica number n, then consider the analytical continuation to real values and only at the end of the computation take the limit n → 0.

The starting point is then the computation of the replicated partition function, which in this specific case becomes:

Z n = i,(ij) dN a i dαij exp   - (ij) (αij -µ/S) 2 2σ 2 /S -βH({N a i })   (A6)
where the overline denotes the average over the disorder, i.e. the average over the Gaussian variables αij. To perform the computation -that will involve quadratic terms in the product of the species abundances -we introduce the overlap matrix Q ab (with diagonal value Qaa) and the external field Ha, where (a, b) represent two replicas of the same system, i.e.:

Q ab = 1 S S i=1 N a i N b i , (A7) Ha = 1 S S i=1 N a i . (A8)
We can thus rewrite the free energy in the replica space as:

F = - 1 βn ln a,(a<b)
dQ ab dQaadHa e SA(Q ab ,Qaa,Ha) (A9)

where the action reads

A(Q ab , Qaa, Ha) = -ρ 2 σ 2 β 2 a<b Q 2 ab 2 + -ρ 2 σ 2 β 2 a Q 2 aa 4 + ρµβ a H 2 a 2 + 1 S i ln Zi . (A10)
In turn, the partition function Zi is

Zi = a dN a i exp (-βH eff ({N a }i)) , (A11)
which depends on the effective Hamiltonian:

H eff ({N a }i) = -βρ 2 σ 2 a<b N a i N b i Q ab + -βρ 2 σ 2 a (N a i ) 2 Qaa 2 + a ρµHaN a i + + Vi(N a i ) + (T -λ) ln N a i . (A12)
The simplest scenario in the panorama of all possible replica techniques corresponds to the replica symmetric (RS) computation, which turns out to be correct as long as the freeenergy landscape is characterized by one single equilibrium state. Any permutation of the replica indices does not affect the matrix structure. In order words, within the RS Ansatz the permutation symmetry of the replicated Hamiltonian is respected.

The overlap matrix is thus parametrized by two values: the self-overlap between replicas inside the same state, q d , and the inter-state overlap, q0. The external field is assumed to be uniform, ∀a.

Q ab = q0 if a = b Qaa = q d if a = b Ha = h ∀a (A13)
The action A then becomes

A(q d , q0, h) = -ρ 2 σ 2 β 2 n(n -1) 4 q 2 0 -ρ 2 σ 2 β 2 n 4 q 2 d + + ρµβ n 2 h 2 + 1 S i ln Zi (A14)
where the partition function is integrated over N a i with an effective Hamiltonian that depends now on the parameters (q d , q0, h). Replica indices are nevertheless still coupled. At this stage, the replica trick comes into play allowing us to decouple replicas by the introduction of an auxiliary Gaussian variable z, with zero mean and unit variance, which makes the expression of the partition function of the form:

Zi = +∞ -∞ dzi √ 2π e -z 2 i /2 n a=1 dN a i e -β a H RS (N a i ,z i ) , (A15)
which is written in terms of the RS Hamiltonian:

HRS(Ni, zi) = -ρ 2 σ 2 β(q d -q0) N 2 i 2 + (ρµh -ziρ √ q0σ)Ni+ + Vi(Ni) + (T -λ) ln Ni = = N 2 i 2 ρ -ρ 2 σ 2 β(q d -q0) + + (ρµh -ziρσ √ q0 -ρ) Ni + (T -λ) ln Ni . (A16)
where we recall that Vi

(Ni) = -ρNi 1 -N i 2
and the immigration rate λ → 0 + .

In the thermodynamic limit, we can safely resort to the Laplace method and evaluate the integral by saddlepoint approximation. From the maximization of the action A(q d , q0, h), we get the corresponding for (q d , q0, h), which, after considering the analytical continuation n → 0, read:

q d = Dz ∞ Nc dN e -βH RS (q 0 ,q d ,h,z) N 2 ∞ Nc dN e -βH RS (q 0 ,q d ,h,z) = N 2 , q0 = Dz ∞ Nc dN e -βH RS (q 0 ,q d ,h,z) N ∞ Nc dN e -βH RS (q 0 ,q d ,h,z) 2 = N 2 , h = Dz ∞ Nc e -βH RS (q 0 ,q d ,h,z) N ∞ Nc dN e -βH RS (q 0 ,q d ,h,z) = N . (A17)
where the calligraphic notation stands for the Gaussian integral Dz ≡ dz √ 2π e -z 2 /2 . More precisely, the most internal average corresponds to the standard average over the Boltzmann measure, while the external one corresponds to the average over the quenched disorder. A worth noticing aspect concerns the correct choice of the extremes for the integration over N . The integral cannot be extended over the interval [0, ∞) and a more attentive analysis is needed because of the term T ln N in the Hamiltonian. The logarithmic term contributes to tilting the quadratic potential and providing a negative, divergent trend for values of the abundances very close to zero. There are two alternatives: either play with the immigration rate λ or carefully select the lower cut-off in the integral over the species abundances, Nc (for more details, see Appendices E and F). We set the immigration rate to zero from the beginning of the computation. Conversely, we probe the optimal value of the cut-off Nc, eventually set to 10 -2 , as detailed also in E 4. Lower cut-off values cannot be selected as they would result in a non-optimal matching between our theoretical analysis and numerical simulations within a Dynamical Mean-Field Theory formalism.

The solution of Eqs. (A17) is obtained by implementing an iterative algorithm. The iteration stops when the relative error between the value of the updated parameter and that at the previous time step is smaller than a given precision value, . We study the convergence of the equations and their stability as a function of the demographic noise, whose amplitude corresponds to the physical temperature T . The analysis at zero temperature has been recently performed in [20]. Our work provides then a full and more comprehensive perspective. We derive the complete phase diagram starting from a high demographic noise initial condition -corresponding to the high-temperature phase -up to a very low demographic noise (corresponding to the zero-temperature limit).

We report below the resulting values at different inverse temperatures β ≡ 1/T obtained via the numerical integration for µ = 10, σ = 1 and a cut-off Nc = 10 -2 .

Upon increasing β, the diagonal and off-diagonal value of the overlap matrix tend to become degenerate with (q d -q0) → 0 (see Table I) and the stability matrix -which is defined by the second derivative of the free energy with respect to the overlap -develops zero modes, precisely at β ≈ 58. The RS solution is then marginally stable. In this regime, to correctly incorporate and describe the complexity of the landscape, a more structured Ansatz, so-called 1RSB, must be introduced.

β q d q0 h 0.

One-step replica symmetry breaking ansatz

In the presence of low demographic noise the RS solution is no longer appropriate to describe the thermodynamics of the system, which can be more correctly identified by a one-step replica symmetry breaking (1RSB) Ansatz. This instability is intrinsically related to the emergence of multiple minima in the free-energy landscape, each of them associated with a different equilibrium configuration. According to the 1RSB approximation, the n replicas are now divided into n/m different groups of m replicas, with 0 ≤ m ≤ 1 and n/m integer. The overlap matrix assumes now three different values

Q ab =      q d if a = b q1 if a, b ∈ B l q0 if a, b ∈ B l (A18)
where B l stands for the group of replicas in the same block. This Ansatz yields the definition of an overlap matrix with m -1 off-diagonal elements equal to q1 and n -m elements equal to q0. A new parameter, m, comes into play which is usually denoted as breaking point parameter in the replica jargon. The external field is instead assumed to be uniform for all replicas a.

Ha = h , ∀a . (A19)
Exactly as in Eq. (A9), the free energy within the 1RSB Ansatz reads:

F 1RSB = - 1 βn ln dq d dq1dq0dh e SA(q d ,q 1 ,q 0 ,h) (A20)
We are now able to write the resulting expression for the action A:

A(q d , q1, q0, h) = -ρ 2 σ 2 β 2 n 4 (n -m)q 2 0 + (m -1)q 2 1 + q 2 d + + ρµβ n 2 h 2 + 1 S i ln Zi , (A21)
where the replicated partition function Zi turns out to be

Zi = dzi √ 2π e -z 2 i /2 n/m a B =1 dta B ,i √ 2π a∈a B dN a i × exp   - t 2 a B ,i 2 -β a∈a B H1RSB(N a i , zi, ta B ,i)   (A22)
hence, according to Eq. (A16), we have:

H1RSB(Ni, zi, ta B ,i) = -ρ 2 σ 2 β(q d -q1) N 2 i 2 + + (ρµh -ta B ,iρσ √ q1 -q0 -ziρσ √ q0)Ni+ + Vi(Ni) + T ln Ni . (A23)
To decouple replica indices according to a 1RSB solution, we have thus introduced a double Gaussian integration leading to two different kinds of averages: i) over the single replica; ii) over a group of replicas in the same block of size m. The two averages correspond respectively to

• 1r = ∞ Nc dN exp [-βH1RSB(N, z, ta B )] • ∞ Nc dN exp [-βH1RSB(N, z, ta B )] , (A24)
• m-r = dta B √ 2π e - t 2 a B 2 ∞ Nc dN exp [-βH1RSB(N, z, ta B )] m • dta B √ 2π e - t 2 a B 2 ∞ Nc dN exp [-βH1RSB(N, z, ta B )] m (A25)
The values of the overlap parameters q d ,q1, q0 and the external field h are obtained by saddle-point approximation of the action, leading to the following compact expressions

q d = N 2 1r m-r q1 = N 2 1r m-r q0 = N 1r 2 m-r h = N 1r m-r (A26)
that can be solved iteratively up to convergence. The initial condition is chosen close to the RS solution but with q1 > q0.

As the overlap parameters stand for the degree of similarity between two replicas in the same block or in different blocks, the found solution is thermodynamically relevant only if q d > q1 > q0. Indeed, in the very low-temperature limit, q d → q1, similarly to the analysis we have performed for the RS solution (see A 1).

The determination of the breaking parameter m requires more attention. We can follow different routes depending on whether one aims to find the m-value that optimizes the expression of free energy or the one that identifies the aging solution in a dynamical framework. In this second case, the parameter m is chosen in such a way that the found solution is marginal: one optimizes over the parameters in Eq. (A26) while selecting m accordingly to the marginal stability condition. We will devote the next Section to a detailed explanation of what marginal stability means and what kind of consequences it yields for the system.

Meaning of the replicon eigenvalue of the stability matrix

To investigate the stability of the different phases, we introduce the Hessian matrix of the free energy, which allows us to study the harmonic fluctuations with respect to δQ ab . Thanks to symmetry group properties of the replica space, the diagonalization of the stability matrix can be expressed in terms of three different sectors. Following [67], we define the three eigenvalues: the longitudinal, λL, the anomalous, λA, and the replicon, λR. We are specifically interested in the computation of the replicon mode as it is responsible for possible RSB effects. By contrast, a zero longitudinal mode can give information in terms of spinodal points describing how a state opens up along an unstable direction and originates then a saddle.

The detection of a vanishing replicon mode from the hightemperature (one single equilibrium) phase is related to the appearance of marginal states. This feature is extremely important because of its intimate connection with out-ofequilibrium aging dynamics.

We consider then the variation of the RS action with respect to the overlap matrix

A(Q ab , Qaa, Ha) = -ρ 2 σ 2 β 2 a<b Q 2 ab 2 -ρ 2 σ 2 β 2 a Q 2 aa 4 + + ρµβ a H 2 a 2 + 1 S i ln Zi , (A27)
where the partition function

Zi = a dN a i exp β 2 ρ 2 σ 2 2 a<b Q ab N a i N b i + β 2 ρ 2 σ 2 a (N a i ) 2 Qaa 2 -ρβµ a N a i H a -βVi(N a i ) -ln N a i (A28)
which implies to the second order

M abcd ≡ - ∂ 2 A ∂Q ab ∂Q cd = = β 2 ρ 2 σ 2 δ (ab),(cd) -(β 2 ρ 2 σ 2 ) N a N b , N c N d c (A29)
where the subscript • c denotes the connected part of the correlator. When evaluated at the saddle point and in the limit n → 0, the stability matrix (A29) can be decomposed as a function of three different correlators

M abcd =M ab,ab δacδ bc + δ ad δ bc 2 + +M ab,ac δac + δ bd + δ ad + δ bc 4 + +M ab,cd (A30)
from which the projection on the replicon subspace is λR = (βρσ) 2 1 -(βρσ) 2 (M ab,ab -2M ab,ac + M ab,cd ) (A31) The • is performed over the effective Hamiltonian, while the average over the quenched disorder is always denoted as •. The replicon mode gives information about the fluctuations inside the innermost block of the overlap matrix, corresponding to the fluctuations within the same state. According to this structure and the different combinations of the replica indices, three elements must be determined

M ab,ab -2M ab,ac + M ab,cd = (N a ) 2 (N b ) 2 -2 (N a ) 2 N b N c + N a N b N c N d (A32) λ R T Figure 6
. Replicon eigenvalue obtained within the replica symmetric (RS) Ansatz for µ = 10, σ = 1 (in orange) and σ = 0.85 (in green) as a function of temperature. Below the dashed black line corresponding to the zero value, the replicon becomes negative and makes the RS approximation no longer valid. Upon decreasing σ, the transition towards the 1RSB phase occurs at a lower critical temperature.

In the simplest scenario corresponding to the presence of one single equilibrium, the expression for the replicon eigenvalue can be further simplified as:

λR = (βρσ) 2 1 -(βρσ) 2 ( N 2 -N 2 ) 2 , ( A33
)
where the averaged difference describes the fluctuations between the first and second moment of the species abundances within one state, namely between the diagonal value q d and the off-diagonal contribution q0 of the overlap matrix. The difference between the two overlap values can also be interpreted as the response function of the single species to an infinitesimal perturbation. A large, diverging response -which is related to the inverse of the replicon mode -is a further evidence that the system is close to a critical point. We can repeat exactly the same procedure to obtain the replicon in the 1RSB Ansatz. The analysis of harmonic fluctuations proceeds exactly as for the RS case. The only difference now is that -as for Eqs. (A24)-(A25) -we have to distinguish between the intra-state average and the inter-state average, that is over a bunch of replicas in the inner block of size m.

λ 1rsb R = (βρσ) 2 1 -(βρσ) 2 ( N 2 1r -N 2 1r ) 2 m-r , (A34)
The value at which the 1RSB replicon vanishes signals a critical phase transition towards a more structured phase. This so-called Gardner phase turns out to be characterized by a hierarchical organization of the different equilibria, whose analytical solution is well described within a Full RSB Ansatz.

For the first time in an ecological context we are able to prove in a rigorous way the analogy between glassy physics and complex energy landscapes in large ecosystems, with respect to previous studies in this direction [20,22].

Appendix B: Derivation of the complete phase diagram

Our analysis, based on increasingly structured Ansatz, has allowed us to derive the complete phase diagram of the Lotka-Volterra model as a function of the demographic noise. By varying this control parameter, we have highlighted the existence of different phase transitions of increasing complexity. Thus far the only available results have been obtained without demographic noise, exactly at zero temperature [19]. ). We can distinguish three different phases: i) a single equilibrium phase above the blue line; ii) a multiple equilibria regime (1RSB stable phase) between the light blue and the orange lines; iii) a Gardner phase, characterized by a hierarchical organization of the different equilibria in the free-energy landscape.

In correspondence of the blue line in Fig. (7) the landscape structure is no longer identified by a single equilibrium but should be replaced by a two-level hierarchy, which results in the appearance of new equilibrium configurations.

In particular we have realized that with very low demographic noise the system undergoes a Gardner transition to a new marginally stable phase: each amorphous state, say a basin, is fragmented into a fractal structure of sub-basins (metabasin). The internal structure of states in which a basin splits is described by the Full RSB solution of the partition function, according to which the replica symmetry is broken an infinite number of times. Similar to what observed in lowtemperature glasses, such a transition will revolutionize our understanding of ecosystems. The emergence of a infinite hierarchy of scales for the order parameter is also related to an infinite number of time-scales playing a central role in population dynamics, as we will discuss in the next Sections.

What a Full RSB regime precisely yields for biological and ecological settings with the appearance of an abundance of soft modes is particularly fascinating and still an open matter of debate.

Appendix C: Results for the spin-glass model without logarithmic interaction

We consider exactly the same system without the logarithmic interaction in the species abundances. In the limit λ → 0 + , the resulting Hamiltonian is

H = i Vi(Ni) + ij αijNiNj (C1)
where the first term has the same quadratic dependence in the species abundances as in Eq. (A2), while the latter represents the pairwise interacting part of the Hamiltonian. In theoretical physics, Eq. (C1) is an example of spin-glass model where the degrees of freedom can be either discrete variables, to model magnetic spins, or continuous variables. The couplings between spins are chosen randomly, taking both positive and negative values and thus resulting in frustration and nonconvex optimization phenomena. The simplest and most classical example of spin glass is represented by the Sherrington-Kirkpatrick (SK) model [68] where each spin variable interacts with all the others in a fully-connected topology.

Connections with the random replicant model

In theoretical ecology, a similar model was first introduced by Diederich and Opper in the 80s, studying the mean-field dynamical evolution of S randomly interacting species with a deterministic self-interaction [34]. This pioneering model is usually known as replicator model. Later on, a detailed study was also performed in [39] by using the replica formalism in the zero-temperature limit. Even though without the demographic noise, this model turns out to be not particularly interesting in a purely ecological context, it can still provide further insights on the replicator equations in many other interdisciplinary domains [69], such as in the study of Nash equilibria [70] in Game Theory.

Given S species, the Lyapunov function of the replicator model is written as a function of the concentration variables xi in the following way

HR = S i<j=1 Jijxixj -a S i=1 x 2 i (C2)
where the control parameter a limits the growth rate of one single species. The global constraint contributes to nonconvex constraint optimization properties with non-trivial symmetry breaking effects.

The couplings Jij are extracted from a Gaussian distribution with zero mean and variance J 2 , i.e.:

P (Jij) = S πJ 2 exp - SJ 2 ij J 2 (C3)
playing the same role as the αijs in our notation. xi are real variables subject to the global constraint i xi = S, ∀xi ≥ 0 with i = 1, ..., S. One can introduce a Lagrange multiplier γ to enforce the normalization condition over xi, which leads to the following expression for the partition function:

Z = S 0 S i=1 dxi i∞ -i∞ dγ exp -β ij Jijxixj -βa i x 2 i + -γ i (xi -1) . (C4)
As usual, one can introduce the replica trick to average over the quenched variables Jij and thus rewrite Eq. (C4) as

Zn = dγαTrnS exp -S α βaX 2 α + γα(Xα -1) + + S 2 β 2 J 2 4 αβ X 2 α X 2 β (C5)
where X are exactly the same n variables as before that have now lost their dependence on the index i because of the equivalence of all sites. By introducing the overlap parameter Q αβ between Xα and X β , one obtains

Zn = dγαdQ αβ exp   -S αβ Q 2 αβ + S α γα   × × TrnS exp [SL(Q, γ, X)] (C6)
where the action L(Q, γ, X) is

L(Q, γ, X) ≡ βJ αβ Q αβ XαX β -βa α X 2 α - α γαXα .
(C7) The second term, absent in the ordinary SK model, is actually due to the global quadratic constraint on the species concentrations. Using the same Ansatz as in [39] and optimizing over Q αβ and γα with the following choice:

Q αβ = qδ αβ + t γα = γ (C8)
one can get to writing the replica symmetric action

LRS(q, t, γ, X) = -β(a -qJ) α X 2 α + 2z βJt -γ α Xα (C9)
where the auxiliary Gaussian variable z has been introduced to decouple replicas. Because under this transformation the integrals become independent, one can safely forget the dependence on the replica index α and re-scale the expression above by J LRS(q, t, γ, X) = -βJ (a/J -q) X 2 -2z t -γ X (C10) According to our notation and in particular to Eqs. (A13)-(A16) considered without the logarithmic term and the immigration parameter, a direct mapping between the two models can be immediately pointed out. We additionally consider a quadratic potential quadratic potential V (X) = -ρX 1 -X 2 with respect to the original model proposed in [39], where we have replaced N ↔ X for consistency of notation. We can thus conclude

ρ -ρ 2 σ 2 β(q d -q0) ⇐⇒ -βJ(a/J -q) (ρµh -zρσ √ q0 -ρ) ⇐⇒ -2z t -γ . (C11)
The Lotka-Volterra model embeds a quadratic one species potential as a function of N 2 i , which contributes to shifting the linear and quadratic terms of an amount ρ. The parameter γ in the random replicant model is related to the (competitive) interaction term µ, while ã ≡ a/J is associated with the heterogeneity parameter σ.

In [39] it was shown there exists a critical value of the control parameter ã =≡ a J , ãc = 1/ √ 2, below which the system develops non-ergodic features and replica symmetry breaking effects. At higher values, irrespective of the initial conditions, only one equilibrium is possible, whereas in the region characterized by ã < ãc the final configuration is strongly affected by a small perturbation that can allow the system to reach different equilibria. This last phase corresponds to a high-competition scenario between the interacting species. Interestingly, the replicator model can be proven to have a one-to-one mapping with many different problems of interest in optimization and disordered systems.

Exact solution in the replica symmetric case

The following analysis focuses on the Lotka-Volterra Hamiltonian, as reported in Eq. (A4), and considered in the absence of the logarithmic term. We will derive exact expressions for the order parameters to be analyzed as a function of temperature up to T = 0. In this case, the temperature explicitly comes out in the weighting factor exp(-βH) and implies for the replicated partition function:

Z n = i,(ij) dN a i dαij exp   - (ij) (αij -µ/S) 2 2σ 2 /S -βH({Ni})   (C12)
Exactly as for the previous analysis, we have performed both a RS and 1RSB computation. In the former, the saddlepoint equations read:

q d = Dz ∞ 0 dN e -βH RS (q 0 ,q d ,h,z) N 2 ∞
0 dN e -βH RS (q 0 ,q d ,h,z) , q0 = Dz ∞ 0 dN e -βH RS (q 0 ,q d ,h,z) N ∞ 0 dN e -βH RS (q 0 ,q d ,h,z) 2 , h = Dz ∞ 0 e -βH RS (q 0 ,q d ,h,z) N ∞ 0 dN e -βH RS (q 0 ,q d ,h,z) .

(C13)

where now the integral over N is extended over the entire positive axis without fixing any cut-off. Dz denotes, as usual, the Gaussian integration over the auxiliary variable z -which is introduced to decouple replicas -with zero mean and unit variance. The RS Hamiltonian then reads:

HRS(N, z) = -ρ 2 σ 2 β(q d -q0) N 2 2 + (ρµh -zρ √ q0σ)N + V (N ) = N 2 ρ -ρ 2 σ 2 β(q d -q0) + (ρµh -zρσ √ q0 -ρK) N . (C14)
Therefore, Eqs. (C13) can be exactly computed in terms of the error function and its combinations, i.e.

q d = Dz e - βρ(K-hµ+ √ q 0 zσ) 2 2-2(q d -q 0 )βρσ 2 2(K -hµ + √ q0zσ) βρ[1 -(q d -q0)βρσ 2 ]+ e βρ(K-hµ+ √ q 0 zσ) 2 2-2(q d -q 0 )βρσ 2 √ 2π 1 + βρ (K -hµ) 2 + 2 √ q0z(K -hµ)σ + (q0 -q d + q0z 2 )σ 2 × × 1 + Erf βρ(K -hµ + √ q0zσ) √ 2 βρ(1 + q0βρσ 2 -q d βρσ 2 ) 1 √ 2Z , (C15) q0 = Dz e - βρ(K-hµ+ √ q 0 zσ) 2 2-2(q d -q 0 )βρσ 2 e βρ(K-hµ+ √ q 0 zσ) 2 2-2(q d -q 0 )βρσ 2 √ πβρ(K -hµ + √ q0zσ) + 2βρ(1 + q0βρσ 2 -q d βρσ 2 )+ + e βρ(K-hµ+ √ q 0 zσ) 2 2-2(q d -q 0 )βρσ 2 √ πβρ (K -hµ + √ q0zσ) Erf βρ(K -hµ + √ q0zσ) √ 2 βρ(1 + q0βρσ 2 -q d βρσ 2 ) 1 Z 2 , (C16) h = Dz √ πβρ(K -hµ + √ q0zσ) + √ 2e - βρ(K-hµ+ √ q 0 zσ) 2 2-2(q d -q 0 )βρσ 2 βρ(1 + q0βρσ 2 -q d βρσ 2 )+ + √ πβρ(K -hµ + √ q0zσ)Erf βρ(K -hµ + √ q0zσ) √ 2 βρ(1 + q0βρσ 2 -q d βρσ 2 ) 1 Z , (C17)
where the partition function reads:

Z = √ πβρ[1 -(q d -q0)βρσ 2 ] 2 1 + Erf βρ(K -hµ + √ q0zσ) √ 2 βρ(1 + q0βρσ 2 -q d βρσ 2 ) . ( C18
)
Exactly as for the previous Section where we have derived the expression for the replicon eigenvalue, we have here:

λR = (βρσ) 2 1 -(βρσ) 2 ( N 2 -N 2 ) 2 , ( C19
)
where the expressions for the averaged values N 2 and N 2 correspond precisely to q d and q0 obtained in Eqs. (C15)-(C16) according to a RS Ansatz. The overline denotes the average over the quenched disorder, which is technically implemented by integrating over the Gaussian variable z.

Replica symmetry broken case

The computation can be easily generalized to more complex scenarios. In particular, in the 1RSB case the Hamiltonian reads

H1RSB(N, z, t) = N 2 2 ρ -ρ 2 σ 2 β(q d -q1) + +(ρµh -tρσ √ q1 -q0 -zρσ √ q0 -ρK)N (C20)
with the only difference of considering now the double Gaussian integration over z and t that account respectively for the single replica average and the m-block average. Even though the structure of the resulting equations for (q d , q1, q0, h), we are still able to get an exact closed-form derivation. Exactly as before for the most general case, the saddle-point equations read

q d = N 2 1r m-r q1 = N 2 1r m-r q0 = N 1r 2 m-r h = N 1r m-r (C21)
to be solved iteratively.

Derivation of the phase diagram: evidence of a Gardner phase

By fixing the carrying capacity K = 1 and ρ = 1 in Eqs. (C15)-(C16)-(C17), we have explored the phase space as a function of the parameter β = 1/T along with the variation of the mean and the variance µ and σ respectively of the interaction matrix. The equations (C13) have been solved iteratively starting from the high-temperature region up to zero temperature. For high values of σ and T there exists only one equilibrium: the system is thus able to recover its equilibrium configuration even starting from different initial conditions. As shown in Fig. (8), upon decreasing T new phases emerge, leading respectively to an intermediate stable 1RSB phase and a Gardner phase in the very low-temperature limit.

We have thus obtained the complete phase diagram at fixed µ = 10. In Eqs. (C11) we have derived the exact mapping between our model and the Random Replicant Model (RRM), pointing out how the mean and the variance of the interaction matrix are essentially related to the Lagrange multiplier γ and the parameter ã ensuring the global quadratic constraint. In our case the phase diagram does not display any special dependence on the mean interaction µ and can be thus analyzed at fixed value, by changing only the heterogeneity and the demographic noise. However, at variance with the original RRM model, fixing µ would correspond to allowing the sum of the species concentrations to vary, i.e. injecting or removing species in the ecosystem.

Appendix D: Derivation of the complexity

We aim to determine the number of local minima N of the free energy corresponding to the different equilibrium configurations in the landscape structure. The logarithm of the number of minima divided by S defines the complexity Σ(f ) (a.k.a. configurational entropy) of the ecosystem. This provides a crucial information to determine what kind of universality class the system belongs to. In particular, in statistical physics of disordered systems two main universality classes are well-known: i) spin-glass models, in which the number of minima of the free energy is not exponential in the system size and the free energy barriers are sub-extensive; ii) structural glasses, for which the number of free-energy minima is actually exponential in the system size. Then, there exists a finite temperature range, between a static transition temperature Ts and a dynamical transition temperature T d , within which the complexity turns out to be finite.

As we deal with mean-field models, we can resort to a simplified description for which the total partition function of the model can be expressed as the sum of α pure states contributions. More precisely, we can define a generic effective potential as a function of local order parameters whose minima are in correspondence one-to-one with the so-called Thouless-Anderson-Palmer states [START_REF] Mézard | Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications[END_REF][71][72][73]. Provided the number of minima of the free energy is exponential in the system size, the partition function can be expressed as:

Z ∼ α e -βN fα(T ) = df α δ(f -fα(T ))e -βN f = df ρ(f )e -βN f ∼ e N [Σ(f * ,T )-βf * ] n . (D1)
Therefore, the configurational entropy corresponds to

Σ(f, T ) ≡ 1 S ln α δ(f -fα(T )) (D2)
while f * ∈ [f0, f th ] satisfies the extremal condition:

d df [f -(1/β)Σ(f, T )] = 0 (D3) namely dΣ(f, T ) df = 1 T . (D4)
Three different regimes can be generically observed in meanfield models:

• i) in the high-temperature phase (corresponding to high demographic noise), namely for above the dynamical transition temperature, the paramagnetic state dominates the free-energy density for any value [f0, f th ];

• ii) between the dynamical and the statical transitions, there exists a value f * such that the quantity f * -(1/β)Σ(f * ) evaluated in f * coincides with the paramagnetic value of the free-energy density. However, in this second regime, the resulting state is composed by an exponential number of metastables states of individual free-energy density f * . Upon crossing the dynamical temperature, the free energy preserves its analyticity without undergoing any true phase transition;

• iii) For temperatures lower the static transition temperature, the leading contribution is due to the lowest free-energy states with f * = f0 and, as a consequence, Σ(f0) = 0. In this phase, the number of states is subexponential in the system size.

While in the very high-temperature phase, the paramagnetic solution is always present, for T < T d , it disappears and is replaced by a non-trivial combination of states. To explicitly compute the complexity of the system, Σ, and grasping the physics behind it, several techniques have been proposed in the last decades. In the following, we will focus on the so-called real replica method [START_REF] Monasson | [END_REF]. It consists in replicating m times the system and coupling the different copies through an infinitesimally small parameter , which will be sent to zero at the end of the computation after the thermodynamic limit. This attractive coupling naturally breaks the replica symmetry: it constrains the m copies to be in the same metastable state yet remaining uncorrelated within a state. As a consequence, the free energy can be simply written as m times the individual contribution fα. In this case:

Zm ≡ α e -N mβfα(T ) = df e N [Σ(f,T )-βmf ] (D5)
and, by evaluating the integral by saddle-point method, we get a similar expression to (D4)

∂Σ(f, T ) ∂f f * (m,T ) = βm (D6)
where one can immediately notice that the breaking parameter m has a clear counterpart in the study of the complexity.

Fixing the temperature and varying only the parameter m, one can expand the complexity around its minimum and obtain:

Σ(f, T ) ≈ Σ(f0) + a(T )(f -f0) + ... (D7)
which implies that for m = 1 the static transition is the solution of the equation a(T ) = 1 T , whereas at small m the condition is precisely replaced by Eq. (D6).

The free-energy density of a system of m different copies can thus be rewritten as [74] φ(m, β) = -

1 βN ln Zm = min f [βmf -Σ(f )] = =βmf * (m, β) -Σ(f * (m, T )) . (D8)
Generalizing the canonical definition of entropy applied to disordered systems, the complexity is thus defined as the Legendre transform of the free energy averaged over quenched disorder.

From the parametric plot of f * (m, β) and its transform Σ(m, β), one can extract the information about the behavior of the configurational entropy at any temperature and, consequently, of the associated TAP states. Therefore, Σ can be explicitly calculated from φ(m, β) or, thanks to the identity:

φ(m, β) = mF 1RSB (D9)
from a direct computation of the 1RSB free energy

Σ =m 2 d dm βF 1RSB = = -m 2 d dm 1 n
ln dq d dq1dq0dh e SA(q d ,q 1 ,q 0 ,h) .

(D10)

We compute then the derivative of the free energy w.r.t m, which eventually leads to

Σ = m 2 ρ 2 σ 2 β 2 4 (q 2 1 -q 2 0 ) + Dz ln dta B √ 2π e - ta B 2 A(z, ta B ) m -m Dz dta B √ 2π e - ta B 2 A(z, ta B ) m ln A(z, ta B ) dta B √ 2π A(z, ta B ) m (D11)
where we have denoted as A(z, ta B ):

A(z, ta B ) ≡ dN e -βH 1RSB (N,z,ta B ) (D12)

We quantitatively evaluate the above expression within the 1RSB phase to determine the nature of the emerging transition. We find strictly positive values of complexity at finite temperature and fully compatible with the results previously obtained at zero temperature [20]. We manage then to prove that the emergent 1RSB stable phase is actually characterized by a finite complexity, i.e. an exponential number of metastable states. Note also that, because the replicon 1RSB becomes negative below the Gardner transition temperature, the complexity is well-defined only in the region for which the 1RSB can be safely applied, i.e. in the interval [f0, fG], fG being the free-energy density at the Gardner transition. Within this formalism, it is also possible to reproduce the complexity curves at different and fixed m starting from its equilibrium value m * . The self-condition equation ∂F 1RSB /∂m = 0 gives indeed the equilibrium value m * , which corresponds to the lowest free energy density and then to a vanishing complexity. The resulting complexity curve is expected to have an increasing trend for lower values of m up to a maximum point where unstable states start to appear and dominate the thermodynamics.

Appendix E: Comparison theory and numerics

Protocol

For comparing with the theoretical results, we sample the dynamical system presented in Equation [START_REF] May | Theoretical ecology: principles and applications[END_REF]. The input parameters of a sample are the system size S, the interaction matrix parameters (µ, σ), the immigration λ, the strength of the demographic noise T (temperature), and the initial condition distribution P[{Ni(0)}i=1..S]. In order to sample one realization of the ecosystem, we perform the following steps:

1. We sample the S-sized symmetric interaction matrix α, from the Gaussian distribution with scaled parameters (µ, σ).

2. We sample the initial conditions Ni(t = 0) from the distribution P[{Ni(0)}i=1..S]. For instance, we use a factorized uniform distribution in [0, 1]:

P[{Ni(0)}i=1..S] = S i=1 1{Ni(0) ∈ [0, 1]}
where 1{.} is the indicator function.

3. We sample the demographic noise {ηi(t)} t=0..tmax i=1..S

, from the white-noise distribution, with temperature T . 4. Then, all three random contributions (interactions, initial conditions and demographic noise) have been dealt with. We can then integrate deterministically the system, to end up with {Ni(t)} t=0..tmax i=1..S

, where tmax is the temporal extent for the simulation. Actually, the above 3 and 4 points are a bit more involved: the implementation of immigration is detailed in Appendix F, and the exact numerical scheme we used is presented in Appendix G. But for simplicity's sake, let's focus on this framework: we fix parameters (S, µ, σ, λ, T ), we sample the three random contributions, we integrate, and we obtain the species populations over time {Ni(t)} t=0..tmax i=1..S . When we reproduce different sets of data by keeping the same parameters, but sampling different randomness, we obtain {N r i (t)} t=0..tmax i=1..S, r=1..N sample .

Observables

In order to compare with the theory, we need to decide on the observables. So far, there are four sources of statistics in the process: the three random parts (interactions, initial conditions and demographic noise) that we labelled with r = 1..N sample , and the species themselves i = 1..S. In the following, we will denote E[X] the average over all those contributions. For example:

E[N (t)N (t )] = S -1 N -1 sample S i=1 N sample r=1 N r i (t)N r i (t )
It can be shown [48] that if the system is large enough (S 1) and the sampling thorough enough (N sample 1), there is a convergence in law of the process. Mainly, there is a well defined limit (S, N sample → ∞) that we can compare with the theory. To fix ideas, we generically use S ∼ 500 and N sample ∼ 50, and we checked there is no (S, N sample ) dependency at this scale. More precisely, in the S → ∞ limit, the free energy is self-averaging, so results should typically not depend on the realization of the sampling. Here for the numerics, as 1 S < ∞, we still use some averaging over the samples to get cleaner data.

The theory is a thermodynamical one, so we will first assume that if we wait for a big-enough twait, the system will reach a time-translationnally invariant (TTI) state. For instance, the two-time correlation C is a function of the time difference:

∀ t ≥ t > twait, E[N (t)N (t )] = C(t, t ) = C(t -t )
We check this numerically. The waiting-time depends on parameters, mainly (σ, T ). However, if we lie in the replica-symmetric (RS) phase, we can always find the TTI state, for rather small waiting times twait ∼ 10 2 .

All the comparisons we will be making are in this state (t ≥ twait), for the RS phase. We will now use a mapping between thermodynamics properties, and dynamical ones. We use the notations from Equations (A13).

h = E [N (t)] q d = C(0) = E N (t) 2 q0 = lim τ →∞ C(τ ) = lim τ →∞ E [N (t)N (t + τ )] ∼ E [N (t)N (tmax)]
The lhs is predicted by the theory, and the rhs are numerical observables.

Example of numerical results in the RS phase

On figure 9, we show that one time observables such as E [N (t)] or E N (t) 2 converge to a constant value in time. This indicates the reach of a TTI state. It can be confirmed by the collapse of two-time observables such as the correlation E[N (t)N (t )] = C(t, t ), that we plot as C(t -t , t ) for different t on figure 10.

We can see that for t > twait ∼ 20 here, the system is indeed TTI, at least regarding these observables. We then read the values of h = E [N (t)] T T I and q d = E N (t) 2 T T I when they stabilize. And we read q0 = E [N (t)N (tmax)] T T I on the collapse of figure 10.

We can also infer a relevant information from figure 10: the timescale for decorrelation τ decorrel . We are only interested in the scaling of this observable, so we introduce a rough estimate. We approximate τ decorrel by the needed time so that the decorrelation decay is of 70%. Mathematically, τ decorrel is then determined by:

C(τ decorrel ) -C(∞) = 0.3 (C(0) -C(∞))

Match in the RS phase

The results are presented on figure 11: the theory matches beautifully the numerics. The dashed line correspond to the read TTI value of q0. Here, the timescale for decorrelation is around τ decorrel ∼ 10.

On figure 12, we show how the timescale for decorrelation τ decorrel diverges as we approach the 1RSB transition from above in temperature. Data seems to indicate a critical exponent as τ decorrel (T ) ∼ (T -T1RSB) -1/2 .

Rough results in the 1RSB phase

In the 1RSB phase, thermodynamics indicate that the system no longer reaches a TTI state. Instead, it presents aging behaviour: the older the system is, the slower it becomes. In Parameters are (S, µ, σ, λ) = (500, 10, 1, 10 -2 ). We consider the observables (h, qD, q0) as a function of temperature. The orange full line is the theory predictions. Blue crosses are numerical results, error bars are taken with respects to the N sample = 50 different samples of the ecosystem. We found twait ∼ 200 to be enough to observe TTI state in all these values of temperature, except for the last point on the left (T = 2.10 2 ): due to slowing down of the dynamics, we had to increase the extent of the simulation and found twait ∼ 3000. the orange dashed line correspond to the critical temperature at which the theory becomes 1RSB. Indeed, numerically we can't observe TTI state below this temperature, even increasing tmax to 10 7 . the simplest case of aging, there is a good understanding of the correlation decay C(t, t ). At equal time, the correlation is the dynamical one C(t , t ) = q d . Then it decorrelates quickly for t > t to an intermediate plateau C(t, t ) = q1 < q d , as the system explores the neighbouring phase space. Eventually, for t t , the correlation decreases to a final plateau C(t, t ) = q0 < q1. The timescale from the intermediate plateau to the final plateau increases with the age t of the system. These theoretical predictions are compared with numerical results on figure 13.

Appendix F: Mathematical issues for immigration implementation

In this part, we detail the mathematical issue for immigration implementation. This problem is independent on the interactions, so we drop them (α = 0 here). We consider the following one-species Ito-stochastic process:

dN dt = N (1 -N ) + √ 2T N η + λI(N ) (F1)
with white noise η(t)η(t ) = δ(t -t ), and immigration Red dashed line is a simple fit. As we approach the transition, the system becomes slower and slower, and the dynamical timescale diverges. Parameters are (S, µ, σ, λ, T, N sample ) = (2000, 10, 1, 10 -2 , 1/80, 40). We plot the rescaled correlation C(t, t )/C(t , t ) as a function of t -t , for different t . The different curves no longer collapse, there is no TTI state anymore. In dotted black and red lines, we respectively show the predictions for the intermerdiate and final plateau values. They do not coincide exactly with the data, but the trends correspond. function I(N ). Immigration is generically implemented so that the populations do not go too close to 0. In the usual immigration, I(N ) = 1, but we will see that this is problematic.

We want to have a hint at the stationary distribution of population P∞(N ) = P (N, t = ∞). In order to obtain it, we change variables so that the noise becomes additive, and not multiplicative any more. Here the relevant change of variables is s(t) = √ N (t), and Ito's lemma gives:

ds dt = s 2 -s 4 + λI(N ) 2s - T /4 s + 2T /4 η
Then we use Langevin-Boltzmann to read the stationary distribution.

In the usual immigration I(N ) = 1 case, the stationary distribution is always integrable:

P∞(N ) = Z -1 N λ/T -1 exp 1 T (N - N 2
2 )

However, the corresponding effective potential V ef f (N ) behaves repulsively around N = 0 only if λ > T :

V ef f (N ) = -N + N 2 2 -(λ -T ) ln N
Indeed, if we introduce an approximate induced cut-off value Ncut(b) such that P∞(N < Ncut) ∼ e -b 1, the scaling yields Ncut(b) ∼ e -bT /λ , which means that the density is still relevant up to e -bT /λ 1. Basically, this means that demographic noise with usual immigration will not prevent populations from reaching very low values. The usual immigration is not strong enough when facing demographic noise. This is indeed problematic, because whenever we will want to actually compute observables, the integrals will be dominated by the domain N ∼ 0. This is wrong physically (important species should be the high population ones), and difficult numerically (integration is illdefined).

In order to solve this, we can implement stronger immigration such as I(N ) = N -α with α > 0. However, another even simpler physical solution is to impose a hard repulsive boundary condition on the problem: an infinite potential at N = λ. In this case, the same steps can be performed and we obtain the stationary distribution:

P∞(N ) = Z -1 N -1 exp 1 T (N - N 2 2 ) 1{N > λ}
which is well-behaved. This is the solution we chose for both the theory predictions and the numerics. We will now detail in the next section how to integrate this process numerically.

Appendix G: Numerical scheme to sample demographic noise 1. Litterature review Numerical simulations need discrete time. However, when discretizing time with bounded random processes, one often encounters a non-zero probability that during one time-step the system will cross the boundary of the system (for example the N ≥ 0 boundary in our case), and become numerically unstable. We review different solutions that have been proposed to solve this issue, and check how they deal with our Lotka-Volterra (LV) system. A more thorough review can be found in [42].

A first naive way to go around the difficulty is to change variable (sqrt, ln...). But this won't work because if the noise becomes treatable, the deterministic part becomes numerically unstable. Most articles then study the numerical integration of processes such as Ṅ = α + βN + √ σN η.

[75] proposes Balanced Implicit Method: they implement a clever discretization scheme so that the boundaries (positivity) are respected. The scheme amounts to Euler's for small time step. It needs a small regularization. It does not work for LV, because it needs very small regularization parameter and time-step to give good results. This is too heavy numerically. [76] derives the exact Fokker Planck solution of a simpler system. But sampling is inefficient (rejection method). [77] builds on this method by improving the sampling method, but this still isn't satisfactory. Eventually, [78] improves again the method, by exactly solving (Fokker-Planck) the full process. The sampling is clever, with Poisson variables. They also indicate a way to solve more elaborate processes, which we will detail in the following section. Our strategy is heavily based on [78].

Our implementation

The idea from [78] is to separate the process into solvable ones. More precisely, we want to solve:

Ṅi = √ Niηi -N 2 i -Ni αijNj -1 .
..+hardW all(λ) (G1) where hardW all(λ) implements the hard wall boundary at N = λ. We will discretize time with a timestep dt, and further subdivise it into three timesteps dt = dt/3. We consider that only one part of the process is active during a subtimestep dt . So the final scheme is the following:

1. From [76], we know how to sample efficiently the demo- There are a lot of different combinations of this kind of schemes. We tried some, and chose this one after a lot of checks on simpler models for which we know the distributions at all times.

Figure 1 .

 1 Figure 1. Phase diagram showing the strength of the demographic noise, T , as a function of the degree of heterogeneity, σ, at fixed µ = 10 and selected value of the cutoff, Nc = 10 -2 .Upon decreasing the noise three different phases can be detected: i) a single equilibrium phase; ii) a multiple equilibria regime between the light blue and the orange lines; iii) a Gardner phase, which turns out to be characterized by a hierarchical organization of the equilibria in the free-energy landscape.

q 1 q 0 Figure 2 .

 02 Figure 2. Zoom on a pictorial landscape. The parameters q0 and q1 denote the size of the largest and the innermost basins respectively within the two-level structure of the 1RSB phase.

Figure 4 .

 4 Figure 4. Decorrelation time as a function of (T -T1RSB) in logarithmic scale. The blue points correspond to numerical data, while the dashed red line is a fit. The decay of the decorrelation time in (T -T1RSB) occurs with an exponent ≈ -0.5.

Figure 7 .

 7 Figure 7. Two-dimensional phase diagram of the Lotka-Volterra model in the presence of demographic noise as a function of the heterogeneity parameter σ (in the presence of the logarithmic term and with a selected value of the cutoff Nc = 10 -2). We can distinguish three different phases: i) a single equilibrium phase above the blue line; ii) a multiple equilibria regime (1RSB stable phase) between the light blue and the orange lines; iii) a Gardner phase, characterized by a hierarchical organization of the different equilibria in the free-energy landscape.

Figure 8 .

 8 Figure 8. Two-dimensional phase diagram showing the instability lines respectively for the one equilibrium (RS) phase and the multiple equilibria phase (1RSB) as a function of σ (without the logarithmic term). Between the two lines (in light blue and orange respectively) a stable 1RSB phase persists. The red dot corresponds to the analytically predicted value in the zero-temperature limit allowing for the estimation of the critical value σc = 1/ √ 2 ≈ 0.707.

Figure 9 .

 9 Figure 9. RS one time observables converge in time. Parameters are (S, µ, σ, λ, T ) = (500, 10, 1, 10 -2 , 10 -1 ). This data comes from only one sample of the process, with discrete timestep dt = 10 -1 . The dashed lines correspond to the read TTI value of h and q d .

Figure 11 .

 11 Figure 11. Comparison with the theory in the RS phase. Parameters are (S, µ, σ, λ) = (500, 10, 1, 10 -2 ). We consider the observables (h, qD, q0) as a function of temperature. The orange full line is the theory predictions. Blue crosses are numerical results, error bars are taken with respects to the N sample = 50 different samples of the ecosystem. We found twait ∼ 200 to be enough to observe TTI state in all these values of temperature, except for the last point on the left (T = 2.10 2 ): due to slowing down of the dynamics, we had to increase the extent of the simulation and found twait ∼ 3000. the orange dashed line correspond to the critical temperature at which the theory becomes 1RSB. Indeed, numerically we can't observe TTI state below this temperature, even increasing tmax to 10 7 .

Figure 12 .

 12 Figure 12. Critical slowing down of the dynamics. We plot the decorrelation time τ decorrel (T ) as a function of T -T1RSB, in a loglog scale. Blues crosses come from the same numerical data as figure 11. Red dashed line is a simple fit. As we approach the transition, the system becomes slower and slower, and the dynamical timescale diverges.

  graphic noise only Ñi(t + dt ) = Gamma Poisson[ Ni(t) T dt ] T dt This corresponds to a process which only feels the demographic noise Ṅi = √ Ni ηi during [0, dt ]. We respectively used the notation Poisson[ω] (Gamma[ω]) for random Poisson (Gamma) variables, with parameter ω. 2. Treating immigration as a reflecting wall. The particle wishes to go to Ñi but bounces on the wall. Ni(t + dt ) = λ + | Ñi(t + dt ) -λ| 3. During [dt , 2dt ], only integrate the blue process Ṅi = -N 2 i : Ni(t + 2dt ) = Ni(t + dt ) 1 + dt Ni(t + dt ) 4. During [2dt , 3dt ], only integrate the pink process Ṅi = -Ni ( αijNj -1): Ni(t + 3dt ) = Ni(t + 2dt ) exp dt 1 -αijNj(t)

  ). This data comes from only one sample of the process, with discrete timestep dt = 10 -1 . The dashed lines correspond to the read TTI value of h and q d .

	correlation [N(t) N(t )]	0.00 0.02 0.04 0.06 0.08 0.10 0.12	10 1	10 0	10 1 time t t	10 2	10 3 t = 100 t = 0 t = 1 t = 10 t = 500 q 0

Figure 10

. RS correlation E[N (t)N (t )] = C(t, t ), plotted as a function of t -t for different t . This data is from the same sample as figure

9

. We see that, up to fluctuations, the correlation collapse as a function of t -t for t > twait ∼ 20.

Acknowledgments -We acknowledge stimulating discussions with G. Bunin and G. Parisi on this subject. This work was supported by the Simons Foundation Grant on Cracking the Glass Problem (# 454935 Giulio Biroli).

Issues of our implementation

After careful tests on simpler models, we used this scheme to compare with the theory. Initially we used a hardwall immigration at λ = 10 -3 . The agreement was quite good for second degree observables (q d , q0), but not for h. This is due to the numerical scheme. Indeed, if T is quite high (T λ), the sampling of the demographic noise sends many O(1) species close to 0, then they bounce on the wall and end up at N = 2λ. Therefore there is an induced concentration of species at N = 2λ. Because of the 2λ peak, there is a subsampling of the O(1) populations.

In order to reduce this issue, we use a higher λ = 10 -2 in the final results that are shown on figure 11. We reckon the slight discrepancy at high temperature between theory and numerics comes from this issue. A solution is still under investigation in Appendix H . We are aware that the method is still in development. However, it is already enough at the moment to beautifully confirm the theory.

Appendix H: Ongoing investigations 1. χ4

A cleaner numerical test for the transition RS to 1RSB would be the divergence of the χ4 correlation. So far, we do not have enough data to present clean results, but in principle this observation should not be too difficult.

Improve the numerical scheme

In the current numerical scheme, we first sample pure demographic noise then implement the hard wall immigration. When doing this, a lot of trajectories do bounce on the wall, which lowers the accuracy of the scheme. A way to solve this would be to directly solve the Fokker-Planck equation associated to the whole process "demographic noise + hard wall". We reckon this can be done adapting the proof from [78].