Nonlinear Analogue of the May-Wigner Instability Transition

Model

Consider a system of $N$ coupled non-linear autonomous ODEs of the form

where $\mu>0$ and the components $f_{i}({\textbf{x}})$ of the vector field $\textbf{f}=(f_{1},\ldots,f_{N})$ are zero mean random functions of the state vector $\textbf{x}=(x_{1},\ldots,x_{N})$. To put this model in the context of the discussion above, if $\mathbf{x}_{e}$ is an equilibrium of (2), i.e., if $-\mu\mathbf{x}_{e}+\mathbf{f}(\mathbf{x}_{e})=\mathbf{0}$ then, in the immediate neighbourhood of $\mathbf{x}_{e}$ system (2) reduces to May’s model (1) with $\mathbf{y}=\mathbf{x}-\mathbf{x}_{e}$ and $J_{jk}=(\partial f_{j}/\partial x_{k})(\mathbf{x}_{e})$.

The non-linear system (2) may have multiple equilibria whose number and locations depend on the realisation of the random field $\textbf{f}(\textbf{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(\mathbf{x})$ such that $\mathbf{f}=-\nabla V$. In this case, system (2) can be rewritten as $d\mathbf{x}/dt=-\nabla L$, with $L(\mathbf{x})=\mu{|{\bf x}|^{2}}/{2}+V({\bf x})$ being the associated Lyapunov function describing the effective landscape. In the domain of $L$, the state vector $\mathbf{x}(t)$ moves in the direction of the steepest descent, i.e., perpendicular to the level surfaces $L(\mathbf{x})=h$ towards ever smaller values of $h$. This provides a useful geometric intuition. The term $\mu{|{\bf x}|^{2}}/{2}$ represents the globally confining parabolic potential, i.e., a deep well on the surface of $L(\mathbf{x})$, which does not allow $\mathbf{x}$ to escape to infinity. At the same time the random potential $V(\mathbf{x})$ may generate many local minima of $L(\textbf{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 ${\mathbf{x}}=0$. In the opposite case of relatively weak confining term, the disorder-dominated 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. .

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 .

To allow for a suitable level of generality we therefore suggest to choose the $N-$dimensional vector field $\mathbf{f}(\mathbf{x})$ as a sum of ‘gradient’ and non-gradient (‘solenoidal’) contributions:

where we require the matrix $A({\bf x})$ to be antisymmetric: $A_{ij}=-A_{ji}$. 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., 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({\bf x})$ the scalar potential and the matrix $A({\bf x})$ the vector potential. The normalising factor $1/\sqrt{N}$ in front of the sum on the right-hand side in (3) ensures that the transversal and longitudinal parts of $\mathbf{f}({\bf 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({\bf x})$ and the components $A_{ij}({\bf 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:

Here the angular brackets $\langle...\rangle$ stand for the ensemble average over all realisations of $V(\mathbf{x})$ and $A(\mathbf{x})$, and $\delta_{in}$ is the Kronecker delta: $\delta_{in}=1$ if $i=n$ and zero otherwise.

For simplicity, we also assume that the functions $\Gamma_{V}(r)$ and $\Gamma_{A}(r)$ do not depend on $N$. This implies

where the ‘radial spectral’ densities $\gamma_{\sigma}(s)\geq 0$ have finite total mass: $\int_{0}^{\infty}\gamma_{\sigma}(s)ds<\infty$. We normalize these densities by requiring that $\Gamma^{\prime\prime}_{\sigma}(0)=\int_{0}^{\infty}s^{2}\gamma_{\sigma}(s)ds=1$. The ratio

is a dimensionless measure of the relative strengths of the longitudinal and transversal components of $\mathbf{f}(\textbf{x})$: if $\tau=0$ then $\mathbf{f}(\textbf{x})$ is divergence free and if $\tau=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 ${\cal N}_{tot}$ of solutions of the simultaneous equations

Certainly, finding ${\cal N}_{tot}$ is a good starting point of any phase portrait analysis.

Had we restricted ourselves to the gradient-descent flows, ${\cal N}_{tot}$ would simply count the number of stationary points (minima, maxima, or saddle-points) on the surface of the Lyapunov function $L({\bf x})$. The problem of counting and classifying stationary points of high-dimensional random energy landscapes of various types attracted considerable interest in recent years . In particular, works 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(\mathbf{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 $m_{c}=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. , counts solutions of simultaneous algebraic equations. Under our assumptions (homogeneity, isotropy and Gaussianity of $V$ and $A$), this formula yields the ensemble average of ${\cal N}_{tot}$ in terms of that of the modulus of the spectral determinant of the Jacobian matrix $(J_{ij})_{i,j=1}^{N}$, $J_{ij}={\partial f_{i}}/{\partial x_{j}}$ (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 $\epsilon_{N}=(1-\tau)/N$. Thus, to leading order in the limit $N\to\infty$,

where $X_{ij}$, $i,j=1,\ldots,N$ are zero mean Gaussians with

and $\xi$ is a standard Gaussian, $\xi\sim N(0,1)$, which is statistically independent of $X=(X_{ij})$. Note that for the divergence free fields $\mathbf{f}(\mathbf{x})$ (i.e., if $\tau=0$) the entries of $J$ are statistically independent in the limit $N\to\infty$, exactly as in May’s model. On the other side, if $\mathbf{f}(\mathbf{x})$ has a longitudinal component ($\tau>0$) then this implies positive correlation between the pairs of matrix entries of $J$ symmetric about the main diagonal: $\left\langle X_{ij}X_{ji}\right\rangle=\tau$ if $i\not=j$. Such distributions of the community matrix has also been used in the neighbourhood stability analysis of model ecosystems . Finally, in the limiting case of curl free fields ($\tau=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:

where $x=\sqrt{N}(m+t\sqrt{\tau})$ and the angle brackets $\langle\ldots\rangle_{X_{N}}$ stand for averaging over the real elliptic ensemble of random $N\times N$ matrices $X$ defined in (9), see also (20). This one-parameter family of random matrices interpolates between the Gaussian Orthogonal Ensemble of real symmetric matrices (GOE, $\tau=1$) and real Ginibre ensemble of fully asymmetric matrices (rGinE, $\tau=0$), see for discussions. Both rGinE and its one-parameter extension (9) have enjoyed considerable interest in the literature in recent years .

The matrix $X$ is asymmetric (unless $\tau=1$) and can have real as well as complex eigenvalues. The latter come in complex-conjugate pairs. Their density, in the limit $N\to\infty$, is constant inside the ellipse with the main half-axis $\sqrt{N}(1\pm\tau)$ and vanishes sharply outside . The corresponding theorem is known as the Elliptic Law and its validity extends beyond the Gaussian matrix distributions . However, in the context of our investigation it is the density of real eigenvalues of $X$ that appears to be most relevant.

Denote by $\rho_{N}^{(r)}(x)$ the density of real eigenvalues of $N\times N$ matrices $X$ (9) averaged over all realisations of $X$. It is convenient to normalize $\rho_{N}^{(r)}(x)$ in such a way that $\int_{\alpha}^{\beta}\rho_{N}(x)\,dx$ gives the average number of real eigenvalues of $X$ in the interval $[\alpha,\beta]$. A crucial observation is that $\rho_{N}^{(r)}(x)$ is directly related to the averaged value of the modulus of the determinant that appears in (10). Namely,

where ${\cal C}_{N}(\tau)=2\sqrt{1+\tau}\,(N-1)!/(N-2)!!$ and $\rho^{(r)}_{N+1}(x)$ is the average density of real eigenvalues of matrices $X$ of size $(N+1)\times(N+1)$. For the limiting case $\tau=0$ this relation appeared originally in , and it can be extended to any $\tau\in[0,1)$ without much difficulty (see SI for a derivation of (11) following the approach of ). In the limiting case of real symmetric matrices $\tau=1$, all eigenvalues of $X$ are real and relation (11) is also valid .

Combining (10) and (11) and changing the variable of integration from $t$ to $\lambda=m+t\sqrt{\tau}$, one can express $\left\langle{\cal N}_{tot}\right\rangle$ 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)\times(N+1)$:

where $S(\lambda)=\frac{(\lambda-m)^{2}}{2\tau}-\frac{\lambda^{2}}{2(1+\tau)}$ and ${\cal K}_{N}(\tau)={N^{\frac{-N+1}{2}}{\cal C}_{N}(\tau)}/\!{\sqrt{\tau}}$. The importance of this relation is due to the fact that $\rho_{N}^{(r)}(x)$ is known in closed form in terms of Hermite polynomials . This allows us to carry out an asymptotic evaluation of the integral in (12) and calculate $\left\langle{\cal N}_{tot}\right\rangle$ in the limit $N\to\infty$. The key finding that emerges from this calculation is that $\left\langle{\cal N}_{tot}\right\rangle$ changes drastically around $m=1$. If $m>1$ then

On the other hand, if $0<m<1$ then, to leading order in the limit $N\to\infty$,

where $\Sigma_{tot}(m)=\frac{1}{2}(m^{2}-1)-\ln{m}>0$ for all $0<m<1$. Therefore, if $m<1$ then $\left\langle{\cal N}_{tot}\right\rangle$ grows exponentially with $N$. The factor in front of the exponential in (14) is given by $\gamma_{\tau}=\sqrt{{2(1+\tau)}/{(1-\tau)}}$ as long as $\tau<1$. The gradient limit $\tau=1$ can be approached by scaling $\tau$ with $N$. Setting $\tau=1-\frac{u^{2}}{N},\,0\leq u<\infty$, one obtains $\gamma_{\tau}=4\sqrt{\frac{N}{\pi}}\,\int_{0}^{\sqrt{1-m^{2}}}\,e^{-u^{2}p^{2}}dp$. This regime describes a weakly non-gradient flow. The corresponding regime for ensembles of asymmetric matrices was discovered long ago .

Close to the phase transition point $m=1$ the complexity exponent vanishes quadratically, $\Sigma_{tot}=(1-m)^{2}$ as $m\to 1$, implying that the width of the transition region around $m=1$ is $1/\sqrt{N}$. According to the general lore of phase transitions, for large but finite $N$ there exists a ‘critical regime’ $m=1+\kappa 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 $\rho_{N}^{(r)}(x)$ in the vicinity of the ‘spectral edge’ $x=(1+\tau)\sqrt{N}$, see Materials and Methods. After rescaling $\lambda$, $\lambda=1+\tau+\frac{\zeta\sqrt{1-\tau^{2}}}{\sqrt{N}}$, the density $\rho_{N}^{(r)}(\lambda\sqrt{N})$ converges to $\frac{1}{\sqrt{1-\tau^{2}}}\rho_{edge}^{(r)}(\zeta)$ in the limit $N\to\infty$, where :

with $\operatorname{erf}(x)=1-\operatorname{erfc}(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-t^{2}}\,dt$. In terms of $\rho_{edge}^{(r)}(\zeta)$, the critical crossover profile is given by

where $c_{\tau}=\sqrt{\tau/(1-\tau)}$. The right-hand side of (16) interpolates smoothly between the two regimes (13) and (14), when parameter $\kappa$ runs from $\kappa=-\infty$ to $\kappa=+\infty$.

Although our investigation is concerned with the ensemble average of the number of equilibria ${\cal N}_{tot}$, we expect that in the limit $N\to\infty$ the deviations of ${\cal N}_{tot}$ from its average $\langle{\cal N}_{tot}\rangle$ 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 $\mathbf{f}(\mathbf{x})$, the system (2) will almost certainly have at least one stationary point, see for the relevant results about the maxima of homogeneous Gaussian fields. Therefore, ${\cal N}_{tot}\geq 1$ and the established convergence of $\langle{\cal N}_{tot}\rangle$ to 1 in the limit $N\to\infty$ actually implies that the probability for ${\cal N}_{tot}$ to take other values than one is close to zero for large $N$. The problem of estimating the deviation of ${\cal N}_{tot}$ from its average value in the opposite regime $0<m<1$ is much harder and is an open and interesting questionIn this context we would like to mention the recent work of Subag who proved that the deviations of ${\cal N}_{tot}$ from $\langle{\cal N}_{tot}\rangle$ in the spherical $p$-spin model are negligible in the limit of large system size. Though that model is different from ours, it is not dissimilar to the gradient limit of $\tau=1$ of our model , for instance, the average number of equilibria grows exponentially with $N$ . Thus one might hope to adopt the technique of to our model. Another relevant reference is ..

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 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\gg 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 . The model considered in is essentially of the form (2) but with the particular choice of $f_{i}=\sum_{j}J_{ij}S(x_{j})$ where $S$ is an odd sigmoid function representing the synaptic nonlinearity and $J_{ij}$ 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 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 $\Sigma_{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 $\mathbf{f}(\mathbf{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, $\langle{\cal N}_{st}\rangle$, over all realisations of the vector field $\mathbf{f}(\mathbf{x})$ in terms the random matrix integral (cf.(10)):

where $\chi_{x}(X)=1$, if all $N$ eigenvalues of matrix $X$ have real parts less than the spectral parameter $x=\sqrt{N}(m+t\sqrt{\tau})$, and $\chi_{x}(X)=0$ otherwise. In the limiting case of a purely gradient dynamics $\tau=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 , with the latter being a well-studied object in the random matrix theory, see e.g. and references therein. This observation can then be used to evaluate $\left\langle{\cal N}_{st}\right\rangle$ for $N\gg 1$. One finds that $\langle{\cal N}_{st}\rangle\to 1$ if $m>1$, whilst if $0<m<1$ then, to leading order in $N$, $\langle{\cal N}_{st}\rangle\propto e^{N\Sigma_{st}}$, with $0<\Sigma_{st}<\Sigma_{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\leq\tau<1$, it would be natural to expect a further phase transition in the plane $(m,\tau)$ such that below a certain number $\tau_{c}(m)$ stable equilibria are no longer exponentially abundant in the limit $N\to\infty$ (i.e. $\Sigma_{st}(m,\tau)\to 0$), with further implications for the global dynamics. Unfortunately, for a fixed $0\leq\tau<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,\tau,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\to\infty$ 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.

Furthermore, at every spatial point $\mathbf{x}$ the vector $\mathbf{f}(\mathbf{x})$ is Gaussian with uncorrelated and identically distributed components,

Therefore $\langle e^{i\mathbf{k}\cdot\mathbf{f}(\mathbf{x})}\rangle=e^{-\sigma^{2}|\mathbf{k}|^{2}/2}$, and evaluating the integral on the right-hand side in (19) one arrives at (7).

Real elliptic matrices and asymptotics of $\langle{\cal N}\rangle_{tot}$ The joint probability density function ${\cal P}_{N}(X)$ of the matrix entries in the elliptic ensemble of real Gaussian random matrices $X$ of size $N\times N$ is given by

where ${\cal Z}_{N}$ is the normalisation constant and $\tau\in[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 $\rho_{N}^{(r)}(x)$ in the elliptic ensemble (20) is known in closed form in terms of Hermite polynomials, see . Assuming for simplicity that $N+1$ is even, one has $\rho_{N+1}^{(r)}(x)=\rho_{N+1}^{(r),1}(x)+\rho_{N+1}^{(r),2}(x)$ where

Here $\psi^{(\tau)}_{k}(x)=e^{-\frac{x^{2}}{2(1+\tau)}}h^{(\tau)}_{k}(x)$ and $h^{(\tau)}_{k}(x)$ are rescaled Hermite polynomials, $h^{(\tau)}_{k}(x)=\frac{1}{\sqrt{\pi}}\int_{-\infty}^{\infty}e^{-t^{2}}\left(x+it\sqrt{2\tau}\right)^{k}\,dt$. This, together with the integral (12) allow one to evaluate $\langle{\cal N}\rangle_{tot}$ in the limit $N\to\infty$. We shall sketch the corresponding evaluation below.

The asymptotics of $\rho_{N}^{(r)}(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 , and outside the support it can also be readily extracted using (21)-(22). In particular, in the bulk, i.e., for $|x|<(1+\tau)\sqrt{N}$, the contribution of (21) to $\rho_{N}^{(r)}(x)$ is dominant and, to leading order in $N$,

At the same time for $|x|>(1+\tau)\sqrt{N}$ both (21) and (22) yield exponentially small contributions to $\rho_{N+1}^{(r)}(x)$, with (22) being dominant. Our evaluation yields

The form of (12) suggest the application of the Laplace method. One easily finds that $S(\lambda)$ has a minimum at $\lambda_{*}=m(1+\tau)$ which belongs to the domain $|\lambda|<1+\tau$ as long as $0<m<1$.Thus, applying the Laplace method and taking into account the asymptotic formula

one arrives at the asymptotic expression (14) for $\langle{\cal N}\rangle_{tot}$ in the parameter range $0<m<1$. For $m>1$ the saddle-point occurs in the domain $\lambda>1+\tau$ so that the analysis requires search for the minimum of $S(\lambda)+\Psi(\lambda)$. After straighforward algebra we find $\frac{d}{d\lambda}\left(S(\lambda)+\Psi(\lambda)\right)=\frac{1}{2\tau}(\lambda+\sqrt{\lambda^{2}-4\tau})-\frac{m}{\tau}$ which is equal to zero at $\lambda=\lambda_{*}=m+\frac{\tau}{m}>1+\tau$. One also verifies that this is a point of minimum for $S(\lambda)+\Psi(\lambda)$ and a further simple calculation then yields $S(\lambda_{*})+\Psi(\lambda_{*})=-\ln{\left(m/\sqrt{\tau}\right)}$. Calculating the saddle-point contribution then yields (13).

The above asymptotic analysis assumes that $0\leq\tau<1$. Let us now discuss the modifications required to study the scaling regime of weakly non-gradient flow $\tau\to 1$ for $0<m<1$. We only need to evaluate the leading contribution to $\rho_{N+1}^{(r)}$ given by (21). By making use of the above integral representation for the scaled Hermite polynomials $h_{k}^{(\tau)}(x)$ and applying the scaling $\tau=1-\frac{u^{2}}{N}$, we can write

where $\Phi_{N}(a)=e^{-a}\sum_{k=0}^{N}{a^{k}}/{k!}$. Recalling that the limit of $\Phi_{N}(a)$ as $N\to\infty$ is 1 if $0<a<1$ and 0 if $a>1$, one obtains

Substituting this expression into the integrand of (12) and evaluating the integral in the limit $N\to\infty$ (hence, $\tau\to 1$) by the Laplace method then yields $\langle{\cal N}\rangle_{tot}$ in the weakly non-gradient regime.

Finally, our calculation of the profile of $\langle{\cal N}\rangle_{tot}$ in the transitional region $m=1+\kappa 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, $\lambda=1+\tau+\frac{\zeta\sqrt{1-\tau^{2}}}{\sqrt{N}}$, where we have, to the leading order in $N$,

This together with (26) and (15) converts (12) to (16).

References

Housholder reflections and partial triangularization of real matrices

The key idea of is based on employing Householder reflections described by matrices

and consider the Hausholder reflection $P_{\mathbf{v}}$ built from the above $\mathbf{v}$ according to (27). Then it is easy to check that $P_{\mathbf{v}}\mathbf{x}=-\mathbf{e}_{1}$. This implies that for any nonzero vector $\mathbf{x}$ there exists a Housholder reflection such that $P_{\mathbf{v}}\mathbf{x}=k\mathbf{e}_{1}$, with $k=-|\mathbf{x}|$.

Let $\lambda$ be a real eigenvalue of the $N\times N$ matrix $A^{(N)}$ with real entries $A^{(N)}_{ij}$, i.e. $A^{(N)}\mathbf{x}=\lambda\mathbf{x}$ for some column $N-$vector $\mathbf{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)}=\prod_{i,j}^{N}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 $\mathbf{w}$, $1$ parameter for $\lambda$, and the remaining $N-1$ parameters for representing the matrix $P$. A convenient parametrization for $P$ comes from employing (28) for the vector $\mathbf{v}$, which shows that the last $N-1$ components of that vector $(v_{2},\ldots,v_{N})^{T}\equiv\mathbf{q}$ can be used as independent variables, whereas normalization fixes the first component. Writing $\mathbf{v}^{T}=\left(\sqrt{1-\mathbf{q}^{T}\mathbf{q}},\mathbf{q}\right)$ with $|\mathbf{q}|<1$ and employing (27) yields an explicit parametrization :

The problem therefore amounts to calculating the Jacobian of the transformation $A^{(N)}\to(\lambda,\mathbf{w},\mathbf{q},A^{(N-1)})$. To that end we start with differentiating (29) which gives

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 $(PdP)$ 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

where the last factor symbolically denotes the part of the Jacobian coming from the transformation $d\mathbf{b}\to d\mathbf{q}$ described in (Housholder reflections and partial triangularization of real matrices). 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 $X_{N}$ of size $N\times N$ whose entries have covariance $\left\langle X_{ij}X_{nm}\right\rangle=\delta_{in}\delta_{jm}+\tau\delta_{jn}\delta_{im}$ is given by

where ${\cal Z}_{N}$ is the corresponding normalization factor

Our goal is to find the JPD of variables $\lambda,\mathbf{w},\mathbf{q},X_{N-1}$ used to perform the partial triangulation of $X_{N}$ via (29), with the role of $A^{(N)}$ played now by $X_{N}$. We have

Taking into account the change-of-measure formula (Housholder reflections and partial triangularization of real matrices) we see that the corresponding JPD can be written as

By definition, the density of real eigenvalues $\rho^{(r)}_{N}(\lambda)$ is obtained by integrating the above JPD over variables $|\mathbf{q}|<1,-\infty<\mathbf{w}<\infty$ and finally over $X_{N-1}$. After performing the integrals over $\mathbf{q}$ and $\mathbf{w}$ we immediately arrive at the relation

where ${\cal C}_{N-1}(\tau)$ 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: $\lambda\to x$ and $N\to N+1$.