Eigenvalue spectra of finely structured random matrices

Random matrix theory allows for the deduction of stability criteria for complex systems using only a summary knowledge of the statistics of the interactions between components. As such, results like the well-known elliptical law are applicable in a myriad of different contexts. However, it is often assumed that all components of the complex system in question are statistically equivalent, which is unrealistic in many applications. Here, we introduce the concept of a finely-structured random matrix. These are random matrices with element-specific statistics, which can be used to model systems in which the individual components are statistically distinct. By supposing that the degree of `fine structure' in the matrix is small, we arrive at a succinct `modified' elliptical law. We demonstrate the direct applicability of our results to the niche and cascade models in theoretical ecology, as well as a model of a neural network, and a directed network with arbitrary degree distribution. The simple closed form of our central results allow us to draw broad qualitative conclusions about the effect of fine structure on stability.

The broad applicability of RMT is in part due to two key properties shared by many random matrix ensembles: the self averaging property and universality.Many ensembles are 'self averaging', meaning that the eigenvalue spectrum of a single realization approaches a deterministic (non-random) limit with increasing matrix size.The term 'universality' [39,40] indicates that this limit does not depend on the intricacies of the specific distribution from which the matrix elements are drawn, and only depends on a small number of its characteristics such as the first and second moments.Hence, in many cases we can do away with the specifics of a given system and focus instead on its statistical properties.
One central result of RMT is the celebrated elliptical law [29,41], which applies in a universal manner to a range of self-averaging random matrix ensembles, can be summarized as follows.If the elements of a large matrix are drawn in diagonally opposite pairs, independently and from the same distribution, then the eigenvalues will be confined to an ellipse in the complex plane, with the possibility of a single outlier.The shape of this ellipse and the location of the outlier eigenvalue are dependent only on the mean and variance of the elements, and on the correlation coefficient between diagonally opposite elements.If this large matrix represents the Jacobian matrix of some dynamical system, then the elliptical law tells us that we do not need to know the precise values of all the matrix elements in order to determine the stability of the system.Instead, we need only know the statistical properties of the interactions between the various components of the system.This observation allows for very general qualitative deductions about which characteristics of a complex system permit stability [32].
However, for the elliptical law to be applicable, it is necessary that all transpose pairs of elements be statistically equivalent.In recent applications, particularly in ecology and neural networks, this assumption breaks down.There is often a hierarchy, or ordering, to the elements which cannot be ignored [25,[42][43][44][45][46][47][48][49][50][51].As a result, many random matrix models have been developed to describe the various structures specific to the system in question, which result in modifications to the usual elliptical law.In this work, we aim to tie together a large subset of these models.
To this end, we introduce the concept of a 'finelystructured' random matrix (FSRM), which has elements with statistics that vary depending on the position in the matrix.By considering this fine structure as a perturbation to an elliptical random matrix (i.e. a random matrix for which the elliptic law applies [29,41,52]), we derive a correction to the elliptical law that depends only on a few summary statistics of the model.A similar correction for the associated outlier eigenvalue can also be determined.
The results we derive are simple and compact, so they are easily applied to many previously studied random matrix models.Indeed, we show that a wide variety of models can be understood as 'finely structured' upon careful analysis, and we demonstrate that our formu-arXiv:2311.02006v1[cond-mat.dis-nn]3 Nov 2023 lae accurately predict the eigenvalue spectra (and consequently the stability) of these systems.Ultimately, we use the modified elliptic law to draw qualitative conclusions about which aspects of fine structure affect the stability of disordered systems.
The rest of this work is laid out as follows.First, in Section II, we define more precisely what we mean by a finely-structured random matrix.In Section III we outline the derivation of the equations that determine the support of the spectrum of a general FSRM, and we use a perturbative approach to derive the modified elliptical law and associated outlier eigenvalue [Eqs.(10) and (13)], which approximate this support.Then, in Sections V to VIII we demonstrate how the modified elliptical law can be applied to a number of random matrix models.Specifically, in Section V we provide an exact solution for the cascade model from theoretical ecology [42,53], and show that this solution reduces to the modified elliptical law when the level of fine structure is small.In Section VI, we approximate the eigenvalue spectra of elliptical matrices with dense and random network structure.In Section VII, we apply our results to a toy neural network model inspired by Ref. [25].Finally, in Section VIII we show how the 'niche model' [32,54] can be understood as a finely-structured system and predict ecosystem stability in this model.We conclude by discussing possible extensions of the methods and results put forward in this work.

II. FINELY-STRUCTURED RANDOM MATRICES
Our focus is a particular ensemble of finely-structured random matrices (FSRMs), which includes many existing generalizations to the elliptical ensemble of Refs.[29,41] as special cases.Specifically, we consider FSRMs J of size N × N , whose elements have the following statistics where angular brackets ⟨•⟩ J , Var J (•) and Cov J (•) denote the average, variance and covariance over realizations of the random matrix J respectively.We note that if the higher moments of J ij drop off sufficiently quickly with N , then the results we obtain will be universal in the sense that we only need to know the moments in Eqs.(1), and not the full joint probability distribution of the elements of J, to find the eigenvalue spectrum [see Section S1 of the Supplemental Material (SM)].
We emphasize that the coefficients u ij , s ij and t ij are not random variables, but are instead fixed parameters of the model.We assume that all of u ij , s ij and t ij are of order N 0 .The factors of N in Eqs.(1) then ensure that we obtain a sensible large-N limit for the spectrum of J [22].The definitions of s and t also imply that we must have s ij ≥ 0 and t ij = t ji .
Similar random matrix ensembles have been studied in a number of different contexts [25,50,53,[55][56][57][58][59][60][61][62].However, as far as we are aware, none of these previous studies allow for non-zero u and t simultaneously, nor do they derive the central result of this paper, the modified elliptical law and associated outlier eigenvalue that we present in Section III.
The FSRMs that we study are a generalization of the elliptical ensemble in the following sense.If the statistics of the (i, j)-element do not depend on i and j (written u ij = µ, s ij = σ 2 and t ij = γσ 2 ), then the random matrix J is an elliptical random matrix of the type in [52].In this case, we say that the statistics of the FSRM J have no fine structure.In the more general case where the statistics of J do depend on position in the matrix, we find that the majority of the eigenvalues are still confined to a bulk region, which is generally not an ellipse, and there are O(N 0 ) outlier eigenvalues due to the non-zero mean values u ij .An example is shown in Fig. 1 (a).We see a 'star' shaped bulk region containing the majority of the eigenvalues, as well as four isolated outlier eigenvalues.
As is common in the study of disordered systems, we assume that the spectra of the large random matrices that we study are self-averaging [22].That is, we will assume that the eigenvalue spectrum, averaged over realizations of J, is equivalent to the spectrum of a single realization of the matrix for large N .We verify this assumption in panel (a) of Figs. 1 to 4, in which our analytical predictions, which describe ensemble averaged properties, are compared to single realizations of the FSRM models of interest.
Similarly to the random matrices studied in Refs.[37,52,63], we write the spectrum of J as the sum of the bulk spectrum and the isolated outlier eigenvalues as follows where λ i [J] are the eigenvalues of J indexed by i, ρ bulk (z) is the continuous bulk part of the spectrum at z = x + iy, and λ (k) outlier are the outlier eigenvalues, indexed by k.For a general random matrix J, the spectrum ρ J (z) can be derived from the disorder-averaged resolvent matrix (see Refs. [29,52,64]), defined by where I is the N × N identity matrix.Similarly to Refs.[25,29,37,65,66], the resolvent matrix is diagonal in the limit N → ∞, and we write G i (z) for the diagonal element [G(z)] ii .In Sections S1 and S2 of the SM, we find a set of equations that self-consistently determine the {G i (z)} for the FSRM ensemble defined in Eqs.(1) using well-known replica methods [29,37,67,68].We then describe how these expressions for the {G i (z)} can, in principle, be solved and used to obtain the spectrum of J.
The equations relating the {G i (z)}, ρ J (z) and the statistics of J are complicated.Non-trivial choices of the matrix statistics for which ρ J (z) can be found explicitly are difficult to come by.Further, when an explicit solution exists, it can be rather complex (see Section V, and Section S6 of the SM for examples), and this complexity can obfuscate the qualitative effect of fine structure on stability.
Therefore, we proceed by treating the fine structure as a perturbation, seeking a fine-structure correction to the elliptical law and outlier eigenvalue which can be interpreted qualitatively.We write the statistics in Eqs.(1) as a sum of index-independent and index-dependent parts where the 'zeroth-order', or elliptical, statistics are (in the limit of large N ) and where we call the matrices u (1) , s (1) , t (1) the finestructure parts of u, s, t respectively.We also define the following 'first-order' (or fine-structure) parameters (again to be understood in the large-N limit), jk .
As we will show, these fully characterize the finestructure correction to the elliptical law and outlier eigenvalue when the amount of fine structure is small.We can interpret these new parameters by noticing that R, S, T, U and V are directly related to the row and column sums of the statistics of J.For example, if we define the i-th row sum of t as t i ≡ j t ij /N , then we have Hence, up to a pre-factor, T is the variance of the row sums of the matrix t.
If the fine-structure contributions are all zero, then the spectrum of J is simply the well-known elliptical law (see Section S6 A 1 of the SM) [29,52].That is, the majority of the eigenvalues are contained in the ellipse and, provided |µ| > σ, there is a single outlier eigenvalue located at (see [52,64,69]),

B. Modified elliptical law
The fine-structure correction to the support of the spectrum of J can be computed by introducing a small perturbation parameter ϵ which measures the extent to which the statistics of J are finely structured.We assume that all elements of u (1) , s (1) and t (1) are proportional to this small quantity ϵ so that the fine structure parameters in Eqs. ( 6) are proportional to ϵ 2 .We then find that the leading-order correction to the elliptical law is not of order ϵ, but is in fact of order ϵ 2 .Here, we only present the final result.For a full and detailed derivation of all results up to and including Eq. ( 12), we refer the reader to Sections S3 and S4 of the SM.All relationships derived in this section are understood to be accurate up to second order in ϵ.
By supposing that the diagonal elements of the resolvent G i (z) have the following small-ϵ expansion we can find an expression for the trace of the resolvent i G i (z)/N which is accurate to second order in ϵ.Using standard methods (see e.g.Ref. [29]), the trace of the resolvent can then be related to the eigenvalue density in the complex plane, and we can thus determine the support of the eigenvalue spectrum.
The support of the bulk spectrum can be expressed in terms of the elliptical parameters σ 2 , γ, as well as the fine-structure parameters R, S and T .In the end, we find the following deformed ellipse with When S, T, R → 0, we recover the usual elliptical law in Eq. ( 7).In Figs. 2 to 4, rather than plot Eq. ( 10) directly by solving for x or y, we use the following parametric representation For a demonstration of the equivalence of this parametric form to Eq. ( 10), see Section S4 A 1 of the SM.The quantity a determines the location of the rightmost edge of the bulk region λ edge , and hence stability.More precisely, we have C. Outlier eigenvalue Within a calculation to second order in ϵ, we find at most one outlier eigenvalue in the spectrum of an FSRM.The outlier can be expressed in terms of the elliptical parameters µ, σ 2 , γ, as well as the fine-structure parameters T, U and V as follows (see Section S4 B in the SM) Similarly to how an elliptical matrix has no outlier if |µ| ≤ σ, the finely structured matrix under consideration only has an outlier at Eq. ( 13) if the following condition holds In the case of an equality in Eq. ( 14), we have λ outlier = λ edge , and the outlier eigenvalue is 'absorbed' into the bulk spectrum.Note that for T, U, V → 0, we recover the outlier for elliptical matrices in Eq. ( 8).Eqs. ( 10), ( 12) and ( 13) are the central results of this paper.Given an arbitrary finely-structured random matrix, Eqs. ( 10) and ( 13) provide an explicit and direct approximation for the support of the spectrum of that matrix.The fine-structure correction to the leading eigenvalue of a general FSRM is given by the maximum of λ edge and λ outlier .We can therefore make general statements about the stability of systems for which the Jacobian matrix is finely-structured using Eqs.( 12) and ( 13) (when the strength of the fine structure is small).

D. Stability
In this section, we use Eqs.( 12) and ( 13) to understand in more detail what kinds of fine structure promote (in)stability.That is, if the matrix J − 1 1 is the Jacobian matrix of a dynamical system linearized about its fixed point (as it is for the models in Sections V to VIII), we see that if any of the eigenvalues of J have a real part that exceeds 1, the fixed point is unstable to perturbations.Thus, by understanding how the fine structure affects the rightmost eigenvalue of J, we gain insight into the effect of fine structure on stability.All mathematical details in this section can be found in Section S5 of the SM.
In our analysis we use the following matrix semi-norm [70] which satisfies ∥A∥ ≥ 0, and where ∥A∥ = 0 if and only if the row sums of A are all equal to zero.We also define the following quantities where s (1)T denotes the transpose of the matrix s (1) , and similarly for u (1)T As ∥ • ∥ is a matrix semi-norm, we can interpret S s as quantifying the amount of fine structure present in the matrix (s + s T )/2, the symmetric part of s.We also interpret S a as quantifying the amount of fine structure in the antisymmetric part of s, and similarly for U s and U a .Because T = ∥t (1) ∥ 2 /σ 4 [see Eqs.(6)] and t (1)T = t (1) , we interpret the parameter T as directly quantifying the amount of (symmetric) fine structure in the matrix t.
The fine-structure parameters R, S, U and V are related to S s , S a , U s , U a and T via The coefficients ρ 1 and ρ 2 , defined in Eqs.(S89) in the SM, are numbers that vary in the range [−1, 1].Their precise values are immaterial to the arguments that follow.Substituting Eqs.(17) into Eq.( 12), the following deductions of the effect of fine-structure on the bulk edge λ edge are drawn.We find that increased antisymmetric fine structure in the matrix s (as measured by S a ) always decreases the value of λ edge .That is If γ ≥ −1/3, then symmetric fine structure in the matrices t and s increases the value of λ edge If γ < −1/3, then symmetric fine structure in the matrices t and s can be stabilizing or destabilizing.We can make similar deductions for the outlier eigenvalue by substituting Eqs.(17) into Eq.( 13).Using Eq. ( 14), we see that Eq. ( 13) only corresponds to an outlier eigenvalue to the right of the bulk spectrum if µ/σ > 1 + δ, where δ goes to zero as the fine-structure correction goes to zero.Hence, to leading order in the fine structure, the factor of (1 − γσ 2 /µ 2 ) in Eq. ( 13) is positive.Using this observation, we find that antisymmetric fine structure in the matrix u decreases the outlier eigenvalue and the parameters U s and T increase the outlier eigenvalue, for any value of γ We can therefore draw broad qualitative conclusions on the effect of fine structure on stability.Equations (18) and (20) tell us that antisymmetric fine structure in the statistics of an FSRM is always a stabilizing influence, and Eqs.(19) and (21) tell us that symmetry in the fine structure of the statistics of a FSRM largely promotes instability.Indeed, for symmetric fine structure to stabilize a given system, it is necessary, but not sufficient, that γ lies in the range [−1, −1/3] and that the inequality in Eq. ( 14) is violated.

IV. SUMMARY OF THE EXAMPLES IN THE FOLLOWING SECTIONS
The remainder of this paper is concerned with the application of the results for the fine-structure corrections to the eigenvalue spectrum presented thus far [Eqs.(10) and (13)].The purpose of these examples is to make clear what constitutes a finely-structured random matrix in context, and to demonstrate the breadth of the applicability of the modified elliptical law.
In the first example (Section V) we study a generalization of the cascade model from theoretical ecology [42,71].This is a simple example of an FSRM ensemble, for which there is an exact solution, which we derive in Section S6 B of the SM.We evaluate the fine-structure correction to the elliptical law for this model using the modified elliptical law, and verify that this correction produces the same result as an expansion of the exact solution for small fine structure.An instance of this model, and our solution, is shown in Fig. 1.
In a second example (Section VI), for which there is no simple exact solution, we show that directed complex networks with dense random network structure can also be treated as FSRMs.We find that the amount of fine structure present depends on three statistical properties of the network, and we derive the fine-structure correction to the elliptical law, which is valid when these statistics are small.The results can be seen in Fig. 2. We also note that our fine-structure corrections reproduce the exact result of Ref. [72] when the network has no undirected links, and the formulae reproduce the results of [65] when the network has no directed links.
The third example (Section VII) considers a square grid of interacting neurons arranged in physical space [25].The interaction between any two neurons depends continuously on the physical distance between them in the grid.We use Eqs.(10) and (13) to compute the finestructure correction to the spectrum of random matrices with statistics modelling the interacting neurons, arriving at a simple closed-form expression, where previous works had to find the spectrum of similar models numerically.Our approximation is compared with numerical results in this case in Fig. 3.
In the final example (Section VIII), we consider a more sophisticated model in theoretical ecology.We use our central results to provide a simple approximation for the eigenvalue spectrum of a random matrix constructed according to the niche model [32,54], for which precise analytical results have so far been difficult to come by.In the case of the niche model, there is no limit of any model parameter one could take to recover the conventional elliptical law.Despite this, our fine-structure correction still works well (an example is shown in Fig. 4).27) and (28).Each marker is the result of diagonalizing a single 8000 × 8000 instance of the matrix.

V. THE CASCADE MODEL OF COMPLEX ECOSYSTEMS A. Model Definition
Following Robert May [33,34], we consider a set of ordinary differential equations that govern the time evolution of a set of species abundances x i , ẋi = f i (x 1 , x 2 , . . ., x N ).(22) We imagine that this system is linearized about a fixed point such that where δx i is the deviation of x i from its fixed point value.
The diagonal elements of the Jacobian are set to −1 so that the system without interactions would be stable.We wish to understand what kinds of interactions (or Jacobian matrix elements J ij ) permit stability.If the largest real part of any eigenvalue of J is greater than one, then the ecological equilibrium will be unstable.May's original work makes use of the circular law [40,73,74], a special case of the elliptical law with zero mean (µ = 0) and uncorrelated matrix entries (γ = 0).More recent work [75][76][77], has given an ecological interpretation to the parameters µ and γ, associating them respectively with the average strength of interspecies interaction and with the proportion of interspecies interactions that are of predator prey type p (the proportion of links satisfying J ij J ji < 0) through γ = cos(πp) (see also Ref. [78]).Here, and later in Section VIII, we also incorporate the possibility of a hierarchy amongst species [42,78], which can be interpreted as fine structure.
We consider the FSRM ensemble with statistics given by where u, s, t are defined in Eqs.(1) [the structure of these matrices is illustrated in panels (u), (s) and (t) of Fig. 1].If γ 1 = 0, then we recover the cascade model of Refs.[42,71].The triangular structure of u reflects the hierarchy amongst species (see Refs. [42,78] for more in-depth discussion).The generalization presented here allows for t ij to be index-dependent.In writing down Eqs. (24) we have anticipated that µ, σ 2 and γ in Eqs.(24) will coincide with the quantities defined in Eqs.(5).

B. Eigenvalue spectrum
To find the fine-structure correction to the elliptical law due to the statistics in Eqs.(24), we suppose that µ 1 , σ 2  1 and γ 1 are proportional to a small parameter ϵ.We make the following substitutions, If ϵ = 0 then the statistics Eqs.(24) reduce to the elliptical law with parameters µ, σ 2 and γ [see Eqs.(5)].
Let us now find the correction due to the fine structure in the cascade model.The calculation of the finestructure parameters R, S, T, U and V [using Eqs.(6)] is detailed in Section S7 C of the SM.We find R = V = 0, and the non-zero terms are Substituting the corrections from Eq. ( 26) into Eqs.( 10), ( 12) and ( 13), we obtain the fine-structure correction to the spectrum of an FSRM with cascade model statistics.We can see how the addition of fine structure affects stability by looking at the correction to the bulk edge [Eq.(12)] and the fine-structure correction to the outlier eigenvalue [Eq.( 13)] Recalling |γ| ≤ 1, we see that λ edge increases if γ 1 is increased and decreases if σ 1 increases.Further, Eq. ( 14) tells us that Eq. ( 28) only corresponds to an outlier eigenvalue if µ/σ > 1 + O ϵ 2 .Hence, when Eq. ( 28) is valid, the term ϵ 2 (1 − γσ 2 /µ 2 ) is positive (to second order in ϵ) so we can conclude that γ 1 increases, and µ 1 decreases the value of λ outlier .As σ 1 , µ 1 encode the amount of asymmetry in the fine-structured parts of the matrices s and u, and as γ 1 encodes the amount of symmetrical fine structure in the matrix t, our findings here are consistent with the discussion in Section III D. The above stated dependence of stability on the various model parameters is also confirmed in Fig. 1 (b).
For FSRMs with statistics as in Eqs.(24), there exist exact expressions for the support of the bulk region and outlier eigenvalues.The solution is derived in Section S6 B of the SM, and verified in panels (a) and (b) of Fig. 1.In Section S7 of the SM, we also demonstrate that expanding the exact solution in small ϵ leads to the same result as the modified elliptical law, thus verifying Eqs.(10), (12) and ( 13) analytically.

VI. DIRECTED COMPLEX NETWORKS AS FINE STRUCTURE A. Model definition
As a second example, we consider the eigenvalue spectra of matrices that represent weighted and directed complex networks.Here, the network structure is the source of fine structure in the model.
In particular, we consider matrices of the form where A is the adjacency matrix of the network (A ij = 1 if a link exists going from node j to node i and is zero otherwise), K encodes the weights of the edges, which can be positive or negative, and • denotes the element-wise product of matrices.Here, K is an elliptical random matrix with statistics where p = ij A ij /N is the average degree of nodes in the network.The scaling with p of the above statistics ensures a sensible dense limit p → ∞, with p/N held constant as N → ∞.As the statistics of K contain no fine structure, our corrections to the elliptical law and outlier eigenvalue will only depend on properties of the network.
The network is constructed as follows: Each pair of nodes i and j is joined with a directed edge i → j, another directed edge j → i, a reciprocal edge i ↔ j, or no edge with probabilities P j←i , P i←j , P i↔j and P i̸ ↔j respectively (P j←i +P i←j +P i↔j +P i̸ ↔j = 1).These events are taken to be mutually exclusive, so our definitions are such that we cannot obtain an undirected link from two directed links.We also define the reciprocity (the ratio of undirected links to total links) [79][80][81][82] of the network as well as the mean undirected, exclusively-in and exclusively-out degrees of each node The follow from the definitions.Note that even if the underlying network is completely undirected (r = 0), the weighted network A • K can still be , the bulk edge without the fine-structure correction.The top three curves when r = 0 correspond to random matrices with positive degree correlations ρ, the bottom three curves correspond to ensembles with the same degree distributions with the links rearranged so that the degree correlation is negative, demonstrating that increasing ρ leads to destabilization.When r = 1, the degree heterogeneity h 2 controls the size of the fine-structure correction, and it stabilizes/destabilizes depending on the value of Γ.The network has the same dichotomous structure as that used in panel (a) and µ = 0.5, σ = 1 and Γ as indicated.Markers are the result of diagonalizing a single 8000 × 8000 realization of the matrix.We verify Eq. ( 36), the fine-structure correction to the outlier eigenvalue, in Fig. S2 of the SM.
considered directed because of the presence of positive and negative weights (see, for example, Ref. [65]).
In Section S8 A of the SM, we show that J = A • K (the random matrix of interest) has the same spectrum as a fully-connected FSRM with statistics Hence, network structure can be interpreted as a manifestation of fine structure, and we can approximate the spectrum of A • K using the modified elliptical law.

B. Eigenvalue spectrum
To find the fine-structure correction to the spectrum we first compute the elliptical parameters [Eqs.(5)] from the statistics Eqs.(32).The definitions of µ and σ 2 in Eqs.(32) coincide with those in Eqs. ( 5), and we find γ = Γr.
The fine-structure parameters R, S, T, U, V are computed directly using Eqs.( 6), and they depend on three additional statistical properties of the network.These are the variance of the undirected degrees h 2 , otherwise known as the degree heterogeneity [65], the correlation coefficient between exclusively undirected degrees and exclusively directed degrees τ , and the correlation coefficient between exclusively in and out-degrees ρ, otherwise known as the degree correlation coefficient [72].Explicitly, these statistics are Written in terms of h 2 , τ and ρ, the fine-structure param-eters are then The fine-structure correction to the spectrum of J due to the network is found by substituting the values of R, S, T, U and V into Eqs.( 10) and ( 13).This approximation is valid provided h 2 , τ , and ρ are all small.For details of the calculation, see Section S8 A of the SM.
We can analyse the effect of network structure on stability by looking at the fine-structure correction to the bulk edge [see Eq. ( 12)], and the fine-structure correction to the outlier eigenvalue [see Eq. ( 13)], Examining Eqs. ( 35) and ( 36), and recalling that there is no outlier if µ/σ < 1 + δ, where δ goes to zero as the fine structure goes to zero, we find that the network heterogeneity h 2 always increases the value of the outlier eigenvalue, and increases the value of λ edge if Γ > 1 − √ 2 or if r < 0.73 [as can be seen by examining the factor multiplying h 2 in Eq. ( 35)].We also observe that λ edge and λ outlier always increase with increasing correlation coefficient ρ.Given that a greater values of h 2 and ρ both connote greater degrees of symmetric fine structure, we see that our results agree with the general remarks made in Section III D.
In Section S8 A of the SM, we further confirm our results by demonstrating that our fine-structure correction [Eqs.(35) and (36)] reduces to the approximate results derived in Ref. [65] when the underlying network is exclusively directed (r = 1), and reduces to the exact results of Ref. [72] when the underlying network is exclusively undirected (r = 0).

VII. FINE STRUCTURE IN A TOY NEURAL NETWORK MODEL A. Model Definition
In the study of neural networks, firing-rate models are used to investigate the dynamical interaction of the neurons.The activation of the i-th neuron x i is often taken to follow dynamical rules of the type [27,49] where J ij dictates the presence and strength of an effect of neuron j onto the activation of neuron i. Similarly to May's analysis of complex ecosystems, when the leading eigenvalue of J is greater than 1, the state x i = 0 becomes unstable, leading to a qualitative change in the dynamics.
Inspired by a toy model proposed in [25], we consider N = n 2 neurons regularly spaced on a finite square grid (with no periodic boundary conditions), and with correlated interactions.We suppose that the degree to which the interaction strengths J ij and J ji between any two neurons are correlated is a function of the physical distance between them in the grid [83,84].For the purposes of illustration of the method, we consider an FSRM ensemble with statistics that depend in a simple way on the distances between neurons where r i = (x i , y i ) is the position of the i-th neuron in the grid and k is a parameter controlling the scale over which correlations in the interaction of one neuron with other neurons varies.We use a lattice spacing of 1/n, so that the total size of the grid is 1 × 1.
We suppose that the neurons are labelled 1, 2, . . ., n in the top row, n + 1, n + 2, . . ., 2n in the second row, and so on, so that there are n 2 neurons in total.The matrix t therefore has size n 2 × n 2 .An example with n = 9 is shown in Fig. 3 (so the matrix of statistics t ij has size 81 × 81).The block-like structure is due to the specific labelling of the neurons; there is a discontinuous jump in the function |r i − r j | at the end of each row.Ultimately, the labelling of the neurons has no effect on the spectrum of J.
If the range parameter k is zero, then we recover the usual elliptical law, so it is fair to assume that the approximation works if k ≪ 1.However, as Fig. 3 (b) demonstrates, our approximation works well for most values of k, not just small values.This is because the fine-structure correction to the elliptical law is small for all values of k.Assuming that the number of neurons in the grid is large, we can compute the elliptical and fine-structure parameters using Eqs.( 5) and ( 6).The mean and variance are µ and σ 2 respectively, and the correlation coefficient is where the F 1 (k) is an integral over the coordinates of the grid with Details of the derivation of Eq. ( 39), as well as of all other calculations in this section, can be found in Section S8 B of the SM.From Eqs. (38), we see that the matrices u and s have no fine structure.Hence, the parameter T is the only non-zero fine-structure parameter where the function F 2 (k) is also an integral over the spatial coordinates of the grid Plugging the values of µ, σ 2 , γ and T into Eqs.( 10), ( 12) and ( 13) gives the fine-structure correction for FSRMs with statistics as in Eqs.(38).
When k = 0, the statistics in Eqs.(38) reduce to those of an elliptical random matrix.Therefore, the amount of fine structure present in the system is guaranteed to be small if k is small.Expanding Eqs. ( 12) and ( 13) in powers of k gives us some insight into the (de)stabilizing effect of the spatial dependence of interactions between neurons We have divided by the zeroth order factors in the fine structure to highlight the relative effect of the fine structure on the spectrum.Our approximation for the bulk and outlier are confirmed in Fig. 3.In particular, panel (b) reveals that the fine structure present in the statistics Eqs.(38) is always destabilizing, which is consistent with the general condition in Eqs.(19).

VIII. THE NICHE MODEL OF COMPLEX ECOSYSTEMS A. Model definition
Our final example builds on the same basic randommatrix framework of Robert May, as described in Section V, but involves a more intricate interaction network.
The niche model is primarily a model of ecological network structure [54], which approximates many statistical features of real ecological networks [36,85].However, due to its complexity, it is a daunting task to find an exact analytical expression for the eigenvalue spectra of the associated interaction matrices.Here we show how Eqs. ( 10) and ( 13) deliver a simple approximate solution to this problem.
First, we describe the steps for the construction of a niche model network.To construct a realization of the network, we draw three random numbers for each of the N species in the system, the 'niche value' η i , 'niche range' d i , and 'niche centre' c i [54,86], whose distributions are described below.These quantities determine where in the hierarchy a species lies, how many species it interacts with, and where in the hierarchy the species it interacts with are, respectively.The network structure is encoded in the matrix E, with Following Ref. [54], the niche value, range and centre of species i are defined by the following steps: 1.The niche value η i is a uniform random variable between 0 and 1. Species are re-labelled so that Species with high niche value are at the top of the hierarchy, and species with low niche value are at the bottom.
2. The niche range is d i = η i X i , where X i is a beta distributed random variable sampled independently for each species from the distribution p( The value of the parameter d is the same for all species, and it is equal to the average niche range By definition, the niche centre of a species is always smaller than that species' niche value, reflecting the idea that species generally 'look down' the hierarchy for food.
The structure of E, and therefore of the niche model network, varies solely due to the average niche range d.
A single realization of E is illustrated in Fig. 5(a).Given a realization of the network E, a random matrix constructed according to the niche model has elements , where K L ij are independent random variables with mean µ L /N and variance σ 2 L /N , and where K U ij are also independent random variables, with mean µ U /N and variance σ 2 U /N .The elements of K L and K U are also correlated via For a fixed E, the matrix J is therefore an FSRM, with statistics Elements of the adjacency matrix A ij are equal to one if either of E ij or E ji are equal to one.Otherwise, A ij = 0. We also assume that µ L > µ U , which implies that if E ij = 1 and E ji = 0, then species i benefits more from the relationship than species j on average.If E ij = E ji = 1, then on average species i and species j benefit an equal amount from their mutual relationship.See Refs.[54,86] for discussion and motivation of the model.
In the following section we will compute the finestructure correction using Eqs.(45), that is, we compute the fine-structure parameters for a fixed instance of the network.However, as we will see, the fine-structure correction depends only on a few summary statistics of {E ij }.

B. Eigenvalue spectrum
We now proceed to calculate the zeroth-order parameters in the fine structure [using Eqs. ( 5)].They are The fine-structure correction [Eqs.(6)] can be expressed in terms of two additional properties of the underlying niche model network: the degree heterogeneity of the niche model network with adjacency matrix A, and a property encoding the amount of asymmetry in the matrix E We also introduce the following parameters, defined in terms of µ L , µ U , σ L , σ U , which measure the amount of asymmetry in the interspecies interaction strengths and variances In terms of the aforementioned quantities, the finestructure parameters read [see Section S8 C for details of the calculation] Therefore, the fine-structure correction to the elliptical law and associated outlier eigenvalue only depend on three characteristics of the niche model network: d, h 2 and f , and not on any particular instance of E. If σ 1 = µ 1 = 0, then the niche model is a special case of the FSRMs in Section VI where the adjacency matrix of the underlying network is symmetric (r = 1).In this case Eqs.(34) and Eqs.(49) give the same results.
Plugging the elliptical [Eqs.(46)] and fine-structure [Eqs.(49)] parameters into Eqs.(10) and (13) gives the fine-structure correction to the support of the spectrum of a random matrix constructed according to the niche model.In particular, we can analyze the combined effect of the network and hierarchical interactions on stability.We find the following correction to the bulk edge and the correction to the outlier eigenvalue is We note that the quantity f is always negative.Recalling the restriction |γ| ≤ 1, and that there is no outlier if µ/σ ≤ 1 + δ, where δ vanishes as the fine structure vanishes, we see that increased asymmetry in both the interaction strengths (i.e.increased µ 1 and σ 2 1 ) and network structure (i.e. more negative f ) reduce λ edge and λ outlier .Hence, increased asymmetry is always stabilizing.The effect of degree heterogeneity is more complicated.If there are no asymmetrical interactions (µ 1 = σ 1 = 0), then larger degree heterogeneity always increases the value of λ outlier and increases λ edge provided γ > 1 − √ 2. If there are asymmetrical interactions, then increased degree heterogeneity can be stabilizing or destabilizing, as it increases the amount of symmetrical and asymmetrical fine structure in the matrices u and s.
As the network structure ultimately depends only on the average niche range d, the parameters h 2 and f must be functions of d only.In Section S8 C of the SM, we find approximate expressions for h 2 and f , valid when d is small, but not so small as to make the network sparse Equations ( 52) are compared to the measured values [obtained from generating E according to rules (1) to ( 4)] in Fig. 5.We find good numerical agreement for all possible values of d, not just small ones.The validity of our results, as well as our assertions on the effect of µ 1 on stability, are verified in Fig. 4 for varying strength of hierarchy between species, where we see that increasing µ 1 is indeed stabilizing.The validity of our assertions on the effect of the degree heterogeneity of the network are verified in Fig. S3 of the SM.The performance of our approximations is surprisingly good considering that the size of the fine-structure correction can be as large as ≈ 70% of the uncorrected value.Figure 5 (b) also reveals that the degree heterogeneity is always non-zero, and has a minimum value of approximately 0.1.Hence, any combination of the model parameters leads to a non-vanishing fine-structure correction.

IX. CONCLUSIONS
The elliptical law carries information about the stability of a densely-connected disordered system in which all components are statistically equivalent.However, in many cases, we have some knowledge about the system's structure and want to know how this feature affects its stability.In this work, we have provided an explicit formula for a large subset of such models that answers this question, so long as the effect of the additional structure is small.Further, our modified elliptical law reveals that antisymmetry in the statistics (as opposed to the elements themselves) of an FSRM is stabilizing, and symmetry is usually destabilizing.We have seen how this observation manifests itself in a number of examples.For instance, by interpreting hierarchy in an ecosystem as antisymmetric fine structure in the statistics of the associated randommatrix models, we are able to conclude that hierarchy is a stabilizing factor.Also, by interpreting the degree heterogeneity of a network as symmetric fine structure, we can conclude that degree heterogeneity is usually a destabilizing factor.
There are several possible avenues for future work.We have treated fine structure as a perturbation to the elliptic law, but there are other aspects which could be treated similarly.For example, the corrections to the elliptic law due to a sparse network structure [87], or additional correlations [88,89] have been addressed in this way.Combining such corrections using a perturbative approach is interesting prospect for future work.Alternatively, we could consider the fine-structure corrections to the stability criteria of some non-linear dynamical system with an FSRM as its interaction matrix, such as the generalized Lotka-Volterra equations [78].
A further interesting possibility for future work would be to drop the restriction that all mean values of the matrix elements, ⟨J ij ⟩ J , scale with 1/N , instead allowing for some to be of order 1, perhaps by imposing conditions similar to those in Ref. [90].We could then consider systems where agents interact with some agents very strongly, and with others randomly.For example, many ecological models include self-regulation of species (see, e.g., Refs.[37,43,48,91,92]), which could be incorporated into an FSRM for which the diagonal elements of u are of order 1, rather than of order 1/N .
Finally, we remark that one can compute all parameters necessary for the modified elliptical law with just one instance of the FSRM J.This opens up the possibility of applications of the modified elliptical law to situations with real-world data, in a similar spirit to uses of the standard elliptical law [93].In Section S1, we compute the disorder average of the resolvent matrix of the N × N finely-structured random matrix (FSRM) J [Eq. (3) in the main text] by means of a saddle point approximation for large N of the so-called eigenvalue potential [29].Given that the disorder-averaged resolvent matrix is diagonal, we then derive a set of self-consistent equations that relate the diagonal elements of the resolvent matrix, the statistics of the FSRM ensemble, and the support of the spectrum of J in Section S2.

ACKNOWLEDGMENTS
In Section S3, we expand the self-consistent equations derived in Section S2 to second order in a small parameter ϵ, which measures the amount of fine structure present in the system.We then solve these simplified equations explicitly in Section S4, deriving our modified elliptical law and outlier eigenvalue as consequences [Eqs.(10) and ( 13) in the main text].
In Section S5, we use the modified elliptical law and outlier eigenvalue to justify the claims made in Section III D of the main text.In particular, we demonstrate that antisymmetric fine structure in the statistics of a FSRM generally stabilizes the system, and symmetric fine structure largely destabilizes the system unless γ < −1/3.If γ < −1/3 then symmetric fine structure can stabilize or destabilize the system.
In Section S6, we analyze a number of FSRM ensembles for which the self-consistent equations for the support of the spectrum can be solved exactly.Specifically, in Section S6 A we present an exact solution in the case where the fine-structure parameter T = 0.This class of FSRM models includes the elliptical law, the cascade model of Ref. [42], and the directed networks model of Ref. [72] as special cases.In Section S6 B we find the support of the spectrum for random matrices constructed according to the cascade model as presented in Section V of the main text.
In Section S7, we demonstrate that the modified elliptical law gives the same result as a direct expansion of the exact result in the case of the cascade model of Section V in the main text, verifying our small-fine-structure correction analytically.
Finally, in Section S8 we detail the computation of the elliptical and fine-structure parameters for the models presented in Sections VI to VIII of the main text.Specifically, we compute the parameters µ, σ, γ, R, S, T, U and V as defined in Eqs. ( 5) and ( 6) of the main text for a directed network, a neural network model and for the niche model from theoretical ecology.

A. The eigenvalue potential
To find the spectrum of an N × N finely-structured matrix J with element specific statistics, we calculate the disorder averaged resolvent matrix, given by Eq. ( 3) in the main text.It is well-known that the density of the bulk region can be derived from the disorder averaged resolvent of J via [see Ref. [94] for detailed derivations of Eqs.(S1) to (S3)] where ∂ z = (∂ x + i∂ y )/2 and Tr denotes the trace.
From Eq. (S1), we see that the eigenvalue density is non-zero when Tr G(z) is a non-analytic function of z.Hence, the support of ρ bulk (z) is the boundary of the region in the complex plane for which Tr G(z) is an analytic function.
Following [29], the problem of computing Tr G(z) can be tackled by considering the eigenvalue potential where z denotes the complex conjugate of z.The potential is related to the resolvent matrix through In the large-N limit, we will compute Φ(z, z) using the replica trick together with a saddle point approximation [29,67,89,95,96].The result is ultimately a set of self-consistent equations which determine the resolvent matrix G(z), as well as spectral density ρ bulk (z).
For the calculation of the eigenvalue potential it is convenient to write the FSRM J [with statistics given by Eqs.(1) in the main text] as the sum of a mean and a fluctuating part where the mean zero random variables w ij satisfy The correlation parameters γ ij are related to the covariance parameters t ij , defined in Eqs.(1) in the main text, via Similarly to [29,37,67], we evaluate the eigenvalue potential Φ(z, z) using replicas.This involves using the identity ln x = lim n→0 (x n − 1)/n, but assuming that n is an integer.This reduces the problem of computing the disorder average of a logarithm to the problem of computing the disorder average of the 'replicas' ⟨x n ⟩ J .
In previous works, it has been shown that the replicas decouple [29,37,65,67], that is, that ⟨x n ⟩ J = ⟨x⟩ n J .Hence, if the replicas decouple, then the average of the logarithm is the logarithm of the average.In this work, we assume that the replicas decouple, so that We will also follow [52,[63][64][65] and assume that if the rank of the matrix u is finite as N increases, then it has an O(1/N ) contribution to the bulk spectrum and can therefore be set to zero for the derivation of the bulk spectrum.
Whilst the elements u ij turn out to have no effect on the bulk spectrum, we will see in section Section S2 B that they are important in determining the location of possible outlier eigenvalues.

B. Disorder average
To carry out the average of the determinant in Eq. (S6) over realizations of J, we first express it as the following Gaussian integral over variables p i=1,...,N where the second equality follows from a Hubbard Stratonovich transformation [95,96].Substituting Eq. (S4), the ensemble average of the disordered part of J can be computed by taylor expanding the exponent in powers of 1/N .One can then re-exponentiate the result, assuming that higher moments decay sufficiently fast with the system size.
Focusing on the parts of the exponent which contain the disordered terms (terms containing the random variable w ij ) we have Now that we have averaged over realizations of J, our expression for the eigenvalue potential Φ(z, z) reads We can perform the integral in Eq. (S9) for large N by transforming the integrand into a form amenable to a saddle point approximation.The transformation involves introducing the following parameters Following [29,65,69], we neglect terms involving (p i ) 2 , (q i ) 2 , p i q i and p i q i as these do not ultimately contribute to the integral to leading order in N .We introduce P i , Q i , R i , R * i into Eq.(S9) by inserting Dirac deltas in their complex exponential form where we have ignored the pre-factors involved in expressing the deltas as complex exponentials because they only contribute a constant factor to the eigenvalue potential, which does not affect the spectrum.
With the complex exponentials inserted, our expression for the eigenvalue potential now reads The integral in the final term is a standard Gaussian integral, it can be explicitly evaluated As the integrand in Eq. ( S12) is now of order exp(N ), we can evaluate the integrals over the parameters and their hatted counterparts with a saddle point approximation.For the hatted parameters, we arrive at the following saddle point equations and for the unhatted parameters Using Eqs.(S14), we find Eliminating the hatted parameters from Eqs. (S14) and (S15), we arrive at.
Similarly to Refs.[37,65], we can relate the saddle point values of R i , R * i to the diagonal elements of the resolvent matrix G i (z).First, we suppose that the eigenvalue potential depended not just on one value of z, but on a diagonal matrix z of complex variables with elements Following the derivation of Eq. (S12) with this alternative definition of the eigenvalue potential, the only difference is that we replace z → z i and z → z i throughout.Comparing the partial derivative of the alternative potential with respect to z i in Eqs.(S12) and (S21) gives As the eigenvalue potential of the FSRM J could be determined from this alternative eigenvalue potential with z i = z, we conclude that the relations in Eqs.(S22) hold for our original eigenvalue potential.
Eliminating R i and R * i in Eqs.(S19) and (S20) for G i , we arrive at the following system of equations for the diagonal elements of the resolvent matrix with Hence, to find the spectrum of the matrix J, one must first simultaneously solve Eqs.(S17), (S18) and (S23) for the diagonal resolvent elements G i .The spectral density can then be computed according to Eq. (S1).

S2. SELF-CONSISTENT EQUATIONS FOR THE SUPPORT OF THE SPECTRUM
In this section we show how Eqs.(S17), (S18) and (S23) simplify if we are only interested in knowing the boundary of the support of the bulk spectrum, and not the density of eigenvalues inside.We then derive a similar condition which determines the outlier eigenvalues.

A. Boundary of the bulk spectrum
There are two solutions to Eqs. (S17), (S18) and (S23), corresponding to Hence, the diagonal elements of the resolvent are functions of z only, and not of z.G i is therefore an analytic function and the spectral density is zero.Therefore, the solution for which Q i = 0 corresponds to the region of the complex plane outside the bulk spectrum.Conversely, if Q i ̸ = 0, the trace of the resolvent matrix is a non-analytic function as Eq.(S23) contains both G i and its complex conjugate.Therefore, the eigenvalue density is non-zero, and we are in the bulk region of the spectrum of J.
Consequently, to find the support of the bulk spectrum we seek the set of points λ bulk where the two solutions to Eqs. (S17), (S18) and (S23) meet.The solutions are distinguished by whether Q i is positive or zero, so we write Q i → δQ i for small positive δ > 0 and expand to leading order in δ, finding Both Eqs.(S26) and (S27) are valid outside the bulk region, including on the boundary, and Eq. ( S25) is valid only on the bulk boundary.Equation (S27) determines G i (z) outside the bulk region. Because We can also manipulate Eq. (S26) into an eigenvalue equation.By dividing Eq. (S26) by the positive quantity p ≡ 1 N i P i and defining B ′ i = P i /p, we obtain where 1 is an N × N matrix of ones.It is a known fact that if any one of a positive matrices elements increases, then so too does the PF eigenvalue of that matrix [97][98][99].Therefore, the RHS of Eq. ( S30) is an increasing function of p.
It has a minimum when p = 0 with value λ PF sG(z)G(z)/N and a maximum value of ∞.We can therefore combine Eqs.(S29) and (S30) into a single condition.Points z outside the bulk spectrum all satisfy with equality only on the boundary of the bulk spectrum.

B. Outliers
So far, we have derived conditions determining the diagonal elements of the resolvent matrix outside the bulk region, as well as the boundary of the bulk region itself.In our derivation of these conditions, the matrix u played no part and was effectively set to 0 [37,63].Outlier eigenvalues are isolated eigenvalues of the FSRM J which lie outside the bulk spectrum.Writing J as sum of its deterministic and random part as in Eq. (S4), the outlier eigenvalues λ where J 0 is equal to the random matrix J with the mean values u set to zero, outlier is outside the bulk region, it is not an eigenvalue of J 0 .Therefore, the matrix λ (k) outlier I − J 0 is invertible, and its inverse is the resolvent [G(z)] ij = δ ij G i (z), where G i (z) are the diagonal elements of the resolvent matrix of J.We note that S8 in this step we have assumed that the resolvent is self-averaging property of the FSRM assemble to equate a single realization of the resolvent of J with its disorder average.That is, we have assumed that Multiplying Eq. (S32) by det[G], we obtain Hence, to find outlier eigenvalues, we first solve Eq. (S27) for the diagonal elements of the resolvent matrix G i (z), which is valid for all z outside the bulk of the eigenvalue spectrum.Then, we substitute G i (z) into Eq.(S34) and enumerate all solutions.Finally, we must check that the solutions to Eq. (S34) actually lie outside the bulk spectrum, which we can do with Eq. (S31).

C. Summary
We now summarize and repeat the relevant equations derived in Sections S1 and S2.These equations self-consistently determine the boundary of the bulk spectrum and outlier eigenvalues of a FSRM with statistics given by Eqs. ( 1) in the main text, we repeat them here First, we find the diagonal elements of the resolvent matrix (valid for complex numbers z outside the bulk spectrum) from the following equation The boundary of the bulk spectrum comprises the set of solutions λ bulk to the following equation where [G(z)] ij = δ ij G i (z) is the diagonal resolvent matrix, the overline indicates element wise complex conjugation, and λ PF [M] is the Perron-Frobenius eigenvalue of the matrix M with positive entries.The outlier eigenvalues are given by the simultaneous solutions λ (k) outlier of the following equation and inequality For explicit examples of Eqs.(S36) to (S39), including a derivation of the standard elliptical law, see Section S6.

S3. APPROXIMATE SELF-CONSISTENT EQUATIONS FOR THE SUPPORT OF THE SPECTRUM A. Approximate statistics of an FSRM
In this section, we derive an explicit formula for the support of the eigenvalue spectrum of a general FSRM when there is a small amount of fine structure.Specifically, we suppose that the statistics of an FSRM J ij , given by Eqs.(4) in Eqs.(S45) and (S46) can both be written as eigenvalue equations, we write them both in the equivalent matrix-vector form where B is the right PF eigenvector of the positive matrix 1 N |G(λ bulk )| 2 σ 2 1 + ϵs (1) , and where C is an eigenvector of the matrix 1  N |G(z)| 2 σ 2 1 + ϵs (1) .Both vectors B and C correspond to an eigenvalue of 1.We note that Eq. ( S25) is essentially the same as Eq.(S48), but is written in terms of the left eigenvector Q.
To proceed with the derivation of the fine-structure corrections to the elliptical law and outlier eigenvalue, we suppose that the unknown variables G i , z, B i and C i can be expanded in powers of ϵ i , (S52) In the remainder of this section, we will expand Eqs.(S44), (S48) and (S49) to second order in ϵ.We will derive our modified elliptical law and outlier eigenvalue [Eqs.(10) and ( 13) in the main text] as the simultaneous solution to these expanded equations.
In the remainder of this section we stop writing the explicit z dependence for the resolvent, and will frequently use the abbreviations i /N and similarly for B i and C i .

B. Second order expansion of Eq. (S36)
On substituting Eqs.(S50) and (S51) into Eq.(S44), we expand to second order in ϵ and equate terms.We proceed order by order in ϵ, substituting the solution for the zeroth order equations in to the first order equations and so on.To order ϵ 0 we find Hence, the zeroth order approximation to the diagonal elements of the resolvent, G i , have no i dependence and G (0) i = G 0 .We may therefore write z 0 as Equating terms of order ϵ, we obtain By summing over the index i, we find We also find the following useful expression for G i , which we will use repeatedly throughout this section Equating terms of order ϵ 2 gives where we have used Eq.(S58) to go from the first to the second line.Again, summing over the index i, we find Finally, we compute z = z 0 + ϵz 1 + ϵ 2 z 2 and arrive at the following compact expression for G, valid to second order in ϵ Hence, to second order in ϵ, we can replace the N equations in Eq. (S44) for each diagonal element of the resolvent G i (z) with a single equation for the trace G(z).
C. Second order expansion of Eq. (S37) We follow the same procedure as in the preceeding section for expanding Eq. (S48) in small ϵ, by substituting Eqs.(S50) and (S52) into Eq.(S48) and equating term by term in powers of ϵ.We also note that, from Section S3 B, we know To zeroth order we have As the elements B i must all be positive, by the Perron-Frobenius theorem the only solution to Eq. (S62) requires that σ 2 |G 0 | 2 = 1 and B (0) i = B 0 for each i.Equating terms of order ϵ, we find .

(S63)
On summing over the index i, we see that G 1 = 0. Finally, the second order equation is Summing over the index i and substituting Eqs.(S58) and (S63) gives where we have written φ for the argument of the summed resolvent G = |G|e iφ .If we now compute the sum Hence, the values of the resolvent on the boundary of the bulk spectrum can be parametrized in terms of its complex argument φ as Finally, we remark that one could derive Eq. (S47) similarly to the expansion of Eq. (S45) that we have performed here.The expansion of Eq. (S47) gives the following condition on outlier eigenvalues where φ (k) is the argument of the outlier eigenvalue λ D. Second order expansion of Eq. (S38) Expanding Eq. (S49) is very similar to the expansion of Eq. (S48) detailed in Section S3 C. The main difference between Eq. (S48) and Eq.(S49) is that, in Eq. (S48), the elements B i must all be positive, and so the Perron-Frobenius theorem guarantees that there is only one solution.In contrast, the elements C i in Eq. (S49) do not need to be positive, and there are N possible eigenvectors C for each possible outlier eigenvalue λ outlier .However, we will see that only one of these solutions satisfies the further condition given in Eq. (S68), so in fact there is only one outlier eigenvalue.
Substituting Eqs.(S50) and (S53) into Eq.(S49), we have To zeroth order, we have One solution to Eq. (S70) is to take C = C 0 and uG 0 = 1.There are an additional N − 1 solutions corresponding to a diverging value of G 0 and the N − 1 linearly independent vectors C (0) such that i C (0) i = 0.However, for small ϵ we know from Section S3 C that |G i (λ bulk )| is O(ϵ 0 ).Therefore, as outlier eigenvalues must satisfy , only the first solution to Eq. (S70) has a chance of corresponding to an outlier eigenvalue when ϵ is small.For the remainder of this section we are only interested in the first solution.
Equating terms of order ϵ in the expansion of Eq. (S69), we find i , summing over the index i, we find G 1 = 0. Finally, equating terms of order ϵ 2 gives (1) i .
On substituting Eqs.(S58) and (S71) into Eq.(S72) and summing over the index i, we arrive at the following expression for the resolvent, evaluated at the outlier eigenvalue To summarize, in this section we have expanded Eqs.(S36) to (S39) to second order in the degree of fine structure (as measured by the small parameter ϵ).In doing so, we have found a set of self-consistent equations for determining the trace of the resolvent matrix [Eq.(S61)], the boundary of the bulk spectrum [Eq.(S66)], and the outlier eigenvalue [Eq.(S73)].We have also determined a criterion for testing whether a complex number is outside the bulk spectrum [Eq.(S68)].This criterion allows us to verify whether the outlier eigenvalue, as determined from Eq. (S73), corresponds to an outlier eigenvalue.

S4. MODIFIED ELLIPTICAL LAW
To find the support of the spectrum of a general FSRM exactly, we must simultaneously solve Eqs.(S36) to (S39) for the complex numbers λ bulk comprising the support of the bulk spectrum of J, and for the isolated outlier eigenvalues λ (k) outlier outside the bulk spectrum.In Section S3 we have seen that if the level of fine structure is small (as measured by the small parameter ϵ), then the support of the spectrum is determined from the simultaneous solution of Eqs.(S61), (S66), (S68) and (S73).In this section we will show how the solution of these equations leads to our modified elliptical law [Eq.(10) in the main text] and to the associated outlier eigenvalue [Eq.(13) in the main text].

A. Support of the bulk spectrum
To find the support of the bulk spectrum, we must simultaneously solve Eqs.(S61) and (S66) to second order in ϵ.Writing the resolvent in polar form G(z) = |G(z)|e iφ(z) , we substitute Eq. (S66) into Eq.(S61) to obtain To second order in ϵ, Eq. ( S74) is equivalent to the modified ellipse in Eq. ( 10) of the main text.To show this, we write λ bulk (φ) = x(φ) + iy(φ) and take the real and imaginary parts of Eq. ( S74) Note that, if we set ϵ = 0, then we recover the usual elliptical law To find the order ϵ 2 correction to the zeroth order ellipse, it is helpful to compute the following terms involving x 2 and y 2 x(φ) σ(1 + γ) By considering the product of Eqs.(S77), we find Substituting Eq. (S78) into Eqs.(S77) and adding up the resulting equations, we find Finally, we use the binomial identity (1 + Aϵ 2 /2) −2 = 1 − Aϵ 2 + O ϵ 4 to move the terms of order ϵ 2 on the left-hand side of Eq. (S79) to the denominator, giving us Eq. ( 10) from the main text In the main text, we drop the factors of ϵ in our statement of the modified elliptical law, absorbing them into the fine-structure parameters T, S and R.

Parametric form for the boundary of the bulk spectrum
In the main text we present the parametric equations Eqs.(11) as a practical way to plot the support of the bulk spectrum [Eq.(S80)] of the FSRM J.The proof is simple, taking we compute (to second order in ϵ) x a Hence, the expressions are equivalent to second order in ϵ.

B. Outlier eigenvalue
As discussed in Section S3 D, whilst a general FSRM can have multiple outlier eigenvalues, a FSRM with a small amount of fine structure only has a single outlier eigenvalue, much like in the usual elliptical law [52].The location of the outlier eigenvalue [Eq.( 13) in the main text] is obtained by substituting Eq. (S73) into Eq.(S61) and expanding to second order in ϵ, we arrive at Eq. ( 13) in the main text We note that Eq. (S83) is only valid for certain values of the model parameters.To check that λ outlier lies outside the bulk we must use Eq.(S68), the fine-structure correction to Eq. (S39).Noting that G(λ outlier ) is real, this condition is which is equivalent to Eq. ( 14) in the main text.

S5. EFFECT OF FINE STRUCTURE ON STABILITY
In Section III D, we make a number of conclusions about the effect of symmetry and antisymmetry in the structure of the statistics of the FSRM J. Central to our analysis is the following function which we claim is a matrix semi-norm.That is, the function ∥ • ∥ satisfies the triangle inequality, absolute homogeneity and positivity The properties Eqs.(S86) all follow from the observation that ∥A∥ = |A1|/ √ N 3 , where 1 is a vector of ones and | ˙| is the usual vector norm |a| = i a 2 i .We also note that ∥A∥ = 0 if only if all column sums of A are equal to zero, which is straightforward to verify from the definition.
B. Writing R, S, U , and V in terms of ρ1, ρ2, Ss, Sa, Us, Ua and T .
C. When fine structure is stabilizing and destabilizing Stability is determined from the rightmost eigenvalue of the spectrum of a FSRM J.The rightmost eigenvalue is either equal to the edge of the bulk spectrum λ edge , or is equal to the outlier eigenvalue λ outlier , for which the fine-structure corrections are given in Eqs. ( 12) and ( 13) in the main text.Written in terms of S s , S a , U s , U a , T, ρ 1 , ρ 2 , they are Recall that |γ| ≤ 1, and that for an outlier to emerge from the bulk we must have µ ≥ σ.Taking derivatives with respect to S a and U a , we have hence, as discussed in Section III D the main text, antisymmetric fine structure in the statistics of a FSRM is stabilizing.
The effect of symmetric fine structure on stability is determined by the remainder of the fine-structure corrections in Eq. (S91) if λ s edge , λ s outlier is positive/negative, then we say that symmetric fine structure in the statistics of the FSRM J is destabilizing/stabilizing. Noticing that Eqs.(S92) are quadratics in the variables √ T and √ S s , completing the square for the variables √ S s and √ U s respectively gives The first term in both expressions is always positive.Therefore, sufficient conditions for the positivity of λ s edge and λ s outlier can be found by analyzing when the second term in both expressions changes sign.As T, S s ≥ 0 and |γ| ≤ 1, the second term in the curly braces on the RHS of Eq. (S93) is positive The LHS of Eq. (S95) takes values between 0 and 1, and the RHS takes values between 3/4 and ∞.Therefore, if the LHS is smaller than 3/4, or if the RHS is greater than 1, then Eq. (S95) is automatically satisfied.Solving for ρ 1 and γ in these two cases, we find that symmetric fine structure in the statistics of a FSRM J is destabilizing if γ ≥ −1/3 or if ρ 1 ≥ − √ 3/2.In the main text, we only quote the first of these conditions.
The situation is more straightforward for the outlier eigenvalue.As T, U s ≥ 0, |γ| ≤ 1 and µ ≥ σ for an outlier to emerge to the right of the bulk, the second term in Eq. (S94) is positive provided The LHS of Eq. (S96) takes values between 0 and 1, and the RHS takes values between 1 and ∞.Therefore, λ s outlier is always positive, and we conclude that symmetric fine structure in the statistics of a FSRM is always destabilizing.
If we were to expand the statistics in small ϵ as we did to derive the fine-structure correction to the elliptical law and outlier eigenvalue, then we recover the same prediction as would be obtained by using the modified elliptical law directly.We first set To find the bulk spectrum and outlier eigenvalues to second order in ϵ, we need the eigenvalues of u and s to second order in ϵ.It is relatively straightforward (the derivations are very similar to those in Section S3) to show that, to second order in ϵ, these eigenvalues are where S and U are as defined in Eq. (S43).Substituting into Eqs.(S100) and (S102), and noting that if T = 0 then both R and V are also equal to 0, we recover the result that would be obtained from the modified elliptical law and outlier eigenvalue.

The elliptical law
An elliptical matrix ensemble is defined by the statistics These are the statistics of a FSRM with no fine structure.Therefore, it is clear that the fine-structure parameter T = 0, and we can use the exact solution detailed in Section S6 A, of course we could also set all fine-structure parameters to zero in the modified elliptic law in Eqs.(10) and ( 13) of the main texts, the result is the same.We need to know the PF eigenvalue of s/N and all eigenvalues of u/N , they are Only the k = 1 outlier of u/N can correspond to an outlier eigenvalue as it is the only eigenvalue of u/N which satisfies |λ k [u/N ]| > λ PF [s/N ].Substituting into Eqs.(S100) and (S102), we find expressions which agree with the modified elliptical law and outlier in Eqs.(10) and ( 13) of the main text.
2. The cascade model of Ref. [42] The cascade model of Ref. [42] is equivalent (albeit with a different parametrization) to the cascade model we present in Section V with the model parameter γ 1 set to zero.The statistics are As the matrix t is a constant matrix, the fine-structure parameter T = 0. Hence, we can again use the exact solution detailed in Section S6 A, for which we need the PF eigenvalue of s/N and all eigenvalues of u/N .In Section S6 B, we show that they are given by where Θ(x) = 1 if x > 0 and is 0 otherwise.Substituting into Eqs.(S100) and (S102) gives the exact support for the bulk spectrum, as well as the location of the outlier eigenvalues for the cascade model.In the case of the boundary of the bulk spectrum, the results here are the same as those obtained in [42,71].As far we are aware, the exact location of the outlier eigenvalues has not been presented before.

Directed networks
In section Section VI in the main text, and in Section S8 A, we claim that our fine-structure approximation for the spectrum of network with directed and undirected links [see Section VI of the main text] reduces to the exact solution when the underlying network contains only exclusively directed links.The FSRM of interest has statistics (obtained by setting P i↔j = 0 in Eqs.(33) of the main text) where P i←j is the probability of a link existing in the network going from node i to node j.As t ij = 0, clearly the fine-structure parameter T vanishes, and therefore we can use Eqs.(S100) and (S102) to find the support of the spectrum.Hence, we only need to know the eigenvalues of P i←j to determine the spectrum, or equivalently, we need to know the eigenvalues of the adjacency matrix A ij .The eigenvalues of undirected and directed networks are considered in Refs.[100,101].Provided the maximum degree of a node in the network is sufficiently small, the leading eigenvalue is equal to where we recall the definition of ρ from Eqs. (33) in the main text.All other eigenvalues are negligible (again, assuming the maximum degree is sufficiently small), hence, by Eq. (S99), there is only one outlier eigenvalue.Substituting into Eqs.(S100) and (S102), we find the bulk boundary and outlier eigenvalue Eqs. (S111) are equivalent to Eq. ( 71) in [72].

Circulant variances and correlations
We now turn to an example FSRM ensemble for which the matrix t contains fine structure, but for which the value of T is nonetheless equal to zero, we consider an FSRM for which the matrices t and s are circulant.That is, we suppose that t ij = t(i − j mod N ) and s ij = s(i − j mod N ).In this case the PF eigenvalue of the matrix s ij is simply the sum of any one of its rows 1 N j s ij .and (S36) with u, s, and t as in Eqs.(24) in the main text.We then expand the exact solution in powers of ϵ and verify that it agrees with the fine-structure correction we find in Section V of the main text.

Solving for the resolvent
To find the support of a FSRM with statistics as in Eqs.(S115), we must solve the self-consistent equations Eqs.(S29), (S34) and (S36).We will start by solving Eq. (S36) for the diagonal elements of the resolvent matrix G i (z).For large N , the sums in Eq. (S36) approach integrals over the continuous functions Eq. (S116) can be exactly solved for G(α, z) with the following manipulations.First, we differentiate Eq. (S116) with respect to α [noting that the distributional derivative of sign(α) is 2δ(α)], the result reads Hence, the resolvent G(α, z) can be found by inverting Eq. (S118) for A(z).Luckily, we do not need to do this to find the spectrum of J, as Eqs.(S29) and (S34) can be manipulated into expressions involving the function A(z) only with relative ease.To find the bulk boundary and any outlier eigenvalues, we must solve Eqs.(S29) and (S34) with the resolvent G(α, z) determined from Eq. (S118) and u and s as in Eqs.(24).In our solution of both equations, we will use the following useful information, which we prove in Section S6 B 5. For a matrix [A] ij = a + a 1 sign(i − j) and any diagonal matrix D, the following holds where the eigenvalues of A are given by where q k runs through the even integers if a > a 1 and runs through the odd integers otherwise.
Hence, we can sum Eqs.(S138) to get Using the binomial identity (1 + Aϵ 2 /2) −2 = 1 − Aϵ 2 + O ϵ 4 to move the terms of order ϵ 2 on the left-hand side of Eq. (S140) to the denominator, we arrive at an expression in the form of the modified elliptical law

Outliers
The exact kth outlier eigenvalues λ (k) outlier for the cascade model with statistics as described in Eqs.(24) of the main text are given by Eqs.(S126) and (S127).We repeat the results here with the small parameter ϵ included where the complex number A k is a solution to As we are assuming that ϵ is small, we see that q ℓ k , which is even if ϵµ 1 ≤ µ and is odd otherwise, must run through the even integers.In fact, q ℓ k must be equal to 0. For all other values of q ℓ k , the eigenvalue λ k [u] is of order ϵ and therefore violates the inequality . Explicitly, we see this by expanding λ k [u] [given by Eq. (S120)] to second order in ϵ, finding Hence, as λ PF [s] = O(1), we need only solve Eq. (S143) for A 1 , and not for each A k .
Assuming that the unknown quantity A 1 can be written A 1 = A 0 + ϵA 1 + ϵ 2 A 2 , we expand Eq. (S143) and equate powers of ϵ, finding On substituting A 1 = ϵγ 1 σ 2 /µ 2 into Eq.(S142) and expanding to second order in ϵ, we arrive at Eq. ( 13) in the main text C. fine-structure correction from the modified elliptical Law We will now reproduce Eq. (S141) using our modified elliptical law, thereby verifying its correctness.Computation of the zeroth order fine-structure parameters from Eqs. (S132) is straightforward, the zeroth order parameters are simply µ, σ 2 and γ.
where the covariance of K ij and K ji is zero for the directed distribution because one of K ij or K ji must be zero if there is a directed link j ← i or i ← j.We have excluded the possibility of a link i ← j and j ← i by definition of the network.
To derive the statistics of A • K over the whole ensemble, rather than just over the distribution of non-zero elements, we show that the statistics of a random matrix constructed according to Eq. (S150) and Eqs.(S151) have the same statistics as a fully connected FSRM J with statistics [Eqs.(32) in the main text] To show this, we consider the eigenvalue potential associated with the matrix A • K with statistics described in Eq. (S150) and Eqs.(S151), it is given by exp Focusing only on the parts which contain factors of A ij K ij for the disorder average, we can write the average over realizations of A • K in terms of the averages ⟨ Assuming that the average node degree p is large, we can expand the exponentials inside the averages in powers of 1/p and average the terms in the expansion.We have exp (q i p j + q i p j ) 2 + Γσ 2 (p i q j + p i q j )(p j q i + p j q i ), (S155) for the disorder average over reciprocated links and exp for the disorder average over directed links.Substituting these results back into Eq.(S154) and using the relation 1 + A/p ≈ exp(A/p), which is valid for large p, we can re-exponentiate the averaged statistics to obtain exp {P i↔j + P i←j }(q i q j + q i q j ) 2 + Γσ 2 P i↔j (q i q j + q i q j )(q j q i + q j q i )   . (S157) Finally, we recognize that one would obtain the exact same expression for the eigenvalue potential from a fully connected FSRM with statistics as in Eq. (S152).Note that it would also be straightforward to modify this argument to allow for the matrix K to be an FSRM, rather than an elliptical matrix.
Test of the fine-structure correction to the outlier of a directed network [Eq.(36)].Parameters for all curves are µ = 1.5, σ = 1.0, γ = 0.5, r = 0.4 and p = 0.1.The parameters p in , p out and p ↔ are set to 1 or 0 as indicated.Curves are generated from Eq. ( 13) in the main text and markers are the maximum eigenvalue of random matrices constructed according the description in Section S8 A with the statistics of non-zero elements given by Eqs.(S151) and on a network with degree distributions given by Eqs.(S164).Markers are the average of 10 realizations of the matrix with N = 4000.

Test of the outlier condition with dichotomous degree distributions
To test our fine-structure correction to the spectrum of a directed network we use a network with dichotomous degree distributions.Specifically, the network has degree distributions (S164) We can compute the relevant statistics for the fine-structure correction, they are On substituting into Eqs.( 10) and ( 13) in the main text, we find the fine-structure correction to the spectrum of random matrices with a dichotomous degree distribution.
In the main text, we claim that if the underlying network is either completely directed (P i←j = 0) or completely undirected (P i↔j = 0), then our results reduce to known cases [65,72].In the following two sections we demonstrate this.
main text] We emphasize that the statistics given are those of a random matrix constructed according to the niche model for a fixed instance of the matrix E. However, as we will show, the elliptical and fine-structure parameters computed from Eqs. (S174) ultimately depend on statistical properties of E, not on any particular instance.
To find the fine-structure correction, we first plug the statistics Eqs.(S174) into Eqs.(S41) to find expressions for the elliptical and fine-structure parameters in terms of µ L , µ U , σ L , σ U , Γ, and the matrix E. We will then show that these parameters do not depend on the specific instance of E that we use, but on the statistics of the elements of E, which are ultimately functions of the average niche radius d only.By approximating these statistics, we finally find expressions for the elliptical and fine-structure parameters in terms of the model parameters µ L , µ U , σ L , σ U , Γ, and d.

fine-structure parameters
We can find the elliptical parameters in a straightforward manner by substituting the statistics Eqs.(S174) into Eqs.(S41) [Eqs.(5) in the main text], which gives where d = N −2 ij E ij .We can also compute the fine-structure parameters similarly, by plugging the statistics into Eq.(S43) [Eqs.(6) in the main text].For example, we can compute S as follows Using the identity , we arrive at where we have defined the following quantities (which are statistical properties of the matrix E), (S178)

S33
The remaining fine-structure parameters can be calculated similarly, and eventually we end up with Whilst we have derived the elliptical and fine-structure parameters for a specific realization of the matrix E, Eqs.(S179) reveal that our correction only depends on the network structure through d, h 2 , and f , which, for large N , are properties of the ensemble from which E is drawn, and not of any particular instance of the network.That is, for large N , we have Hence, d, h 2 , and f are non-random functions of d only, as that is the only parameter used to construct realizations of the network.If we can determine these functions, then we have the fine-structure correction to the spectrum of a random matrix constructed according to the niche model.

Approximating the statistics of the network
First, we observe that d, h 2 , and f can all be written in terms of the following quantities, which are simpler We also repeat the definition of E ij from Section VIII A in the main text with the dependence of E ij on the niche position η i , range d i , and center c i (which are all random variables, also defined in Section VIII A of the main text) made explicit For the definitions of the random variables η i , d i , and c i , see Section VIII A of the main text.
It turns out that E in i is much easier to compute than E out i , so we start with E in i .As η j is a uniformly distributed random variable on the interval [0, 1], and as both of c i − d i /2 and c i + d i /2 are, by definition, confined to the interval [0, 1], the sum j E ij /N is simply equal to the length of the interval [c i − d i /2, c i + d i /2] when N is large.That is, we have As an immediate consequence, we also have d = i d i /N = d.The expression for E in i E is compared to numerical simulations in Fig. S4(b).
On the other hand, E out i is more difficult to evaluate, so we resort to an approximation.We write where and where δ i is the error in the approximation.The error depends on the higher central moments of η i , d i , and c i .Whilst δ i itself is generally not small [see Fig. S4(c)], we will see that the sum i δ i /N vanishes for large N .
The average of η i , d i , and c i over realizations of E are straightforward to compute using their definitions [see Section VIII A of the main text].The variables η i are uniformly distributed random variables on the interval [0, 1], and they are ordered so that η 1 < η 2 < • • • < η N .The fact that these variables are ordered implies that ⟨η i ⟩ E = i/N and that the variance of η i is sub-leading in N [102,103].That is, each η i is equal to its average, up to a term sub-leading in N .The niche range of each species is given by d i = Xη i , where X is a beta distributed random variable drawn independently for each species with PDF p X (x) = 1 2r − 1 (1 − x) 1 2r −2 .The average value of d i over realizations of E is therefore The distribution of the niche center of species i, c i , is i , where c (1) i is a uniform random variable on the interval [d i /2, η i ] and c (2) i is uniform on the interval [d i /2, 1 − d i /2].To compute ⟨c i ⟩ E , we have to average over both the distribution of d i and of c i , this can be done using the law of total expectation.Recalling d i = Xη i , we have where H(x) = max(0, x).
By inspection, we can simplify ⟨c i ⟩ E with the following approximation The expressions in Eqs.(S188) and (S189) are compared with one another in Fig. S4(a).Overall, we find that there is a maximum error of ∼ 5%.The two expression coincide with each other if any of the following conditions is fulfilled: η i < 2/3, η i = 1, d = 1/2, or d = 0. We will use the approximation in Eq. (S189) to compute the fine-structure parameters.
Substituting for ⟨η i ⟩ E , ⟨d i ⟩ E , and our approximation to ⟨c i ⟩ E [⟨η i ⟩ E = i/N , Eqs. (S186) and (S189)] for the summand in Eq. (S184), we arrive at  To determine E out i from E approx ij , we note that, in the large N limit, the summations in Eqs.(S181) can be approximated by integrals over the variable α ≡ i/N from 0 to 1. Defining E approx (i/N, j/N ) = E approx ij , and similar for E ij , E in i , E out i , and the error term δ i , we have Noting that Let us now use E in (α) and E out (α) [Eqs.(S183) and (S191)] to find expressions for h 2 and f in terms of the average niche range d only.Neglecting terms proportional to δ(α), we have On carrying out the integrals, we find (S193) To obtain Eqs.(52) in the main text from Eq. (S193), we take the order [2/2] Padé approximant of h 2 around d = 0, which incurs a maximum relative error of approximately 0.004 when d = 1/2, or approximately 5% of the value of h 2 .
On substituting h 2 and f into Eqs.(S179), we have successfully found an explicit fine-structure correction for FSRMs constructed according to the niche model, in terms of only the model parameters µ L , µ U , σ L , σ U , Γ and d.We compare our approximation h 2 and f in Fig. 5(b) and (c) in the main text.The fine-structure correction to the outlier eigenvalue is verified in Fig. 4(b) and the correction to the bulk of the spectrum is verified in Fig. S3.

FIG. 1 .
FIG. 1. Support of the spectrum of FSRMs constructed according to the cascade model.(u), (s) and (t): heatmaps of the matrices u, s and t as defined in Eqs.(24).(a): The bulk and outliers of a single 8000 × 8000 cascade model FSRM with µ = 2.5, µ1 = 3.5, σ = 0.5, σ1 = 0.1, γ = 0, γ1 = 1.The blue solid lines and circles are the exact analytical prediction for the boundary of the bulk and outliers [given in Section S6 B of the SM].(b): The leading eigenvalue λ1[J] of a cascade model FSRM with γ = 0.2 and all other model parameters as indicated.We also divide the leading eigenvalue by λ1[J0], the leading eigenvalue of a corresponding elliptical random matrix with statistics µ, σ and γ to isolate the effect of fine structure.For the blue, orange and red curves (squares and triangles), λ1[J0] = σ(1 + γ) because the bulk edge is the leading eigenvalue.For the green curve (circles), λ1[J0] = µ + γσ 2 /µ because the outlier is the leading eigenvalue.Solid lines are the exact solution [Section S6 B] and dashed lines are the fine-structure corrections to the leading eigenvalue given in Eqs.(27) and(28).Each marker is the result of diagonalizing a single 8000 × 8000 instance of the matrix.

FIG. 2 .
FIG.2.The spectrum of a directed network with statistics as in Eqs.(32).(u),(s), (t): Heatmaps of the statistics of J = A • K for a single realization of the network A with N = 40, coloured squares indicate the presence of a link and white squares indicate no link.(a) Eigenvalues of a single 8000 × 8000 realization of J. Blue solid line and circle are the fine-structure correction to the eigenvalue spectrum of a directed network with µ = 2, σ = 1, Γ = 0.7, orange dashed line and triangle show the elliptical law with no fine-structure correction.Half of the nodes have expected degrees (k ↔ , k in , k out ) = (325, 300, 300) and the other half have expected degrees(75,100,100).(b) Ratio of the bulk edge eigenvalue of J to σ(1 + Γr), the bulk edge without the fine-structure correction.The top three curves when r = 0 correspond to random matrices with positive degree correlations ρ, the bottom three curves correspond to ensembles with the same degree distributions with the links rearranged so that the degree correlation is negative, demonstrating that increasing ρ leads to destabilization.When r = 1, the degree heterogeneity h 2 controls the size of the fine-structure correction, and it stabilizes/destabilizes depending on the value of Γ.The network has the same dichotomous structure as that used in panel (a) and µ = 0.5, σ = 1 and Γ as indicated.Markers are the result of diagonalizing a single 8000 × 8000 realization of the matrix.We verify Eq. (36), the fine-structure correction to the outlier eigenvalue, in Fig.S2of the SM.

FIG. 3 .
FIG.3.Spectrum of a FSRM with statistics as in Eqs.(38).(u), (s) and (t): Heatmaps of the matrices u, s and t in Eqs.(38).(a) Eigenvalues of a single 8000 × 8000 realization of the matrix.Blue solid line and circle are the fine-structure correction to the eigenvalue spectrum [Eqs.(10) and (13)], orange dashed line and triangle show the elliptical law with no fine-structure correction.Parameters are µ = 1.1, σ = 1, Γ = 1 and k = 3.(b) Blue curves (squares): ratio of the fine-structure prediction for the bulk edge of J to the bulk edge of an elliptical matrix with no fine structure.Orange curves (triangles): as for the blue curves but for the outlier eigenvalue.Solid curves indicate the fine-structure correction [Eqs.(10) and (13)], dashed lines indicate the small-k correction [Eqs.(44)].Blue (squares) markers are the result of diagonalizing a single 8000×8000 realization of the matrix, orange (triangles) markers are the result of diagonalizing 100 instances of 4000 × 4000 realizations of the matrix and averaging.Parameters are µ = 1.1, σ = 1 and Γ = 1, with k as indicated.

FIG. 4 .
FIG. 4. Spectrum of an FSRM constructed according to the niche model [Eqs.(45)].(u), (s), (t): Heatmaps of the statistics u, s, t with N = 50.Colour indicates the strength and presence of a link, and white indicates no link in the network.(a) Predictions from Eqs. (10) and (13) are shown in blue (solid line and circle), and the elliptical approximation to the spectrum is shown in orange (dashed line and triangle).Grey markers are the result of exact diagonalization of a 12000 × 12000 realization of the matrix with µL = 3, µU = 5, σL = 0.6, σU = 1.4,γ = −0.6 and d = 0.05.(b) The ratio of the maximum eigenvalue of a niche model matrix to that of an elliptical matrix with no fine-structure correction.Lines indicate the result of the finestructure correction [Eqs.(10) and (13)], markers are computed from the average of 10 realizations of 4000 × 4000 matrices.The remaining system parameters are µL + µU = 8, σL = 1.2, σU = 0.8, d = 0.1, and the parameters Γ and µL − µU are as indicated in the figure.

FIG. 5 .
FIG. 5. Structure and statistics of the matrix E in the niche model.(a) A single realization of the matrix E with N = 80 and d = 0.15, constructed according to rules (1) to (4) in Section VIII A. Panel (b): Solid line shows the degree heterogeneity (h 2 ) as given in Eqs.(52), symbols are from simulations.Panel (c): Solid line shows the quantity f in Eqs.(47), and the symbols are again from computer generated networks.Simulation data in (b) and (c) are from a single realization of E with N = 2000.
) Now, because the vector B ′ has all positive entries [see Eqs.(S10)], it must be the right PF eigenvector of the positive matrix with elements [G(z)G(z)(p1 + s)] ij /N = |G i (z)| 2 (p + s ij )/N for some positive number p, with eigenvalue 1. Writing λ PF [M] for the PF eigenvalue of a matrix M with all entries positive, we can write Eqs.(S25) and (S26) as

45 FIG
FIG. S4.(a) Comparison of the exact (dashed lines) expression for ⟨ci⟩ E [Eq.(S188)] and its approximation (solid lines) [Eq.(S189)].The functions coincide when d = 0 or d = 1/2 and only differ in a small range for intermediate values of d, shown zoomed in as indicated in (a1).Markers are the result of averaging 10000 realizations of the niche center ci.(b) and (c) Comparison of numerical (faded wiggley lines) values of E in i E and E out i E [Eqs.(S181), averaged over 10 realizations] and their approximations [(a) solid straight lines, (b) solid piecewise straight lines] as given in Eqs.(S183) and (S191).
Writing R, S, U , and V in terms of ρ 1 , ρ 2 , S s , S a , U s , U a and T .contains additional details of the calculations presented in the main text, as well as additional figures and examples of FSRM eigenvalue spectra.