Properties of equilibria and glassy phases of the random Lotka-Volterra model with demographic noise
PProperties of equilibria and glassy phases of the random Lotka-Volterra modelwith demographic noise
Ada Altieri, Félix Roy,
2, 1
Chiara Cammarota,
3, 4 and Giulio Biroli Laboratoire de Physique de l’École normale supérieure, ENS, Université PSL,CNRS, Sorbonne Université, Université de Paris F-75005 Paris, France Institut de physique théorique, Université Paris Saclay, CEA, CNRS, F-91191 Gif-sur-Yvette, France Department of Mathematics, King’s College London, Strand London WC2R 2LS, United Kingdom Dip. Fisica, Università "Sapienza", Piazzale A. Moro 2, I-00185, Rome, Italy
In this letter we study a reference model in theoretical ecology, the disordered Lotka-Volterramodel for ecological communities, in the presence of finite demographic noise. Our theoreticalanalysis, which takes advantage of a mapping to an equilibrium disordered system, proves thatfor sufficiently heterogeneous interactions and low demographic noise the system displays a multipleequilibria phase, which we fully characterize. In particular, we show that in this phase the number ofstable equilibria is exponential in the number of species. Upon further decreasing the demographicnoise, we unveil a
Gardner transition to a marginally stable phase, similar to that observed injamming of amorphous materials. We confirm and complement our analytical results by numericalsimulations. Furthermore, we extend their relevance by showing that they hold for others interactingrandom dynamical systems, such as the Random Replicant Model. Finally, we discuss their extensionto the case of asymmetric couplings.
Introduction – Lotka-Volterra equations describing thedynamics of interacting species are a key ingredient fortheoretical studies in ecology, genetics, evolution andeconomy [1–6]. Cases in which the number of species isvery large are becoming of general interest in disparatefields, such as in ecology and biology, e.g. for bacte-ria communities [7, 8], and economy where many agentstrade and interact simultaneously both in financial mar-kets and in complex economic systems [9, 10].The theoretical framework used in the past for a smallnumber of species is mainly based on the theory of dy-namical systems [11–16]. When the number of ordi-nary differential equations associated with the Lolta-Volterra (LV) model becomes very large, i.e. for manyspecies, methods based on statistical physics become ide-ally suited. Indeed, several authors have recently inves-tigated different aspects of community ecology, such asproperties of equilibria, endogeneous dynamical fluctu-ations, biodiversity, using ideas and concepts rooted instatistical physics of disordered systems [5, 17–27]. Simi-lar investigations have been also performed for economicsystems [28]. The complexity of dealing with a large num-ber of interacting species can actually become a welcomenew ingredient both conceptually and methodologically.In fact, different collective behaviours can emerge. Asit happens in physics, such phases are not tied to thespecific model they come from, hence they can be char-acterized and characterize systems in a generic way [29].From this perspective, it is natural to ask which kind ofdifferent collective behaviors arise from LV models in thelimit of many interacting species and what are their mainproperties [19, 20]. These questions, which have startedto attract a lot of attention recently, tie in with the anal-ysis of the properties of equilibria [30–32].Here we focus on the disordered Lotka-Volterra modelof many interacting species, which is a representativemodel of well-mixed community ecology [33], and can be mapped or related to models used in evolutionary gametheory and for economic systems [34–38]. We considerthe case of symmetric interactions and small immigra-tion and work out the phase diagram as a function ofthe degree of heterogeneity in the interactions and of thestrength of the demographic noise. Compared to previ-ous works [5, 19, 20, 39] adding demographic noise notonly allows us to obtain a more general picture, but alsoto fully characterize the phases and connect their proper-ties to the ones of equilibria. In particular, we shall showthat the number of stable equilibria in the LV model isexponential in the system size and their organization inconfiguration space follows general principles found formodels of mean-field spin-glasses. Our findings, whichare obtained for symmetric interactions, provide a use-ful starting point to analyze the non-symmetric case, aswe shall demonstrate by drawing general conclusions onproperties of equilibria in the case of small asymmetry.Henceforth we focus on the disordered Lotka-Volterramodel for ecological communities [5, 19] defined by theequations: dN i dt = N i − N i − (cid:88) j, ( j (cid:54) = 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 zeromean and covariance (cid:104) η i ( t ) η j ( t (cid:48) ) (cid:105) = 2 T N i ( t ) δ ij δ ( t − t (cid:48) ) (we follow Ito’s convention). This noise term allows us toinclude the effect of demographic noise in a continuoussetting [40–42]; the larger is the global population thesmaller is the strength, T , of the demographic noise. Im-migration from the mainland is modeled by a reflectingwall for the dynamics at N i = λ , since this is more prac-tical for simulations than the usual way of adding a λ inthe RHS of Eq. (1) (see the Appendix for more details). a r X i v : . [ c ond - m a t . d i s - nn ] S e p The elements of the interaction matrix α ij are indepen-dent and identically distributed variables such that:mean [ α ij ] = µ/S var [ α ij ] = σ /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 dis-tribution: P ( { N i } ) = exp (cid:18) − H ( { N i } T (cid:19) (3)where H = − (cid:88) i (cid:18) N i − N i (cid:19) + (cid:88) i Figure 1. Phase diagram showing the strength of the demo-graphic noise, T , as a function of the degree of heterogeneity, σ , at fixed µ = 10 and selected value of the cutoff, N c = 10 − .Upon decreasing the noise three different phases can be de-tected: i) a single equilibrium phase; ii) a multiple equilib-ria regime between the light blue and the orange lines; iii)a Gardner phase, which turns out to be characterized by ahierarchical organization of the equilibria in the free-energylandscape. transition line corresponding to the blue curve in Fig.1: λ R = ( βσ ) (cid:104) − ( βσ ) ( (cid:104) N i (cid:105) − (cid:104) N i (cid:105) ) (cid:105) = 0 , (5)where β = 1 /T . The average (cid:104)·(cid:105) is the thermodynamicsaverage taken over the effective Hamiltonian (4), while · denotes the average over the quenched disorder associ-ated to the random interactions ( i is a dummy index sincestatistically all species are equivalent after average overthe interactions). Physically, the condition above canbe shown to correspond to a diverging response function[19, 20], and is a signature of the system being at the edgeof stability, namely at a critical point in the parameterspace.Below the blue curve there exist multiple states—which one is reached dynamically depends on the initialcondition. Such states correspond to dynamically fluctu-ating equilibria that are stable to perturbations and thathave typically an overlap in configuration space given by q = 1 S (cid:88) i (cid:104) N i (cid:105) α (cid:104) N i (cid:105) β , (6)where α and β denote the average within two genericstates α and β . One can similarly define the intra-stateoverlap q = S (cid:80) i (cid:104) N i (cid:105) α . See Fig. 2 for a pictorial repre-sentation of these two quantities and the organization ofequilibria in phase space. This is (in the replica jargon)the so-called one-step replica symmetry breaking phase( RSB). In order to characterize the properties of thisphase of the LV-model, we have computed the number q q Figure 2. Zoom on a pictorial landscape. The parameters q and q denote the size of the largest and the innermost basinsrespectively within the two-level structure of the RSB phase. of states, and hence of equilibria, using methods devel-oped for structural glasses [44]. More specifically, wehave computed the complexity Σ (see the Appendix fordetails), which is defined as the logarithm of the numberof equilibria with a given free-energy density f normal-ized by the number of species S . This allow us to showthat 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, theheterogeneity in the interactions becomes even more im-portant and a second phase transition takes place. Inorder to locate it, we repeat exactly the same procedureas for the single equilibrium phase but now within oneof the typical states with a given free-energy f [45]. Thecomputation is more involved (it corresponds to analyzethe stability of the RSB Ansatz) and leads to the con-dition: λ = ( βσ ) (cid:104) − ( βσ ) (cid:104) ( (cid:104) N (cid:105) r − (cid:104) N (cid:105) r ) (cid:105) m -r (cid:105) = 0 (7)where the two different averages correspond to i) theintra-state average (cid:104)·(cid:105) r and the inter-state average, (cid:104)·(cid:105) m -r . All technical details of the calculation will followin the Appendix. The critical temperature that resultsfrom the equation above leads to the orange line in Fig.1. Crossing this line results in a fragmentation of eachstate into a fractal structure of sub-basins [46] (see thelandscape on the bottom in Fig. 1): each state becomesa meta-basin that contains many equilibria, all of themmarginally stable, i.e. poised at the edge of stability [20],and organized in configuration space in a hierarchicalway, as it was discovered for mean-field spin glasses [43].This phase, which is called Gardner , plays an importantrole in the physics of jamming and amorphous materi-als [47]. Our results unveil its relevance in theoreticalecology by showing that it describes the organization ofequilibria in the symmetric disordered LV-model at lowenough demographic noise and for highly heterogeneouscouplings.We now present numerical simulation results that con-firm and complement our analytical study. We numer-ically integrate the stochastic equation Eq. (1) using aspecifically designed method (see the Appendix for de-tails). The initial abundance for each species is drawnindependently in [0 , . time t t c o rr e l a t i o n [ N ( t ) N ( t )] t = 0 t = 1 t = 10 t = 100 t = 500 q Figure 3. Correlation function C ( t, t (cid:48) ) as a function of t − t (cid:48) inthe one-equilibrium (Replica Symmetric) phase. The curvesclearly collapse to the theoretically predicted value, q , fortimes t (cid:48) > t wait . There are three sources of randomness for a given sam-ple: the interactions, the initial conditions and the demo-graphic noise. In the following, we obtain numerically theaverage correlation function defined by E [ N ( t ) N ( t (cid:48) )] = 1 SN sample S (cid:88) i =1 N sample (cid:88) r =1 N ri ( t ) N ri ( t (cid:48) ) (8)where E [ X ] stands for the average over all those sourcesof randomness. If the system size is sufficiently large( S (cid:29) ) as well as the sampling set, with N sample (cid:29) , itcan be shown that the stochastic process converges in law[48]. We generically choose S ∼ and N sample ∼ and eventually verify there is no ( S, N sample ) -dependencyat this scale. We find that in the high-temperature phasea time-translationally invariant (TTI) state is reachedafter a finite time-scale t wait : ∀ t ≥ t (cid:48) > t wait E [ N ( t ) N ( t (cid:48) )] = C ( t, t (cid:48) ) (cid:39) C ( t − t (cid:48) ) . (9)This convergence to a TTI regime is shown in Fig. 3.The long-time limit of C ( t − t (cid:48) ) is the overlap between twogeneric configurations belonging to the single equilibriumstate: the dashed line in Fig. 3 is the analytical predic-tion for lim t − t (cid:48) →∞ C ( t − t (cid:48) ) which is in perfect agreementwith the numerics. We have also checked that this agree-ment holds upon varying T and for other observables, theresults are reported in the Appendix (see Fig. 11). Fromthe time-dependence of C ( t − t (cid:48) ) one can estimate thetypical time-scale characterizing dynamical fluctuationswithin the single equilibrium phase. Formally, we define τ decorrel by the identity: C ( τ decorrel ) − C ( ∞ )( C (0) − C ( ∞ )) = 0 . (10)In Fig. 4, we plot τ decorrel as a function of ( T − T RSB ) ,where T is the critical value of T at which the singleequilibrium phase becomes unstable (blue line in Fig. 1).We find that the thermodynamic instability is accompa-nied by a dynamical transition at which τ decorrel divergesas a power law with an exponent close to . , see Fig. 4.For small demographic noise, i.e when T is below theblue line of Fig. 1, previous results on the dynamics ofmean-field spin glasses [49–51] suggest that the LV-modelshould never reach an equilibrium stationary state, andinstead it should display aging [52, 53]. In fact, one ex-pects that among the very many equilibria the dynam-ics starting from high-temperature-like initial conditionsfalls in the basin of attraction of the most numerous andmarginally stable equilibria, and display aging behaviour.This is indeed what we report in Fig. 5 which shows thatthe longer is the age of the system, t (cid:48) , the longer it takesto decorrelate. The landscape interpretation of this phe-nomenon is that the system approaches at long times apart of configuration space with many marginally stableequilibria. This leads to aging because the longer is thetime, the smaller is the fraction of unstable directionsto move, hence the slowing down of the dynamics, butthe exploration never stops and eventually the systemnever settles down in any equilibrium [54–56]. The twodashed lines in Fig. 5 correspond to our analytical pre-diction for the intra-state and the inter-state overlaps ofthe marginally stable equilibria. The agreement is satis-factory but larger times would be needed to fully confirmit. Temperature T T RSB d e c o rr e l ( T ) -0.56 slope Figure 4. Decorrelation time as a function of ( T − T ) inlogarithmic scale. The blue points correspond to numericaldata, while the dashed red line is a fit. The decay of thedecorrelation time in ( T − T ) occurs with an exponent ≈ − . . Our characterization of the phases and the dynamicsof the LV-model has important consequences on relatedsystems, in particular on the so-called random replicantmodels (RRMs) that consist of an ensemble of replicantsevolving according to random interactions. Given their time t-t' r e s c a l e d c o rr e l a t i o n C ( t , t ) / C ( t , t ) t'=1e+02t'=1e+03t'=1e+04t'=1e+05t'=1e+06 q / q d q / q d Figure 5. Rescaled correlation as a function of ( t − t (cid:48) ) , fordifferent t (cid:48) , showing aging dynamics. The dashed black andred lines correspond respectively to the theoretical predictionsfor q and q both rescaled by the analytical prediction for C ( t (cid:48) t (cid:48) ) (called q d in the Appendix). numerous applications in biology, optimization problems[34, 57] as well as evolutionary game theory [58, 59],RRMs still attract great theoretical interest. The RRM,which was introduced in [34] and further studied in [39], isremarkably similar to the disordered LV-model we stud-ied. In the case of symmetric interactions, one can sim-ilarly map the problem onto an equilibrium statisticalphysics one with the following Hamiltonian: H R = − S (cid:88) i Cracking the Glass Problem ( Appendix A: Thermodynamics and replicaformalism The evolution of the species abundances N i in theecosystem (with i = 1 , ..., S ) is regulated by the followingdynamical equation: dN i dt = − N i ∇ N i V i ( N i ) + (cid:88) j, ( j (cid:54) = i ) α ij N j + (cid:112) N i η i ( t ) + λ (A1)where η i ( t ) is a white noise with covariance (cid:104) η i ( t ) η j ( t (cid:48) ) (cid:105) = 2 T δ ij δ ( t − t (cid:48) ) and T is the tempera-ture. In the presence of a (demographic) noise – anintrinsic population randomness due to birth, death andunpredictable interaction events – Eq. (A1) representsa generalized Langevin equation with a one-speciesquadratic potential V i ( N i ) : V i ( N i ) = − ρ i (cid:18) K i N i − N i (cid:19) , (A2)which allows us to precisely recover the well-knownLotka-Volterra equations [64, 65]. The adimensional pa-rameter ρ i = r i /K i denotes the ratio between the growthrate, r i , and the carrying capacity, K i . For simplicity, wewill assume no species dependence on these parametersand set r i = 1 , K i = 1 in the following.For a first-order analysis, the interaction matrix α ij isassumed to be symmetric: its elements are i.i.d. Gaussianvariables with mean and variance respectively:mean [ α ij ] = µ/S var [ α ij ] = σ /S . (A3)We can also suppose to add a small degree of asymmetryand perform a perturbative expansion. We should ex-pect, even for arbitrary small asymmetry, that the criti-cal (multiple equilibria) phase described in the main textwill be suppressed, as proven long ago in spin-glass andneural network contexts [61, 62].In the case of symmetric interactions, the dynamicalequation (A1) admits an invariant probability distribu-tion in terms of an Hamiltonian operator H , as can beproven by writing the Fokker-Planck equation of the cor-responding stochastic process. More precisely, we cansafely define an Hamiltonian H and solve the problemexactly within the replica formalism [43], as we will showin detail in the next Section. We thus define H = (cid:88) i V i ( N i ) + (cid:88) i 1. Replica symmetric ansatz To deal with the thermodynamics of disordered sys-tems and to give a precise characterization of possiblephase transitions associated with the emergence of com-plex collective behaviours we can resort to the replicamethod [43, 66]. It was originally introduced to studyspin-glass models but has now become a cornerstone for avast class of complex systems. The introduction of repli-cas allows us to handle quantities in which the disorderplays a key role and eventually obtain the free energy bymeans of the identity − βF = lim n → ln Z n n . (A5)In principle, the free energy should depend on the spe-cific realization of the disorder. However, in the thermo-dynamic limit we can safely ignore this issue since thefree energy will converge to a unique value, thanks to itsself-averaging property. Operatively, one should computethe quantity on the r.h.s for integer values of the replicanumber n , then consider the analytical continuation toreal values and only at the end of the computation takethe limit n → .The starting point is then the computation of the repli-cated partition function, which in this specific case be-comes: Z n = (cid:90) (cid:89) i, ( ij ) dN ai dα ij exp − (cid:88) ( ij ) ( α ij − µ/S ) σ /S − βH ( { N ai } ) (A6)where the overline denotes the average over the disorder, i.e. the average over the Gaussian variables α ij . To perform thecomputation – that will involve quadratic terms in the prod-uct of the species abundances – we introduce the overlap ma-trix Q ab (with diagonal value Q aa ) and the external field H a ,where ( a, b ) represent two replicas of the same system, i.e. : Q ab = 1 S S (cid:88) i =1 N ai N bi , (A7) H a = 1 S S (cid:88) i =1 N ai . (A8)We can thus rewrite the free energy in the replica space as: F = − βn ln (cid:90) (cid:89) a, ( a
2. One-step replica symmetry breaking ansatz In the presence of low demographic noise the RS solution isno longer appropriate to describe the thermodynamics of thesystem, which can be more correctly identified by a one-stepreplica symmetry breaking ( RSB) Ansatz. This instabilityis intrinsically related to the emergence of multiple minimain the free-energy landscape, each of them associated with adifferent equilibrium configuration. According to the RSBapproximation, the n replicas are now divided into n/m dif-ferent groups of m replicas, with ≤ m ≤ and n/m integer.The overlap matrix assumes now three different values Q ab = q d if a = bq if a, b ∈ B l q if a, b (cid:54)∈ 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 − off-diagonal elements equal to q and n − m elementsequal to q . A new parameter, m , comes into play whichis usually denoted as breaking point parameter in the replicajargon. The external field is instead assumed to be uniformfor all replicas a . H a = h , ∀ a . (A19)Exactly as in Eq. (A9), the free energy within the 1RSBAnsatz reads: F RSB = − βn ln (cid:90) dq d dq dq dh e S A ( q d ,q ,q ,h ) (A20)We are now able to write the resulting expression for theaction A : A ( q d , q , q , h ) = − ρ σ β n (cid:2) ( n − m ) q + ( m − q + q d (cid:3) ++ ρµβ n h + 1 S (cid:88) i ln Z i , (A21) where the replicated partition function Z i turns out to be Z i = (cid:90) dz i √ π e − z i / n/m (cid:89) a B =1 (cid:40)(cid:90) dt a B ,i √ π (cid:89) a ∈ a B dN ai × exp − t a B ,i − β (cid:88) a ∈ a B H RSB ( N ai , z i , t a B ,i ) (cid:41) (A22)hence, according to Eq. (A16), we have: H RSB ( N i , z i , t a B ,i ) = − ρ σ β ( q d − q ) N i ρµh − t a B ,i ρσ √ q − q − z i ρσ √ q ) N i ++ V i ( N i ) + T ln N i . (A23)To decouple replica indices according to a RSB solution, wehave thus introduced a double Gaussian integration leadingto two different kinds of averages: i) over the single replica;ii) over a group of replicas in the same block of size m . Thetwo averages correspond respectively to (cid:104)·(cid:105) r = (cid:82) ∞ N c dN exp [ − βH RSB ( N, z, t a B )] · (cid:82) ∞ N c dN exp [ − βH RSB ( N, z, t a B )] , (A24) (cid:104)·(cid:105) m -r = (cid:82) dt aB √ π e − t aB (cid:16)(cid:82) ∞ N c dN exp [ − βH RSB ( N, z, t a B )] (cid:17) m · (cid:82) dt aB √ π e − t aB (cid:16)(cid:82) ∞ N c dN exp [ − βH RSB ( N, z, t a B )] (cid:17) m (A25)The values of the overlap parameters q d , q , q and the ex-ternal field h are obtained by saddle-point approximation ofthe action, leading to the following compact expressions q d = (cid:104)(cid:104) N (cid:105) r (cid:105) m -r q = (cid:104)(cid:104) N (cid:105) r (cid:105) m -r q = (cid:104)(cid:104) N (cid:105) r (cid:105) m -r h = (cid:104)(cid:104) N (cid:105) r (cid:105) m -r (A26)that can be solved iteratively up to convergence. The initialcondition is chosen close to the RS solution but with q > q .As the overlap parameters stand for the degree of similaritybetween two replicas in the same block or in different blocks,the found solution is thermodynamically relevant only if q d >q > q . Indeed, in the very low-temperature limit, q d → q ,similarly to the analysis we have performed for the RS solution(see A 1).The determination of the breaking parameter m requiresmore attention. We can follow different routes depending onwhether one aims to find the m -value that optimizes the ex-pression of free energy or the one that identifies the agingsolution in a dynamical framework. In this second case, theparameter m is chosen in such a way that the found solu-tion is marginal: one optimizes over the parameters in Eq.(A26) while selecting m accordingly to the marginal stabilitycondition. We will devote the next Section to a detailed ex-planation of what marginal stability means and what kind ofconsequences it yields for the system. 3. Meaning of the replicon eigenvalue of thestability matrix To investigate the stability of the different phases, we in-troduce 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 expressedin terms of three different sectors. Following [67], we definethe three eigenvalues: the longitudinal , λ L , the anomalous , λ A , and the replicon , λ R . We are specifically interested in thecomputation of the replicon mode as it is responsible for pos-sible RSB effects. By contrast, a zero longitudinal mode cangive information in terms of spinodal points describing howa state opens up along an unstable direction and originatesthen a saddle.The detection of a vanishing replicon mode from the high-temperature (one single equilibrium) phase is related to theappearance of marginal states . This feature is extremelyimportant because of its intimate connection with out-of-equilibrium aging dynamics.We consider then the variation of the RS action with re-spect to the overlap matrix A ( Q ab , Q aa , H a ) = − ρ σ β (cid:88) a
Gardner phase turns out to be characterized by ahierarchical organization of the different equilibria, whose an-alytical solution is well described within a Full RSB Ansatz. For the first time in an ecological context we are able to provein a rigorous way the analogy between glassy physics and com-plex energy landscapes in large ecosystems, with respect toprevious studies in this direction [20, 22]. Appendix B: Derivation of the complete phasediagram Our analysis, based on increasingly structured Ansatz, hasallowed us to derive the complete phase diagram of the Lotka-Volterra model as a function of the demographic noise. Byvarying this control parameter, we have highlighted the exis-tence of different phase transitions of increasing complexity.Thus far the only available results have been obtained withoutdemographic noise, exactly at zero temperature [19]. σ T Single equilibrium phase multiple equilibria phaseGardner phase (cid:1)(cid:2)(cid:1)(cid:2)(cid:3)(cid:2)(cid:2)(cid:4)(cid:1)(cid:2)(cid:3)(cid:2)(cid:5)(cid:1)(cid:2)(cid:3)(cid:2)(cid:5)(cid:4)(cid:1)(cid:2)(cid:3)(cid:2)(cid:6)(cid:1)(cid:2)(cid:3)(cid:2)(cid:6)(cid:4)(cid:1)(cid:2)(cid:3)(cid:2)(cid:7)(cid:1)(cid:2)(cid:3)(cid:2)(cid:7)(cid:4)(cid:1)(cid:2)(cid:3)(cid:2)(cid:8)(cid:1)(cid:2)(cid:3)(cid:9)(cid:4) (cid:1)(cid:2)(cid:3)(cid:10) (cid:1)(cid:2)(cid:3)(cid:10)(cid:4) (cid:1)(cid:2)(cid:3)(cid:11) (cid:1)(cid:2)(cid:3)(cid:11)(cid:4) (cid:1)(cid:5) (cid:1)(cid:5)(cid:3)(cid:2)(cid:4) (cid:1)(cid:5)(cid:3)(cid:5) (cid:1)(cid:5)(cid:3)(cid:5)(cid:4) (cid:1)(cid:5)(cid:3)(cid:6)(cid:5)(cid:1)(cid:12)(cid:13)(cid:14)(cid:15)(cid:16)(cid:15)(cid:17)(cid:18)(cid:15)(cid:14)(cid:19)(cid:1)(cid:1)(cid:15)(cid:20)(cid:21)(cid:22)(cid:23)(cid:17)(cid:15)(cid:16)(cid:15)(cid:22)(cid:24)(cid:1)(cid:16)(cid:15)(cid:20)(cid:12)(cid:25)(cid:14)(cid:16)(cid:22)(cid:15)(cid:26)(cid:16)(cid:12)(cid:1)(cid:12)(cid:13)(cid:14)(cid:15)(cid:16)(cid:15)(cid:17)(cid:18)(cid:15)(cid:23)(cid:1)(cid:15)(cid:20)(cid:21)(cid:22)(cid:23)(cid:17)(cid:15)(cid:16)(cid:15)(cid:22)(cid:24)(cid:1)(cid:16)(cid:15)(cid:20)(cid:12) Figure 7. Two-dimensional phase diagram of the Lotka-Volterra model in the presence of demographic noise as afunction of the heterogeneity parameter σ (in the presenceof the logarithmic term and with a selected value of the cut-off N c = 10 − ). We can distinguish three different phases: i)a single equilibrium phase above the blue line; ii) a multipleequilibria regime ( RSB stable phase) between the light blueand the orange lines; iii) a Gardner phase, characterized bya hierarchical organization of the different equilibria in thefree-energy landscape.In correspondence of the blue line in Fig. (7) the landscapestructure is no longer identified by a single equilibrium butshould be replaced by a two-level hierarchy, which results inthe appearance of new equilibrium configurations.In particular we have realized that with very low demo-graphic noise the system undergoes a Gardner transition toa new marginally stable phase: each amorphous state, saya basin, is fragmented into a fractal structure of sub-basins( metabasin ). The internal structure of states in which a basinsplits is described by the Full RSB solution of the partitionfunction, according to which the replica symmetry is brokenan infinite number of times. Similar to what observed in low-temperature glasses, such a transition will revolutionize ourunderstanding of ecosystems. The emergence of a infinite hi-erarchy of scales for the order parameter is also related to an infinite number of time-scales playing a central role in popu-lation dynamics, as we will discuss in the next Sections.What a Full RSB regime precisely yields for biological andecological settings with the appearance of an abundance ofsoft modes is particularly fascinating and still an open matterof debate. Appendix C: Results for the spin-glass modelwithout logarithmic interaction We consider exactly the same system without the loga-rithmic interaction in the species abundances. In the limit λ → + , the resulting Hamiltonian is H = (cid:88) i V i ( N i ) + (cid:88) ij α ij N i N j (C1)where the first term has the same quadratic dependence in thespecies abundances as in Eq. (A2), while the latter representsthe pairwise interacting part of the Hamiltonian. In theoreti-cal physics, Eq. (C1) is an example of spin-glass model wherethe degrees of freedom can be either discrete variables, tomodel magnetic spins, or continuous variables. The couplingsbetween spins are chosen randomly, taking both positive andnegative values and thus resulting in frustration and non-convex optimization phenomena. The simplest and most clas-sical example of spin glass is represented by the Sherrington-Kirkpatrick (SK) model [68] where each spin variable interactswith all the others in a fully-connected topology. 1. Connections with the random replicant model In theoretical ecology, a similar model was first introducedby Diederich and Opper in the 80s, studying the mean-fielddynamical evolution of S randomly interacting species witha deterministic self-interaction [34]. This pioneering modelis usually known as replicator model . Later on, a detailedstudy was also performed in [39] by using the replica formal-ism in the zero-temperature limit. Even though without thedemographic noise, this model turns out to be not partic-ularly interesting in a purely ecological context, it can stillprovide further insights on the replicator equations in manyother interdisciplinary domains [69], such as in the study ofNash equilibria [70] in Game Theory.Given S species, the Lyapunov function of the replicatormodel is written as a function of the concentration variables x i in the following way H R = S (cid:88) i 2. Exact solution in the replica symmetric case The following analysis focuses on the Lotka-Volterra Hamil-tonian, as reported in Eq. (A4), and considered in the absenceof the logarithmic term. We will derive exact expressions forthe order parameters to be analyzed as a function of temper-ature up to T = 0 . In this case, the temperature explicitlycomes out in the weighting factor exp( − βH ) and implies forthe replicated partition function: Z n = (cid:90) (cid:89) i, ( ij ) dN ai dα ij exp − (cid:88) ( ij ) ( α ij − µ/S ) σ /S − βH ( { N i } ) (C12) Exactly as for the previous analysis, we have performedboth a RS and RSB computation. In the former, the saddle-point equations read: q d = (cid:90) D z (cid:32) (cid:82) ∞ dNe − βH RS ( q ,q d ,h,z ) N (cid:82) ∞ dNe − βH RS ( q ,q d ,h,z ) (cid:33) ,q = (cid:90) D z (cid:32) (cid:82) ∞ dNe − βH RS ( q ,q d ,h,z ) N (cid:82) ∞ dNe − βH RS ( q ,q d ,h,z ) (cid:33) ,h = (cid:90) D z (cid:82) ∞ e − βH RS ( q ,q d ,h,z ) N (cid:82) ∞ dNe − βH RS ( q ,q d ,h,z ) . (C13)where now the integral over N is extended over the entirepositive axis without fixing any cut-off. D z denotes, as usual,the Gaussian integration over the auxiliary variable z – whichis introduced to decouple replicas – with zero mean and unitvariance. The RS Hamiltonian then reads: H RS ( N, z ) = − ρ σ β ( q d − q ) N ρµh − zρ √ q σ ) N + V ( N )= N (cid:2) ρ − ρ σ β ( q d − q ) (cid:3) + ( ρµh − zρσ √ q − ρK ) N . (C14)Therefore, Eqs. (C13) can be exactly computed in terms ofthe error function and its combinations, i.e. q d = (cid:90) D z (cid:40)(cid:34) e − βρ ( K − hµ + √ q zσ )22 − qd − q βρσ (cid:32) K − hµ + √ q zσ ) (cid:112) βρ [1 − ( q d − q ) βρσ ]+ e βρ ( K − hµ + √ q zσ )22 − qd − q βρσ √ π (cid:20) βρ (cid:2) ( K − hµ ) + 2 √ q z ( K − hµ ) σ + ( q − q d + q z ) σ (cid:3)(cid:21) ×× (cid:20) Erf (cid:32) βρ ( K − hµ + √ q zσ ) √ (cid:112) βρ (1 + q βρσ − q d βρσ ) (cid:33)(cid:21)(cid:33)(cid:35) √ Z (cid:41) , (C15) q = (cid:90) D z (cid:40)(cid:34) e − βρ ( K − hµ + √ q zσ )22 − qd − q βρσ (cid:32) e βρ ( K − hµ + √ q zσ )22 − qd − q βρσ √ πβρ ( K − hµ + √ q zσ ) + (cid:112) βρ (1 + q βρσ − q d βρσ )++ e βρ ( K − hµ + √ q zσ )22 − qd − q βρσ √ πβρ ( K − hµ + √ q zσ ) Erf (cid:34) βρ ( K − hµ + √ q zσ ) √ (cid:112) βρ (1 + q βρσ − q d βρσ ) (cid:35)(cid:33)(cid:35) Z (cid:41) , (C16) h = (cid:90) D z (cid:20) √ πβρ ( K − hµ + √ q zσ ) + √ e − βρ ( K − hµ + √ q zσ )22 − qd − q βρσ (cid:112) βρ (1 + q βρσ − q d βρσ )++ √ πβρ ( K − hµ + √ q zσ ) Erf (cid:20) βρ ( K − hµ + √ q zσ ) √ (cid:112) βρ (1 + q βρσ − q d βρσ ) (cid:21)(cid:21) Z , (C17)where the partition function reads: Z = √ πβρ [1 − ( q d − q ) βρσ ] (cid:20) Erf (cid:32) βρ ( K − hµ + √ q zσ ) √ (cid:112) βρ (1 + q βρσ − q d βρσ ) (cid:33)(cid:21) . (C18)Exactly as for the previous Section where we have derived theexpression for the replicon eigenvalue, we have here: λ R = ( βρσ ) (cid:104) − ( βρσ ) ( (cid:104) N (cid:105) − (cid:104) N (cid:105) ) (cid:105) , (C19)where the expressions for the averaged values (cid:104) N (cid:105) and (cid:104) N (cid:105) correspond precisely to q d and q obtained in Eqs. (C15)-(C16) according to a RS Ansatz. The overline denotes theaverage over the quenched disorder, which is technically im-plemented by integrating over the Gaussian variable z . 3. Replica symmetry broken case The computation can be easily generalized to more complexscenarios. In particular, in the RSB case the Hamiltonianreads H RSB ( N, z, t ) = N (cid:2) ρ − ρ σ β ( q d − q ) (cid:3) ++( ρµh − tρσ √ q − q − zρσ √ q − ρK ) N (C20)with the only difference of considering now the double Gaus-sian integration over z and t that account respectively for thesingle replica average and the m -block average. Even thoughthe structure of the resulting equations for ( q d , q , q , h ) , weare still able to get an exact closed-form derivation. Exactlyas before for the most general case, the saddle-point equationsread q d = (cid:104)(cid:104) N (cid:105) r (cid:105) m -r q = (cid:104)(cid:104) N (cid:105) r (cid:105) m -r q = (cid:104)(cid:104) N (cid:105) r (cid:105) m -r h = (cid:104)(cid:104) N (cid:105) r (cid:105) m -r (C21)to be solved iteratively. 4. Derivation of the phase diagram: evidence of aGardner phase By fixing the carrying capacity K = 1 and ρ = 1 in Eqs.(C15)-(C16)-(C17), we have explored the phase space as afunction of the parameter β = 1 /T along with the variationof the mean and the variance µ and σ respectively of theinteraction matrix. The equations (C13) have been solvediteratively starting from the high-temperature region up tozero temperature. For high values of σ and T there existsonly one equilibrium: the system is thus able to recover its equilibrium configuration even starting from different initialconditions. As shown in Fig. (8), upon decreasing T newphases emerge, leading respectively to an intermediate stable RSB phase and a Gardner phase in the very low-temperaturelimit.We have thus obtained the complete phase diagram at fixed µ = 10 . T σ (cid:1)(cid:2)(cid:1)(cid:2)(cid:3)(cid:4)(cid:1)(cid:2)(cid:3)(cid:5)(cid:1)(cid:2)(cid:3)(cid:6)(cid:1)(cid:2)(cid:3)(cid:7)(cid:1)(cid:8)(cid:1)(cid:8)(cid:3)(cid:4)(cid:1)(cid:8)(cid:3)(cid:5)(cid:1)(cid:8)(cid:3)(cid:6)(cid:1)(cid:2)(cid:3)(cid:5) (cid:1)(cid:2)(cid:3)(cid:6) (cid:1)(cid:2)(cid:3)(cid:7) (cid:1)(cid:8) (cid:1)(cid:8)(cid:3)(cid:4) (cid:1)(cid:8)(cid:3)(cid:5) (cid:1)(cid:8)(cid:3)(cid:6) (cid:1)(cid:2)(cid:3)(cid:4)(cid:5)(cid:6)(cid:7)(cid:8)(cid:9)(cid:4)(cid:10)(cid:4)(cid:7)(cid:11)(cid:12)(cid:1)(cid:2)(cid:13)(cid:3)(cid:4)(cid:5)(cid:6)(cid:7)(cid:8)(cid:9)(cid:4)(cid:10)(cid:4)(cid:7)(cid:11)(cid:14)(cid:15)(cid:16)(cid:3)(cid:17)(cid:18)(cid:19)(cid:20)(cid:4)(cid:21)(cid:7)(cid:4)(cid:22)(cid:5) RS phase 1RSB stable phaseFull RSB phase Figure 8. Two-dimensional phase diagram showing the insta-bility lines respectively for the one equilibrium (RS) phaseand the multiple equilibria phase ( RSB) as a function of σ (without the logarithmic term). Between the two lines (inlight blue and orange respectively) a stable RSB phase per-sists. The red dot corresponds to the analytically predictedvalue in the zero-temperature limit allowing for the estimationof the critical value σ c = 1 / √ ≈ . .In Eqs. (C11) we have derived the exact mapping betweenour model and the Random Replicant Model (RRM), pointingout how the mean and the variance of the interaction matrixare essentially related to the Lagrange multiplier ˜ γ and theparameter ˜ a ensuring the global quadratic constraint. In ourcase the phase diagram does not display any special depen-dence on the mean interaction µ and can be thus analyzedat fixed value, by changing only the heterogeneity and thedemographic noise. However, at variance with the originalRRM model, fixing µ would correspond to allowing the sumof the species concentrations to vary, i.e. injecting or remov-ing species in the ecosystem. Appendix D: Derivation of the complexity We aim to determine the number of local minima N of thefree energy corresponding to the different equilibrium config-urations in the landscape structure. The logarithm of thenumber of minima divided by S defines the complexity Σ( f ) (a.k.a. configurational entropy ) of the ecosystem. This pro-vides a crucial information to determine what kind of univer-sality class the system belongs to. In particular, in statisticalphysics of disordered systems two main universality classesare well-known: i) spin-glass models, in which the number ofminima of the free energy is not exponential in the system sizeand the free energy barriers are sub-extensive; ii) structuralglasses, for which the number of free-energy minima is actu-ally exponential in the system size. Then, there exists a finitetemperature range, between a static transition temperature T s and a dynamical transition temperature T d , within whichthe complexity turns out to be finite.As we deal with mean-field models, we can resort to a sim-plified description for which the total partition function of themodel can be expressed as the sum of α pure states contri-butions. More precisely, we can define a generic effective po-tential as a function of local order parameters whose minimaare in correspondence one-to-one with the so-called Thouless-Anderson-Palmer states [43, 71–73]. Provided the number ofminima of the free energy is exponential in the system size,the partition function can be expressed as: Z ∼ (cid:88) α e − βNf α ( T ) = (cid:90) df (cid:88) α δ ( f − f α ( T )) e − βNf = (cid:90) dfρ ( f ) e − βNf ∼ e N [ Σ( f ∗ ,T ) − βf ∗ ] n . (D1)Therefore, the configurational entropy corresponds to Σ( f, T ) ≡ S ln (cid:88) α δ ( f − f α ( T )) (D2)while f ∗ ∈ [ f , f th ] satisfies the extremal condition: ddf [ f − (1 /β )Σ( f, T )] = 0 (D3)namely d Σ( f, T ) df = 1 T . (D4)Three different regimes can be generically observed in mean-field models: • i) in the high-temperature phase (corresponding to highdemographic noise), namely for above the dynamicaltransition temperature, the paramagnetic state domi-nates the free-energy density for any value [ f , 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 paramag-netic value of the free-energy density. However, in thissecond regime, the resulting state is composed by anexponential number of metastables states of individualfree-energy density f ∗ . Upon crossing the dynamicaltemperature, the free energy preserves its analyticitywithout undergoing any true phase transition; • iii) For temperatures lower the static transition tem-perature, the leading contribution is due to the lowest free-energy states with f ∗ = f and, as a consequence, Σ( f ) = 0 . In this phase, the number of states is sub-exponential in the system size.While in the very high-temperature phase, the paramagneticsolution is always present, for T < T d , it disappears and isreplaced by a non-trivial combination of states. To explic-itly compute the complexity of the system, Σ , and graspingthe physics behind it, several techniques have been proposedin the last decades. In the following, we will focus on theso-called real replica method [44]. It consists in replicating m times the system and coupling the different copies through aninfinitesimally small parameter (cid:15) , which will be sent to zeroat the end of the computation after the thermodynamic limit.This attractive coupling naturally breaks the replica symme-try: it constrains the m copies to be in the same metastablestate yet remaining uncorrelated within a state. As a conse-quence, the free energy can be simply written as m times theindividual contribution f α . In this case: Z m ≡ (cid:88) α e − Nmβf α ( T ) = (cid:90) dfe N [Σ( f,T ) − βmf ] (D5)and, by evaluating the integral by saddle-point method, weget a similar expression to (D4) ∂ Σ( f, T ) ∂f (cid:12)(cid:12)(cid:12)(cid:12) f ∗ ( m,T ) = βm (D6)where one can immediately notice that the breaking parame-ter 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 ob-tain: Σ( f, T ) ≈ Σ( f ) + a ( T )( f − f ) + ... (D7)which implies that for m = 1 the static transition is the so-lution of the equation a ( T ) = T , whereas at small m thecondition is precisely replaced by Eq. (D6).The free-energy density of a system of m different copiescan thus be rewritten as [74] φ ( m, β ) = − βN ln Z m = min f [ βmf − Σ( f )] == βmf ∗ ( m, β ) − Σ( f ∗ ( m, T )) . (D8)Generalizing the canonical definition of entropy applied to dis-ordered systems, the complexity is thus defined as the Leg-endre transform of the free energy averaged over quencheddisorder.From the parametric plot of f ∗ ( m, β ) and its transform Σ( m, β ) , one can extract the information about the behaviorof the configurational entropy at any temperature and, con-sequently, of the associated TAP states. Therefore, Σ can beexplicitly calculated from φ ( m, β ) or, thanks to the identity: φ ( m, β ) = mF RSB (D9)from a direct computation of the RSB free energy Σ = m ddm (cid:16) βF (cid:17) == − m ddm (cid:18) n ln (cid:90) dq d dq dq dh e S A ( q d ,q ,q ,h ) (cid:19) . (D10)We compute then the derivative of the free energy w.r.t m ,which eventually leads to Σ = m ρ σ β q − q ) + (cid:90) D z ln (cid:20)(cid:90) dt a B √ π e − taB A ( z, t a B ) m (cid:21) − m (cid:90) D z (cid:82) dt aB √ π e − taB A ( z, t a B ) m ln A ( z, t a B ) (cid:82) dt aB √ π A ( z, t a B ) m (D11)where we have denoted as A ( z, t a B ) : A ( z, t a B ) ≡ (cid:90) dNe − βH RSB ( N,z,t aB ) (D12)We quantitatively evaluate the above expression within the RSB phase to determine the nature of the emerging transi-tion. We find strictly positive values of complexity at finitetemperature and fully compatible with the results previouslyobtained at zero temperature [20]. We manage then to provethat the emergent RSB stable phase is actually character-ized by a finite complexity, i.e. an exponential number ofmetastable states. Note also that, because the replicon RSBbecomes negative below the Gardner transition temperature,the complexity is well-defined only in the region for which the RSB can be safely applied, i.e. in the interval [ f , f G ] , f G being the free-energy density at the Gardner transition.Within this formalism, it is also possible to reproducethe complexity curves at different and fixed m startingfrom its equilibrium value m ∗ . The self-condition equation ∂F RSB /∂m = 0 gives indeed the equilibrium value m ∗ , whichcorresponds to the lowest free energy density and then to avanishing complexity. The resulting complexity curve is ex-pected to have an increasing trend for lower values of m up toa maximum point where unstable states start to appear anddominate the thermodynamics. Appendix E: Comparison theory and numerics1. Protocol For comparing with the theoretical results, we sample thedynamical system presented in Equation (1). The input pa-rameters of a sample are the system size S , the interactionmatrix parameters ( µ, σ ) , the immigration λ , the strength ofthe demographic noise T (temperature), and the initial con-dition distribution P [ { N i (0) } i =1 ..S ] . In order to sample onerealization 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 N i ( t = 0) from thedistribution P [ { N i (0) } i =1 ..S ] . For instance, we use afactorized uniform distribution in [0 , : P [ { N i (0) } i =1 ..S ] = S (cid:89) i =1 { N i (0) ∈ [0 , } where { . } is the indicator function.3. We sample the demographic noise { η i ( t ) } t =0 ..t max i =1 ..S , fromthe white-noise distribution, with temperature T .4. Then, all three random contributions (interactions, ini-tial conditions and demographic noise) have been dealtwith. We can then integrate deterministically the sys-tem, to end up with { N i ( t ) } t =0 ..t max i =1 ..S , where t max is thetemporal 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 Ap-pendix G. But for simplicity’s sake, let’s focus on this frame-work: we fix parameters ( S, µ, σ, λ, T ) , we sample the threerandom contributions, we integrate, and we obtain the speciespopulations over time { N i ( t ) } t =0 ..t max i =1 ..S .When we reproduce different sets of data by keeping thesame parameters, but sampling different randomness, we ob-tain { N ri ( t ) } t =0 ..t max i =1 ..S, r =1 ..N sample . 2. Observables In order to compare with the theory, we need to decideon the observables. So far, there are four sources of statis-tics in the process: the three random parts (interactions, ini-tial conditions and demographic noise) that we labelled with r = 1 ..N sample , and the species themselves i = 1 ..S . In thefollowing, we will denote E [ X ] the average over all those con-tributions. For example: E [ N ( t ) N ( t (cid:48) )] = S − N − sample S (cid:88) i =1 N sample (cid:88) r =1 N ri ( t ) N ri ( t (cid:48) ) It can be shown [48] that if the system is large enough( S (cid:29) ) and the sampling thorough enough ( N sample (cid:29) ),there is a convergence in law of the process. Mainly, there isa well defined limit ( S, N sample → ∞ ) that we can comparewith the theory. To fix ideas, we generically use S ∼ and N sample ∼ , 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 typicallynot depend on the realization of the sampling. Here for thenumerics, as (cid:28) S < ∞ , we still use some averaging over thesamples to get cleaner data.The theory is a thermodynamical one, so we will first as-sume that if we wait for a big-enough t wait , the system willreach a time-translationnally invariant (TTI) state. For in-stance, the two-time correlation C is a function of the timedifference: ∀ t ≥ t (cid:48) > t wait , E [ N ( t ) N ( t (cid:48) )] = C ( t, t (cid:48) ) = C ( t − t (cid:48) ) We check this numerically. The waiting-time depends onthe parameters, mainly ( σ, T ) . However, if we lie in thereplica-symmetric (RS) phase, we can always find the TTIstate, for rather small waiting times t wait ∼ .All the comparisons we will be making are in this state( t ≥ t wait ) , for the RS phase. We will now use a mappingbetween thermodynamics properties, and dynamical ones. Weuse the notations from Equations (A13). h = E [ N ( t )] q d = C (0) = E (cid:2) N ( t ) (cid:3) q = lim τ →∞ C ( τ ) = lim τ →∞ E [ N ( t ) N ( t + τ )] ∼ E [ N ( t ) N ( t max )] The lhs is predicted by the theory, and the rhs are numer-ical observables. 3. Example of numerical results in the RS phase On figure 9, we show that one time observables such as E [ N ( t )] or E (cid:2) N ( t ) (cid:3) converge to a constant value in time.This indicates the reach of a TTI state. It can be confirmedby the collapse of two-time observables such as the correla-tion E [ N ( t ) N ( t (cid:48) )] = C ( t, t (cid:48) ) , that we plot as C ( t − t (cid:48) , t (cid:48) ) fordifferent t (cid:48) on figure 10.We can see that for t > t wait ∼ here, the system isindeed TTI, at least regarding these observables. We thenread the values of h = E [ N ( t )] TTI and q d = E (cid:2) N ( t ) (cid:3) TTI when they stabilize. And we read q = E [ N ( t ) N ( t max )] TTI on the collapse of figure 10.We can also infer a relevant information from figure 10: thetimescale for decorrelation τ decorrel . We are only interestedin the scaling of this observable, so we introduce a roughestimate. We approximate τ decorrel by the needed time so thatthe decorrelation decay is of . Mathematically, τ decorrel is then determined by: C ( τ decorrel ) − C ( ∞ ) = 0 . C (0) − C ( ∞ )) 4. Match in the RS phase The results are presented on figure 11: the theory matchesbeautifully the numerics. time o n e t i m e o b s e r v a b l e s [ N ( t )] h [ N ( t ) ] q d Figure 9. RS one time observables converge in time. Pa-rameters are ( S, µ, σ, λ, T ) = (500 , , , − , − ) . Thisdata comes from only one sample of the process, with dis-crete timestep dt = 10 − . The dashed lines correspond to theread TTI value of h and q d . time t t c o rr e l a t i o n [ N ( t ) N ( t )] t = 0 t = 1 t = 10 t = 100 t = 500 q Figure 10. RS correlation E [ N ( t ) N ( t (cid:48) )] = C ( t, t (cid:48) ) , plottedas a function of t − t (cid:48) for different t (cid:48) . This data is from thesame sample as figure 9. We see that, up to fluctuations, thecorrelation collapse as a function of t − t (cid:48) for t (cid:48) > t wait ∼ .The dashed line correspond to the read TTI value of q . Here,the timescale for decorrelation is around τ decorrel ∼ .On figure 12, we show how the timescale for decorrelation τ decorrel diverges as we approach the 1RSB transition fromabove in temperature. Data seems to indicate a critical ex-ponent as τ decorrel ( T ) ∼ ( T − T RSB ) − / . 5. Rough results in the 1RSB phase In the 1RSB phase, thermodynamics indicate that the sys-tem no longer reaches a TTI state. Instead, it presents agingbehaviour: the older the system is, the slower it becomes. In h q d Temperature q T RSB theorynumerics Figure 11. Comparison with the theory in the RS phase.Parameters are ( S, µ, σ, λ ) = (500 , , , − ) . We considerthe observables ( h, q D , q ) as a function of temperature. Theorange full line is the theory predictions. Blue crosses arenumerical results, error bars are taken with respects to the N sample = 50 different samples of the ecosystem. We found t wait ∼ to be enough to observe TTI state in all thesevalues of temperature, except for the last point on the left( T = 2 . ): due to slowing down of the dynamics, we had toincrease the extent of the simulation and found t wait ∼ .the orange dashed line correspond to the critical temperatureat which the theory becomes 1RSB. Indeed, numerically wecan’t observe TTI state below this temperature, even increas-ing t max to .the simplest case of aging, there is a good understanding ofthe correlation decay C ( t, t (cid:48) ) . At equal time, the correlation isthe dynamical one C ( t (cid:48) , t (cid:48) ) = q d . Then it decorrelates quicklyfor t > t (cid:48) to an intermediate plateau C ( t, t (cid:48) ) = q < q d , asthe system explores the neighbouring phase space. Eventu-ally, for t (cid:29) t (cid:48) , the correlation decreases to a final plateau C ( t, t (cid:48) ) = q < q . The timescale from the intermediateplateau to the final plateau increases with the age t (cid:48) of thesystem. These theoretical predictions are compared with nu-merical results on figure 13. Appendix F: Mathematical issues for immigrationimplementation In this part, we detail the mathematical issue for immi-gration implementation. This problem is independent on theinteractions, so we drop them ( α = 0 here). We consider thefollowing one-species Ito-stochastic process: dNdt = N (1 − N ) + √ T N η + λI ( N ) (F1)with white noise (cid:104) η ( t ) η ( t (cid:48) ) (cid:105) = δ ( t − t (cid:48) ) , and immigration Temperature T T RSB d e c o rr e l ( T ) -0.56 slope Figure 12. Critical slowing down of the dynamics. We plot thedecorrelation time τ decorrel ( T ) as a function of T − T RSB , in aloglog scale. Blues crosses come from the same numerical dataas figure 11. Red dashed line is a simple fit. As we approachthe transition, the system becomes slower and slower, and thedynamical timescale diverges. time t-t' r e s c a l e d c o rr e l a t i o n C ( t , t ) / C ( t , t ) t'=1e+02t'=1e+03t'=1e+04t'=1e+05t'=1e+06 q / q d q / q d Figure 13. 1RSB aging. Parameters are ( S, µ, σ, λ, T, N sample ) = (2000 , , , − , / , . Weplot the rescaled correlation C ( t, t (cid:48) ) /C ( t (cid:48) , t (cid:48) ) as a function of t − t (cid:48) , for different t (cid:48) . 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 intermerdiateand final plateau values. They do not coincide exactly withthe data, but the trends correspond.function I ( N ) . Immigration is generically implemented sothat the populations do not go too close to 0. In the usualimmigration, I ( N ) = 1 , but we will see that this is problem-atic.We want to have a hint at the stationary distribution ofpopulation P ∞ ( N ) = P ( N, t = ∞ ) . In order to obtain it, wechange variables so that the noise becomes additive, and notmultiplicative any more. Here the relevant change of variablesis s ( t ) = √ N ( t ) , and Ito’s lemma gives: dsdt = s − s + λI ( N )2 s − T / s + (cid:112) T / η Then we use Langevin-Boltzmann to read the stationarydistribution.In the usual immigration I ( N ) = 1 case, the stationarydistribution is always integrable: P ∞ ( N ) = Z − N λ/T − exp 1 T ( N − N However, the corresponding effective potential V eff ( N ) be-haves repulsively around N = 0 only if λ > T : V eff ( N ) = − N + N − ( λ − T ) ln N Indeed, if we introduce an approximate induced cut-offvalue N cut ( b ) such that P ∞ ( N < N cut ) ∼ e − b (cid:28) , the scal-ing yields N cut ( b ) ∼ e − bT/λ , which means that the density isstill relevant up to e − bT/λ (cid:28) .Basically, this means that demographic noise with usual im-migration will not prevent populations from reaching very lowvalues. The usual immigration is not strong enough when fac-ing demographic noise. This is indeed problematic, becausewhenever we will want to actually compute observables, theintegrals will be dominated by the domain N ∼ . Thisis wrong physically (important species should be the highpopulation ones), and difficult numerically (integration is ill-defined).In order to solve this, we can implement stronger immi-gration such as I ( N ) = N − α with α > . However, anothereven simpler physical solution is to impose a hard repulsiveboundary condition on the problem: an infinite potential at N = λ . In this case, the same steps can be performed and weobtain the stationary distribution: P ∞ ( N ) = Z − N − exp (cid:20) T ( N − N (cid:21) { N > λ } which is well-behaved. This is the solution we chose forboth the theory predictions and the numerics. We will nowdetail in the next section how to integrate this process nu-merically. Appendix G: Numerical scheme to sampledemographic noise1. Litterature review Numerical simulations need discrete time. However,when discretizing time with bounded random processes,one often encounters a non-zero probability that during onetime-step the system will cross the boundary of the system(for example the N ≥ boundary in our case), and becomenumerically unstable. We review different solutions thathave been proposed to solve this issue, and check how theydeal with our Lotka-Volterra (LV) system. A more thoroughreview can be found in [42].A first naive way to go around the difficulty is to changevariable (sqrt, ln...). But this won’t work because if the noise becomes treatable, the deterministic part becomesnumerically unstable. Most articles then study the numericalintegration of processes such as ˙ N = α + βN + √ σN η .[75] proposes Balanced Implicit Method: they implementa clever discretization scheme so that the boundaries (pos-itivity) are respected. The scheme amounts to Euler’s forsmall time step. It needs a small regularization. It doesnot work for LV, because it needs very small regularizationparameter and time-step to give good results. This is tooheavy numerically.[76] derives the exact Fokker Planck solution of a simplersystem. But sampling is inefficient (rejection method). [77]builds on this method by improving the sampling method, butthis still isn’t satisfactory. Eventually, [78] improves again themethod, by exactly solving (Fokker-Planck) the full process.The sampling is clever, with Poisson variables. They alsoindicate a way to solve more elaborate processes, which wewill detail in the following section. Our strategy is heavilybased on [78]. 2. Our implementation The idea from [78] is to separate the process into solvableones. More precisely, we want to solve: ˙ N i = √ N i η i − N i − N i (cid:16)(cid:88) α ij N j − (cid:17) ... + hardW all ( λ ) (G1)where hardW all ( λ ) implements the hard wall boundary at N = λ . We will discretize time with a timestep dt , and furthersubdivise it into three timesteps dt (cid:48) = dt/ . We consider thatonly one part of the process is active during a subtimestep dt (cid:48) . So the final scheme is the following:1. From [76], we know how to sample efficiently the demo-graphic noise only ˜ N i ( t + dt (cid:48) ) = Gamma (cid:20) Poisson[ N i ( t ) T dt (cid:48) ] (cid:21) T dt (cid:48) This corresponds to a process which only feels the de-mographic noise ˙ N i = √ N i η i during [0 , dt (cid:48) ] . We re-spectively used the notation Poisson[ ω ] ( Gamma[ ω ] ) forrandom Poisson (Gamma) variables, with parameter ω .2. Treating immigration as a reflecting wall. The particlewishes to go to ˜ N i but bounces on the wall. N i ( t + dt (cid:48) ) = λ + | ˜ N i ( t + dt (cid:48) ) − λ | 3. During [ dt (cid:48) , dt (cid:48) ] , only integrate the blue process ˙ N i = − N i : N i ( t + 2 dt (cid:48) ) = N i ( t + dt (cid:48) )1 + dt (cid:48) N i ( t + dt (cid:48) ) 4. During [2 dt (cid:48) , dt (cid:48) ] , only integrate the pink process ˙ N i = − N i ( (cid:80) α ij N j − : N i ( t + 3 dt (cid:48) ) = N i ( t + 2 dt (cid:48) ) exp dt (cid:48) (cid:16) − (cid:88) α ij N j ( t ) (cid:17) There are a lot of different combinations of this kind ofschemes. We tried some, and chose this one after a lot ofchecks on simpler models for which we know the distributionsat all times. 3. Issues of our implementation After careful tests on simpler models, we used this schemeto compare with the theory. Initially we used a hardwallimmigration at λ = 10 − . The agreement was quite goodfor second degree observables ( q d , q ), but not for h . Thisis due to the numerical scheme. Indeed, if T is quite high( T (cid:29) λ ), the sampling of the demographic noise sends many O (1) species close to 0, then they bounce on the wall and endup at N = 2 λ . Therefore there is an induced concentrationof species at N = 2 λ . Because of the λ peak, there is asubsampling of the O (1) populations.In order to reduce this issue, we use a higher λ = 10 − in the final results that are shown on figure 11. We reckonthe slight discrepancy at high temperature between theoryand numerics comes from this issue. A solution is still underinvestigation in Appendix H . We are aware that the methodis still in development. However, it is already enough at themoment to beautifully confirm the theory. Appendix H: Ongoing investigations1. χ A cleaner numerical test for the transition RS to 1RSBwould be the divergence of the χ correlation. So far, we donot have enough data to present clean results, but in principlethis observation should not be too difficult. 2. Improve the numerical scheme In the current numerical scheme, we first sample pure de-mographic 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 thiswould be to directly solve the Fokker-Planck equation associ-ated to the whole process "demographic noise + hard wall".We reckon this can be done adapting the proof from [78].[1] R. May, A. R. McLean, et al. , Theoretical ecology: prin-ciples and applications (Oxford University Press on De-mand, 2007).[2] K. Faust and J. Raes, Nature Reviews Microbiology ,538 (2012).[3] V. Bucci and J. B. Xavier, Journal of molecular biology , 3907 (2014).[4] R. M. Goodwin, Chaotic economic dynamics (OxfordUniversity Press, 2003).[5] D. A. Kessler and N. M. Shnerb, Physical Review E ,042705 (2015).[6] D. S. Maynard, Z. R. Miller, and S. Allesina, NatureEcology & Evolution , 91 (2020).[7] J. A. Chandler, J. M. Lang, S. Bhatnagar, J. A. Eisen,and A. Kopp, PLoS genet , e1002272 (2011).[8] J. Lloyd-Price, A. Mahurkar, G. Rahnavard, J. Crabtree,J. Orvis, A. B. Hall, A. Brady, H. H. Creasy, C. Mc-Cracken, M. G. Giglio, et al. , Nature , 61 (2017).[9] S. F. Risk, OECD Economic Outlook (2012).[10] S. Thurner, IFP/FGS Working Paper (2011).[11] L. Barreira and C. Valls, in Dynamical Systems (Springer, 2013) pp. 57–86.[12] R. MacArthur, Theoretical population biology , 1(1970).[13] D. Tilman, Resource competition and community struc-ture (Princeton university press, 1982).[14] S. Ruan, in Delay differential equations and applications (Springer, 2006) pp. 477–517.[15] J. Vano, J. Wildenberg, M. Anderson, J. Noel, andJ. Sprott, Nonlinearity , 2391 (2006).[16] J. H. van Opheusden, L. Hemerik, M. van Opheusden,and W. van der Werf, SpringerPlus , 474 (2015).[17] C. K. Fisher and P. Mehta, Proceedings of the NationalAcademy of Sciences , 13111 (2014).[18] C. A. Serván, J. A. Capitán, J. Grilli, K. E. Morri-son, and S. Allesina, Nature ecology & evolution , 1237(2018).[19] G. Bunin, Physical Review E , 042414 (2017). [20] G. Biroli, G. Bunin, and C. Cammarota, New Journalof Physics , 083051 (2018).[21] M. Tikhonov and R. Monasson, Physical review letters , 048103 (2017).[22] A. Altieri and S. Franz, Physical Review E , 010401(2019).[23] M. T. Pearce, A. Agarwala, and D. S. Fisher, Proceed-ings of the National Academy of Sciences (2020).[24] F. Roy, M. Barbier, G. Biroli, G. Bunin, et al. , PLOSComputational Biology , 1 (2020).[25] R. Marsland, W. Cui, and P. Mehta, Scientific reports , 1 (2020).[26] L. Sidhom and T. Galla, Physical Review E , 032101(2020).[27] I. Dalmedigos and G. Bunin, arXiv preprintarXiv:2002.04358 (2020).[28] J. Moran and J.-P. Bouchaud, Physical Review E ,032307 (2019).[29] Along a similar vein, models of liquids and crystals thatare used in physics are disparate, approximate and ofteninaccurate with respect to real system. Yet, the proper-ties of the phases that arise from their studies provide aprecise and quantitative description of the phases foundin nature.[30] Y. V. Fyodorov and B. A. Khoruzhenko, Proceedings ofthe National Academy of Sciences , 6827 (2016).[31] Y. V. Fyodorov and P. Le Doussal, Journal of Physics A:Mathematical and Theoretical , 474002 (2018).[32] Y. V. Fyodorov, G. Ben Arous, and B. A. Khoruzhenko,arXiv preprint arXiv:2008.00690 (2020).[33] M. Barbier, J.-F. Arnoldi, G. Bunin, and M. Loreau,Proceedings of the National Academy of Sciences ,2156 (2018).[34] S. Diederich and M. Opper, Physical Review A , 4333(1989).[35] T. Galla and J. D. Farmer, Proceedings of the NationalAcademy of Sciences , 1232 (2013).[36] J. B. Sanders, J. D. Farmer, and T. Galla, Scientificreports , 1 (2018). [37] S. Solomon et al. , Advances in Complex Systems (ACS) , 301 (2000).[38] J. Moran and J.-P. Bouchaud, Physical Review E ,032307 (2019).[39] P. Biscari and G. Parisi, Journal of Physics A: Mathe-matical and General , 4697 (1995).[40] G. Domokos and I. Scheuring, Journal of Theoretical Bi-ology , 535 (2004).[41] T. Rogers, A. J. McKane, and A. G. Rossberg, EPL(Europhysics Letters) , 40008 (2012).[42] H. Weissmann, N. M. Shnerb, and D. A. Kessler, Phys-ical Review E , 022131 (2018).[43] M. Mézard, G. Parisi, and M. Virasoro, Spin glass the-ory and beyond: An Introduction to the Replica Methodand Its Applications , Vol. 9 (World Scientific PublishingCompany, 1987).[44] R. Monasson, Physical review letters , 2847 (1995).[45] We focus on the ones giving the leading contribution tothe partition function; considering a different value wouldjust slightly shift the transition line but keeps qualita-tively unaltered the conclusions.[46] P. Charbonneau, J. Kurchan, G. Parisi, P. Urbani, andF. Zamponi, Nature communications , 1 (2014).[47] L. Berthier, G. Biroli, P. Charbonneau, E. I. Corwin,S. Franz, and F. Zamponi, The Journal of chemicalphysics , 010901 (2019).[48] G. Ben Arous, A. Dembo, and A. Guionnet, ProbabilityTheory and Related Fields , 619 (2006).[49] H. Sompolinsky and A. Zippelius, Physical Review B ,6860 (1982).[50] S. Franz and M. Mézard, Physica A: Statistical Mechan-ics and its Applications , 48 (1994).[51] L. F. Cugliandolo and J. Kurchan, Physical Review Let-ters , 173 (1993).[52] L. F. Cugliandolo, in Slow Relaxations and nonequilib-rium dynamics in condensed matter (Springer, 2003) pp.367–521.[53] G. Biroli, Journal of Statistical Mechanics: Theory andExperiment , P05014 (2005).[54] J. Kurchan and L. Laloux, Journal of Physics A: Math-ematical and General , 1929 (1996).[55] L. F. Cugliandolo, J. Kurchan, and L. Peliti, PhysicalReview E , 3898 (1997).[56] G. Parisi, Proceedings of the National Academy of Sci-ences , 7948 (2006). [57] W. Mende, The predator-prey model: do we live in aVolterra world? (Akademie-Verlag, 1986).[58] J. M. Smith and J. M. M. Smith, Evolution and the The-ory of Games (Cambridge university press, 1982).[59] K. Sigmund et al. , Evolutionary Game Dynamics: Amer-ican Mathematical Society Short Course, January 4-5,2011, New Orleans, Louisiana , Vol. 69 (American Math-ematical Soc., 2011).[60] L. Berthier and J. Kurchan, Nature Physics , 310(2013).[61] J. Hertz, G. Grinstein, and S. Solla, in AIP ConferenceProceedings , Vol. 151 (American Institute of Physics,1986) pp. 212–218.[62] J. Hertz, G. Grinstein, and S. Solla, in Heidelberg collo-quium on glassy dynamics (Springer, 1987) pp. 538–546.[63] A. Altieri, in Jamming and Glass Transitions (Springer,2019) pp. 133–152.[64] A. J. Lotka, Proceedings of the National Academy of Sci-ences , 410 (1920).[65] V. Volterra, Variazioni e fluttuazioni del numerod’individui in specie animali conviventi (C. Ferrari,1927).[66] G. Parisi, Physical Review Letters , 1946 (1983).[67] A. Bray and M. Moore, Journal of Physics C: Solid StatePhysics , 79 (1979).[68] S. Kirkpatrick and D. Sherrington, Phys. Rev. Lett ,1792 (1975).[69] R. Cressman and Y. Tao, Proceedings of the NationalAcademy of Sciences , 10810 (2014).[70] J. Berg and M. Weigt, EPL (Europhysics Letters) ,129 (1999).[71] D. J. Thouless, P. W. Anderson, and R. G. Palmer,Philosophical Magazine , 593 (1977).[72] A. Crisanti and H.-J. Sommers, Journal de Physique I ,805 (1995).[73] A. Cavagna, I. Giardina, and G. Parisi, Physical ReviewB , 11251 (1998).[74] To be more definite, the exact protocol requires the com-putation of φ ( m, β ) = − βN lim n → ∂ n ( Z m ) n , reproducing asystem of ( m × n ) copies and eventually taking the limit n → .[75] G. N. Milstein, E. Platen, and H. Schurz, SIAM Journalon Numerical Analysis , 1010 (1998).[76] L. Pechenik and H. Levine, Physical Review E , 3893(1999).[77] E. Moro, Physical Review E , 045102 (2004).[78] I. Dornic, H. Chate, and M. A. Munoz, Physical ReviewLetters94