Local sign stability and its implications for spectra of sparse random graphs and stability of ecosystems

We study the spectral properties of sparse random graphs with different topologies and type of interactions, and their implications on the stability of complex systems, with particular attention to ecosystems. Specifically, we focus on the behaviour of the leading eigenvalue in different type of random matrices (including interaction matrices and Jacobian-like matrices), relevant for the assessment of different types of dynamical stability. By comparing numerical results on Erdős–Rényi and Husimi graphs with sign-antisymmetric interactions or mixed sign patterns, we propose a sufficient criterion, called strong local sign stability, for stability not to be affected by system size, as traditionally implied by the complexity-stability trade-off in conventional models of random matrices. The criterion requires sign-antisymmetric or unidirectional interactions and a local structure of the graph such that the number of cycles of finite length do not increase with the system size. Note that the last requirement is stronger than the classical local tree-like condition, which we associate to the less stringent definition of local sign stability, also defined in the paper. In addition, for strong local sign stable graphs which show stability to linear perturbations irrespectively of system size, we observe that the leading eigenvalue can undergo a transition from being real to acquiring a nonnull imaginary part, which implies a dynamical transition from nonoscillatory to oscillatory linear response to perturbations. Lastly, we ascertain the discontinuous nature of this transition.


Introduction
Understanding the stability of dynamical systems is a fundamental question in various fields of science, ranging from ecology [1,2] and economics [3,4] to neuroscience [5] and chemistry [6,7]. In many cases, the stability analysis of a dynamical system can be reduced to a spectral problem involving a matrix, as discussed in Refs. [8,9]. Therefore, there has been significant interest in understanding how the statistical properties of matrix elements impact the spectral properties of the matrix, which in turn can shed light on the stability of the underlying dynamical system.
As early as the 1970s, using random matrices May has studied, for instance, the stability of fully connected ecosystems [10]. Although this fueled significant interest [11], it is only recently that the influence of sparse network structure on dynamical stability has been studied. Indeed, following pioneering work on the spectra of symmetric Erdős-Rényi graphs [12,13,14,15], recent papers studied the spectra of random, directed graphs [16,17,18,19,20,21,22] and the spectra of random graphs with antagonistic, mutualistic, or competitive interactions [23].
A surprising finding of these more recent works is that the spectra of sparse random graphs are strongly affected by the sign patterns of their matrix entries, i.e., whether sign(A ij A ji ) = 0 (unidirectional interactions), sign(A ij A ji ) = −1 (sign-antisymmetric interactions), or sign(A ij A ji ) = 1 (sign-symmetric interactions). Notably, the spectra of (infinitely large) sparse random graphs are confined to a region in the complex plane with bounded real part when sign(A ij A ji ) ∈ {0, −1} for all pairs i, j of nodes [19,23], whereas the spectra of (infinitely large) sparse random graphs encompass the full real axis when there exists a finite proportion of links with sign(A ij A ji ) = 1 [23].
In the present paper, we provide a simple, sufficient criterion for the finiteness of the real part of the leading eigenvalue, which is the eigenvalue with the largest real part, of an infinitely large, sparse, random matrix. To this aim, we rely on the concept of sign stability.
Sign stability appeared first in studies on qualitative economics [3,4,24] in the 1960s when economists were studying the impact of qualitative properties of interaction matrices on the stability of economic systems, such as, the sign of the elements in the interaction matrices. The importance of sign stability was soon realised for ecology [25,26,27,28] and later it was also considered in chemistry [29]. A matrix is sign stable if for any choice of the absolute values of the non zero elements its leading eigenvalue is also negative, and therefore bounded from above. In order for a matrix to be sign stable, it must satisfy specific constraints on its topology and sign pattern [3,30]. Interestingly, tree graphs admit sign stable structures, for instance, directed tree graphs and antagonistic tree graphs are sign stable. On the other hand, in general, if cycles are present the sign stability property is lost.
As sparse random graphs, e.g., sparse, Erdős-Rényi graphs, contain cycles, they are not sign stable. However, sparse random graphs are locally tree-like [31,32,33], and as tree graphs admit sign stable structures, we say that random graphs are locally sign stable if their interaction patterns match those of sign stable graphs. We then conjecture that the leading eigenvalue of infinite sparse random graphs has finite a real part if the random graphs are locally sign stable. Hence, we argue that local sign stability allows us to predict the stability of large, sparse network structures and hence extends sign stability to sparse random graphs.
One important aspect of sign stability, which also applies to local sign stability, is that it refers to all matrices with the same topology and sign pattern, independently from the absolute value of their nonzero elements. Sign stability is therefore a particularly robust type of stability, as it characterizes an infinitely large set of matrices. This aspect is particularly relevant in ecological applications for at least two reasons. First, most of the time, as discussed in the following section, to determine whether M is stable we need to determine its leading eigenvalue, which requires knowledge of the matrix entries M ij . Unfortunately, in applications it is often the case that only partial information about the matrix M is available. For example, in the context of ecology, it is relatively easy to determine both the foodweb of the trophic interactions between species, i.e., whether M ij = 0 or M ij = 0, and the type of the interactions, inter alia, predator-prey (sign(M ij M ji ) = −1), mutualistic (sign(M ij ) = sign(M ji ) = 1) or competitive (sign(M ij ) = sign(M ji ) = −1) interactions. On the other hand, it is significantly more difficult to determine the strengths |M ij | of the trophic interactions between species [2,34]. This raises the question whether stability can be determined from the sign pattern of the entries of the matrix M. Second, different kinds of ecosystem stability (linear stability, structural stability, feasibility) can be studied by looking at the properties of different matrices obtained from the matrix A of inter-species interactions without sign or topological alterations. In these cases, as explained in more details in the following section, if the topological properties and the sign pattern of the interaction network grant its local sign stability, the ecosystem can be declared at once feasible and stable, both with respect to small fluctuations and small changes in the external conditions.
A second interesting problem for sparse random graphs that we address is whether the leading eigenvalue is real-valued or whether it has a nonzero imaginary part. As it will be recalled, the presence of pairs of conjugate complex leading eigenvalues has important consequences on the dynamical behaviour of the models associated as it gives rise to oscillatory dynamics in the vicinity of the fixed point with frequency of oscillations inversely proportional to the absolute value of the imaginary part of the leading eigenvalue. This aspect is especially relevant for locally sign stable random graphs, as their leading eigenvalue is finite. In particular for these cases, we will discuss how depending on the choice of model's parameter both situations can arise and we will describe the transition between the two corresponding dynamical phases.
The paper is structured as follows: we first review in Sec. 2 the problem of stability in ecology, including a general discussion on the dependence of stability on the system size, and we specify the random matrix models that we study in this paper. In Sec. 3, we review the properties of sign stability with particular attention to its occurrence in tree graphs, before presenting the main results of this paper in the Secs. 4 and 5. In Sec. 4, we identify local sign stability as an important feature for the spectra of sparse random matrices and the stability of complex systems. In particular, in this section we present the paper's main conjecture on sign stability of sparse random graphs, and subsequently, we present several numerical test of the proposed conjecture through direct diagonalisation of the models of sparse random matrices, as defined before in Sec. 2. In Sec. 5, we determine the imaginary part of the leading eigenvalue for several type of sparse, random graphs, and in particular we identify a transition from a regime where the leading eigenvalue is real to a regime where the leading eigenvalue come in a pair of complex eigenvalues. We end the paper with a discussion in Sec. 6, and a few Appendices with technical details.

Stability in Ecology and Model Setup
In this section, we introduce the setup of this paper. In Secs. 2.1 and 2.2, we review the relation between, on one hand, the spectral properties of matrices and, on the other hand, the stability of linear dynamical systems and (nonlinear) ecosystems, respectively. The reader not needing the basic mathematical background or not interested in the ecological applications can skip these two sections. In Sec. 2.3, we review the concepts of absolute and size-dependent stability, which play an important role in this paper. Lastly, in Sec. 2.4, we define the random matrix models that we study in this paper, and in particular in Sec. 2.4.3, we discuss the canonical model parameters that we use.

Stability in linear dynamical systems
Let x(t) ∈ R N , with t ≥ 0 a time index, denote the evolution in time of the state of a system consisting of N components. The simplest model for a dynamical system of N interacting components is given by a linear differential equation of the form where M ∈ R N ×N is an arbitrary matrix. The asymptotic state x ∞ = lim t→∞ x(t) is determined by the eigenvalues λ i (M) of the matrix M. When all the λ i (M) have negative real parts, then x ∞ = 0 [8], and we say that the matrix M is stable. On the other hand, when there exists at least one eigenvalue with a positive real part, then x ∞ does not exist, as the norm of x(t) diverges for large t, and we say that M is unstable. If we order the eigenvalues such that In the intermediate regime for which we say that the matrix is marginally stable. Note that for marginally stable systems x(t) may still diverge as a function of t if the matrix M has degenerate eigenvalues [35,36]. Additionally, the transient dynamical behaviour after external perturbations is revealed by the imaginary part of the leading eigenvalue (λ 1 (M)), where (·) denotes the imaginary part of a complex number. If (λ 1 (M)) = 0, then the transient is nonoscillatory, whereas a nonzero imaginary part implies oscillatory behaviour of x(t) in the vicinity of the origin. The absolute value | (λ 1 (M))| determines the frequency of oscillations of the slowest mode when the system is stable, and of the fastest unstable mode when the system is unstable.
Since in this paper we consider large complex systems, following Ref. [23], we introduce here two variants of linear stability in the limit of large N . Consider a sequence of matrices M N growing in size N ∈ N. In this case, we can distinguish two classes of matrix sequences, viz., those for which the real part of the leading eigenvalue converges to a finite value, i.e., lim and those for which lim Take as an example of the former the nondirected star graph (M N ) ij = δ i,1 (1 − δ j,1 )+δ j,1 (1−δ i,1 ) (with λ 1 = √ N − 1) and as an example of the latter, the directed star graph (M N ) ij = δ i,1 (1 − δ j,1 ) (with λ 1 = 0), where δ a,b is the Kronecker delta function. In the former case, there exists a finite d > 0 such that where 1 N stands for the identity matrix of size N ×N , and hence the sequence M N −d1 N is characterized by linear absolute stability. In the latter case, any constant shift d renders the ensemble stable up to a certain size N , such that for for N > N , and we speak of linear size-dependent stability.
Although linear systems simplify significantly the dynamics of complex systems, they can be insightful for the study of complex systems, such as, ecological systems [11], neural networks [37,38], chemical interaction networks [39], and economic models [40], whenever the interest is to understand the transient dynamics of systems in the vicinity of a stable fixed point by linearizing the system of dynamical equations around the fixed point. We discuss the connection to more general non linear dynamics in more detail in the following section, where other types of stability properties for dynamical systems are reviewed on an example of an ecological model.

Stability in Ecology
The possibility to predict and control the fate of ecosystems is of immediate concern for our lives, which strongly depend upon them. Therefore the concept of their stability in theoretical ecology has been investigated for decades leading to the classification of different types of stability of potential practical relevance.
We give a brief overview of the different notions of stability studied in the literature. Famously modelled by dynamical systems, ranging from simple one-or two-species population evolution [41,42] to more recent studies about multi-species interactions [43], ecosystem's stability with respect to small perturbations around putative fixed points has been investigated at length with linear stability analyses, [10,11], accompanied by considerations on their global or local stability. These approaches assume the existence of at least one equilibrium of ecological relevance, i.e., with nonnegative species' abundances. However it has been pointed out that the conditions for the existence of ecologically meaningful equilibria, called feasibility, are far from trivial [44,45]. Note that in feasible equilibria several species of the original pool can be extinct, and therefore do not enter in the final composition of the ecosystem. Yet, conditions for its (un-)invadability, corresponding with linear stability with respect to potential immigration of species, must be explicitly looked at. A different kind of stability is represented by structural stability, which refers to the sensitivity of (feasible and linearly stable) equilibria (and of their stability and invadability) to changes in the ecological parameters. The present paper discusses in the same general context feasibility, structural and linear stability of model of ecosystems defined by several types of (sparse) interaction graphs.
The well-known complexity-stability trade-off in models of ecological systems generally affects their feasibility [44,45], linear stability [10,46,47,11] and structural stability [48,49,50,51,52]. Specifically, many models of ecosystems assembled from a pool of S interacting species show a dependence of their stability properties on S. We will refer in short to this situation by calling it size-dependent stability. Conversely, we define absolute stability the stability property of a system when it is not affected by its system size. A formal definition of size-dependent stability in linear systems can be found in the previous section. Note that in many cases absolutely stable models can be trivially constructed from models with size-dependent stability by simply rescaling the inter-species interactions by S. However, in absence of direct biological evidences of very weak inter-species interactions, such rescaling can appear as an artefact. Therefore in the following we will call absolutely stable only those models for which stationary points of the dynamics can be made available and stabilized without rescaling the interactions by the system size.
The natural consequence of size-dependent stability is to severely constraint the possibility for these models to account for the emergence of a large biodiversity. We will discuss how ecosystems whose interactions are structured according to sparse, locally tree-like, graphs, in some special cases, can benefit from absolute stability. Such property will affect for different reasons all the different kinds of stability and feasibility, previously reviewed, allowing rich biodiversity to emerge.
To illustrate the mathematical implications of the different concepts of stability and the ways to address questions about them, we refer to a generalised Lotka-Volterra model where the i-th species' abundance N i ∈ R + , with i ∈ {1, 2, . . . , S}, obeys the following dynamical equation where the {α ij } are the entries of the interaction matrix A and where ζ is the immigration rate (to be sent to zero before extracting the results, but useful to avoid considering ecosystems with trivial extinctions, and therefore to grant uninvadability). The other parameters of the model, appearing in the first self-regulation term on the right hand side, are the growth rates in isolation r i and the carrying capacities K i /r i . Depending on the choice of the graph structure, determined by the non zero α ij and of their sign and strength, different type of ecological models can be obtained and studied, from unstructured ecosystems to hierarchical food-webs, from predator-prey (antagonistic) types of interactions to mutualistic or competitive ones. In this family of models, feasibility requires the existence of fixed points of the dynamics, therefore the existence of at least one non trivial meaningful solution N * , with elements {N * i }, to the set of equations f i ({N j }) = − lim ζ→0 ζ/N i . All extinct species i will have N * i = 0 and f i ({N * j }) < 0, while surviving ones will be characterised by N * i = 0 and f i ({N * j }) = 0, therefore The existence of a solution to the last equation is granted by the invertibility of the matrix B with elements restricted to the S * surviving species: For B to be invertible, it is needed that none of its eigenvalue is null or, alternatively, that, in the infinitely large S limit, the continuous part of the spectrum does not include the origin of the complex plane and none of the isolated eigenvalues is null. Feasibility also requires that all elements in N * are non negative. Explicit characterization of the probability to observe a feasible equilibrium are determined for ecosystems on dense unstructured graphs [45], some special ecologically inspired structure of the graph [48,53], and can studied numerically in some more contexts, but no result is known in general. Naturally, the requirement of having nonnegative abundances is more stringent than the condition for the existence of a solution to Eq. (10), yet in the case previously studied the failing of the first condition closely anticipate the breaking of the second [45]. Following this observation and in absence of a general rule able to asses full-fledged feasibility, we will consider the condition for the existence of a non trivial N * as a good proxy for feasibility.
Interestingly, the matrix B is also directly relevant for structural stability, defined as the stability of the abundances of surviving species, N * , to small perturbations of the ecological parameters and consequently small perturbation of The perturbed equations for the N * i then read which can be derived with respect to ξ k revealing that the susceptibility of N * i to little variations of f i is directly determined by the inverse of B: Again a singular behaviour emerges when the spectrum of B contains the origin of the complex plane hinting to a large susceptibility of the solution of N * to ecological parameters. Finally, the classical information on linear stability (stability with respect to dynamical fluctuations as induced by demographic noise, for instance), or Lyapunov stability concerning the domains of attraction of fixed point of the dynamics, is obtained by linearizing the system of dynamical equations around the fixed point N * hence therein evaluating the Jacobian J, a.k.a. the community matrix, with elements Note that contributions to the Jacobian coming from extinct species is diagonal and negative. The non trivial part comes from the S * surviving species and it gives rise to the S * × S * matrix with a non trivial stripy structure where the elements in each row are all rescaled by the same factor N * i . As discussed in the section about linear dynamical systems, linear stability requires that the real part of the leading eigenvalue λ 1 (J) is negative, and a nonzero imaginary part gives rise to oscillatory dynamics in the vicinity of the fixed point with frequency of oscillations inversely proportional to | (λ 1 (J)))|.
The generalised Lotka-Volterra model discussed in this section provides an example of the structure of the matrices of interest when focusing on different facets of the stability of ecological systems. In these structures the interaction matrix A always plays an important role on the elements outside the diagonal, while the self regulation mechanism represented by the carrying capacities contributes to the non trivial diagonal. Moreover, in the Jacobian, each row is multiplied by the abundance of the corresponding species at the fixed point. Note that different examples of single species self-regulation mechanisms contained in the definition of f i , such as those including the so called Allee effect [54] for instance, can lead to less straightforward connections between feasibility, structural and linear stability.
A nowadays widespread approach to model ecosystems with large number of species is to account for the large variety of self-regulation and interaction mechanisms by introducing random parameters, so that the matrices A, B, J may be represented by random matrices. In the simple generalised Lotka-Volterra model considered above, this choice would require to introduce at least one probability distribution p(α) for the amplitudes of inter-species interactions α ij . When referring specifically to feasibility or structural stability determined by B, the simplest setting would imply assuming uniform growth rates across different species and unitary carrying capacities so that B Id = A+1. However, more generally it can be important to include in B the contribution of non trivial diagonal terms of a diagonal matrix D, which are extracted from a second distribution p D (d) to describe the variability of carrying capacities B = A + D. Recall that the important stability trait of B or B Id is whether the spectrum contains or not the origin, which is answered by checking that the smallest real eigenvalue is positive. Under this perspective it is completely equivalent to check whether the largest real eigenvalue of −B or −B Id , or equivalently (as long as the spectra of A is symmetric around the origin) of A − 1 or of A − D, is negative. When focusing on J, instead, the variability of the stationary abundances N * i becomes a more relevant factor in the structure of the matrix. For simplicity, the carrying capacities are then set to the unity, p D (d) represents the distribution of the abundances, placed on the elements of a diagonal matrix D, and we look at J = − DB Id , or equivalently at J =D(A − 1), checking also in this case that the real part of its leading eigenvalue is negative.
Other key ingredients for model selection are the choice of the sign of interactions α ij and α ji and the graph structure of inter-species interactions. The first aspect is related to which type of ecological behaviour determines the inter-species interactions: mutualistic (both positive), competitive (both negative), or predator-prey (of opposite sign). The first two cases will be also called sign-symmetric and the last sign-antisymmetric. The interaction graph structure that can be considered spans from fully connected graphs to several types of sparse graphs. In this work we focus on the influence of short and long cycles on sparse graphs, therefore we will discuss and compare the results obtained on tree graphs (no cycles), Erdős-Rényi graphs (typical cycles of O(log(S * ))), and pure Husimi trees (cycles of fixed, short, length).
All the ecologically motivated characteristics of the random matrices B Id , B and J highlighted and discussed in this section will be encoded in the different types of random matrix models introduced in the next section. In that context the ecological notation is abandoned in favour of a more general random matrix notation where N is the size of the matrix instead of S or S * , and a random diagonal matrix called D can either represent the diagonal matrix of inverse carrying capacities, or contain the elements of the vector of abundances N * .

Absolute stability and size-dependent stability for interaction-like and Jacobian-like ensembles
For all the ensembles of random matrices introduced in the previous section, and defined more generally in the next section, we are interested in behaviour of their spectra when the matrix size is large, i.e., N 1. The stability of the corresponding dynamical system is assured when the spectrum of the associated interaction-like matrix does not include the origin of the complex plane and when the spectrum of the corresponding Jacobian-like matrix has leading eigenvalue with negative real part.
As discussed in the previous sections, we distinguish between models whose stability properties are not affected by their system size N , which we call absolutely stable models, and models whose stability is lost for N larger than a finite size N , and hence their stability is size-dependent.
Absolute stability, in terms of feasibility, linear, and structural stability, is granted when the real part of all eigenvalues of the corresponding relevant matrix is negative for all N , as discussed in Sec. 2.2. In other words, we require that [λ i ] < 0, ∀i ∈ {1, 2, . . . , N } and ∀N ∈ N. Therefore, a necessary condition for absolute stability is that the spectra have a real part bounded from above in the large size limit, which implies that [λ i ] < a, ∀i ∈ {1, 2, . . . , N } and ∀N , where a is a finite constant. In such settings, absolute stability is obtained whenever, for all N , the matrix is equipped with a diagonal that has elements that are smaller than some finite −d < 0, such that all eigenvalues have negative real part.
Note that the stability of densely connected models does not depend on their sign pattern. Hence, dense matrices are either size-dependent, or absolutely stable when their matrix entries are properly rescaled by N , see, e.g., Refs, [11,55,56]. Instead, we focus in this paper on sparsely connected models for which, interestingly, the sign pattern of interactions determines whether the matrix is absolutely stable or exhibits size-dependent stability [23]. Indeed, for large, sparse, random matrices the existence of a finite upper bound for the real part of all eigenvalues may arise in specific settings without the need of any ad-hoc global rescaling by system size, contrarily to the case of dense models.

Random Matrix Models
We define the random matrix models that we study in this paper.
2.4.1. The model structure. We distinguish two type of random matrices defined on sparse random graphs, namely, interaction-like matrices B, and Jacobian-like matrices J. Both matrices B and J are obtained from an interaction matrix A that is the adjacency matrix of a weighted graph, and which specifies the network of interactions between the system constituents.
Interaction-like matrices are the sum of the interaction matrix A and a diagonal matrix D, i.e., where the entries D i ≥ 0 are independent and identically distributed random variables drawn from a distribution p D (d) with d ∈ R + , and the interaction matrix A will be defined in the next subsection. In the special case of D i = 1, for all i, we get what we call the shifted interaction matrix where 1 is the identity matrix. Instead, Jacobian-like matrices J are defined as the product of B Id with D, viz., where the entries D i ≥ 0 are as before independent and identically distributed random variables drawn from a distribution p D (d) with d ∈ R + .

The interaction matrix A. The random matrices A have elements
where C ij ∈ {0, 1} are the entries of the adjacency matrix of a nonweighted, nondirected, random graph, and the α ij ∈ R are the weights of the edges of the graph, which represent the strengths of the interactions. The weights α ij and α ji are pairs of random variables extracted from a probability distribution p α (u, l) that is symmetric under the exchange of its arguments, viz., where and p(x) is a probability distribution supported on R + ; notice that we are particularly interested in models with unbounded support as their norm diverges in the infinite size limit, which is important for the findings in this paper. The constants π S and π O determine the (anti)correlation between the sign of α ij and the sign of α ji (their absolute values are uncorrelated). For π O = 1 elements opposite to the main diagonal of A have opposite signs, i.e., α ij α ji < 0. We speak of sign antisymmetric or antagonistic interactions. When π O = 0, then elements opposite to the main diagonal have the same sign, i.e., α ij α ji > 0. In this case, we speak of sign-symmetric interactions. The elements are positive if π S = 1 and negative if π S = 0, sometimes referred to as mutualistic and competitive interactions, respectively [55]. For intermediate values of π O ∈ (0, 1), we speak of a mixture models, as it contains a mixture of sign-antisymmetric and signsymmetric interactions.
For the adjacency matrix C ij we focus in this paper on two models. One is a random graph model that is locally tree-like, i.e., it has a small number of cycles of small length. The second model is deterministic and has many cycles of small length. In this way, we will be able to address the effect of cycles on our results. The models considered are: • Erdős-Rényi graphs: There are two closely related variants of the Erdős-Rényi (ER) random graph model [57,58]. In the first model, a graph is chosen uniformly at random from the collection of all graphs which have N nodes and M edges. In the second model, the number of nodes N is fixed and each edge connecting two of them exists with a probability q, which is fixed and independent from every other edge. The Erdős-Rényi graphs we use are built according to the second model, that we denote as G(N, q), setting q = c/(N − 1) where c is it is the average number of edges on a single node, also called connectivity of the node. In the limit N → ∞, with c fixed, Erdős-Rényi graphs are locally tree-like graph in the sense that with probability one the finite neighborhood of a randomly selected node is a tree and the typical cycles length grows like log N [31, 32, 33]. • Husimi trees: Husimi trees are connected graphs for which no edge lies on more than one cycle [59]. Loosely said, Husimi trees are trees built out of edges and cycles, such as, triangles, quadrilaterals, pentagons, etc. Husimi trees were introduced by Harary and Uhlenbeck [59], who recognised this graph structure in Husimi's virial expansion of the equation of state of a nonideal gas [60], and whose terminology we adapt in this paper. If Husimi trees are built out of one type of cycle, then one speaks of pure Husimi trees, as in panel c) of Fig.1, while otherwise they are mixed Husimi trees. If the cycles are triangles, then one speaks of a Husimi cactus. When the tree structure is regular with coordination number c, a pure Husimi tree with cycles of length can be defined by using the notation (c, )-pure Husimi tree. For instance, a (c, 1)-pure Husimi tree is a Cayley tree or a Bethe lattice.
As a recap, the models studied across this work can be identified from the matrix structure (interaction like B, shifted interaction B Id , Jacobian like J), the choice of the distributions p and p d , the interaction sign choice (antagonistic, mixture), the graph structure encoded in C (tree, G(N, q) Erdős-Rényi graphs, (c, )-pure Husimi tree).

Canonical model parameters
Here we list the parameters that we use in the numerical results shown in the following sections. Any variations on these will be reported in the figures captions.
As anticipated, we deal with three matrix structures, viz., shifted interaction In numerical examples, we need to specify the adjacency matrix C, the distribution p α of weights (α ij , α ji ), and the distribution p D of diagonal entries D i .
We consider two ensembles of adjacency matrices C, viz., sparse, Erdős-Rényi graphs and Husimi trees. Erdős-Rényi graphs are drawn from the G(N, q) model with q = c/(N − 1). Since we are interested in sparse graphs, the connectivity is kept fixed at c = 2, and thus does not scale with the graph size N . The Husimi trees that we employ are (4, 4)-pure Husimi tree, and hence all cycles have length = 4, as shown in panel (b) of Fig. 2 for the case of a Husimi tree dressed with sign-antisymmetric interactions.
For the probability distribution p of the absolute values |α ij | of the off-diagonal matrix entries α ij , which appears in Eq. (21), we use a truncated Gaussian distribution. In particular, we truncate a Gaussian distribution with mean µ G = 1.0 and variance σ G = 0.6 so that it is supported on the positive part of the real line, viz., where G µ G ,σ G (x) is a Gaussian distribution with mean µ G and variance σ G , erf(x) is the error function, and ϑ(x) is the Heaviside function.
Notice that the first two moments of the truncated Gaussian distribution take the expressions and where · TG denotes the average with respect to the truncated Gaussian distribution. From Eqs. (21), (24), and (25), the first two cumulants of the distribution p α (u, l) of the off-diagonal element pairs follow readily as and where · denotes the average with respect to p α (u, l).
The matrix sign pattern is set by the choices of π O and π S as defined in Eq. (22). In our work we focus on antagonistic models with π O = 1 in which the interactions are sign-antisymmetric and mixture models with π O = 0.9 and π S = 0.5 characterised by a majority of sign-antisymmetric interactions and a smaller portion of sign-symmetric ones.
In numerical examples we are using for p D a uniform distribution supported on [d min , d max ], i.e., where d min and d max are, respectively, the minimum and maximum values of the uniform distribution. Notice that where µ D and σ D are the mean and standard deviation of p D , respectively. We choose µ D > 0 and σ D < µ D / √ 3 such that p D (d) is supported on a subset of the positive real axis and, thereby, all the D i are positive. More specifically, in Sec. 4 we set µ D = 1 and σ D = 0.5, whereas in Sec. 5 we set µ D = 1 and σ D takes values equally spaced between 0 and 0.30.
We diagonalise matrices with linear algebra routines of the Numpy submodule linalg of Python3.

Sign stable matrices
We review sign stability of matrices M ∈ R N ×N , which plays an important role in this paper, as we extend this concept to random graphs in the next section. As anticipated in the Introduction, sign stability refers to all matrices with the same topology and sign pattern, independently from the value of their non zero elements. Therefore, sign stability was introduced in studies on qualitative economics in the 1960s [4,3,24], and found a decade later applications in, among others, qualitative ecology [25,26,27,28], and chemistry [29].
We say that a matrix M is sign stable if any matrix M with the same sign pattern is stable. Note that sign stability is a stronger condition than stability, as it requires that the real part of the eigenvalues of all matrices M in the equivalence class are negative, where sign(x) ∈ {−, 0, +}. Note that matrices M in the equivalence class M(M) can be generated from M through where y ij > 0. Also, note that M(M 1 ) = M(M 2 ), for all M 2 ∈ M(M 1 ). On first sight, sign stability may appear as a too strong condition to be useful. However, as will become soon evident, there exist several, interesting examples of equivalence classes M that are sign stable. Moreover, sign stable matrices have a tree structure, which will make them important for the spectral theory of random graphs that we discuss in the next section. In what follows we discuss necessary and sufficient conditions for the sign stability of the equivalence class M generated by the matrix M.
In the case when M ii < 0 for all i, sufficient and necessary conditions for sign stability have been derived by Maybee and Quirk [3]. These are (see Theorem 3 in [3]): and for all m ≥ 3, Condition (33) implies that edges are either directed or nondirected with antagonistic weights, and condition Eq. (34) states that there are no directed cycles of length m ≥ 3.
In the case when M ii ≤ 0, the conditions Eqs. (33) and (34) are sufficient and necessary conditions for marginal, sign stability. In Refs. [26,30,61] the marginal stable case has been studied in more detail, in particular to determine the conditions for which lim t→∞ x(t) diverges.
We will not discuss the derivation of the conditions Eqs. (33) and (34), as these can be found in detail in Refs. [3,30]. However, we discuss a few notable examples of sign stable matrices, and show explicitly that they are sign stable. In particular, we consider the adjacency matrices of • Oriented tree graphs: these are the adjacency matrices of tree graphs with unidirectional edges, as sketched in Panel (a) of Fig. 1. The unidirectionality of the edges implies and hence Eq. (33), and the absence of cycles in a tree graph implies condition Eq. (34). The sign stability of oriented tree graphs follows readily from the following two facts: (i) all eigenvalues of oriented tree graphs are equal to zero [19,62], which is a direct consequence of the Coefficients Theorem for Directed Graphs, see Appendix A; (ii) the oriented tree property is preserved in the equivalence class M generated by an oriented tree graph. • Husimi trees built out of unidirectional feed-forward cycles: these are Husimi trees built out of motifs that are feedforward cycles, as illustrated in Panel (c) of Fig. 1.
Adjacency matrices of such graphs are sign stable for exactly the same reason as the adjacency matrices of oriented tree graphs: (i) λ j = 0 for all values of j, see Appendix A; (ii) the orientedness and feedforward structure are preserved in the equivalence class M generated by Husimi trees built out of unidirectional feed-forward cycles.
• Nondirected tree graphs with antagonistic interactions: this is a tree graph with predator-prey edges, i.e., M ij > 0 ⇔ M ji < 0 for all i = j, as illustrated in Panel (b) of Fig. 1. The condition (33) ensues from the antagonistic nature of the edges, and the condition (34) from the absence of cycles in a tree graph implies condition. The sign stability of nondirected tree graphs with antagonistic interactions follows from the fact that (i) tree graphs with antagonistic interactions and negative diagonal elements have eigenvalues with negative real part, as we show in Appendix B and Appendix C; (ii) the antagonistic nature of the interactions is preserved in the equivalence class M generated by an antagonistic tree graph.
Note that not all sign stable matrices are tree graphs or Husimi trees. Let's comment on a couple of simple notable examples to illustrate how general is the concept of sign stability. An upper-diagonal matrix with negative entries on the diagonal is a sign stable matrix on any graph structure, even fully connected, as all eigenvalues correspond to the element on the diagonal. On the other hand, antisymmetric matrices for which the entries satisfy M ij = −M ji and M ii = 0, are not sign-(marginally)stable, although having all imaginary eigenvalues, similarly to antagonistic tree graphs with zero diagonal entries (see Appendix B). The reason is that, at variance with antagonistic tree graphs with zero diagonal entries, most matrices in the corresponding equivalence class M of antisymmetric matrices are not marginally stable, because they are not antisymmetric and their eigenvalues can have nonnull (positive and negative) real parts.
A last comment on the extension of sign stable ensembles is that adding negative terms to the diagonal of matrices in a sign stable random matrix ensemble does not affect their stability. This leads to an asymmetry in the spectra of sign stable matrices, as the real part of the eigenvalues is bounded from above, but can extend towards infinity on the negative real axis.

Implications of local sign stability on the spectra of random graphs
In this section, we identify a general criterion, which we call local sign stability, for the finiteness of the average of the real part of the leading eigenvalue of infinitely large, sparse, random graphs, i.e., lim N →∞ [λ 1 ] < ∞; such random matrix models can be made absolutely stable through a constant shift of the diagonal entries, as we discussed in Sec. 2.3.
To introduce local sign stability, we first review the results of Ref. [23], which studied the asymptotic behaviour of the real part of the leading eigenvalue for Erdős-Rényi graphs. This paper finds that for Erdős-Rényi graphs with sign-antisymmetric weights the average, real part of the leading eigenvalue converges to a finite limit as a function of N . This result came as a surprise as the norm of the associated adjacency matrix diverges in the infinite size limit. Moreover, it is sufficient to decorate the Erdős-Rényi graph with a finite fraction of sign-symmetric weights to have a [λ 1 ] that diverges as a function of N , as expected for an Erdős-Rényi ensemble with diverging norm. In what follows, we aim to identify the property that causes the finiteness of the real part of the leading eigenvalue in Erdős-Rényi graphs with sign-antisymmetric weights, and we aim to extend this property so that it can be applied to other random graph ensembles.
To understand the distinction between Erdős-Rényi graphs with sign-antisymmetric and sign-symmetric weights, we build on the locally tree-like property of Erdős-Rényi graphs, see e.g. Refs. [31,32,33]. We say that a random graph is locally tree-like when for large enough N with almost certainty the finite neighbourhood of a randomly selected node is a tree. As discussed in Sec. 3, trees with sign-antisymmetric weights are sign stable, whereas this property does not hold for trees with sign-symmetric weights. However, in general, random graphs have simple cycles and hence, given the condition in Eq.(34), are not sign stable. As a consequence, we can not rely directly on the concept of sign stability introduced in the previous section to understand the distinction between sign-antisymmetric and sign-symmetric Erős-Rényi graphs. This leads us to introduce the weaker condition of local sign stability.
Definition (Local Sign Stability). Let M N be a sequence of adjacency matrices of weighted graphs of N nodes. Let M (i) N (d) be the adjacency matrix of the subgraph generated by a uniformly and randomly selected node i and all of its nodes located within a distance d of i. We say that M N is locally sign stable if the probability that M Since locally sign stable matrices are not sign stable, they can have a leading eigenvalue with a positive real part. Nevertheless, since finite neighbourhoods of large graphs are with almost certainty sign stable, sign stability is broken by either cycles of length ln(N ), which diverge for large enough size N , or by a small number of cycles of finite length. Consequently, we conjecture that the real part of the leading eigenvalue can be positive, but will not grow indefinitely when the matrix size N increases, as the contribution of cycles is not significant enough to make the leading eigenvalue diverge.
We summarise this as follows.
Conjecture (Local sign stability (LSS) conjecture). Let M N represent a sequence of adjacency matrices of random graphs of size N × N , with N ∈ N. Excluding pathological cases (as discussed below), it holds that: In what follows, we make a few comments about the LSS conjecture, before testing its validity on numerical examples.
The LSS conjecture holds (trivially) when the matrix norm M N is bounded. Indeed, the matrix norm is always larger or equal than the real part of the leading eigenvalue, i.e., M N ≥ (λ 1 (M N )) [63]. The interesting cases for the LSS conjecture are those for which lim N →∞ M N = ∞, e.g., for weighted Erdős-Rényi graphs.
The LSS conjecture characterises a sufficient condition, but not a necessary condition to have finite [λ 1 ]. Indeed, for antisymmetric matrices M N it holds that [λ 1 ] = 0, nevertheless they are not locally sign stable.
The LSS conjecture excludes a number of pathological cases, on which we comment now. An example of a pathological case is when the entries M ij with large absolute value, |M ij | 1, concentrate on cycles of finite length, such as, nondirected edges (cycle of length 2) or triangles of nondirected edges (cycle of length 3). For example, suppose that the weights of the nonzero entries of M N are drawn from a distribution with unbounded support, so that the largest value M max = max i,j {M ij : i, j ∈ V } diverges as a function of N . In addition, assume that the two largest matrix entries of the M ij are located on one edge, say M 12 = M max and M 21 = max i,j {M ij : i, j ∈ V } \ {M 12 }. Then in the limit of large N , the leading eigenvalue λ 1 ≈ √ M 12 M 21 , and thus the leading eigenvalue diverges as a function of N , even though M N may be sign stable. Hence, in order for the conjecture to hold, we require that large values of the random weights M ij do not concentrate on (are uncorrelated with) cycles in the graph. Note that in general this will be the case, unless otherwise explicitly imposed, as cycles of finite length are rare in random graphs. Another pathological case is when one node in the graph is incident to a divergent number of small cycles, such as, nondirected edges (cycles of length 2) or triangles of nondirected edges (cycles of length 3). For example, consider the case where N j=1 ϑ(sign(M 1j M j1 ) − 1/2) ∈ O(N ), such that node 1 is incident to a logarithmically increasing number of cycles of length 2 with interactions of the same sign. Then the leading eigenvalue λ 1 ∈ O √ N , even though the graph is locally sign stable. Since we do not have a mathematically rigorous argument for the LSS conjecture, we will rely on known results from the literature on spectra of graphs and numerical tests.
First, let us consider known results on spectra of random graphs. Random directed graphs are locally sign stable, as locally they are trees with unidirectional links, and as shown in Refs. [19,20,21], their leading eigenvalue is finite, which is in correspondence with the LSS conjecture. Another example of a sign stable random graphs are Erdős-Rényi graphs with sign antisymmetric weights, and Ref. [23] shows that the real part of the leading eigenvalue is in this case finite as well. On the other hand, symmetric Erdős-Rényi graphs do not satisfy local sign stability condition Eq. (33) and their spectra contain the whole real axis in the infinitely large limit, see e.g. Refs. [12,13,64,65].
Since examples already present in literature are limited, we now employ new numerical simulations to further test the LSS conjecture. In particular we break local sign stability in two different ways, by removing one of two fundamental ingredients of sign stability: sign-antisymmetric entries and the locally tree-like structure. We compare spectral results for antagonistic (i.e., purely sign-antisymmetric) Erdős-Rényi graphs, which according to our definition are locally sign stable, with two matrix ensembles that are not locally sign stable, namely, mixture Erdős-Rényi graphs (i.e., locally tree-like but with a small fraction of sign-symmetric links) and antagonistic Husimi trees (i.e., keeping the sign-antisymmetric links, while dropping the locally tree-like structure). Hence, the mixture ensemble has a network topology that is favourable, but its interaction pattern breaks local sign stability, while the antagonistic Husimi tree has a favourable interaction pattern, but a network topology that breaks local sign stability.
Since different stability criteria rely on different type of matrices, we consider three cases, as introduced in Sec. 2, namely, shifted interaction matrices B Id in Sec. 4.1, interaction-like matrices B in Sec. 4.2, and Jacobian-like matrices J in Sec. 4.3. Note that the case of shifted interaction matrices as was also studied in Ref. [23], but here we also study the spectra of antagonistic Husimi trees that has not been considered in Ref. [23]. Hence, here we go beyond previous works by breaking local sign stability through both the sign of the weights J ij and through the network topology C ij , and in addition, we go beyond previous works by considering matrices with diagonal disorder and Jacobian-like matrices.

Shifted interaction matrices B Id
First, we test the validity of the LSS conjecture on shifted interaction matrices, B Id = A − 1, as defined in Sec. 2.4. The B Id matrices are the simplest class of matrices that we consider, as the diagonal entries are constant, and hence we can focus on the contribution of the off-diagonal matrix entries to the spectrum. Figure 3 plots the average value of the real part of the leading eigenvalue, (λ 1 (B Id )) , as a function of N for the three ensembles under study, i.e., antagonistic Erdős-Rényi graphs, mixture Erdős-Rényi graphs, and antagonistic Husimi trees. These results, obtained by numerically diagonalising matrices, corroborate the conjecture that local sign stability sets the asymptotic behaviour of the leading eigenvalue. In fact, Fig. 3 shows a qualitative difference in the behaviour of the real part of the leading eigenvalue as a function of the matrix size N : for mixture Erdős-Rényi graphs and antagonistic Husimi trees, [λ 1 ] increases monotonically as a function of N , while for antagonistic Erdős-Rényi graphs, [λ 1 ] quickly converges to a finite value.  So far, we have considered how local sign stability affects the leading eigenvalue of B Id . Instead now, we investigate the effect of local sign stability on the full spectra of matrices B Id , which are plotted in Fig. 4. Figure 4 shows a qualitative difference between, on one hand, the spectra of antagonistic Erdős-Rényi graphs (red), and on the other hand, the spectra of mixture Erdős-Rényi graphs (blue) and antagonistic Husimi trees (orange). Indeed, in the latter two cases the spectrum develops long tails on the real axis, while in the former case the tails are absent.
Analysing how the spectra evolve as a function of the matrix size N , we have found that the tails on the real axis elongate as the matrix size increases, populating larger and larger portions of the real axis (results not shown), which is in agreement with the results on the divergence of the leading eigenvalue in Fig. 3. On the other hand, the antagonistic Erdős-Rényi graph has a spectrum that remains confined in a part of the complex plane that has finite width along the real axis, even when the matrix size increases.
Focusing on the imaginary parts of the spectra, we have found that for all three ensembles under study the spectra grow vertically as a function of N , covering an ever larger portion of the imaginary axis. Notice that the latter result is expected as the matrix norm diverges as a function of N , and M ≥ |λ 1 (M)|. Hence for the antagonistic, Erdős-Rényi ensemble the divergence of the norm materialises exclusively into the growth of the tails of the spectrum parallel to the imaginary axis.
In light of what we outlined in Sec. 2, the observed qualitative difference in the width of the spectrum on the real axis between ensembles that are locally sign stable and those that are not, indicates that local sign stability may be an important characteristic of stability in ecological models. In fact, both structural stability and feasibility require that the origin of the complex plane is not part of the spectrum of interaction-like matrices, and therefore, are not compatible with spectra that exhibit tails covering the whole real axis. On the contrary, the spectrum of antagonistic Erdős-Rényi graphs contains a finite portion of the real axis and, therefore, as explained in Sec.2, the origin of the complex plane can be excluded from the spectrum after a finite shift of the diagonal entries.
Another interesting feature that we observe in Fig. 3 is a, so-called, reentrance effect in the spectrum of antagonistic Erdős-Rényi graphs. The reentrance effect implies that the width of the spectrum is small for eigenvalues with (λ) ≈ 0. Increasing (λ), the width of the spectrum increases until it reaches a maximum at (λ) = (λ 1 ), after which the width of the spectrum decreases again to vanish at large values of (λ). As a consequence, the leading eigenvalue of antagonistic Erdős-Rényi has typically a finite imaginary part, i.e., (λ 1 ) = 0, and hence the leading eigenvalue comes in pairs with its complex conjugate. This reentrance effect, which was already observed in Ref. [23], will be discussed in-depth in Sec. 5.

Interaction-like matrices B
In the previous Section, we have tested the LSS conjecture for interaction matrices with a trivial diagonal. In the present section, we test this conjecture for interaction matrices with fluctuating (negative) diagonal entries, i.e., B = A − D. Therefore, to do not break the condition of local sign stability, the D i s are drawn independently from a distribution p D supported on a subset of the positive real axis. In particular we choose for simplicity a uniform distribution p D (d) supported on [d min , d max ], with d min > 0, even though the main results we obtain on the divergence of the leading eigenvalue also holds for more general distributions p D as long as it is supported on R + . Figure 5 plots (λ 1 (B)) as a function of the matrix size N in the three cases considered before in Fig. 3, albeit now with a distribution p D that has a nonzero variance. The results of Fig. 5 are in correspondence with those of Fig. 3, further establishing the connection between local sign stability and the asymptotic finiteness of the leading eigenvalue. Indeed, for mixture, Erdős-Rényi graphs and antagonistic, Husimi tree graphs the real part of the leading eigenvalue is steadily growing with N , whereas for antagonistic, Erdős-Rényi graphs it converges to a finite value (as a function of N ).
A more detailed look at Fig. 5 reveals that for small values of N < 10 2 , the average leading eigenvalue, (λ 1 ) , of the antagonistic Erdős-Rényi graph increases as a function of N , before it eventually saturates at its asymptotic value for N 10 2 . The transient behaviour of (λ 1 (B)) at small values of N is different from the fast convergence of (λ 1 (B)) in Fig. 3. Moreover, according to Fig. 5 the asymptotic value is approximately equal to −d min , the largest possible value of the diagonal entries; notice that to consider finite size effects, Fig. 5 shows in fact −D min (green stars), with Hence, for the parameters chosen in Fig. 5 the real part of the leading eigenvalue of a large antagonistic Erdős-Rényi graph is determined by the largest possible diagonal element −d min .
Note that −d min does not always determine the leading eigenvalue of antagonistic Erdős-Rényi graphs. For example, let us consider the limiting case of a trivial diagonal with no disorder, i.e., p D (d) = δ(d−1), as discussed in the previous Sec. 4.1. The results of Fig.3 show that (λ 1 ) is larger than −d min = −1, and hence its value is not directly related to d min . Thus, depending on the model parameters, the asymptotic behaviour of the leading eigenvalue of antagonistic Erdős-Rényi graphs is either set to −d min or it is determined by a complex interplay of various parameters. In Sec. 5, we will study the transition between these two regimes in more detail. Nevertheless, we emphasize that in both cases the real part of the leading eigenvalue converges to a finite value, marking a qualitative difference with respect to mixture Erdős-Rényi graphs and antagonistic Husimi trees.   Figure 6 shows the full spectra of the matrices considered in Fig. 5. Comparing the spectra in Fig. 6 with those in Fig. 4, we observe again tails of eigenvalues on the real axis for the mixture Erdős-Rényi and antagonistic (4, 4)-pure Husimi tree ensemble. Note that in Fig. 6 we also observe a tail on the real axis in the spectrum of the antagonistic Erdős-Rényi graph. However, in the latter case, the tail of eigenvalues does not grow indefinitely as a function of N , and instead it is confined to the interval [−d max , −d min ], in agreement with the results in Fig. 5. The tails on the real axis show that also for disordered diagonal entries, local sign stability is related to structural stability and feasibility and consequently local sign stability is expected to be relevant for the stability of ecosystems.
Another difference between Figs. 4 and 6 is that in the latter we do not observe a reentrance effect in the spectrum of the antagonistic Erdős-Rényi graph. We stress however that this is due to the choice of model parameters, and that in general interaction-like matrices can also exhibit reentrance effects. In Sec. 5, we will discuss in detail how reentrance effects appear in models with diagonal disorder, and how they depend on the model parameters.

Jacobian-like matrices J
Lastly, we test the LSS conjecture, i.e., Eq.(36), on Jacobian-like matrices, J = DB Id , that have a distinctive stripy structure, which is relevant for linear stability analysis in ecology as explained in Sec. 2. In order to ascertain that J is locally sign stable, we extract the entries of the diagonal matrix D from a uniform distribution p D (d) supported on [d min , d max ], with d min > 0. Further details on the various parameters can be found in Section 2.4.3. Figure 7 depicts the average, real part of the leading eigenvalue, (λ 1 (J))) , as a function of the matrix size N in the three canonical models of interest, mirroring the analysis in Figs. 3 and 5. The numerical results corroborate the LSS conjecture, viz., (λ 1 (J)) diverges as a function of N for mixture, Erdős-Rényi graphs and antagonistic, Husimi trees, while (λ 1 (J)) rapidly converges to a finite value for antagonistic, Erdős-Rényi graphs. In addition, in agreement with the results in Fig. 5, Fig. 7 shows that for antagonistic Erdős-Rényi graphs the average leading eigenvalue saturates at a value that is approximately equal to −D min (after a transient regime for small values of N ).
In Appendix D we refine the results of Fig. 7 for (λ 1 (J)) as a function of N , by considering the limit d min → 0. In this limit, the antagonistic, Husimi tree exhibits strong transient effects, which we call the Husimi plateau. Nevertheless, in the d min → 0 regime, the results remain in correspondence with Fig. 7 and the LSS conjecture.
The results in Fig. 7 have interesting implications for the linear stability of ecosystems. Recalling the classical linear stability condition (λ 1 (J)) < 0, the results in Fig. 7 imply that system size N is not an important parameter for the stability of systems defined on locally sign stable graphs, while system size N is an important parameter in the general case of models defined on graphs with sign-symmetric interactions or with short cycles, as stability is only attained for small enough values of N .   Figure 8 shows the spectra of Jacobian-like matrices for the three canonical models under study. The Jacobian-like spectra have an arrow-like shape, which resembles those already observed for the dense version of Jacobian-like matrices [66,67]. However, importantly in the sparse case a clear distinction is observed between, on one hand, antagonistic Erdős-Rényi graphs, and on the other hand, mixture Erdős-Rényi graphs and antagonistic Husimi trees. The latter two exhibit long tails on the real axis that increase with system size, while the former does not exhibit such tails. Hence again, the divergence of the real part of the leading eigenvalue is due to tails that develop on the real axis of the spectra of sparse random graphs, and such tails are absent in locally sign stable ensembles.
Note that the spectrum of the antagonistic, Erdős-Rényi graph in the left panel of Fig. 8 does not exhibit a reentrance effect, similar to the spectrum of the interactionlike matrix in the left panel of Fig. 6, but different from the spectrum of the shifted interaction matrix in the left panel of Fig. 4. We stress that this is due to the choice of model parameters, and in fact Jacobian-like matrices can also exhibit reentrance effects. In the next section we will investigate how reentrance in the spectra of Jacobian-like matrices is governed by an interplay between diagonal disorder and network structure.

Discontinuous transition in the imaginary part of the leading eigenvalue
As shown in the left Panel of Fig. 4, the boundary of the spectrum of antagonistic Erdős-Rényi graphs exhibits a reentrance in correspondence of the real axis. In this case, the leading eigenvalue comes in a pair of complex conjugate values with finite imaginary part. From a dynamical systems point of view, the reentrance effect is interesting, as the imaginary part of the leading eigenvalue determines the frequency of oscillations of the slowest mode of relaxation towards the fixed point. Hence, when the leading eigenvalue has a nonzero imaginary part, then the leading, relaxation mode is oscillatory, while for leading eigenvalues that are real the leading mode is nonoscillatory. As a consequence, a transition from a phase in which the leading eigenvalue comes in a pair of two conjugate, complex values to a phase in which it is real corresponds, from a dynamical perspective, to a transition from an oscillatory to a nonoscillatory relaxation dynamics. For this reason we call it a dynamical transition.
Note that in general reentrance effects can be interesting both for interaction-like matrices and Jacobian-like matrices. In the latter, they identify oscillatory dynamics in nonlinear systems of the generalised Lotka-Volterra type in the vicinity of a fixed point. In the former, they may simply identify oscillatory dynamics of corresponding linear systems. However, for the generalised Lotka-Volterra model discussed in this paper, the dynamical behaviour is partially accessible by studying the Jacobian-like matrix only, and the reentrance effect visible in the spectrum of interaction-like matrix does not have any dynamical consequences. On the other hand, the emergence of a reentrance in the spectrum of the interaction-like matrix for the generalised Lotka-Volterra model can still represent an interesting piece of information in situations where the system is feasible/structurally stable, because the spectrum does not include the origin, although (complex) eigenvalues happen to have real parts of both positive and negative sign. In such cases the sign of the leading eigenvalue cannot directly determine structural stability or feasibility, as otherwise naively expected.
Antagonistic Erdős-Rényi graphs with constant diagonal entries, i.e., shifted interaction matrices, exhibit a dynamical transition as a function of the mean degree c, as shown in Ref. [23]. Indeed, for large c the boundary of the spectrum resembles the elliptic law, and thereby the leading eigenvalue is typically real in the large size limit. On the other hand, for small c the boundary of the spectrum shows a reentrance effect, as shown in the left Panel of Fig. 4, and the leading eigenvalue is typically complex. In addition, Ref. [23] shows that the dynamical transition is a continuous transition in the following sense: the imaginary part (λ 1 ) of the typical value λ 1 of the leading eigenvalue, for which we provide a mathematical definition later, equals zero at the transition, and therefore the frequency of the oscillations near the transition is small.
In the present Section, we also study the dynamical transition in the leading eigenvalue, but now for the interaction-like and Jacobian-like matrices with signantisymmetric weights. Note that in Figs. 6 and 8 we have not observed reentrance effects for interaction-like and Jacobian-like matrices, and consequently their leading eigenvalue is real. Instead, in Fig. 4 we have observed a reentrance in the spectra of shifted interaction matrices with sign antisymmetric weights, and consequently in this case the leading eigenvalue comes in a pair of two conjugate, nonreal eigenvalues. Since the latter random matrix ensemble can be obtained from the former ensembles in the limit of small σ D , a transition, possibly related to the reentrance of the spectra in correspondence of the real axis, is to be expected at intermediate values of σ D .
Interestingly, the dynamical transition we find in this section for interaction-like and Jacobian-like matrices occurs while the spectrum is still reentrant in correspondence of the real axis, and hence it features a discontinuous jump in the imaginary part (λ 1 ) of the leading eigenvalue, at variance with the nature of the transition studied in Ref. [23]. This Section is devoted to a careful study of this transition. In particular, we further discuss the definition of a suitable control parameter for the transition, both in the interaction-like and Jacobian-like matrices, involving σ D and other relevant parameters of the model. For both cases, in Sec. 5.1 we locate the transition point through a finite size scaling analysis and in Sec. 5.2, by some additional finite size scaling studies reported in Appendix E, we also show that the transition takes place with a discontinuous jump.

Locating the transition point
As anticipated by the results of the previous section, a change in the strength of the diagonal disorder, as quantified by σ D , can have a direct impact on the imaginary part of the leading eigenvalue of interaction-like and Jacobian-like matrices, determining a qualitative change in the relaxation dynamics of a corresponding dynamical system.
Concerning other model's parameter, in the interaction-like case, µ D is only responsible of a global shift of the spectrum and cannot affect the imaginary part of the leading eigenvalue. A global rescaling of all matrix elements also affects trivially the spectrum. First non trivial changes in the spectrum emerge when σ D changes with respect to the scale of the off-diagonal elements, represented by their variance σ, and therefore the relevant control parameter for the dynamical transition must be σ D /σ ‡.
For Jacobian-like matrices the situation is more involved as their off-diagonal elements are the result of the product of pairs of random variables D i , A ij , with probability distribution p D and p. In this case, a global rescaling of the D i gives also a trivial global rescaling of all the elements of the matrix. Any other modification of the parameters of p D and p induces a non trivial modification of the probability distribution of the off-diagonal elements, which cannot be exactly recast in terms of the variation of a simple control parameter. However, it is possible to derive a rough estimate of the scale of the off-diagonal elements via their variance σ O = σ µ 2 D + σ 2 D , which suggests that the relevant control parameter is approximately given by the ratio We analyse in Fig. 9 the probability P (λ 1 ∈ R) that the leading eigenvalue is real as a function of the control parameters defined respectively for interaction-like and Jacobian-like matrices of antagonistic Erdős-Rényi graphs. For small values of s B and s J , P (λ 1 ∈ R) ≈ 0.1, and hence with high probability the leading eigenvalue has a nonzero imaginary part. On the other hand, for large values of s B and s J , P (λ 1 ∈ R) ≈ 1, and hence the leading eigenvalue is typically real. Because of the aforementioned dynamical significance of (λ 1 ), we call the former the oscillatory phase and the latter the nonoscillatory phase. ‡ Additional non trivial modifications of the spectra are introduced by changes in the relative importance of µ G and σ G for fixed σ, as well as in general every change in the distribution of the diagonal and off-diagonal elements. Conversely, here we consider cases in which changes in σ are only obtained by changing both µ G and σ G , with µ G /σ G fixed. To show that the transition from an oscillatory to a nonoscillatory phase is a proper phase transition, we also plot in Fig. 9 the probability P (λ 1 ∈ R) for different system sizes N . Notably, the transition becomes sharper as the system size increases, and the curves for different values of N intersect at one single point, which we denote by s d , indicating where the dynamical transition takes place. We also show results obtained at different values of σ, σ D , µ D , which all collapse when plotted as function of s B and s J . In particular, this evidence confirms that s J can be effectively used as relevant control parameter for the transition §.
Note that in the oscillatory phase P (λ 1 ∈ R) ≈ 0.1 for large values of N , and hence there is a small nonzero probability that the leading eigenvalue is real. Hence, the imaginary part of the leading eigenvalue is not a self-averaging quantity, which is in correspondence with the numerical results from Ref. [23]. § Be warned that for Jacobian-like matrix, s J can be safely considered the relevant control parameter only as long as the modifications in the model parameters σ, µ D , σ D do not give rise to a significant change in the shape of the distribution of the off-diagonal elements when the transition takes place, as it is the case for the two series of data for the largest system size in the right panel of Fig.. 9.

Characterising the discontinuity of the transition
In this Section, we determine whether the dynamical transition occurs with a continuous or discontinuous variation of the imaginary part of the leading eigenvalue. Figure 10 displays the distribution p (λ 1 ) (x) of the imaginary parts of the leading eigenvalue obtained by generating a large sample of interaction-like matrices B (left Panel) and Jacobian-like matrices J (right Panel); notice that for each pair of conjugate leading eigenvalues with nonzero imaginary part we only focus on the one with (λ 1 ) > 0. From Fig. 10, we observe that the distribution p (λ 1 ) (x) consists of two parts, viz., , and the dashed vertical lines denote the mode of these fitted distributions, which provide the estimates for (λ 1 ). The delta peaks in zero denote the probability P (λ 1 ∈ R) that an eigenvalue is real. a delta distribution at zero carrying the fraction P (λ 1 ∈ R) of matrix realisations that have a real-valued leading eigenvalue, and a continuous distribution corresponding with eigenvalues that have nonzero imaginary part: We define (λ 1 ) := argmax x∈R +p (λ 1 ) (x), as the typical value of the imaginary part of the leading eigenvalue.
As the control parameter changes from s < s d , to s ≈ s d and to s > s d (obtained for σ D = 0.05, σ D = 0.10 and σ D = 0.15 in correspondence of σ G = 0.6 and µ G = 1.0) it can be observed that the maximum of the continuous part of the distribution, indicated by the vertical dashed line, is mostly independent on the control parameter both for interaction-like and Jacobian-like matrices. The change in the control parameters mainly affects the weight P (λ 1 ∈ R) carried by the two parts of the distribution p (λ 1 ) (x), as it is usual with discontinuous phase transitions and at variance with the continuous dynamical transition found in Ref. [23] for shifted interaction matrices as a function of c. Since the numerical results in Fig. 10 are obtained at finite N , we perform in Appendix E a detailed, finite size scaling analysis of the distribution p (λ 1 ) (x), which confirms that the discontinuity in (λ 1 ) at the transition is to be expected also in the large N limit.
Following the approach implemented in [23], to get an estimate of the typical value, (λ 1 ), of the portion of the distribution corresponding no nonzero imaginary part we take the mode of a γ-distribution γ(x; α, β) = β α x α−1 e −βx Γ(α) fitted on the histogram of (λ 1 ) > 0. Here Γ(α) is the gamma function with parameters α, β ∈ R + real and positive.
The typical value (λ 1 ) for different sizes N is plotted in Fig. 11 as a function of the control parameter for the interaction-like and Jacobian-like case and with the same settings as in Fig. 9. As it can be observed again, (λ 1 ) is an almost constant function of the control parameter and its value at the transition point is positive until it vanishes at s d , showing the discontinuous nature of the transition and the associated finite size effects. For further evidence, the lighter markers show that (λ 1 ) is nonzero also in the nonoscillatory phase (the dynamical transition point is indicated by the vertical dotted line) when the leading eigenvalue is typically real. Finally, notice that while the location of the transition is only controlled by s B or s J , respectively, the value of (λ 1 ) depends on the particular choice of the model parameters. This value reflects the global rescaling of the matrix as discussed at the beginning of Sec. 5.1.
Lastly, to develop a better understanding about the mechanism of the discontinuity in the dynamical transition of interaction-like and Jacobian-like matrices, we investigate their spectra in Figs. 12 and 13 at different values of the control parameters. Let us first discuss the spectra of antagonistic, interaction-like matrices in Fig. 12. We observe that the spectrum contains two parts, viz., a cloud of eigenvalues that have a nonzero imaginary part and a segment of real-valued eigenvalues. In particular, the width of the segment of real eigenvalues is approximately equal to [−D max , −D min ], and hence the segment width increases as a function of σ D (or as a function of s B , for fixed σ). The discontinuous nature observed for the dynamical transition is originated from the competition between the width of the cloud of complex eigenvalues and the width of the segment on the real axis. It can be observed that for s B < s d , the leading eigenvalue belongs to the cloud of complex eigenvalues, and since the shape of this cloud is reentrant the imaginary part of the leading eigenvalue in this regime is nonzero. On the other hand, for s B > s d , the leading eigenvalue belongs to the segment of real eigenvalues, The values of (λ 1 ) are estimated as shown in Fig. 10, viz., as the mode of a γ-distribution fitted to a histogram built out of 1000 matrix realisations. The error bars are computed by standard error propagation from the fitting parameters. The parameters employed for empty and filled markers are the same as those used in Fig. 9, respectively. The vertical dotted line depicts the transition value s d . Note that the markers after the transition are denoted in a lighter shade because they represent data which are merely the result of finite size effects: in the nonoscillatory phase the number of complex leading eigenvalues decreases with system size. After the transition the eigenvalues are instead typically real, as represented by the solid red line ( (λ 1 ) = 0). Its left dashed continuation shows rare real leading eigenvalues also present in the oscillatory phase. and hence has null imaginary part. Since the eigenvalue cloud is still reentrant at the point s B = s d when the segment width overtakes the width of the cloud, the transition is discontinuous. This behaviour is different from the continuous transition driven by the connectivity c at s B = 0, as discussed in Ref. [23]. In this case, the transition in the imaginary part of the leading eigenvalue is the result of a gradual reshaping of the spectrum resulting in the progressive disappearance of the reentrance in correspondence of the real axis at large c.
For the spectra of Jacobian matrices, shown in Fig. 13, the qualitative picture is similar to what we have discussed for Fig. 12 for interaction-like matrices, viz., the spectrum consists of two parts, one being a cloud of complex eigenvalues reentrant in correspondence of the real axis and the other is a segment of real eigenvalues approximately supported on [−D max , D min ]. Again, a discontinuous transition on the imaginary part of the leading eigenvalue takes place as the cloud of eigenvalue is reentrant when the width of the segment overtakes the width of the cloud of eigenvalues.
In both cases, in the nonoscillatory phase, (λ 1 ) ≈ −D min , as the support of the real eigenvalues is well approximated by [−D max , −D min ], which is consistent with the results in Figs. 5 and 7. Moreover, in both cases, it emerges that the reentrance of the cloud of complex eigenvalue in correspondence of the real axis also becomes slightly less pronounced when σ D increases, which hints to the possibility that, for specific settings, the disappearance of the reentrance could take place before the largest real eigenvalues become the leading one. In such situation a continuous transition of the imaginary part of the leading eigenvalue is to be expected, instead of the discontinuous transition observed here. In  characteristic arrow shape, which has also been observed before for dense matrices, see Ref. [66,67]. This different behaviour can produce a qualitative difference in the phase diagram containing discontinuous and continuous dynamical transitions of interactionlike and Jacobian-like matrices, which will be also interesting to study.

Discussion
In this paper, we have focused on two interesting features of the spectra of sparse random graphs. First, we have analyzed how the leading eigenvalue of random graphs depends on the system size, taking into account the sign pattern of matrix entries and the graph topology. Second, we have studied how the imaginary part of the leading eigenvalue transitions from zero to a nonzero value as a function of the model parameters in the large system limit. These features are specific to sparse network topologies and do not appear in highly connected graphs.
More specifically, this paper presents a simple, general criterion for predicting the asymptotic behavior of the leading eigenvalue of infinitely large, sparse, random graphs, based on the concept of local sign stability. Although examples of the importance of the sign pattern and network topology on the stability of some matrices ensembles have been highlighted previously in specific cases (Refs. [19,23]), we extend the result by conjecturing that infinitely large graphs that are locally sign stable have a leading eigenvalue with finite real part. We also check our conjecture with numerical tests on sparse random graphs with different topologies and sign patterns, and for different types of matrices such as adjacency matrices and Jacobian-like matrices.
Aside the numerical evidence, we can also provide an intuitive explanation for the relation between local sign stability and the finiteness of the real part of the leading eigenvalue within the context of linear, dynamical system defined on a locally sign stable graph. When a perturbation is applied to a random node its neighbours react. However the finite neighbourhood is, by definition, stable, therefore the perturbation will be damped. This scenario could, in principle, be overturned by cycles of length O(log(N )) that break sign stability, but since they are long we do not expect their effect to be strong enough to destabilise the system in all cases. The situation is different for large, random graphs that are not locally sign stable. In such graphs there exist nodes that are not locally stable, and hence perturbations applied to these nodes will always locally grow and destabilize the system.
One outstanding feature of local sign stability is that it is not affected by the absolute values of the weights, and thus it refers to the set of graphs with the same network topology and sign pattern. As a consequence, whenever the interaction matrix A of an ecosystem is locally sign stable, then also the corresponding interaction-like B and Jacobian-like J matrices are locally sign stable, potentially ensuring at once feasibility, structural stability and linear stability. This aspect is particularly interesting in the ecological context where, as anticipated in the Introduction, it is difficult to quantify the strength of interactions between species [2,34]. In this regard, Ref. [28] has shown that for a general predator-prey model defined on a tree, the Jacobian matrix evaluated on an isolated, feasible equilibrium point has eigenvalues with negative real parts when the Jacobian matrix has negative diagonal elements, and hence the equilibrium point is linearly stable for such systems. The concept of local sign stability that we have introduced in the present paper generalises the results in Ref. [28], as it applies to a broader class of systems, including predator-prey ecosystems defined on random graphs which are locally tree-like.
Our results challenge the classical notion that complex, dynamical systems exhibit a trade-off between size and stability [10,11,34,46,71]. Instead, we demonstrate that this trade-off depends on the sign pattern and topology of the graph. More specifically, Ref. [23] showed that the real part of the leading eigenvalue of infinitely large, antagonistic Erdős-Rényi graphs is finite. The present work generalises this result by extending it to various matrix structures relevant in ecology and by identifying the common feature underlying the finiteness of the real part of the leading eigenvalue.
Additionally, note that the stability properties of sparse, locally tree-like graphs are much more sensitive to the sign pattern, and therefore to the type of the interactions, than dense models. Indeed, in dense graphs the sign pattern does not alter how the leading eigenvalue scales with system size [11,55,56]. Systematic studies of the stability properties of sparse systems are far from being achieved, although potentially relevant for many applications of dynamical systems. In particular new interest is arising on the rich phenomenology of sparse ecosystems [72]. We believe that our paper gives an important contribution by grasping a general criterion for stability of sparse ecosystems.
Concerning applications, our findings on the importance of local sign stability in enhancing stability are consistent with empirical observations on real food webs, which are graphs that represent predator-prey interactions in ecological systems. According to our results, locally tree-like structures and antagonistic interactions stabilize large ecosystems, and hence, these are the structures we expect to observe. Empirical studies have shown that food webs are indeed locally tree-like, with a number of cycles that is similar to those found in locally tree-like Erdős-Rényi graphs [73]. Notably, other networks such as social and technological networks have a significantly larger number of cycles, a feature unique to food webs [73]. In addition, Ref. [74] found that the weights of long cycles in real food webs are typically smaller than in random matrices, further underlying the importance of locally tree-like structures for large ecosystems.
The second main result presented in this paper concerns the imaginary part of the leading eigenvalue (λ 1 ) of antagonistic Erdős-Rényi graphs. The spectra of antagonistic Erdős-Rényi interaction-like and Jacobian-like matrices display a transition driven by the strength of the diagonal disorder σ D from a low σ D -phase in which the leading eigenvalue typically has nonzero imaginary part to a large σ D -phase in which it is real. In the first phase the spectrum is characterised by a reentrance effect around the real axis and this is the cause for it to have a pair of complex conjugate leading eigenvalues.
This reentrance effect is specific to sparse, low connectivity graphs with signantisymmetric weights. In fact, as soon as we move away from the low connectivity limit the spectra tend, under fairly general conditions, to have an elliptic shape. It is worth noting that the spectral reentrance observed in sparse random graphs is a qualitatively new feature, which has only been previously observed in the recent work [23]. That paper identifies a transition in the spectra of antagonistic Erdős-Rényi graphs with zero diagonal entries driven by the connectivity c, from a large c phase where (λ 1 ) = 0 to a small c phase where typically (λ 1 ) = 0, and which also displayed a reentrance effect. The present paper generalizes this result by extending it to more elaborate matrix structures. In particular, we have shown that the reentrance effect persists even in the presence of a disordered diagonal or a stripy, Jacobian-like, structure, but only as long as their disorder is not too large.
It is important to stress once more that while the transition described in [23] is continuous, the one discussed in the present paper is discontinuous. In the former case, the reentrance disappears progressively as the connectivity increases. On the other hand, in the latter case, the spectrum is composed of a cloud of complex eigenvalues and of a segment lying on the real axis. In the low σ D phase, the leading eigenvalue is determined by the cloud and is therefore complex. As σ D increases, both of these components gradually change: the cloud reshapes, reducing the reentrance effect, while the real segment elongates. At the transition, the right tip of the real segment reaches the convex hull of the cloud, and the leading eigenvalue starts to be real. This overtaking occurs before the reentrance has completely disappeared, and therefore we observe a jump in the leading eigenvalue imaginary part, giving rise to a discontinuous transition. At the same time it is possible that, for specific settings, the reentrance effect could disappear completely before (or together with) the overtaking by the segment on the real axis thus making the transition continuous. A more detailed study of the phase diagram describing both the reentrance effect and the transition in (λ 1 ) is left for future work. Moreover, since this reentrance effect is peculiar to sparse graphs, we expect the network topology to be relevant and therefore it would be interesting to investigate its impact.
The transition in (λ 1 ) discussed in this paper, especially with regard to Jacobianlike matrices, is of interest in the context of dynamical systems. In this framework, the imaginary part of the Jacobian leading eigenvalue determines the oscillation frequency of the slowest mode of relaxation towards the related fixed point. Accordingly, the response to a perturbation around a fixed point is oscillatory when the leading eigenvalue has imaginary part different from zero, nonoscillatory otherwise. Our result on the transition indicates that the dynamical response of a nonlinear system defined on antagonistic Erdős-Rényi is oscillatory only if the strength of the diagonal disorder σ D is small compared to the off-diagonal one σ O .
In order to appreciate the actual implication of this dynamical transition it would be interesting to derive a phase diagram describing the dynamical behaviors of antagonistic systems defined on sparse graphs. Indeed, in the context of dynamical systems with nonsymmetrical interactions the spectra of Jacobian matrices can play a role only in the case in which the dynamics is attracted by fixed points while chaos and limit cycles are not predominant. In dense ecological models such a phase diagram has been derived and it identifies regions with multiple attractors where solutions are often chaotic [75,76,77,78]. An analogous phase diagram for sparse ecosystems with predator-prey interactions is not known and its derivation is an interesting endeavour left for future work.
Appendix A. All the eigenvalues of the adjacency matrices of weighted, oriented graphs without directed cycles are equal to zero Let M ∈ R N ×N represent the adjacency matrix of a weighted, directed graph. We assume that M ii = 0, so that there are no self-links. If M ij = 0, then the directed edge from i to j is absent, while if M ij ∈ R \ 0, then there exists a directed edge from i to j weighted by the value of M ij .
We say that a directed graph is oriented if all of its edges are unidirectional, i.e., for all pairs of indices i, j. Additionally, we say that a graph has no directed cycles when Eq. (34) holds. Note that cycles that are nondirected, as for example the feedforward cycles in Panel (c) of Fig. 1, are allowed. The characteristic polynomial of an adjacency matrix of an oriented graph without directed cycles is given by and consequently all the eigenvalues of M are equal to zero, i.e., for all i ∈ {1, 2, . . . , N }. The Eq. (A.2) is a specific case of the so-called Coefficients Theorem for directed graphs, which we revisit here. First, we present the Coefficients Theorem for unweighted graphs, i.e., for M ij ∈ {0, 1}, which is Theorem 1.2 in Ref. [62], and then we present the Coefficients Theorem for weighted graphs, i.e., for M ij ∈ R. Before stating the Coefficients Theorem, we need the following definition: A linear directed graph is a directed graph for which it holds that all vertices have an indegree and outdegree equal to one. Hence, linear directed graphs are composed out of one or more directed cycles.
Theorem 1 (Coefficients Theorem for directed graphs (Milic [79], Sachs [80], and Spialter [81])). Let be the characteristic polynomial of an arbitrary directed graph G. Then a n = L∈Ln (−1) p(L) , with n = 1, 2, . . . , N, where L n is the set of all linear directed subgraphs L of G with exactly n vertices; p(L) denotes the number of connected components of L (i.e., the number of directed cycles of which L is composed).
Theorem 2 (Coefficients theorem for weighted directed graphs (Devadas Acharya where E(L) is the set of edges in the linear directed subgraph L.
In the case that G has no directed cycles, all coefficients a i in Theorem 1 or Theorem 2 are equal to zero, and Eq. (A.2) follows as a corollary of the Coefficients Theorem for directed graphs. for all i = 1, 2, . . . , N . As interactions are antagonistic, it holds that for all pairs (i, j) for which either M ij = 0 or M ji = 0. The tree condition implies that Eq. (34) holds. We assume that M ii = 0. The arguments we present are adapted from Ref. [3], albeit applied to the case of antagonistic, tree graphs.
First we define a general class of, so-called, strictly quasi-antisymmetric matrices, and we show that these matrices have purely imaginary eigenvalues. Second, we show that antagonistic trees are strictly quasi-antisymmetric.

Appendix B.1. Eigenvalues of strictly quasi-antisymmetric are imaginary
We say that a matrix Q ∈ R N ×N is strictly quasi-antisymmetric when where η is a symmetric, positive definite matrix; notice that Refs. [83,84] define strictly quasi-symmetric matrices, which are related to PT-symmetric matrices in quantum mechanics [85]. Strictly quasi-antisymmetric matrices have imaginary eigenvalues as they are similar to an antisymmetric matrix K. Indeed, since η is positive and symmetric, it is a diagonalisable matrix with positive eigenvalues. We define √ η as the square root of η that is positive definite, which is the (unique) symmetric, matrix that has eigenvalues that are equal to the positive square roots of the eigenvalues of η. Consequently, we may define the matrix which has real-valued entries as √ η has real-valued entries. The matrix K is antisymmetric, as where in the last step we have used that √ η is a symmetric matrix. Hence, since Q is similar to K, both matrices share the same eigenvalues [63], and since K is an antisymmetric with real-valued entries, the eigenvalues of Q are purely imaginary.
Appendix B.2. Adjacency matrices of trees with antagonistic weights are strictly-quasi-antisymmetric We show that the adjacency matrices M of antagonistic trees are strictly quasiantisymmetric i.e., they satisfy Eq. (B.3). To this aim, following Ref. [3], we explicitly construct the matrix η .
We select a random node in the graph, which we call the root node, and we label it as i 0 . Subsequently, we consider (i) the set of nodes V 1 that are neighbours of i 0 ; (ii) the set of nodes V 2 that are neighbours of nodes in V 1 excluding the root node; (iii) the set of nodes V 3 that are neighbours of V 2 and are not part of V 1 , etc.. Eventually we obtain a partitioning V = {i 0 } ∪ V 1 ∪ V 2 ∪ . . . V of the set of vertices of the tree graph G = (V, E) associated with M, where is the depth of the tree rooted at i 0 .
The matrix η is a diagonal matrix with elements on the diagonal defined as follows. We set the matrix entry associated to the root node to unity, and determine the other nodes through a recursion. In particular, we set for i ∈ V k and j ∈ V k+1 . Since M ij M ji < 0, the elements η ii > 0, and η is a symmetric, positive definite matrix. The procedure of constructing η is sketched in Fig. B1. Notice that the value of each element η ii is determined by its path to the root node i 0 , and this path is unique, as the graph is a tree. Next, we show that so that M is strictly quasi-antisymmetric. Components wise, the right-hand side of Eq. (B.8) reads If i = j, or i = j and the nodes are not each other's neighbours, then (B.10) Figure B1: Illustration of the partitioning of nodes in a tree with root node i 0 . The sets V 1 , V 2 , and V 3 contains nodes that are a distance 1, 2, are 3 separated from the root node i 0 .
On the other hand, when i ∈ V k and j ∈ V k+1 , then Eq. (B.7) applies, and we obtain where M is the adjacency matrix of a weighted, antagonistic tree as in Appendix B, and D is a diagonal matrix with nonnegative diagonal entries, i.e., [D] ii = D ii ≥ 0. Hence, the distinction with the matrix M in Appendix B is that the diagonal entries of M can be nonzero, and therefore for clarity we added the prime. We show that for this ensemble all eigenvalues have nonpositive real parts, i.e., (λ i (M )) ≤ 0 (C.2) for all i = 1, 2, . . . , N . The derivation consists of two parts. First, in Appendix C.1, we show that Eq. (C.2) holds for matrices K built from subtracting a nonnegative diagonal matrix to an antisymmetric matrix. Second, in Appendix C.2, we show that M is related to a matrix K by a similarity transformation, and hence they share the same eigenvalues. We end this appendix with a related result: in Appendix C.3 we show that the width spectrum of a diagonal matrix shrinks when we add to it an antisymmetric matrix.
Appendix C.1. Spectra of antisymmetric matrices with nonpositive diagonal entries Let us consider matrices of the form where K is an antisymmetric matrix (instead of the adjacency matrix of an antagonistic tree in Eq. (C.1)), and where D is a nonnegative, diagonal matrix. We show that for matrices of this form, (λ i (K )) ≤ 0 (C.4) for all i = 1, 2, . . . , N . Indeed, in this case, for any vector z ∈ C N , it holds that Since M is related to K by a similarity transformation, they share the same eigenvalues, and since K has nonpositive eigenvalues, as we we have shown in Sec. Appendix C.1, also M has nonpositive eigenvalues.
Using a similar line of reasoning it follows that (λ 1 (M )) ≤ −D min , (C. 18) where M is now of the form Eq. (C.1) with M the adjacency matrix of a weighted, antagonistic tree and D a diagonal matrix with diagonal entries that can be negative and positive.

Appendix D. Husimi Plateau in Jacobian-like matrices
We refine the results in Fig. 7 by analysing (λ 1 (J)) as a function of N in the limit of d min ≈ 0, where d min is the smallest value of d that belongs to the support set of p D (d). In this limit, the leading eigenvalue of antagonistic, Husimi trees exhibits strong transient effects as a function of N , and hence it is important to carefully extrapolate the results to large N .
To study the influence of a value d min ≈ 0 on the leading eigenvalue, we extract the diagonal entries D i from a distribution p D (d) = p a (d) that is plotted in Fig. D1. As illustrated by Fig. D1, the distribution p a (d) develops a peak around zero, i.e., d ≈ 0, whose weight increases a function a, which is the main reason why we use p a (d) and not the uniform distribution considered before in Fig. 7 where p HN (d) is a (narrow) half-normal distribution centered at zero, obtained by setting µ = 0 and σ = 0.05 on the right-hand side of Eq. (23), and p U (d) is a uniform distribution with a support that is not touching zero, and we denote its mean and standard deviation by µ U and σ U , respectively. We set the parameters µ U and σ U such that, as a varies, the first two cumulants of p a (d) are identical to those of p D (d) used in Secs. 4.2 and 4.3.  Figure D1: The distribution p a (d), as defined in Appendix D, for three values of the parameter: a = 0 (left), a = 0.1 (center) and a = 0.2 (right). Figure D2 plots (λ 1 (J)) as a function of N for antagonistic Husimi trees with diagonal elements drawn from p a (d) and for three values of a. In the left panel, the support of the diagonal distribution p a (d) does not include d = 0, and hence the functional behaviour of (λ 1 (J)) is analogous to the one shown in Fig.7. On the other hand, in the middle and right panels p a (d) develops a peak around d = 0, and consequently the monotonic increase of (λ 1 (J)) as a function of N slows down significantly at intermediate values of N , leading to the appearance of a plateau. Note that the plateau is a transient effect, as for large enough values of N the steady increase of (λ 1 (J)) continues. Observe that increasing a only widens the plateau from the left side, making it appear at smaller sizes of fN . Hence, also when p a (d) is peaked around d = 0, (λ 1 (J)) diverges as a function of N , albeit with a strong, transient, plateau effect.
We end with some final remarks. Due to the Husimi plateau, we should carefully assess finite size effects in Husimi trees. In particular, we could wrongly conclude that (λ 1 ) converges to a finite value when not considering large enough values of N . The Husimi plateau only occurs in Jacobian-like matrices when the support of the diagonal distribution contains d = 0, and hence we conjecture that it is related to the stripy structure of Jacobian-like matrices. In this Appendix, we determine the effect of a finite size N on the distribution p (λ 1 ) (x) plotted in Fig. 10. As we show, increasing the system size N , the discontinuity of the transition becomes more pronounced. Figure E1 plots the distribution p (λ 1 ) (x) for different sizes N when the control parameter s is roughly at the transition s d (in particular we set s 0.08, in correspondence with σ D = 0.10, σ G = 0.6 and µ G = 1.0). As shown in Fig. E1, the typical value (λ 1 ) of the distribution p (λ 1 ) (x) is, up to numerical accuracy, independent of N , confirming the discontinuous nature of the transition found in Fig. 10. Interestingly, as the size N increases, the γ-distribution  Figure E1: Histograms of the imaginary part of the leading eigenvalues of interactionlike (left panels) and Jacobian-like (right panels) matrices of antagonistic Erdös-Rényi graphs. The data shown corresponds with the parameter value σ D = 0.10 at which the dynamical transition takes place in Fig. 10, albeit now for different matrix sizes N as indicated in the legend. The data is plotted identically as in Fig. 10. Additionally, the inset shows the histograms for the range [0, 0.45] on the x-axis. For reasons of clarity, the figures at the bottom show the same data as the figures at the top, except for the absence of the histograms for N = 250 and N = 1000.
peaked at (λ 1 ) gets more narrow, indicating that in the infinite size limitp (λ 1 ) (x), as defined in Eq. (39), p (λ 1 ) (x) possibly converges to the delta distribution δ(x − (λ 1 )). Moreover, focusing on the behaviour of the nonreal histogram near the zero (see the insets in the top row of Fig.E1), we see find that the bin closest to the real delta peak decreases as a function of N , while the second bin increases as a function of N .