Analytic solution of the resolvent equations for heterogeneous random graphs: spectral and localization properties

The adjacency matrix of random graphs stores the interactions between the constituents of large complex systems [1] ranging from physical and biological to social and technological systems. The empirical spectral density of the adjacency matrix and the localization of its eigenvectors are important to understand algorithms for node centrality [2, 3] and community detection [4, 5], as well as the interplay between the structure of networks and dynamical processes on them. In fact, the leading eigenpair of the adjacency matrix governs the spreading of diseases [6, 7], the synchronization transition [8, 9], and the linear stability of large complex systems [10–13]. In condensed matter physics, models defined on random graphs represent mean-field versions of finite-dimensional lattices which mimic the effects of finite coordination number. The spin-glass transition and Anderson localization have been intensively investigated on random graph structures over the past years [14, 15], and they continue to attract a lot of interest [16–20].

The spectral and localization properties of the random adjacency matrix are determined by the resolvent matrix G . The average of its diagonal elements G_ii yields the empirical spectral density, while the average of $|G_{ii}|^2$ gives the inverse participation ratio (IPR) [21–23], which characterizes the volume of the eigenvectors. The full probability density of G_ii fulfils a system of distributional equations, derived in a series of fundamental works [24–26] using the cavity and the replica methods of spin-glass theory (see [27] for a review of these techniques). The resolvent distributional equations are exact on locally tree-like random graphs [28] and they provide a solid framework to investigate the spectral properties of sparse and heterogeneous random graphs. Heterogeneity broadly refers to local fluctuations in the graph structure, such as randomness in the degrees or in the interaction strengths between the nodes, while sparseness means that the average degree is finite (the degree counts the number of edges attached to a node). The numerical solutions of the resolvent equations have led to a profusion of results for the spectral and localization properties of random graphs [22, 29, 30] with different topological features, including short loops [31], modularity [32], and degree-degree correlations [33]. Analogous resolvent equations describe the spectral properties of stochastic matrices on random graphs [23, 34].

Despite the resolvent distributional equations have led to a tremendous progress in the field, they admit analytic solutions only for sparse regular graphs [28, 31], whose local structure is homogeneous, and for high-connectivity random graphs [25], where the mean degree is infinitely large and the graph becomes homogeneous on account of the law of large numbers. In each case, the distribution of G_ii is a Dirac-δ and the spectral density follows from a simple algebraic equation for the average resolvent. In a recent paper [35], the resolvent equations for the configuration model of random graphs with a geometric degree distribution have been studied in the high-connectivity limit [1, 36, 37]. In this case, the average resolvent fulfils a transcendental equation and the spectral density diverges at centre of the spectrum [35]. These analytic findings are interesting for at least two reasons. First, they imply that the spectral density of high-connectivity random graphs is not generally given by the Wigner law of random matrix theory [38], but it explicitly depends on the degree distribution. Indeed, as rigorously proven in the Wigner universality only holds for degree distributions that become sharply peaked in the high-connectivity limit. Second, these findings hint the existence of a rich and nontrivial family of analytic solutions of the resolvent equations, sandwiched between the sparse and the dense regime, in which the average degree is large but the network heterogeneities are still relevant for the spectral properties. The analytic results in [35] are limited, however, to a geometric degree distribution.

In this paper we generalize the results of [35] and we extract the analytic solution of the resolvent equations for random graphs with arbitrary degree distributions in the high-connectivity limit. To our knowledge, this is the first example of a full analytic solution for the probability density of G_ii for undirected random graphs with an heterogeneous structure. We show that the spectral density, the IPR, and the distribution of the local density of states (LDOSs) are fully determined by the choice of the degree distribution. We present explicit results for a negative binomial degree distribution, in which the variance of the degrees is controlled by a single parameter $0 \lt \alpha \lt \infty$ that enables to interpolate between homogeneous ($\alpha \rightarrow \infty$) and strongly heterogeneous graphs ($\alpha \rightarrow 0$). In this way, we are able to thoroughly investigate the impact of degree fluctuations on the spectral and localization properties. We show that the spectral density undergoes a transition as the degree fluctuations become stronger. For α > 1, the spectral density is a regular function, whereas it displays either a power-law or a logarithmic divergence at the zero eigenvalue provided $\alpha \in (0,1)$ or α = 1, respectively. In the regime of weak degree fluctuations ($1 \ll \alpha \lt \infty$), the spectral density has a finite support and large complex systems interacting through the underlying adjacency matrix can be found in a (linearly) stable state. From the analytic results for the IPR and the distribution of the LDOSs, we show that all eigenvectors of the adjacency matrix are extended for any α. In particular, the distribution of the LDOSs at the zero eigenvalue exhibits a power-law tail with exponent α + 1, which emphasizes the prominent role of the degree fluctuations and clarifies the singular behaviour of the spectral density. We support our theoretical findings by comparing them with numerical diagonalizations of large adjacency matrices and, as a byproduct, we show that the adjacency matrix of a single graph instance can be decomposed for large average degree as a product between a Gaussian random matrix [38] and the square root of the degree matrix. Such decomposition provides a straightforward way to sample the adjacency matrix of the configuration model of networks in the high-connectivity limit.

The paper is organized as follows. In the next section we introduce the random graph model and the resolvent equations for its adjacency matrix. In section 3 we derive the distributional equations and the analytic expression for the probability density of the resolvent in the high-connectivity limit using the law of large numbers. In section 4 we discuss explicit results for the spectral density, the IPR, and the distribution of the LDOSs in the case of a negative binomial degree distribution. We present a summary and a discussion of our results in section 5, and we provide a more rigorous derivation of the analytic solution of the resolvent distributional equations in the appendix.

We consider a simple and undirected random graph with N nodes. The graph structure is specified by the set of binary random variables $\{c_{ij} \}$ ($i,j = 1,\dots,N$), in which $c_{ij} = c_{ji} = 1$ if there is an undirected edge between nodes i and j (i ≠ j), and $c_{ij} = 0$ otherwise. In addition, we associate a symmetric coupling strength $J_{ij} = J_{ji} \in \mathbb R$ to each edge $i\leftrightarrow j$. The degree $k_i = \sum_{j = 1}^N c_{ij}$ of a node i gives the number of nodes attached to i, and the degree distribution

$$\begin{equation} p_k = \lim_{N \rightarrow \infty} \frac{1}{N}\sum_{i=1}^N\delta_{k,k_i} \end{equation} \tag{ 1 }$$

gives the fraction of nodes with degree k. The average degree stands as

$$\begin{equation} c = \sum_{k=0}^\infty k p_k. \end{equation} \tag{ 2 }$$

We study the spectral properties of the N × N adjacency random matrix A , with elements

$$\begin{equation} A_{ij} = c_{ij}J_{ij}, \end{equation} \tag{ 3 }$$

where the coupling strengths J_ij are, apart from the symmetry constraint $J_{ij} = J_{ji}$, independently and identically distributed random variables drawn from a distribution p_J with mean zero and standard deviation $J/\sqrt{c}$. The nonzero values of $\{c_{ij} \}$ are randomly assigned following the configuration model of networks [36, 37, 39], in which a single graph instance is uniformly chosen at random from the set of all random graphs with a given degree sequence $k_1,\ldots,k_N$ sampled from p_k . The configuration model allows us to fix the degree distribution from the outset and study its impact on the spectral properties.

The adjacency matrix A has a complete set $\{\boldsymbol{v}_\mu\}_{\mu = 1,\ldots,N}$ of orthonormal eigenvectors that fulfil

$$\begin{equation} \boldsymbol{A} \boldsymbol{v}_\mu = \lambda_\mu \boldsymbol{v}_\mu, \end{equation} \tag{ 4 }$$

with $\{\lambda_\mu\}_{\mu = 1,\ldots,N}$ the set of eigenvalues. The empirical spectral density of A reads

$$\begin{equation} \rho(\lambda) = \lim_{N \rightarrow \infty} \frac{1}{N}\sum_{\mu=1}^N \delta(\lambda - \lambda_\mu). \end{equation} \tag{ 5 }$$

The IPR of an eigenvector $\boldsymbol{v}_\mu$ is defined as

$$\begin{equation} Y_{\mu} = \sum_{i=1}^N (v_{\mu,i})^4, \end{equation} \tag{ 6 }$$

where $v_{\mu,i}$ is the ith component of $\boldsymbol{v}_\mu$. The IPR distinguishes between localized and extended eigenvectors in the large N limit. The components of a localized eigenvector are nonzero on a finite number of nodes and the corresponding IPR is of order $\mathcal{O}(N^0)$, whereas an extended eigenvector is spread over a finite fraction of nodes and the IPR vanishes as $\mathcal{O}(1/N)$ in the large N limit. Since the extent of the eigenvectors typically depends on the corresponding eigenvalue, it is sensible to introduce the eigenvalue-dependent IPR [21–23, 40]

$$\begin{equation} \mathcal{I}(\lambda) = \lim_{N \rightarrow \infty} \frac{ \sum_{\mu=1}^N\delta(\lambda-\lambda_\mu)Y_\mu}{ \sum_{\mu=1}^N\delta(\lambda-\lambda_\mu)}, \end{equation} \tag{ 7 }$$

which is the average of Y_µ over all eigenvectors in an infinitesimal spectral window around λ.

The spectral properties of A follow from the resolvent matrix

$$\begin{align} \boldsymbol{G}(z) = (\boldsymbol{I} z-\boldsymbol{A})^{-1}, \end{align} \tag{ 8 }$$

where I is the N × N identity matrix and $z = \lambda- i\epsilon$ lies in the lower complex half-plane. The diagonal elements of G determine the spectral density and the eigenvalue-dependent IPR for a finite regularizer $\epsilon$ > 0 according to [40]

$$\begin{align} & \rho_{\epsilon}(\lambda)= \frac{1}{\pi} \lim_{N \rightarrow \infty} \frac{1}{N} \sum_{i=1}^N \mathrm{Im}G_{ii}(z), \end{align} \tag{ 9 }$$

$$\begin{align} &\mathcal{I}_{\epsilon}(\lambda) = \frac{\epsilon}{\pi\rho_{\epsilon}(\lambda)} \lim_{N \rightarrow \infty} \frac{1}{N} \sum_{i=1}^N |G_{ii}(z)|^2. \end{align} \tag{ 10 }$$

Working with finite $\epsilon$ amounts to replace the Dirac-δ distributions appearing in equations (5) and (7) by Cauchy distributions with a scale parameter $\epsilon$ [25, 26]. The spectral observables $\rho(\lambda)$ and $\mathcal{I}(\lambda)$ are reconstructed by taking the limit $\epsilon \rightarrow 0^{+}$ in equations (9) and (10).

By introducing the joint probability density of the real and imaginary parts of $G_{ii}(z)$,

$$\begin{equation} \mathcal{P}_{z}(g) = \lim_{N\to \infty}\frac{1}{N}\sum_{i=1}^N \delta\left[g- G_{ii}(z)\right], \end{equation} \tag{ 11 }$$

the spectral density and the IPR can be written as

$$\begin{align} & \rho_{\epsilon}(\lambda)= \frac{1}{\pi} \mathrm{Im}\Bigg[\int_{{\mathbb H}^+} dg \mathcal{P}_{z}(g) g\Bigg], \end{align} \tag{ 12 }$$

$$\begin{align} &\mathcal{I}_{\epsilon}(\lambda) = \frac{\epsilon}{\pi\rho_{\epsilon}(\lambda)}\int_{{\mathbb H}^+} dg \mathcal{P}_{z}(g) |g|^2. \end{align} \tag{ 13 }$$

The symbol $\mathbb{H}^+$ represents the complex upper half-plane and we have introduced the shorthand notation $dg = d \mathrm{Re} g \, d \mathrm{Im} g$.

We see that $\rho_{\epsilon}(\lambda)$ and $\mathcal{I}_{\epsilon}(\lambda)$ are determined by the moments of $\mathcal{P}_z(g)$. In the limit $N \rightarrow \infty$, the local structure around a randomly chosen node of a graph drawn from the configuration model converges to a tree [28], and the probability of finding a short loop in a finite neighbourhood of the node in question goes to zero. This property implies that the resolvent diagonal elements for a single graph instance fulfil the equations [33]

$$\begin{equation} G_{ii}(z) = \frac{1}{z - \sum_{j\in \partial_i}J_{ij}^2 G_{jj}^{(i)}(z)} \quad (i=1,\dots,N), \end{equation} \tag{ 14 }$$

where $\partial_i$ is the set of nodes adjacent to i. The complex variable $G_{jj}^{(i)}$ is the jth-diagonal element of the resolvent on a graph in which node $i \in \partial_j$ and all its incident edges have been deleted [41]. The variables $\{G_{jj}^{(i)} \}$ are determined from the fixed-point solutions of the so-called cavity equations

$$\begin{equation} G_{jj}^{(i)}(z) = \frac{1}{z - \sum_{\ell \in \partial_j\backslash i}J_{j\ell}^2 G_{\ell\ell}^{(j)}(z) } \quad i \in \partial_j, \end{equation} \tag{ 15 }$$

where $\partial_j \backslash i$ is the set of nodes connected to j excluding i. The total number of cavity variables $\{G_{jj}^{(i)} \}$, defined on the edges of the graph, is $\sum_{i = 1}^N k_i$. The solutions of equation (15) lead to approximations for the resolvent diagonal elements of single graph instances when N is large. In the limit $N \rightarrow \infty$, equations (14) and (15) become exact and it is more convenient to work with the distributions of $G_{ii}(z)$ and $G_{jj}^{(i)}(z)$. Given that both sides of equation (14) are equal in distribution, the probability density $\mathcal{P}_z(g)$ is determined from

$$\begin{equation} \mathcal{P}_{z}(g) = \sum_{k=0}^\infty p_k \int_{{{\mathbb{H}}^+}} \bigg[ \prod_{\ell=1}^{k}d g_\ell \mathcal{Q}_{z}(g_\ell) \bigg] \int_{\mathbb{R}} \bigg[\prod_{\ell=1}^{k} d {J_{\ell}} p_J(J_\ell)\bigg] \delta\Bigg (g - \frac{1}{z - \sum_{\ell=1}^{k}J_{\ell}^2 g_\ell}\Bigg), \end{equation} \tag{ 16 }$$

where $\mathcal{Q}_{z}(g)$ is the probability density of the cavity variables $G_{jj}^{(i)}(z)$, defined as

$$\begin{equation} \mathcal{Q}_{z}(g) = \lim_{N \rightarrow \infty} \frac{\sum_{i,j=1}^N c_{ij} \delta\left[g- G_{ii}^{(j)}(z)\right] }{\sum_{i,j=1}^N c_{ij} }, \end{equation} \tag{ 17 }$$

which solves the self-consistent equation

$$\begin{equation} \mathcal{Q}_z(g) = \sum_{k=1}^\infty \frac{k}{c}p_k \int_{\mathbb{H}^{+}} \bigg[\prod_{\ell=1}^{k-1}d g_\ell \mathcal{Q}_z(g_\ell) \bigg] \int_{\mathbb{R}} \bigg[\prod_{\ell=1}^{k-1} d {J_{\ell}} p_J(J_\ell)\bigg] \delta\Bigg (g - \frac{1}{z - \sum_{\ell=1}^{k-1}J_{\ell}^2 g_\ell}\Bigg). \end{equation} \tag{ 18 }$$

Equations (16) and (18) are the distributional version of the resolvent equations. Once we solve equation (18) and find a fixed-point solution for $\mathcal{Q}_z(g)$, the probability density $\mathcal{P}_{z}(g)$ of the diagonal elements of the resolvent follows from equation (16). As we will consider the high-connectivity limit $c \rightarrow \infty$, it is interesting to introduce the joint probability density $\mathcal{W}_{z}(s)$ of the complex variable

$$\begin{equation} S(z) \stackrel{d}{=} \sum_{\ell=1}^{k}J_{\ell}^2 g_\ell, \end{equation} \tag{ 19 }$$

which consists of a sum of independent and identically distributed random variables. The distribution $\mathcal{P}_{z}(g)$ is written in terms of $\mathcal{W}_{z}(s)$ as

$$\begin{equation} \mathcal{P}_{z}(g) = \int_{\mathbb{H}^{+}} d s \mathcal{W}_{z}(s) \delta\Bigg (g - \frac{1}{z - s}\Bigg). \end{equation} \tag{ 20 }$$

In the context of tight-binding models for the diffusion of an electron on a graph [14, 42], S(z) is known as the self-energy and its distribution $\mathcal{W}_{z}(s)$ plays an important role in the study of the Anderson localization transition. We note from equation (20) that the moments of the resolvent diagonal elements $G_{ii}(z)$ are determined by the distribution $\mathcal{W}_{z}(s)$ of the self-energy. For instance, equations (12) and (13) become

$$\begin{align} & \rho_{\epsilon}(\lambda)= \frac{1}{\pi} \mathrm{Im}\Bigg[\int_{{\mathbb H}^+} ds \frac{\mathcal{W}_z(s)}{z-s}\Bigg], \end{align} \tag{ 21 }$$

$$\begin{align} &\mathcal{I}_{\epsilon}(\lambda) = \frac{\epsilon}{\pi\rho_{\epsilon}(\lambda)}\int_{{\mathbb H}^+} ds \frac{\mathcal{W}_z(s)}{|z-s|^2}. \end{align} \tag{ 22 }$$

The exact equations (16) and (18), albeit having a complicated structure, represent a major step in our understanding of the spectral properties of random graphs, since they can be solved numerically using a Monte-Carlo iterative method called population dynamics [25–27]. Below we present an analytic solution of these equations for $c \rightarrow \infty$ and arbitrary degree distributions.

In this section we present a straightforward approach, based on the law of large numbers, to solve the resolvent equations and determine the spectral and localization properties of random graphs in the high-connectivity limit $c \rightarrow \infty$. In the appendix, we discuss a more rigorous derivation based on the characteristic functions of the probability densities $\mathcal{P}_{z}(g)$ and $\mathcal{Q}_{z}(g)$.

Our starting point are the cavity equations (14) and (15) for a single graph instance, expressed in terms of the self-energy

$$\begin{equation} S_{i}(z) = \sum_{j \in \partial_i} J_{i j }^2 G_{j j}^{(i)}(z) \end{equation} \tag{ 23 }$$

as follows

$$\begin{equation} G_{ii}(z) = \frac{1}{z - S_{i}(z)}, \qquad G_{jj}^{(i)}(z) = \frac{1}{z - S_{j}(z) + J_{ij}^2 G_{ii}^{(j)}(z) }. \end{equation} \tag{ 24 }$$

In the high-connectivity limit $c \rightarrow \infty$, $S_{i}(z)$ is a sum of a large and random number k_i of independent and identically distributed random variables. By the law of large numbers, $S_{i}(z)$ is asymptotically given by

$$\begin{equation} S_{i}(z) \stackrel{c \rightarrow \infty}{\longrightarrow} \kappa_i J^{2} \langle G \rangle , \end{equation} \tag{ 25 }$$

where $\kappa_i = k_i/c$ is the rescaled degree of node i and the expectation value of $G_{ii}^{(j)}(z)$ is defined as

$$\begin{equation} \langle G \rangle = \int_{\mathbb H^+} dg \mathcal{Q}_{z}(g) g. \end{equation} \tag{ 26 }$$

Therefore, the spatial fluctuations of $S_{i}(z)$ are solely governed by κ_i in the limit $c \rightarrow \infty$. By assuming that the empirical distribution of $\kappa_1,\dots,\kappa_N$ converges to $\nu(\kappa)$ as $c \rightarrow \infty$,

$$\begin{equation} \nu(\kappa) = \lim_{c\to\infty} \sum_{k=0}^\infty p_k\delta\Big(\kappa - \frac{k}{c}\Big), \end{equation} \tag{ 27 }$$

the probability density $\mathcal{W}_z(s)$ of $S_{i}(z)$ is obtained by the change of variables set by equation (25), namely

$$\begin{equation} \mathcal{W}_z(s) = \frac{1}{J^2 \mathrm{Im} \langle G \rangle } \nu \left( \frac{\mathrm{Im} s}{J^2 \mathrm{Im} \langle G \rangle } \right) \delta \left[\mathrm{Re} s - \frac{\mathrm{Re} \langle G \rangle }{\mathrm{Im} \langle G \rangle } \mathrm{Im} s \right]. \end{equation} \tag{ 28 }$$

We note that $\mathcal{W}_z(s)$ depends itself on the first moment $\langle G \rangle$ of $G_{ii}^{(j)}(z)$. This is computed by substituting the large c behaviour of $G_{ii}^{(j)}(z)$,

$$\begin{equation} G_{ii}^{(j)}(z) \stackrel{c \rightarrow \infty}{\longrightarrow} \frac{1}{z - \kappa_i J\,^{2} \langle G \rangle }, \end{equation} \tag{ 29 }$$

in the definition of the probability density $\mathcal{Q}_z(g)$ of $G_{ii}^{(j)}(z)$, equation (17), leading to the self-consistent equation

$$\begin{equation} \mathcal{Q}_z(g) = \int_{0}^{\infty} d \kappa \, \nu (\kappa) \, \kappa \, \delta \left( g - \frac{1}{z - \kappa J^2 \langle G \rangle} \right). \end{equation} \tag{ 30 }$$

The fixed-point equation for the first moment $\langle G \rangle$ readily follows from the above expression

$$\begin{equation} \langle G \rangle = \int_{0}^{\infty} d \kappa \frac{ \nu (\kappa) \kappa }{z - \kappa J^2 \langle G \rangle}. \end{equation} \tag{ 31 }$$

Equation (28) is one of our main analytic results, since the distribution $\mathcal{W}_z(s)$ of the self-energies determines all moments of the diagonal elements of the resolvent. Despite the fact we have considered the high-connectivity limit $c \rightarrow \infty$, the distribution $\mathcal{W}_z(s)$ retains information about the degree fluctuations through $\nu(\kappa)$. When the tail of the degree distribution p_k decays fast enough, such as in the case of regular and Erdös–Rényi random graphs [39], the rescaled degree distribution is given by $\nu(\kappa) = \delta(\kappa-1)$, and $\mathcal{W}_z(s)$ reduces to

$$\begin{equation} \mathcal{W}_z(s) = \delta \left( \mathrm{Im} s - J^2 \mathrm{Im} \langle G \rangle \right) \delta \left( \mathrm{Re} s - J^2 \mathrm{Re} \langle G \rangle \right). \end{equation} \tag{ 32 }$$

The above class of solutions describes homogeneous random graphs [35], in which the the self-energy $S_i(z)$ is equal to its mean value $J^2 \langle G \rangle$ and the spectral density is given by the Wigner law

$$\begin{equation} \rho_{\mathrm{w}} (\lambda) = \frac{\sqrt{4 J^2 - \lambda^2 }}{2 \pi J^2} \, \boldsymbol{1}_{(-2J ,2J)}(\lambda), \end{equation} \tag{ 33 }$$

where $\boldsymbol{1}_{\mathcal{A}}(x)$ denotes the indicator function, i.e. $\boldsymbol{1}_{\mathcal{A}}(x) = 1$ if $x \in \mathcal{A}$, and $\boldsymbol{1}_{\mathcal{A}}(x) = 0$ otherwise.

Let us derive some consequences of equation (28). By inserting equation (28) in equations (21) and (22), we obtain

$$\begin{equation} \rho_{\epsilon}(\lambda) = \frac{1}{\pi} \mathrm{Im}\Bigg[ \int_{0}^{\infty} d\kappa\frac{\nu(\kappa)}{z- \kappa J^2 \langle G \rangle }\Bigg] \end{equation} \tag{ 34 }$$

and

$$\begin{equation} \mathcal{I}_{\epsilon} (\lambda) = \frac{\epsilon}{\pi \rho_{\epsilon}(\lambda)} \int_{0}^{\infty} d\kappa\frac{\nu(\kappa)}{|z- \kappa J^2 \langle G \rangle|^2}. \end{equation} \tag{ 35 }$$

We can also derive the analytic expression for the joint distribution $\mathcal{P}_z(g)$ of the diagonal elements of the resolvent. Let G(z) and S(z) be independent complex random variables distributed according to $\mathcal{P}_z(g)$ and $\mathcal{W}_z(s)$, respectively, then equation (24) entails

$$\begin{equation} G(z) \stackrel{d}{=} \frac{1}{z-S(z)}. \end{equation} \tag{ 36 }$$

By making a two-dimensional change of variables and using equation (28), we find

$$\begin{equation} \mathcal{P}_z(g) = \frac{1}{J^2 |g|^4 \mathrm{Im} \langle G \rangle} \nu \left( \frac{ \mathrm{Im} \left( z - g^{-1} \right) }{J^2 \mathrm{Im} \langle G \rangle } \right) \delta \left[ \mathrm{Re} \left( z - g^{-1} \right) - \frac{\mathrm{Re} \langle G \rangle }{\mathrm{Im} \langle G \rangle } \mathrm{Im} \left( z - g^{-1} \right) \right]. \end{equation} \tag{ 37 }$$

The above equation determines how the distribution of the diagonal part of the resolvent depends on the distribution $\nu(\kappa)$ of rescaled degrees. The object $\mathcal{P}_z(g)$ contains much information about the spectral and localization properties of the adjacency matrix [20, 23, 43]. For instance, by marginalizing $\mathcal{P}_z(g)$ with respect to $\mathrm{Re} g$, we can calculate the empirical distribution of $\{\mathrm{Im} G_{ii} \}_{i = 1,\dots,N}$, which is essentially the (regularized) LDOSs [20, 40, 43]

$$\begin{equation} \rho_i(z) = \frac{1}{\pi} \mathrm{Im} G_{ii}(z) = \frac{1}{\pi} \sum_{\mu=1}^N \frac{\epsilon |v_{\mu,i} |^2 }{\left( \lambda -\lambda_{\mu} \right)^2 + \epsilon^2 } \quad (i=1,\dots,N). \end{equation} \tag{ 38 }$$

In the limit $\epsilon \rightarrow 0^{+}$, the empirical distribution of $y_i = \mathrm{Im} G_{ii}(z)$ ($i = 1,\dots,N$) characterizes the spatial fluctuations of the eigenvector amplitudes $|v_{\mu,i}|^2$ corresponding to the eigenvalues around λ. By integrating equation (37) over $\mathrm{Re} g$, we find the expression for $|\lambda|\gt0$

$$\begin{align} \begin{gathered} {P_z}(y) = \left\{ {{\omega _ + }(y){\mkern 1mu} \nu \left( {\frac{{y{{\left[ {x_ + ^2(y) + {y^2}} \right]}^{ - 1}} - }}{{{J^2}{\mkern 1mu} {\text{Im}}\langle G\rangle }}} \right)} \right. \hfill \\ \left. {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + {\omega _ - }(y){\mkern 1mu} \nu \left( {\frac{{y{{\left[ {x_ - ^2(y) + {y^2}} \right]}^{ - 1}} - }}{{{J^2}{\mkern 1mu} {\text{Im}}\langle G\rangle }}} \right)} \right\}{{\mathbf{1}}_{(0,{y_{\text{e}}})}}(y), \hfill \\ \end{gathered} \end{align} \tag{ 39 }$$

where the support of $P_z(y)$ is determined by

$$\begin{equation} y_{\mathrm{e}} = \frac{|\mathrm{Re} \langle G \rangle | + \sqrt{\left( \mathrm{Re} \langle G \rangle \right)^2 + \left( \mathrm{Im} \langle G \rangle \right)^2}} {2 \left( |\lambda| \mathrm{Im} \langle G \rangle + \epsilon \, |\mathrm{Re} \langle G \rangle| \right) }. \end{equation} \tag{ 40 }$$

The functions $\omega_{\pm}(y)$ and $x_{\pm}(y)$ are defined as

$$\begin{equation} \omega_{\pm}(y) = \frac{1}{J^{2} \, \mathrm{Im} \langle G \rangle \Big{|} x_{\pm}^2(y) - y^2 + 2 \frac{\mathrm{Re} \langle G \rangle }{\mathrm{Im} \langle G \rangle } \, y \, x_{\pm}(y) \Big{|} } \end{equation} \tag{ 41 }$$

and

$$\begin{equation} x_{\pm}(y) = \frac{1 \pm \sqrt{1 - 4 \left( \lambda + \frac{\mathrm{Re} \langle G \rangle }{\mathrm{Im} \langle G \rangle } \epsilon \right) \left[ \left( \lambda + \frac{\mathrm{Re} \langle G \rangle }{\mathrm{Im} \langle G \rangle } \epsilon \right) y^2 - \frac{\mathrm{Re} \langle G \rangle }{\mathrm{Im} \langle G \rangle } y \right]}} {2 \left( \lambda + \frac{\mathrm{Re} \langle G \rangle }{\mathrm{Im} \langle G \rangle } \epsilon \right) }. \end{equation} \tag{ 42 }$$

Equation (39) shows that $P_z(y)$ has a finite support for $|\lambda| \gt0$ regardless of the shape of the distribution ν. By setting λ = 0 in equation (37) and then integrating over $\mathrm{Re} g$, we obtain

$$\begin{equation} P_z(y) = \frac{1}{J^2 \mathrm{Im} \langle G \rangle y^2 } \nu \left( \frac{- \epsilon + y^{-1} }{J^2 \mathrm{Im} \langle G \rangle } \right) \quad (\lambda=0). \end{equation} \tag{ 43 }$$

In contrast to equation (39), the distribution $P_z(y)$ at λ = 0 can have an infinite support depending on the choice of the rescaled degree distribution ν.

The analytic expressions of the previous section are valid for any distribution of rescaled degrees $k_i/c$ that converges to $\nu(\kappa)$ as $c \rightarrow \infty$. Here we discuss results for random graphs with the negative binomial degree distribution [44]

$$\begin{equation} p_k^{\mathrm{(b)}} = \frac{\Gamma (\alpha + k)}{\Gamma(\alpha)}\frac{1}{k!}\left(\frac{c}{\alpha}\right)^k \frac{1}{(1 + \frac{c}{\alpha})^{\alpha +k}}, \end{equation} \tag{ 44 }$$

where $\Gamma(x)$ is the Gamma function and $0 \lt \alpha \lt \infty$ is a continuous parameter. The variance $\sigma_{b}^{2}$ of $p_k^{\mathrm{(b)}}$ is related to α as follows

$$\begin{align} \sigma_{b}^{2} = c + \frac{c^2}{\alpha}. \end{align} \tag{ 45 }$$

In a previous work [35], we have shown that spectral density of the configuration model does not converge to the Wigner law if the relative variance of the degree distribution does not vanish as $c \rightarrow \infty$. The negative binomial degree distribution provides a controllable way to investigate the effect of degree fluctuations on the spectral and localization properties of random graphs by varying a single parameter. In fact, given that

$$\begin{equation} \lim_{c \rightarrow \infty} \frac{\sigma_{b}^{2}}{c^2} = \frac{1}{\alpha}, \end{equation} \tag{ 46 }$$

by changing α we are able to explore the entire range of degree fluctuations for $c \rightarrow \infty$. The limit $\alpha \rightarrow \infty$ corresponds to homogeneous random graphs, whose spectral properties are governed by random matrix theory [38], whereas the limit $\alpha \rightarrow 0$ characterizes random graphs with strongly heterogeneous degrees. The geometric degree distribution is recovered for α = 1 [35]. Inserting equation (44) in equation (27), we obtain the analytic form of $\nu(\kappa)$

$$\begin{align} \nu_{b} (\kappa) = \frac{\alpha^\alpha \kappa^{\alpha - 1} e^{-\alpha\kappa}}{\Gamma(\alpha)}. \end{align} \tag{ 47 }$$

The above expression is the only input to the general formulae of the previous section, from which we can derive several analytic results as a function of α.

4.1. Spectral density

In the high-connectivity limit, the empirical distribution $\nu(\kappa)$ of rescaled degrees determines the spectral density $\rho(\lambda)$, as the latter is given by the free multiplicative convolution of $\nu(\kappa)$ with the Wigner law $\rho_{\mathrm{w}}(\lambda)$. This rigorous result, proven in [45], essentially means that the adjacency matrix can be decomposed in the limit $c \rightarrow \infty$ as a product of X and D , where D is the degree matrix with elements $D_{ij} = \kappa_i \delta_{ij}$, and X is a random matrix in which the diagonal entries are zero and the off-diagonal elements are independent random variables drawn from a Gaussian distribution with mean zero and variance $J^2/N$. The product $\boldsymbol{X}\boldsymbol{D}$, however, is non-Hermitian and its eigenvalues could be complex numbers. Fortunately, D is a positive operator and $\boldsymbol{D}^{1/2} \boldsymbol{X} \boldsymbol{D}^{1/2}$ is Hermitian, with the same moments as $\boldsymbol{X}\boldsymbol{D}$, which allows us to rewrite the adjacency matrix as

$$\begin{align} \boldsymbol{A} = \boldsymbol{D}^{1/2} \boldsymbol{X} \boldsymbol{D}^{1/2}. \end{align} \tag{ 48 }$$

This interesting decomposition, which follows as a natural interpretation of the main theorem in [45], can be used to study the spectral properties of random graphs with a prescribed degree distribution and large c without having to run sophisticated algorithms to sample graphs from the configuration model. This is precisely the strategy we adopt below, i.e. we compare our theoretical findings with numerical results obtained from diagonalizing equation (48).

Let us determine the spectral density of random graphs with a negative binomial degree distribution. Substituting equation (47) in equations (31) and (34), and evaluating the integrals over κ, we obtain

$$\begin{equation} \rho_{\epsilon}(\lambda) = \frac{1}{\pi} \mathrm{Im}\Bigg[\frac{z^2 + \gamma^2J^2}{z\gamma^2J^2}\Bigg], \end{equation} \tag{ 49 }$$

where the dimensionless variable $\gamma \in \mathbb{C}$, defined in terms of $\langle G \rangle$ as

$$\begin{equation} \gamma = \frac{z}{J^2 \langle G \rangle}, \end{equation} \tag{ 50 }$$

solves the transcendental equation

$$\begin{equation} \gamma^2 J^2= \frac{z^2}{\left(-\alpha \gamma e^{-\gamma} \right)^{\alpha} \Gamma (1-\alpha,-\alpha \gamma) - 1}, \end{equation} \tag{ 51 }$$

with $\Gamma(a,\xi)$ ($a \in \mathbb R$ and $\xi \in \mathbb C$) denoting the incomplete Gamma function. The solution of the fixed-point equation (51) yields the regularized spectral density (49) for any $0 \lt \alpha \lt \infty$. We recall that the strength of the degree fluctuations is controlled only by α (see equations (45) and (46)). By setting α = 1 in equation (51), we recover the equations for the spectral density of random graphs with a geometric degree distribution [35].

Figure 1 compares the regularized spectral density $\rho_{\epsilon}(\lambda)$ computed from the solutions of equation (51) with numerical results for the eigenvalues obtained from diagonalizing the adjacency matrix of equation (48). The agreement between our theoretical findings and numerical diagonalization results is excellent. In particular, we note from figure 1 that degree fluctuations modify the tails of the spectral density as well as its behaviour around λ = 0.

Zoom In Zoom Out Reset image size

We have shown in a previous work [35] that $\rho(\lambda)$ has a logarithmic divergence at λ = 0 for α = 1. In order to understand how this singular behaviour depends on α, we need to extract the functional form of $\gamma = \gamma(z)$ as $|z| \to 0$. We follow [35] and make the assumption

$$\begin{equation} \gamma(z) = \frac{\beta_1}{J} z + \frac{\beta_2(\alpha,z)}{J^2 }z^2, \end{equation} \tag{ 52 }$$

where the coefficient β₁ is independent of z and $\beta_2(\alpha,z)$ satisfies $\lim_{|z|\to 0} z^2\,\beta_2(\alpha,z) = 0$. Inserting the above ansatz in equation (51) and expanding the result up to $\mathcal{O}(z^2)$, one finds that β₁ and $\beta_2(\alpha,z)$ are given by

$$\begin{equation} \beta_1 = - i \end{equation} \tag{ 53 }$$

and

$$\begin{equation} \beta_2(\alpha,z) = -\frac{1}{2}\Bigg[i^{2\alpha}\alpha^\alpha \bigg(\frac{\beta_1}{J}z\bigg)^{\alpha-1}\Gamma(1-\alpha) + \frac{\alpha}{1-\alpha}\Bigg] \end{equation} \tag{ 54 }$$

in the regime $\alpha \in (0,1)$.

The last step is to substitute equation (52) in equation (49) and compute the limit $\epsilon \to 0^+,$ which leads to the power-law divergence

$$\begin{equation} \rho(\lambda) = \frac{1}{\pi J} \Bigg[ \alpha^{\alpha} \sin\bigg(\frac{\pi}{2}\alpha\bigg)\Gamma(1-\alpha) \left( \frac{|\lambda|}{J} \right)^{\alpha-1} - \frac{\alpha}{(1-\alpha)}\Bigg] \quad (0 \lt \alpha \lt 1) \end{equation} \tag{ 55 }$$

for $|\lambda|\to 0$. By taking the limit $\alpha \rightarrow 1$ in equation (55), we recover the logarithmic divergence obtained in [35]

$$\begin{equation} \rho(\lambda) = -\frac{1}{\pi J}\Bigg[E + \log \left( \frac{|\lambda|}{J} \right) \Bigg] \quad (\alpha=1), \end{equation} \tag{ 56 }$$

with E representing the Euler–Mascheroni constant. Figure 2 compares equation (55) with numerical solutions of equations (49) and (51) for $|\lambda| \ll 1$. The numerical results deviate from the analytic expression for values of λ below a certain threshold $|\lambda_*|$. As $\epsilon$ decreases, $|\lambda_*|$ shifts towards smaller values, confirming that the discrepancy between the numerical data and equation (55) is due to the finite values of $\epsilon$ used in the numerical solutions.

Zoom In Zoom Out Reset image size

In the homogeneous limit $\alpha \rightarrow \infty$, the variance of the rescaled degree distribution $\nu(\kappa)$ vanishes and we expect to recover the Wigner law. By using the functional relation [46]

$$\begin{equation} \Gamma (1-\alpha,-\alpha \gamma) = - \alpha \Gamma (-\alpha,-\alpha \gamma) + \left(- \alpha \gamma \right)^{-\alpha} e^{\alpha \gamma} \end{equation} \tag{ 57 }$$

and the asymptotic formula [47]

$$\begin{equation} - \alpha \left( - \alpha \gamma \right)^{\alpha} e^{-\alpha \gamma} \Gamma (-\alpha,-\alpha \gamma) = \frac{1}{\left( \gamma -1 \right)} + \frac{\gamma}{\alpha \left( \gamma -1 \right)^3} + \mathcal{O}(\alpha^{-2}) \quad (\alpha \gg 1), \end{equation} \tag{ 58 }$$

we derive from equation (51) an approximate equation for γ

$$\begin{equation} \gamma^2 J^2 = z^2 (\gamma -1) - \frac{z^2 \gamma }{\alpha (\gamma -1) }. \end{equation} \tag{ 59 }$$

In the limit $\alpha \rightarrow \infty$, the above expression reduces to a quadratic equation, whose solution yields the Wigner law (see equation (33)).

Here we do not derive the analytic expression for the spectral density $\rho(\lambda)$ that arises from solving the cubic equation (59), but we characterize the support of $\rho(\lambda)$, which plays a pivotal role for the stability of complex systems [13]. In general, when the largest eigenvalue of the adjacency matrix A is finite, there exists a regime of model parameters where the stationary states of a large complex system coupled through A are linearly stable. Thus, complex systems interacting through the symmetric random matrix of equation (48) can be in a stable state in the limit $\alpha \rightarrow \infty$, in view of the finite support of the Wigner law. An interesting question here is whether the support of $\rho(\lambda)$ remains finite when a small amount of heterogeneity is introduced ($1 \ll \alpha \lt \infty$). In order to resolve this issue, we study the discriminant $\Delta(\lambda)$ of the cubic equation (59) in the limit $\epsilon \rightarrow 0^{+}$. If $\Delta(\lambda) \gt 0$, then equation (59) has only real roots and $\rho(\lambda) = 0$, whereas if $\Delta(\lambda) \lt 0$, then equation (59) admits a pair of complex-conjugate solutions, yielding $\rho(\lambda)\gt0$. As shown in figure 3, the discriminant is zero at a certain value $\lambda = \lambda_{\mathrm{b}}$, which implies that $\rho(\lambda)$ has a finite support. The spectral edge $\lambda_{\mathrm{b}}$ of $\rho(\lambda)$ consistently approaches the value $\lambda_{\mathrm{b}} = 2 J$ of the Wigner law as α increases.

Zoom In Zoom Out Reset image size

4.2. Eigenvector localization and the distribution of the LDOSs

In this section we analyse the effect of degree fluctuations on the IPR and on the LDOSs for a negative binomial degree distribution. Substituting equation (47) in equation (35) and calculating the integral over κ, we obtain the regularized IPR around an eigenvalue λ

$$\begin{align} & \mathcal{I}_{\epsilon} (\lambda) = \frac{\epsilon}{\pi \rho_\epsilon(\lambda)}\mathrm{Im}\Bigg\{ \frac{\alpha^\alpha \gamma}{z \left[ \mathrm{Im} \left( z/\gamma \right) \right]^{\alpha-1} \left[ \epsilon+ \gamma \, \mathrm{Im} \left( z/\gamma \right) \right]} \nonumber\\ & \quad \quad\quad\ \times \Bigg[ \epsilon^{\alpha-1}\exp\left( \frac{\epsilon \, \alpha}{\mathrm{Im} \left( z/\gamma \right)}\right) \Gamma\left(1-\alpha,\frac{\epsilon \, \alpha}{\mathrm{Im} \left(z/\gamma\right)}\right) - \left[ -\gamma \mathrm{Im} \left( z/\gamma \right) \right]^{\alpha-1} e^{-\alpha\gamma} \Gamma \left( 1-\alpha,-\alpha\gamma \right) \Bigg]\Bigg\},\nonumber\\ \end{align} \tag{ 60 }$$

where γ fulfils equation (51) and the regularized spectral density $\rho_\epsilon(\lambda)$ is given by equation (49).

It is well-established that the eigenvectors of random graphs with finite c become localized in the tails of the spectral density due to the existence of hubs in the graph structure [6, 22, 30]. It is natural to ask whether such localized states survive for $c \rightarrow \infty$ in the presence of degree fluctuations. Figure 4 shows the spectral density and the IPR derived from equations (49) and (60) as a function of $\epsilon$ for large values of $|\lambda|$. In the regime $\epsilon \rightarrow 0^+$, the spectral density $\rho_{\epsilon}(\lambda)$ converges to a finite limit and the IPR vanishes as $\mathcal{I}_{\epsilon}(\lambda) \propto \epsilon$. The same picture holds for other values of α, which demonstrates that all eigenvectors corresponding to nonzero eigenvalues are extended.

Zoom In Zoom Out Reset image size

The distribution $P_z(y)$ of the LDOS probes the spatial fluctuations of the eigenvectors and it gives important information about localization phenomena. In the limit $\epsilon \rightarrow 0^{+}$, the distribution $P_z(y)$ within the localized phase typically exhibits a singularity at $y \simeq \epsilon$, due to the extensive number of sites at which $\mathrm{Im} G_{ii} \simeq \epsilon$ [14]. Differently from that, our results show that $P_z(y)$ converges to a regular, $\epsilon$-independent function in the limit $\epsilon \rightarrow 0^{+}$, highlighting the extended nature of the eigenvectors. Figure 5(a) compares equation (39) with numerical results obtained from the solutions of equations (16) and (18) using the population dynamics algorithm [32] for α = 1 and large values of c. The agreement between theoretical and numerical results is excellent over the central portion of the distribution. The discrepancy close to $y = y_{\mathrm{e}}$ in figure 5(a) is due to strong finite-connectivity effects, since the convergence of the numerical results to the asymptotic behaviour for $c \rightarrow \infty$ is extremely slow.

Zoom In Zoom Out Reset image size

Figure 5(b) shows that $P_z(y)$ diverges at the edge $y = y_{\mathrm{e}}$ for any value of α. Moreover, the distribution $P_z(y)$ develops an additional power-law singularity at y = 0 when α < 1, which is a genuine effect of strong degree fluctuations and a direct consequence of the shape of ν (see equation (47)). In the limit $\alpha \rightarrow \infty$, the graph becomes homogeneous and $P_z(y)$ converges to a Dirac-δ distribution centered at $y_{\mathrm{e}} = \pi \rho_{\mathrm{w}}(\lambda)$ ($\epsilon \rightarrow 0^{+}$), where $\rho_{\mathrm{w}}(\lambda)$ is given by equation (33).

With the aim of clarifying the singular behaviour of the spectral density (see figure 2), we turn our attention to the statistics of the LDOS at λ = 0. Substituting equation (47) in equation (43), setting $\epsilon$ = 0, and then noting that $\mathrm{Im} \langle G \rangle = 1/J$ for $z \rightarrow 0$ (see equations (50) and (52)), we obtain the simple analytic result

$$\begin{equation} P_{0}(y) = \frac{\alpha^{\alpha}}{\Gamma(\alpha) J^{\alpha}} \frac{\exp{\left(- \frac{\alpha}{J y} \right) }}{y^{\alpha+1}}, \end{equation} \tag{ 61 }$$

which reveals the unbounded character of the LDOS fluctuations at λ = 0. Figure (6)(a) confirms the exactness of expression (61) by comparing this equation with results obtained by numerically solving equations (16) and (18) for α = 1 and c = 400. It is interesting to contrast $P_{0}(y)$ with the distribution of the LDOS in the extended phase of regular random graphs with on-site random potentials [20, 48, 49]. While for regular random graphs with on-site disorder the distribution of the LDOS decays exponentially fast beyond a certain scale [48, 49], the power-law tail of equation (61) implies that the qth moment $\overline{y^q} = \int_{0}^{\infty} d y \, y^q \, P_{0}(y)$ diverges for $\alpha \leqslant q$, whereas

$$\begin{equation} \overline{y^q} = \frac{\alpha^q \Gamma \left( \alpha - q \right)}{J^q \Gamma(\alpha) } \end{equation} \tag{ 62 }$$

for $\alpha \gt q$. Figure 6(b) shows that $\overline{y}$ does not coincide with the most probable value of the distribution $P_{0}(y)$ due to its skewed shape. As α increases and the graph becomes more homogeneous, the distribution $P_{0}(y/\overline{y})$ gradually becomes more symmetric and concentrated around its mean value. For $1 \ll \alpha \lt \infty$, $P_{0}(y/\overline{y})$ is a Gaussian distribution with variance $\mathcal{O}(1/\alpha)$, and it ultimately converges to $P_{0}(y/\overline{y}) = \delta(y/\overline{y} -1)$ in the homogeneous limit $\alpha \rightarrow \infty$. We remark that the distribution of the LDOS characterizing the master operator of the sparse Barrat–Mézard trap model also exhibits a power-law tail at specific eigenvalues [23].

Zoom In Zoom Out Reset image size

The resolvent distributional equations for the spectral properties of heterogeneous random graphs do not have analytic solutions for finite mean degree c. In the limit $c \rightarrow \infty$, such equations admit a trivial solution, typical of random graphs with a homogeneous structure, in which the resolvent elements are all equal to their mean value. Here we have shown how to distill a nontrivial analytic solution of the resolvent distributional equations, valid in the high-connectivity limit, which explicitly depends on the shape of the degree distribution. This solution enables to perform a thorough analysis of the impact of degree heterogeneities on the spectral and localization properties of the adjacency matrix.

We have presented several results for the spectral and localization properties of random graphs with a negative binomial degree distribution, in which the network heterogeneity, measured by the relative variance of the degree distribution (see equation (46)), is governed by a single parameter $\alpha \in (0,\infty)$. When the degree fluctuations are sufficiently strong ($0 \lt \alpha \leqslant 1$), the spectral density $\rho(\lambda)$ diverges at the zero eigenvalue λ = 0. More specifically, the function $\rho(\lambda)$ exhibits either a logarithmic or a power-law singularity if α = 1 or $\alpha \in (0,1)$, respectively. In addition, we have shown that $\rho(\lambda)$ has a finite support in the regime of weak degree fluctuations ($1 \ll \alpha \lt \infty$), which implies that large complex systems coupled through highly connected random graphs can be found in a linearly stable state [13], at least when the variance of the degree distribution is small enough. An interesting open question is whether $\rho(\lambda)$ becomes unbounded below a critical value of α, or if the support of $\rho(\lambda)$ remains always finite.

We have shown that the IPR vanishes for nonzero eigenvalues and the corresponding eigenvectors are extended for any amount of degree fluctuations. We point out that this picture is not in conflict with recent results [18, 50] that show the existence of localized eigenvectors in the tails of the spectral density of critical random graph models. In fact, our results for the absence of localization hold for $c = \mathcal{O}(N^{a})$ (a < 1) [17], while in critical random graphs the mean degree scales as $c = \mathcal{O}(\ln N)$.

In order to further examine the nature of the eigenvectors and the singular behaviour of the spectral density, we have computed analytically the distribution of the LDOS, which quantifies the spatial fluctuations of the eigenvector amplitudes throughout the graph (see equation (38)). The distribution of the LDOS attains a nonsingular, $\epsilon$-independent limit as $\epsilon \rightarrow 0^{+}$, confirming the absence of localized eigenvectors in the high-connectivity limit [14]. The importance of degree fluctuations is more evident at the zero eigenvalue λ = 0, where the distribution of the LDOS exhibits a power-law tail with exponent α + 1 (see equation (61)). In particular, the divergence of the mean value of the LDOS at λ = 0 explains the singular behaviour of $\rho(\lambda)$ for $\alpha \leqslant 1$.

It is interesting to compare our analytic expression for the distribution of the LDOS at λ = 0 with the analogous result for the extended phase of sparse regular random graphs with on-site disorder [48, 49]. In the latter class of models, the distribution of the LDOS decays exponentially fast and all its moments are finite, whereas in the present model the qth moment diverges for $\alpha \leqslant q$. In particular, the second moment of the LDOS for regular random graphs only diverges as the critical point for the Anderson transition is approached from the delocalized phase by increasing the strength of the diagonal disorder [43, 51, 52]. In an analogous way, the second moment of the LDOS in the present model is finite for α > 2 and it diverges for $\alpha \leqslant 2$, which seems to suggest that highly-connected random graphs with strongly fluctuating degrees lie in a critical regime [40, 51]. Besides constituting an interesting benchmark to study how degree heterogeneities affect the spectral properties of networks, our analytic findings open the possibility to investigate how the interplay between on-site disorder and fluctuations in the network topology modify the Anderson localization transition.

Overall, our results uncover an interesting high-connectivity regime in which the resolvent equations admit exact and nontrivial solutions that incorporate heterogeneous features of the network topology. Thus, it would be interesting to generalize the techniques developed in this work to solve the resolvent equations for the adjacency matrix of directed random graphs [41, 53] and networks with loops [31], as well as the analogous equations for the Laplacian matrix on graphs [54]. We also expect that the resolvent equations of heavy-tailed random matrices [55, 56] can be studied with our approach, since such ensembles can be mapped on a sparse random-matrix ensemble through the introduction of a cutoff that distinguishes between strong and weak matrix elements [22]. Finally, it would be interesting to solve the resolvent equations for the master operator of the sparse Barrat–Mézard model [23] in the high-connectivity limit, which could lead to analytic results for the low-temperature localization properties. Work along these lines is under way.

J D S acknowledges a fellowship from CNPq/Brazil. F L M thanks London Mathematical Laboratory and CNPq/Brazil for financial support.

All data that support the findings of this study are included within the article (and any supplementary files).

In this appendix we present a more formal derivation of equations (28) for the probability density $\mathcal{W}_z(s)$, from which all subsequent results for the spectral and localization properties follow. By inspecting equations (16) and (19), we note that $\mathcal{W}_z(s)$ can be written as

$$\begin{equation} \mathcal{W}_z(s) = \sum_{k=0}^\infty p_k \int_{\mathbb{H}^+} \left[\prod_{\ell=1}^{k}d g_\ell \mathcal{Q}_{z}(g_\ell) \right] \int_{\mathbb{R}} \left[ \prod_{\ell=1}^{k} d {J_{\ell}} p_J(J_\ell) \right] \delta\Bigg ({s - \sum_{\ell=1}^{k}J_{\ell}^2 g_\ell}\Bigg), \end{equation} \tag{ A1 }$$

with $d g = d \mathrm{Re} g \, d \mathrm{Im} g$. In a similar fashion, one can introduce the distribution $\mathrm{W}_{z}(s)$ associated to $\mathcal{Q}_{z}(g)$. The average of an arbitrary function f(G) of the cavity resolvent G distributed according to $\mathcal{Q}_{z}(g)$,

$$\begin{equation} \langle f(G) \rangle = \int_{{{\mathbb{H}}^+}} dg \mathcal{Q}_{z}(g) f(g), \end{equation} \tag{ A2 }$$

is recast in the form

$$\begin{equation} \langle f(G) \rangle = \int_{{{\mathbb{H}}^+}} ds \mathrm{W}_{z}(s) f \left( \frac{1}{z-s} \right) \quad d s = d \mathrm{Re} s \, d \mathrm{Im} s, \end{equation} \tag{ A3 }$$

where the expression for $\mathrm{W}_{z}(s)$ is inferred from equation (18)

$$\begin{equation} \mathrm{W}_{z}(s) = \sum_{k=1}^\infty \frac{k}{c}p_k \int_{{{\mathbb{H}}^+}} \left[\prod_{\ell=1}^{k-1}d g_\ell \mathcal{Q}_{z}(g_\ell) \right] \int_{\mathbb{R}} \left[ \prod_{\ell=1}^{k-1} d {J_{\ell}} p_J(J_\ell) \right] \delta\Bigg ({s - \sum_{\ell=1}^{k-1}J_{\ell}^2 g_\ell}\Bigg). \end{equation} \tag{ A4 }$$

The quantity $\mathrm{W}_{z}(s)$ is the probability density of the random variable defined in equation (19) with the replacement $k \rightarrow k-1$. In particular, it follows from equation (A2) that the average resolvent $\langle G \rangle$ on the cavity graph is given by

$$\begin{align} \langle G \rangle = \int_{{\mathbb{H}}^+} ds \frac{\mathrm{W}_{z}(s)}{z-s}. \end{align} \tag{ A5 }$$

The distributions $\mathcal{W}_{z}(s)$ and $\mathrm{W}_{z}(s)$ fully determine the spectral properties of the adjacency matrix.

Our aim is to calculate the joint distributions $\mathcal{W}_{z}(s)$ and $\mathrm{W}_{z}(s)$ for $c \rightarrow \infty$. Given that $\mathcal{W}_{z}(s)$ and $\mathrm{W}_{z}(s)$ are distributions of sums of independent and identically distributed random variables, it is natural to work with the characteristic functions of such distributions. Let $\mathcal{V}(u,v)$ and $\mathrm{V}(u,v)$ be the characteristic functions of, respectively, $\mathcal{W}_{z}(s)$ and $\mathrm{W}_{z}(s)$, defined as

$$\begin{align} \mathcal{V}(u,v) &= \int_{{\mathbb{H}}^+} ds \mathcal{W}_{z}(s) \exp{\left(- i u \mathrm{Re} s - i v \mathrm{Im} s \right)}, \end{align} \tag{ A6 }$$

$$\begin{align} \mathrm{V}(u,v) &= \int_{{\mathbb{H}}^+} ds \mathrm{W}_{z}(s) \exp{\left(- i u \mathrm{Re} s - i v \mathrm{Im} s \right)}. \end{align} \tag{ A7 }$$

Inserting equations (A1) and (A4) in the above expressions, we obtain

$$\begin{align} &\mathcal{V}(u,v) = \sum_{k=0}^\infty p_k \exp \big[k\mathcal{S}_c(u,v)\big], \end{align} \tag{ A8 }$$

$$\begin{align} &\mathrm{V}(u,v) = \sum_{k=1}^\infty \frac{k}{c}p_k \exp \big[(k-1)\mathcal{S}_c(u,v)\big], \end{align} \tag{ A9 }$$

with

$$\begin{equation} \mathcal{S}_c(u,v) = \ln \Bigg[\int_{{\mathbb{H}}^+} d g \mathcal{Q}_{z}(g) \int_{-\infty}^{+\infty}d J p_J(J) \exp\Big(- i u J^2 \mathrm{Re} g - i v J^2 \mathrm{Im} g \Big) \Bigg]. \end{equation} \tag{ A10 }$$

Since the second moment of the coupling strengths is of $\mathcal{O}(1/c)$, the leading term of the above equation for $c \gg 1$ is given by

$$\begin{equation} \mathcal{S}_c(u,v) = - i u \frac{J^2}{c} \mathrm{Re} \langle G \rangle - i v \frac{J^2}{c} \mathrm{Im} \langle G \rangle , \end{equation} \tag{ A11 }$$

where we assumed that $\langle G \rangle$ attains a well-defined limit for $c \to \infty$. The substitution of the above expression for $\mathcal{S}_c(u,v)$ in equations (A8) and (A9) leads to the following equations for $c \rightarrow \infty$

$$\begin{align} \mathcal{V}(u,v) = \int_{0}^{\infty} d \kappa \, \nu(\kappa) \exp{\left( - i u \kappa J^2 \mathrm{Re} \langle G \rangle - i v \kappa J^2 \mathrm{Im} \langle G \rangle \right)}, \end{align} \tag{ A12 }$$

$$\begin{align} \mathrm{V}(u,v) = \int_{0}^{\infty} d \kappa \, \kappa \, \nu(\kappa) \exp{\left( - i u \kappa J^2 \mathrm{Re} \langle G \rangle - i v \kappa J^2 \mathrm{Im} \langle G \rangle \right)}, \end{align} \tag{ A13 }$$

where the probability distribution $\nu(\kappa)$ of the rescaled degrees is defined in equation (27). Performing the inverse Fourier transform of $\mathcal{V}(u,v)$ and $\mathrm{V}(u,v)$, we get

$$\begin{align} &\mathcal{W}_{z}(s) = \int_{0}^{\infty} d\kappa\,\nu(\kappa) \delta(s - \kappa J^2 \langle G\rangle), \end{align} \tag{ A14 }$$

$$\begin{align} &\mathrm{W}_{z}(s) = \int_{0}^{\infty} d\kappa\,\kappa\,\nu(\kappa) \delta(s - \kappa J^2 \langle G\rangle). \end{align} \tag{ A15 }$$

Equation (A14) means that the complex random variable S, distributed according to $\mathcal{W}_{z}(s)$, is equal in distribution to the random variable $\kappa J^2 \langle G\rangle$. Thus, given $\nu(\kappa)$, equation (28) follows by making a change of variables. The self-consistent equation for $\langle G \rangle$, equation (31), is readily obtained by inserting equation (A15) in equation (A5). This completes the calculation of $\mathcal{W}_{z}(s)$.