Many-species ecological fluctuations as a jump process from the brink of extinction

Highly-diverse ecosystems exhibit a broad distribution of population sizes and species turnover, where species at high and low abundances are exchanged over time. We show that these two features generically emerge in the fluctuating phase of many-variable model ecosystems with disordered species interactions, when species are supported by migration from outside the system at a small rate. We show that these and other phenomena can be understood through the existence of a scaling regime in the limit of small migration, in which large fluctuations and long timescales emerge. We construct an exact analytical theory for this asymptotic regime, that provides scaling predictions on timescales and abundance distributions that are verified exactly in simulations. In this regime, a clear separation emerges between rare and abundant species at any given time, despite species moving back and forth between the rare and abundant subsets. The number of abundant species is found to lie strictly below a well-known stability bound, maintaining the system away from marginality. At the same time, other measures of diversity, which also include some of the rare species, go above this bound. In the asymptotic limit where the migration rate goes to zero, trajectories of individual species abundances are described by non-Markovian jump-diffusion processes, which proceeds as follows: A rare species remains so for some time, then experiences a jump in population sizes after which it becomes abundant (a species turnover event) and later sees its population size gradually decreasing again until rare, due to the competition with other species. The asymmetry of abundance trajectories under time-reversal is maintained at small but finite migration rate. These features may serve as fingerprints of endogenous fluctuations in highly-diverse ecosystems.


I. INTRODUCTION
In ecological communities, interactions between species can drive changes in population sizes.For few-species communities, experiments find dynamics including stable equilibria, periodic oscillations and chaos [1][2][3], which are explained in terms of dynamical models of interacting populations, such as Lotka-Volterra or resourcecompetition models [4].Yet many natural ecosystems, from microbes in a grain of soil to plants in a rainforest, can be staggeringly diverse.Nevertheless, the dynamics of such highly-diverse communities are far less understood.
Observations on highly diverse systems show that the distribution of species abundances (population sizes) at a given time is often very broad, with many species at very low population size [5,6].Time fluctuations in abundances can be very large, with "blooms" and significant species turnover, where the species at high and low abundances are exchanged over time [7][8][9].
Theoretically, dynamical models of interacting populations with many variables can be notoriously challenging to analyze.They are parameterized by very many parameters describing the interactions between species, which are unknown and often unrealistic to obtain from measurements.This has prompted a change of paradigm (following similar ideas in physics and other fields), replacing unknown parameters by randomly sampled ones [10,11], and looking for typical and universal properties of the many-variable systems.Both physicists and ecologists aim to classify the different broad behaviors and the robust features of each, formalized within physics as "phases".An important contribution of statistical physics is the ability to provide mathematical frameworks, giving systematic answers to questions on these robust properties.This work provides such a framework for one such challenging phase.
Two distinct phases that have recently attracted much attention, are a phase where the abundances of different species reach a fixed point, and another where they fluctuate indefinitely [12,13].These distinct behaviors have been observed in controlled experiments where properties of microbial communities are varied [14], highlighting the power of robust theoretical predictions when applied to ecological phenomena.
Most of the theoretical work so far has been devoted to the fixed point phase, and many of its properties are well understood, including the abundance distribution, and limits on the fixed point stability that signal the transition to the fluctuating phase [12,13,15].Predictions obtained theoretically for the fixed-point phase stand in contrast with empirically observed broad abundance distributions, and are also unable to account for situations where large abundance fluctuations are observed.The dynamical phase, which is the focus of this work, holds promise of addressing these limitations.
Much less has been known about the dynamical phase.For well-mixed systems (no explicit space), that are coupled to the outside by a migration of all species, simulations have shown that the system reached a stationary chaotic state [16].Similar results have been obtained from simulations of coupled spatial locations [17,18].However, for well-mixed systems in the absence of migration, a dynamical slowdown is observed, along with large population fluctuations [16,18].Analytical results for many-species fluctuating dynamics have been derived when interactions are fully anti-symmetric, in which case a stationary chaotic state is reached even without migration [18].Yet, this state is sensitive to the anti-symmetry that is not expected to hold generally in nature [18].A few-species model featuring dynamical slowdown and large fluctuations is the three-species rock-paper-scissors dynamics, cycling between and ever-closer to three unstable fixed points (a heteroclinic orbit) [19].This elegant model serves as an instructive analogy for many-species dynamics [18,20,21], yet is only of limited relevance to many-species properties such as diversity, abundance distributions, and stability.A many-species exactly-solvable dynamical toy model that features dynamical slowdown was introduced in [21].Yet due to its special structure, this model did not include key features of ecological systems.In particular, questions relating diversity and linear stability that appear in many other systems cannot be addressed, and it also did not include migration that interrupts the dynamical slowdown process.
In this work, we provide a systematic analytical framework for the dynamical phase, for the Lotka-Volterra model with randomly sampled interactions, in the limit of many species, and when migration rates are small.In this phase, fluctuations in population sizes are caused solely by interactions between species, without changes in the environment.We show that it is precisely in the limit of low migration, where fluctuations in population sizes get slower and larger, that many of the striking signatures of this phase emerge.These features, which are now listed, include commonly observed traits of highdiversity natural ecosystems, such as a broad abundance distribution and species turnover, see points 1 and 2 below.In addition, they shed new light on the definition and measurement of diversity when species turnover is involved, see point 4, and uncover new phenomena that could be used as fingerprints of endogenous fluctuations, see point 5.
1.At any given time, abundances are broadly distributed, with many rare species with population size close to the minimal value set by the migration and many species with small population sizes but yet much larger than the minimal one.We show analytically that over this intermediate range, the abundance distribution scales as a power law with exponent −1 and characterize the corrections to this power law behavior when the migration rate is finite.Broad distributions have been measured in natural ecosystems [6,22], but their form and origin have been debated [5,23].In contrast, in the fixed point phase, the number of species in this intermediate range is small and the abundances are not broadly distributed [13].
2. The dynamics exhibit species turnover, where the species at high abundance are exchanged over time with species that were at low abundances.In particular, we show that there are always species at low abundances that are able to grow.This is in contrast to the fixed point phase, where rare species cannot invade.

3.
A long timescale emerges, possibly extending over many generations: as the migration gets lower, temporal changes in both abundances and growth rates become slower.This timescale scales as the absolute value of the logarithm of the migration rate, a prediction that might be directly tested in controlled experiments.Understanding the various timescales involved is key to understanding ecological dynamics [24] and this work identifies a new collective mechanism through which long timescales are robustly self-generated by species interactions in highly-diverse model ecosystems.The combination of species that can invade and long timescales allows even for species that are very rare to reach high abundance later.Hence, rare species at a given time can be important in the future, and the distinction of rare versus abundant species is not fixed in time.Due to the long timescale, the species at high abundance lie close to a fixed point, which would have been stable in the absence of the other species.
4. How do many species coexist in highly-diverse ecosystems is one of the key questions in ecology.Turnover events and broad abundance distributions raise questions as to how one might even define diversity.We show that while the list of abundant species changes in time, the fraction above any given threshold abundance fluctuates only by a little.Furthermore, in the low migration limit, the number of species with high and intermediate abundance is well-defined, namely insensitive to the precise choice of the thresholds.Finally, we show that the number of high-abundance species, which lie close to a fixed point, is strictly below the stability bound (known as the "May bound" [11]).Equivalently, this fixed point is fully stable (as opposed to marginally-stable).In contrast, the number of high-abundance and intermediateabundance species together is not constrained by the stability bound and exceeds it.Therefore measures of species-richness that also probe the power law region of the species abundance distribution are not constrained by May's stability bound.In fact, most of the species whose population sizes are well above the minimal value set by migration, have small abundances.
5. Lastly, we show that "bloom" dynamics, where a species grows from rare to abundant, until eventually going back to rare, are strongly asymmetric under time-reversal: The trajectory of the population size starts with a quick increase and decreases back gradually.This is an interesting signature of purely endogenous fluctuations in high-diversity ecosystems that could be looked for in empirical time series.
In addition to the above predictions, we also consider isolated systems (namely, without migration from the outside) and show that there, the timescale and size of the abundance fluctuations will continue to grow in time indefinitely.This will inevitably lead to extinctions of many species, once the finite size of the populations is taken into account.The core of our argument is built on identifying the appropriate transformations of time and abundances, for which all dynamical properties (including steady-state distributions and two-time correlations) collapse for different values of migration, when migration is low.These scaling relations are verified exactly in collapse of simulation data.We obtain a well-defined stochastic process for these transformed variables in the limit of small but positive migration rate.The resulting picture has features that set it apart from generic many-variable chaotic dynamics.First, abundances follow non-Markovian jumpdiffusion dynamics.Second, even though the dynamics are slow, the system is maintained strictly away from marginality, in contrast with glassy systems [25,26].Last, the emergence of the long timescale results from dynamical slowdown (aging) that would have continued indefinitely at zero migration.Yet, this slowdown is not due to the existence of a rough landscape.We trace all the unique phenomenology back to the possibility of extinctions in the absence of migration.In other words, this is a consequence of the multiplicative nature of the dynamics, which are constrained to positive population sizes with absorbing boundaries at zero population sizes.

II. MODEL DEFINITION, THE TWO PHASES
We start with the Lotka-Volterra system of equations for i = 1 . . .S where S is the total number of species.The variables N i represent the abundances (population sizes) of the different species, and so N i ≥ 0 at the initial time and is guaranteed to remain so throughout.λ ≥ 0 represents migration from an external source, which for simplicity is taken to be the same for all i.Eq. ( 1) is a standard rescaling [13,[27][28][29] of the Lotka- For simplicity, we take all r i = r.Eq. ( 1) is obtained by rescaling time by r, and introducing the nondimensional parameters ). Below we study the limit λ ≪ 1, but assume that absolute population sizes n i , whose minimal size is around D i /r, are always large enough so that demographic stochasticity can be neglected.The large parameter |ln λ| will play an important role in the following, and since the absolute population sizes k i can be large (e.g., up to 10 10 bacterial cells in one milliliter), (C-F) The corresponding distributions for N and ln N , denoted by P (N ) and P (ln N ) respectively.The distributions contain three parts (see E,F): a top part, which remains O λ 0 when λ → 0 + , containing S * top species; an intermediate part, at λ ≪ N ≪ 1; and a low part at O(λ).The intermediate part contains a finite fraction of the species in the dynamical phase, but not in the fixed point phase.(G) The three definitions of species richness, normalized by var(α) so that 1 is the fixed-point stability bound, as a function of the interaction variability σ.For σ < σc a fixed point is reached, all three definitions coincide and lie below the stability bound, and agree with known theory (solid line).At σ = σc the stability bound is reached.Beyond it, there are persistent fluctuations, and the three definitions no longer coincide: S * top is lower than the stability bound, while S * inter and S * growth are above it.The diversities are obtained by solving the rescaled dynamics defined in Sec.IV, and agree quantitatively with careful analysis of the full equations of motion (see Sec. VII).Simulation parameters for all figures are given in Appendix A 6.
|ln λ| can be reasonably large in ecologically relevant set-tings.
The limit of large S is relevant to many natural highdiversity ecosystems, with hundreds to thousands of species, in communities from microbes to trees [5,30].We assume that α is a random matrix with Gaussian entries, and for mathematical convenience we carry most of the analysis in the case where the interaction coefficients are all sampled independently from each other such that ⟨α ij ⟩ = µ/S and ⟨α ij α kl ⟩ − ⟨α ij ⟩⟨α kl ⟩ = σ 2 δ ik δ kl /S.By using combinations of numerical simulations and analytical calculations, we later show that the qualitative picture presented in this paper is nonetheless robust to the addition of correlations (positive or negative) of the interaction coefficients within pairs of species, see Sec.VII .
The system exhibits different dynamical behaviors depending on the parameters [13].When µ > 0 and the heterogeneity of interactions is smaller than a critical value σ < σ c the dynamics reach a fixed point, see Fig. 1(A).In it, some species are absent and others remain present.This definition is straightforward when λ = 0, where a fixed point dN i /dt = 0 in Eq. ( 1) implies either N i = 0 that are the absent species, or N i > 0 which are the present species.At the stable fixed point reached, the extinct species have negative growth rates Ṅi /N i < 0 and so cannot invade, and the subset of present species is linearly stable.For small λ > 0, the absent species are now at values N i of order λ and the present species are unaffected by the small λ.The requirements for a stable fixed point remain the same.
Above σ > σ c , the system evolves indefinitely, see Fig. 1(B), without ever settling at a stable fixed point (in fact, such stable, uninvadable fixed points do not exist [31]).

III. PHENOMENOLOGY
We now describe key features of the dynamics when σ ≥ σ c and λ ≪ 1, starting with the Species Abundance Distribution (SAD), and then turn to the dynamics.

A. Species abundance and stability
The species abundance distribution P (N ) in the dynamically-fluctuating phase is shown in Fig. 1(D,F).It spans many orders of magnitude, ranging from O(1) population sizes to order O(λ) populations sizes.This is indeed the order of magnitude of the minimal value allowed by migration, which we call the migration floor in the following.Compared to the equilibrium situation, Fig. 1(C,E), there are many more species with abundances in the intermediate range λ ≪ N ≪ 1 and the fraction of species there remains finite even for small λ, see Fig. 1(F).In this range, P (N ) appears to be approximately a power law.The value of P (ln N ) changes slowly, so the distribution of the abundances, P (N ), is expected to behave roughly as N −1 [18,32], see Fig. 1(C).In Sec.
V B, we show that in the λ → 0 + limit the power law is indeed exactly N −1 , and refine this picture with precise corrections to this N −1 behavior.
A common way to define the number of present species (the "species richness") is by those whose abundance lies above some value.We consider two definitions based on this idea, and an additional criterion based on the invasion growth rate.The three proposed definitions of species richness are as follows: 1. S * top is the number of species belonging to the right peak in Fig. 1(F).Their abundance N i remains finite even for small λ.These are the top or abundant species.

S *
inter includes the number of "intermediate" species.It is all species except those belonging to the left peak in Fig. 1(F).They satisfy N i ≫ λ.

A third definition, S *
growth , can be obtained through invasion experiments.If the population size of species i is set to a value that is small but wellabove the migration floor, λ ≪ N new i ≪ 1, and keeping all other N j unchanged, species i will grow with Ṅi /N i invasion = g i > 0. g i is called the "invasion growth rate" (see, e.g.[33,34]).For the dynamics Eq. ( 1), growth counts the number of species with positive growth rate.
At the fixed point phase, all three definitions coincide, S * top = S * inter = S * growth = S * , when λ is small.In the fluctuating phase, given the presence of species in the range λ ≪ N ≪ 1, one might worry that species richness S * top is not well-defined, in that it relies on an arbitrary cutoff on the abundances to decide which species are "present" or "absent", with a similar concern for S * inter .In the limit λ → 0 + , we show below that the peaks in see Fig. 1(F) become narrow compared to |ln λ|, and therefore the definitions for S * top , S * inter become sharp.Focusing on S * top , we argue in Sec.VII that for λ small but reasonable for ecological applications, these quantities can be measured in practice, and are not far from their asymptotic values.
The values of S * top , S * inter , S * growth are plotted in Fig. 1(G).They are compared to a bound on species richness coming from linear stability (horizontal dashed line): If a subset of S * abundant species is at a fixed point, meaning that the abundances satisfy 1 − N i − j(̸ =i) α ij N j = 0, this fixed point will typically be linearly stable (ignoring other species) if S * var(α) ≤ 1 [11,13].In the fixed point phase, it is known that S * var(α) < 1 up to the transition [13]; thus, the stability bound is not saturated in this phase.As shown in Fig. 1(G), after crossing the transition, we have again S * top var(α) < 1, meaning that the richness of abundant species lies strictly below the stability bound.In Sec.V A, we will further show that the subset of abundant species lies at any time in the vicinity of a fixed point, which is thus linearly stable.This is perhaps surprising, compared to the symmetric case where the stability bound is saturated [35], and given that the dynamics are slow, as shown below.We return to this result in Sec.V A.
In contrast, S * inter var(α) and S * growth var(α) continue to grow above the bound.This is not in contradiction to the stability bound, since S * inter and S * growth count species with intermediate abundances λ ≪ N i ≪ 1, and some of them have positive growth rates, so do not satisfy the condition Ṅi / The proportion of species with intermediate abundances (λ ≪ N ≪ 1) among those with population size well above the migration floor (N ≫ λ) grows with the standard deviation of the interactions std(α), as can be seen in Fig. 1(G).

B. Dynamics and timescales
We now describe the long-time phenomenology of the dynamics, when λ > 0. The transient regime, see Fig. 2(A), is later discussed in Sec.VI together with the closely related dynamics of isolated systems, for which λ = 0.When λ > 0 and at long times, the species abundances fluctuate forever and their autocor- Crucially, the dynamics become slow when λ ≪ 1 and feature an emergent statistical invariance between realizations at different λ under the rescaling of time t → t/| ln λ|.As we show in Fig. 2(C), the autocorrelation functions for different values of λ, but identical values of σ and µ, indeed collapse to a single master curve when plotted against t/ |ln λ|, namely C λ (t, t + | ln λ|s) → Ĉ(s) when λ → 0 + .Thus, a unique timescale τ ∼ | ln λ| characterizes the autocorrelation function.We defer the proof of this result to Sec.IV where we derive the rescaled dynamics, that are also solved numerically to obtain the function Ĉ(s), see Eqs. (5,7,8).
Timescales of order | ln λ| are expected to appear, since it takes a time t ∼ | ln λ| for a population to grow from the migration floor λ to O(1) population size under exponential growth with finite growth rate.Importantly, we find that there is no shorter timescale τ ≪ | ln λ| relevant for describing the effective dynamics of a single species.Indeed, the master curve Ĉ(s) is regular at s → 0 + .Recall that λ is non-dimensional, and that | ln λ| can be quite large in ecologically relevant contexts, see Section II, in which case fluctuations become correlated over many generation times.

IV. RESCALED DYNAMICS
In Sec.III we described the phenomenology of the long-time fluctuating dynamics when 0 < λ ≪ 1.We have seen that the species abundances N i (t) fluctuate over a long timescale of order | ln λ|.Furthermore, the log-abundances ln N i dynamically explore values from O(ln λ) to O(1).These observations motivate defining rescaled variables s ≡ t/| ln λ| and z i ≡ ln (N i ) / |ln λ|.z i (s) turns out to follow a well-defined stochastic process when λ → 0 + (that no longer includes any λ dependence).It is a single-variable process describing the probability of trajectories of a single species within the large system.Subsection IV A is devoted to the mathematical derivation of this process.Subsection IV B describes its properties, and can be read independently of subsection IV A. Later sections discuss its implications to the species abundance distribution and species diversity.

A. Derivation
Our starting point is the effective single-species stochastic dynamics for N i (t), previously obtained in the limit where the number of species is very large and for any value of λ [16].For the sake of completeness, we briefly outline the steps leading to these dynamics.In Eq. ( 1), the dynamics of N i (t) are driven by the influence of all other species through g i (t) = 1 − j α ij N j (t), which involves the sum of many weakly-correlated contributions.Dynamical Mean-Field Theory (DMFT), valid when S → ∞, shows that the g i (t) are identically distributed Gaussian processes, and are independent for i ̸ = j.This implies that population sizes N i (t) behave as independent realizations of the single-variable stochastic process, where the subscript i has been dropped.The first two moments of the Gaussian process g(t) obey self-consistent closure relations, that relate the input noise g(t) to the output N (t), Here the angular brackets ⟨.⟩ denote an average over the realizations of g(t) (and the initial conditions N (0) which are irrelevant at large times).The derivation of Eq. ( 2) follows a standard procedure [36][37][38][39] which was applied to the Lotka-Volterra equations in [16].In the long-time limit where ⟨N (t)⟩ → ⟨N ⟩, the entire dynamics is controlled by the two-time correlation C(t, t ′ ) ≡ ⟨N (t)N (t ′ )⟩, since the Gaussian process g(t) is completely characterized by its correlations and mean, given in (3,4).
In order to study the behavior of these equations when λ ≪ 1, we introduce z ≡ ln(N )/| ln λ| and s ≡ t/| ln λ| so that Eq. ( 2) becomes The non-linear terms become impenetrable boundaries when λ → 0 + , since Hence the process z(s) is confined between −1 and 0 when λ → 0 + .The confinement originates from the migration term and the self-regulation term, proportional to N 2 i , in Eq. ( 1).The effective noise g(s) is not able to push z(s) to outside of the confining region, because its mean and variance are finite provided population sizes do not blow up, see Eqs. (3,4).Thus z(s) obeys in the low migration limit, where W (z) and W (z + 1) account for the confining boundaries at z = 0 and z = −1.The autocorrelation function of g(s) is proportional to the master function Ĉ(s) introduced in Sec.III, since ⟨g(s)g(s We now derive the evolution of the abundance N (s).Beyond the fact that N (s) is the main quantity of interest, such an evolution is also necessary to derive the self-consistent equations (3,4) in the the λ → 0 + limit, and obtain a closure of the DMFT equations in terms of the process z(s).When z(s) < 0 it is clear from the definition N (s) = exp (|ln λ| z(s)) that N (s) = 0 in that limit.However, this relation appears ambiguous in the double limit z → 0 and λ → 0 + .To remove the ambiguity we use the impenetrability condition at the z = 0 boundary, namely W (z(s)) = g(s) when z(s) = 0. Since by definition of z and Eq. ( 5), W (z(s)) = N (s), we obtain the relation where Θ(z) is the Heaviside function with the convention Θ(0) = 1.The existence of this impenetrability condition rests on the facts that (i) g(s) does not have a white noise component, which follows from the fact that its mean and variance remain finite when λ → 0 + and (ii) that Ĉ(s) has a well-defined limit when λ → 0 + , which agrees with numerical simulations (see Fig. 2

B. Properties of the limit dynamics
Equation ( 5) describes the effective evolution of z i (s) for any single species within the many-variable system.It is driven by a Gaussian noise g i (s) and confined to values −1 ≤ z i ≤ 0. The mean and variance of the Gaussian noise g i (s) are determined self-consistently through Eqs.(3,4), which become Eqs.(7,8) in the limit λ → 0 + .The noise g i (s) can be interpreted within the original many-species species dynamics Eq. ( 1) as the effective growth rate set by all the other species, g i (s) = 1 − j α ij N j (s), which can indeed be shown to be Gaussian distributed in the limit of many species.When −1 < z i (s) < 0, meaning λ ≪ N i ≪ 1, the dynamics read z ′ i (s) = Ṅi (t)/N i (t) = g i (s), resulting in exponential growth or decay of the population sizes under the effective growth rate set by all the other species.
In the limit λ → 0 + , the dynamics of the population size N i (s) is related to that of z i (s) by Eq. ( 6).When z i (s) < 0, Eq. ( 6) yields N i (s) = 0 which naturally follows from taking the limit λ → 0 + in the definition N i = exp(| ln λ|z i ).Yet Eq. ( 6) goes further, relating N i and z i in the double limit where both λ → 0 + and z i = 0.The interpretation of Eq. ( 6) is that the species with z i = 0, meaning with a finite population size N i , are in a slowly-changing fixed point.Indeed, Eq. ( 6) can also be obtained from the many-body dynamics using the slowness of the dynamics discussed in Sec.III when λ → 0 + : by Eq. ( 1), slow changes Ṅi (t) ≃ 0 require Equations (3,4) guarantee that g(s) has finite variance and mean, so that g(s) can not contain a white noise component.Also, by solving these rescaled DMFT equations following the procedure detailed in Appendix A 6, we find that Ĉ′ (0 + ) is finite, thus confirming the absence of fast-time scale in the Lotka-Volterra dynamics, see Fig. 8.We further show that Ĉ′ (0 + ) ̸ = 0, meaning that g(s) is rough, namely nowhere differentiable (like in Brownian motion [40]).The fact that Ĉ′ (0 + ) ̸ = 0 follows from Eq. ( 10) below, see the discussion there.An example run of the rescaled dynamics is shown in Fig. 3.As can be seen, −1 ≤ z ≤ 0, and z spends finite time intervals at the boundaries, which is a consequence of the finite memory of the noise g(s).Indeed, if g > 0 and z = 0 at a given time, it will remain so for a finite amount of time.
Lastly, we discuss important features of the dynamics of the population sizes N (s) seen in Fig. 3. First, N (s) is rough since N = g when z = 0, while z(s) is more smooth since it is an time integral of g(s).Second, the limiting process N (s) is not continuous in time and features jumps from 0 to a finite value, after which N (s) continuously reaches 0. The jumps represent species that grow from very small abundances at a finite growth rate, and thus their time to reach N = g from a small fixed value (that doesn't depend on λ, say N = 10 −5 ) is finite in time t, and so vanishes in the rescaled time s.These jumps are precisely species turnover events, that drive the change in the composition of the abundant species.

A. Diversities revisited
In Sec.III A, we proposed three definitions for diversity S * top , S * inter , S * growth .Denote ϕ top ≡ S * top /S and similarly for ϕ inter , ϕ growth .All these quantities are welldefined in the limit λ → 0 + .Indeed, they take simple forms in terms of the limiting process described in Sec.IV: We now return to the discussion of ϕ top and its relation to the stability bound, see Sec.III A. As followed from the rescaled dynamics, the abundant species counted in ϕ top are approximately at a fixed point, while all other species have negligible abundances, and so do not affect this fixed point.Thus, one indeed expects the bound S * top var(α) = σ 2 ϕ top ≤ 1 to hold, as is clear in Fig. 1(G).Indeed, a fixed point with S * top coexisting species and S * top var(α) > 1 would be typically linearly unstable to perturbations [11,12,35].
A natural question is: Is the stability bound saturated, i.e. σ 2 ϕ top = 1, resulting in fixed points of abundant species that are near marginal stability?One could perhaps expect marginal stability, as the rate at which lowabundance species are added to the subset of abundant ones is slow, "gently" perturbing the fixed points, and so would perhaps allow S * top to increase up to the stability bound.In addition, marginality is reached in Lotka-Volterra dynamics [35], Eq. ( 1), with symmetric interactions (α ij = α ji ).And more generally, slow dynamics are in many cases associated with marginality (see, e.g., [41]).
Yet, as was shown in Sec.III A, see Fig. 1(G), σ 2 ϕ top < 1 so the stability bound is not saturated.Con-sequently, as λ → 0 + , the high-abundance species lie at any time in the vicinity of a fixed point that is linearly stable to perturbations applied to those species.Note that this fixed point changes over time, as it is destabilized by the growth of species from rare.
To obtain the bound σ 2 ϕ top < 1, we prove that while species turnover is a slow process, the jumps when going from rare to abundant, seen clearly in Fig. 3, are sufficient to significantly perturb the subset of abundant species and prevent it from reaching marginal stability.We prove that by deriving an exact relation that links diversity with temporal fluctuations, defined as follows.
Let N jump be the size of these jumps.Let G be the rate of incoming species weighted by N 2 jump : that is, the sum of N 2 jump over all jumps taking place in a unit of rescaled time s in the many-species dynamics, and divided by the number of species S. The relation, derived in Appendix A 3, reads: Here Ĉ (s) is the autocorrelation of N (s) in the rescaled time defined in Sec.III.Ĉ′ (0 + ) is finite, see Fig. 2(C) and Fig. 8 in Appendix A 1. Eq. ( 10) then limits ϕ top to be below the stability bound; indeed, marginal stability 1 = σ 2 ϕ top is only possible if G = 0, namely species do not perform jumps, in contradiction with the dynamics in Sec.IV.Additionally, we note that G > 0 implies Ĉ′ (0 + ) > 0, thus proving that the trajectories N i (s) are rough when N i (s) ̸ = 0. Put differently, the introduction of one new species leads to the removal of others, with the average number of removed species growing as one approaches the marginal diversity.The balance, in which one species is removed for each one introduced, sets S * top , that only reaches some fraction of the bound 1/var(α).In dynamics that reach a fixed point, the requirement that all species involved have positive abundance is known as feasibility [42], and it is what limits the diversity in the fixed point phase, σ < σ c in Fig. 2(F).Eq. 10 can be thought of as an extension of the requirement to dynamics, a form of "dynamical feasibility".

B. Species Abundance Distribution revisited
We now return to the species abundance distribution P (N ).As mentioned in Sec.III A, P (N ) behaves roughly as 1/N in the intermediate range λ ≪ N ≪ 1.Using the rescaled dynamics, we can refine this statement.The dynamics of z i (s) (Sec.IV) spend a finite fraction of the time at the boundaries z = −1, 0. This translates to two delta-peak contributions in P (z) at these values.In addition, there is a regular contribution for −1 < z < 0. Together, this reads  where h(z) is a smooth function with 0 −1 h(z) dz = 1.ϕ top , ϕ inter were defined in Eq. ( 9).Fig. 4(B) shows the collapse of P (z = ln N/| ln λ|) as λ → 0 + to this form.
Eq. ( 11) sets the form of the abundance distribution in the intermediate range λ ≪ N ≪ 1. Changing variables from z to N , we get This refines the 1/N dependence with an additional, slowly-varying correction.As we discuss in Sec.VIII, this correction can appear to change the power law exponent of the species abundance distribution when λ is finite, and only parts of the entire distribution are sampled.The distribution of top species can be inferred from P [g|z], the distribution of the growth rate g conditioned on the value of the rescaled abundance z.For N ≥ 0, we get Fig. 4(A) shows the convergence of the distribution P (N ) to this limiting distribution as λ → 0 + .The results in the limit λ → 0 + were obtained by solving numerically the rescaled DMFT equations, see Eqs. (5,7,8).The limiting distribution P (N ) deviates from the truncated Gaussian SAD obtained in the fixed point phase.

A. Dynamics of an isolated system
Isolated systems are characterized by zero migration rate, λ = 0, which is a singular limit of the Lotka-Volterra system of equations with a large number of species, in the chaotic phase.Indeed, the timescale | ln λ| which characterizes these dynamics at finite λ, diverges when λ → 0. For λ = 0, the dynamics are not time-translation invariant but forever slow down in time, as evidenced in a linear growth of the correlation time as a function of the elapsed time, a behavior known in physics as 'aging'.Formally, C λ=0 (t, t + tt ′ ) → Ĉ(t ′ ) when t is large, see Fig. 5(A) for the collapse of results from numerical simulations.Again, the master curve Ĉ(t ′ ), which depends on the parameters σ and µ, is regular at t ′ → 0 + .This behavior can be understood as follows: When λ = 0 the lowest values of ln N reached after time t are of order ln N ∼ −t.If a species changes to positive growth rate at this time, it will therefore take another time t for its population size to be O (1).This sets the correlation time.A similar mechanism was found in another, exactly-solvable, model [21].
The proof of the existence of this aging regime is given in App.A 2. Using a reasoning similar to that employed in Sec.IV, we show that the transformed variables s ≡ ln t, z ≡ ln(N )/t obey a well-defined set of DMFT equations which become time-translation invariant when t → ∞.These transformations reflect a growth of both timescales and log-fluctuations with the elapsed time.The resulting process z(s) is different from that of Sec.IV, and is described in detail in Appendix A 2.
An important consequence of this result is that the collective deterministic dynamics in an isolated system drives the population size of any species arbitrarily close to 0 as time grows.Considering that actual populations sizes are finite, integer numbers, this process will inevitably lead to the extinction of many species, subsequently leading to an arrest of the fluctuations, as suggested in [17,18].

B. Transient dynamics at finite λ
When migration is present, this dynamical slowdown provides a mechanism by which the correlation time grows until time t transient ∼ | ln λ|, where the correlation time reaches the value | ln λ| discussed above.This is indeed the time it takes for a finite fraction of the species in the community to reach the migration floor, when starting with all species with population sizes of O(1), and before which λ can be safely set to zero.
We verify in simulations that the transient dynamics is characterized by a linear growth of timescale with the elapsed time, which is interrupted at a time t transient ∼ | ln λ|.For that, we measure the time constant τ 2 (λ, t) it takes for the autocorrelation function C λ (t + τ, t) − C λ (∞, t) to reach a fraction e −1 of its τ = 0 value, starting at initial time t and with migration rate λ.Fig. 5(C) shows that the growth and saturation of the timescale follow the scaling relation where f is a smooth function with f (x) ∝ x at small x, encoding the slowdown of the dynamics in the transient regime, and f (x) approaching a constant as x → ∞, encoding the time-translation invariant behavior of the long-time dynamics, with correlation timescale | ln λ|.

VII. ROBUSTNESS OF THE PREDICTIONS
The theory is built around two limits, S → ∞ and λ → 0 + , and assumes that the interaction coefficients α ij are sampled independently.In this last part of this work, we assess the robustness of our predictions when these assumptions are relaxed.We find that the main qualitative features discussed in this work are robust against changes in model definition, and relevant even at reasonable values of migration rate and number of species.

A. Finite migration rate λ and number of species S
We start by discussing how the key quantities, ⟨N ⟩, N 2 and ϕ top vary with S, λ, when measured from abundance data gathered in numerical simulations of the original many-species dynamics Eq. ( 1).Regarding the moments ⟨N ⟩ , N 2 , we find that the dependence on S and λ is very weak for ⟨N ⟩ and weak for N 2 , respectively within roughly 1% and 10% in the inspected range of parameter, see Fig. 11 in Appendix A 5.
To measure ϕ top at finite λ, we test three options.The first measure, ϕ N top (S, λ, ϵ), is simply a threshold on the abundance: counting all the species with N i ≥ ϵ at a given time, for some chosen ϵ.Unsurprisingly, ϕ N top (S, λ, ϵ) is more sensitive to λ than ⟨N ⟩ or N 2 .For moderate values of S and λ, we find that the asymptotic value ϕ top nevertheless provides a reasonable estimate of ϕ N top (S, λ, ϵ), see Fig. 6(A, B), though the discrepancy increases with the strength of the interactions, see Fig. 12(B).This last point suggests that disentangling species at high and intermediate abundances becomes harder when the scale of the interactions, and thus the amplitude of the fluctuations, increases.The second measure, denoted by ϕ g top (S, λ, ϵ), refines the first one and requires N i ≥ ϵ but also g i > 0. We find that both measures converge to the asymptotic ϕ top as 1/S in S. Yet the convergence with λ is significantly different, with ϕ N top converging as |ln λ| −1/2 , and ϕ g top faster, as |ln λ| −1 , see Fig. 6(B, D).This highlights the relevance of the growth rate g in determining "top species", namely members of fixed points.We find that the measure ϕ N top (S, λ, ϵ) is often above the stability bound, in contrast to the asymptotic ϕ top which is always below it.Due to the faster convergence, ϕ g top (S, λ, ϵ) can be either above or below this bound, depending on the parameters, see Fig. 6(C).The origin of these convergence rates is discussed in Appendix A 4.
Finally, we consider a situation where one is able to manipulate the system, by removing certain species and continuing the dynamics.Namely, after a long time so that transients have passed, we kill all species with N i /g i < 1/2.These correspond to species with positive growth rate that are in the midst of their jump, below halfway, as well as those with negative growth rate.(Recall that, asymptotically as λ → 0 + , abundant species are those for which that N i /g i = 1.)Then, we run again the dynamics, and we find that the remaining species reach a stable equilibrium.The properties of the equilibrium obtained in this way are strikingly similar to those predicted by the asymptotic theory, even for reasonable λ.The asymptotic diversity and the asymptotic distribution in Eq. ( 13) can be almost exactly reproduced, see Appendix A 6.

B. Correlation between the matrix elements
In addition to taking the asymptotic limits in S, λ, the model above assumes that interactions are statistically asymmetric with vanishing correlation coefficient corr(α ij , α ji ) = 0. We show that two of our main qualitative results-how the timescale grows and that top diversity is below the stability bound-also hold when this assumption is relaxed.Allowing for more symmetric or anti-symmetric interactions, we take a correlation coefficient corr(α ij , α ji ) = γ with −1 < γ < 1.
The growth of the timescale τ ∼ τ λ ≡ |ln λ| is clearly seen in Fig. 10(A,B) in Appendix A 4, the equivalent ) is within 5% of the value predicted by the rescaled DMFT, and well below the stability bound.For all reasonable values of the migration rate, (here already when λ ≳ 10 −40 ), the measured top diversity is significantly above the stability bound.ϕ N top (S → ∞, λ, ϵ) was inferred from the finite S scaling, panel A. (C,D) Same as panels A,B, but for ϕ g top (S, λ, ϵ).The convergence to the asymptotic ϕtop is more rapid; in particular, ϕ g top can be either below or above the stability bound.
of Fig. 2(C).This is expected for the same reason as when γ = 0: the time for a population to grow from λ to N ∼ 1 scales as |ln λ|.As to the diversity ϕ top , we conjecture that the slowness of the dynamics at λ ≪ 1 still introduces a clear partition between nearly-extinct species with ln N ∼ ln λ and abundant species with N = O(1).For the abundant species, the long timescale implies that they are near a fixed point.In Appendix A 4, we generalize the relation between fluctuations and diversity from Sec. V A, Eq. ( 10), to give As for γ = 0, we expect that the growth of timescale when λ → 0 + goes hand-in-hand with jumps in the rescaled dynamics of the abundances, meaning G > 0. This implies 1 + 1 − 4γϕ top σ 2 2 − 4ϕ top σ 2 > 0, so that ϕ top lies strictly below the linear stability bound [13].This is confirmed by numerical simulations, see Fig. 10(C,D) in Appendix A 4.

VIII. DISCUSSION
We begin the discussion by summarizing key predictions presented above.We identified several signatures of the many-species Lotka-Volterra dynamics, when the abundnaces does not reach a fixed point and migration rates are small.They may serve as footprints of endogenously-driven fluctuations in experimental or natural situations.
First, we predict that endogenously-driven fluctuations in time would lead to broadly distributed abundance distributions, in contrast with the fixed point distributions.For very low migration, the abundance distribution at intermediate values (meaning small population sizes well above the migration floor) are predicted to behave as P (N ) ∼ 1/N .At any finite λ, there are slowly varying corrections to this 1/N behavior, see Eq. ( 12).As P (N ) is broad, one can always define a slowly-varying 'local' power law ν(z) by the slope of ln P (N ) versus ln N at given z = ln N/ |ln λ|.It has corrections of order , which might explain deviations from ν = −1 in observed abundance distributions [6].These corrections vary with z and are arbitrarily large over the entire range −1 < z < 0, so no unique exponent other than −1 can be defined over the entire range of N .
Second, we predict that a single timescale controls the dynamics, predicted to grow as | ln λ| when lowering the migration rate.This could be tested in controlled experiments, for example via the abundance autocorrelation.
Moving on to diversity, different definitions of diversity give different results.For example, the number of species that can grow from rare at a given time, is generally different from the number of species that have high abundance, in contrast with the situation at a fixed point.We show that the number of species with high, intermediate and low abundance is well defined, namely insensitive to the precise threshold above (or below) which the number of species is counted.We show that the number of species at high abundance is significantly below the May bound (their fraction of the total number of species is below the fraction allowed by the bound).This makes the community of high-abundance species a stable equilibrium of the dynamics if the other species are removed.The distance to the stability bound increases with the strength of the interactions σ, see Fig. 1(G), and also Fig. 12(B) in Appendix A 4. In contrast, when including the intermediate ones (in the power law regime of P (N )), the total number goes above the bound.The same is true when including all species in the pool that may invade.It is this last fact that drives species turnover: there are always species that grow from rare to replace the ones at high abundance.
Last, a key property of the dynamics at low migration is the existence of jumps from rare in the dynamics of population sizes, see Fig. 3.For finite λ, this manifests itself in a strong asymmetry of "blooms", namely trajectories where the abundance of a species increases from rare before returning there.This can be clearly observed in time series at finite λ, see Fig. 7, and would be very interesting to look for in experimentally-measured time series.
In conclusion, the dynamically-fluctuating phase of high-diversity ecological communities is a promising direction to explain key features of natural high-diversity ecosystems.We offered a list of additional predictions expected in this phase.It would be interesting to further investigate the robustness of these features upon modifying the structure of the interaction matrix.Along this line, we note that the chaotic dynamics appear at finite S also without any beneficial interactions (α ij < 0) between species, and the analysis above is expected to hold.We further note that the limit of large σ, µ in our framework connects to other asymptotic limits [18,43].A recent work on the strongly interacting case [44] suggests from numerics that the qualitative picture of the present study may extend to that regime, in particular the existence of a growing timescale, and the fact that dynamics evolve in the vicinity of fixed points.A different and very interesting direction for future research is understanding how these results extent to spatially-extended metacommunities [17,18], beyond a constant migration from an unspecified species pool, as assumed here.This question is pertinent, given that the present work shows that chaotic fluctuations cannot generically be sustained in isolated high-diversity systems, due to extinctions.Trajectories of N (t) while at high-abundance display a clear asymmetry in time, with a rapid initial increase and more gradual decrease.This is conspicuous even for migration rates that are not very small (here λ = 10 −6 ).It is the finite-λ counterpart of the asymptotic behavior at λ → 0 + , featuring sharp jumps from zero to positive N , see Fig. 3. Shown are trajectories that go above N = 10 −2 at some time tin and stay above it until time tin + ∆t, and reach N ≥ 0.1 at some intermediate time.Here ∆t ≃ 5.8 |ln λ| that is the most common length of such trajectories.In light grey are 22 example trajectories, and the thick blue line shows the average over many such trajectories.Note however that in the λ = 0 case, the process z(s) is confined in the negative half-line by a harmonic potential and not by a hard boundary, see Eq. (A2).
Under the condition that g(s) has time-translation invariant correlations, Eq. (A2) manifestly predicts that z(s) reach at long time a time-translation invariant state.The closure equations Eqs.(A4)-(A5) then self-consistently show the validity of the time-translation invariant ansatz in rescaled time s, eventually showing the validity of the aging scaling given in Eq. (A1).In Fig. 9, we show the large-time convergence of the distributions P (N ) and P (z) to those predicted by the λ = 0 rescaled dynamics presented here.

Diversity limited by fluctuations
The two cases λ → 0 + and λ = 0 share a crucial property: the long-time dynamics are very slow, reflecting the fact that the system evolves in the vicinity of feasible fixed points (all population sizes are ≥ 0).More precisely, at any given time s, the system is close to a fixed point comprised of some abundant species with O(1) population sizes (corresponding to the fraction ϕ top of species with z(s) = 0), and some rare species (corresponding to the fraction 1 − ϕ top of species with z(s) < 0).Species turnover happens because these fixed points are invadable, meaning that some nearly-extinct species have positive growth rate.Furthermore, we see that the subset of abundant species, when taken alone, is linearly stable and not marginal.Indeed, the fraction ϕ top of top species does not saturate the stability bound with ϕ top σ 2 < 1, see Fig. 1(G).Despite the fact that abundant species are found in the vicinity of a stable fixed point, the dynamics exhibit abundance fluctuations due to the continuous flux of incoming species from the pool of nearly extinct ones.
Here we derive a relation between temporal fluctuations and the observed diversity valid both in the aging (λ = 0) and chaotic (for λ → 0 + ) regimes of the many-body Lotka-Volterra system of equations that then establishes the linear stability of the subset of abundant species ϕ top σ 2 < 1, see Eq. (10).We take advantage of the slowness of the dynamics to generalize the calculation of the fluctuations induced by a random perturbation to a fixed point [13].Between the times s and s + ds, new species become abundant and induce a perturbation on the species already abundant at time s.We then relate the fluctuations of their population sizes to the amplitude of the effective perturbing field and conclude by using the DMFT closure relations.Henceforth we use the notations Θ = Θ(z(s)) and for any quantity x(s) we write x(s) = x and δx = x(s + ds) − x(s).Using N (s) = g(s)Θ(z(s)), we obtain δN = (Θ + δΘ)δg + gδΘ .(A6) We assume that over short times intervals ds, the changes in g(s) scale to leading order as δg ∼ √ ds, as for Brownian motion.Furthermore, δΘ ∈ {−1, 0, 1} and its moments scale as O(ds).To prove the latter, we evaluate the mean number of species with δΘ = −1, namely the mean number of species going from z < 0 to z = 0 between s and s + ds (which is also the mean number of species going from z = 0 to z < 0 in the same time interval).Denoting where we used the above result to obtain the last equality.We note that a generic time-translation invariant Gaussian process can be generated from the Langevin equation, with η(s) a Gaussian white noise and J(s) a suitably chosen memory kernel that enforces ⟨ξ(s)ξ(s ′ )⟩ = C(s − s ′ ), i.e. in Fourier space (with J(s < 0) = 0) This means that the O( √ ds) increments of δg are statistically independent from the previous history of the system.We can therefore write Hence we obtain, using lim ds→0 In the right-hand side we recover the quantity G introduced in Sec.V A of the main text and the above equation reduces to Eq. (10).Note that we can express G as In terms of the many-body dynamics, the above equation becomes

Generalization for γ ̸ = 0
We now consider the case γ ̸ = 0 in the limit λ → 0 + .We generalize the previous argument connecting fluctuations, diversity and species turnover based on the many-body equations of motion and derive Eq. ( 14).We assume that the previous scaling for the timescale of the correlation matrix holds, which we check numerically, see Fig. 10.In other words, at any time s = t/| ln λ|, the system is close to a fixed point, meaning that some species (a fraction ϕ top of Lastly, for the surviving species, we have Based on the DMFT analysis of the γ = 0 case, we assume that the O( √ ds) perturbing fields h j are statistically independent of the state of the system at time s and that to leading order in ds, h j , h k are uncorrelated for j ̸ = k.Therefore we obtain, We furthermore assume that to compute the trace we can take α * to be sampled from the same ensemble as the full interaction matrix α (albeit with a smaller size), neglecting the correlations induced by the dynamical selection of the community of surviving species.Relying on the cavity method, the validity of this approximation has been argued in the fixed point phase [13,45] and more generally at the level of the average spectral density, see [35] (the possible existence of two outlying eigenvalues [45] is sub-extensive in the trace).We thus have, Together, we obtain Hence, we recover Eq. ( 14) of the main text and where I(s, s + ds) is the subset of species experiencing a jump from rare during the interval [s, s + ds].For γ = 0, Eq. (A9) reduces to Eq. (A7) derived within the DMFT formalism.Additional figures and details for the Disscussion section

Convergence with S and λ
Fig. 11, presents the convergence of ⟨N ⟩ , N 2 , showing that it is quantitatively quite robust to changes in λ, S. One implication is that other definitions of diversity besides than species richness, will be quite robust.For example, defining diversity using the inverse Simpson index is relatively robust in S, λ, due to the robustness of ⟨N ⟩ , N 2 .As briefly discussed in Sec.VIII, to measure ϕ top at finite S and λ, we test two options.The first is simply a threshold on the abundance: counting all the species with N i ≥ ϵ at a given time, for some chosen ϵ ≪ 1.Let ϕ N top (S, λ, ϵ) be the diversity measured this way.The second measure, denoted by ϕ g top (S, λ, ϵ), utilizes the invasion growth rate g i .In it, in addition to N i ≥ ϵ we also requires g i > 0. For both measures, we find a similar form for the dependence on λ, S, ϵ.For ϕ N top , where A = P (N = 0 + ), c 1,2 are constants, and the exponent β N = 1/2.The expression for ϕ g top is of the same form, only with β g = 1.This form of convergence can be understood as follows.We find that the asymptotic distribution P (z) diverges as (−z) −1/2 when z → 0 − .This can be traced back to the fact that ż = 0 when z(s) leaves the confining wall.Therefore, by defining a cutoff N ≥ ϵ, or equivalently z ≥ ln ϵ/| ln λ|, the error made in sampling the distribution P (z) for z < 0 scales as | ln ϵ/ ln λ|, hence the result in Eq. (A10).However, when conditioned on g > 0, the distribution P(g > 0, z) is found to be finite when z → 0 − stemming from the fact that ż > 0 when z(s) reaches the wall from below.Thus, by defining a cutoff N ≥ ϵ, or equivalently z ≥ ln ϵ/| ln λ|, the error made in sampling the distribution P(g > 0, z) for z < 0 only scales as | ln ϵ/ ln λ|, hence the exponent β g = 1.These properties of the steady-state distribution P(g, z) seem to be robust features of persistent random walkers confined by hard obstacles and have been discussed in other contexts [46][47][48].Since in practice |ln λ| would not be a very large number, the difference in the convergence in λ is important.Quantitatively, ϕ N top , ϕ g top are typically larger than ϕ top , and the convergence in λ of ϕ g top is indeed much faster, again highlighting the relevance of g.We find that the measure ϕ N top (S, λ, ϵ) is often above the stability bound, in contrast to the asymptotic ϕ top which is always below it.ϕ g top (S, λ, ϵ) can be either above or below this bound, depending on the parameters, see Fig. 6.

Identifying top species by selection and subsequent dynamics
The separation of top species from the rest is guaranteed for λ → 0 + .Yet as discussed in Sec.VII, for finite S, λ this separation is not very clear when looking at the abundances, see P (N ) in Fig. 4(A).By using insights gleaned from the theory, we have shown, see Sec.VII and Fig. 6, that this separation can be much better defined even away , and removing species that are about to grow and disrupt the system, here the orange trajectory, and those with a negative invasion growth rate (by removing species i if Ni/gi < 1/2), and then continuing the dynamics (dashed lines).The outcome is an equilibrium.(B) The distribution of abundances thus obtained is close to the asymptotic distribution obtained from the rescaled dynamics.It is much closer than the bare abundance distribution P (N ), and is also very different from the distribution in the equilibrium phase.S = 5000, σ = 1.8, µ = 10, λ = 10 −10 .
from the asymptotic limit, by also considering g i , the invasion growth rate defined in Sec.III A.Here we show, that if we can manipulate the system, by removing certain judiciously-chosen species and rerunning the dynamics, the asymptotic distribution P (N ) can be remarkably well reproduced at reasonable S, λ.Beyond potential applications to experiments, this seems to suggest that the separation between abundant species and the others is still present, even when it seems blurred with other measures.The lack of clear separation is most pronounced at not-very-small N (say N ∼ 0.1).There, one finds: (1) abundant species, (2) species that grow from rare and are about to disrupt the abundant species, and (3) species with negative growth rate that are leaving the abundant subset.The abundant species (group 1) are characterized by appreciable N , and being at a fixed point, g ≃ N .They occupy an equilibrium that is fully stable if not disrupted by species that grow from rare (group 2).These can have significant g, with still N small compared to g.The idea here is to remove species (2) and ( 3), which we do by removing all species that have N i /g i < 1/2, corresponding to species with positive growth rate that are in the midst of their jump, below halfway, and those with negative growth rate in Fig. 3,7.Then, we run again the dynamics, and we find that the remaining species reach a stable equilibrium.The properties of the equilibrium obtained in this way are remarkably similar to those predicted by the asymptotic theory, even for reasonable λ.Fig. 13 shows the obtained abundance distribution P (N ).It is similar to the asymptotic distribution, and in contrast quite different from the instantaneous unfiltered P (N ) from Fig. 4(A).Furthermore, the resulting distribution is distinct from those in the fixed point phase (when σ < √ 2), as shown by a comparison with the predicted distribution there, a truncated Gaussian.The asymptotic species richness ϕ top = 0.28 is recovered within 1%; this is lower than the stability bound, at ϕ = 0.31.
Numerical methods Here we detail the numerical procedures used in producing the figures.They are of two types: (1) Full many-variable simulations of Eq. ( 1).

Figure 1 .
Figure 1.Species richness and species abundance distributions at fixed points (left column) and a persistently fluctuating state (right column).(A) Dynamics reaching a fixed point.(B) Persistent dynamics, where the abundances fluctuate indefinitely in the range λ ≲ Ni ≲ 1.(C-F)The corresponding distributions for N and ln N , denoted by P (N ) and P (ln N ) respectively.The distributions contain three parts (see E,F): a top part, which remains O λ 0 when λ → 0 + , containing S * top species; an intermediate part, at λ ≪ N ≪ 1; and a low part at O(λ).The intermediate part contains a finite fraction of the species in the dynamical phase, but not in the fixed point phase.(G) The three definitions of species richness, normalized by var(α) so that 1 is the fixed-point stability bound, as a function of the interaction variability σ.For σ < σc a fixed point is reached, all three definitions coincide and lie below the stability bound, and agree with known theory (solid line).At σ = σc the stability bound is reached.Beyond it, there are persistent fluctuations, and the three definitions no longer coincide: S * top is lower than the stability bound, while S * inter and S * growth are above it.The diversities are obtained by solving the rescaled dynamics defined in Sec.IV, and agree quantitatively with careful analysis of the full equations of motion (see Sec. VII).Simulation parameters for all figures are given in Appendix A 6.

Figure 2 .
Figure 2. Species dynamics and correlations.(A)Species are initialized with order-one values.We follow ln Ni(t).In the transient regime, the changes in ln Ni(t) are comprised of three elements: downward motion, roughly in straight line (maintaining its slope for a large part of the decay between largest and smallest values, corresponding to exponential decay of Ni); upward motion, roughly in a straight line; and slow changes at high population values, Ni(t) ≃ O(1).The dynamics slow down, with the time spent in each of these elements increasing with the time since the start of the dynamics, see Sec.VI.As the excursions grow longer, the average value of ln Ni(t) decreases linearly in t.This transient regime ends when a finite fraction of species reach the migration floor Ni(t) ≃ O(λ), after a time of order |ln λ|.At long times, after the transient regime is over, ln Ni(t) performs the three dynamical elements described above, and a forth one, slowly changing around the migration floor.The vertical red dashed line marks the crossover between these two regimes.(B) At long times, timescales do not grow in time anymore, as seen in the correlation function that depends only on time differences, C(t, t ′ ) = C(t − t ′ ).(C, D) This timescale is proportional to |ln λ|.This is demonstrated by the data collapse when plotting C against t/ |ln λ| presented in panel(C), compared with the data without rescaling shown in (D).

Figure 3 .
Figure 3.The rescaled dynamics.Trajectories of N and z ≡ ln(N )/ |ln λ|, as a function of the rescaled time s ≡ t/ |ln λ|, in an example run of the limiting rescaled dynamics when λ → 0 + .

Figure 4 .
Figure 4. Collapse of species abundance distributions.Numerically measured distributions of N (A) and z = ln N/| ln λ| (B) converge to the distributions predicted by the rescaled process as λ → 0 + . 0

Figure 5 .
Figure 5. Growth of timescales in dynamics without migration.(A,B) When λ = 0, the collapse of C(t, t + τ ) as a function of τ /t demonstrates the linear growth of timescales with the elapsed time.Data without rescaling, as a function only of τ , is shown for comparison in (B).(C,D) Crossover at finite λ from a transient regime exhibiting the λ = 0 phenomenology to the long-time behavior, as identified by the time constant τ2 of the correlation function.After an initial transient of finite time duration, τ2 grows linearly with the elapsed time and independently of λ, a trademark of the λ = 0 dynamics.The crossover to the long-time behavior happens around times proportional to | ln λ|, after which τ2 stabilizes to a value proportional to |ln λ|.This scaling of the crossover is manifest in the data collapse in (C), with the data without rescaling displayed in (D) for comparison.τ2(t, λ) is defined as the time when [C(t, t + τ2) − C(t, ∞)] / [C(t, t) − C(t, ∞)] = e −1 .The value of C(t, ∞) is estimated through an exponential fit of the function C(t, t + τ ) as a function of τ .

Figure 6 .
Figure 6.Average top diversity measures, ϕ N top (S, λ, ϵ), ϕ g top (S, λ, ϵ), as a function of S and λ. ϕ N top is the fraction of species with Ni > ϵ = 10 −3 ; ϕ g top also requires gi > 0. (A) ϕ N top (S, λ, ϵ) as a function of S converges as ϕ N top (S → ∞, λ, ϵ) + #/S, at fixed λ.Data points at S = 500 are excluded from the linear fit.(B) Large S value of the top diversity, ϕN top (S → ∞, λ, ϵ), converges as |ln λ| −1/2 to its asymptotic value ϕ N top (S → ∞, λ → 0 + , ϵ).The asymptotic value ϕ N top (S → ∞, λ → 0 + , ϵ) is within 5% of the value predicted by the rescaled DMFT, and well below the stability bound.For all reasonable values of the migration rate, (here already when λ ≳ 10 −40 ), the measured top diversity is significantly above the stability bound.ϕ N top (S → ∞, λ, ϵ) was inferred from the finite S scaling, panel A. (C,D) Same as panels A,B, but for ϕ g top (S, λ, ϵ).The convergence to the asymptotic ϕtop is more rapid; in particular, ϕ g top can be either below or above the stability bound.

Figure 7 .
Figure 7. Asymmetry of blooms under time-reversal.Trajectories of N (t) while at high-abundance display a clear asymmetry in time, with a rapid initial increase and more gradual decrease.This is conspicuous even for migration rates that are not very small (here λ = 10 −6 ).It is the finite-λ counterpart of the asymptotic behavior at λ → 0 + , featuring sharp jumps from zero to positive N , see Fig.3.Shown are trajectories that go above N = 10 −2 at some time tin and stay above it until time tin + ∆t, and reach N ≥ 0.1 at some intermediate time.Here ∆t ≃ 5.8 |ln λ| that is the most common length of such trajectories.In light grey are 22 example trajectories, and the thick blue line shows the average over many such trajectories.

Figure 9 .
Figure 9. Collapse of species abundance distributions.Numerically measured distributions of N (A) and z = ln N/t (B) converge to the distributions predicted by the rescaled process as t → ∞.

ds→0 1
Sds i∈I(s,s+ds) N i (s + ds) 2 , where I(s, s + ds) is the subset of species experiencing a jump from rare during the interval [s, s + ds].Following Eq. (A8), the coefficient G can also expressed in terms of the steady-state distribution P[g, z] G = +∞ 0 dg P[g, z = 0 − ]g 3 .

Figure 11 .Figure 12 .
Figure 11.Behavior of the first and second moments of the population size.⟨N ⟩ and N 2 are robust predictions that weakly depend on the number of species S and the migration rate λ.In the inspected range of parameters, the variations in ⟨N ⟩ are of the order of 1% and those in N 2 of the order of 10%.Parameters: σ = 1.8, µ = 10, average over 80 realizations.

Figure 13 .
Figure 13.Identifying the top species by combined selection and dynamics.(A) Top species are identified by stopping the dynamics (solid lines), and removing species that are about to grow and disrupt the system, here the orange trajectory, and those with a negative invasion growth rate (by removing species i if Ni/gi < 1/2), and then continuing the dynamics (dashed lines).The outcome is an equilibrium.(B) The distribution of abundances thus obtained is close to the asymptotic distribution obtained from the rescaled dynamics.It is much closer than the bare abundance distribution P (N ), and is also very different from the distribution in the equilibrium phase.S = 5000, σ = 1.8, µ = 10, λ = 10 −10 .