Anatomy of the eigenstates distribution: a quest for a genuine multifractality

Motivated by a series of recent works, an interest in multifractal phases has risen as they are believed to be present in the Many-Body Localized (MBL) phase and are of high demand in quantum annealing and machine learning. Inspired by the success of the RosenzweigPorter (RP) model with Gaussian-distributed hopping elements, several RP-like ensembles with the fat-tailed distributed hopping terms have been proposed, with claims that they host the desired multifractal phase. In the present work, we develop a general (graphical) approach allowing a self-consistent analytical calculation of fractal dimensions for a generic RP model and investigate what features of the RP Hamiltonians can be responsible for the multifractal phase emergence. We conclude that the only feature contributing to a genuine multifractality is the on-site energies' distribution, meaning that no random matrix model with a statistically homogeneous distribution of diagonal disorder and uncorrelated off-diagonal terms can host a multifractal phase.


Introduction
A study of eigenstates of random Hamiltonians is crucial for understanding the spectral and transport properties of various systems.For example, basis-invariant ensembles like the Gaussian Orthogonal Ensemble (GOE) or Gaussian Unitary Ensemble (GUE) [1] have basisinvariant distributions of eigenstates, meaning that, on average, all their components are of the same order.This leads to the ergodicity of eigenstates and, thus, makes the system conducting.On the other hand, the 1d Anderson model's eigenstates are exponentially localized around their maxima, meaning that any local perturbation to the system only affects a finite number of such eigenstates, and the system is an insulator.
The examples of ergodic and localized states are just the opposite ends of the broad spectrum of what is available.For example, systems at the critical points of the Anderson transition between the ergodic and localized phases are known to host the so-called fractal, or even multifractal, eigenstates [2].In a thermodynamic limit of the system size going to infinity, such critical eigenstates occupy an infinite number of sites, being still a measure zero of the total system size.The difference between fractal and multifractal eigenstates is even more subtle than that between fractal and localized ones.This is the difference in the probability density functions (PDF) of the eigenstate coefficients in these two cases.Fractal states are charac-terized by a single power-law tail in PDF which usually provides the only energy/time scale in the system in addition to the standard ones: the bandwidth and the mean level spacing.Here, by an additional energy/time scale, we mean the Thouless energy/time, like in [3,4], which appears to be the only parameter, needed to collapse the power-law distributions of wave-function coefficients.Multifractal states, in turn, have multiple (running) power-law exponents at different scales.This means that the above distribution of eigenstate coefficients cannot be approximated by a single power-law, but needs many power-law-like pieces at different wave-function amplitudes.The simplest example of such a multifractal distribution would be a log-normal one, known to be the signature of weak multifractality [2].In this respect, the multifractal distribution, being the most general type of smooth distribution, provides the set of energy scales and, thus, such systems usually have more rich time evolution and hierarchical local spectral structure, needed for the quantum-algorithm speed-up applications [5] and faster learning of artificial neural networks [6].For these applications, the robust multifractal phases of matter are of crucial importance.
As another motivation, we consider a series of recent works [7][8][9][10], providing numerical evidence of the Hilbert space multifractality in the entire phase of many-body localization (MBL) in disordered interacting quantum systems.Since the direct study of MBL is challenging numerically and analytically, it is of particular interest and high demand to develop some proxy models that mimic the necessary properties of the MBL systems but are easier to tackle.One of such proxies is an Anderson model on random regular graphs (RRG) [11,12], possessing a hierarchical structure, similar to the many-body Hilbert space.It was believed to host multifractal states [13][14][15][16][17][18][19][20][21], even before similar claims about the actual many-body systems have been made [12].To simplify the problem even more, a random-matrix proxy to the RRG was proposed, namely the log-normal Rosenzweig-Porter model (LN-RP) [18,22].Due to the multifractal and heavy-tailed nature of the log-normal distribution, it was believed to host a genuine multifractal phase.However, the suggested mapping implied that the RRG corresponds to the such parameters of the LN-RP models, where the multifractal phase disappears at the tricritical Anderson transition point, meaning that the observed RRG multifractality can be just a finite-size effect, see also [23][24][25][26][27][28].
Nevertheless, the very existence of the multifractal phase in the log-normal Rosenzweig-Porter model (unlike the fractal one in the Gaussian RP [29]) has never been proven mathematically and was based on the numerical evidence and the heuristic argument predicting multifractality for any RP model with heavy-tailed distribution, e.g., a Lévy-RP model [30,31].Given that the Gaussian Rosenzweig-Porter model [32][33][34][35][36][37], eight years after the discovery of the fractal phase there [38], is still the only analytically tractable model [29,30,39,40] hosting an entire phase of fractal states, it makes sense to look for an analytical approach to the other RP models, like Lévy- [30,31] or LN-RP [3,18,22], as such attempt does not look hopeless.In this paper, we show that it is indeed possible.
The paper is organized as follows.In Sec. 2, we define our main object of study, the spectrum of fractal dimensions, and discuss how it allows us to distinguish the multifractality from the fractality.Section 3 introduces a graphical approach to deal with the spectra of fractal dimensions (SFDs) of different independent random variables.Methodologically, this section contains the paper's main results.In Sec. 4, we demonstrate the application of the developed machinery to the Gaussian Rosenzweig-Porter model and compare our result to the previously known ones.Section 5 shows a predictive power of the above-developed method by applying it to the Lévy-RP model and demonstrates that its local density of states (LDOS) has a fractal, but not multifractal distribution.In Sec. 6, we do the same for the log-normal RP model and, again, claim the absence of LDOS multifractality.In Sec. 7, we analyze if the lack of LDOS multifractality implies the absence of eigenstates multifractality and conclude that, for models with RP-like LDOS SFDs, it does.Finally, in Sec. 8, we prove that no RP-like model, i.e., a model with a regular on-site disorder, no correlations between hopping elements, and no spatial structure, can host a multifractal phase.

Multifractality and the spectrum of fractal dimensions
For a non-negative random variable X (e.g., the eigenstate amplitude |ψ(i)|2 at site i), the spectrum of fractal dimensions (SFD) f X (α; N ) parameterizes the probability density function p X (x) as p X (x)dx = p X (N −α ) ln(N )N −α dα = ln(N )N f X (α;N )−1 dα . ( For the wave-function amplitudes, N is the system size, which is considered to be large.In other words, we focus mainly on the large-N limit f X (α) = lim N →∞ f X (α; N ) if it is not stated explicitly otherwise.The intuition behind the spectrum of fractal dimensions is as follows: in the set of N independent samples of X , we will find around N f X (α;N ) samples of the order of X ∼ N −α .Note that the parameter N can, in principle, be any number, but, in our physical applications, we will associate it with the system size.
The SFD contains all the necessary information to extract the eigenstate fractal dimensions D q and the critical exponents τ q = (q − 1)D q [2].In the large-N limit, τ q is given by the Legendre transform of f (α) in the saddle-point approximation: (2) Here the first equality is the definition of τ q via the logarithm log N of base N of the inverse participation ratio IPR q = i |ψ(i)| 2q , with the sum approximated by the averaging over the probability distribution (1), noted by 〈. ..〉.In the r.h.s. of (2), the Legendre transform is given by the parameter α q , minimizing qα − f (α).For a smooth f (α), α q is defined as the solution of the equation f ′ (α q ) = q. 1 A pictorial representation of this minimization is shown in Fig. 1.
From this result, one can readily identify at least two important values of α: α 0 and α 1 .The value α 0 corresponds to the maximum of the SFD, f (α 0 ) = max α { f (α)} = 1. 2 Using the intuitive meaning of SFD mentioned above, one can deduce that, in the sufficiently large sample, the realizations with x ∼ N −α 0 will prevail over any other values of x, making this value typical for the ensemble.On the other hand, the mean value of x corresponds to α 1 .For example, the SFDs of normalized wave functions with 1 = N i=1 |ψ(i)| 2 = N |ψ(i)| 2 always lie below the line f (α) = α and have at least one common point with this line in the thermodynamic limit.Finally, this interpretation allows us to give a mathematically strict definition of the support set dimension D [41].To see this, let us calculate D 1 as a limit of D q for q → 1 via the L'Hôpital's rule: (3) Thus, since α 1 is responsible for the normalization, D 1 = α 1 represents the actual, coherent with its physical meaning, support set dimension.
As it has been just shown, in a generic case, the different fractal dimensions are given by the different points of the SFD.However, the SFD can also have discontinuities in its first derivative.In this case, each α, corresponding to one of such discontinuities, contributes to D q Figure 1: Graphical representation of the relation between the spectrum of fractal dimensions f (α) and the critical exponents τ q , Eq. ( 2).Here, the solid red line shows the SFD of a certain distribution, while the tangent dashed lines of different colors with the slopes 1 (green), 2/3 (orange), and 1/3 (blue) correspond to Legendre transform from the r.h.s of (2).The label "D 1 " corresponds to the value of α = α q→1 , responsible for the normalization N N −α ∝ N 0 and, hence, represents the actual, correctly defined, support set dimension [41], see the main text for more details.
in an entire range of different q's, meaning that the same fraction of the eigenstate components contributes to different fractal dimensions and D q stays constant in the corresponding range of q's.The eigenstates are called multifractal [2] if D q is not constant for integer positive q or, in other words, it f (α) does have finite values at α < α 1 . 3Otherwise, we will talk about fractality.
Note that the concept considered here is very similar in spirit to the large-deviation theory (see, e.g., the review [42]), where all the physical objects, like the partition or generating functions, probability distributions of extensive variables are considered in the same way.The main difference is in the large scaling parameter: in the large-deviation theory it is usually time, while in our case it is given by the logarithm of the matrix size ln N .For more examples of such analogies, one can also look into [43] and references therein.

Graphical algebra
Working with probabilities, we often need to calculate probability distributions of composite random variables, i.e., a PDF of a sum or a product.Such quantities can be obtained via proper integral convolutions or using characteristic functions.However, the calculation difficulty grows rapidly with the complexity of the composite random variable's expression.Working with fractal spectra brings a whole new life to such operations due to the possibility of approximate integration using the Laplace method.In this section, we derive simple pictorial rules allowing the calculation of composite random variables' SFDs on the fly.
Below, we will consider the following operations: an exponentiation of a random variable (Sec.3.1), a sum of two independent random variables (i.r.v.s) (Sec.3.2), an "ensemble mix-

𝑛
Figure 2: Pictorial representation of random variable exponentiation.In this case, n is assumed to be greater than one.Here and below, we omit labeling of the axis, always assuming the horizontal axis to represent α and the vertical axis to represent f (α).The dashed horizontal line marks the level f (α) = 1, the solid horizontal line corresponds to f (α) = 0, and the labels "α 0 " and "nα 0 " correspond to the typical values of α before and after the exponentiation.
ture" of several independent random variables with different distributions (Sec.3.3), and a product of two i.r.v.s (Sec.3.4).For the exponentiation and the ensemble mixture, the derivations of the graphical rules are simple regardless of the random-variable (r.v.) distributions.
The derivations for the rest two operations, become more straightforward under the assumption that their SFDs are convex functions.However, this does not limit the applicability of the results since sums and products of any two i.r.v.s can always be decomposed into a superposition of the operations involving i.r.v.s with convex SFDs only.

Raising a random variable to a power
We start by considering the simplest case, namely, by expressing a spectrum of fractal dimensions f X n (α) of X n in terms of the spectrum of fractal dimensions f X (α) of X .Since, by definition p X (N −α ) ∝ N f X (α)+α−1 and, by variables substitute, p X n (x) = n −1 x 1/n−1 p X (x 1/n ), we get that and, hence, Pictorially, this is just an n-times stretching of the spectrum of fractal dimensions in a horizontal direction with a fixed point α = 0, see Fig. 2. For the negative n, this is also a reflection with respect to the vertical axis.

Sum of two independent random variables
The main point for the sum rule, written in terms of SFDs, is that the expression In other words, as soon as x < y, we can always neglect x completely, and vice versa.
Without loss of generality, let us assume that α 0 (X ) < α 0 (Y ), i.e. that typically Y is negligible.Then there are three cases: α < α 0 (X ), α 0 (X ) < α < α 0 (Y ), and α > α 0 (Y ).let us consider them one by one.The case α < α 0 (X ) is simple: since N −α is bigger than the typical values of both X and Y , here, both P(α X < α Y |α X = α) and P(α Y < α X |α Y = α) are of the order of unity as the typical values α 0 (X ), α 0 (Y ) are in the interval α X ,Y > α.The case α > α 0 (Y ) > α 0 (X ) is not much different: since now the typical values of X and Y are both bigger than N −α , here, the corresponding conditional probabilities are denoted 4 by N f X ,Y (α)−1 .And, as both f X (α) and f Y (α) in this region are smaller than one, this result shows a suppression of zeros which is quite logical: the more positive terms we add, the fewer probable we get something small.Finally, in the intermediate case, α 0 (X ) < α < α 0 (Y ), the probability P(α X < α Y |α X = α) is of the order of one, while the probability In the last equality we used the fact that f Y (α)−1 < 0 at α < α 0 (Y ).A pictorial representation of this result is shown in Fig. 3.

Weighted ensemble mixture of several independent random variables
Suppose that we are interested in the SFD of a random variable X with a probability density function proportional to a weighted sum, j w j = 1, of other probability density functions: Such a definition corresponds to the notion of conditional probability and the chain rule.Indeed, going from the sum to an integral and making substitutions w i → p(i)d i and p X i (x) → p(x|i), we arrive at the probably more familiar expression Now, the transformation of the PDF to the SFD in the saddle-point approximation gives the relation where f (i) is the SFD corresponding to the probability density function p(i).In the previous section, we already saw a similarly looking relation (7) for two discrete i = X , Y values with the same weights.It was a particular case of this more general mix rule. 5his is the mix rule that allows us to consider explicitly only r.v.s with convex SFDs without losing generality, as it was mentioned earlier in the introduction to this chapter.Indeed, it is easy to notice that any SFD can be defined as a mix of convex SFDs, while both sum and product of two i.r.v.s, being bi-linear operations, can be represented as the mixes of sums and products of the convex SFDs.
: Graphical representation of the weighted mix rule.In this particular case, the resulting distribution consists of the 'blue' distribution with the weight N 1 and the 'orange' distribution with the weight N ω .Notice how the orange distribution in the r.h.s was shifted vertically according to its weight.After making such shifts for all involved SFDs, the resulting SFD is obtained by a simple envelope.
Figure 5: Graphical representation of a product of the two i.r.v.s in terms of the weighted mix.In this figure, the 'blue' SFD plays the role of weights, while the 'orange' SFD plays the role of the shifted distribution.Interchanging the roles of the two will not alter the final result.In the rightmost panel, there are many replications of the orange curve, shifted such that their maxima (marked by the black points) always lie on the blue curve.The envelope of this construction corresponds to the desired product.

Product of two independent random variables
To calculate the SFD of the product of two i.r.v.s X ∼ N −ξ and Y ∼ N −η , let us consider the corresponding convolution in terms of ξ and η: Calculating this integral in the saddle-point approximation, we arrive at the expression suspiciously similar to the one from the mix rule: And this is not a surprise: the product can be understood purely in terms of the previously introduced mix rule (Sec.3.3).To see this, notice that multiplication by a constant N −s results in a horizontal shift of the multiplied SFD by the distance s.Thus, the multiplication by an arbitrarily distributed random variable X = N −α can be now considered as a superposition of different horizontal shifts α (multiplication by a constant) with different vertical shifts f (α)−1 (the constant's weights), restoring the same mix rule, see Fig. 5.

Extensive sums of i.i.d random variables
In Sec.3.2, we have already seen how to calculate a spectrum of fractal dimensions of a sum of two independent random variables.It can be easily generalized to a sum of any finite number of r.v.s.At the same time, we can imagine that an extensive sum of r.v.s has to have something to do with its typical value, which does not follow from the finite sum rule, e.g., a central limit theorem predicts a square-root dependence of the typical value with the number of terms in the sum.To make the required generalizations, let us consider a non-negative random variable X defined by its spectrum of fractal dimensions f X (α) and a random variable S defined as a sum of N β independent copies of the random variable X : Our task now is to calculate f S (α).
We start by focusing on such values of x = N −α that appear in any typical sample of size N β .Specifically, let us consider α such that its typical number of realizations one can estimate a probability to find exactly N g i −1+β δα terms of the order of N −α i in a particular realization of the sum as exp{− N g i − N f X (α i ) 2 /2N f X (α i )+1−β δα} which is double-exponentially small in g i , provided g i ̸ = f X (α i ).Thus, neglecting the double-exponentially small contributions, one can write a typical value of the sum as [22,44] meaning that α 0 (S) = min Ω {α − f X (α) + 1 − β}.In addition, one can deduce that f S (α > α 0 (S)) = −∞, which is just a consequence of the 'zeros-suppression' effect mentioned in Sec.3.2 in its extreme case.At the same time, the values of α corresponding to f X (α) + β < 1 with α < α 0 (S) are unlikely to be represented in a typical realization of the sum even by a single term.However, in case of such an unlikely event, the whole sum will be determined by this very contribution.Thus, these rare events should be handled with the mix rule (Sec.3.3): a probability density for such an event to contribute equals to N β p α(X ) (α) ∼ N f X (α i )+β−1 ≪ 1, resulting in the following graphical rules for obtaining the SFD of an extensive sum: 2. Draw a unit-slope line having exactly one common point with f X (α) + β in the area where f X (α) + β ≥ 1.
3. A point α 0 where this unit-slope line crosses the horizontal level f = 1 is our new typical value, according to (16).There is nothing to the right of this point due to the zerossuppression effect.

4.
To the left of this point f S (α) just equals f X (α) + β.
As one can see from this construction, all tails with a slope greater than unity eventually die as β → ∞, in agreement with the standard central limit theorem.In addition, in agreement with the generalized central limit theorem for the one-sided stable distributions, in order for distributions to be stable, the tails of their PDFs should decrease not faster than ∝ s −2 , which is exactly what they do, see Sec. 5 for more details.For such distributions, the unit-slope line touches f X (α) + β exactly at f = 1 for any β > 0, thus, f S (α) never develops discontinuities, and the tail never dies.
We have just derived the rule for extensive summation in the graphical language.For the derivation of the same result in a more traditional mathematical fashion, see Appendix A.

Gaussian Rosenzweig-Porter model
The Gaussian Rosenzweig-Porter ensemble is essentially just an ensemble of N × N Gaussian orthogonal (or unitary) random matrices with broken rotational symmetry: where the elements of H GOE/GU E are i.i.d.Gaussian r.v.s, with zero mean and unit variance.In this section, we show how the graphical rules defined above can help self-consistently calculate the spectrum of fractal dimensions of the local density of states of the model from the first principles.We do it using the cavity method in its diagonal approximation [39,45].The idea of the cavity method is to use an exact expression relating diagonal elements of an N ×N Green's function , where H and H (i) differ by a single site i (H (i) has a small "cavity"): here, z = ϵ − iη is a complex-valued parameter with a small imaginary part η to ensure the existence of G(z).For η > 0, one can define a local density of states ν i (ϵ) as Following [45,46], we choose η = N β δ ϵ with 0 < β < D 1 and δ ϵ being a mean level spacing in the corresponding part of the spectrum.Such a choice allows us to get meaningful physical results for any system in a delocalized phase.Unless otherwise noted, we consider bulk with The idea of the diagonal approximation is to say that, in a thermodynamic limit, the sum in the denominator of ( 18) is dominated by the diagonal elements of the reduced Green's function [45,47]: The approximation is considered valid, provided the system is in a non-ergodic phase. 6However, the SFD's graphical algebra introduced above can not be directly applied to this expression as it also contains complex variables.To proceed, we need to either generalize our graphical algebra to complex random variables or write the local density of states explicitly staying in the real domain.While the former approach is also possible (see Appendix D), for simplicity, here we employ the latter one: In this expression, we neglected the real part of the diagonal self-energy j j (ϵ) compared to the on-site disorder amplitude 7 and omitted η compared to the broadening Γ i (ϵ).These simplifications, again, restrict the applicability of the following results to the non-ergodic delocalized part of the phase diagram.Now, the problem of calculating the spectrum of fractal dimensions of the local density of states ν i (ϵ) appears as a self-consistent problem [31,44,45,47], and below, we show how to solve it using the SFD graphical algebra introduced earlier.
First, let us consider the broadening Γ i (ϵ).Since the broadening can be expressed as an extensive sum, its SFD possesses all the characteristic properties of extensive sums, allowing us to guess the SFD almost completely.Indeed, since this is an extensive sum, there can be nothing to the right of α = α 0 (Γ i ) corresponding to the broadening's typical value, i.e., f Γ i (α > α 0 (Γ i )) = −∞.On the other hand, the terms entering the sum do not correspond to any fat-tailed distribution: V i j is Gaussian by definition, and the LDOS is bounded from above by inverse amplitude of the on-site disorder, which is of the order of N 0 .This means that there can be nothing also to the left of α 0 (Γ i ), i.e., f Γ i (α < α 0 (Γ i )) = −∞.As a result, its SFD should look like The only unknown parameter left is the value of this typical value Γ i ∼ N −c , which we parameterized by α 0 (Γ i ) = c.
The site index i is used in the previous paragraph in Γ i and ν i to specify the random quantities, corresponding to this site.Due to the statistical homogeneity, assumed by the selfconsistent cavity formulation of the problem and respected by the model under consideration, further, we omit this unnecessary (double) indexing and write just Γ , or f Γ (α), or, similarly, p ν (x), where possible.Therefore, the index-free symbols Γ , ν, etc., should be considered as shortcuts for "level broadening", "local density of states", etc.Now, having the shape of the broadening fixed, we can calculate the distribution of ν, substitute it to the definition of Γ , and find c self-consistently.To perform the first step, one should ensure that all terms entering the expression for ν are independent -which is not true because Γ enters the expression two times.However, since, from the SFDs' point of view, Γ is just a constant, we can still perform this step without introducing any further complications.
The calculation is depicted below: ← stretching, see Sec. 3.1; (24) As one can see from a comparison with the result obtained in [38], the shape of the SFD we have just found is correct.Now, let us write the self-consistency equation: ← an SFD of the normal distribution; (28) In the last equation, (31), we demonstrated the elevation of the highest blue point with the small dashed line.Further we will use the same notation.Comparing (31) to (22), one can find c = γ − 1, in complete agreement with the previously known results.The ergodic and Anderson transitions then follow from the equations c = 0 (Γ ∝ N 0 , thus, all ν i are roughly of the same amplitude) and c = 1 (Γ ∝ δ ϵ and, thus, the normalization of ν is contributed by a finite fraction of large components), giving, as expected, While the result ( 27) is not new, its graphical calculation already provides some insights.For example, the orange linear region with the slope 1/2 between α = ±c comes from the fact that our diagonal matrix elements are homogeneous in the bulk of the distribution and, thus, this should be a general feature of many similar random matrix models. 8The only contribution from the off-diagonal elements to the final result (27) was in the form of fixing the value of c, which is reflected by the different colors we use for different contributions.The shape itself was controlled only by the statistics of the on-site disorder and, in particular, by the fact that p ϵ i (x) is finite and non-zero as x → ϵ.
While due to the introduced above restriction Γ ≫ η ≫ δ ϵ , we cannot apply our method in the localized phase, γ > 2, let us see how the Anderson localization formally looks like from our graphical self-consistency equation's point of view.When γ approaches 2 from below, the blue points at α = 2γ − 2 − c and α = γ − c from (31) approach each other until they coincide for γ = γ AT = 2.For γ > 2, the same picture looks like which already contradicts our initial guess (22) for the level broadening.Nevertheless, one can try to use (33) as a new guess, which results in a self-consistency equation c = γ + c − 2. The unknown c drops out of this equation, leaving us with the only point this construction can hold for, γ = 2. And, while this attempt does not lead us to a correct solution, it hints at an important conclusion about the localized phase: the level broadening for γ > 2 no longer has the form (22), originated from the collective contribution of many sum's terms |V i j | 2 ν j .Instead, as one may assume from the previously known results [38], it should also contain an individual contribution from the localized wave function's maximum to preserve the LDOS normalization condition.

Lévy Rosenzweig-Porter model
Inspired by the success of the Rosenzweig-Porter model with normally distributed off-diagonal amplitudes and by the distribution of the effective matrix elements in the Hilbert space, derived for many-body disordered models [49], a similar RP ensemble with the fat-tailed Lévy-distributed amplitudes was proposed [30,31].The main motivation was that this Lévy Rosenzweig-Porter should host the desired multifractal phase.This section reviews the statement and provides our own analysis of the proposed model.
The Lévy Rosenzweig-Porter ensemble is the ensemble of random matrices (17), with the uncorrelated diagonal on-site energies ϵ i , distributed according to some narrow sizeindependent distribution like the normal or box distribution, and with the i.i.d.off-diagonal elements V i j of typical amplitudes N −γ/2 , distributed according to the following PDF with the parameter µ: Such a polynomial decay of the PDF tails makes the off-diagonal elements' distribution heavytailed for µ < 2 (variance is undefined) and fat-tailed for µ < 1 (mean of the absolute value is undefined).The SFDs of the distributions with such polynomial tails look like These heavy-and fat-tail properties are precisely what led the authors of [31] to an elegant argument supporting the existence of the multifractal phase in the Lévy-RP models.The argument is based on estimating the fractal dimensions D 1 and D ∞ and concludes that they are not equal, meaning the wave functions are multifractal.
Below, we calculate a fractal spectrum of the bulk eigenstates of the Lévy Rosenzweig-Porter (Lévy-RP) model, analogously to how we did it for the Gaussian Rosenzweig-Porter model in Sec. 4. As we are mostly interested in the (multi)fractal phase, which is expected [31] to exist for 1 < µ < 2 and 2/µ < γ < 2, this is exactly the parameter range we consider.Thus, we start by defining the fractal spectrum of |V i j | 2 and trying to guess the SFD of the broadening Γ .Rescaling (35) according to the exponentiation rule from Sec. 3.1 and using α 0 = γ/2, we get A product of |V i j | 2 and the LDOS ν j is at least as heavy-tailed as |V i j | 2 .Taking also into account that ν is normalizable, and, i.e., bounded, we conclude that the extensive sum (37) where the right-wing tails disappear, as usual, due to the effect of zero suppression in extensive sums.The only unknown here, like in the Gaussian RP case from Sec. 4, is the scaling exponent c of the typical value of Γ .
Next, to obtain the SFD of ν, we use the mix rule from Sec. 3.3.Indeed, taking for each fixed Γ -value the corresponding conditional distribution (27) and the distribution of 'weights', given by (37), we get where the new 'zeros' part to the right of c originates from the rare realizations of The label with the left arrow demonstrates that f ν (c + 0) = 1 − µc.To finish the calculations, we write the self-consistency equation: As one can see from (38), the SFD of the LDOS in the Lévy-RP model, similarly to the Gaussian RP case, corresponds to such f (α) that f (α < α 1 ) = −∞ and, hence, it is just fractal, not multifractal.Moreover, calculating the fractal dimensions D 1 and D ∞ using this SFD 9 and the self-consistent value of γ e f f , we get which coincides with the results for D 1 from [31], but not with the ones for D ∞ .The fractal phase with 0 < D 1 < 1, as in the Gaussian RP case from Sec. 4, gives a way to the ergodic one as soon as c = 0, provided while the localized phase appears when c = 1, giving As one can see from these equations, the purpose of introducing γ e f f was in preserving the analogy with the Gaussian RP model, as the fractal dimension takes a universal form D q>1/2 = 2 − γ e f f as well as the corresponding phase diagram (32).At this point, we are ready to draw a phase diagram of the Lévy-RP model according to its LDOS SFD.To do that, in addition to what we already did, we need to explore the regions µ > 2 and µ < 1.The former case of µ > 2 reproduces the results of the Gaussian RP: It follows from the fact that the hopping distribution ceases to be heavy-tailed, the slope µ/2 from (39) becomes larger than one, and, as a consequence, the extensive sum from (40) produces the same expression for c * as in the Gaussian RP model, giving the ergodic and the Anderson localization transitions at γ E T = 1 and γ AT = 2 for any µ > 2. The latter case of µ < 1 is a bit more interesting: similarly to the situation with γ > 2 from the end of Sec. 4, c drops out of the corresponding self-consistency equation, giving the Anderson localization transition line as γ AT = 2/µ.However, the support set dimensions (42) on this line are now not zero but one, meaning that the line corresponds also to the ergodic transition.The resulting phase diagram is given in Fig. 7.
We already said that, due to the heavy-tailed distribution of the off-diagonal elements, the convergence of numerical simulations to the thermodynamic limit of the Lévy-RP model is very slow, see also [44].However, an attempt to verify our analytical prediction can be seen in Fig. 8: an extrapolation to infinite sizes [12,38,50] is shown by black dots, and our theoretical prediction is by the thick red line(s).We have used different kinds and directions of extrapolation to diminish the finite-size effects in the most efficient manner, please see the caption of Fig. 8 for more details.
The meaning of the two red lines, dashed and solid, to the right of the typical value, α > α 0 = 0.55, is the following: the dashed line is plotted according to (38), while the solid line shows the actual behavior of the SFD in this region supported by the extrapolation of our numerical results.One can see that the analytical calculations, described above in Eq. (38), are confirmed by the numerical simulations for all α < α 0 below the typical value.The main difference between numerics and analytical predictions appears for α > α 0 .Here one should mention that for the Lévy-RP model this issue does not affect the values of D q even for negative q as in both cases of the behavior f (α > α 0 ) the fractal dimension diverges at q < −µ/2.Nevertheless, let us consider the reason why the calculations above failed to capture the behavior at α > α 0 .This should originate from one of the approximations we made on the way.Here, we consider them one by one and point out which is the most crucial.Based on the estimation from Appendix B, we can conclude that the diagonal cavity approximation holds in the parameters range 2/µ < γ < 2, which we defined at the beginning of Sec. 5.This range corresponds to the non-ergodic extended phase and, thus, it is not this approximation, which is responsible for the inconsistency between (38) and simulation in Fig. 8.
As a result, the only other approximation which might be crucial is the one we made following [31]: it was throwing away the real part of the self-energy compared to the diagonal disorder.This approximation seems reasonable until the typical value of Re Σ i is less than that of the on-site disorder and the self-energy distribution is less heavy-tailed than the on-site disorder.Indeed, provided the conditions above hold, the SFD of |ϵ i + Re Σ i | will be identical to that of |ϵ i |.This is not the case for the Lévy-RP model.
Note that, strictly speaking, to prove the above statement, one should generalize the notion of SFD from the distributions on the positive axis to the ones on the entire real axis (see Appendix C for details of this approach).However, here we provide the simpler reasoning.Indeed, first, the shape of SFD for the variable as one of the terms will dominate the sum and determine its sign.And second, the linear decay of the SFD with the unit slope at α > α 0 (| Re Σ i |) follows from the finite value of the PDF for |ϵ i + Re Σ i | at zero (analogously to the Porter-Thomas distribution in GOE/GUE ensembles), guaranteed, in turn, by the mutual independence of ϵ i and Σ i .As a result, since the PDF of |ϵ i | is also finite at zero, one can neglect the difference between the SFDs of Returning back to the Lévy-RP after this comment, we have to conclude that the heavytailed hopping distribution violates the second condition because it is more heavy-tailed than the on-site disorder distribution.This prevents us from neglecting the real part of self-energy and is the source of the discrepancy at α > α 0 .
More concretely, the above observation forces us to consider both real and imaginary parts of the self-energy or, equivalently, of Green's function G ii , making it necessary to generalize our SFD algebra of independent r.v.s to that of two correlated ones or, equivalently, to the complex-  ⃝-2 ⃝ and 3 ⃝-6 ⃝: in the former case, the straight lines are plotted through the points with the same derivative q = f ′ ν (α) for q = −2, −1, 0, 2.5, 3, and the lines' intersections are highlighted by the black points, while, in the latter case, the ansatz f ν (α(N ); N ) ∼ f ν (α; ∞) + A/ ln(N ), α(N ) ∼ α + B/ ln(N ) is applied along the specified directions: at fixed f (α), A = 0, for 3 ⃝, 5 ⃝, at fixed α, B = 0 for 4 ⃝, and at a certain ratio A/B for 6 ⃝.
valued random variables.This broad and difficult topic is covered extensively in Appendix D, but the take-home message from that generalization is encouraging: the real part of the selfenergy can only contribute to the atypically small values of LDOS, α > α 0 .This fact can also be understood intuitively: while we increase the denominator of ( 21) keeping Γ i fixed, we decrease the value of the LDOS.The smallest value of the LDOS we get in such a way is its typical value, which is realized by a typical value of |ϵ − ϵ i | ∼ O(1).Finally, as Re Σ i starts to dominate when it is larger than O(1), it only contributes to the SFD part to the right of α 0 .Hence, we can continue neglecting it, provided we are only interested in D q with q > 0, given by α ≤ α 0 .
In order to compare our results, Eq. ( 38) and Fig. 8, with the previous claims, let us remind the results from [31].In [31], the authors have considered in detail the miniband structure of the Lévy-RP model, and, both with the phenomenological arguments and the cavity method, derived the fractal dimension D 1 , the results for which have been consistent from both methods as well as been numerically confirmed.The derivation of D ∞ has been done only within the phenomenological arguments with the assumption that it is governed by a typical miniband broadening.It is well-known that the numerical studies of D q for large q are usually very difficult due to the exponentially smaller statistics and very strong finite-size effects.Here, Figure 9: A spectrum of fractal dimensions of intensities of hopping matrix elements in the log-normal RP model.
we avoid having any phenomenological arguments and calculate the spectrum of fractal dimensions directly from the cavity method.The results show indeed the controversy to the phenomenological arguments of [31] for q > 1, including the typicality argument of the miniband broadening for q → ∞.

Log-normal Rosenzweig-Porter model
Another candidate to host a genuinely multifractal phase, circulating in the literature, is the socalled log-normal Rosenzweig-Porter model [3,18,22].Its definition is analogous to other RP models, but the distribution of hopping matrix elements is defined as a real-valued log-normal distribution, with the parameters scaling with the system size N : Proceeding according to (1), we find that its SFD is parabolic, and hence truly multifractal, see Fig. 9. Indeed, if we want a truly multifractal wave function, why not to start from a truly multifractal hopping distribution?As we usually do at the beginning of the graphical calculations, let us define the range of parameters we consider and guess the first step's SFD.Defining the hopping SFD by the relation 4a), let us start from considering 1 < γ < 2 and small a = pγ > 0. The 'smallness' of a will be defined later.But, considering that, for a → 0, the hopping distribution approaches a narrow distribution with all moments well-defined, it is reasonable to assume that the LN-RP LDOS in this regime will be close to the Gaussian RP LDOS.With this idea in mind, let us start from the Gaussian RP LDOS SFD (27) and multiply it by the squared hopping element, with the log-normal distribution: (46) One can see how the orange segment with the slope 1/2 appears smoothly in the blue parabola between the points A and B, where the parabola's slope is also equal to 1/2.Thus, the values of A, B, and ????????????? calculate the level broadening The gray dot located at X = γ− c − 2a marks the point on the SFD (46) where its slope is equal to one.The tangent line with a unit slope, touching this point, crosses the level From the shape of the SFD for the level broadening, we see that the resulting self-consistent LDOS SFD has the RP-like shape for α < γ e f f .At the same time, for α > γ e f f it should show a very low probability of zeros.This is the only place where the parabolic input reveals itself.Formally speaking, this SFD is already multifractal for D q with negative q, but from the physical application perspective and according to the definition we gave at the end of Sec. 2, we consider it as a trivial fractal phase.From Eq. ( 48), one can find the part of the phase diagram at γ < 2, Fig. 10.Indeed, the case of c = 0 corresponds to the level broadening of the order of the entire bandwidth which should correspond to the ergodic transition.This gives In the r.h.s.we have used the standard substitution a = pγ.let us now consider the question of the parameter's a smallness.The geometric construction pictured in (47) holds only until X > c * .This implies the restriction γ < 2. Given that, we can now draw a part of our model's phase diagram, Fig. 10.
Next, let us try to move to higher values of γ.To do that, let us start with some γ in the fractal phase, 1 + a < γ < 2, and gradually increase it until we reach the unexplored territory.
When the gray dot from Eq. ( 47) crosses the horizontal dashed line f (α) = 1, we can no longer determine c * as a crossing point between the unit-slope tangent line and the level f (α) = 1.Instead, we have to use the following construction: From this we find that c * = γ − c − 2 a(1 − c), and the self-consistent solution for c is now Since the part of the SFD for the level broadening to the left of c does not have even a single point with a slope less than or equal to 1/2, the resulting SFD for the LDOS in the LN-RP model again has the shape of the LDOS SFD for the Gaussian RP, except for the part to the right of the typical value α = c.Note that after the substitution of a = pγ Eq. ( 51) agrees well with (49) in [3] and (C3) in [31] for τ * ≡ α(γ, p) ≡ c.The square root in the self-consistent expression (51) for c immediately tells us that the result is valid only until γ ≤ 2 + a/2.From the geometrical point of view, the line γ = 2 + a/2 corresponds to the situation when the orange segment of ( 46) touches the level f (α) = 0. Notice a similarity with the case γ > 2, discussed at the end of Sec. 4 in the context of the Gaussian RP model: the change in geometry happening for γ > 2 + a/2 calls to modify the self-consistency equation once more, but, if one tries to do it, one would quickly realize that it is impossible as c drops out of the equation, leaving us with just the expression for this borderline itself.Similarly to the discussion at the end of Sec. 4, this disappearance of c from the self-consistency equation hints that one of the basic conditions cannot be satisfied, namely, the wave-function normalization condition, leading to the emergence of the localization peak f ν (α = −1) = 0 and signaling the Anderson transition.Thus, the expression for the Anderson transition from the fractal phase is given by where the latter expression is straightforwardly derived from γ = 2+ a/2 after the substitution a = pγ.Especially intriguing is that, like in the Lévy-RP model at γ = 2/µ and µ < 1, the fractal dimensions on the Anderson transition lines are finite, implying a discontinuity of these quantities.Given that, in the LN-RP case, the line is the borderline of the fractal phase and, hence, the borderline of our method's applicability region, we cannot rule out the possibility of having a fine-tuned multifractality right on this line.
To finish the LN-RP diagram exploration, let us find the line corresponding to the ergodicity transition at γ > 2. To do that, we solve the equation c = 0 and get the latter form, again, is the known expression obtained from the former one by the substitution a = pγ.As a result, summarizing Eqs. ( 49), (52), and (53), we obtain the full phase diagram for the LN-RP model, see Fig. 11, confirming the previous results from [3,18,22,31].It has a tricritical point analogous to the Lévy-RP case.According to [18], the case of RRG corresponds to the bisector a = γ (p = 1) in this phase diagram, which crosses the above tricritical But how is this possible that even the model with a multifractal distribution of its hopping elements failed to host a multifractal phase?One of the possible answers to this question lies in the observable we studied: recall that the local density of states, being an average over an extensive number of wave functions, does not necessarily reproduce all the details of the individual wave function distributions.Thus, as it was shown, e.g., for the Anderson model on the Cayley tree in [19,20], the absence of multifractality in the LDOS SFD does not eliminate the possibility of having the multifractal statistics of the eigenstates.
In order to examine this opportunity of having the wave-function multifractality together with the fractality of the LDOS, we consider the relation between their SFDs in the next section.

Relation between LDOS and eigenstate distributions
We are unaware of any method as powerful as the self-consistent cavity method but for individual wave functions ψ n (i) instead of the local density of states ν i (ϵ).And, while we cannot calculate the corresponding SFD f |ψ| 2 (α) directly, we can infer restrictions for this function implied by the shape of f ν (α).Surprisingly, this analysis can be performed for any random Hamiltonian, and the result of this section goes beyond the Rosenzweig-Porter family.
By definition (19), the local density of states ν i (ϵ) is proportional to the average of the squared wave functions' amplitudes |ψ n (i)| 2 over the energy window η around the energy ϵ.This energy window for the averaging is controlled by the Lorentzian function meaning that its tails decrease as (ϵ − E n ) −2 .First, for simplicity, let us consider a box kernel instead of the Lorentzian one and define the box LDOS νi (ϵ) as Later, we return back to the standard Lorenzian LDOS.From now on, let us fix a specific site i of our system.In the RP family's models, all sites are statistically equivalent, but, in general, they do not have to be.Each realization of a random Hamiltonian H will then produce a single realization of νi (ϵ) = νi (ϵ; H) and of the order of to the value of νi (ϵ; H), Eq. (54).Our first task is to relate the distribution of |ψ n (i) 2 | at |ϵ − E n | < η to the distribution of νi (ϵ).For this, we consider 0 < β ≪ 1 in order to assume that the wave functions in this small energy window are statistically equivalent. 10he values of |ψ n (i; H)| 2 from each realization of H may or may not be correlated.Hence, since our graphical language was developed only for independent random variables, we cannot directly apply the generalized central limit theorem from Sec. 3.5 to this case of Eq. ( 54).Having said that, let us consider the whole ensemble of H and a set Ωα consisting of all |ψ n (i; H)| 2 contributing to νi (ϵ; H) such that νi (ϵ; H) = N −α : The unconditional distribution p |ψ| 2 (x) of |ψ n (i)| 2 from our energy window can be obtained from the conditional distribution Ωα by a probability chain rule as Transposing from the probability density functions to the corresponding spectra of fractal dimensions, we recover the mix rule from Sec. 3.3: This formula directly relates the SFD of the eigenvalues f |ψ| 2 (α) to the SFD of the box LDOS f ν(ξ).Now, let us consider f |ψ| 2 (α|N −α ∈ Ωξ ) and infer what it may look like.First, from the definition of Ωξ , Eq. ( 55), we know that |ψ| 2 from our conditional distribution p |ψ| 2 (x|x ∈ Ωξ ) cannot exceed N −ξ δ ϵ = N −ξ−1 , see Eq. (54).Indeed, otherwise, such |ψ| 2 would give rise to ν > N −ξ at least for some realizations of H related to Ωξ .Second, we know that the point Otherwise, we would get ν < N −ξ for some H contributing to our conditional distribution.Hence, the conditional SFD is constrained to have the following form: Here, the blue part of the graph is the necessary part: f |ψ| 2 (α|N −α ∈ Ωξ ) has to be equal to 1 at α = 1 + ξ and to −∞ at α < 1 + ξ.In the interval α > 1 + ξ, our conditional SFD can have any possible shape allowed by our definition of SFD, f |ψ| 2 (α|N −α ∈ Ωξ ) ≤ 1.Several examples of that are depicted in (58) by the thin lines of different colors.Hence, applying the mix rule to this manifold of conditional SFDs, we get that the unconditioned SFD f |ψ| 2 (α) can differ from the corresponding SFD f ν(α) of the box LDOS, Eq. ( 54), only in the region α > α 0 (ν).
To the left of this point, these two SFDs must coincide. 11inally, let us return back to the proper LDOS with the Lorentzian kernel.It differs from the box LDOS, which we have just examined, because it can be dominated, in principle, by the wave functions outside the energy window ϵ ± η.This possibility lifts the restriction for the point {1 + ξ, 1} to belong to f |ψ| 2 (α|N −α ∈ Ωξ ) for all ξ except ξ = α 1 (ν) − 1 ≡ D 1 (ν) − 1.Indeed, the latter is just a consequence of the wave-function normalization condition.Thus, we arrive at the following conclusion about the relation between the distributions of the eigenfunctions and the local density of states: For the log-normal and Lévy Rosenzweig-Porter models, it means that D q (|ψ| 2 ) = D q (ν) for q ≥ 1/2.This conclusion follows from the fact that, for both of these models, For most practical applications, where only q > 1/2 matters, these models show only fractal, but not multifractal properties.

Absence of multifractality in Rosenzweig-Porter models
As one may have already guessed from the previously considered models, the finding of a multifractal phase hosted by an RP model is far from trivial.In this section, we will prove that it is, in fact, impossible.Before proceeding to the proof itself, let us focus on an important limitation of our method.Indeed, as one can guess, the above-developed method is insensitive to the changes in the PDF that do not affect the SFD.In this sense, all the sparse graph models and the standard Anderson models on the lattices cannot be described by this method, as the localization transition and the fractality are not governed by the scaling with the system size N .
As a remarkable example, let us consider a random Hamiltonian H = H RP +A ER , where H RP is a Gaussian RP model Hamiltonian from (17), and A ER is an adjacency matrix of a random Erdos-Rényi graph, with a fluctuating finite number of non-zero hopping terms of the order of unity.The hopping of this 'Erdos-Rényi-RP model' corresponds to the SFD hereafter, let us assume γ > 1 + c with a certain positive c > 0. Note that the additional blue point at the origin in the above plot corresponds to a finite number of "neighbors" for any site of this model, connected to it by the hopping of the order one.This is given in addition to the all-to-all RP-like hopping of the scaling with the system size amplitude N −γ/2 , see, e.g., [51] for the case of correlated hopping of this kind.Performing the calculations of the SFD for ν, one can straightforwardly show that for γ > 1 + c any SFD curve, respecting Mirlin-Fyodorov symmetry [52], f ν (α) = α + f ν (−α), which is finite only on the support |α| < c, i.e., f ν (|α| > c) = −∞, satisfies the self-consistent cavity equation for such hopping.For example, one can take the parabolic one What's the catch?As we have mentioned above, for the Anderson model on sparse random graphs, the hopping scaling is not the only thing that matters for the localization and ergodicity breaking.In principle, different on-site energy distributions with the same SFD will lead to different LDOS SFDs, see, e.g., [17] for the RRG.Analogously, the different lattice dimensionalities of the standard Anderson model drastically change the localization diagram [53,54].Thus, from the perspective of Laplace's method in its leading order, the problem is ill-defined for such models, as it is not the SFD of the off-diagonal matrix elements and N -scaling, but PDFs and prefactors that resolve the localization phase diagram.This leads to the solution's ambiguity within the above graphical method and its inapplicability to such problems.In the following paragraphs, we focus on the models, where a multifractal segment in the LDOS SFD necessarily originates from the hopping SFD.This assumption implies that the resulting solution is solely determined by the SFDs of hopping terms and on-site energies.
That being said, let us prove that conventional RP-like models, i.e., the models differing from the Gaussian RP only by the distribution of the i.i.d.uncorrelated hopping elements without the Erdos-Rényi component, cannot host any multifractal phase.To do that, we are going to exploit a similar anatomical approach to what we used in Sec.7: we assume the solution found, trace back its features to the input distributions, and conclude if such a selfconsistent solution can actually exist or not.
So, for concreteness, let us start with the shape of the LDOS SFD given, e.g., by ( 61).As will be shortly seen, this particular choice does not affect the argument.Being a part of the iteration procedure, this shape ought to originate from mixing together different Gaussian-RPlike Lorenzian shapes of LDOS ( 27) corresponding to different fixed values of Γ : please see the orange dashed line in (62) for Γ ∼ N −ξ and recall, e.g., how we obtained (38).Schematically, this inheritance can be illustrated by ν : here, the blue part of the LDOS SFD corresponds to the blue part of the broadening SFD, and, since, due to the Mirlin-Fyodorov symmetry, the magenta part of f ν (α) is controlled by the same blue part of f Γ (α), the thin differently-colored curves of the SFD for the broadening demonstrate the unimportance of f Γ (α < 0) as it can only affect f ν (α > c).The fact that the blue region of the broadening SFD produces the identical region on the LDOS SFD is due to the derivative of f Γ (α) in this region being smaller than 1/2.Indeed, as soon as it becomes larger than 1/2, the corresponding contribution becomes subdominant with respect to the contribution of independent on-site energies, leading to the orange straight lines ∝ α/2 of the Poisson distribution we have already used to.
In its turn, the broadening SFD is obtained from the SFD of |V i j | 2 ν j by the extensive summation according to Sec. 3.5.Because of the zeros suppression effect, its tail, f Γ (α < α 0 (Γ )), may only originate from the tail of In our case, it must look like here, as usual, the solid horizontal axis marks the level f (α) = 0. What happens to the right of α = c is, again, unimportant, as long as c stays the typical value of the extensive sum To answer this question, we need to consider the product rule from Sec. 3.4.According to this rule, we know that, if the SFD of the product contains a point with some derivative, a corresponding point with the same derivative must be present on at least one multipliers' SFD.Moreover, if both multipliers contain points with the same derivative, the corresponding point on the product's SFD has a well-defined value and position, e.g., in our case, From ( 62) and (63), we know that, for any point from the blue region, α q (|V | 2 ν) = α q (ν), and leaving us with the requirement for f |V | 2 (α) to pass through the point {0, 0} and to have a discontinuous first derivative at this point.But it means that the multifractal shape of f ν (α) we assumed in our example originated not from the hopping or on-site energies SFDs but from something more subtle like specific PDFs or prefactors, i.e., it takes us beyond the conventional RP models' family and, thus, proves the absence of LDOS multifractality inside it. 12Finally, we use the results of Sec. 7 to conclude that eigenstates of RP-like models, following the LDOS, also cannot be multifractal from the point of view of fractal dimensions D q with q ≥ 1/2.

Conclusions and discussion
To sum up, in this paper, we have developed a powerful graphical method to calculate the multifractal properties, such as the spectrum of fractal dimensions and self-energies of the local density of states, in various models of the Rosenzweig-Porter type.We consider the Gaussian, Lévy, and log-normal Rosenzweig-Porter models and, with our method, easily reproduce their phase diagrams and fractal dimensions D 1 .
In addition to that, we have calculated the entire spectra of fractal dimensions f (α) for the local density of states and found its relation to the one for the eigenstate coefficients.As a result, we have explicitly shown that, in all such models, the phase diagram may suggest the only type of the non-ergodic extended phase, which is fractal, but not multifractal, if we are speaking about the positive integer orders q > 0 of the fractal dimension D q .The only formally multifractal part in f (α) we have found is beyond the point with the tangent slope 1/2, α > α 1/2 .This corresponds to the small or even negative moments q < 1/2 of the fractal dimensions D q .
Having these calculations, we have managed to track back the origin of all parts of the spectrum of fractal dimensions and concluded that the uncorrelated models with i.i.d.hopping terms and conventional Poisson disorder can host only fractal phases.This statement can be, in principle, generalized to the Erdös-Rényi graphs with random hopping amplitudes and finite fraction p ∼ O(1) of non-zero edges.Statistically non-homogeneous distributions of the on-site disorder (like in [44]) may lead to a zero fraction of multifractal states, but the question if it can lead to the formation of an entire multifractal phase remains open.Another possibility for creating genuine multifractality is adding hopping-term [55] and on-site disorder correlations [56,57].In principle, such correlations can be taken into account by combining the ideas of the Sherman-Morrison formula [46] with the graphical method, developed in this work.
The absence of multifractality, though, does not exclude non-trivial and anomalously slow dynamics in Rosenzweig-Porter models, see, e.g., [3,30], which has a direct application to many-body disordered systems close to the MBL transition.The generalization of the developed graphical method for the effects of correlations of the local density of states may become a good way to map such many-body systems to their random-matrix proxies.In this direction, one of the most prominent things is to focus on the frozen dynamical phase, suggested in [3], where the return probability of the wave-packet spreading can be stuck after some finite-time evolution.
Another direction to look at with the developed method is to consider a generic multifractal Rosenzweig-Porter model, obeying the RRG symmetry (see Eq. ( 6) in [3]) and focus on the origin of the tricritical point, found in the phase diagrams of Lévy-and log-normal RP models.
Finally, since our approach relies on Laplace's method of approximate integration, it does not only lead to a self-consistent solution in the thermodynamic limit N → ∞ (which has been done in this paper) but also may allow calculating sub-leading orders of the finite-size scaling for N ≫ 1.In this case, the cavity equation is not supposed to be solved self-consistently but to be viewed as a generator of the RG flow going to the known self-consistent fixed point.Among others, this point of view provides a way to analyze how to minimize the finite-size effects and speed up the convergence of numerical methods.

𝛽 𝛼 − (𝛼
Figure 12: A given spectrum of fractal dimensions f (α) (in red) and one of the possible realizations of the coarse-grained empirical SFD histogram g i .The dashed green line with the unit slope shows the value of s resulting from this particular histogram realization.The fact that this particular realization of s is smaller than N −ξ follows from that the blue histogram lies below the dashed orange line with unit slope α − ξ.The plotted configuration illustrates the case when α−ξ > f X (α)−1+β for all α > ξ.The opposite case with α − ξ < f X (α) − 1 + β for all α > ξ differs by the fact that, instead of only one bin being overfilled, it requires multiple bins being simultaneously underfilled, which is much less probable and leads to the (double-)exponential suppression of zeros, f S (α > α 0 ) = −∞.
we can define a probability that s is less than N −ξ as a probability that all g i lie below the line α − ξ, see Fig. 12.So now, we should calculate the likelihood of the empirical counts n i to be described by a vector n, and then sum up the probabilities of all realizations respecting the condition above: In principle, this approach could work, but in practice, it probably does not because of the total count conservation condition i n i = N β .Indeed, given that the probability to have x ∝ N −α i is p i ∝ N f X (α i )−1 ∆α, the probability P(n| i n i = N β ) itself is easy to write down: However, if we want to sum it over different n, we would have a hard time doing that because of the factor N β ! and inability to sum up all factors p n i i /n i !independently.That is why, following the ideas of classical statistical mechanics, we now introduce a variable "number of particles" Σ = Σ i n i and a narrow distribution over different values of Σ peaked at Σ = N β .Choosing this distribution to be the Poisson distribution we cancel the nasty factorial and leave ourselves with the much easier expression to operate with: here, we introduced ν i = p i N β as the expected filling of the i'th bin.For sufficiently large N β , the Poisson weight (A.6) becomes asymptotically close to a Gaussian one with mean N β and the standard deviation N β/2 , meaning that the probabilities of having n such that Σ i n i = N β ′ with β ′ ̸ = β are double-exponentially suppressed in the thermodynamic limit.Hence, being interested only in the large-N β scaling behavior of the distribution of S, we can use (A.7) instead of (A.5) and write P(s < N −ξ ) from (A.4) in a large-N β limit as where M i (ξ) = N α i −ξ ∆α is the largest possible count in the bin i, obeying the condition is the upper incomplete Gamma functions, and is a probability of all "prohibited" for the given ξ events.The only thing left to do now is to carefully write down the asymptotic expressions of the Gamma functions and go to the thermodynamic limit.First, consider the case when α − ξ > f X (α) − 1 + β for all α > ξ or, in other words, when M i (ξ) ≫ ν i for all i|M i (ξ) > 0, see, again, Fig. 12.In this limit, we can write where γ(1 + M i (ξ), ν i ) is the lower incomplete Gamma function.In the given limit, all terms in the sum are exponentially small: Moreover, since ∆α does not scale with N and the number of terms stays the same as N → ∞, and because of this increasingly small M i (ξ) −M i (ξ) factor, each subsequent term is exponentially smaller than the previous one.In other words, where ν(1) = N f X (α i )−1+β ∆α is the expected filling at the point where M i (ξ) = 1, or, in other words, at α i = ξ − log N (∆α) → ξ.Since in this case f X (ξ) − 1 + β < 0, ν(1) ≪ 1, and Taking into account also the factor e −N β p(ξ) , we get and, thus, Now, consider the second case, i.e., when the line α − ξ has one or more intersections with f X (α) − 1 + β in the region α > ξ.Because in this case, there are such i at M i (ξ) > 0 that ν i ≫ M i (ξ), the probability P(s < N −ξ ) becomes exponentially small.As written at the end of the caption to Fig. 12, this can be understood in simple entropic terms.Alternatively, it can be extracted from (A.8) and (A.9), using the Gamma function asymptotic, similarly to how we have done it for the intersection-free case.Moreover, as ξ and f X (ξ) enter all the expressions as the exponents of N , this smallness is, in fact, double exponential, meaning that P(s < N −ξ ) ∝ exp N −ξ .Thus, without any involved calculations, we can say that, in this case, f S (α) → −∞ as N → ∞.This result manifests the absence of zeros for extensive sums of non-negative random variables, which agrees with the suppression of zeros reported for finite sums in Sec 3.2.

B Applicability of the diagonal cavity approximation to the Lévy-RP model
The diagonal cavity approximation consists in substituting the exact self-energy j,k̸ =i V i j G (i) jk V ki with the diagonal cavity self-energy j̸ =i G (i) j j |V i j | 2 .This approximation relies on the assumption that, for non-ergodic states, typical off-diagonal terms of the Green's function G are negligible compared to the diagonal ones, see, e.g., [45,47].Such assumption is indeed quite reasonable, and it can be seen from the following simple picture: if all eigenstates are non-ergodic and, hence, each of them occupies N D 1 sites, which is a measure zero of all sites for D 1 < 1, then the probability to find 〈i|n〉 〈n| j〉 of the order of N −D 1 is of the order of N −2(1−D 1 ) ≪ 1.In other words, with the probability N 0 this off-diagonal projector element is much smaller than the local density of states on the support set.However, since there are many more off-diagonal elements than diagonal ones, let us proceed to more rigorous reasoning.
To obtain the (asymptotic) equation ( 20), all we need to do is to neglect, for some reason, the off-diagonal terms in the exact block matrix inversion formula (18).This simplification was used in numerous papers, including [31,45,[58][59][60][61].However, it is not so simple to justify the simplification in each particular case.Following [59] and [58], one may write the analogous to (18) expression for the off-diagonal Green's function elements and try to estimate the off-diagonal contribution based on the generalized central limit theorem.Assume the off-diagonal contribution is negligible: for V i j following a distribution with the PDF p(x) decreasing like ∼ x −1−µ and a typical value equal to N −γ/2 , we get and V i j G (i) Comparing the terms in the r.h.s. of the last equation, we get that the assumption is true for γµ > 2 and false otherwise. 13Still, the diagonal approximation is sometimes used even for γµ < 2, as, e.g., in [31].We do not know any analytical justification for that.

C Graphical algebra of (SFD-)symmetrically distributed r.v.s
Before trying to generalize the rules from Sec. 3 to the case of complex random variables, let us consider a more straightforward case, namely, the case of a real random variable supported on the entire real axis.The problem with this generalization is that α ∈ only describes an absolute value of the random variable |x| = N −α .Thus, ideally, we would need some additional construction to describe the sign-related aspects.For example, one could have written the PDF p(x) of this random variable as w + p + (x) + w − p − (x), where p ζ (x), ζ = ±1, equals zero for ζx < 0 correspondingly, and assign conventional SFDs f ± (α) to each of the one-sided probability density functions p ± (x) separately.However, it turns out that, for the purpose of the present paper, we only need to consider symmetric or, rather, SFD-symmetric distributions.
After restricting ourselves to this case, we only need to adjust the meaning of f (α), its properties, and the rules discussed above.For example, while all even moments of symmetrically distributed random variables can still be found through their relation to tangent lines discussed around Fig. 1, all odd moments in terms of the symmetric SFD are now ill-defined.The exponentiation rule from Sec. 3.1 remains the same, with the only remark stating that even and odd powers now correspond to one-sided and SFD-symmetric distributions, respectively.The product rule from Sec. 3.4 and the ensemble-mixture rule from Sec. 3.3 stay exactly the same, and only the sum rule from Sec. 3.2 requires a substantial modification. 14This modification is given below.
Consider two SFD-symmetric random variables X and Y .Their sum is also SFD-symmetric and, hence, to draw f X +Y (α), it is enough to know the probability density of the absolute value of the sum, which can be expressed as The SFD corresponding to p |X |+|Y | (s) can be calculated using the original sum rule from Sec. 3.2, the SFD corresponding to the weighted sum of different PDFs is given by the mix rule from Sec. 3.3, so the only thing missing here is a rule for the subtraction of non-negative random variables.let us fill this gap.To do that, we write where Then, we substitute |s| = N −α and χ = N −ξ and go from the PDFs to the SFDs: This 'subtraction rule' is important by itself as it provides a way for further generalizations of the notion of the SFD to not only SFD-symmetric distributions.Finally, getting back to (C.1) and assuming, as earlier, that α 0 (X ) < α 0 (Y ) and, for simplicity, that both f X (α) and f Y (α) are convex, we get using the mix from Sec. 3.3 with equal weights ∝ N 0 .As one can see, the new contribution may prevent the suppression of zeros.This is, indeed, expected, considering that the zeros suppression effect originated from summing up strictly non-negative r.v.s.For example, adding together two normally distributed variables X and Y with f X ,Y (α) = (1−α) θ (α), where θ (α) is −∞ for negative α and 1 for positive ones, we get back another normally distributed one as we should: Here, we used (1), given as p X +Y (|s| = N −α ) ≡ N f X +Y (α)+α−1 .

D Graphical algebra of complex-valued Green's functions' SFDs
As one may guess, the task of generalizing the notion of the spectra of fractal dimensions to complex random variables is quite challenging.However, since we do not pursue the goal of being general but the goal of applying the algebra to our cavity Green's functions, we will generalize only those operations we actually need.let us see what the operations are: 1. Sum of independent real and complex random variables.This operation is necessary to compute the denominator of the r.h.s. of (20) and add together the real on-site energy and the complex self-energy.
2. Inverse of a random complex variable.This operation helps to get the Green's function from the sum of the on-site energy and the self-energy.
3. Product of independent real and complex random variables.This one is needed to compute individual terms of the sum corresponding to the (diagonal) cavity self-energy.
4. Sum of the extensive number of i.i.d.complex random variables (c.r.v.s).This is to compute the (diagonal) self-energy from the cavity Green's function and hopping amplitudes.
5. Taking real/imaginary parts of the c.r.v.s.This one is needed to finally get the distribution of the desired LDOS from the distribution of the complex-valued Green's function.
Restricting our interest to this set of operations, let us now define what we will mean by the SFD of a complex-valued random Green's function.The observation allowing us to simplify the definition significantly is the fact that, in the middle of the spectrum, the real part of Green's function is distributed (SFD-)symmetrically, see Appendix C for the definition of SFDsymmetric functions.Because of that, and because the imaginary part always has a one-sided distribution, we can define the desired SFD Because this generalization forces us to consider not planar diagrams but 2D surfaces embedded into 3D space, it requires some time to get used to the notation and find solid ground.To speed up this process, consider a couple of examples.
2D Gaussian distribution By definition, the probability density of the two-dimensional Gaussian distribution is a product of two independent one-dimensional PDFs.Assuming the first one has width N −α R0 and the second one N −α I 0 , we get the following SFD expression for the corresponding complex distribution: Here, as earlier, θ (x) is −∞ for negative and 1 for positive values.Picturing 2D surfaces on 2D paper clearly and concisely is generally challenging.To do that, we use polygon mesh: assuming our surfaces to be polyhedra and to have no curved regions, we plot different facets with different colors, mark edges and vertices with straight lines, arrows, and points, and specify edges' slopes with text labels.The 2D Gaussian distribution we defined above is depicted in Fig. 13a.While the contour plot of the surface would be more general and could depict both curved and polyhedral surfaces, it leads to, in some sense, a "raster" image requiring the specification of as many contours as possible for the complete representation of the surface.In contrast, the notation we propose relies, of course, on the absence of multifractality, but can describe each facet of a fractal SFD by just two vectors.In addition, it appears to be very convenient from the practical point of view: the aforementioned operations in this notation look easier than in any other we considered.
Fully-correlated 2D Gaussian distribution Next, let us see how correlations between our random variable's real and imaginary parts affect the picture.For simplicity, consider the extreme case with Re X = Im X resulting in the Fig. 13b.The behavior along the only edge, which is shown as both −α R and −α I , should be understood as follows: the labels are the coordinate-dependent parts of the parameterization of the SFD value along the edge, and, since it can be parameterized using both α R and α I , both ways of labeling are possible and equivalent.In this particular case of the edge directed at π/4-angle to the axes, the change of the parameterization is particularly simple and leaves the functional dependence intact.
(a) The SFD (D.2) of a c.r.v. with real and imaginary parts being independent Gaussian r.v.s.
The SFD of a complex variable X + iX N α R 0 −α I 0 with normally distributed X .
Figure 13: Graphical representations of the SFDs of c.r.v.s with (a) uncorrelated and (b) fully correlated real and imaginary parts.The solid point denotes the location of the typical (absolute) value of the random variable, while the encircled number gives the corresponding SFD value.The labels ∝ −α R and ∝ −α I show how the SFD behaves along the edges.Notice that the diagonal edge on panel (b) has both labels close to it.This is only to demonstrate their equivalence as any edge having non-zero angles to both axes can be labeled using both α R and α I , up to the author's choice.Below, we will choose one of the two ways of labeling for each specific situation separately.The colorful zone represents a facet with a finite value of the SFD, while the white zone represents the area with f (α R , α I ) = −∞.In addition, the arrows show the direction of the SFD's decay along the edges.This does not provide any additional information but makes it easier to read the diagrams.The value of the SFD on the facet is schematically shown by the color gradient and is only present here for the aesthetic beauty.

D.1 Real and imaginary parts of the complex r.v.s.
Imagine we have obtained the SFD for our complex-valued Green's function.How to extract the information about the local density of states from it?To do that, we need to calculate the SFD f Im X (α) of the marginal distribution of the Green's function's imaginary part.From the definitions (1) and (D.1), we find that 3) the expression for f Re X (α) is analogous.Pictorially, drawing a marginal SFD from the joint SFD is an orthographic projection along the corresponding axis, see, e.g., Fig. 14.

D.2 Sum of the real and complex r.v.s.
A PDF p X (x I , x R ) of any distribution can be written as p Im X (x I )p Re X (x R |x I ), where p Im X (x I ) is a PDF of the marginal distribution of Im X and p Re X (x R |x I ) is a PDF of the conditional distribution of Re X given that Im X = x I .Noticing that r.v.X + Y has the same p Im X (x I ) as X provided Y ∈ , we reduce the task of calculating the distribution of X + Y to the task of calculating the distributions of Y + Re X for different fixed values of Im X which can be done using the already familiar rules from Sec. 3.2 and Appendix C. Thus, the corresponding graphical algorithm is: 1. Split the 2D surface f X (α R , α I ) of the r.v.X's SFD into individual 1d slices with fixed α I .
2. Shift each of the slices vertically such that their maxima become 1.After this normalization, the slices can be viewed as the SFDs f Re X (α R |α I ) of the corresponding conditional distributions.
3. Calculate the SFD of the sum Re X + Y under the condition, posed by each of these slices, using the summation and subtraction rules from Sec. 3.2 and Appendix C.
4. Assemble the 2D surface from the resulting 1d slices restoring their initial maximal heights by the appropriate vertical shifts given by the SFD f Im X (α).
To illustrate the scheme, consider a sum of the complex fully correlated Gaussian r.v.X from Fig. 13b and the real one-sided Lévy r.v.Y from (37) with c > α R0 .The normalized slices with fixed imaginary α I , f Re X (α R |α I ), of the two-dimensional SFD of X together with the corresponding vertical shifts given by the marginal SFD f Im X (α I ) are Notice that the slices are associated with symmetric conditional distributions, but the results of Sec.C cannot be applied directly because the Lévy distribution is not symmetric.However, keeping in mind (C.4), the remark directly after it, and considering the SFD of X + Y as the equal-weight mix of a sum and a difference between the Lévy and the half-normal distributions, we get the following result 15 for the f Re X (α R |α I ) of the individual slices for SFD of the sum Finally, assembling all the slices together with the heights defined by f Im X (α I ), we get Here and further, we omit the specification of vertical and horizontal axes, assuming them always to be α R and α I .For brevity, we also omit the specification of the π/4-angle.From now on, we will always use the same solid-line arc notation to specify such angles.The different colors of the facets mark their relation to the original r.v.s X and Y , allowing, as earlier, to trace back easily the features we see in the resulting distribution.The black-and-white dashed line connecting the vertices marks the discontinuity between the facets originating from the second case in (D.5), and the colors of the vertices mark the facets they belong to: the vertex at {α R 0 , α I 0 } belongs only to the orange facet, while the other vertex belongs to both.Consequently, the two arrows and two non-equivalent labels below and above the discontinuous edge show the slopes for the orange and the blue facets, respectively.

D.3 Product of the real and complex r.v.s
To start with, consider a generic complex r.v.X and multiply it by a constant C = N −c .The resulting two-dimensional SFD of C X is just a shift of the original one: A product of X and a generic real random variable Y can be viewed as a weighted ensemble mixture of such shifts, with the weights controlled by the SFD of Y .Moreover, since all the shifts happen along the lines α R − α I = s = const, one can view the product operation similarly to the previously considered sum operation, i.e., as applying the same independent real operations to different slices characterized by different values of s.
To better understand the operation, consider a product of X sampled from the uncorrelated 2D Gaussian distribution from Fig. 13a and Y sampled, again, from the Lévy distribution (37).For s < s 0 = α R 0 − α I 0 , the normalized fixed-s slices of the complex SFD of X , parameterized by α R , together with the Lévy distribution's SFD, look like for s > s 0 , the diagrams are the same up to the substitution R → I.Then, by multiplying these two SFDs according to the usual rules from Sec. 3.4 and assembling different slices together according to the corresponding weights of each slice, we get The correctness of this result can be again checked by multiplying Y separately by (independent) Re X and Im X and assembling the complex-valued r.v.from its (independent) real and imaginary parts.

D.4 Inversion of the complex random variable
In general, the graphical exponentiation of complex r.v.s is a mess because of the mixing of different phases.However, to calculate cavity Green's function, we only need to learn how to get an SFD of 1/X from the SFD of X , which is doable.The rule for the inversion operation can be deduced from the following two observations: 1. Inversion preserves the argument of c.r.v.mod π.This means that if the point of the original 2D plane corresponded to the tangent N −α I /N −α R , its image will correspond to the same tangent.In other words, the quantity α I − α R conserves for all points under the inversion operation.This leads to the conclusion that all actions again occur along the π/4 slices, as with the product discussed above.
2. Inversion operation inverses the absolute value, which is, in our case, dominated either by N −α R or N −α I .Thus, if the original point lies above the line α I = α R , its α R changes the sign.The same happens with the α I for the points below the line.
The resulting recipe for drawing the SFD of a complex random variable's inverse is shown in Fig. 15.In Fig. 15b, we illustrated how the angles arctan(2) and arctan(1/2) may arise on the diagrams from the inversion operation.Because they will often appear in real calculation, we introduce a special notation for them: we specify the former with a double-line arc and the latter with a dashed-line arc.Below, we omit the text labels and only use this notation.

D.5 Generalized central limit theorem for c.r.v.s
Finally, let us derive the rule for extensive summation of i.i.d.complex r.v.s.For simplicity, we do it in a graphical fashion similarly to how we did it in Sec.3.5, and not in a mathematical style of Appendix A. Thus, we reevaluate the arguments surrounding (16).What this equation says is that each set of terms x i ∝ N −α i , considered individually, sums up to a definite quantity of the order of N −α i N f X (α i )+1−β δα, where N f X (α)+1−β δα is the expected number of such terms.The final result is obtained then as a sum of all such contributions from different α and equals   to the dominant one coming from α i such that f ′ X (α i ) = 1 (saddle point).Now, we are going to use exactly the same logic here: our starting point for the derivation is the analysis of extensive sums of r.v.s defined by narrow SFDs concentrated around a single point on the {α R , α I } plane.In the case of the one-sided distribution with an SFD concentrated around α = γ, an SFD of a sum of N β such i.i.d.random variables would be just another similar SFD concentrated around γ − β.In the case of general SFD-symmetric distribution, the result would be equal to a subtraction of the absolute values of the extensive sums of positive and negative terms, leading to an SFD with the same typical amplitude 16 but also with the linear part ∝ −α to the right of γ − β.In the case we are discussing now, the answer would depend on the nature of real and imaginary parts of our complex random variable.

D.5.1 One-quarter complex distribution.
First, assume both real and imaginary parts of our random summand correspond to one-sided distributions.This situation results in the following expression: Indeed, since, in this case, the sum of identical terms is just a multiplication by N β , it leads to a simple shift of the initial concentrated SFD.And, since the shift is again happening solely along the lines with α R − α I = s = const, we can, in principle, treat the extensive sum of i.i.d.c.r.v.s similarly to how we approached the product and the inverse, i.e., by applying the corresponding real operation to individual s-slices and properly assembling the results together.However, the rules for this 'proper assembling' are quite tricky, and one should be careful while summing up the contributions from different slices.To see why, consider an SFD with not one but two peaks.If they correspond to different values of s and, hence, do not lie on the same slice, the resulting extensive sum may look like On the r.h.s.diagram, the white circles show separate contributions from each peak, and the black circle shows the actual typical value of the sum.Notice that the real and imaginary parts of this value come from different peaks' contributions: its α R is equal to min{α R 0 1 , α R 0 2 } − β, and its α I is equal to min{α I 0 1 , α I 0 2 }−β.In other words, the sum has the same binary nature as the sum of two non-negative r.v.s from Sec. 3.2.This nature is most evident from the expression For summands with SFDs having more than two finite-valued slices, the generalization is similar to the one for the sum rule from Sec. 3.2 to any N -independent number of terms.In other words, using the same coarse-graining argument for the slices-enumerating parameter s as we exploited in Appendix A for α, we arrive at the conclusion that all extensive sums of i.i.d.c.r.v.s must have SFDs topologically similar to here, the double-dot-dashed line {α R (s), α I (s)} represents the manifold of typical contributions from different s-slices, and the right angle shown by the red dashed line is fixed by the two points on the manifold, the ones with α R,I 0 = min s α R,I (s).Consequently, the black point {α R 0 , α I 0 } represents a typical value of the corresponding extensive sum, and the rest of the finite-valued sum's SFD, if any, must lie in a dashed quarter-plane region.This is what the zeros suppression effect looks like on the complex plane.Now, after clarifying the situation with the typical value of the extensive sum, let us find out what happens to the distribution's tails, i.e., how rare events contribute to the sum, and how the dashed quarter-plane region from (D.12) must look like in any particular case.For this purpose, consider i.i.d.summands X with the SFD f X (α R , α I ), an arbitrary point {α R , α I } from the dashed region of (D.12), and calculate the probability for the extensive sum S to take the corresponding value.In contrast to the real-valued case from Sec. 3.5 when the only significant contribution to f S (α < α 0 ) came from rare terms of the order of N −α , in the complex-valued situation, there is another worth-mentioning possibility to get the extensive sum's value of the order of N −α R + iN −α I : namely, apart from having one very rare term of this order, one may consider having two less rare terms of the orders N −α R + iN − αI and N − αR + iN −α I , with αR,I > α R,I .As a result, the correct expression describing the tails of the SFD f The first option under the outer MAX operation in the r.h.s represents the direct contribution from the point {α R , α I }, while the second option represents the aforementioned two-point contribution.As the terms corresponding to the points are independent, the probability density for them to occur in the same sum is a product of their individual probability densities, hence the sum of the corresponding SFDs minus one.
As one can guess, the outer MAX operation in (D.13) may result in an additional facet on the sum's SFD compared to the SFD of the summand.An example of this behavior is shown in Here, the upper-right blue points with the heights B = 1 + β and the small dashed line inbetween them schematically represent the SFD of the summand X elevated by the value of β as we usually pictured at similar diagrams corresponding to the real-valued random variables, see Fig. 6.Continuing the analogy, we used the double-dot-dashed line to mark the points where the slice-wise unit-slope lines, having exactly one common point with f X (α R , α I ) + β in the area where f X (α R , α I ) + β ≥ 1, cross the plane f (α R , α I ) = 1.As one can see, the positions and heights of the resulting SFD's vertices then follow from this geometrical construction.
Thus, for the extensive sum of c.r.v.s with a generic PDF different from zero only in a single quarter of the complex plane, the graphical rules can be formulated as follows: 1.For each slice s = α R − α I parameterized by either α R or α I , apply the extensive summation rule according to the algorithm from Sec. 3.5.
2. Draw the curve {α R (s), α I (s)} corresponding to the typical contributions from each slice and determine the position of the typical value of the whole sum as {min s α R (s), min s α I (s)}.
3. Complete the tails of the SFD with the use of (D.13).

D.5.2 Two-quarter complex distribution.
A case when the PDF of a distribution is non-zero not in one but in two quarters of the complex plane is the most relevant for our physical application.Indeed, since the associated with the local density of states imaginary part of Green's function is always greater than or equal to zero while its real part in the bulk has an SFD-symmetric distribution, the distribution of the complex Green's function falls exactly into the category of two-quarter complex distributions.So now, without loss of generality, assume the half-plane corresponding to a non-zero PDF to be the upper half-plane.Then, as at the beginning of App.D.5.1, consider the trial SFD concentrated around a single value of {α R , α I }, meaning that all of the terms in the sum are assumed to have the real and imaginary parts of the order of N −α R 0 and N −α I 0 correspondingly.Since the imaginary part of our trial random variable is non-negative, the imaginary part of the sum is also non-negative, and the zeros suppression effect for the imaginary part takes place as usual, meaning that the SFD of the extensive sum of N β such terms cannot have finite values away from the line α I = α I 0 −β.However, since the real part of the individual terms is assumed to be SFD-symmetric, the real part of the sum does not show the same zeros suppression.To see this, it is enough to separate the sum's terms into two partial sums consisting of the terms with positive and negative real parts, sum up each of the partial sums separately, and subtract one from another according to (C.4).As a result, instead of (D.9), we arrive at The case of the distributions supported on the whole complex plane can be considered similar to the ones described above.However, since our applications do not touch this topic, we will leave this generalization as an exercise for the reader.

D.6 Full Gaussian RP cavity LDOS SFD calculation
Now, after developing all the necessary techniques, let us try to self-consistently calculate a full SFD of the Gaussian RP Green's function.As we usually did in the cases of one-sided real LDOS distributions, let us start from an educated guess.In this case, let us guess the cavity self-energy Σ i = j̸ =i |V i j | 2 G j j , with Σ t y p ∝ N −c R + iN −c I .What we know about this quantity is that, due to the absence of the tails in its distribution and the zeros suppression effect, its imaginary part has an SFD concentrated around N −c I (let us pretend we do not know that c I = γ − 1).We can also assume that, since Re Σ i in the middle of the spectrum definitely has a non-zero finite PDF at Re Σ i = 0, its marginal SFD is of the Gaussian type, like (23) or (C.6).These observations, after considering the rule depicted in Fig. 14, result in the guess Σ : The rest of the calculations can be summarized as follows.After adding the independent and normally distributed on-site disorder to the above self-energy according to App.D.4, we get Then, after inverting this result according to App.D.4, we get In this parameterization, it becomes clear that the typical contributions from the whole edge to the extensive sum from (D.22) have the same real part, α R = γ − 1.
From this calculation, we see that the typical values of Re Σ i and Im Σ i have the same scaling and the result (D.20) gives the same LDOS as we obtained earlier in Sec. 4 neglecting Im Σ i .But probably the most useful observation here is the fact that the interior of the colored area from (D.21) had no impact on the distribution of the self-energy.This is, of course, just a trivial consequence of the complex version of the zeros suppression effect, but, still, it will have a significant consequence for some of our considerations in App.E. In addition, one can see a curious fact that the marginal SFD of the real part of the RP Green's function (D.20) is not trivial and reminds the SFD we obtained for the Lévy-RP in Sec. 5.

E Full Lévy-RP cavity LDOS SFD calculation
Finally, let us self-consistently calculate a full SFD of the Lévy-RP Green's function and, hopefully, obtain the desired local density of states in agreement with our numerical experiments, Fig. 8.However, instead of guessing the self-energy distribution for the Lévy-RP model, let us exploit a different approach: let us start from some arbitrary SFD and iterate the diagonal cavity Green's function expression (20) until the iterations converge.A choice of the same self-energy SFD as in the Gaussian RP case looks like a good option, so let us start from (D.18).As one can see, the exterior of the colored area and the slopes of its borders are identical to those from the Gaussian RP (D.20).And, since, for µ > 1, the slope of the white arrow parameterized by α I is always smaller than µ/2, the interior of the area will not contribute to the essential parts of the subsequent calculations, similar to how it happened in the previous cases.These two facts allow us to claim that the subsequent iterations with the use of Green's function G ii from (E.5) will lead to the exactly the same result as we obtained using (D.20).The SFD of G ii |V i j | 2 will differ from the result (E.2) only by the insignificant changes in the green area, and the self-energy will then fully coincide with the result (E.3), making our solution self-consistent, and the result (E.5) the desired answer.Thus, we can stop here and write the self-consistency equations based on (E.3).From that geometrical construction, one finds that c I = (γ − 2/µ)/(2 − 2/µ), in agreement with the value following from (40).In its turn, c R = γ − 2/µ is not equal to and even smaller than c I , in contrast to the Gaussian RP case and (D. here, the colors of the lines, again, reflect the colors of the facets whose edges contributed to the result, revealing its origins.The difference with the result (38) is in the value of f ν (c I + 0).

Figure 3 :
Figure 3: Pictorial representation of a sum of two i.r.v.s with convex SFDs.

ΣFigure 6 :
Figure 6: Graphical representation of the generalized central limit theorem.The summation result is shown in red.

Figure 7 :
Figure 7: Localization phase diagram for the Lévy Rosenzweig-Porter model, according to its LDOS distribution.

Figure 8 :
Figure 8: Spectrum of fractal dimensions of the LDOS of the Lévy-RP model for γ = 1.6, µ = 1.8, γ e f f = 1.55, and several system sizes, including the thermodynamic limit theoretical prediction f t ν and extrapolation of numerical results f n ν .A thick dashed red line shows the wrong prediction f ν (c + 0) = 1 − µc caused by us neglecting the real part of the self-energy.The correct prediction, f ν (c + 0) = 2 − µγ/2, is shown by the solid red line and derived in Appendix E. For the extrapolation, we used two different techniques for 1⃝-2 ⃝ and 3 ⃝-6 ⃝: in the former case, the straight lines are plotted through the points with the same derivative q = f ′ ν (α) for q = −2, −1, 0, 2.5, 3, and the lines' intersections are highlighted by the black points, while, in the latter case, the ansatz f ν (α(N ); N ) ∼ f ν (α; ∞) + A/ ln(N ), α(N ) ∼ α + B/ ln(N ) is applied along the specified directions: at fixed f (α), A = 0, for 3 ⃝, 5 ⃝, at fixed α, B = 0 for 4 ⃝, and at a certain ratio A/B for 6 ⃝.
and γ + c.Next, we   =  1 ergodic f r a c t a l 1 2

Figure 10 :
Figure 10: A part of the phase diagram for the LN-RP LDOS below γ = 2.The question marks signify the parameter range which we have not yet described.

2 𝑎Figure 11 :
Figure 11: A full LN-RP LDOS phase diagram.The dark blue line shows where the expression γ = 2 a is applicable, and the red point marks the tricritical point.

Figure 14 :
Figure 14: Orthographic projections of the SFD from Fig. 13a.The colors of the projections indicate edges contributing to these projections.Notice that the projections would be the same also for the SFD from Fig. 13b.
Inversion-invariant shapes, origins and images are shown by the same color.
A general case of inversion, an origin and an image are shown by different colors.

Figure 15 :
Figure 15: Graphical representations of the behavior of the SFD of a c.r.v.under inversion.The points connected by the double-sided arrows are the source and the image of each other.The red line shows a set of points invariant under the inversion operation.Therefore this red line separates each of the dashed double-sided arrows into two equal parts.The white and grey sectors of the complex plane are the source and the image of each other.
∆α R,I = α R,I 2 − α R,I 1 , ξ = ∆α R /(∆α R + ∆α I ) is a constant consistent with the zero slope between the encircled 'ones' at the points {α R i , α I i }, and, for brevity, we assume all colored arrows of the same color within the same diagram to have the same labels.This latter convention appears to be very useful on dense diagrams with many different facets.A detailed diagram corresponding to the same sum is given by (D.15) the same partitioning of the sum, instead of (D.10), we obtain another diagram with the same position of the peak but with the additional linear part ∝ −α R to the right of it.Finally, for the extensive sum of i.i.d.c.r.v.s having a generic two-quarter distribution, after the same partitioning and summing up the partial sums according to App.D.5.1, the result can be obtained as a mix of all contributions from all possible pairs of points corresponding to the two partial sums, meaning that the only modification we need to introduce to the rule from App. D.5.1 is the addition of the linear parts ∝ −α R to the right of each point.As an example, the extensive sum from (D.14) but with an SFD-symmetrized17 real part would look like N β i=1 = (D.17)Here we use the same notion of ξ from (D.14).

/ 2 1 − 2 1
!(D.20) one can see how the fracture of the edge from (D.19) happens on the line α R = α I , in agreement with the rules pictured in Fig. 15.Also, in (D.20), we label the edge from {0, c I } to {−c I , −c I } using the variable α I to simplify taking the imaginary part of the Green's function, as, in this parameterization, it is immediately clear that, according to App.D.1, the LDOS SFD calculated as Im G ii with G ii from (D.20) is equal to (27), as it should.Then, after multiplying (D.20) by |V i j | 2 from (29) according to App.D.3, we obtainG ii |V i j | 2 :where we changed the parameterization of the edge with a red arrow from α I in (D.20) to α R in (D.21) to simplify the drawing of the double-dot-dashed contour line in the final expression for the self-energy SFDΣ i = j̸ =i G j j |V i j | 2 : Then, we naturally arrive at (D.20) and, after multiplying it with |V i j | 2 defined by |V i j | 2 : used different colors for different facets to specify where these facets originated from with respect to (E.1).From what we got, one can already see that our initial guess for the self-energy was incorrect.However, this was a part of the plan, so let us continue.The second-iteration guess for the self-energy obtained from (E.2) looks like Σ :(E.3)Here, the small triangular blue facet originates from the tails of (E.1) because of (D.13), and the yellow facets appear because of the SFD-symmetric nature of the Green's function's real part, see App.D.5.2.The absence of any green facets signals that, again, as in the Gaussian RP case considered in App.D.6, the (green) interior of the angle from (E.2) formed by the black arrows does not contribute to the self-energy distribution.Continuing our iterations using (E.3) as the self-energy distribution, and adding it to the diagonal disorder defined by (23) according to App.D.2, we getϵ i + Σ i : (E.4)here, the blue and the yellow facets originate from the corresponding facets in (E.3), while the orange facet originates from the on-site energies distribution(23).Notice the black-and-white dashed line connecting the vertices {0, c I } and {0, 0} and the filling colors of the vertices.As it was already introduced in App.D.2, this notation depicts discontinuity between the blue and orange facets.Finally, by inverting (E.4) according to App.D.4, we get the complex Green's functionG ii = (ϵ i + Σ i ) −1 : (E.5) 22).Finally, let us extract the real-valued LDOS SFD from this complex Green's function distribution (E.5).It can be done using the projection rule from App. D.1 that givesIm G ii = ν i : ∝ α/2 −c I ← 2 − µγ/2 ∝ −µα/2 c I (E.6)