Antagonistic interactions can stabilise fixed points in heterogeneous linear dynamical systems

We analyse the stability of large, linear dynamical systems of variables that interact through a fully connected random matrix and have inhomogeneous growth rates. We show that in the absence of correlations between the coupling strengths, a system with interactions is always less stable than a system without interactions. Contrarily to the uncorrelated case, interactions that are antagonistic, i.e., characterised by negative correlations, can stabilise linear dynamical systems. In particular, when the strength of the interactions is not too strong, systems with antagonistic interactions are more stable than systems without interactions. These results are obtained with an exact theory for the spectral properties of fully connected random matrices with diagonal disorder.


Introduction
We consider a dynamical system described by n variables x j ∈ that are labeled by indices j = {1, 2, . . ., n} = [n] and where t ∈ + is the time index.The evolution in time of the variables x j (t) is described by a set of randomly coupled, linear differential equations of the form where the A jk are the entries of a random matrix A of dimension n × n.The fixed point ⃗ x = 0 of the set of Eqs.(1) is stable when all eigenvalues of A have negative real parts.On the other hand, if there exists at least one eigenvalue with a positive real part, then the fixed point is unstable.
Differential equations of the form Eq. ( 1) appear in linear stability analyses of complex systems described by nonlinear differential equations of the form ∂ t ⃗ y(t) = ⃗ f ( ⃗ y(t)) where ⃗ y = ( y 1 , y 2 , . . ., y n ).For example, in theoretical ecology ecosystems are modelled with Lotka-Volterra equations, where the variable ⃗ y quantify the population abundances of the different species in the population [1].Other examples are models for neural networks, for which ⃗ y represents the neuronal firing rates or the membrane potentials [2][3][4], and models of economies [5], for which ⃗ y represents economic variables such as the prices of goods.If the differential system determined by ⃗ f admits a fixed point, defined as ⃗ f ( ⃗ y * ) = 0, then the dynamics of ⃗ x = ⃗ y − ⃗ y * near the fixed point is given by Eq. ( 1), where A is the Jacobian of ⃗ f .The linear stability of a complex system that settles in a fixed point state is thus determined by the real part of the leading eigenvalue λ 1 , which is defined as an eigenvalue of the Jacobian matrix A that has the largest real part.
To study complex systems, Wigner [6], Dyson [7] and May [8], among others, suggested to study random matrices A, and the task at hand is then to determine the real part of the leading eigenvalue as a function of the parameters that define the random matrix ensemble.Although one should be careful in drawing conclusions about the dynamics of nonlinear systems from the study of randomly coupled linear differential equations, random matrix theory has the advantage of providing analytical insights about the influence of interactions on linear stability.In fact, linear stability analyses with random matrix theory have been used to study the onset of chaos in random neural networks [2][3][4], the stability of ecosystems modelled by Lotka-Volterra equations with random interactions [8][9][10][11][12][13][14], economies [15], or gene regulatory networks [16].Moreover, although traditionally random matrix models are fully connected, recently exact results have been derived for the stability of linear models defined on complex networks [17][18][19].
So far, stability analyses for randomly coupled, linear dynamical systems have mainly focused on matrices A with diagonal entries, which we also call the growth rates, that are fixed to a constant value d, i.e., [8] where δ j,k is the Kronecker delta function, and where the coupling strengths J jk are random variables drawn from a certain distribution.Following Refs.[9,20], we consider the case where the pairs of random variables (J jk , J k j ) are independent and identically (i.i.d.) distributed random variables drawn from a distribution with 〈J i j 〉 = 0 , 〈J 2 i j 〉 = σ 2 , and where the variance σ 2 of the entries J i j quantifies the strength of the interactions, and τ ∈ [−1, 1] is the Pearson correlation coefficient between the variables J jk and J k j .The sign of the parameter τ is important in theoretical ecology as it determines the nature of the trophic interactions between two species.If the interactions are on average competitive (J i j < 0 and J ji < 0) or mutualistic (J i j > 0 and J ji > 0), then τ > 0. On the other hand, if the interactions are on average antagonistic (J i j > 0 and J ji < 0 or J i j < 0 and J ji > 0), then τ < 0 [9][10][11]19].
In theoretical ecology, antagonistic interactions are also called predator-prey interactions as they describe trophic interactions between two species for which one predates on the other.The leading eigenvalue of random matrices of the form (2) is given by [8,[21][22][23] Re(λ 1 ) = σ(1 It follows from Eq. ( 4) that in the case of homogeneous relaxation rates d < 0 is required for a linear system to be stable.Hence, when the diagonal entries of A are fixed to a constant value d, then interactions J jk always destabilise fixed points in large dynamical system.In the model given by Eq. ( 2) it holds that in the absence of interactions (J i j = 0) either all variables are stable (when d < 0) or all variables are unstable (when d > 0).In this paper, we relax this condition and consider random matrix models with growth rates A j j = D j that fluctuate from one variable to the other.In the symmetric case (τ = 1), such random matrices are called deformed Wigner matrices [24][25][26] and in this case a functional equation that determines the spectral distribution in the limit of large n has been derived by Pastur in Ref. [24].Another case that has been studied in the literature is when A is the adjacency matrix of a random directed graph with diagonal disorder [17,18,27], which corresponds in the dense limit with τ = 0 [28], and for which a simple equation for the boundary of the spectrum as a function of the distribution of diagonal matrix entries has been derived.
On the other hand, in the present paper we focus on the case of heterogeneous relaxation rates D j with negative τ, which is, as discussed, of particular interest for ecology.This case has been studied signficantly less in the literature, a notable exception being Ref. [29].Here, apart from completing the theory of Ref. [29] by deriving analytical results for eigenvalue outliers, which are important when considering the leading eigenvalue, we also show that for negative τ the leading eigenvalue can be negative, even if a finite fraction of the relaxation rates D j are positive.The latter finding, not discussed in Ref. [29], implies that antagonistic interactions can stabilise linear systems when the interactions are neither too weak nor too strong, even if a finite fraction of the variables are unstable in isolation, and this constitutes the main result of this paper.
The paper is organised as follows.In Sec. 2 we define the model that we study, which is a fully connected random matrix with diagonal disorder.In Sec. 3, we discuss the cavity method, which is a method from theoretical physics that we use to study the model in the limit of infinitely large random matrices.In Sec. 4, we present the main results for the boundary for the spectrum of fully connected matrices with diagonal disorder and in Sec. 5 we present analytical results for the eigenvalue outlier.In Sec. 6 we use the obtained theoretical results to derive phase diagrams for the linear stability of fixed points.We end the paper with a discussion in Sec. 7. The paper also contains a few appendices where we present details about the mathematical derivations.

Fully connected random matrices with diagonal disorder
We consider the random matrix model where the (off-diagonal) pairs (J i j , J ji ) are i.i.d.random variables drawn from a joint distribution p J 1 ,J 2 with moments as specified in the Eqs.(3), where the diagonal elements D j are i.i.d.random variables drawn from a distribution p D , and where µ ∈ is a constant shift of the matrix elements.Note that without loss of generality we have set 〈J i j 〉 = 0, as a nonzero average value can be incorporated into the parameter µ.
As will become clear later, just as is the case for the circular law [30,31], in the limit of n ≫ 1 the boundary of the spectrum of A is a deterministic curve in the complex plane that depends on the distribution p J 1 ,J 2 of (J i j , J ji ) only through its first two moments given in Eq. ( 3), and hence we will not need to specify p J 1 ,J 2 .On the other hand, the boundary of the spectrum of A depends in a nontrivial way on the distribution p D , and therefore it will be interesting to study the effect that the shape of p D has on the leading eigenvalue.Due to the constant shift µ, the spectrum may also contain a single (deterministic) eigenvalue outlier in the limit of large n ≫ 1 [17,27,32,33].
In the special case when p D (x) = δ(x − d) we recover the model given by Eq. ( 4).A more interesting case is when the growth rates D j are heterogeneous, and arguably the most simple model for heterogeneous growth rates considers that the D j can take two possible values, yielding a bimodal distribution with d − < 0, d + > 0, and p ∈ [0, 1], and where δ(x − d) denotes the Dirac delta distribution.
In this example, a fraction (1 − p) of variables x j are unstable in the absence of interactions (σ 2 = 0).We also consider cases where p D is a continuous distribution.One example of a continuous distribution is the uniform distribution defined on an interval Since the uniform distribution is supported on a bounded set, we will also consider an example for which p D has unbounded support, namely, we will consider the Gaussian distribution with zero mean and unit variance.
The main question we address in this paper is whether the interaction variables J i j can stabilise a linear dynamical system even when a finite fraction of variables are unstable in the absence of interactions, i.e., a finite fraction of species i ∈ [n] have a positive growth rate D i .In other words, we ask whether it is possible to have Re λ 1 < 0 even when there exists a value d > 0 such that p D (d) > 0.

Cavity method for the empirical spectral distribution of infinitely large matrices
We determine the leading eigenvalue λ 1 in the case of µ = 0, when the spectrum of A has no outliers in the limit of n → ∞.
In this case, the leading eigenvalue of the adjacency matrices A is determined by the empirical spectral distribution ρ of the eigenvalues λ j of A, defined by for all z = x + i y ∈ .The spectral distribution determines the leading eigenvalue of the continuous part of the spectrum through Equation (10) holds as long as the spectrum of A does not have eigenvalue outliers [17,27], which for the model defined in Sec. 2 is the case as long as µ = 0 [17,27].
The convergence in Eq. ( 9) should be understood as weak convergence [31], which implies that the average of any bounded and continuous function f (z) defined on the complex plane converges in the limit of large n to dz ρ(z) f (z).Also, we can drop the average in the righthand side of Eq. ( 9) as the spectral distribution converges almost surely and weakly to ρ [31], and hence also the leading eigenvalue λ 1 as defined in Eq. ( 10) is a deterministic variable for large values of n.
The limiting distribution ρ of random matrix models as defined in Sec. 2 have been studied before in several special cases.Notably, for the symmetric case with τ = 1 Pastur derived a functional equation that determines ρ [24].Recently, the symmetric case was revisited in [26], and in that reference also the large deviations of λ 1 were computed in the case when the matrix entries J i j are drawn from a Gaussian distribution; note that large deviations are not universal and depend on the statistics of (J i j , J ji ) as determined by the distribution p J 1 ,J 2 .In the case when τ = 0 and p D is a bimodal distribution the spectral distribution ρ has been determined in Refs.[34,35] and the τ = 0 case for general p D has been considered in [28].For random directed graphs with a prescribed distribution of indegrees and outdegrees, which corresponds with the case τ = 0 in the limit of large mean degrees, a simple equation was derived for the boundary of the spectrum in Refs.[17,18,27].Lastly, Ref. [29] obtained analytical results for the spectrum when τ < 0 and µ = 0.
We determine the spectral density ρ(z) from the resolvent of the matrix A, which can be determined with the cavity method [36,37].The resolvent is defined as where 1 n is the identity matrix of size n.The spectral distribution can be expressed in terms of the diagonal elements of the resolvent by [34] ρ For non-Hermitian matrices, the eigenvalues are in general complex-valued, and therefore in the limit of n → ∞ we cannot get ρ(z) from TrG(z) [23].To overcome this, we use the Hermitization method [34] that considers the enlarged 2n × 2n matrix where we have introduced a regulator η that keeps all quantities well-defined in the limit of large n, where A T is the transpose of the matrix A, and where z * is the complex conjugate of z.The inverse of the matrix H is where In the limit η → 0, we obtain where 0 n is the matrix with zero entries.Hence, combining Eqs. ( 12) and ( 16), we find that where Tr 21 is a block trace over the diagonal of the lower-left block of H −1 .Defining the jk-th block of the generalized resolvent as the spectral distribution (17) can be written as [37] ρ(x, y) = lim where In Appendix A, we use the cavity method to derive a selfconsistent equation for the matrix g at fixed η in the limit of n ≫ 1, viz., where 〈. . .〉 D denotes the average over the distribution p D .
Note that to derive (21) we have determined g at finite values of η in the limit of large n, and afterwards we take the limit of η → 0. Hence, we interchange the two limits in Eq. ( 19), which is not evident as the leading order, correction terms in Eq. ( 16) at large values of n and small values of η intertwine the two limits.Demonstrating that these two limits can be interchanged constitutes the main challenge in rigorous approaches to non-Hermitian random matrix theory, see e.g.Refs.[31,[38][39][40][41][42].This involves bounding the rate at which the least singular value of z1 n −A converges to zero for large values of n, as the correction terms in ( 16) depend on the inverse of the matrices I l and I r .In this paper, we use the theoretical physics approach, i.e., we exchange the two limits in good faith and then corroborate theoretical results with direct diagonalisation results.In the next section, we use the Eq. ( 21) together with Eq. ( 19) to determine the boundary of the support set of ρ in the complex plane.

Boundary of the spectrum
The support set of ρ(z) is defined as where • denotes the closure of a set.From Eq. ( 10) it follows that the support set determines the leading eigenvalue whenever the spectrum does not contain eigenvalue outliers [27], which for the model defined in Sec. 2 is the case as long as µ = 0.The support set S follows from the solutions to the Eqs.( 19)-( 21).The Eq. ( 21) admits two types of solutions [17,19].First, there is the trivial solution for which g 11 = g 22 = 0 and ∂ z * g 21 = 0, yielding a distribution ρ = 0 for z / ∈ S. Second, there is the nontrivial solution for which g 11 > 0 and g 22 > 0 and ∂ z * g 21 ̸ = 0, yielding the probability distribution ρ > 0 for z ∈ S.
Although the trivial solution solves the set of Eqs.(21) for any value of z and for η = 0, it is only for z / ∈ S that the trivial solution is relevant.Indeed, when z ∈ S the trivial solution is unstable with respect to infinitesimal small perturbations, and hence the regulator η > 0 guarantees that the spectral distribution for z ∈ S is determined by the nontrivial solution.As a consequence, the boundary of the support set S follows from a linear stability analysis of the Eqs.( 21) around the trivial solution [19].Expanding the Eqs.(21) in small values of g 11 > 0 and g 22 > 0, we obtain that for all values of z ∈ S it holds that and the boundary of the support set is given by Note that, in general, Eq. ( 24) is coupled with the Eq. ( 21) and therefore these equations have to be solved together.
In what follows, we first analyse the Eqs.( 21)and ( 24) in two limiting cases, and then we discuss the general case.

Symmetric matrices with
For symmetric random matrices the Eq. ( 21) reduces to a functional equation for the resolvent of a Wigner matrix with diagonal disorder derived originally by Pastur in Ref. [24].Indeed, in this case g 22 = g 11 = 0 for all z with nonzero imaginary part so that for all z / ∈ , which is identical to Equation (1.6) in Ref. [24].Since g 21 is the Stieltjes transform of the spectral distribution defined on the real line, we can use the Sokhotski-Plemelj inversion formula (see e.g.Chapter 8 of [43]) to obtain the spectral distribution.Note that the delta distribution δ( y) on the right hand-side of Eq. ( 26) specifies that the eigenvalues of A are real, and hence the distribution ρ defined on the complex plane equals zero for all values y ̸ = 0.
4.2 Uncorrelated interaction variables J i j and J ji (τ = 0) In the absence of correlations between J i j and J ji , the Eq. ( 24) decouples from the Eq. ( 21).Therefore, the "τ = 0"-case is mathematically simpler to solve than the "τ ̸ = 0"-case.The boundary of the support set S is determined by the values of λ ∈ that solve the equation which is closely related to the results obtained for the boundary of spectra of random directed graphs in Refs.[17,18,27] and to those of perturbed random matrices with uncorrelated matrix entries [28].Equation ( 27) implies that for τ = 0 the leading eigenvalue satisfies In other words, in the absence of correlations between the interaction variables J i j and J ji , interactions always increase the real part of the leading eigenvalue and have thus a destabilising effect on system stability.Let us analyse the boundary of the spectrum and the leading eigenvalue for a couple of examples.As shown in Appendix B, when p D (x) is the uniform distribution supported on the interval [d − , d + ], then the boundary of the support set S is given by values of (x, y) that solve a result that was also obtained in [29].For (d + − d − )/σ 2 ≪ 1, we recover the celebrated circular law [31,44], while for (d + − d − )/σ 2 ≈ 1 the formula Eq. ( 29) expresses a deformed circular law replacing the constant radius σ by y (d In Fig. 1(a) we have plotted the curve Eq. ( 29) for the case d + = 1 and d − = −1 and we show that this theoretical results is well corroborated by the spectrum obtained from numerically diagonalising a matrix.From Eq. ( 29) it follows that the leading eigenvalue is given by Eq. ( 30) reveals that Re(λ 1 ) > d + for any value of σ, and hence the interactions make the system less stable.For d + = d − = d we recover the formula Eq. (4), and in the limit of large As a second example we consider the case when p D is a Gaussian distribution with zero mean and unit variance.In Fig. 1(c), we compare the solution to Eq. ( 27) with the spectrum of a random matrix drawn from the ensemble defined in Sec. 2. In this case the spectrum S in the limit n → ∞ contains the whole real axis, contrarily to the case where p D is a uniform distribution (compare Fig. 1(a) with Fig. 1(c)).The distinction between the two cases follows  32) and (33).
from the fact that p D is supported on a compact interval in the uniform case, while it is supported on the whole real axis in the Gaussian case.Indeed, Eq. ( 27) implies that in the former case the spectrum S is a finite subset of the complex plane, while in the latter case it contains the real axis.Consequently, for a compactly supported distribution p D the leading eigenvalue converges to a finite value as a function of n, while for a distribution p D that is supported on the real axis the leading eigenvalue diverges.The rate of divergence as a function of n of the average of the leading eigenvalue, 〈λ 1 〉, is determined by the scaling of the maximum value of the diagonal entries D i as a function of n.Since the maximum of n i.i.d.random variables drawn from a Gaussian distribution with zero mean and unit variance scales as log(n) (see Theorem 1.5.3 in Ref. [45]), it holds that when p D is Gaussian, where D max = max {D 1 , D 2 , . . ., D n } and where O(•) is the big O notation.

The case of generic correlations between J i j and J ji (τ ∈ [−1, 1])
We consider now the case of nonzero correlations between the interaction variables J i j and J ji .
In this case, it is more difficult to find the values of z that solve the Eq. ( 24), as contrarily to the τ = 0 case Eq. ( 24) is coupled with Eq. ( 21).Nevertheless, we can simplify the Eqs.( 24) and ( 21) by using generic properties of H and A.
Using that H is Hermitian, which is implied by the definition Eq. ( 18), we obtain that g 12 = g * 21 , Im(g 11 ) = 0 and Im(g 22 ) = 0.In addition, since A and A T have the the same statistical properties, we can set g 11 = g 22 .Also, since we are interested in the boundary of the continuous part of the spectrum, which is located at the edge between the trivial and the nontrivial solutions, we can set g 11 = g 22 = 0, as this is satisfied for the trivial solution.Furthermore, we make the ansatz that Im(g 12 ) is independent of the distribution p D , and therefore Im(g 12 ) = y/σ 2 (τ − 1), which is the solution when p D (x) = δ(x).In addition, using that g 11 g 22 = 0, we can express the Eq. ( 21) as and Eq. ( 24) reads We could not simplify these equations further, and hence we will obtain the boundary of the spectrum by solving the Eqs.(32)(33).
In Figs. 2 and 3, we corroborate the boundary of the spectrum, obtained from solving the Eqs.(32)(33), with numerical results for the eigenvalues of matrices of finite size, obtained with numerical diagonalisation routines.We show the boundary of the spectrum for the case of the bimodal distribution p D given by Eq. ( 6). Figure 2 compares two spectra with the same σ but different values of τ, whereas Fig. 3 considers one negative value of τ and observes how the spectrum changes as a function of σ.Note that the the real part of the leading eigenvalue Re(λ 1 ) decreases as a function of τ.
The leading eigenvalue is obtained by solving Eqs. ( 32)-( 33) at y = 0.For bimodal p D we obtain a quartic equation in x and we identify the largest real-valued solution of this quartic equation with Re(λ 1 ).We have obtained an analytical expression for Re(λ 1 ) as a function of the system parameters, which we omit here as it is a long mathematical formula without clear use.However, it can be found in the Supplemental Material of Ref. [29].
For uniform p D the Eqs.( 32)-( 33) can be solved explicitly as shown in Appendix C. Remarkably, in this case we obtain a simple, analytical expression for the leading eigenvalue, viz., where d + > d − ∈ .One readily verifies that for τ = 0 Eq. ( 34) reduces to Eq. ( 30), for τ = 1 Eq.( 34) recovers the result in Ref. [26] for the case of symmetric matrices with entries drawn from a Gaussian distribution, and for σ = 1 it is equivalent to a formula that appeared in the Supplemental Material of [29].Since the sign of the second term of Eq. ( 34) is equal to the sign of τ, the leading eigenvalue λ 1 decreases as a function of negative values of τ.
In Fig. 4 we compare Eq. ( 34) with numerical results of the leading eigenvalue obtained through the direct diagonalisation of matrices of finite size n.The numerics corroborate well the analytical results that are valid for infinitely large n.We make a few interesting observations from Fig. 4: (i) for τ = 0, the leading eigenvalue is a monotonically increasing function of the interaction strength σ implying a continuous increase of the width of the spectrum as a function of σ; (ii) for τ = −0.8, the leading eigenvalue is a nonmonotonic function of σ.Initially, for small values of σ, the width of the spectrum decreases as a function of σ, while for large enough values of σ the width of the spectrum increases as a function of σ; (iii) for τ = −1, the leading eigenvalue is monotonically decreasing.In this case, the width of the spectrum decreases continuously as a function of σ and converges for large σ to a vertical spectrum centered on the mean value of p D .

Eigenvalue outlier
Now, we determine the leading eigenvalue when µ ̸ = 0.Even though, the continuous part of the spectrum is not affected by µ (see Appendix E), the spectrum may have an eigenvalue outlier, which can be the leading eigenvalue; this is illustrated in Panel (b) of Fig. 1.Hence, for µ ̸ = 0 the leading eigenvalue  where λ isol is the eigenvalue outlier, if it exists, and λ c 1 is the leading eigenvalue of the continuous part of the spectrum, as defined by Eq. ( 10) with λ 1 replaced by λ c 1 .In what follows, we determine λ c 1 and λ isol .The leading eigenvalue λ c 1 of the support set S is, in the limit of large n, independent of µ.Indeed, as shown in Appendix E, the boundary of the support set S solves the Eqs.( 21) and (24), just as was the case for µ = 0.
To determine the eigenvalue outlier we follow the theory for eigenvalue outliers of random matrices, as developed in Ref. [17], which is also based on the cavity method, albeit works in a different way as in this approach recursion relations are derived for entries of right eigenvectors, instead of the entries of the resolvent.Following this approach, we show in the Appendix F that the eigenvalue outliers λ isol of A in the limit for large n solve where g 21 is the trivial solution to Eq. ( 21), i.e., for g 11 = g 22 = 0.For τ = 0, which leads to the equation for the outlier λ isol .
In general for τ ̸ = 0, we do not have an explicit expression for g 21 , and hence we express λ isol in terms of the functional inverse f 21 of g 21 , namely, where z = f 21 (g 21 (z)) .( 40) Using Eq. ( 21) for g 22 = g 11 = 0, corresponding to the trivial solution, we obtain the following selfconsistent equation for g 21 (z), If we can rewrite Eq. ( 41) as z = f 21 (g 21 ), then we readily obtain the functional inverse f 21 of g 21 .Now, we determine f 21 for specific distributions of the diagonal disorder D. First, we consider the case when p D is uniform, as in Eq. (7).It then holds for real values of z ∈ that Im(g 21 ) = 0, and from which it follows that and thus Inserting the expression Eq. ( 44) into Eq.( 39), we obtain which is an explicit analytical expression for the outlier when D is uniformly distributed.Equation (45) shows that for negative τ the eigenvalue outlier decreases monotonically as a function of σ, which is different from the nonmonotonic behaviour λ c 1 as a function of σ, see Eq. (34).In Fig. 5, we plot λ 1 = max λ c 1 , λ isol as a function of σ.For small values of σ, the leading eigenvalue λ 1 decreases rapidly as a function of σ, as λ 1 is an outlier, until λ isol = λ c 1 , at which point the outlier stops existing and λ 1 is located at the boundary of S. Also for the case of bimodal p D , as given by Eq. ( 6), we can obtain an explicit expression for λ isol .Following the same steps as for uniform disorder, we get Again, as in the case for uniform disorder, the outlier decreases monotonically as a function of σ when τ < 0. Comparing Eqs. ( 45) with (46), we make the following interesting observation.Both equations take the form where λ (0) isol is the corresponding eigenvalue outlier for τ = 0 solving Since Eq. ( 47) holds for both uniform and bimodal p D , we conjecture that Eq. ( 47) holds for general p D .The Eq. ( 47) is a convenient result as λ (0) isol can be obtained easily from solving Eq. (48).Notice that for constant diagonal matrix entries, i.e., D i = d for all i ∈ [n], Eq. ( 47) is consistent with Theorem 2.4 of Ref. [33] for the eigenvalue outliers of finite rank perturbations of elliptic random matrices.34).Black lines show the maximum between Eq. ( 45), for the eigenvalue outlier, and Eq. ( 34), for the leading eigenvalue of S. Each marker represents the largest eigenvalue of one matrix realisation of size n = 3000.

Stability of linear dynamical systems
We discuss the implications of the spectral results obtained in the previous two sections for the stability of linear systems of the form given by Eq. (1).

Uncorrelated interactions destabilise dynamical systems
For τ = 0 it holds that Re(λ 1 ) ≥ d + for all values of σ [see Eq. ( 28)], which has a couple of interesting implications for the stability of linear dynamical systems.First, a linear dynamical system with τ = 0 cannot be stable if the support of p D covers the positive axis.Second, interactions J i j destabilise linear dynamical systems as λ 1 is an increasing function of σ (see also Fig. 3).Third, if the support of p D covers the whole real line, then the leading eigenvalue λ 1 diverges as a function of n.In the latter case we obtain a tradeoff between diversity, as measured by n, and stability, as measured by Re(λ 1 ) [8,46].Indeed, when p D has unbounded support, then for any realisation of the system parameters σ, τ, and p D , there will exist a value n * so that with large probability Re(λ 1 ) > 0 when n > n * .In Ref. [19], the latter scenario is referred to as size-dependent stability, as the system size n is an important parameter in determining system stability.

Antagonistic interactions can render dynamical systems stable
In the case of negative τ values the interactions J i j can stabilise linear dynamical systems when they are neither too strong nor too weak.To understand how this works, consider linear systems A for which there exist values x ∈ + with p D (x) > 0, such that the system is unstable in the absence of interactions.As illustrated in Fig. 3, adding antagonistic interactions to a linear system can retract the real part Re(λ 1 ) of the leading eigenvalue and make it negative.This example demonstrates that unlike the uncorrelated case with τ = 0, interactions can contribute to the stability of a system when τ < 0. However, as shown in Figs. 4 and 5, for large values of the interaction strength σ the leading eigenvalue increases as a function of σ, and hence antagonistic interactions stabilise linear dynamical systems as long as they are neither too strong nor too weak.Fig. 6 draws the lines of marginal stability, corresponding with Re(λ 1 ) = 0, in the (σ, τ) plane for µ ≤ 0 and for homogeneous growth rates p D (x) = δ(x − D) (dotted line), for a bimodal distributions p D (dashed line), and for a uniform distribution p D (solid line).In these cases, the leading eigenvalue λ 1 is located at the boundary of the support set S, such that, For all cases we have set 〈D〉 = −1, so that we see the effect of fluctuations in D on system stability.Note that for the dotted line d + = −1, whereas for the dashed and solid lines d + = 0.1.As a consequence, for the dotted line a stable region exists when τ = 0 and σ is small enough, while for the other cases there is no stable region when τ = 0. Interestingly, for negative values of τ and for interaction strengths σ that are neither too weak nor too strong, there exists a stable region with Re(λ 1 ) < 0. This region exists even though d + > 0 (solid and dashed lines).On the other hand, for τ = 0 a stable region can only exist when d + < 0, which is the case of the dotted line with homogeneous rates.
Figure 6 shows that Re(λ 1 ) is independent of p D for large values of σ and fixed 〈D〉.We explore this universal behaviour in more depth.Expanding the expression of Re(λ 1 ) for Eq. ( 34) in large values of σ we obtain Identifying the mean and variance of the uniform distribution p D in Eq. ( 49), we can write where 〈〈D 2 〉〉 represents the variance of the diagonal elements.If D is a deterministic variable with zero variance, then we recover the Eq. ( 4).This suggests that when the interactions are strong enough, only the first moment of the diagonal elements is important, rather than the distribution of their elements.Although the relation Eq. ( 50) is derived for the uniform case, numerical evidence shows that it also holds for the bimodal case, and therefore we conjecture that it holds for arbitrary p D distributions.Demonstrating the validity of the Eq. ( 50) beyond the uniform case would be an interesting extension of the present work.Figure 7 draws the lines of marginal stability in the (σ, τ), similar to Fig. 6, albeit with µ > 0. In this case, the leading eigenvalue λ 1 is either the eigenvalue outlier, λ 1 = λ isol , when it exists, or is the leading eigenvalue of the support set S, i.e., λ 1 = λ c 1 .The eigenvalue outlier exists when σ is small enough, and in this regime the system stability increases as a function of σ.For fixed τ, at a value σ = σ * the eigenvalue outlier merges into S, i.e., λ isol = λ c 1 , and for σ > σ * the leading eigenvalue is λ c 1 .We can clearly notice this transition in Fig. 7 due to the cusp that appears in the lines of marginal stability.Hence, for σ > σ * , the lines of marginal stability in Fig. 7 are identical to those in Fig. 6.Comparing Panels (a) and (b) in Fig. 7, we notice that increasing the parameter µ reduces system stability.This can be understood from the fact that the eigenvalue outlier λ isol increases as a function of µ, while the boundary of the continuous spectrum is independent of µ.Note that negative values of µ have no influence on system stability as in this case the eigenvalue outlier, when it exists, is located on the real axis at the negative side of S.

Discussion
We have obtained exact results for the leading eigenvalue of random matrices of the form Eq. ( 5), where the pairs (J i j , J ji ) are i.i.d.random variables drawn from a joint distribution with moments as given in Eq. ( 3), and where the diagonal elements D i are i.i.d.random variables drawn from a distribution p D .If the Pearson correlation coefficient τ = 0, then the boundary of the spectrum solves the Eq. ( 27), which implies that λ 1 ≥ d + [see Eq. ( 28)].Hence, in this case interactions render dynamical systems less stable, irrespective of the form of p J 1 ,J 2 and p D .On the other hand, if the Pearson correlation coefficient τ between the pairs (J i j , J ji ) is negative and the variance of the distribution p D is nonzero, then λ 1 can exhibit a nonmonotonic behaviour as a function of the strength σ of the off-diagonal matrix elements, as illustrated in Figs. 4 and 5.As a consequence, antagonistic interactions that are neither too strong nor too weak can stabilise linear dynamical systems when the diagonal entries D i are heterogeneous, see Fig. 6.
The results in Figs. 4 and 5 can also be understood perturbatively.Indeed, a perturbative expansion of the eigenvalues of A in the parameter σ starting from the diagonal case with σ = 0 leads to the expression (see Appendix D) For the leading eigenvalue j = 1 it holds that D 1 − D i ≥ 0 for all values i ∈ [n], and hence the second term is negative whenever J i j J ji < 0, leading to the initial decrease of the leading eigenvalue λ 1 in Fig. 4. For larger values of σ we need to consider the higher order terms in the perturbative expansion of λ 1 , which are in general positive leading to the nonmonotonic behaviour of λ 1 in Fig. 4. We discuss various interesting questions for future research, both from the mathematical and ecological point of view.
For distributions p D that are uniform and bimodal, we have obtained the analytical expressions Eqs.(46) and (45) for the eigenvalue outlier, which led us to conjecture Eq. ( 47) for general distributions p D .However, we have no proof of this general expression for the eigenvalue outlier, and hence it will be interesting to construct a proof of this simple, generic formula in future work.
For the case of a uniform distribution p D of the diagonal elements we have obtained an analytical expression for λ 1 , which is given either by Eq. (34) or Eq. ( 45), depending on the value of µ, and in the case of τ = 0 we have obtained a closed form expression for the boundary of the support of the spectral density ρ, which is given by Eq. ( 29) and which also holds for uniform p D .The peculiar solvability of the uniform disorder case is consistent with results obtained recently in Ref. [26] for symmetric matrices (τ = 1).Reference [26] shows that when D j = a + b j/n, with a and b arbitrary constants, and when the entries J i j are complex-valued and Gaussian distributed, then an explicit expression for the joint distribution of eigenvalues can be obtained.Based on the results in the present paper, one may speculate that these results are extendable to the case of τ = 0, which will be interesting to explore in future work.
An important extension of the present work are models of the form where C i j is now the adjacency matrix of a random graph.The case of random directed graphs with a prescribed degree distribution p K in ,K out of indegrees and outdegrees has been solved in Refs.[17,18,27].This is the sparse equivalent of the "τ = 0"-case, and in fact the oriented and locally tree-like structure of random directed graphs leads to a decoupling similar to those of Eqs. ( 21) and (24) in the "τ = 0"-case.For this reason, random directed graphs are analytically tractable, and Refs.[17,18,27] derived for the boundary of the spectrum an equation similar to Eq. ( 27), but with a prefactor that is given by (σ 2 + µ 2 ) and analogously, Refs.[17,18,27] derived for the eigenvalue outliers an equation similar to Eq. ( 54), viz., The case of antagonistic interactions (τ < 0) is considerably more challenging as one needs to know the distribution of the diagonal entries [G] ii of the resolvent, which is non trivial in the sparse case, see Ref. [19].Nevertheless, Ref. [19] analysed the antagonistic case without diagonal disorder and found that systems with antagonistic interactions are significantly more stable than systems with mutualistic and competitive interactions (in fact, in the limit n → ∞ they are infinitely more stable).Ref. [19] did however not study the effect of diagonal disorder on system stability.
In an ecological setting, the Jacobian matrix of a set of randomly coupled Lotka-Volterra equations have a specific structure, viz., all elements in a row are multiplied by the population abundance so that The spectra of such random matrix ensembles have been studied in Refs.[48,49], and it would be interesting to study the stabilising effect of antagonistic interactions in this setup.
The question of stability is also relevant for the study of experimental systems, see e.g.Refs.[50][51][52], and the matrices studied in the present paper are null models for real-world systems.However, as discussed in detail in Ref. [53], most ecological data on foodwebs is qualitative, and obtaining quantative data in particular on the Jacobian, is challenging.
Other interesting applications of the theory developed in the present paper are the study of Turing patterns that are governed by randomly coupled chemical reactions, which in Fourier space involves a random matrix with diagonal disorder [54].
Let us end the paper with a word of caution when using the present results to understand the dynamics of nonlinear systems.In a linear system there exist one fixed point, i.e., ⃗ x = 0, and the system parameters will not affect the existence and uniqueness of this fixed point.However, in nonlinear systems this is in general not the case, see e.g.[55,56], and therefore system stability can also be affected by bifurcations that eliminate fixed points.Moreover, for certain problems, such as stability in ecosystems, the fixed point has to be feasible, which for ecosystems implies that all entries of the fixed point are nonnegative [1], and this also constitutes an interesting random matrix theory problem [57].Another issue is that A is the Jacobian matrix, which is in general different from the interaction matrix.Nevertheless, studies in, among others, ecology [12][13][14] and neuroscience [2,3], show that nonlinear systems do exhibit regimes with one unique stationary fixed point and random matrix theory can provide insights on system stability in this regime.In the ecological context for symmetric interactions J i j = J ji , Refs.[12] shows that May's stability argument, albeit in the symmetric setting, applies when the number of extinct species is correctly taking into account, and the corresponding spectrum of the Hessian is described by random matrix theory.Moreover, note that for symmetric interactions replica theory can be used to determine the leading eigenvalue of the Hessian, see Refs.[12,14,58], while for nonsymmetric interactions this is not possible.
When preparing the manuscript, we became aware of the preprint [59] that also studies the spectral properties of matrices of the type defined in Sec. 2. However, the paper [59] discusses the case of τ > 0 for which interactions further destabilise fixed points, whereas we were interested in the potentially stabilising effect of interactions for τ < 0.

A Derivation of the generalised resolvent equation (21) using the Schur complement formula
We derive the Eq. ( 21) using two useful properties.First, we use that permutation of a matrix and matrix inversion are two commutable operations.Indeed, let P be the orthogonal matrix that represents an arbitrary permutation of the integers {1, 2, . . ., n}, then We use this property to perform the permutation [18, 37] where H is the matrix of permuted entries of H.The effect of this permutation is to bundle together the elements of H that depend on pairs of entries (A i , A ji ).Second, we use the Schur inversion formula for the inverse of a 2×2 block matrix [37,60], where s d = a − bd −1 c and s a = d − ca −1 b are the Schur complements of the blocks d and a, respectively.If we choose for a the upper diagonal 2 × 2 block of H, then we obtain where G (1)   jℓ is defined as in Eq. ( 18), but for a matrix A (1) obtained by deleting the 1-th row and 1-th column of A. Permuting the entries of the matrix so that the 1-th row and 1-th column are swapped with the i-th row and i-th column, and using again Eq. ( 56), we obtain the analogous formula Note that the sum in Eq. ( 60) that runs over the indices ℓ and j contains a very large number of terms in the limit of n ≫ 1.Assuming that the law of large numbers applies to this sum -which can be verified to be the case, see Sec.
with u, v ∈ [n] \ {i} and u ̸ = v.In Eq. ( 61) we have used that the pairs (J i j , J ji ) are identically and independently distributed random variables.Moreover, we have used that 〈[G (i) uu ] 12 〉 and 〈[G (i)  uv ] 12 〉 are independent of i, u and v, as in the definition of the random matrix model for A all indices are equivalent.Since 〈J iu 〉 = 0, the second term in Eq. ( 61) equals zero, and using that σ 2 = 〈J 2 ik 〉 and τσ 2 = 〈J ik J ki 〉, we obtain where o(•) is the small o notation.Taking the ensemble average of Eq. ( 62) and, in the limit of large n, identifying for all i ∈ [n], and for all u ∈ [n] \ {i}, we obtain Eq. (21).Equation (64) follows from the fact that A (i) is drawn from the same ensemble as A, except that n → n − 1.
B Derivation of Eq. ( 29) for the boundary of the spectrum when τ = 0 and p D is uniform Equation ( 27) for the uniform distribution Eq. ( 7) gives where we have used z = x + i y.Using the formula for the indefinite integral of 1/(x 2 + a 2 ) when a ̸ = 0, we obtain that for y ̸ = 0 C Derivation of Eq. ( 34) for the leading eigenvalue λ 1 in the case of uniformly distributed diagonal elements We derive Eq. ( 34) for the leading eigenvalue λ 1 when p D is the uniform distribution given by Eq. ( 7).

D Perturbation theory for the leading eigenvalue
We use perturbation theory to understand the functional behaviour in Fig. 4 of Re(λ 1 ) as a function of σ.
Let D + σJ, where D is a diagonal matrix, σ a small parameter, and J an arbitrary σindependent matrix with zero-valued diagonal entries.Let λ For σ = 0, it holds that D 1 = D max and thus the denonimator in Eq. ( 78) is positive.From this it follows that when J i j J ji < 0, λ 1 initially decreases as a function of σ.However, when σ is larger, then the O(σ 3 ) becomes relevant, which provides the nonmonotonic behaviour in Fig. 6.

E Boundary of the spectrum when µ ̸ = 0
We show that the support set S of ρ(z) is, in the limit of large n, independent of µ, and hence in this limit the boundary of S solves the Eqs.( 21) and (24).Following the derivation of Appendix A, we obtain, instead of the self-consistent Eq. (60) for G ii , the self-consistent equation with u, v ∈ [n] \ {i} and u ̸ = v.Using that σ 2 = 〈J 2 ik 〉 and τσ 2 = 〈J ik J ki 〉, Eqs. ( 80) and (80) yield Eq. ( 62), and consequently, in the limit of large n the quantity G ii that determines the resolvent of A is, neglecting subleading order terms in n, independent of µ.

Figure 1 :
Figure 1: Spectra of three random matrices A as defined in Eq. (5) for the uncorrelated case τ = 0 and with diagonal elements that are independently drawn from a uniform distribution [Panel (a) and Panel(b)] or a Gaussian distribution [Panel (c)].Markers denote the eigenvalues of a random matrix of size n = 3000 and with offdiagonal elements A i j = J i j + µ/n, where the J i j are drawn independently from a Gaussian distribution with zero mean and unit variance and where µ is as given in the subfigure captions.The red solid line denotes the solution to Eq. (27), which provides boundary of the support set S in the limit of infinitely large n.The eigenvalue outlier is indicated by an arrow in Panel (b).Panel (a) and (b) show the analytical solution Eq. (29) and Panel (c) is obtained by numerically solving Eq. (27).

Figure 2 :
Figure 2: Comparison between the spectra of two random matrices A with two different values of τ.Eigenvalues plotted are for two matrices of size n = 3000 whose diagonal elements are drawn from the bimodal distribution Eq.(6) with d − = −1, p = 0.9 and d + = 0.1, and whose off-diagonal entries are drawn from a normal distribution with zero mean µ = 0, variance σ 2 /n = 1/n, and τ = 0 [Panel (a)] or τ = −0.7 [Panel (b)].The red line denotes the solution to the Eqs.(32) and(33).

Figure 3 :
Figure 3: Comparison between the spectra of random matrices A with different values of the interaction strength σ.Eigenvalues plotted are for three matrices of size n = 3000 whose off-diagonal elements (J i j , J ji ) are drawn from a joint Gaussian distribution with zero mean, a Pearson correlation coefficient τ = −0.7,and a variance σ 2 /n as indicated.The diagonal elements follow a bimodal distribution with parameters p = 0.9, d − = −1, d + = 0.1, and µ = 0.The red line denotes the solution to the Eqs.(32) and(33).

Figure 4 :
Figure4: Effect of the interaction strength σ on the real part of the leading eigenvalue λ 1 when µ = 0 and for τ = 0 (triangle, dotted), τ = −0.8(circle, solid) and τ = −1 (diamond, dashed).Lines show the Eq.(34).Markers are numerical results obtained for random matrices A with diagonal elements D j that are drawn independently from a uniform p D supported on the interval [d − , d + ] = [−1, 0.1] and with pairs of offdiagonal elements (J i j , J ji ) that are drawn independently from a normal distribution with mean 0, variance σ 2 /n, and Pearson correlation coefficient τ as provided.Each marker represents the largest eigenvalue of one matrix realisation of size n = 7000.

Figure 5 :
Figure 5: Effect of the interaction strength σ on the real part of the leading eigenvalue λ 1 for random matrices with µ = 2 and all other parameters the same as in Fig. 4. Similar to Fig. 4, solid lines/circles correspond with τ = −0.8 and dashed lines/triangles with τ = −1.Gray lines show Eq.(34).Black lines show the maximum between Eq. (45), for the eigenvalue outlier, and Eq.(34), for the leading eigenvalue of S. Each marker represents the largest eigenvalue of one matrix realisation of size n = 3000.

Figure 6 :
Figure 6: Phase diagram for the stability of linear dynamical systems with antagonistic interactions when µ ≤ 0. Lines denote values of (τ, σ) of marginal stability, i.e.Re(λ 1 ) = 0, separating a stable region with Re(λ 1 ) < 0 (below the lines) from an unstable region with Re(λ 1 ) > 0 (above the lines).Results shown are for the random matrix model defined in Sec. 2 and for various distributions p D with the same mean 〈D〉 = −1.The solid line represents a uniform disorder on the interval I = [−2.1,0.1]; the dashed line represents a bimodal disorder with parameters p = 0.5, d − = −2.1,d + = 0.1; and the dotted line represents the case where all diagonal elements take the value −1 with no disorder.

Figure 7 :
Figure 7: Phase diagram for the stability of linear dynamical systems with antagonistic interactions and for µ > 0. Lines show parameter values for which the system is marginally stable (Re[λ 1 ] = 0).The parameters are the same as in Fig. 6, except for µ, which is set to µ = 1 in Panel (a) and µ = 5 in Panel (b).