Nonlinear Analogue of the May-Wigner Instability Transition

We study a system of $N\gg 1$ degrees of freedom coupled via a smooth homogeneous Gaussian vector field with both gradient and divergence-free components. In the absence of coupling, the system is exponentially relaxing to an equilibrium with rate $\mu$. We show that, while increasing the ratio of the coupling strength to the relaxation rate, the system experiences an abrupt transition from a topologically trivial phase portrait with a single equilibrium into a topologically non-trivial regime characterised by an exponential number of equilibria, the vast majority of which are expected to be unstable. It is suggested that this picture provides a global view on the nature of the May-Wigner instability transition originally discovered by local linear stability analysis.

We study a system of N 1 degrees of freedom coupled via a smooth homogeneous Gaussian vector field with both gradient and divergence-free components. In the absence of coupling, the system is exponentially relaxing to an equilibrium with rate µ. We show that, while increasing the ratio of the coupling strength to the relaxation rate, the system experiences an abrupt transition from a topologically trivial phase portrait with a single equilibrium into a topologically non-trivial regime characterised by an exponential number of equilibria, the vast majority of which are expected to be unstable. It is suggested that this picture provides a global view on the nature of the May-Wigner instability transition originally discovered by local linear stability analysis.
complex systems | equilibrium | model ecosystems | random matrices | W ill diversity make a food chain more or less stable?
The prevailing view in the mid-twentieth century was that diverse ecosystems have greater resilience to recover from events displacing the system from equilibrium and hence are more stable. This 'ecological intuition' was challenged by Robert May in 1972 [1]. At that time, computer simulations suggested that large complex systems assembled at random might become unstable as the system complexity increases [2]. May's 1972 paper complemented that work with an analytic investigation of the neighbourhood stability of a model ecosystem whereby N species at equilibrium are subject to random interactions.
The time evolution of large complex systems, of which model ecosystems is one example, is often described within the general mathematical framework of coupled first-order nonlinear ordinary differential equations (ODEs). In the context of generic systems, the Hartman-Grobner theorem then asserts that the neighbourhood stability of a typical equilibrium can be studied by replacing the non-linear interaction functions near the equilibrium with their linear approximations. It is along these lines that May suggested to look at the linear model J jk y k , j = 1, . . . , N , [ 1 ] to study the stability of large complex systems. Here J = (J jk ) is the coupling matrix and µ > 0. In the absence of interactions, i.e., when all J jk = 0, system (1) is self-regulating: if disturbed from the equilibrium y1 = y2 = . . . = yN = 0 it returns back with some characteristic relaxation time set by µ. In an ecological context yj(t) is interpreted as the variation about the equilibrium value, yj = 0, in the population density of species j at time t. The element J jk of the coupling matrix J, which is known as the community matrix in ecology, measures the per capita effect of species k on species j at the presumed equilibrium. Generically, the community matrix is asymmetric, J jk = J kj . For complex multi-species systems information about the interaction between species is rarely available at the level of detail sufficient for the exact computation of the community matrix and a subsequent stability analysis. Instead, May considered an ensemble of community matrices J assembled at random, whereby the matrix elements J jk are sampled from a probability distribution with zero mean and a prescribed variance α 2 . This is similar to the approach taken by Wigner in his description of statistics of energy levels of heavy nuclei via eigenvalues of random matrices, which proved to be very fruitful [3]. A detailed review of May's model in the light of recent advances in random matrix theory can be found in [4].
The linear system (1) is stable if, and only if, all the eigenvalues of J have real parts less than µ. Invoking Wigner's arguments for studying eigenvalues of large random matrices, May claimed that for large N the largest real part of the eigenvalues of J is typically α √ N . Obviously, the model's stability is then controlled by the ratio m = µ/(α √ N ). For N large, system (1) will almost certainly be stable if m > 1 and unstable if m < 1, with a sharp transition between the two types of behaviour with changing either µ, α or N . In particular, for fixed µ, α system (1) will almost certainly become unstable for N sufficiently large.
Despite the simplistic character of May's model [5], his pioneering work gave rise to a long standing 'stability versus diversity' debate, which is not fully settled yet [4,6,7,8,9], and played a fundamental role in theoretical ecology by prompting ecologists to think about special features of real multi-species ecosystems that help such systems to remain stable. Variations of May's model are still being discussed nowadays in the context of neighbourhood stability, see [4] and references therein.
One obvious limitation of the neighbourhood stability analysis is that it gives no insight into the model behaviour beyond the instability threshold. Hence May's model has only limited bearing on the dynamics of populations operating outof-equilibrium. An instability does not necessarily imply lack of persistence: populations could coexist thanks to limit cycles or chaotic attractors, which typically originate from unstable equilibrium points. Important questions posing serious challenge then relate to classification of equilibria by stability, studying their basins of attraction, and other features of global dynamics [4]. Over the last years, as the computing power grew, non-linear models have increasingly been used to investigate population dynamics on the global scale by means of numerical integration of the corresponding system of ODEs [10,11,12,13]. Although such investigations captured a rich variety of types of behaviour such as fold bifurcations when points of equilibrium merge/annihilate [14], or various types of chaotic dynamics [15,16], they provide little analytic insight and are limited to small to medium-sized systems.
In this paper we attempt to investigate the generic properties of the global dynamics of large complex multi-species systems by retaining in our model only the bare essentialsnonlinearity and stochasticity. Much in the spirit of May's original approach, the model we propose is simple enough to allow for an analytic treatment yet at the same time is rich enough to exhibit a non-trivial behaviour. In particular, our model captures an instability transition of the May-Wigner type, but now on the global scale. It also sheds additional light on the nature of this transition by relating it to an exponential explosion in the number of equilibria. Interestingly, despite the nonlinear setting of the problem the random matrix ideas again play a central role in our analysis.
Similar to the May's linear model our toy model is likely to have rather limited practical significance for quantitative description of real ecosystems but it might provide insight into the generic qualitative features of such systems and beyond, e.g. machine learning [17] or financial ecosystems [18,19]. The idea of destabilisation by interaction is of relevance far beyond the mathematical ecology context as applications of systems of many coupled non-linear ODEs are vast (e.g. complex gene regulatory networks [20], neural networks [21,22], or random catalytic reaction networks [23]). where µ > 0 and the components fi(x) of the vector field f = (f1, . . . , fN ) are zero mean random functions of the state vector x = (x1, . . . , xN ). To put this model in the context of the discussion above, if xe is an equilibrium of (2), i.e., if −µxe + f (xe) = 0 then, in the immediate neighbourhood of xe system (2) reduces to May's model (1) with y = x − xe and J jk = (∂fj/∂x k )(xe). The non-linear system (2) may have multiple equilibria whose number and locations depend on the realisation of the random field f(x). To visualise the global picture, it is helpful to consider first a special case of a gradient-descent flow, characterised by the existence of a potential function V (x) such that f = −∇V . In this case, system (2) can be rewritten as dx/dt = −∇L, with L(x) = µ|x| 2 /2 + V (x) being the associated Lyapunov function describing the effective landscape. In the domain of L, the state vector x(t) moves in the direction of the steepest descent, i.e., perpendicular to the level surfaces L(x) = h towards ever smaller values of h. This provides a useful geometric intuition. The term µ|x| 2 /2 represents the globally confining parabolic potential, i.e., a deep well on the surface of L(x), which does not allow x to escape to infinity. At the same time the random potential V (x) may generate many local minima of L(x) (shallow wells) which will play the role of attractors for our dynamical system. Moreover, if the confining term is strong enough then the full landscape will only be a small perturbation of the parabolic well, typically with a single stable equilibrium located close to x = 0. In the opposite case of relatively weak confining term, the disorderdominated landscape will be characterised by a complicated random topology with many points of equilibria, both stable and unstable. Note that in physics, complicated energy landscapes is a generic feature of glassy systems with intriguingly slow long-time relaxation and non-equilibrium dynamics, see e.g. [24].

Model
The above picture of a gradient-descent flow is however only a very special case since the generic systems of ODEs (2) are not gradient. The latter point can easily be understood in the context of model ecosystems. For, by linearising a gradient flow in a vicinity of any equilibrium, one always obtains a symmetric community matrix, whilst the community matrices of model ecosystems are in general asymmetric. Note also a discussion of an interplay between non-gradient dynamics in random environment and glassy behaviour in [25].
To allow for a suitable level of generality we therefore suggest to choose the N −dimensional vector field f (x) as a sum of 'gradient' and non-gradient ('solenoidal') contributions: where we require the matrix A(x) to be antisymmetric: Aij = −Aji. The meaning of this decomposition is that vector fields can be generically divided into a conservative irrotational component, sometimes called 'longitudinal', whose gradient connects the attractors or repellers and a solenoidal curl field, also called 'transversal'. As discussed in, e.g., [26] such a representation is closely related to the so-called Hodge decomposition of differential forms and generalises the well-known Helmholtz decomposition of the 3−dimensional vector fields into curl-free and divergence-free parts to higher dimensions. Correspondingly, we will call V (x) the scalar potential and the matrix A(x) the vector potential. The normalising factor 1/ √ N in front of the sum on the right-hand side in (3) ensures that the transversal and longitudinal parts of f (x) are of the same order of magnitude for large N .
Finally, to make the model as simple as possible and amenable to a rigorous and detailed mathematical analysis we choose the scalar potential V (x) and the components Aij(x), i < j, of the vector potential to be statistically independent, zero mean Gaussian random fields, with smooth realisations and the additional assumptions of homogeneity (translational invariance) and isotropy reflected in the covariance structure: Aij(x)Anm(y) = a 2 ΓA |x − y| 2 (δinδjm − δimδjn) . [ 5 ] Here the angular brackets ... stand for the ensemble average over all realisations of V (x) and A(x), and δin is the Kronecker delta: δin = 1 if i = n and zero otherwise. For simplicity, we also assume that the functions ΓV (r) and ΓA(r) do not depend on N . This implies [27] where the 'radial spectral' densities γσ(s) ≥ 0 have finite total mass: ∞ 0 γσ(s)ds < ∞. We normalize these densities by is a dimensionless measure of the relative strengths of the longitudinal and transversal components of f (x): if τ = 0 then f (x) is divergence free and if τ = 1 it is curl free.

Results
Determining and classifying all points of equilibria of a dynamical system with many degrees of freedom is a well-known formidable analytical and computational problem. In this paper we shall focus our investigation on the simplest, yet informative characteristic of system (2) by counting its total number of equilibria, that is the total number Ntot of solutions of the simultaneous equations − µxi + fi(x1, . . . , xN ) = 0, i = 1, . . . , N.
[ 6 ] Certainly, finding Ntot is a good starting point of any phase portrait analysis. Had we restricted ourselves to the gradient-descent flows, Ntot would simply count the number of stationary points (minima, maxima, or saddle-points) on the surface of the Lyapunov function L(x). The problem of counting and classifying stationary points of high-dimensional random energy landscapes of various types attracted considerable interest in recent years [28,29,30,31,32,33]. In particular, works [28,29] study such energy landscapes generated by a potential equivalent to the above Lyapunov function. One of the main conclusions of that study is that for N large the topology of the Lyapunov function changes drastically with decrease of the strength of the confining term relative to that of the interaction term in L(x). The change manifests itself in the emergence of multitude of equilibria, exponential in number. Such a transition is intimately connected to the spin-glass like restructuring of the Boltzmann-Gibbs measure induced by the Lyapunov function when the latter is treated as an effective energy landscape.
We shall prove below that for N large the general autonomous system (2)-(3) exhibits a similar drastic change in the total number of equilibria when the control parameter drops below the threshold value mc = 1. As in the case of gradient systems, the proof involves the Kac-Rice formula as a starting point. However, performing the subsequent steps requires quite different mathematical techniques due to the asymmetry of the Jacobian matrix for non-gradient systems. The Kac-Rice formula, see e.g. [34], counts solutions of simultaneous algebraic equations. Under our assumptions (homogeneity, isotropy and Gaussianity of V and A), this formula yields the ensemble average of Ntot in terms of that of the modulus of the spectral determinant of the Jacobian matrix (Jij) N i,j=1 , Jij = ∂fi/∂xj (see Materials and Methods): thus bringing the original non-linear problem into the realms of the random matrix theory. The probability (ensemble) distribution of the matrix J can easily be determined in closed form. Indeed, the matrix entries of J are zero mean Gaussian variables and their covariance structure, at spatial point x, can be obtained from (4)-(5) by differentiation: where N = (1 − τ )/N . Thus, to leading order in the limit N → ∞, where Xij, i, j = 1, . . . , N are zero mean Gaussians with and ξ is a standard Gaussian, ξ ∼ N (0, 1), which is statistically independent of X = (Xij). Note that for the divergence free fields f (x) (i.e., if τ = 0) the entries of J are statistically independent in the limit N → ∞, exactly as in May's model. On the other side, if f (x) has a longitudinal component (τ > 0) then this implies positive correlation between the pairs of matrix entries of J symmetric about the main diagonal: XijXji = τ if i = j. Such distributions of the community matrix has also been used in the neighbourhood stability analysis of model ecosystems [8]. Finally, in the limiting case of curl free fields (τ = 1), the matrix J is real symmetric.
The representation (8) comes in handy as it allows one to express (7) as a random matrix integral: ) and the angle brackets . . . X N stand for averaging over the real elliptic ensemble of random N × N matrices X defined in (9), see also (20). This oneparameter family of random matrices interpolates between the Gaussian Orthogonal Ensemble of real symmetric matrices (GOE, τ = 1) and real Ginibre ensemble of fully asymmetric matrices (rGinE, τ = 0), see [35] for discussions. Both rGinE and its one-parameter extension (9) have enjoyed considerable interest in the literature in recent years [36,37,38,39,40].
The matrix X is asymmetric (unless τ = 1) and can have real as well as complex eigenvalues. The latter come in complex-conjugate pairs. Their density, in the limit N → ∞, is constant inside the ellipse with the main half-axis √ N (1±τ ) and vanishes sharply outside [41,39,35]. The corresponding theorem is known as the Elliptic Law and its validity extends beyond the Gaussian matrix distributions [42,43]. However, in the context of our investigation it is the density of real eigenvalues of X that appears to be most relevant.
Denote by ρ (r) N (x) the density of real eigenvalues of N × N matrices X (9) averaged over all realisations of X. It is convenient to normalize ρ  N (x) is directly related to the averaged value of the modulus of the determinant that appears in (10). Namely, where CN (τ ) = 2 √ 1 + τ (N − 1)!/(N − 2)!! and ρ (r) N +1 (x) is the average density of real eigenvalues of matrices X of size (N + 1) × (N + 1). For the limiting case τ = 0 this relation appeared originally in [44], and it can be extended to any τ ∈ [0, 1) without much difficulty (see SI for a derivation of (11) following the approach of [45] ). In the limiting case of real symmetric matrices τ = 1, all eigenvalues of X are real and relation (11) is also valid [28].
Combining (10) and (11) and changing the variable of integration from t to λ = m + t √ τ , one can express Ntot for system (2) with N degrees of freedom in terms of the density of real eigenvalues in the elliptic ensemble of random matrices (9) of size (N + 1) × (N + 1): where S(λ) The importance of this relation is due to the fact that ρ (r) N (x) is known in closed form in terms of Hermite polynomials [39]. This allows us to carry out an asymptotic evaluation of the integral in (12) and calculate Ntot in the limit N → ∞. The key finding that emerges from this calculation is that Ntot changes drastically around m = 1. If m > 1 then lim N →∞ Ntot = 1.
Close to the phase transition point m = 1 the complexity exponent vanishes quadratically, Σtot = (1 − m) 2 as m → 1, implying that the width of the transition region around m = 1 is 1/ √ N . According to the general lore of phase transitions, for large but finite N there exists a 'critical regime' m = 1 + κN −1/2 where the number of equilibria changes smoothly between the two 'phases' (13) and (14). A quick inspection of (12) shows that the corresponding crossover profile is determined by the profile of ρ After rescaling λ, λ with erf(x) = 1−erfc(x) = 2 √ π x 0 e −t 2 dt. In terms of ρ (r) edge (ζ), the critical crossover profile is given by where cτ = τ /(1 − τ ). The right-hand side of (16) interpolates smoothly between the two regimes (13) and (14), when parameter κ runs from κ = −∞ to κ = +∞. Although our investigation is concerned with the ensemble average of the number of equilibria Ntot, we expect that in the limit N → ∞ the deviations of Ntot from its average Ntot are relatively small. This is certainly the case above the critical threshold, for m > 1. For, under some additional technical assumptions on the decay of correlations for f (x), the system (2) will almost certainly have at least one stationary point, see [48] for the relevant results about the maxima of homogeneous Gaussian fields. Therefore, Ntot ≥ 1 and the established convergence of Ntot to 1 in the limit N → ∞ actually implies that the probability for Ntot to take other values than one is close to zero for large N . The problem of estimating the deviation of Ntot from its average value in the opposite regime 0 < m < 1 is much harder and is an open and interesting question 1 .

Discussion
In this paper we introduced a model describing generic large complex systems and examined the dependence of the total number of equilibria in such systems on the system complexity as measured by the number of degrees of freedom and the interaction strength. The inspiration for our work came from May's pioneering study [1] of the neighbourhood stability of large model ecosystems. Our outlook is complementary to that of May's in that it adopts a global point of view which is not limited to the neighbourhood of the presumed equilibrium.
In the context of model ecosystems our analysis is applicable to complex multi-species communities in which each kind of species on its own becomes extinct and thus interaction between species is key to persistence of the community. The key feature of our analysis is that in the presence of interactions, as the complexity increases, there is an abrupt change from a simple set of equilibria (typically, a single equilibrium for large number of species N 1) to a complex set of equilibria, with their total number growing exponentially with N . In the latter regime, we expect the stable equilibria to be only a tiny proportion of all the multitude of equilibria, see discussion below, which is indicative of long relaxation times and transient non-equilibrium behaviour.
We expect this sharp transition in the phase portrait to be shared by other systems of randomly coupled autonomous ODE's with large number of degrees of freedom. To that end, it is appropriate to mention that very recently a similar 'explosion in complexity' was reported in a model of neural network consisting of randomly interconnected neural units [22]. The model considered in [22] is essentially of the form (2) but with the particular choice of fi = j JijS(xj) where S is an odd sigmoid function representing the synaptic nonlinearity and Jij are independent centred Gaussian variables representing the synaptic connectivity between neuron i and j. Although being Gaussian, the corresponding (non-gradient) vector field is not homogeneous and thus seems rather different from our choice and not easily amenable to a rigorous analysis. Nevertheless, a shrewd semi-heuristic analysis of [22] revealed that close to the critical coupling threshold the two models actually display very similar behaviour, described essentially by the same exponential growth in the total number of equilibria with rate Σtot(m). This fact points towards considerable universality of the transition from (13) to (14) and suggests that the crossover function (16) may be universal as well.
Our model captures an abrupt change in the dynamics of large complex systems on the macroscopic scale. At the same time zooming in to classify each and every equilibrium point into locally stable or unstable seems a hard task. For, although linearising the field f (x) around a given equilibrium is fairly straightforward, with the outcome being the Jacobian matrix (8), conditioning by the positions of equilibria and taking into account all eventualities is a highly non-trivial task. Given the stochastic setup of our model the question about stability of individual equilibria may be even the wrong question to ask, whereas addressing the statistics of the number of stable equilibria seems very appropriate.
Arguments similar to those in the previous section yield the ensemble average of the total number of stable equilibria, Nst , over all realisations of the vector field f (x) in terms the random matrix integral (cf.(10)): [ 17 ] where χx(X) = 1, if all N eigenvalues of matrix X have real parts less than the spectral parameter x = √ N (m+t √ τ ), and χx(X) = 0 otherwise. In the limiting case of a purely gradient dynamics τ = 1, the rescaled Jacobian matrix X is real symmetric with all N eigenvalues real. In this case the above integral can be related to the probability density of the maximal eigenvalue of the GOE matrix [29,30], with the latter being a well-studied object in the random matrix theory, see e.g. [50] and references therein. This observation can then be used to evaluate Nst for N 1. One finds [29] that Nst → 1 if m > 1, whilst if 0 < m < 1 then, to leading order in N , Nst ∝ e N Σ st , with 0 < Σst < Σtot. Thus, in the case of purely gradient dynamics, as the complexity increases, large nonlinear autonomous systems assembled at random undergo an abrupt change from a typical phase portrait with a single stable equilibrium to a phase portrait dominated by an exponential number of unstable equilibria with an admixture of a smaller, but still exponential in N , number of stable equilibria.
It was suggested to us by J.-P. Bouchaud that in the general case of non-gradient dynamics 0 ≤ τ < 1, it would be natural to expect a further phase transition in the plane (m, τ ) such that below a certain number τc(m) stable equilibria are no longer exponentially abundant in the limit N → ∞ (i.e. Σst(m, τ ) → 0), with further implications for the global dynamics. Unfortunately, for a fixed 0 ≤ τ < 1 only vanishing fraction of order of N −1/2 of eigenvalues of X remain real, and the relation of the integral in Eq. (17) to statistics of the largest real eigenvalue in the elliptic ensemble seems to be lost. This fact prevented us so far from reliable counting of stable equilibria for the general case of non-gradient flows. In principle, for given values of parameters N, τ, m one may attempt to evaluate the ensemble average in the integral in Eq. (17) numerically, and then to evaluate numerically the integral itself. Although such a procedure seems tractable, its actual implementation with sufficient precision is not straightforward, especially in the limit N → ∞ due to the exponentially large values involved. Clarification of the status of the picture suggested by J.-P. Bouchaud and related studies remain an important outstanding issue and is left for a future work.

Materials and Methods
Kac-Rice Formula The expected number Ntot of simultaneous solutions to the system of equations (6) in R N is given by the formula, see e.g. [34], [ 18 ] where δ(x) is the multivariate Dirac δ-function, dx N is the volume element in R N and J ij (x) = ∂f i /∂x j are matrix elements of the Jacobian matrix J = (J ij ). By our assumptions (4) -(5) the random field f (x) is homogeneous and isotropic. For such fields samples of f and J taken in one and the same spatial point x are uncorrelated, f l · ∂f i /∂x j = 0 for all i, j, l. This is well known and can be checked by straightforward differentiation. In addition, the field f is Gaussian, hence the f (x) and J(x) are actually statistically independent. This simplifies the evaluation of the integral in (18) considerably. Indeed, the statistical average in (18) factorizes and, and since | det (−µδ ij + J ij (x)) | does not vary with x, the integrand effectively reduces to . [ 19 ] Furthermore, at every spatial point x the vector f (x) is Gaussian with uncorrelated and identically distributed components, Therefore e ik·f (x) = e −σ 2 |k| 2 /2 , and evaluating the integral on the right-hand side in (19) one arrives at (7).
Real elliptic matrices and asymptotics of N tot The joint probability density function P N (X) of the matrix entries in the elliptic ensemble of real Gaussian random matrices X of size N × N is given by Tr XX T − τ X 2 , [ 20 ] where Z N is the normalisation constant and τ ∈ [0, 1). It is straightforward to verify that the covariance of matrix entries X ij is given by the expression specified in (9). The mean density of real eigenvalues of ρ (r) N (x) in the elliptic ensemble (20) is known in closed form in terms of Hermite polynomials, see [39]. Assuming for simplicity that N + 1 is even, one has ρ (r) [ 21 ] and ρ (r),2 N −1 (u) du. [ 22 ] Here ψ (τ ) This, together with the integral (12) allow one to evaluate N tot in the limit N → ∞. We shall sketch the corresponding evaluation below.
The asymptotics of ρ (r) N (x) in the bulk and at the edge of the support of the distribution of real eigenvalues in the real elliptic ensemble were found in [39], and outside the support it can also be readily extracted using (21)- (22). In particular, in the bulk, i.e., for |x| < (1 + τ ) √ N , the contribution of (21) to ρ (r) N (x) is dominant and, to leading order in N , [ 23 ] At the same time for |x| > (1 + τ ) (21) and (22) The form of (12) suggest the application of the Laplace method. One easily finds that S(λ) has a minimum at λ * = m(1+τ ) which belongs to the domain |λ| < 1+τ as long as 0 < m < 1.Thus, applying the Laplace method and taking into account the asymptotic formula [ 26 ] one arrives at the asymptotic expression (14) for N tot in the parameter range 0 < m < 1. For m > 1 the saddle-point occurs in the domain λ > 1 + τ so that the analysis requires search for the minimum of S(λ) + Ψ(λ). After straighforward algebra we find d is equal to zero at λ = λ * = m + τ m > 1 + τ . One also verifies that this is a point of minimum for S(λ) + Ψ(λ) and a further simple calculation then yields S(λ * ) + Ψ(λ * ) = − ln m/ √ τ . Calculating the saddle-point contribution then yields (13).
The above asymptotic analysis assumes that 0 ≤ τ < 1. Let us now discuss the modifications required to study the scaling regime of weakly non-gradient flow τ → 1 for 0 < m < 1. We only need to evaluate the leading contribution to ρ (r) N +1 given by (21). By making use of the above integral representation for the scaled Hermite polynomials h (τ ) k (x) and applying the scaling τ = 1 − u 2 N , we can write ρ (r),1 with I N (λ) given by where Φ N (a) = e −a N k=0 a k /k!. Recalling that the limit of Φ N (a) as N → ∞ is 1 if 0 < a < 1 and 0 if a > 1, one obtains ρ (r) Substituting this expression into the integrand of (12) and evaluating the integral in the limit N → ∞ (hence, τ → 1) by the Laplace method then yields N tot in the weakly non-gradient regime.
Finally, our calculation of the profile of N tot in the transitional region m = 1 + κN −1/2 uses the fact that in such a regime the main contribution to the integral (12) comes from the neighbourhood of the spectral edge, , where we have, to the leading order in N , This together with (26) and (15) converts (12) to (16).

Appendix: Supporting Information
In this appendix we express the mean density of real eigenvalues in the real elliptic ensemble, see Equation (35) below, of N × N matrices XN in terms of the modulus of the spectral determinant of real elliptic matrices of size (N − 1) × (N − 1). Our derivation is based on the method suggested in [44], however our Jacobian computation differs from the one given in [44] and may be of interest on its own. For a similar calculation in the context of complex matrices see [45].

Housholder reflections and partial triangularization of real matrices
The key idea of [44] is based on employing Householder reflections described by matrices where IN is N × N identity matrix, v is column vector in R N of unit length. Its transpose, v T , is row vector, so that the Kronecker product v ⊗ v T is a matrix. Pv describes the reflection about the hyperplane with normal v and passing through the origin. Obviously, Pv is is symmetric and orthogonal: P T v = Pv and P 2 v = IN . Let e1 be the first vector of the standard basis in R N , i.e., e T 1 = (1, 0, . . . , 0). For any unit vector x, |x| = 1 define v = x + e1 2(1 + x1) . [ 28 ] and consider the Hausholder reflection Pv built from the above v according to (27). Then it is easy to check that Pvx = −e1. This implies that for any nonzero vector x there exists a Housholder reflection such that Pvx = ke1, with k = −|x|.
Let λ be a real eigenvalue of the N × N matrix A (N ) with real entries A (N ) ij , i.e. A (N ) x = λx for some column N −vector x of unit length. Our goal is to demonstrate that it is always possible to represent that matrix as Considering now the volume element dA (N ) = N i,j dA (N ) ij our next goal to write it down in terms of N 2 independent variables parametrizing the right-hand side of (29), that is (N −1) 2 variables parametrizing A (N −1) , N − 1 components of w, 1 parameter for λ, and the remaining N − 1 parameters for representing the matrix P . A convenient parametrization for P comes from employing (28) for the vector v, which shows that the last N − 1 components of that vector (v2, . . . , vN ) T ≡ q can be used as independent variables, whereas normalization fixes the first component. Writing v T = 1 − q T q, q with |q| < 1 and employing (27) yields an explicit parametrization : The problem therefore amounts to calculating the Jacobian of the transformation A (N ) → (λ, w, q, A (N −1) ). To that end we start with differentiating (29) which gives dA (N ) = P (P dP ), where we employed that dP P = −P dP and used the notation [A, B] = AB − BA for the matrix commutator. A direct calculation using (30) shows that the matrix (P dP ) can be symbolically written as where the expression for dF is immaterial for our goals, and Substituting (31) into the expression for dA (N ) we find that From this expression we easily read off the required Jacobian to be given by Jacobian = | det (λIN−1 − A (N −1) )| ∂b ∂q , [ 33 ] where the last factor symbolically denotes the part of the Jacobian coming from the transformation db → dq described in (32). A straightforward calculation shows that so that finally we arrive at the change-of-measure formula Elliptic Ensemble of Gaussian Random Matrices. The Joint Probability Density (JPD) of the Elliptic Ensemble of Gaussian random matrices XN of size N × N whose entries have covariance XijXnm = δinδjm + τ δjnδim is given by Tr(XN X T N ) − τ Tr(X 2 N ) , [ 35 ] where ZN is the corresponding normalization factor ZN = 2 N/2 π N (N +1)/2 (1 + τ ) N (N +1)/4 (1 − τ ) N (N −1)/4 .
Our goal is to find the JPD of variables λ, w, q, XN−1 used to perform the partial triangulation of XN via (29), with the role of A (N ) played now by XN . We have Tr XN X T N = λ 2 +w T w+Tr XN−1X T N −1 , Tr X 2 N = λ 2 +Tr X 2 N −1 . Taking into account the change-of-measure formula (34) we see that the corresponding JPD can be written as P(λ, w, q, XN−1) = 2 N (1 − q T q) By definition, the density of real eigenvalues ρ (r) N (λ) is obtained by integrating the above JPD over variables |q| < 1, −∞ < w < ∞ and finally over XN−1. After performing the integrals over q and w we immediately arrive at the relation where CN−1(τ ) is a certain constant which can be found explicitly. This is precisely equivalent to Equation (13) in the main text with the obvious change of notations: λ → x and N → N + 1.