Marginally Stable Equilibria in Critical Ecosystems

In this work we study the stability of the equilibria reached by ecosystems formed by a large number of species. The model we focus on are Lotka-Volterra equations with symmetric random interactions. Our theoretical analysis, confirmed by our numerical studies, shows that for strong and heterogeneous interactions the system displays multiple equilibria which are all marginally stable. This property allows us to obtain general identities between diversity and single species responses, which generalize and saturate May's bound. By connecting the model to systems studied in condensed matter physics, we show that the multiple equilibria regime is analogous to a critical spin-glass phase. This relation provides a new perspective as to why many systems in several different fields appear to be poised at the edge of stability and also suggests new experimental ways to probe marginal stability.

Many complex systems in Nature organize in states that are poised just at the edge of stability. The growing evidence comes from physics [1], biology [2], ecology [3], neuroscience [4,5] and economy [6]. One important common trait of all examples is that they are formed by strongly interacting units-species, neurons, agents and particles depending on the situation. The possible explanations of such phenomenon are varied. They include the need for flexibility and adaptiveness to time-varying conditions [2,7], balance between functionality and stability [7], self-organized criticality [8], self-organized instability [9], and continuous constraints satisfaction [1]. Here we address this problem focusing on generalized Lotka-Volterra (LV) equations. They provide a simple and general setting to study assemblies of interacting degrees of freedom; as such they are used in several fields [10][11][12][13]. In particular, they provide a canonical model for ecosystems, with growing connections to systems across biology [10,11,14]. The study of stability of equilibria and their properties using LV equations and generalizations has become a very active research subject. Several important results were obtained recently; in particular general techniques to count the number of equilibria and their properties have been developed [15], and criticality and glassiness have been found to be emergent properties of ecosystems [16][17][18]. Our approach unifies these different perspectives and, by a mapping to condensed matter systems, reveal their generality beyond LV models. Henceforth, in order to describe it, we shall use the terminology employed in theoretical ecology.
In the model we consider, an ecological community is assembled from a pool of available species. We focus on the case relevant for the examples cited above, and in many other situations, when the number of species is large. Since the detailed parameters of all interactions are not known in the majority of cases, and in any case not all details are expected to matter [19], we follow the long tradition pioneered by May in ecology [20] and Wigner in physics [21], and sample the interactions randomly. However, we go beyond May's classical work since randomness is here introduced at the level of interactions between all possible species, while the community self-organizes by choosing which species are present. In other words, the number and identity of the species that are present in the community is selected dynamically [22,23]. Understanding the emergent stability of the equilibria reached dynamically and its dependence on the external parameters is the main purpose of this work. We find, in agreement with [16,18,24], that when the interactions are weak or highly uniform, only one equilibrium is present and is determined mainly by selfregulation within each species. For stronger and more heterogeneous interactions, multiple equilibria emerge. Our main result is that when this happens, all possible states of the system are close to be marginally stable for large number of species and this determines the diversity of the ecosystem, see Fig. 1. Marginal stability has several important consequences, in particular it leads to extreme susceptibility to small perturbations. This situation is referred as "critical" in the physics literature [25]. May famously suggested that complexity and interactions limit the stability of ecosystems [20]. Our results provide a complementary perspective: complex ecological communities reduce dynamically their instability through a reduction of the possible number of surviving species, i.e. diversity, and eventually reach a marginally stable state saturating May's bound. Since this phenomenon stems from a dynamical process, it holds for a broad range of system parameters. It is robust against a range of variations in the model, including different functional forms of responses and interactions, as well as noise. Although in many physical cases criticality arXiv:1710.03606v1 [cond-mat.stat-mech] 10  emerges only at phase transitions, i.e. for very special values of the parameters, there also exist critical phases of matter which instead cover a wide portion of their phase diagram. By relating the LV model to systems studied in condensed matter physics, the multiple-equilibria regime is shown to be akin to a critical spin-glass phase. This suggests that the applicability of our results goes well beyond the LV model we consider and it offers a possible explanation of why so many different systems are found at the edge of stability: they are in a critical marginally stable phase. It also makes clear that this result, while general, is expected to have a well-defined regime of validity, as we shall explain at the end of this work. Finally it makes predictions on some distinctive features of the dynamical behaviour of ecological systems at criticality, which will be interesting to test.
The Lotka-Volterra model we focus on is defined as follows. There are S species in the regional pool, whose abundance is N i ≥ 0. The dynamical equations read (1) where r i is the intrinsic growth rate of species i, and K i is the carrying capacity. It corresponds to the equilibrium abundance to which species i would self-regulate in absence of interaction. For sake of clarity, in the following we focus on the case where r i and K i are constants (set equal to one by rescaling the other parameters). Later, we shall consider the effect of variability in r i and K i , and also different functional responses by replacing The interaction between species is encoded in the matrix α ij . We also add a small (infinitesimal) immigration rate λ to ensure that all invadable species exist 1 . Finally, η i (t) is a white noise with variance 2ω 2 , and √ N i captures the scaling of demographic noise 2 . We consider a symmetric interaction matrix α ij = α ji , corresponding to competitive (or weakly mutualistic) interactions; we will discuss in the conclusion the effect of asymmetry. Except for this constraint, no additional structure (such as trophic levels or space) is included, and the entries α ij are taken to be independent identically distributed random variables. Note that, as already anticipated above, the assumption on the randomness is done at the level of the pool and not of the community. The random variable α ij can be drawn from any distribution without long tails, all that matters are its mean and variance. It turns out that the parameters that play a role in the final theory are the average number of links, C, per site and the first two moments of α ij though the combination µ = Cmean[α ij ] and σ 2 = Cvar[α ij ]. For the sake of clarity, we now focus on the case C = S in which all species interact, extensions are discussed at the end of this paper. The LV eqs. can be rewritten in a way that makes their relationship with stochastic equations studied in physics more transparent: Without noise these eqs. admit a Lyapounov function In presence of noise, eqs. (2) are generalized Langevin equations. In the SI we show that they represent equilibrium dynamics of a thermal system with temperature T = ω 2 and characterized by the following effective Hamiltonian, or energy, H = −L + i T log N i . As a consequence, the long-time stationary probability distribution is the Boltzmann law: P = e −H/T /Z, where the partition function Z guarantees the normalization. This result reveals that understanding the equilibria and the dynamics associated with the LV eqs. (1) can be exactly reformulated as a problem of statistical physics of thermal disordered systems, in which the N i represent the degrees of freedom interacting via random couplings α ij and subjected to individual Therefore, we can deduce properties of the equilibria reached dynamically by a thermodynamical analysis. In particular, without noise, i.e. at zero temperature, the equilibria correspond to the minima of the energy function, i.e. to the ground state and the metastable states.
In order to obtain the full solution of this model we shall use tools developed in statistical physics of disordered systems, including replica computations, discussed in detail in the SI. The phase diagram without noise, i.e. at zero temperature, and for small migration, λ → 0 + , has been obtained in [24]. We reproduce it in the SI for completeness (Fig.  6) and show that it coincides, as it should, with the one obtained by the replica method. One finds that when interaction strengths are all identical, all species coexist in the community when µ > 0 4 . By increasing the variability in the interactions, more and more species are driven out of the community. They are characterized by N i = 0 for λ → 0 + (we will call them "extinct" henceforth). The equilibrium reached dynamically is stable under perturbations, that is by changing and there is a gap in the spectrum of the corresponding stability matrix 5 [24], see Fig. 1(A). Note that in this regime the final community composition is unique, independent of the assembly history, e.g. , the initial conditions for the dynamics. This picture persists up to σ c = 1/ √ 2. By increasing randomness in the interactions above σ c a transition to another regime, which is sharp for large systems, takes place. Our purpose is to study the phase reached when crossing the transition. Henceforth we continue to focus on the zero-noise case, the effect of the demographic noise is discussed at the end of the paper. In physics terms, the single equilibrium regime corresponds to a "paramagnetic phase" where the zerotemperature values of the degrees of freedom N i are mainly fixed by the external potential V i (N i ) and a meanfield anti-ferromagnetic (competitive) interaction. By increasing the randomness in the α ij s the system undergoes a zero temperature phase transition toward a spin-glass phase, characterized by many local minima of the energy (or maxima of the Lyapounov) function and, hence, multiple equilibria. We have used the replica method to study it (see SI) and found that the regime with multiple equilibria corresponds, technically, to a full replica symmetry breaking solution. On the basis of all previous analysis of mean-field 6 spin-glasses [26,27], we can then make general statements about the regime with multiple equilibria. First, it is characterised by a large number of equilibria. These equilibria are minima of the energy, separated by regions with higher energies that form what are called barriers. The lowest equilibria are typically separated by barriers that diverge in the large S limit, while the higher ones by barriers of order of one, i.e. that do not scale with S [28]. Second, and central to our discussion, all minima display a stability matrix characterized by arbitrary small eigenvalues for large S, i.e. minima are marginally stable and characterized by flat directions in the energy landscape at quadratic order. The ground state has this property and, naturally, the higher energy local minima too. We now explain the main findings of our thermodynamic analysis and relate them to random matrix theory results, see SI for details. The two main observables we focus on are: (i) where the star indicates that only non-extinct species are considered and (ii) φ ≡ S * /S, which is the fraction of species present in the community, called diversity in what follows. For identical interaction strengths, i.e. at σ = 0, all species coexist (φ = 1) and is constant across species and increasing. Concomitantly, the diversity decreases. As found in [24], at σ c = 1/ √ 2 the system undergoes a sharp transition from the single to the multiple equilibria regime. This corresponds physically to a phase transition to the spinglass phase. In this phase we find that for all equilibria and all species i, where φ c is the value of the diversity at the transition. These identities hold throughout the entire multiple equilibria phase and they are consequences of the criticality of the spin-glass phase. We also find that for σ > σ c the number of equilibria, i.e. of energy landscape minima is exponential 7 in S. Figs. 2 and 3 confirm our analytical predictions by showing respectively numerical results for and φ corresponding to a given equilibrium reached dynamically. 7 The number of minima scales as e hS where h goes to zero at the transition from one to multiple equilibria. One cannot do general statements on its order of magnitude, since it depends on the external parameters and on the model, in particular the functional response f (N ). When comparing to numerical and experimental results, it is important to keep in mind that it can be small, as we found for instance for the standard LV model. In consequence even for rather large S ∼ 100 the number of equilibria may be modest, see SI.
The second identity in (3) corresponds to a saturated form of May's original bound 8 , var(α ij ) = 1/(4S * ) = 1/(4C * ), where C * is the average number of interactions per surviving species [20]. In order to reveal this connection with random matrix theory we focus on the S * × S * stability matrix M * associated to a given equilibrium, defined by the relation ( Using the fixed point equation corresponding to (1) it is easy to check that In this equation the indices i, j have to be reduced to the surviving species since extinct species remain so if one adds an infinitesimal perturbation ξ j , and do not contribute to the stability of the equilibrium reached dynamically 9 . Following procedures developed for meanfield spin-glasses [29], one can show (see SI) that the spectrum of M * ij is identical for large S to the one of a S * ×S * matrix with independent identically distributed Gaussian off-diagonal entries having the same first and second moment of α ij . This is by no means trivial since the equilibrium reached dynamically, and hence the identity of the surviving species, depend on α ij and induce correlations in the off-diagonal elements of M * ij [30]. The relation with random matrices implies that the eigenvalue density of the stability matrix is a Wigner semi-circle with support , as we indeed find numerically, see Fig. 1. Moreover, this directly connects (3), which holds in the entire spin-glass phase, to marginal stability. We have therefore recovered May's original stability bound but in a saturated form: the number of surviving species, Sφ, is exactly the one guaranteeing that the system is poised at the edge of stability, similarly to what was proposed in the self-organized instability scenario [9].
Let us now discuss extensions and the range of validity of our results. We have verified that our conclusions on the multiple equilibria regime continue to hold for several different convex functional responses f (N ), the standard f (N ) = N (1 − N ) being only an example among others, and with variability in the values of r i and K i . This is a direct consequence of the properties of the spin-glass phase to which the multiple equilibria regime is related to. In these more general cases, the critical character of this phase is encoded in the following identity valid for the average of the square of the single species response 8 The prefactor 4 comes from the symmetry of α ij in the present model. 9 In the limit of small migration, λ → 0 + , extinct species are those that cannot invade: Adding an infinitesimal ξ i does not change this property and thus the species remains extinct. in the whole multiple equilibria regime: Note that eq. (5) reduces to the first identity in eq. (3) when single species response are identical, as previously discussed. We show in Fig. 4 the numerical results obtained for f (N ) = 1 − N − (N − 1) 2 /4 confirming this prediction: the RHS of eq. (5) is less than one in the single equilibrium phase, it reaches one at the transition and then remains stuck to this value in the whole multiple equilibria phase. As before, the criticality of the spin-glass phase implies marginal stability. Indeed, similarly to the standard LV case, the spectrum of the stability matrix M * ij = V (N * i )δ ij + α ij can be shown to be identical for large S to that of a S * × S * matrix with independent identically distributed Gaussian off-diagonal entries having the same first and second moment of α ij , and independent identically distributed diagonal entries with the same statistics of V (N * i ). As we show in the SI, the condition that the left edge of the eigenvalues density touches zero for this class of random matrices is φσ 2 1 ii turns out to be identical to eq. (5). Therefore we obtain that the multiple equilibria regime is indeed generically characterized by marginal stability and, by doing so, we derive a new generalized version of May's bound (eq. (5)). Remarkably, these properties hold despite the fact that for this general class of random matrices the density of eigenvalues is no longer a shifted semi-circle and the singularity at the left edge is not necessarily a square root 10 , as shown numerically in Fig. 1 for the f (N ) = 1 − N − 3/4(N − 1) 2 case (which has f (N ) > 0 for some values of N , corresponding to an Allee effect). Note that for general f (N ) the phase diagram is modified compared to the standard LV model. For instance for f (N ) = 1 − N − 3/4(N − 1) 2 the single equilibrium phase is absent even for infinitesimally small interactions. In conclusion, our investigations show that the class of f (N ) leading to the marginally stable multiple-equilibria phase, i.e. the phase boundaries of the critical spin-glass phase in the physics terminology, is quite large. Determining its boundaries is an important and interesting task that we leave for future studies. Based on previous results on mean-field glassy systems [32], it is possible that the property that the multiple 10 The singularity at the edge depends on the statistics of the diagonal components M * ii , i.e. of V (N * i ). It is a square root if the distribution of the random variable V (N * i ) approaches the left edge of the support slower than linearly, in the other cases it may be a square root or instead inherit the singularity of V (N * i ) at the left edge [31]. equilibria probed by the system are marginally stable is robust, even though the detailed properties of the landscape might be different depending on the shape of f (N ) 11 . Another extension of our work worthy of future analysis concerns the role of the interactions network. As long as the connectivity C per species is large and the underlying structure rather homogeneous, e.g. no fat tails in the distribution of the local connectivity, the mean-field approach we developed is a very good approximation. In consequence, our results are expected to hold also in this more general setting where C 1 but C is not necessarily equal to S. As a first interesting extension one could consider LV models on random regular graphs [33].
The properties (and the existence) of the multiequilibria phase continue to hold also in presence of small noise. On the basis of previous studies on mean-field spin-glasses [27,28], we can state several general facts.
In presence of small noise the system moves around between multiple dynamical states. These correspond to low temperature spin-glass states associated with the local minima discussed above. Only the low energy minima are able to trap the dynamics of the ecosystem over long periods of time, while the ones higher in energy are instead separated by small energy barriers; the transitions between them are so frequent that their identities as separate states disappear even for small noise. The stability of these dynamical states can described by the matrix M ij , which is a generalization of the stability matrix and is defined as the (matrix) inverse of ∂ Ni ∂ξj ( · denotes the average over the noise). As the stability matrix in absence of noise, M ij is positive definite and has arbitrary small eigenvalues for large S, thus leading to marginal stability. Physically, this is a consequence of the fact that the spin-glass phase is not destroyed by small thermal fluctuations and is critical over a finite range of temperatures [26]. In summary, by mapping the LV models to thermal disordered systems and studying their thermodynamics, we find that marginal stability is a property of all communities that are reached dynamically by an ecosystem in the multiple equilibria phase. Eq. (5), combined with our random matrix analysis, relates this property to the single species response. This is the main result of our work: it represents an exact statement of May's stability bound [20], with three notable differences: (1) it follows from an exact analysis of the communities reached dynamically rather than from a-posteriori assumptions on the stability equations, it is thus a property of the emergent community; (2) it is saturated, with an equality rather than an inequality; (3) it is more general, allowing to incorporate non-linear f (N ).
Having established that marginal stability is a generic property of the multiple equilibria phase, we now discuss some of its consequences and propose measurable tests. The most striking effects are expected to appear in dynamical phenomena. Again, previous results on dynamics of mean-field spin glasses provide useful guidelines [27]. In particular, starting the LV dynamics from random initial conditions one expects slow relaxations toward the minima and history dependence for large S. Both phenomena are tightly linked to marginal stability which results in flat directions in configurations space. The response to perturbations is also expected to be very unusual: marginal stability should lead to strong and wildly fluctuating non-linear responses [34] and avalanches of extinctions and invasions [1] 12 . Working out the relevance of this phenomena in various ecological contexts certainly merits future research. The results we found for symmetric interactions have also important consequences for cases where α ij are asymmetric. Indeed, given that the multiple equilibria are marginally stable, we expect that adding asymmetry leads immediately to a chaotic behavior in which the system moves among the different regions of configurations space corresponding to the vestige of those equilibria [36][37][38].
Throughout this paper, we stressed that the unusual properties of the multiple equilibria phase are related to the criticality of the corresponding spin-glass phase. In the following, we show that this relationship also suggests new ways to test for marginal stability. Criticality corresponds to a state in which the microscopic degrees of freedom are all strongly correlated, which naturally leads to singular responses. The properties of the stability matrix in the multiple equilibria phase are one facet of this phenomenon; diverging fluctuations are another. Note, however, that simple fluctuations such as do not capture criticality (the overbar denotes the average across the species). One has to focus on other kinds of correlations. See SI for detailed computations on this point and the following. As a matter of fact, in the LV model we considered, which can be mapped onto thermal dynamics, diverging responses and fluctuations are exactly related by the fluctuation-dissipation relation ∂ Ni . On this basis, and following previous work on glassy systems [39], we propose to probe criticality, or the collective nature of the equilibria and more generally of the dynamical states by looking at a fourth-order correlation function, called χ 4 (t, t ), which reads: where δN i (t) = (N i (t)− N i (t) )/ C(t, t). This function allows one to measure to what extent species are dynamically correlated and is therefore a way to quantitatively test criticality and marginal stability. In the LV model and for σ < σ c the dynamics becomes stationary after a short transient. In this case χ 4 (t, t ) depends on t − t only. Its long time limit (t − t 1) is equal to: whereas for t = t and small temperatures χ 4 (t, t) is equal to twice χ SG . As shown in Fig. 5, where we check these analytical predictions by numerical simulations for S = 400, χ 4 (t − t ) is a decreasing 13 function of t − t . In the condensed matter theory context χ SG is known as the spin-glass susceptibility and is known to diverge in the entire spin-glass phase. We indeed recover this result and connect it to marginal stability since where ρ(λ) is the density of eigenvalues of the matrix M ij . The divergence of χ SG as S 1/3 comes from the fact that the minimum eigenvalue of M scales as S −2/3 [40].
In the inset of Fig. 5 we show the behavior of χ SG obtained analytically in the large S limit. Note that if the singularity of ρ(λ) is milder than a square root, as it is the case for example for f (N ) = 1−N −3/4(N −1) 2 , then one needs to consider high-order moments. The bottom line of this discussion is that if data on the time-dependence of abundances is available, the function χ 4 (t, t ) allows one to measure to what extent species dynamics are correlated and test directly for criticality and marginal stability. Even though in the simple LV case, measuring χ SG would be sufficient for that purpose, measuring the time dependent four-point function χ 4 (t, t ) is the way to go in order to obtain information in more general cases, which may be neither stationary, nor related to thermal equilibrium dynamics.
In conclusion, our analysis of Lotka-Volterra equations in the limit of large species offers an explanation of why many systems in Nature are poised at the edge of stability: we have shown that when the parameters of an ecosystem cross the limit of stability, the system dynamically self-adapts to remain exactly marginally stable. It does so reducing the number of species in such a way to saturate May's bound, which therefore emerges as a result of a dynamical process. This leads to a whole critical phase with multiple marginally stable equilibria, which is expected to be present for several different models and to display highly non trivial dynamical behaviors. Its consequences can be relevant and important in many fields [3,6,7,41,42]. The noisy simulation in Fig. 5 was run with a simple Ito discretization of the stochastic differential equation: where A i are the noiseless terms of dN i /dt and B i = 2T N i . The simulation was run with T = 10 −5 .

Spin glasses
Many of the theoretical methods used in this work were originally developed to study disordered systems in physics, and in particular spin-glasses. Here we make some brief comments on such systems and how they relate to the present problem. The behavior described in this work requires variability in the interaction strengths α ij (as measured by their standard deviation σ). In physics, systems where interactions between the constituents exhibit analogous variations are known as "disordered systems". In particular, a spin glasses is systems where magnetic interactions vary between the different pairs of atoms. Models of such systems, starting with [43], traditionally use binary variables to model the state of each magnetic spin σ ∈ {−1, 1}, while the complex interactions were modeled as i.i.d. Gaussian random variables. For example, the first spin glass model [43] is the Edwards and Anderson model where the Hamiltonian (the energy of a state) is given by where only terms involving pairs of nearest neighbours (i, j) are included in the sum, and J i,j are randomly sampled with zero mean. The first model to be solved is its fully connected analog, the Sherrington Kirkpatrick (SK) model where all pairs of spins are taken to interact [44,45]. A number of other spin glass models were considered along the years [26,[46][47][48]. Among the most interesting ones, we mention those where interactions involve p > 2 variables at a time, called p-spin models [46] It turned out that spin glass models show interesting new phases and phase transitions and can be classified in different universality classes which correspond to different macroscopic behaviours. The techniques developed for the solution of spin glasses are in general useful to describe systems with high level of frustration [49] (i.e. absence of optimal solutions for certain instances of the couplings). For this reason they are widely applied nowadays in different fields including condensed matter (magnetic systems, supercooled liquids), biology (protein folding, neural networks), social sciences, economics, computer sciences (optimization theory, machine learning). The major obstacle that had to be tackled while dealing with the solution of spin glass models is represented by the task of performing disorder averages of quantities of physical interest to extract information on the macroscopic behaviour of the system. The information about equilibrium is contained in the so-called free energy of the system where the overline represents the average over difference instances of the disordered couplings, the inner sum runs over all the possible configurations of spins, and the argument of the logarithm is commonly called partition function The operation of taking the average usually is reduced to perform a Gaussian integration. This would have required little computation effort indeed had the logarithm not been on the way. Yet its presence cannot be neglected nor the operation of taking the average simply performed on the logarithm's argument (the last procedures is called annealed calculation but does not lead to the correct solution in the interesting regimes). To keep the order of the operations and yet end up with an analytically tractable problem the so-called replica trick was introduced [26]. It amounts to use the fact that log x = lim n→0 log x n n , which can be easily verified. Hence the free energy can be written in terms of the so called replicated partition function Z n as where the power n can be interpreted, before the limit n → 0 is taken, as we were focusing on n independent copies of the same system in presence of a unique sample of random couplings. Averages over different realizations of the disorder are in this form straightforward. The copies in the spin glass jargon are usually called replicas.
To give an intuition of the physical meaning of the results that can be obtained within replica computations we must remember that frustrated systems are usually characterized by a multi minima structure of the energy, or any equivalent cost function that might be of interest. This arrangement of minima is uniquely associated to any instance of the random couplings. The role of replicas is the one of revealing the main features of this multi minima structure by independently probing different minima. In fact one of the most important piece of information that comes out from a replica computation is the average width of equilibrium minima and the average distance between pairs of them, or more in general the hierarchical arrangements of minima in the space of configurations. All this is contained in the structure of overlap Q a,b (or similarity, which accounts for the inverse of distance in the space of configurations) between pairs of replicas a, b: where AR denotes the measure over all replicas with the Boltzmann weight and after averaging over disorder (RHS of the first equation above). Along a typical replica computation a new conceptual obstacle arises when the free energy is rewritten in terms of an integral over all the possible choices of the overlap matrix. In the thermodynamic limit (N → ∞) the Laplace method (or saddle point method) can be applied to evaluate the integral but it requires to find the overlap matrix that maximizes the argument of the integral. This operation requires the introduction of a good ansatz for the overlap matrix. The currently used scheme was proposed by Parisi [50,51] and subsequently proved to be the one providing the correct saddle point result [52]. It is called Replica Symmetry Breaking (RSB) scheme and will be discussed in more details in the following sections.
Depending on the number of steps of breaking of the replica symmetry required to get a meaningful solution [53] (i.e. stable in the replica space), we could end up working with Replica Symmetric (RS) scheme, one step Replica Symmetry Breaking (1RSB), ∞ steps Replica Symmetry Breaking (FRSB), just to mention the most relevant ones. This differentiation allows to classify spin glass models and characterize the features of their relevant phases and phase transition. It turns out that the ecological model we consider in this work, at large values of σ, is characterized by a FRSB solution. The FRSB solution represents, as stressed in the main text, a critical phase. From the technical point of view, it is marginally stable, meaning that within the Laplace method it corresponds to an extremum with vanishing small eigenvalues.

Mapping to a thermal disordered system
We assume the starting equation to be with α i,j a symmetric matrix, and η i a white noise: η i (t)η j (t ) = 2T δ ij δ(t − t ). Lotka-Volterra (LV) classical system of equations for interacting species can be obtained by setting f i (N i ) = K i − N i . By scaling the variables N i by K i = µ K : N i = N i /µ k and consequently also the relevant parameters and also the function f i ( N i ) = f i (N i )/µ K (which is true in the LV case) we get an identical dynamical equation in terms of the new variables and the old α ij , and r i . Without loss of generality then we set µ K = 1 and forget about all the tildas. We also prefer to define an interaction matrix θ ij = r i α ij /K i = ρα ij . In doing this we assume that the ratio ρ = r i /K i is i-independent so that α ij and θ ij are symmetric at the same time.
The dynamical equations then read We want to show that this equation admits an invariant probability distribution in terms of a Hamiltonian H. To do so we derive the corresponding Fokker-Planck equation [54]. We consider a generic observable O({N j }) and the time derivative of its average over the thermal noise d dt O ({N j }) . This derivative will obtain us the time derivative of the distribution of our variables originated by the thermal noise itself: Adopting the Ito convention [54], the LHS of the equation (8) corresponds to In the Ito prescription variables are not correlated with noise at the same time so, also using equation (7), the previous equation becomes where Re-writing now the average over the noise as an average over P ({N j }, t) and integrating by parts we get By comparison with equation (8) we can now write the Fokker-Planck equation as The equilibrium distribution must satisfy then which is obtained by imposing ∀i By asking that this equilibrium distribution is also of the form with Z the usual normalizing partition function, we obtain and hence that so finally Note that having λ > T (even when T → 0) is a fundamental element to have a regularized P ({N j }) at small N j . In conclusion the original dynamical equations describe the dynamical evolution of a system whose thermodynamics is determined by the Hamiltonian just obtained.

Replica computation
We use a replica approach to analyze the thermodynamics of the system characterized by the Hamiltonian (18). Recall that θ ij are assumed to be Gaussian distributed with mean ρµ/S and variance ρ 2 σ 2 /S. We evaluate the free-energy of the system by applying the replica trick to perform sample averages We then evaluate the replicated partition function as where the average over the disorder contained in V remains to be done, and which, through standard replica manipulations [26], becomes In the last equation the effective on-site partition function Z i = a dN a i exp (−βH ef f ({N a } i )) is obtained from an effective Hamiltonian where the order parameters satisfy the self-consistent equations

Replica Symmetric Solution
Replica symmetry is correct when only a single equilibrium (or minimum) in the (free-)energy landscape governs the thermodynamic behavior. Formally, one assumes: The free-energy expression becomes in this case with an effective Hamiltonian for Z i To decouple replicas we exploit standard properties of Gaussian integrals and we get By maximizing the action S at the exponent of (24) we get the following saddle point (SP) equations on the introduced parameters where b, c denotes replica indices and p, r, are the powers of the abundance we are considering. In the n → 0 limit the formula above can be expressed in terms of thermal averages · 1R over single species and single replica with Hamiltonian H RS , and the disorder average · representing the average over the disorder contained in V (N ) and the Gaussian integral over z with mean zero and unit variance where Hence we can write The average over a single replica coincides with the standard average first over the Boltzmann weight and then over the quenched disorder. It is possible to reduce the latter to the former because the model is mean-field. The resulting equations have a clear interpretation: each species is subjected to its own potential V and two extra terms due to the overall mean-field interaction with the rest of the system. Two non fluctuating terms, one quadratic and the other linear, plus a fluctuating linear term proportional to z. The latter can make the minimum of the overall potential zero (extinction) or larger than zero (survival).
Note that for fluctuating K i , the LV potential is and the average in the 1replica computation above is performed over the Gaussian distributed K i with mean µ K = 1 and variance σ K .

Zero Temperature Limit
In the zero temperature limit q 0 → q D at the same pace as T so it is useful to define the variable ∆q = ρβ(q D − q 0 ), where ρ has been inserted just for convenience in writing. The equations of the three parameters in this limit are hence expressed in terms of h, q 0 , and ∆q. In this limit the thermal averages over the 1replica measure are evaluated by saddle point-method at N * , which is the positive minimum of the Hamiltonian H RS (N ) when it exists, or zero. Hence the SP equations become where θ(x) is the Heaviside function, θ(x) = 1 for x > 0 and zero otherwise. Note that in the case of last equation on ∆q = N 2 − N 2 we had to Taylor expand H RS for small T separately in the case of extinction N * = 0 and survival N * > 0. The LV case is particularly simple since N * reads Until now K and z were two separate Gaussian variables (with averages 1 and 0 and variances σ K and 1, respectively) over which we are averaging. Combining together these two variables into z with 0 average and variance 1 we get with value of the random variable z corresponding to extinction −∆ = − 1−µh √ q0σ 2 +σ 2 K . The expression for q 0 , h, and ∆q are hence immediately obtained as being and with These equations coincide with the one obtained by the cavity method in [24].

One step replica symmetry breaking equation
In the multiple minima phase the RS solution is unstable, as already checked in [24]. This implies the existence of multiple equilibria. In order to characterize these equilibria one has to study what kind of RSB solution emerges. In the following we consider the 1RSB solution: the n replica are divided into n/m groups and Q ab = q 1 for a = b both in the same group, Q ab = q 0 for a, b in different groups and Q aa = q D and H a = h. Once introduced this ansatz the computation is similar to the RS one. The free-energy expression is in this case with an effective Hamiltonian, H ef f ({N a } i ), for Z i : To decouple replicas we exploit standard properties of Gaussian integrals and we get By maximizing the action A at the exponent of (43) in the n → 0 limit we get the following SP equations on the introduced parameters and H a(a B ) (N a(a B ) , z, z a B B ). These averages in the n → 0 limit can be expressed in terms of thermal averages · 1R over single species and single replica with Hamiltonian H 1RSB (N, z, averages · mR over the Gaussian variable z B with additional weight given by and averages · V representing the average over the disorder contained in V (N ) and the Gaussian integral over z with mean zero and unit variance. Using all this we have Hence we can write

1RSB Zero Temperature Limit
Also in the case we have to considered rescaled variable in the limit T → 0: ρ(q D − q 1 )β = ∆q ∼ O(1), and the scaling of the replica breaking order parameter m is such that βm remains of the order of one. Hence, in the following we will introduce the notation βm = m and keep m ∼ O(1). In this limit, similarly to the RS case, the SP equations read as follows: Previously though the expressions were even simpler because we had to compute averages of the kind  (N  *  (z, z B ), z, z B )] when N (z, z B ) * = 0. For the same reason the normalization constant is non trivial and must be evaluated. The LV choice of V allows for simple explicit expression of the equations. In particular as in the LV RS case, we combine the Gaussian variables z and K into z and for every z B we get with the new value of the random variable z B corresponding to extinction For every given z B the additional weight involving when N * is non null. Hence the normalization constant is With this in mind and defining we can finally write the 1RSB self consistence equations as follows .
Everywhere we could determine also the m given by a SP equation, which satisfies the following condition .
What we do is instead to use m as a parameter through which we can select minima of the 1RSB structure at different energy levels. This allows to compute the number of minima with a given energy, using m as a parameter conjugated to the energy [48]. The logarithm of the number of minima divided by S is called configurational entropy. It is proportional to the derivative of the free energy with respect to m [48]. Note that, by definition of m, the configurational entropy of minima corresponding to the equilibrium in the 1RSB phase is null.
In the n → 0 and T → 0 limit the free-energy reads and the configurational entropy is .

Instability of the 1RSB phase and marginality condition for the FRSB phase
We now study the (in)stability of the 1RSB phase and, more generally, obtain the condition for the stability of RSB phases. To obtain the stability condition we consider a generic k − RSB phase and study fluctuations δQ ab only inside the inner blocks of the Parisi matrix. This is the so-called replicon eigenvalue and corresponds physically to fluctuations within a state. As we shall discuss in the next section, this is directly related to the Hessian (stability matrix) around one equilibrium. We call L the action that has to be extremized at the saddle-point and we study its Hessian with respect to the fluctuations described above: where the average is done over the effective Hamiltonian. All replica indices belong to the same block of size m × m, one of the inner ones. A single inner block is analogous to a replica m × m symmetric matrix so the corresponding Hessian matrix can be diagonalized rather easily (see Almeida and Thouless). Within the same block there are three independent matrix elements depending whether some replica index are the same: The replicon eigenvalue (see [53]) is λ = P − 2Q + R with degeneracy m(m − 1)/3. The condition for stability is λ ≥ 0. Marginal stability corresponds to λ = 0.
To evaluate the replicon eigenvalue in the 1RSB phase we need to consider a single block a B and evaluate Hence the replicon eigenvalue can be expressed in a more transparent way as where the second moment of N 2 within one single state (or equilibrium) appears. Using the fluctuation dissipation relation one can rewrite the previous equation in term of single species responses: In the FRSB phase the replicon is exactly zero, this is related to the criticality of the phase [26], thus one obtain the equation: (σρ) 2 ∂N ∂ξ 2 = 1 which encodes the marginality condition at finite temperature. In the small T limit this computation is analogous to the one performed for ∆q. At the end one gets the simpler expression which leads to eq. (5) of the main text. For the usual Lotka-Volterra case in which V (N ) is quadratic one gets H 1RSB (N * ) = ρ(1 − σ 2 ∆q) which does not depend on N * . Thus the replicon eigenvalue and consequently the marginality condition is particularly simple: As spotted earlier, the expression of ∆q in terms of correlation function is very similar. In the LV case the similarity becomes even closer: .
Using together the equation on ∆q and the marginality condition one finds two simpler appealing expressions: where φ = θ(N * ) . The first equation is the the limit of stability given by the May Bound: the fraction of surviving species in any equilibrium should be such that the Wigner semi-circle touches zero. The second is a general result valid in the marginal phase. These analytical predictions have been tested in Fig.2 and Fig.3 of the main text.

Phase diagram and numerical solution of the mean-field equations
The replica symmetric phase was already studied in [24]. Our results agree with the previous one. In particular we find three phases, see Fig. 6. By increasing σ for µ > 0 the single equilibrium phase becomes unstable toward the multiple equilibria (spin-glass) phase when its replicon eigenvalue vanishes. The instability toward the unbounded growth phase is signalled by a concomitant divergence of N and N 2 − N 2 . Note that the transition line to the unbounded growth phase was determined within the RS ansatz so it is only an approximation for µ > 0. We have checked by numerical simulations that it is actually a good approximation. Crossing the transition toward the multiple equilibria phase one finds that the RS phase becomes unstable and one has to break replica symmetry. We have found that also the 1RSB solution is unstable even though much less than the RS one, see Fig. 7 where the replicon eigenvalue is plotted for the RS and the 1RSB phases for the standard LV model with r i = K i = 1 for µ = 2 as a function of σ. We didn't look for 2RSB solutions and directly assumed that the stable phase if the FRSB one as found generically in spin-glass models [26]. We validated this assumption by comparison with numerical simulations that show marginal stability in the multiple equilibria phase, a property valid only for the FRSB phase. Note that, although unstable, the 1RSB provides a very good approximation as we have checked by comparison with numerical simulations. For example, in Fig. 8 and 9 we show 1−σ 2 ∆q and φσ 2 for the standard LV model with r i = K i = 1 for µ = 2. These two quantities have respectively to stick to the values 2 and 1/4 as discussed in the main text and found by numerical simulations. As shown the 1RSB is already a very good approximation of the correct results, corresponding to the FRSB phase. We have also computed the configurational entropy. Given that the 1RSB is unstable, we cannot determine even approximatively the most numerous equilibria [26]. The values we found for the configurational entropy as a function of energy within the 1RSB ansatz for the standard LV model with r i = K i = 1 for µ = 2 and σ = 0.88 are very small, in the range 10 −3 − 10 −4 . It would be interesting (but also quite involved) to obtain the correct result within a FRSB computation. Anyhow, it is important to keep in mind that the number of equilibria in realistic situations can be modest depending on the value of the configurational entropy and the total number of species.  9. φσ 2 as a function of σ for the standard LV model with ri = Ki = 1 for µ = 2 from the RS and 1RSB solution. For σ > σc both phases are unstable but the 1RSB result is very close to the correct one corresponding to φσ 2 = 1/4.

Random Matrix Analysis
As explained in the text, the stability of a given equilibrium is governed by the S * × S * stability matrix M * ij which is defined by the equation (M * ) −1 ij = ∂N * i ∂ξ * j and reads In order to study its spectral properties we focus on the resolvent, defined as G(λ) = (S * ) −1 Tr(λ1 − M * ) −1 and add an infinitesimal negative imaginary part to λ. Following exactly the same procedure developed for mean-field spin-glasses in [55], one can construct a perturbative expansion in α ij , sum the leading contributions for S → ∞ and obtain the equation: valid only for indices corresponding to surviving species. By summing over the surviving species one finds where the average is over the distribution of the V (N * i )s. This equation allows one to study the density of eigenvalues ρ(λ) of M * thanks to the relation ImG(λ) = πρ(λ). Since eqs. (64,65) are also the equations satisfied by the resolvent of a random matrixM * with independent identically distributed Gaussian off-diagonal entries having the same first and second moment of α ij , and independent identically distributed diagonal entries with the same statistics of V (N * i ), we conclude that M * andM * are equivalent as far as the average spectrum is concerned (a relation that we checked explicitly by numerics). A marginally stable equilibrium is characterized by arbitrary small eigenvalues of its stability matrix, i.e. it is such the left edge of the support of ρ(λ) is zero. This implies that ImG(λ) becomes arbitrary small for λ → 0. In consequence, close to λ = 0, we can develop the self-consistent equation on the resolvent as: The condition for having a non-zero imaginary part is that when collecting all terms on the RHS the coefficient on the linear term in ImG is positive for λ > 0 and vanishes at λ = 0 14 . This leads to the equation Using relation (64), and replacing 1 V (N * i )+σ 2 φRG(λ) by (M * ) −1 ii in the identity above, we obtain the equation for marginal stability quoted in the text: Dynamical four-point correlation function χ4(t, t ) In the following we derive the analytical results quoted in the main text on χ 4 (t, t ) in the limit of small noise. First, we define the S * × S * matrix A as Let's call also N * i the abundance of the surviving species in the limit of zero noise. For small noise their abundances have fluctuations of the order √ T around the zero-noise value, whereas instead abundances of species with N * i = 0 have fluctuations of the order T . For small noise and at leading order, one can do a quadratic approximation of the Hamiltonian and find for the surviving species: The definitions of C(t, t ) and χ 4 read The correlation at equal time is related (for small noise) to α * via Note that the species characterized by zero abundance in the zero noise limit do not contribute at leading order in T since they would give a contribution O(T 2 ). At long times the correlation function vanishes This also happens for the correlation between different species: And for the correlations In the standard LV parameterization, the center is at b = 1 and a = 2σ √ φ, so one finds and χ 4 (t, t) = 2 As discussed in the main text and shown in Fig.5, C(t, t) is featureless at the transition while χ 4 (t, t) diverges.