Second largest Eigenpair Statistics for Sparse Graphs

We develop a formalism to compute the statistics of the second largest eigenpair of weighted sparse graphs with $N\gg 1$ nodes, finite mean connectivity and bounded maximal degree, in cases where the top eigenpair statistics is known. The problem can be cast in terms of optimisation of a quadratic form on the sphere with a fictitious temperature, after a suitable deflation of the original matrix model. We use the cavity and replica methods to find the solution in terms of self-consistent equations for auxiliary probability density functions, which can be solved by an improved population dynamics algorithm enforcing eigenvector orthogonality on-the-fly. The analytical results are in perfect agreement with numerical diagonalisation of large (weighted) adjacency matrices, focussing on the cases of random regular and Erd\H{o}s-R\'enyi graphs. We further analyse the case of sparse Markov transition matrices for unbiased random walks, whose second largest eigenpair describes the non-equilibrium mode with the largest relaxation time. We also show that the population dynamics algorithm with population size $N_P$ does not actually capture the thermodynamic limit $N\to\infty$ as commonly assumed: the accuracy of the population dynamics algorithm has a strongly non-monotonic behaviour as a function of $N_P$, thus implying that an optimal size $N_P^\star=N_P^\star(N)$ must be chosen to best reproduce the results from numerical diagonalisation of graphs of finite size $N$.


Introduction
The second largest eigenvalue and the associated second eigenvector of a N × N matrix J is of great significance in many areas of science, with plenty of applications. In coding theory, the Hamming distance of a binary linear code can be expressed as a function of the second largest eigenvalue of the coset graph associated to the code [1]. In biology, it has been shown in [2] that the second largest eigenvalue of cancer metabolic networks describes the speed of cancer processes. In the context of clustering methods based on the adjacency matrix of a graph, the second eigenvector encodes inter-cluster connectivity, complementing the information about intra-cluster connectivity included in the top eigenvector [3,4]. Moreover, in Principal Component Analysis, the second eigenvector of the covariance matrix of standardised data represents the direction that accounts for the second largest source of variability within the dataset [5,6].
The second largest eigenvalue plays a pivotal role in the study of complex systems and graph theory, representing topological features of the graphs [7]. If the spectral gap, i.e. the distance between the largest and second largest eigenvalue, is large, then the graph has good connectivity and expansion properties [8]. Therefore, many results have been derived about bounds for the second largest eigenvalue (see e.g. [9,10]). In particular, bipartite regular graphs with very wide spectral gap are called expanders (magnifiers if not bipartite) and have been widely studied since the seminal work of Alon [11]. To shed light on the expansion properties of regular graphs, specific bounds have been derived for their second largest eigenvalue (see e.g. [1] and [12]).
The knowledge of the spectral gap is essential for random walks on undirected graphs, which are substantially equivalent to finite time-reversible Markov chains, as pointed out by Lovasz in his survey [13]. Indeed, up to log-factors, the inverse spectral gap of the transition matrix represents the mixing rate of the Markov chain, i.e. how fast the state probability vector of a Markov chain approaches the limiting stationary distribution [14], given by the top right eigenvector of the transition matrix. The inverse of absolute value of second largest eigenvalue of the transition matrix denotes the largest relaxation time or mixing time, and the corresponding eigenvector describes the nonequilibrium mode with the slowest decay rate. The second largest eigenpair of Markov transition matrices also plays an important role in all processes that are described by means of random walks on graphs, such as out of equilibrium dynamics of glassy systems (see e.g. [15,16]) and search algorithms such as Google PageRank [17].
In our analysis, we will be dealing with sparse symmetric random matrices, i.e. weighted adjacency matrices of undirected graphs. We focus on the case of high sparsity, i.e. when the probability of two nodes being connected is p = c/N , with c being the constant mean degree of nodes. In this sparse limit, numerical studies have shown that most of the eigenvectors of a random regular graph, as well as almost-eigenvectors ‡ [18], follow a Gaussian distribution [19], whereas Erdős-Rényi eigenvectors are localised especially for low values of c. The statistics of the first eigenvector components for very sparse symmetric random matrices has been first considered in the seminal work by Kabashima and collaborators [20] and subsequently in a more systematic way in our previous work [21].
Following the framework developed in [21], we look at the second largest eigenpair problem as the top eigenpair problem for a deflated version of the original sparse matrix (see discussion in Section 2). We will be implementing a Statistical Mechanics formulation of the top eigenpair problem of the deflated matrix, using both the cavity (Section 3) and replica (Section 4) methods in a unified way.
Both the replica and cavity methods from the physics of disordered systems have been employed in the realm of random matrix theory for a long time. The replica method, traditionally used in the physics of spin glasses [22], was first introduced in the context of random matrices by Edwards and Jones [23] to compute the average spectral density of random matrices defined in terms of the joint probability density function of their entries. Later on, the same approach proved useful to derive the spectral density of Erdős-Rényi adjacency matrices as the solution of a (nearly intractable) integral equation in the seminal paper of Bray and Rodgers [24]. Later, approximation schemes such as the single defect approximation (SDA) and the effective medium approximation (EMA) [25,26] were developed. An exact alternative approach was introduced in [27]: starting from Bray-Rodgers replica-symmetric setup [24], the functional order parameters of the theory are expressed as continuous superpositions of Gaussians with fluctuating variances, as suggested by earlier solutions of models for finitely coordinated harmonically coupled systems [28]. This formulation gives rise to non-linear integral equations for the probability densities of such variances, which can be efficiently solved by a population dynamics algorithm. We will follow a similar approach in Section 4.
The cavity method [29] has been employed in the study of disordered systems and sparse random matrices as a more direct alternative to replicas. It is exact for highly sparse tree structures [30]. As shown in [31], one of the advantages of the cavity method is that it provides the spectral density for very large single instances of sparse random graphs. Both methods, known to lead to the same results for the spectral density [32], recover the Kesten-McKay law for the spectra of random regular graphs [33,34], the Marčenko-Pastur law and the Wigner's semicircle law respectively for sparse covariance matrices, and for Erdős-Rényi adjacency matrices in the large c limit [27,31]. Likewise, the spectral density of sparse Markov matrices [35,36], graphs with modular [37] and small-world [38] structure, and with topological constraints [39] have been obtained. Also, the spectral density in the complex plane of sparse non-Hermitian matrices has been considered in [40][41][42][43].
As in [21], we will provide a cavity single-instance derivation for our problem. Generalising the single-instance results in the thermodynamic limit, we will show that even for the second eigenpair problem the cavity method leads to the same stochastic recursions obtained from the replica treatment. The crucial difference between the present work and [21] is the presence of the orthogonality condition between the top and second eigenvectors in the set of final recursion equations. The population dynamics algorithm employed to solve these recursions, complemented by a wise implementation of the orthogonality constraint, allows us to characterise the distributions of the cavity fields in the thermodynamic limit, and to disentangle the individual contributions of different degrees to the second eigenvector's entries. The plan of the paper is as follows. In Section 2, we will formulate the problem in terms of deflation and provide the main starting point. In Section 3, we will describe the cavity approach to the problem, first for the single instance case, and then in the thermodynamic limit, highlighting the role of the orthogonality constraint (in 3.3). In Section 4, we formulate a general replica approach to the second largest eigenpair problem for a general degree distribution, first focussing on the second largest eigenvalue problem (in 4.1), then on the density of corresponding eigenvector's components (in 4.2), and again emphasising the role of the orthogonality condition (in 4.3). In Section 5 we focus on the case of the random regular graph: we analytically show how the solution for the top eigenpair of the deflated adjacency matrix gets modified as the deflation parameter is changed. In Section 6, we specialise our results to the case of Markov transition matrices representing random walks on graphs. In Section 7, we provide the details of the population dynamics algorithm, focussing on how the extra orthogonality constraint is implemented. We also provide convincing evidence that -at odds with what is commonly believed -the algorithm with finite population size N P does not actually capture the thermodynamic limit N → ∞, in that there is a non-trivial relation between the size N of the adjacency matrix being diagonalised, and the size N P of the population one should ideally use to numerically compute its spectral properties. More precisely, the accuracy -measured with different metrics -with which the population dynamics algorithm reproduces numerical diagonalisation of matrices (graphs) of size N has a strongly non-monotonic behaviour as a function of N P , thus implying that an optimal size N P = N P (N ) must be chosen to best reproduce the diagonalisation results. Finally, in Section 8 we offer a summary of results.

Formulation of the problem
Given a real symmetric matrix J = (J ij ) and its top eigenpair (λ 1 , u), we define a deflated matrixJ(x) = (J ij (x)) bỹ In the present paper we will be mostly concerned with sparse matrices, where J ij = c ij K ij are the i.i.d. entries of the sparse symmetric matrix J. The top eigenvector u of J is normalised to N . The dense matrix uu T /N represents the projector onto the top eigenspace of the original matrix J. The entries of the original matrix J are defined in terms of the connectivity matrix c ij ∈ {0, 1}, i.e. the adjacency matrix of the underlying graph, and the random variables K ij encoding bond weights. Within our formalism, we are able to handle any kind of highly sparse degree connectivity -where the mean node degree k = c is a finite constant that does not scale with N (entailing c/N → 0 as N → ∞).
The bond weights K ij will be i.i.d. random variables drawn from a parent pdf p(K) with bounded support. This setting is sufficient to ensure that the largest eigenvalue λ 1 of J will remain of ∼ O(1) for N → ∞.
For any value of x, the matrices J andJ(x) share the same set of eigenvectors (see Section 3.3.2 in [44]). The range of x is [0, λ 1 ], where the boundaries of this range correspond respectively to no deflation (x = 0 ⇒ J =J) and full deflation (x = λ 1 ). When the value of x is larger than the spectral gap -which e.g. happens when x = λ 1 , i.e. the largest eigenvalue of the original matrix J § -then the second largest eigenvalue of J and the corresponding eigenvector become the top eigenpair of the matrixJ. Conversely, the top eigenvector of J, u, is still an eigenvector ofJ, but corresponding to a zero eigenvalue. All other eigenpairs are unchanged.
By setting up a formalism based on the statistical mechanics of disordered systems, we aim to find the average λ 2 J (or typical) value of the second largest eigenvalue λ 2 of J, and the density ρ J,2 (w) = 1 2 ). The second eigenpair statistics of the matrix J is obtained by finding the top eigenpair of the deflated matrixJ(x) when x = λ 1 . Thus, in order to obtain the desired quantities, we analyze the average largest eigenvalue λ 1 J and the density ρJ 1 ) of the deflated matrixJ, where the average · J is taken over the distribution of the matrixJ.
We provide: • the second largest eigenpair statistics λ 2 J and ρ J,2 (w) of the matrix J, i.e. the solution corresponding to the maximum deflation forJ, in the case of a generic connectivity p(k) with bounded maximum degree, found via the cavity method (Section 3). We also offer an equivalent replica derivation for the same problem (Section 4). In this general case, the solution is available via population dynamics simulations (Section 7); • an explicit analytical solution for λ 1 J and ρJ (w) in the specific case of the adjacency matrix of a random regular graph (RRG), showing that the solution requires the deflation parameter x to exceed the spectral gap (Section 5); • the second largest eigenpair statistics of the unbiased random walk Markov transition matrix. In this case, the deflation parameter x is set precisely to 1, i.e. equal to the largest eigenvalue of the Markov transition matrix. Also in this case, an analytical description is provided for the RRG connectivity case (Section 6).
The equations (37), (38), (39) and (40) found within the cavity framework (see Section 3.3) represent the solution of the second largest eigenpair problem in the § In the thermodynamic limit, the value such that full deflation is achieved is actually the average largest eigenvalue λ 1 J of the matrix J. thermodynamic limit, and constitute the main result of this paper. We notice that they are completely equivalent to the equations (117), (118), (119) and (121) found within the replica framework (see Section 4.3).
We will follow the same protocol used in [21]. Focussing on the matrixJ, the problem can be formulated as the optimisation of a quadratic functionĤ(v), according to which v 1 is the vector normalised to N that realises the condition as dictated by the Courant-Fischer definition of eigenvectors. The round brackets (·, ·) indicate the dot product between vectors in R N . It is easy to show thatĤ (v) is bounded and attains its minimum when computed on the top eigenvector. For a fixed matrixJ, the minimum in (3) can be computed by introducing a fictitious canonical ensemble of N -dimensional vectors v at inverse temperature β, whose Gibbs-Boltzmann distribution reads where the delta function enforces normalisation. Clearly, in the low temperature limit β → ∞, only one "state" remains populated, which corresponds to v = v 1 , the top eigenvector of the matrixJ.

Full deflation: cavity method
In this section, we present the cavity derivation of the single instance equation for the second largest eigenpair problem. The formalism shown here differs from that presented in [21]: here we analyse the partition function of the Boltzmann distribution (5), rather than a soft-constrained version of it. This allows us to include hard constraints within the cavity framework. The equations expressing the solution can be easily generalised to the thermodynamic limit case, reproducing the same equations that will be found by the replica formalism in Section 4, which constitute the main results of this work. We consider a random N × N symmetric matrixJ = J ij , with real entries. The matrix entries are defined asJ where J ij = c ij K ij are the entries of the sparse symmetric random matrix J, λ 1 is its largest eigenvalue and u its top eigenvector, normalised to N , which we assume to be known. The vector u will be also referred to as the probe eigenvector. The entries of the original matrix J are defined in terms of the c ij ∈ {0, 1}, which is the connectivity matrix, i.e. the adjacency matrix of the underlying graph, and the K ij , which encode bond weights. The dense matrix 1 N uu T represents the projector onto the first eigenspace of the original matrix J. The matrixJ is also said to be a fully-deflated matrix.

Cavity single-instance derivation for the largest eigenvalue ofJ
In the full deflation setting, given a single instance matrixJ, its largest eigenvalueλ 1 is the second largest eigenvalue λ 2 of the original matrix J. It can be defined as The partition function explicitly reads The square in the exponent can be written as with the identification The definition of the order parameter q is enforced via the integral identity By also employing a Fourier representation of the Dirac delta enforcing the normalisation constraint, the partition function then becomes where and C includes all the pre-factors.
At the saddle point, the following stationarity conditions hold, where the angular brackets indicate averaging w.r.t. the distribution By looking at the saddle point condition (14), in what follows we can identify iq = λ 1 q and define iλ ≡ λ, such that (17) becomes The components v i are found in the β → ∞ limit by the cavity method applied to the distribution (18) (see also Subsection 3.2). Here we will follow the protocol detailed in Section 3.1 of [21], reporting the key steps to make this paper self-contained.
By making a tree-like assumption on the structure of the highly sparse graph encoded in the original matrix J that we deflate, the marginal pdf w.r.t. a certain component i is given by where ∂i denotes the immediate neighbourhood of i. The factorisation over the neighbouring nodes of i is due to the fact that in a tree-like structure the nodes j ∈ ∂i are connected with each other only through i. The distribution P (i) j (v j |λ 1 q, λ) is called marginal cavity distribution: it is the distribution of the components v j defined on the neighbouring nodes of i, in the network in which i has been removed.
In the same way, (see for instance Eq. (11) in [21]), for any j ∈ ∂i the cavity marginal pdf satisfies the self-consistent equation where ∂j\i indicates the set of neighbours of the node j with the exclusion of i.
A Gaussian ansatz provides the solution to the self consistent equation, viz.
where the parameters ω (i) j and h (i) j are called cavity fields. By inserting the ansatz in (20) and performing the Gaussian integrals, the set of self-consistent equations represented by (20) translates into a set of recursions for the cavity fields, Likewise, by means of (21), the marginal distribution P i (v i |λ 1 q, λ) can be written as where In the β → ∞ limit, entailing that the components v i of the second largest eigenvector, representing the ground state of the system with Boltzmann distribution (5), are given by the ratios h i /ω i . A very quick proof of this will be given in Section (3.2). Therefore (15) and (16) become Considering (7) and (12), the largest eigenvalue ofJ (corresponding to the second largest eigenvalue of J) is then given bỹ

Cavity single-instance derivation of the top eigenvector ofJ
We rewrite the Boltzmann distribution (5) to perform explicitly the β → ∞ limit and show that the components of the ground state vector, i.e. the top eigenvector of the deflated matrixJ, are the ratios h i /ω i , as anticipated in (27). We recall that -since we are considering a full deflation -this eigenvector is the one corresponding to the second largest eigenvalue of the original matrix J, defined in (30). Indeed, Eq. (5) can be written as where we have employed (10), (11) and the Fourier representation of the Dirac delta. Using a saddle-point approximation of the integral appearing in (31) for N finite but large, we get where we have again used (14), (15) and (16) and the subsequent identifications iq = λ 1 q, iλ ≡ λ. The normalisation constant Z N defined in (8) can be in turn evaluated via a saddle-point approximation, as already done in (12), leading to where Z N (λ 1 q, λ) is defined in (13). Inserting (33) in (32), we find that the pre-factors and the exponential with square brackets cancel out. The distribution (5) then reduces to (18) . As done in Subsection 3.1, Eq. (18) can be manipulated as in Section 3.1, where ω i and h i are defined by (25) and (26). Computing the β → ∞ limit, we find It can be noticed that the distribution (33) is substantially equivalent to the grand-canonical distribution (Eq. (7) in [21]) which we adopt as the starting point of the cavity treatment in [21]. as expected.
Eq. (35) defines the components of the top eigenvector v of the fully-deflated matrix J in terms of (25) and (26). Because of the full deflation, v also represents the second largest eigenvector of the original matrix J. Moreover, in view of the orthogonality between u and v, it follows that q = 0, viz.
where the components u i and v i are naturally referring to the same node i with degree k i of the network represented by J. To summarise, in the single instance case the solution is given by the cavity recursions (22) and (23) along with the normalisation condition (29) and the orthogonality constraint (36). The value λ = λ 2 represents the second largest eigenvalue, and according to the same mechanism explained in Appendix A of [21], it is the only value that satisfies the normalisation condition (29).

Cavity method: thermodynamic limit and orthogonality condition
Following the reasoning of Section 3.2 in [21], in the limit N → ∞ we can consider the joint probability density of the cavity fields ω j taking values around respectively ω and h, where dπ (ω , h ) ≡ dω dh π (ω , h ), and the average · {K} k−1 is taken over k − 1 independent realisations of the random variable K. The distribution k c p(k) represents the probability that a randomly chosen link points to a node of degree k and c = k , and appears in (37) as cavity fields are related to links. Eq. (37) generalises in the thermodynamic limit the recursions (22) and (23) in the case of full deflation (x = λ 1 ) for which q = 0.
By using the law of large numbers, in the thermodynamic limit the normalisation condition reads whereas the orthogonality constraint becomes Here, ρ J (u|k) is the conditional distribution of the top eigenvector's component of J w.r.t. to the degree k. Similarly, the distribution of the top eigenvector's components of the fully deflated matrixJ, i.e. the second largest eigenvector of J, is obtained in terms of averages w.r.t. the distribution π(ω, h) as We notice that in the equations (38), (39) and (40), the degree distribution p(k) naturally crops up, as they encode properties related to nodes, rather than links. Moreover, Eq. (30) generalises to the thermodynamic limit case, expressing the fact that as found in [21] the parameter λ represents the average largest eigenvalue of the deflated matrix J, i.e. the average second largest eigenvalue of the matrix J, We anticipate that the latter result is equivalent to the average second largest eigenvalue explicitly found by the replica approach in Eq. (104).
Enforcing the orthogonality condition given by (39) is crucial to find the correct solution. The conditional pdf ρ J (u|k) appearing in (39) is given by omitting the k-sum in the expression for the density of the top eigenvector components ρ J (u) (see Eq. (111) in [21]), reported here where π 1 (a, b) indicates the distribution of cavity fields of type a and b for the top eigenpair problem ¶. The integration w.r.t. the conditional distribution ρ J (u|k) in (39) generalises to the thermodynamic limit the fact that both the components u i and v i = h i ω i in (36) refer to the same node i with degree k i . Indeed, by comparing (42) with (40) and (39), we notice that the components of u are still coupled to those of v in (39) through their structure, as they both refer to the same degree k (see for more details Section 7). The replica derivation provides an independent proof of this result in Section 4.3. Therefore, in order to enforce the constraint (39) correctly, we need to impose strict orthogonality on-the-fly, i.e. while the components of the top eigenvector u and the components of the second largest eigenvector v are being evaluated at the same time by averaging w.r.t. respectively π 1 and π, as prescribed by (42) and (40). The way strict orthogonality is imposed is via a correction to the components of v: the details of this procedure and the corresponding algorithm are given in Section 7. We remark that the condition q = 0 holds whenever x exceeds the spectral gap.
To summarise, the equations (37), (38), (39), (40) and (41) represent the solution of the second largest eigenpair problem in the thermodynamic limit and constitute the ¶ In the context of the top eigenpair problem, the cavity field of type a has the role of an inverse cavity variance (similarly to ω for the second largest eigenpair problem), whereas b represents a cavity bias (similarly to h here). See Section 3 in [21]. main result of this paper. This set of equations must be generally solved by a population dynamics algorithm, as detailed in Section 7. We anticipate that they will be completely equivalent to the equations (117), (118), (119), (120) and (121), respectively, found within the replica framework (See Section 4.3).

Full deflation: replica derivation
In this section, we evaluate the average (or typical) value of the largest eigenvalue and the density of top eigenvectors' components of the matrixJ within the replica framework. Our derivation applies to any graph with degree distribution p(k) having finite mean. For weighted adjacency matrices with a Poissonian distribution, we also ask that its support be bounded to ensure that their average largest eigenvalue is finite in the thermodynamic limit.

Typical largest eigenvalue
The u i represents the i-th component of u, the top eigenvector of the original matrix J (normalised to N ) which we assume to be known. We recall that u will be also referred to as the probe eigenvector. Within the framework of the configuration model [38], the joint distribution of the matrix entries J ij is where the distribution P ({c ij }|{k i }) of connectivities {c ij } compatible with a given degree sequence {k i } is given by and the pdf p (K ij ) of bond weights (over a compact support whose upper edge is denoted by ζ) can be kept unspecified until the very end. Our derivation will follow the procedure presented in Appendix B in [21].
Here we fix x = λ 1 J : in this setting, the second largest eigenvalue of J is given in terms of the largest eigenvalue ofJ. This can be computed as the formal limit in terms of the quenched free energy of the model defined in (5).
The partition function explicitly reads we can linearise the square in the exponent of (46) by means of a Hubbard-Stratonovich identity as follows, and therefore the partition function reads (48) The average overJ then reduces to computing the average over J. It is computed using the replica trick as follows where n is initially taken as an integer, and then analytically continued to real values in the vicinity of n = 0. The replicated partition function is Taking the average w.r.t. the joint distribution (44) of matrix entries yields [21,38] exp where the average · K is taken w.r.t. the pdf of the bond weights p(K). A Fourier representation of the Kronecker deltas expressing the degree constraints in (44) has been employed. Employing a Fourier representation of the Dirac delta enforcing the normalisation constraint, the replicated partition function thus becomes where we omit irrelevant proportionality constants.
In order to decouple sites, we introduce the functional order parameter where the symbol v denotes a n-dimensional vector in replica space. We then consider its integrated version [21,38] and enforce the latter definition using the integral identity In terms of the integrated order parameter (54) and its conjugate, the replicated partition function can be written as The multiple integral in the last two lines is the product of N n-dimensional integrals, each related to both k i and u i , i.e. the degree and the eigenvector component of the node i. It can be expressed by means of the law of large numbers in the following way: where Log denotes the principal branch of the complex logarithm, and Each of the φ i integrals can be performed by rewriting the last exponential factor as a power series, viz.
with i = 1, . . . , N . Therefore, by invoking the Law of Large Numbers, the single site integral I (57) can be expressed as where we have used (61) Here, p(k) is the actual degree distribution of the graph and ρ J (u|k) represents the distribution of the top eigenvector's components of the original matrix J conditioned on the degree k. As shown in [21], the variables u i are strongly correlated to the k i so their dependence on the k i must be taken into account.
Therefore, the replicated partition function takes a form amenable to a saddle point evaluation in the large N limit (assuming we can safely exchange the limits n → 0 and N → ∞) where and where we consider k min = 0 henceforth. The stationarity of the action S n w.r.t. variations of ψ andψ requires that the order parameter at the saddle point ψ and its conjugateψ satisfy the following coupled equations which have to be solved together with the stationarity conditions w.r.t. each component λā of λ and zā of z (forā = 1, . . . , n), Apart from the extra averages w.r.t. p(k) and ρ J (u|k), the equations (69) and (70) share some similarities with the saddle-point equations leading to the spectral density of sparse random graphs [24,27] and to those leading to the top eigenpair statistics of sparse symmetric matrices [21]: similarly to the latter case, the harmonic "Hamiltonian" of this problem is real-valued and includes the inverse temperature β. Following [21,27], we will now search for replica-symmetric solutions written in the form of uncountably infinite superpositions of Gaussians with a non-zero mean. As in the case for the top eigenvector, our ansatz will be preserving permutational symmetry between replicas, but not the rotational invariance in replica space, since this symmetry would not produce a physically meaningful result for this problem. where We remark that our replica symmetry assumption has proved to be generally exact in the random matrix context and specifically for the spectral problem of sparse random matrices [23,24,27,45]. Moreover, the representation of the order parameter as a superposition of Gaussian pdfs leads to the correct solution for harmonically coupled systems [28], such as the one described in the present work.
The calculation follows the same path traced in Appendix B of [21]. In (75) and (76), π andπ are auxiliary normalised joint pdfs of the parameters appearing in the Gaussian distributions. The ψ 0 andψ 0 are determined such that the distributions π(ω, h) andπ(ω,ĥ) are normalised.
Expressing the order parameter in this form allows us to perform explicitly the vintegrals in the action S n , eventually leading to simpler coupled stationarity equations for π andπ. The convergence of the v-integrals requires ω >ω and ω > ζ (where ζ is the upper edge of the support of the pdf p(K) of bond weights).
Rewriting the action in terms of π andπ, after performing the v-integrations, and extracting the leading n → 0 contribution the full action now reads with where we have introduced the shorthands We then introduce in (85) the scalar order parameter via the integral representation We intentionally choose to label the scalar order parameter (86) as ψ 0 , since this highlights that the O(1) terms in the action (78) match the terms arising from the saddle-point evaluation of the normalisation constant M, exhibiting the same scaling in N . Indeed, by using the same argument as in (60), we find The stationarity conditions for S M are and entailing that The two conditions above exhibit a gauge invariance [21,38]. Once the same gauge has been fixed for the saddle-point solution of M and the O(1) terms of the action (78) in the numerator, they cancel out so that the action (78) is O(n) as expected.
The stationarity condition w.r.t. λ entails where the average · P is taken with respect to the Gaussian measurē More explicitly, (93) reads We note that the β-dependent term vanishes as β → ∞.
The stationarity condition w.r.t. z entails where the average · P is taken with respect to the Gaussian measure (94). More explicitly, The stationarity condition with respect to variations of π, δS δπ = 0, iŝ Similarly, the stationarity condition with respect to variations ofπ, δS δπ = 0, produces the condition where the brackets · {K} k−1 denote averaging with respect to a collection of k − 1 i.i.d. random variables K, each drawn from the bond weight pdf p(K). We recall that p(k) appearing in (100) is already the actual degree distribution of the graph with finite mean c and bounded maximal degree. Following [21], we relabel the constant terms λ ≡ iλ and q ≡ −iz since they both turn out to be real-valued. We eventually find The parameter λ must be tuned as to enforce the supplementary condition (95) as β → ∞, which reads whereas (97) gives the following condition for q The structure of the action (78) is the same as that found in [21] (see for instance Section 4.1.1 there), except for the term S 4 (z) ≡ S 4 (q). Therefore, building on the same reasoning, the average largest eigenvalue ofJ, i.e. the average second largest eigenvalue of J is given by where λ and q are defined by (102) and (103). Eq. (101), along with the conditions (102) and (103), are typically solved by population dynamics, as shown in Section 7. They represent the generalisation in the large N limit of the single-instance recursions (22) and (23) along with the conditions (28) and (29).

Density of top eigenvector's components
In this section, we provide the derivation for the density of components of the top eigenvector of the matrixJ, in terms of π, λ and q. As in the previous subsection, we consider the deflation parameter x = λ 1 J , and therefore the top eigenvector of the deflated matrixJ corresponds to the second eigenvector of the original matrix J. We will be following the same approach of Section 4.2 in [21]. We will report here the main steps to keep this paper self-contained. In this statistical mechanics framework, the quantityρ is defined such that in the limit β → ∞ it gives the density of the top eigenvector components for a given N × N deflated symmetric random matrixJ. The simple angle brackets ... stands for thermal averaging with respect to the Gibbs-Boltzmann distribution (5) of the system Defining an auxiliary partition function as where δ is a smooth regulariser of the delta function, the quantity (105) can be formally expressed as Averaging now over the matrix ensemble and sending β → ∞ at the very end, the density of the top eigenvector's components is eventually given by the formula equivalent to Eq. (95) in [21].
To compute the average of the logarithm of the auxiliary partition function Z (β) (t,J; w), we employ the replica trick The replicated partition function takes the form where ψ andψ are functional order parameters + . For large N , we employ a saddle point approximation where the starred objects satisfy self-consistency equations where t can be safely set to zero, since the partial derivative ∂ ∂t in (110) only acts on terms containing an explicit dependence on t. The constant M is the same normalisation defined in (85).
The stationarity conditions defining ψ ,ψ , λ and z at the saddle point for t = 0 are identical to those found in Section 4.1. The explicit n-dependence of the action S (β) n ψ ,ψ , λ , z ; t, ; w is again extracted by representing the order parameters ψ andψ as infinite superpositions of Gaussians. The explicit t-dependence appears in the so-called "single-site" term of the action, i.e.
+ We use the same symbols ψ andψ as in Section 4.1.
By making the identifications iλ ≡ λ and q ≡ −iz as before, taking the t-derivative at t = 0 and considering the limits → 0 and β → ∞ as prescribed by (110), we eventually find

The orthogonality condition
Here we show that in case of full deflation, x = λ 1 J , Eq. (103) encodes the orthogonality-on-average condition between w and the probe eigenvector u. Indeed, the orthogonality condition reads where P J (u, w) indicates the joint probability density of the first and second largest eigenvector's components of J and the conditional pdf ρ J,2 (w|u, k) is obtained by (115) without considering the u-integration and the k-sum. The conditional pdf ρ J (u|k) is given by omitting the k-sum in the expression for the density of the top eigenvector components ρ J (u) (42) (see Section 3.3). The comparison between (116) and (103) implies that q = 0. Taking into account the average orthogonality condition q = 0, the equations (101), (102), (103), (104) and (115) simplify to We remark here that all the observations made in Section 3.3 about the orthogonality condition expressed in (119) hold here. We observe that the components of u are still coupled to those of w in (119) through their structure, as they both refer to the same degree k. Therefore, as anticipated in Section 3.3, we need to impose strict orthogonality on-the-fly by correcting the components of w: all the details of this procedure and the corresponding algorithm can be found in Section 7.
Eq. (117), (120) and (121) complemented with the conditions (118) and (119) provide the solution of the second largest eigenpair problem in the large N limit. They are completely equivalent to the thermodynamic limit equations (37), (40) and (41), which -along with the conditions (38) and (39) found within the cavity methodembody the main results of this paper. Eq. (117) is solved via a Population dynamics algorithm, taking into account the conditions (118) and (119), as detailed in Section 7. Figure 1 shows the numerical results in the case of an ER adjacency matrix with c = 4 and k max = 12. We find λ 2 J = 4.463, within a 2% error w.r.t. the value λ 2,∞ = 4.565 obtained by extrapolation from the direct diagonalisation data. The bottom right panel of Figure 1 refers instead to the case of ER weighted adjacency matrix with c = 4 and k max = 12. We consider the case of uniform distribution of bond weights, p(K) = 1/2 for K ∈ [1,3]. In this case, we find λ 2 J = 9.5016, within a 2.5% error w.r.t. the reference value λ 2,∞ = 9.7452 obtained by extrapolation from the direct diagonalisation data. In the plot, we compare the pdf of second largest eigenvector's components obtained via population dynamics with results from the direct diagonalisation of 2000 matrices of size N = 5000. Figure 2 compares the second largest eigenvector's components pdf obtained via population dynamics in the case of ER adjacency matrix with c = 10 and k max = 22 with the results obtained via direct diagonalisation of matrices of increasing size. In this case, we find λ 2 J = 6.656, within a 0.4% error w.r.t. the value λ 2,∞ = 6.658 obtained by extrapolation from the direct diagonalisation data. We observe that there are finite size effects in the distribution of eigenvector components that are significantly stronger than those observed in the eigenvalue problem.

Random regular graphs
For non-weighted adjacency matrices of RRGs, the degree distribution is simply p(s) = δ s,c , and the bond weights distribution is trivially p(K) = δ(K − 1), resulting in a constant probe top eigenvector u, i.e. ρ J (u) = ρ J (u|c) = δ(u−1). The largest eigenvalue λ 1 is non-random and pinned to the value λ 1 = c. The spectral density is given by the and N = 5000 (yellow). As N increases, we notice that the direct diagonalisation curves approach the pdf generated by population dynamics with a fairly large population size, N P = 10 5 . Figure 3),

Kesten-McKay distribution (See
In this section we look at the behaviour of the solution for a generic value of the deflation parameter x in the range [0, c]. Therefore, the value of q is in principle non-zero. We remark that q = 0 holds surely in the case of full deflation, as in Section 3 and 4. For a general value of the deflation parameter x, the equation (101) for π(ω, h), along with the conditions (102) and (103) become respectively and the density of top eigenvector's components of the deflated matrixJ (115) is given for general x by We will show that the solution of the self-consistency equation (123) along with (124), (125) and (126) crucially depends on the value of the deflation parameter x. We recall here that the range of x is [0, c], where the boundaries of this range correspond respectively to no deflation (x = 0) and full deflation (x = c). We anticipate that in the outer regime 0 ≤ x < c−2 √ c − 1 (see Figure 3), the probe eigenvector u = {1, 1, . . . , 1}, i.e. the top eigenvector of the original matrix J, is also the top eigenvector of the deflated matrixJ, with corresponding largest eigenvalue c−x lying outside the bulk of the Kesten-McKay spectrum [33,34]. Conversely, in the bulk regime i.e. when c − 2 √ c − 1 < x ≤ c, the top eigenvector's components density is a standard normal distribution, with corresponding largest eigenvalue 2 √ c − 1. The probe all-one eigenvector u is still an eigenvector ofJ but refers to an eigenvalue c − x < 2 √ c − 1. In other words, we show that the second largest eigenpair of the RRG adjacency matrix is given by λ 2 J = 2 √ c − 1 and ρ J,2 (w) = N (0, 1). Figure 3 explains graphically the outer and bulk regimes.
The abrupt change of the solution (from constant u to normally distributed when x hits the value c − 2 √ c − 1) reflects the fact that the usual peaked ansatz for the RRG case (see [21]) is not valid in the bulk regime c − 2 √ c − 1 < x ≤ c. Therefore, in order to solve the self-consistent equation (123), we choose a "mixed" ansatz of the form for realω andh.
We further show that in the range 0 ≤ x < c − 2 √ c − 1, the solution is provided by a peaked ansatz, i.e. σ 2 = 0 -just like in the case of the largest eigenpair of the original matrix J -whereas in the range c − 2 √ c − 1 ≤ x < c, the variance σ 2 must be finite. Indeed, by inserting (127) into (123) and performing the r.h.s. integrals, we find (128) Comparing (128) with the ansatz (127), we find that the following relations must be satisfiedω From the last condition (131), we can infer that if σ 2 > 0, thenω = √ c − 1, i.e. a finite variance of the distribution of components pinsω to a specific value. Only if σ 2 = 0 , thenω can assume values other than √ c − 1, according to Eq. (129). Inserting the ansatz (127) in the normalisation condition (124) and in the condition (125), we find two extra conditions to fix respectively σ 2 and q, By combining (129), (130) and (133), we find an expression for q in terms ofω and h, which in turn can be inserted into Eq. (130) to givē If we equate the expression in the round brackets in (135) to zero and compare it to (129) after a slight rewriting, we find the value of λ explicitly, From (136), we also find the explicit dependence ofω on x. Indeed, by solving forω we getω which yields a x-dependent real solution only for 0 ≤ x < c − 2 √ c − 1. Only in this regime,ω =ω(x) can assume values other than √ c − 1, entailing from (131) a peaked solution for π. * Conversely, for any x > c − 2 √ c − 1, Eq. (137) would produce a x-dependent complex solutionω(x), which is not acceptable for this problem (recall that ω and h must be real), thus implying

RRG-deflated top eigenvalue: outer regime
From (138), it follows that σ 2 = 0 in the outer regime. From (132) and (133), we thus find When solving (139), we must discard the other possible solution q = 0, since it would not satisfy the normalisation constraint (124). Equipped with this information and also taking into account (127), (136) and the identityh =ω − 1, which follows from (134), the O(n) terms of the action S n in (78)keeping only the leading β → ∞ terms -are expressed as follows x , (143) * We remark that in this regime a finite variance solution for π that pinsω to √ c − 1 is still possible, but yields a higher ground state free energy F J than the peaked solution. Indeed, F J = − N 2 λ 1 J . See Sections 5.1 and 5.3.
Summing up all terms and recalling from (136) that λ = c − x, the action at the saddle point reads which implies from (49) for the average of the largest eigenvalue ofJ the formula Therefore, the deflation with a parameter x in the regime 0 ≤ x < c − 2 √ c − 1 has the effect of decreasing the top eigenvalue c of the original RRG adjacency matrix J by a quantity x, as long as it lies outside the spectral bulk of the Kesten-McKay distribution. In the next subsection, we will show that the corresponding eigenvector is still the top eigenvector of J.

RRG-deflated density of top eigenvector components: outer regime
As found at the beginning of this Section, within the range 0 ≤ x < c − 2 √ c − 1 , the ansatz for π is delta-peaked, since σ 2 = 0. We show that a peaked ansatz of this sort corresponds to the top eigenvector of the matrixJ being all-ones: this means that for 0 ≤ x < c − 2 √ c − 1 the top eigenvector ofJ is exactly the probe eigenvector u. Indeed, by inserting the ansatz (127) in (126) and taking into account (136) and (139), we find but, from (139), implying where the choice of the "+" sign solution is not restrictive.
In conclusion, as long as the largest eigenvalue c − x of the deflated matrixJ lies outside the spectral bulk (i.e. for 0 ≤ x < c − 2 √ c − 1), the corresponding top eigenvector w is equal to the probe eigenvector u = (1, ..., 1) T , i.e. the top eigenvector of J.

RRG top eigenvalue: bulk regime
In this range, we have shown in (138) that the variance σ 2 is positive, giving rise to a mixed "delta-Gaussian" ansatz for π. The parameter σ 2 being positive implies that ω must be pinned to the value √ c − 1. From (129), it follows that λ = 2 √ c − 1. The values of q andh are determined by the normalisation (124) and orthogonality (125) conditions. Indeed, the change in the ansatz corresponds to a change in the structure of the largest eigenvector w ofJ. As shown in Section 4.3, the orthogonality condition reads where ρJ (u|c) = δ(u − 1) is the conditional distribution of the probe eigenvector entries and (126) has been used. Comparing (150) with (125) we infer that q = 0. Moreover, inserting q = 0 in (132) and (133), we can respectively infer that Equipped with this information and also by taking into account (127), the O(n) terms of the action S n in (78) -keeping only the leading β → ∞ terms -are expressed as Summing up all terms and exploiting the identities (129) and (130), the action at the saddle point reads which implies from (49) that the average of the largest eigenvalue ofJ is corresponding to the upper edge of the Kesten-McKay distribution. As expected, the eigenvalue does not depend on the normalisation of the corresponding eigenvector, encoded in σ 2 . Since this result holds for any x in c − 2 √ c − 1 < x ≤ c, including the case of full deflation when x = λ 1 J = c and the first eigenmode u of the original matrix J is associated to a zero eigenvalue, we conclude that the average second largest eigenvalue of the matrix J is

RRG density of top eigenvector components: bulk regime
In this range of values for x, we show that the "delta-Gaussian" ansatz for π(ω, h) leads to a Gaussian-distributed top eigenvector of the matrixJ. Since this result is valid also in case of full deflation, i.e. x = c, we can conclude that the eigenvector corresponding to the second largest eigenvalue of a random regular graph adjacency matrix J is normally distributed . We then identify in x = c − 2 √ c − 1 ⇐⇒ λ = 2 √ c − 1 a transition point for the structure of the distribution of the top eigenvector's components ofJ(x), at which the parameter q changes discontinuously from q = ±1 to 0.
We now evaluate the density of the top eigenvector components in the range Inserting the ansatz (127) in (126) and taking into account (138), (150), (151) and (152), we find We remark that this analytical result is in excellent agreement with the statistics of the second largest eigenvector components of the RRG adjacency matrices found by population dynamics, as shown in Figure 4. Moreover, it is compatible with previous known results about eigenvectors of random regular graphs [18,19].

Sparse random Markov transition matrices
In this section, we apply the deflation formalism to an ensemble of transition matrices W for discrete Markov chains in a N -dimensional state space, in order to characterise the statistics of the second largest eigenpair. This kind of Markov chain can be represented as a random walk on a graph. We remark here that the second largest eigenpair encodes non-equilibrium properties of a Markov process. Indeed, the inverse of the (absolute value) of the second largest eigenvalue represents the slowest relaxation time, whereas We remark that our method cannot provide the eigenvector statistic for x = c − 2 √ c − 1. Indeed, for this specific value of x, the probe eigenvector u is forced to correspond to the eigenvalue 2 √ c − 1, which retains its own eigenvector, thus artificially creating a degeneracy. Our method is based on the assumption of non-degeneracy of eigenvalues, so we are not able to give a result about eigenvectors in this marginal case. the associated second eigenvector is the non-equilibrium mode with the largest relaxation time.
We will then employ a full deflation, by setting x = λ 1 (W ) = 1. The evolution equation for the Markov chain states probability vector at time t, p(t), is given in terms of the matrix W by The transition matrix W is such that W ij ≥ 0 ∀(i, j) and i W ij = 1 ∀j. For an irreducible chain, the top right eigenvector of the matrix W corresponding to the Perron-Frobenius eigenvalue λ 1 = 1 represents the unique equilibrium distribution, i.e. v (1) = p eq . The matrix W is in general not symmetric. However, if the Markov process satisfies a detailed balance condition, i.e. W ij p eq j = W ji p eq i ∀(i, j), it can be symmetrised via a similarity transformation, yielding The symmetrised matrix W S and its deflated versionW S will be the target of our analysis: even though W S is not itself a Markov matrix since the columns normalisation constraint is lost, in view of the detailed balance condition W S has the same (real) spectrum as W , and its top eigenvector u is given in terms of the top right eigenvector of W , p eq , as It is actually well-known that the relation between the eigenvectors of W and those of W S holds in general and is not limited to the case of the top one.
We will consider the case of an unbiased random walk: the matrix W is then defined as where c ij represents the connectivity matrix and k j = i c ij is the degree of node j.
In this case, the top right eigenvector of W is proportional to the vector expressing the degree sequence: for our purposes, we choose the inverse of the mean degree as proportionality constant, i.e. p eq i = k i /(N c). The symmetrised matrix W S is expressed as with its top eigenvector being u where p(k) is the degree distribution of the connectivity matrix {c ij }.
In order to avoid isolated nodes and isolated clusters of nodes, we consider degree distributions with k min ≥ 2 and finite mean degree † †. We will provide a treatment for a generic distribution p(k) with the aforementioned properties and the analytical solution for the random regular connectivity case with degree distribution p(k) = δ k,c .

Replica analysis: second largest eigenpair of Markov transition matrices
We focus on the fully deflated symmetrised version of the Markov matrix W , that is where and u represents the top eigenvector of W S , normalised to N , i.e.
Here, c represents the mean degree, c = k . Our aim is to find the typical largest eigenvalue ofW S , which corresponds to the typical second largest eigenvalue of W S . In the next subsection, we will characterise the distribution of the top eigenvector ofW S , equivalent to the second eigenvector of W S .
We follow the same formalism illustrated in Section 4.1. The partition function reads (167) † † A suitable candidate could be a shifted Poissonian degree distribution with k min = 2, i.e.
with mean degree c =c + 2.
By expressing the delta function in (167) via its Fourier representation and employing the change of variableṽ i √ k i ← v i , the partition function becomes The square in the exponent of (168) can be linearised by a Hubbard-Stratonovich transform as in (47). The resulting partition function, where we rename theṽ i variables as v i to avoid cumbersome notation, reads The average w.r.t. the matrix ensemble ofW S reduces to averaging over the connectivity matrix C = {c ij }. By using the replica trick, we need to compute Henceforth, the derivation will exactly match the steps in Section 4.1. Therefore, we will just report the final equations, corresponding to (117) along with (118) and (119). By taking into account (164) and the existence of k min = 2, we find We remark that in the Markov case a bounded largest degree is not strictly necessary as the spectrum is always bounded. However, we will consider a k max for practical purposes. The self-consistency equation (171) along with the normalisation condition (172) and the orthogonality constraint (173) is solved by a population dynamics algorithm (See Section 7). The RRG connectivity case is analytically tractable, as shown in Section 6.2.
In analogy to Section 4.2, the density of the top eigenvector's component of the matrixW S , corresponding to the second largest eigenvector of W S , is given by where π(ω, h) satisfies the self-consistency equation (171), supplemented by the normalisation condition (172) and the ortogonality condition (173). Figure 5 compares the pdf of the second largest eigenvector's components obtained via population dynamics with results obtained via direct diagonalisation, for the unbiased random walk Markov matrix case with shifted Poisson degree distribution (k min = 2). We study both the low (c 6, left panel) and high (c 12, right panel) connectivity cases. In the c 6 case with k max = 12, we find λ 2 W S = 0.7456, within a 0.7% error w.r.t. the value λ 2,∞ = 0.7504 obtained by extrapolation from the direct diagonalisation data. In the c 12 case with k max = 22, we find λ 2 W S = 0.5530, within a 0.1% error w.r.t. the value λ 2,∞ = 0.5524 obtained by extrapolation from the direct diagonalisation data. As a reference point, the average value of the second largest eigenvalue in the RRG case with the same c is λ 2 (W S ) RRG = 0.5528. We notice that the agreement near the peak of the distribution is slightly worse for the low connectivity case: this is in agreement with the finding that finite-size effects are generally more pronounced for lower c (see also discussion in section 7.4).

Unbiased random walk on a RRG: second largest eigenpair statistics
For a random regular graph, for which p(k) = δ k,c , we note that the matrix W S reduces to implying that all results about the RRG adjacency matrix case stated in Sections 5.3 and 5.4 carry over to this case too, but with all eigenvalues rescaled by 1/c. As expected, λ 1 (W S ) RRG = 1, and the second largest eigenvalue corresponding to a N (0, 1)distributed eigenvector is . The spectral gap for this kind of Markov matrix as a function of c is then g(c)

The curse of orthogonality
With the exception of the unweighted adjacency matrix of a RRG, Eq. (117)supplemented with the conditions (118) and (119) -must be generally solved via a Population Dynamics algorithm, a Monte Carlo technique deeply rooted in the statistical mechanics of spin glasses [46,47]. The algorithm we use bears some similarity with the one employed in [21]. Here, we will highlight the main differences that stem from the presence of the orthogonality condition (119). We recall that the Eqs. (117), (118) and (119) refer to the case of full deflation, where we look at the top eigenpair of the deflated matrixJ (the second largest eigenpair of the matrix J).
Some observations are in order before sketching the algorithm. As we stated in [21], within the population dynamics algorithm the definition of the h variables in Eq. (117) is effectively converted into a stochastic linear update of h values. Its stability can only be achieved for λ = λ 1 J . For any λ > λ 1 J , the variables of type h will shrink to zero, whereas for λ < λ 1 J they will explode in norm. In our scenario, where we consider λ < λ 1 J , the recursion is thus a priori unstable, unless it is otherwise constrained. Therefore, if unconstrained, the population will never spontaneously evolve towards a stable regime, which would at the same time satisfy the conditions (118) and (119).
As anticipated in Section 4.3, this observation entails that the orthogonality condition (119) must be strictly enforced on-the-fly -by imposing a correction to the fields h, which once again have no fixed scale given by their update equation. Enforcing the constraint (119) is equivalent to looking for a self-consistent solution of (117) in a smaller, constrained space. Only once the condition (119) has been enforced, a new stable non-trivial fixed point arises, and the behaviour of the h-variables is similar to that in the top eigenvector case: for any value λ > λ 2 J , the variables h under iteration of the modified population dynamics algorithm shrink to zero, whereas for λ < λ 2 J they will explode in norm. Hence, Eq. (117) -taken together with the condition (119) -admits a stable, hence normalisable solution, such that Eq. (118) is naturally satisfied only for λ = λ 2 J : after the orthogonality correction has been enforced, the procedure we follow is then exactly identical to that used in [21].

The algorithm
Taking into account the observations made in section 7.1, we briefly sketch the algorithm in the case of full deflation.
Two pairs of (coupled) populations with N P members each {(a i , b i )} 1≤i≤N P and {(ω i , h i )} 1≤i≤N P are randomly initialised, taking into account that both a i and ω i must be larger than ζ, the upper edge of the support of the pdf p(K). We typically choose N P = 10 5 or larger. In what follows, the parameter λ is the candidate second largest eigenvalue of J, whereas λ 1 J is the average top largest eigenvalue of J. The first population is employed to solve the top eigenpair problem, and the other to solve the second eigenpair problem; the latter is constrained by results of the former due to the orthogonality constraint.
We therefore first run a short population dynamics simulation following Section 6 in [21] involving only the population {(a i , b i )} 1≤i≤N P to find the solution for the first eigenpair problem and the value λ 1 J . This first simulation acts as an equilibration phase for the fields contributing to the largest eigenpair. Then, for any suitable value of λ ∈ R < λ 1 J , the following steps are iterated until stable populations are obtained: and replace two randomly selected pairs (a j , b j ) and (ω j , h j ) where j ∈ {1, ..., N P } with the pairs a (new) , b (new) and ω (new) , h (new) .
(iv) Compute the components of the top eigenvector u and the candidate second largest eigenvector v. In order to create a sample estimate of the eigenvectors statistics corresponding to the two top eigenvalues, we initialise two empty vectors, where s = k j is exactly the "degree" drawn from p(s) in step (iv)(a) and used to build each component v j in step (iv)(c).
A sweep is completed when all the N P pairs (a j , b j ) and (ω j , h j ) have been updated at least once according to the steps above. The update of the pairs (a, b) is stable, thanks to the prior equilibration phase. The convergence is assessed by looking only at the first moments of the two vectors formed by the N P samples of the pairs (ω, h). The parameter λ is varied according to the behaviour illustrated in Section 7.1: starting from an initial "large" value λ < λ 1 J , it is then progressively decreased until a non trivial distribution for the h is achieved, in correspondence of the value λ = λ 2 J . Indeed, we observe that for any λ > λ 2 J , the h shrink to zero, whereas for any λ < λ 2 J , they blow up in norm.
Some comments are in order: • the condition expressed in (183) is a Gram-Schmidt orthogonalisation, taking place after every microscopic update of the fields; • the correction does not take place for components v j related to s = 0, as both v j and u j are zero; • in step (iv)(c), we can clearly see that the components u j and v j are coupled through their degree and the set of bond weights, as anticipated in Section 4.3. Indeed, for any j, the s i.i.d. realisations of the weights {K} s and the "local neighbourhood" that we dynamically create at every step (c) must be exactly the same for both u j and v j . In other words, both u j and v j must have the same update history.

Potential for simplifications in special cases
The steps (iv) and (v) of the algorithm are computationally heavy. We are able in some cases to simplify them.
• For adjacency matrices of RRGs, where the variables a, b and ω are constant, the correction (183) translates to forcing the mean of the h to be zero after every update. Both steps (iv)-(v) are then replaced by whereh indicates the sample mean of the h population.
• In the E-R case (both weighted and non-weighted), we take advantage of the fact that in the thermodynamic limit there is no statistical distinction between the cavity fields ω and h (respectively a and b) and the denominator and numerator in (182), (respectively in (181)), even in presence of the truncation of the Poissonian degree distribution †. Hence, we can consider just one couple of fields per species to represent a component, so we identify M =N P . Steps (iv) and (v) are then replaced by (iv) Compute eigenvectors u and v as (187) † Provided that the largest degree is reasonably large. The only difference between the distribution π(ω, h) and the distribution of the denominator and numerator of (121) can be observed because of the contribution coming from the largest degree, whose probability to occur is negligible. (v) Compute the correction as 7.4. Is the population dynamics algorithm really capturing the thermodynamic limit?
When no simplification can be used, as in the case of Markov matrices, the population dynamics algorithm can be relatively slow, due to the number of nested updates it requires. In these cases, we have therefore been often forced to consider a population size N P smaller than the values we would have typically wished (N P = O(10 5 ) or more). However, what may appear as a limitation at first sight turned out to be a blessing, in that it made us aware of an interesting interplay between the size N P of the population dynamics, and the size N of the graph whose spectral properties were to be reproduced.
Indeed, we have collected convincing evidence that population dynamics at finite N P does not really capture the thermodynamic limit N → ∞: for a given graph size N 1, there is an optimal size of the population N P = N P (N ) that best captures the spectral properties of that finite-size graph, and the accuracy between "theory" and numerical diagonalisation has a strongly non-monotonic behaviour as a function of N P . Similarly, a population of given size N P reproduces well spectral properties of graphs around a certain optimal size N , but its accuracy rapidly deteriorates if the graph size N is markedly different from N . Of course, the higher N P (e.g. in cases where it is possible to employ N P = O(10 5 ) or larger), the better the large N limit is captured (see e.g. the case in Fig. 2).
This intriguing phenomenon may be related to the existence of loops, which seem to be more relevant in the eigenvector problem than the spectral problem. Indeed, whatever N P is, the cavity fields of type ω and h will have common predecessors within their own species after ∼ ln(N P )/ ln(c − 1) updates. This implies the presence of loops in the population dynamics update history, which lead to correlations between different members of the population. Therefore, the assumption of population elements independently drawn from an ensemble, which underlies (37) (or equivalently (117)) is violated. That assumption in turn implements the notion that loops in the underlying graph that is being described will diverge in the thermodynamic limit.
To quantify this effect, we compare the cumulative distribution function (CDF) of the second eigenvector's components of Markov matrices with Poissonian shifted degree distribution, obtained via population dynamics at various N P , with the result from direct diagonalisation of matrices from the same ensemble at a given size N = 1000for both low and high mean degree.
In Figure 7, we assess the similarity of the two distributions using two figures of merit. The first (left) is the p-value of a 2-sample Kolmogorov-Smirnoff (KS) test: the larger the p-value, the strongest the evidence in favor of the hypothesis that the two distributions are the same. The second (right) is based on the analysis of a so-called quantile-quantile plot (Q-Q plot), which is the scatter plot of the quantiles of the two sets of data. Precisely, we focus on the slope m of the best fit regression line y = mx + b of the Q-Q plot, considered between the first and third quartile (respectively, the 0.25 and 0.75 quantiles), to limit spurious effects coming from the under-sampling of the tails. The slope m is directly proportional to the correlation coefficient between the quantiles of the two distributions, and m = 1 for identical distributions.
The existence of an optimal population size N P for a given graph size N -and the non-monotonic behaviour of the accuracy with N P -is quite evident in the left panels. The optimal value of N P (N ) is consistently identified by both figures of merit. However, the effect is more pronounced in the case of low connectivity (top row of Figure  7) -where finite size effects are indeed stronger -than in the case of high connectivity (bottom row of Figure 7).

Conclusions
In summary, we have developed a formalism to compute the statistics of the second largest eigenvalue and of the components of the corresponding eigenvector for some ensembles of sparse symmetric matrices, i.e. weighted adjacency matrices of graphs with Figure 7: Top panels: low mean degree case, c 6, reference matrix size N = 1000. The left panel shows the base-10 logarithm of the p-value of the KS two-sample test comparing the two empirical cdfs corresponding to different population sizes. We notice that the p-values are all rather low, yet there is a clear maximum value at N P 5000, and the nonmonotonic behaviour is quite pronounced. The right panel shows the slope m of the best-fit regression line of the Q-Q plot between the 25% and 75% quantiles, for various population sizes. The closer m is to 1, the better the agreement. The plot confirms again that the best agreement with our reference distribution is obtained with N P 5000. Bottom panels: high mean degree case, c 12, reference matrix size N = 1000. On the left, we show the p-value of the KS two-sample test against N P in linear y-scale. The curve is much flatter than the low-c case, and the p-values are all significant, suggesting a high level of similarity between the two distributions throughout the full range of N P . On the right, we plot the slope m of the best fit regression line of the Q-Q plot between the 25% and 75% quantiles, for various population sizes. For this figure of merit, we again observe a rather flat value of the slope between N P 2000 and N P 6000, where m 1 (within a 0.2% error). At high c, we indeed observe negligible finite size effects in the direct diagonalisation samples at different sizes N , and this phenomenon seems to be present also in the population dynamics simulations. finite mean connectivity. By assuming that the top eigenpair is known, we show that for a given matrix, computing the second largest eigenpair is equivalent to computing the top eigenpair of a deflated matrix, obtained by subtracting from the original matrix the dense matrix representing a rank-one perturbation proportional to the projector onto its first eigenstate. As in [21], the search for the top eigenpair of the deflated matrix is then transformed into the optimisation of a quadratic Hamiltonian on a sphere: introducing the associated Gibbs-Boltzmann distribution and a fictitious inverse temperature β, the top eigenvector represents the ground state of the system, reached in the limit β → ∞. In order to extract this limit, we have employed two Statistical Mechanics methods, cavity and replicas. We started analysing the case of a single-instance matrix within the cavity framework, introducing a new cavity formulation that allows for the inclusion of hard constraints.
The single-instance cavity method easily leads to recursion equations, which represent the essential ingredient to obtain the solution of the problem in the thermodynamic limit. We also derive the exact same equations from a completely alternative replica derivation, confirming the equivalence of the two methods in the thermodynamic limit. We employed an improved population dynamics algorithm to solve the stochastic recursion (37) complemented by the conditions (38) and (39), (or equivalently (117) along with (118) and (119)). We found that the convergence of the algorithm is driven not only by the largest eigenvalue of the deflated matrix (i.e. the second largest eigenvalue of the original matrix) but, most essentially, by the fact that the orthogonality condition (39) (or equivalently (119)) be correctly enforced. Some ensembles permit simplifications of the algorithm used to enforce orthogonality, which we exploited to speed up convergence.
The simulations show excellent agreement between the theory and the direct diagonalisation of large matrices, and allow us to unpack the contributions to the average density of the second eigenvector's components coming from nodes of different degrees.
Our study clearly demonstrates that -in contrast to beliefs commonly held in the community -population dynamics at finite N P is fundamentally incapable of analysing properties representing the thermodynamic limit behaviour. This discovery is in some sense due to the fact that finite size effects are much stronger for eigenvectors than for eigenvalues (in particular for matrices without random edge weights). That finite population size effects are quantitatively related to finite size effects is, in retrospect, not really surprising, given the clear analogy existing between the emergence of correlations in population values -through loops of common ancestors of population updatesand common ancestors created through loops in random graphs of finite size, in which the scaling of loop lengths with population and graph size follows basically the same logarithmic law.
In the case of the RRG adjacency matrix, we also analytically studied the pdf of the components of the top eigenvector of the deflated matrix as the deflation parameter is continuously changed, showing the abrupt change of the solution as soon as the deflation parameter becomes larger than the spectral gap of the Kesten-McKay distribution.
Lastly, we applied our formalism to sparse Markov matrices representing unbiased random walks on a network, for which the second largest eigenpair plays an important role encoding non-equilibrium properties.