Multifractality and statistical localization in highly heterogeneous random networks

We consider highly heterogeneous random networks with symmetric interactions in the limit of high connectivity. A key feature of this system is that the spectral density of the corresponding ensemble exhibits a divergence within the bulk. We study the structure of the eigenvectors associated with this divergence and find that they are multifractal with the statistics of eigenvector elements matching those of the resolvent entries. The corresponding localization mechanism relies on the statistical properties of the nodes rather than on any spatial structure around a localization centre. This “statistical localization” mechanism is potentially relevant for explaining localization in different models that display singularities in the bulk of the spectrum of eigenvalues.


I. INTRODUCTION
After three decades of an extensive analysis of its main properties, the configuration model has stood its ground as the simplest but realistic theoretical model for complex networks [1,2].A particularly appealing property of this model is that the degree distribution can be freely specified, making it possible to capture a large range of random networks varying from highly homogeneous (all nodes with same degree) to very heterogeneous [3].
In addition to the topological structure of the network, the edges are often weighted, physically representing the strength of an interaction between nodes [4,5].This is encoded in a matrix that contains the weights of all edges in the network.For the theoretical approach, rather than focusing on individual matrices, one reduces them to their defining properties and represents these in a statistical ensemble of weight matrices [6].The interplay between interaction and structure can be studied by analyzing, for instance, the spectral properties of the corresponding ensemble of random matrices (see for instance ref. [7] for a classical work in the case of symmetric sparse matrices or ref. [8] for a more recent perspective on non-Hermitian random matrices).
Analytical methods for computing spectra of configuration model networks can be found in references [9][10][11][12].Most of the theory is based on message-passing algorithms and the corresponding cavity (or resolvent) equations [12,13].Recently, in references [3,14] an analysis of these equations for highly heterogeneous networks, i.e. those with non-vanishing relative variance of the degree distribution (in the high connectivity limit), has been performed.The main result is that in the thermodynamic limit, classical results of Random Matrix Theory do not apply.This includes in particular Wigner's semicircle law to describe the spectrum of eigenvalues.Deviations from this are accompanied also by other striking properties, such as the appearance of singularities in the bulk of the spectral density.
In this paper, we study highly heterogeneous random networks using the negative binomial distribution as the degree distribution for the configuration model, with interactions sampled from a Gaussian with zero mean and a variance that scales inversely with the mean connectivity.As has been shown in reference [14], this setting leads to a family of networks controlled by a single parameter that allows one to cover the range from highly heterogeneous all the way to essentially homogeneous networks.
In reference [14], the resolvent equations were derived and solved for this system.Their solution revealed that in the case of highly heterogeneous graphs, there is a divergence in the spectral density at eigenvalue zero.This can be traced back to the distribution of the local Density of States, which acquires a power-law tail leading to an infinite first moment.
Here we investigate the structure of the eigenvectors associated with the divergence of the spectral density.Our main finding is that they are localized and exhibit "strong multifractal" behavior.We rationalize this by finding an anticorrelation between the degree of the node and the corresponding amplitude of the eigenvector.It turns out that localization occurs in nodes with low degree (relative to the mean), which typically are far apart from each other.We denote the corresponding localization mechanism "statistical localization" as it is driven by the statistical properties of the node(s).This contrasts with the standard case of Anderson localization on random networks [15][16][17], where eigenvectors are localized on nodes around a single localization centre.We conjecture that the new mechanism of statistical localization is relevant also for other types of models that exhibit singularities within the bulk of the spectral density such as the Poisson random graph with Gaussian couplings [7] or the sparse Barrat-Mézard trap model [18,19].

II. HETEROGENEOUS WEIGHTED RANDOM NETWORKS
We provide in this section the essential pieces of information that define the ensemble of random networks that we study.For more details, we refer to the reader arXiv:2311.15808v1[cond-mat.dis-nn]27 Nov 2023 to the original paper [14].Let us consider a simple and undirected network with N nodes.We consider the configuration model with a negative binomial distribution of degrees k, where c is the mean connectivity, Γ(•) is the Gamma function and 0 < γ < ∞ controls the heterogeneity of the distribution1 .Indeed, the relative variance of the degree distribution is given by Each node i is assigned a degree k i drawn randomly from p k , and nodes are then randomly connected following the standard configuration model prescription [2].If there is a link between nodes i and j we set the random variable c ij = 1, otherwise it is set to zero.Additionally we consider interactions J ij between nodes as random variables sampled independently on each link (with J ij = J ji ) from a distribution p J with mean zero and standard deviation J/ √ c.Overall, the weight matrix A has elements In the high-connectivity limit c → ∞, it makes sense to consider the distribution of the rescaled degrees κ = k/c.This distribution ν(κ) is formally defined as follows For p k as given in eq. ( 1), one finds that ν(κ) is a Gamma distribution with shape parameter γ and scale parameter 1/γ, i.e.
The regime γ < 1 characterizes random networks with strongly heterogeneous degrees as the relative variance of the degree distribution is greater than 1 (cf.eq. ( 2)).In the limit , of the weight matrix (eq.( 3)) exhibits a power law divergence with 0 < γ < 1.This behavior can be rationalized via the connection between the resolvent and the spectral density [20], namely where the G ii are the diagonal elements of the resolvent matrix G = ((λ−iϵ)I −A) −1 , with I the identity matrix.As a matter of fact, networks generated by the configuration model exhibit locally a tree-like structure [2], this implies that we can use the cavity method [10,12] to estimate the diagonal elements of the resolvent matrix, also called local Density of States (lDOS).Ref. [14] shows that at λ = 0 the distribution of y i = Im G ii in the limit ϵ → 0 is given by By virtue of the resolvent identity (eq.( 7)), the mean value of the lDOS distribution determines the spectral density.From the power-law tail of the distribution (8) one can then clearly see the origin of the divergence in equation (6).
A final result that we want to introduce here comes from reference [21].This states that in the limit c → ∞, the spectral density of the weight matrix A is given by the free multiplicative convolution [22,23] of ν(κ) with the Wigner semicircle law ρ W (λ).As an implication, one has that the weight matrix can be decomposed as the product of two different matrices, X and D, where the latter is the degree matrix with elements D ij = κ i δ ij and the former is a random (symmetric) matrix with zeros on the diagonal and off-diagonal elements drawn from a Gaussian distribution with mean zero and variance J2 /N .After symmetrization, one can write This decomposition is used in reference [14] to compare the predictions (in particular, for the spectral density) from theory to numerical results using exact diagonalization, without explicitly constructing the configuration model.We will use the same approach to obtain the eigenvectors associated with the modes around λ = 0 numerically.

III. WAVEFUNCTION STATISTICS
We exploit the decomposition in equation ( 9) to generate an interaction matrix A. Then we diagonalize it with the Lanczos algorithm to extract the eigenvectors with associated eigenvalues closest to zero 2 .In Figure 1((a suggests that the distributions have power-law tails for large x, with γ-dependent exponents 3 .In fact, we found that the exponents agree with those for the resolvent, i.e.P (x) ∼ P 0 (x) ∼ x −(1+γ) for sufficiently large x, as shown in Figure 1((b)).This agreement of the statistics of the resolvent and eigenvector entries is consistent with previous results for eigenvectors with multifractal properties, see for instance [25].
The power-law tail motivates the following ansatz for the whole distribution above a scale x min (N ), close to 3 The figure also suggests the existence of a left power-law tail.
This would contribute to the scaling of the exponents τ (q) for negative q (see eq. ( 16)) and equivalently to the multifractal spectrum for γα > 1 (cf.eq. ( 15)).As this piece is not relevant for the understanding of the localization properties, we do not consider it further in our analysis.
the mode of the distribution with b(N ) and x min (N ) functions to be determined.Their scaling can be estimated by using the normalization of P (x) and its first moment [25]; the latter follows from the normalization of the eigenvectors, 1 = These two relations imply b ∼ N γ−1 and x min ∼ N 1−1/γ .As x min sets the scale from which the power-law becomes valid, we identify this as the typical value of the distribution.This scaling is confirmed using results from Exact Diagonalization as can be appreciated in Figure 2. Substituting the scaling for b(N ) in our ansatz (eq.( 10)), we write the distribution P (x) in the generic way with A = O(N 0 ) a normalization constant.The spectrum of fractal dimensions f (α), which is defined as [26,27] where on the r.h.s.x = N 1−α or conversely α = 1 − ln x/ ln N .Substitution of eq. ( 13) into eq.( 14) gives the expression FIG. 3: (a) Scaling of the wavefunction moments as encoded in the exponents τ (q) for N = 2 14 .The dashed lines correspond to the linear piece of the theoretical prediction (17).See numerical details in the appendix.
where the upper cutoff α max = 1 γ corresponds to the lower cutoff x min .Finally, the Legendre transform of equation ( 14) gives the set of exponents characterizing the q-th moment of the distribution Substitution of f (α) (eq.( 15)) into equation ( 16) yields Figure 3 compares this result with Exact Diagonalization.

IV. DISCUSSION
The scaling of the moments I q as given by equation (17) has been found across different models and the corresponding phase has been described variously as "quasilocalized" [27], or localized with multifractal properties [28], or localized with "strong multifractal" behavior [17].The corresponding models are the Gaussian Rosenzweig-Porter, random Levy matrices, and the Anderson model on small-world networks, respectively.References [27,28] point to a power-law (instead of exponential) localization as the origin of this kind of behavior.In contrast, reference [17] (see also [29]) focusses on the existence of a length scale that characterizes the decay of wave functions along the typical branches as responsible for the scaling above.These mechanisms, however, do not seem to be operative for our system.
In order to find the mechanism that governs the behavior of our model, we construct finite instances of the networks generated with the configuration model.Then we investigate the correlation between the squared eigenvector entry (or "mass") and the degree of the node.Our finding is that low-degree nodes, typically leaves (i.e.nodes with degree one), concentrate the mass of the eigenvectors.For mildly heterogenous networks, i.e. those with 1 < γ < 2, the anticorrelation between degree and eigenvector mass is easily visible in a scatter plot (c.f. Figure 4).If we extrapolate this finding to c → ∞, we expect that for a given instance nodes with small relative degree will concentrate the mass of the eigenvector, and what is more relevant, that those nodes may be far apart from each other on the network.Indeed, we find that the dominant nodes from figure 4 lie at distances of the order of the graph diameter from each other.Thus what matters for localization is the identity of each node and its statistical properties, rather than its spatial location in the network.In this sense, we say that the system exhibits "statistical localization" and that this leads to the multifractal behavior encoded in the exponents (17).
To illustrate further the mechanism described above, consider the Bouchad trap model on a sparse random graph [30].The ground state of this system is the Boltzmann distribution, and the set of exponents characterizing the scaling of the moments is of the form (17), with temperature instead of γ as the control parameter (see equation 4.15 in reference [31]).Thus, for finite instances, the nodes that would carry most of the mass of the vector are the ones with the largest energies, corresponding to the deepest traps.Those nodes are spatially uncorrelated and are instead identified by local statistical properties.
At this point, it is worth comparing our results with earlier observations of localization in heterogenous networks.As a representative example of this class we consider the Laplacian on Poisson (Erdös-Rényi) random graphs [7,32] (see also [33]).In both cases, for small |λ|, localization is driven by low degree nodes.However, for our model, this happens in the bulk of the spectrum and is accompanied by a divergence in the spectral density; whereas for the Poisson random graph, localization happens in the (Lifshitz [34]) tail of the spectrum, which is separated from the bulk by a mobility edge.Additionally, while it has been observed that localized states in this tail are centred on geometric defects with abnormally low connectity [32] and that the density of states is dominated by low degree nodes [7], it remains an open question whether those states exhibit multifractal behavior and the existence of multiple localization centres.This would be interesting to study in future work.
In conclusion, we have analysed the wavefunction statistics for highly heterogeneous random graphs from the configuration model defined by a negative binomial degree distribution and in the limit of high connectivity.We have found that the distribution of the eigenvector mass for the modes that contribute to the divergence of the spectral density is a power-law, with the same exponent as the one for the local Density of States.Those eigenvectors exhibit strong multifractality, and for finite graphs are essentially localized in the lowest-degree nodes.On this basis, we introduce the concept of "statistical localization" that is to be contrasted with the standard one of "spatial localization".We leave for future work an extensive analysis of the eigenstate correlations of these modes as has been done for the Anderson model [35] or the Gaussian Rosenzweig-Porter model [27], in order to investigate specific footprints of the statistical localization mechanism.It will also be interesting to carry out an analysis similar to the one in this paper for the sparse Barrat-Mézard trap model [18,19], which like the model considered here exhibits divergences in the bulk of the spectral density.Work in this direction is in progress.Finally, we point out that the phase described in earlier studies as a "frozen phase" has similar phenomenology to the one observed here [36][37][38], and thus may potentially also be viewed as an instance of the notion of statistical localization that we have introduced here.A quantitative analysis along these lines is left as an interesting task for future work.
FIG. 1: (a) Distribution of the (logarithm) of squared wavefunction amplitudes, x i = N |ψ i | 2 , for three different eigenvectors of interaction matrices of size N = 2 14 generated for different heterogeneity parameters γ.(b) Estimation of the power law tail of the distribution P (x) ∼ x −g(γ) , fitted from data for x > x * with x * = 20 x typ well above the mode (∼ x typ ) of the distribution.Error bars show statistics over 128 instances of A, using the 50 eigenvectors with eigenvalues closest to λ = 0 for each instance.

FIG. 4 :
FIG. 4: Correlation among the scaled eigenvector mass, x i = N |ψ i | 2 and the node degree k i for the eigenvector with eigenvalue closest to 0 for an instance constructed via equation (3), for a finite configuration model network with mean connectivity c = 400 and γ = 1.2 and of size N = 8187.(This network is the giant connected component of the original configuration model graph with 2 13 = 8192 nodes.)