Stability-instability transition in tripartite merged ecological networks

Although ecological networks are typically constructed based on a single type of interaction, e.g. trophic interactions in a food web, a more complete picture of ecosystem composition and functioning arises from merging networks of multiple interaction types. In this work, we consider tripartite networks constructed by merging two bipartite networks, one mutualistic and one antagonistic. Taking the interactions within each sub-network to be distributed randomly, we consider the stability of the dynamics of the network based on the spectrum of its community matrix. In the asymptotic limit of a large number of species, we show that the spectrum undergoes an eigenvalue phase transition, which leads to an abrupt destabilisation of the network as the ratio of mutualists to antagonists is increased. We also derive results that show how this transition is manifest in networks of finite size, as well as when disorder is introduced in the segregation of the two interaction types. Our random-matrix results will serve as a baseline for understanding the behaviour of merged networks with more realistic structures and/or more detailed dynamics.


Introduction
A central goal of community ecology is to identify the mechanisms that maintain biodiversity in natural communities. One way to approach this problem is through the lens of ecological networks -abstract representations of interactions between taxa in a community (Montoya et al. 2006;Landi et al. 2018;Delmas et al. 2019). Classically, studies of ecological networks have tended to focus on a single type of interaction, e.g. trophic interactions in a food web (Pimm 1982;Dunne et al. 2002;Dunne 2012) or mutualistic interactions in a plant-pollinator network (Bascompte and Jordano 2013;Valdovinos 2019). Over the last decade or so, however, there has been a growing recognition that further progress requires the construction and study of networks that contain multiple interaction types (Ings et al. 2009;Fontaine et al. 2011;García-Callejas et al. 2018) and even combining ecological with social interactions (Felipe-Lucia et al. 2021).
In this paper we consider the stability of ecological networks based on the tripartite structure of Fig. 1. These networks consist of three guilds (nominally plants, herbivores and mutualists) interacting through two distinct sets of interactions, one antagonistic and the other mutualistic. This tripartite structure represents a minimal example of a merged network (Fontaine et al. 2011;Sauve et al. 2014), and reflects the largescale structure seen in empirical networks such as those discussed by Melián et al. (2009), Pocock et al. (2012), Kéfi et al. (2016), Sauve et al. (2016), Miller et al. (2021), Laha et al. (2022). In the framework of multi-layer networks (Pilosof et al. 2017), Fig. 1 could be viewed as a two-layer multiplex network with antagonists in one layer, mutualists in the other and plants common to both. Two recent surveys analysed structural features and robustness of a number of empirical tripartite networks from the literature, including those with mutualist and antagonist partitions as considered here; Timóteo et al. (2021) found that the importance of a species is positively correlated between the two bipartite subnetworks, and Domínguez-García and Kéfi (2021) found that the robustness (with respect to plant losses) of tripartite mutualist-antagonist networks could be understood in terms of the robustnesses of the two bipartite networks composing them.
The aim of the present contribution is to provide an analysis of the dynamical stability of this network structure when the subnetwork interactions are described by random matrices. The use of random matrices to shed light on ecosystem stability has a long history (May 1972(May , 1974 and acts as a baseline scenario against which more detailed and ecologically-motivated studies can be compared (Allesina and Tang 2012). The stability of networks comprising a mixture of antagonistic and mututalistic interactions has previously been considered from a random-matrix viewpoint (Mougi and Kondoh 2012;Suweis et al. 2014;Sellman et al. 2016). However, these studies have been of (stochastically) homogeneous models, i.e. without the tripartite structure of Fig. 1. The persistence and stability of an ecosystem with this structure was considered in Sauve et al. (2014) using numerical simulations for small networks. In contrast, the focus of our work is on analytic results for ecosystems consisting of a large number of species.
Specifically we investigate the stability of a community matrix with block-structure that reflects Fig. 1 and with random interactions within the blocks. Using the results of  Marčenko and Pastur (1967) and Benaych-Georges and Nadakuditi (2011) we give an account of key spectral features of such matrices in the limit of large network size. We show that the spectrum consists of a bulk component and a pair of eigenvalues of large magnitude. The properties of these latter split the behaviour of the model into two distinct phases. In the antagonist-dominated phase, the large eigenvalues are complex and the stability properties of the model are determined by the bulk. In the mutualistdominated phase, the large eigenvalues are real and serve to destabilise the system. The transition between these two phases represents an eigenvalue phase transition (Baik et al. 2005) and can be driven, for example, by an increase in the relative fraction of mutualists in the ecosystem. Whilst this phase transition strictly takes place in the asymptotic limit, we also derive an expression for the stability-determining eigenvalue valid at finite system size. Finally, we also consider a scenario in which we introduce a degree of disorder into the network of Fig. 1 by exchanging the type of a certain random fraction of interactions. We show that the phase transition can survive in the presence of disorder, but vanishes if the disorder is too strong.

Tripartite network model
We consider a system of N P plants (or producers), N M mutualists and N M herbivores such that the total number of consumers (mutualists and herbivores) is N C = N M +N H . We then define s = N C /N P as the ratio of consumers to plants and r = N M /N C as the ratio of mutualists to consumers. We assume that, close to equilibrium, the dynamics of the populations is described by community matrix K = −D + A, where D describes intraspecific competition and A describes interspecific interactions. As is conventional, we assume D to have elements D i j = dδ i j with δ i j the Kronecker delta symbol and d the competition strength. We then choose K to reflect the tripartite structure of the network in Fig. 1. In the first instance, we take A to have the following block structure in which the species have been ordered plants first, then mutualists and finally herbivores. The block M is of dimension N P × N M and describes the plant-mutualist interactions; block H is of dimension N P × N H and describes plant-herbivore interactions. We take all matrix elements of M and H to be positive such that the nature of the interactions is encoded in the explicit signs shown in Eq. (1). Parameters μ M and μ H are included to describe asymmetries in the two directions of each interaction type. The forefactor N −1/2 P is included for mathematical convenience. We note that a similar matrix block structure was considered by Johnson et al. (2014) but with sign assignments appropriate to exclusively trophic interactions.
This overall structure is then supplemented by a random model for the matrix elements of the blocks. Considering the mutualists, we set M i j = e M i j m i j where the interaction coefficients m i j are treated as independent identical random variables with mean m and variance σ 2 m , and where the factors e M i j describe the links of the plant-mutualist sub-network: species that interact have e M i j = 1; those that do not have e M i j = 0. The mutualist subnetwork connectance is therefore We assume a random structure for the plant-mutualist network with e M i j chosen from a Bernoulli distribution with probability P(e M i j = 1) = c M for all pairs i, j. The first and second cumulants of the matrix elements of M are thus An analogous set-up and notation is employed for the matrix elements of H .

Eigenvalue spectrum
The stability of the equilibrium described by community matrix K is determined by its N C + N P eigenvalues, , which obey K v = v. Writing eigenvector v in terms of three vectors, one for each component of the system, v = (v P , v M , v H ) and taking the block structure of K into account we arrive at the set of three equations Eliminating v M and v H from the first of these, we obtain the eigenvalue equation in which Accordingly, each positive eigenvalue λ > 0 contributes two real eigenvalues to the spectrum of K , whereas a negative eigenvalue λ < 0 contributes an imaginary pair to the spectrum of K . The number of zero eigenvalues of K depends on s, the ratio of consumers to plants. For s > 1, i.e. more consumers than plants, matrix W is of full rank and thus has N P non-zero eigenvalues λ = 0, equating to 2N P non-zero eigenvalues for K . Matrix K then has an additional N C − N P zero eigenvalues. On the other hand, if s < 1, W is rank deficient and has only N C non-zero eigenvalues.
Matrix K then has a total of N P − N C zero eigenvalues. The main features in the spectrum of W can be appreciated from numerical diagonalisation. For concreteness, we draw the nonzero elements of M from a half-normal distribution with probability density function The standard deviation of this distribution is σ m = m √ π/2 − 1. The nonzero elements of H are generated in a similar fashion but using parameter h. Figure 2 shows the numerical spectrum of instances of matrix W across the range of values of the mutualist ratio r with N P and N C fixed. For s > 1 [Figs. 2a and b] the spectrum clearly separates into a compact "bulk" spectrum located around λ = 0, plus up to two large eigenvalues situated outside the bulk with magnitudes strongly dependent on r . The situation for s < 1 [Figs. 2c and d] is similar, except that here the single bulk is replaced by two disjoint lobes with a further collection of eigenvalues of zero. The large eigenvalues are equally apparent in this case.

Eigenvalue phase transition
We now consider the spectrum of W analytically and show that the largest eigenvalue undergoes an eigenvalue phase transition in the N P → ∞ limit with the ratios s and r held fixed. We begin by writing W = B + P with where matrix X has dimensions N P × N C and elements and with T = diag τ 1 , τ 2 , . . . τ N C with The second matrix in the decomposition of W reads and The elements of matrix X are independent random variables each with zero mean and unit variance. Moments beyond the second play no role in our results in the large-N P limit, and thus we take all X i j to be distributed identically. Turning to Eq. (11), we observe that the second term is a summation over a large number of independent variables, and thus gives a contribution of the order N 1/2 C for large N C . Overall, then, this term scales like N −1/2 P and is therefore negligible in comparison with the first term which scales like θ ∼ 1. Thus, for P, we obtain the approximate rank-1 form with J the N P × N P matrix of ones. We first consider the matrix B by itself. Let its ordered eigenvalues be λ 1 (B) ≥ · · · ≥ λ N P (B) and define the probability measure μ B as the limiting empirical eigenvalue Marčenko and Pastur (1967) [see also Silverstein and Bai (1995)], the Stieltjes transform of μ B is given by where ν(τ ) is the N P → ∞ probability distribution function of the values τ 1 , . . . , τ N C . For the case in hand, this reads and thus we obtain and its inverse Again from Marčenko and Pastur (1967), the edges of connected components of spectrum can be obtained as In particular, this gives us b + and b − the upper and lower edges of the bulk part of the spectrum, i.e. the supremum and infimum of the support of μ B . Figure 2 shows the boundaries obtained from Eq. (19) superimposed on the numerical spectra. The boundaries, including b ± , give good approximation to the numerical bulk edges across the entire range of r . Now consider the complete matrix, W = B+P. Applying Theorem 2.1 of Benaych-Georges and Nadakuditi (2011) to the rank-1 perturbation of Eq. (14), we have that in the N P → ∞ limit the uppermost eigenvalue of W almost surely converges as Given that z μ B is a monotonically decreasing function at the spectrum edge and that N P is large, this becomes Thus, in the N P → ∞ limit, the largest eigenvalue of W undergoes an abrupt eigenvalue phase transition (Baik et al. 2005) from a value b + ∼ (N P ) 0 below the transition to a value λ 1 = N P θ that is "macroscopic" in the sense that it scales with the size of the ecosystem. This transition occurs at N P θ = b + but given the different scaling of the two sides of this equation, the transition point is effectively given by θ = 0 in the asymptotic limit. Thus, we see that the phase transition occurs when r = r with the critical mutualist ratio The significance of this is that eigenvalue λ 1 (W) determines the stability of community matrix K and is thus an order parameter for a phase transition in the stability of the model. We note that the quantity N P θ is the row sum of matrix W, as previously discussed by Allesina and Tang (2012) in the context of isolated eigenvalues of mutualist matrices. Considering the lowest eigenvalue, we similarly find This also exhibits a phase transition but this is of little significance from a point of view of the stability. Figure 2 shows the macroscopic value N P θ superimposed on the numerical spectra, where it is seen to be close to numerically-obtained large eigenvalues when they lie significantly outside the bulk. Deviations from this good agreement occur in the region close to where N P θ and b ± cross, and this is a consequence of the finite value of N P in this figure.
Finally we note that there is a secondary, less dramatic stability transition implicit in the above results and visible in Fig. 2. This transition occurs at the point b + = 0. Below this point, the spectrum of W is completely negative [see for example the lefthand edge of Fig. 2a] and thus √ μ M λ will be imaginary for all eigenvalues. The stability of community matrix K will, from Eq. (6), then be given purely by the intraspecific competition term. Above this point, √ μ M λ will be real for some λ and this will then start to reduce the stability. low r to the macroscopic disjoint value at high r . In this section we derive an expression for λ 1 that captures this behaviour at large but finite N P . In doing so, we consider a slightly more general model, namely where M 1 and M 2 are random blocks with relationship left unspecified, and similarly for H 1 and H 2 . This generalisation allows us to discuss the effect of disorder in the next section. Clearly the model of the previous section is recovered by setting With the interaction matrix of Eq. (23), the N P × N P matrix equivalent of Eq. (5) is An approximate account of the spectrum of W can be obtained by considering this matrix in the eigenbasis of its ensemble average W . Doing so allows us to identify the scaling properties of different parts of the matrix and derive an approximate expression valid for N P 1. Details of this calculation are described in the Appendix with the results as follows. The spectrum of W is seen to approximately contain N P − 2 eigenvalues at . This set of eigenvalues is the approximate representation of the bulk in this calculation. More importantly, we also obtain two non-trivial eigenvalues given by in which, similar to Eq. (12), we have and where In the asymptotic limit, Eq. (26) gives λ + → N P θ , which recovers the macroscopic eigenvalues obtained previously, and λ − → φ which then becomes part of the bulk.  Fig. 2. Results for a range of different values of plant number N P are shown. Clear agreement between numerical and analytic results is observed, with the degree of agreement increasing with N P as expected. The only significant deviation between the two at large N P occurs at very small values of r , where λ + overestimates λ 1 . In this regime the scaling arguments leading to Eq. (26) break down.
One additional feature that is revealed by this finite-size analysis is the role of correlations between the strengths of the two directions of the interactions  ] c = 0 and Ω 2 vanishes. The prediction for λ 1 in this case becomes the piecewise linear function λ 1 = max [φ, φ + N P θ ], similar to that found in Sect. 4 in the N P → ∞ limit. Thus the "avoided crossing" that occurs in the correlated case at large but finite N P gives way to an actual crossing when the degree of correlation is zero.

Interaction disorder
As defined by their interactions with plants, the consumers in Fig. 1 are either 100% mutualist or 100% herbivore. In this section we look what happens when we move away from this perfectly ordered scenario and swap the types of a random selection of the species interactions. We consider the same tripartite network as before, but with a probability f M we swap the (++) signs of the mutualistic interactions to the (+−) signs of an antagonistic one. For the antagonistic interactions, we do the opposite with a probability f H . The end result is that the animal species no longer act with well defined roles, but differently across their connected plant species. Mathematically, this is incorporated into the framework of Sect. 5 by selecting the matrix elements of Eq. (23) to be Here, m i j and h i j are random variables as before, m i j and h i j are further independent random variables chosen from the same distributions, and f Using the results of the previous section, we see that, in the asymptotic limit, the largest eigenvalue of the matrix W is λ 1 = max[0, N P θ ] with θ of Eq. (27) in terms of the ensemble averages Results are shown in Fig. 4 as phase diagrams in the plane defined by the mutualist ratio r and the fraction of swapped interactions f M = f H = f . When there is exact symmetry between the interactions strength, i.e. m = h and γ = 1 (Fig. 4 middle), the behaviour observed in the fully coherent case, f = 0, is preserved at finite f , up to a fraction f = 1/2. The transition still occurs at r = 1/2 and the only change is that the size of order parameter in mutualist phase is diminished.
At the point f = 1/2 the dominant roles of the "antagonistic" and "mutualistic" species swap and the two phases reverse accordingly. Away from this exact symmetry we have more complicated behaviour. For m < h (assuming γ = 1) the mutualist phase occupies a diminished area of the phase diagram. For values of f close to 1/2 the system remains in the antagonist phase across the whole range of r and no instability transition takes place. In contrast, for m > h, it is the antagonist phase that is diminished and for f close to 1/2 the system remains in the mutualist phase for all r . Indeed, specifying for simplicity the case of c M = c H , f M = f H = f and γ = 1 (as in Fig. 4), we find the value of r for which θ = 0 to be For a phase transition to occur, we require that 0 < r < 1, which gives as the condition on the interaction and disorder strengths for the existence of a phase transition.

Discussion
We have studied the stability of a community matrix with random elements organised according to the tripartite structure of Fig. 1. This structure can be viewed as the mergence of two bipartite networks with changes in the mutualist ratio r interpolating between them. At r = 0 we have a bipartite predator-prey model, whose community matrix has purely imaginary eigenvalues (Johnson et al. 2014); at r = 1 we have a bipartite mutualist model, whose stability is determined by a single large real eigenvalue (Allesina and Tang 2012). We have shown here that the emergence of this macroscopic eigenvalue as a function of r is abrupt in the asymptotic limit, and that this represents an eigenvalue phase transition that takes place when r reaches the critical value r . The behaviour of the model is thus split into two distinct regimes. For r < r , we obtain an "antagonistic phase" in which the community matrix K can be stabilised by a small (i.e. order N 0 P ) value of the intraspecific competition d. In this regime, changing the number of mutualists by a small amount does not appreciably affect the stability of the system. In contrast, for r > r we obtain a "mutualist phase" in which stability is governed by the macroscopic eigenvalue such that the community matrix requires a large intraspecific competition, d ∼ N 1/2 P , to stabilise it. Moreover, in this phase, small changes in the mutualist number lead to large changes in λ 1 and hence dramatic changes in the stability of the system. Whilst we have framed this discussion in terms of the behaviour as a function of the ratio r , the phase transition can also be driven by changes in other parameters, e.g. mean interaction strengths or connectance. For example, setting θ = 0 for the ordered model of Sect. 2 shows that, with all other parameters fixed, the phase transition will occur when the ratio of the interaction strengths reaches a value A large positive eigenvalue of the community matrix is associated with the growth of populations away from the equilibrium values. The rapidity of growth seen here is a consequence of the positive feedback of "mutual benefaction" associated with the mutualistic interactions (Scheuring 1992). It should be born in mind, however, that the community matrix represents a linearisation of a more complex/detailed dynamical model. Thus, instability should not be interpreted as an uncontrolled growth, but rather as the indication of the shift of the ecosystem away from its equilibrium to a different one, the properties of which are outside the scope of the original model. This new equilibrium may well include fewer species than in the starting community. The scaling in Eq. (1) is chosen such that the bulk spectrum of K converges in the N P → ∞ limit. This follows from the convergence of the spectrum of W with bulk edges b ± ∼ (N P ) 0 . In the spectrum of K , this leaves the scaling of the macroscopic eigenvalue as = −d + √ μ M N P θ ∼ √ N P . A slightly different choice would be to replace the N −1/2 P forefactor out the front of Eq. (1) with N −1 P . This would make extreme eigenvalue of K converge to = −d + √ μ M θ ∼ N P 0 in the limit, whilst the bulk would shrink to a point at = −d. Given the different scalings of the two parts of the spectrum, we might also consider a slightly different model in which this forefactor in Eq. (1) is omitted and we instead choose to scale the matrix-element distribution such that m ∼ N −1 P and σ m ∼ N −1/2 P and similarly for the herbivores. This situation would then be similar to the scaling in e.g. Bunin (2017), Galla (2018). This choice invalidates the assumptions upon which above calculations are based (operator P can no longer be approximated as above; the scaling arguments of Sect. 5 no longer hold). However, because here (and unlike references just cited) the matrix elements are strictly non-negative, m i j , h i j ≥ 0, the only way in which to obtain this scaling m, h, σ m and σ h is with a distribution that is concentrated at values of m extremely close to zero but which nevertheless has a very long tail. This scaling therefore converts most of the interactions in the community matrix into vanishing small ones (with a few very large ones to achieve the required mean) which seems rather at odds with the set-up of the model.
Phase transitions of the type described here have been described in related ecological models. Hogg et al. (1989) discussed an S × S community matrix with random elements distributed homogeneously about a finite mean, and showed that this can give rise to a macroscopic eigenvalue of size √ S (in the scaling convention pursued here). This model was also considered by  and  who described the emergence of a macroscopic eigenvalue as a phase transition from stability to instability as the mean value of the matrix elements is increased. The connection with the antagonist-mutualist system discussed here is that an increased mean interaction strength could arise from an increase in the number of mutualistic interactions. This scenario was discussed by Suweis et al. (2014) in another homogeneous model, introduced by Mougi and Kondoh (2012), that consists only of antagonists and mutualists. In these homogeneous models, the macroscopic eigenvalue can be understood in terms of the row sum of the community matrix K (Allesina and Tang 2012). In contrast, in the structured model discussed here, it is the row sum of the folded matrix W that is important. Sauve et al. (2014) considered the same tripartite topology as discussed here. Their conclusion was that stability was enhanced by the connectance and species diversity of the mutualistic subnetwork but decreased by the connectance and diversity of the antagonistic one. This is opposite to the behaviour described here (where at best increased mutualistic interactions can leave the stability unaltered in the strict asymptotic limit). There are several reasons for this discrepancy. Sauve et al. (2014) consider a relatively small number of species (compared with the limiting behaviour here) described by a more detailed dynamical model with a particular subset of parameters chosen to guarantee that the mutualist subsystem was independently stable. Perhaps even more importantly is that the stability to which they refer is that of a final equilibrium state, obtained through the time-evolution of the model. This will in general have fewer species than the starting community and possess additional structure. It is also clear that non-random structure can affect the stability of such networks, as Sauve et al. (2016) show in their comparison of random and empirical networks.
It seems therefore that an important future direction is to look at the feasibility of dynamical systems with interaction topology similar to that in Fig. 1. This would also allow better connection with the persistence studies of Sauve et al. (2014) and functional extinctions of Sellman et al. (2016). The dynamical cavity method has proved very useful for studying feasibility in homogeneous models (Bunin 2017, Galla (2018), and it remains on open possibility that this approach can be extended to structured models as studied here. Studying related few-species models, Ringel et al. (1996) commented that the destablizing effect of mutualisms "…is more than canceled by an increased chance of feasibility", and it will be interesting to see whether this also applies for large networks of interacting species and, in particular, in the presence of a phase transition.
Finally we note that our analysis can be generalised to a more complex merged networks, reminiscent of that in Pocock et al. (2012), which consist of a central guild of N P species (plants in Fig. 1) that interacts with a set of G other guilds, all of which are non-interacting with one another. Similar to Eq. (23), the interacting part of community matrix will be of the block form where A i and B i are random matrices of dimension N P × N i with N i the number of species in guild i. In contrast to Eq. (23), the interaction signs [e.g. (+,−) for antagonisms] here are given by the matrix elements, rather than being explicit in the block structure. Similar to above, finding the spectrum of this K reduces to the problem of finding the spectrum of W = in the asymptotic limit, whose position relative to the bulk around zero determines the stability. From the interaction signs, a guild of mutualists gives a positive contribution to λ, whilst and antagonist give a negative one. Competitive interactions in which elements of both A i and B i are negative would then also give a positive contribution to λ, thus moving the system in the direction of instability.
Acknowledgements This work was supported by Newcastle University's Institute for Sustainability. We acknowledge helpful discussions with Darren M. Evans.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A Finite-N P expression for the largest eigenvalue
In this appendix we provide details on how the approximate expression for λ 1 in Sect. 5 is obtained. Our starting point is the more general model of Eq. (23) and the corresponding W matrix of Eq. (24). The ensemble average of W evaluates as with φ and θ as in Eq. (25) and Eq. (27) and where 1 and J are the N P × N P unit matrix and matrix of ones, respectively. The eigenvalues of W are φ + N P θ and φ, this latter being (N P −1)-fold degenerate. The eigenvector u (1) belonging to eigenvalue φ + N P θ has elements u (1) In the degenerate subspace, we choose the eigenvector set u (α) : 2 ≤ α ≤ N P with elements These states are normalised, orthogonal to u (1) but not orthogonal to each other. Indeed we have These eigenvectors are easier to work with than orthogonalised ones, and since we work in the N P limit, the effect of the finite overlap u (α) · u (β) is negligible.
Arranging u (α) as columns of matrix U , the full W matrix in this (ensembleaveraged) basis is V ≡ U T WU , with matrix elements in which indices α, β ≥ 2. Considering the scaling of the ensemble averages of these quantities, we find V 11 ∼ N P , V α1 = V 1β = 0, V αα ∼ N 0 P , and V αβ ∼ N −1 P for α = β. Similarly, the variances scale as Var (V 11 ) ∼ N 0 P , Var (V α1 ) ∼ Var V 1β ∼ N 0 P , Var V αβ ∼ N −1 P . Thus, with N P large, we drop V αβ for α = β as being negligible, and replace the diagonal elements with their ensemble averages, as their fluctuations are negligible in comparison. Due to the vanishing of their mean, the elements V α1 and V 1β are left as is. Taken together, this gives the approximation This can be diagonalised exactly. We obtain N P − 2 eigenvalues of value φ and two nontrivial ones. These latter have the form of λ ± from Eq. (26) with Evaluating the leading term in the ensemble average of Ω 2 we find A similar calculation for the variance show that it scales like Var Ω 2 ∼ N −1 P such that the fluctuations about the mean vanish in the limit. We thus approximate Ω 2 with Ω 2 , thus giving the result stated in Eq. (28).