Anatomy of the eigenstates distribution:

a quest for a genuine multifractality

Anton Kutlin1,2$\star$ and
Ivan M. Khaymovich3,4

1 Abdus Salam International Center for Theoretical Physics, Strada Costiera 11, 34151 Trieste, Italy

2 Max Planck Institute for the Physics of Complex Systems,

Nöthnitzer Straße 38, 01187-Dresden, Germany

3 Nordita, Stockholm University and KTH Royal Institute of Technology Hannes Alfvéns väg 12, SE-106 91 Stockholm, Sweden

4 Institute for Physics of Microstructures, Russian Academy of Sciences, 603950 Nizhny Novgorod, GSP-105, Russia

${}^\star$<EMAIL_ADDRESS>
August 27, 2024

§ ABSTRACT

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 Rosenzweig-Porter (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 basis-invariant 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 characterized by a single power-law tail in PDF which provides usually the only energy/time scale in the system in addition to the standard ones: the bandwidth and the mean level spacing.
Multifractal states, in turn, have multiple (running) power-law exponents at different scales. 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 [3] and faster learning of artificial neural networks [4]. For these applications, the robust multifractal phases of matter are of crucial importance.

As another motivation, we consider a series of recent works [5, 6, 7, 8], 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) [9, 10], possessing a hierarchical structure, similar to the many-body Hilbert space.
It was believed to host multifractal states [11, 12, 13, 14, 15, 16, 17, 18, 19], even before similar claims about the actual many-body systems have been made [10]. To simplify the problem even more, a random-matrix proxy to the RRG was proposed, namely the log-normal Rosenzweig-Porter model (LN-RP) [16, 20]. 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 [21, 22, 23, 24, 25, 26].

Nevertheless, the very existence of the multifractal phase in the log-normal Rosenzweig-Porter model (unlike the fractal one in the Gaussian RP [27]) 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 [28, 29]. Given that the Gaussian Rosenzweig-Porter model [30, 31, 32, 33, 34, 35], eight years after
the discovery of the fractal phase there [36], is still the only analytically tractable model [27, 28, 37, 38] hosting an entire phase of fractal states, it makes sense to look for an analytical approach to the other RP models, like Lévy- [28, 29] or LN-RP [16, 20, 39], as such attempt doesn't look hopeless. In this paper, we show that it is indeed possible.

The paper is organized as follows.
In Sec. <ref>, 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 <ref> 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. <ref>, we demonstrate the application of the developed machinery to the Gaussian Rosenzweig-Porter model and compare our result to the previously known ones.
Section <ref> 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. <ref>, we do the same for the log-normal RP model and, again, claim the absence of LDOS multifractality.
In Sec. <ref>, 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. <ref>, 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 $|\psi(i)|^2$ at site $i$), the spectrum of fractal dimensions (SFD) $f_X(\alpha;N)$ parameterizes
the probability density function $p_X(x)$ as

\begin{equation}\label{eq:formal_def_of_sfd}
p_X(x)\rmd x = p_X(N^{-\alpha}) \ln(N) N^{-\alpha} \rmd \alpha = \ln(N) N^{f_X(\alpha;N)-1}\rmd\alpha.
\end{equation}

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(\alpha) = \lim_{N\rightarrow\infty} f_X(\alpha; 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(\alpha;N)}$ samples of the order of $X\sim N^{-\alpha}$.
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 $\tau_q=(q-1)D_q$ [2].

In the large-$N$ limit, $\tau_q$ is given by the Legendre transform of $f(\alpha)$ in the saddle-point approximation:

\begin{equation}\label{eq:tau_q}
\tau_q \stackrel{\text{def}}{=}
-\log_N \IPR_q %\sum_i |\psi(i)|^{2q}
\simeq
-\log_N \lrb{N\mean{|\psi(i)|^{2q}}}
\sim -\log_N \left(
\int \rmd\alpha N^{-\alpha q} N^{f(\alpha)}
\right) \xrightarrow[N\rightarrow\infty]{} q\alpha_q - f(\alpha_q) \ .
\end{equation}

Here the first equality is the definition of $\tau_q$ via the logarithm $\log_N$ of base $N$ of the inverse participation ratio $\IPR_q %\stackrel{\text{def}}{=}
\sum_i |\psi(i)|^{2q}$,
with the sum approximated by the averaging over the probability distribution (<ref>), noted by $\mean{\ldots}$.
In the r.h.s. of (<ref>), the Legendre transform is given by the parameter $\alpha_q$, minimizing $q\alpha-f(\alpha)$.
For a smooth $f(\alpha)$, $\alpha_q$ is defined as the solution of the equation $f'(\alpha_q)=q$ [If the relations provides several solutions for $\alpha_q$, one should pick the one maximizing $f(\alpha_q)-\alpha_q$.]. A pictorial representation of this minimization is shown in Fig. <ref>.

Graphical representation of the relation between the spectrum of fractal dimensions $f(\alpha)$ and the critical exponents $\tau_q$, Eq. (<ref>). 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 (<ref>). The label “$D_1$" corresponds to the value of $\alpha = \alpha_{q\to 1}$, responsible for the normalization $N\mean{N^{-\alpha}} \propto N^0$ and, hence, represents the actual, correctly defined, support set dimension [40], see the main text for more details.

From this result, one can readily identify at least two important values of $\alpha$: $\alpha_0$ and $\alpha_1$. The value $\alpha_0$ corresponds to the maximum of the SFD, $f(\alpha_0)=\max_\alpha\{f(\alpha)\}=1$ [The maximum value of $f(\alpha)$ is always $1$ due to the normalization condition of the probability density function (<ref>).].
Using the intuitive meaning of SFD mentioned above, one can deduce that, in the sufficiently large sample, the realizations with $x\sim N^{-\alpha_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 $\alpha_1$. For example, the SFDs of normalized wave functions with $1=\sum_{i=1}^N|\psi(i)|^2=N\mean{|\psi(i)|^2}$ always lie below the line $f(\alpha)=\alpha$ 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$ [40]. To see this, let's calculate $D_1$ as a limit of $D_q$ for $q\rightarrow 1$ via the L'Hôpital's rule:

\begin{equation}\label{eq:D_1}
D_1 = \lim_{q\rightarrow 1} D_q = - \sum_i |\psi(i)|^{2}\log_N |\psi(i)|^{2} \sim \int \alpha N^{f(\alpha;N)-\alpha} \ln N \rmd\alpha\xrightarrow[N\rightarrow\infty]{}\alpha_1.
\end{equation}

Thus, since $\alpha_1$ is responsible for the normalization, $D_1=\alpha_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 $\alpha$, corresponding to one of such discontinuities, contributes to $D_q$ 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(\alpha)$ does have finite values at $\alpha<\alpha_{1}$ [Here and further we define $\alpha_1$ as the leftmost of all equivalent ones if $f(\alpha)$ has a straight segment of a finite length with the slope $1$.].
Otherwise, we will talk about fractality.

§ 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. <ref>), a sum of two independent random variables (i.r.v.s) (Sec. <ref>), an
“ensemble mixture”
of several independent random variables with different distributions (Sec. <ref>), and a product of two i.r.v.s (Sec. <ref>). 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 doesn't 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}(\alpha)$ of $X^n$ in terms of the spectrum of fractal dimensions $f_X(\alpha)$ of $X$.
Since, by definition $p_X(N^{-\alpha}) \propto N^{f_X(\alpha)+\alpha-1}$
and, by variables substitute, $p_{X^n}(x) = n^{-1} x^{1/n-1}p_X(x^{1/n})$, we get that

\begin{equation}
p_{X^n}(N^{-\alpha}) = N^{f_{X^n}(\alpha)+\alpha-1} \propto N^{(1-1/n)\alpha} N^{f_X(\alpha/n)+\alpha/n-1},
\end{equation}

and, hence,

\begin{equation}
f_{X^n}(\alpha) = f_X(\alpha/n).
\end{equation}

Pictorially, this is just an $n$-times stretching of the spectrum of fractal dimensions in a horizontal direction with a fixed point $\alpha=0$, see Fig. <ref>. For the negative $n$, this is also a reflection with respect to the vertical axis.

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 $\alpha$ and the vertical axis to represent $f(\alpha)$. The dashed horizontal line marks the level $f(\alpha)=1$, the solid horizontal line corresponds to $f(\alpha)=0$, and the labels “$\alpha_0$" and “$n\alpha_0$" correspond to the typical values of $\alpha$ before and after the exponentiation.

§.§ Sum of two independent random variables

The main point for the sum rule, written in terms of SFDs, is that

the expression $x+y=N^{-\alpha_X} + N^{-\alpha_Y}$ is binary: it equals $N^{-\alpha_X}$ if $\alpha_Y>\alpha_X$ and $N^{-\alpha_Y}$ otherwise. In other words, as soon as $x<y$, we can always neglect $x$ completely, and vice versa. If $x+y=N^{-\alpha}$, then $\alpha=\min\{\alpha_X,\alpha_Y\}$. It allows us to write that

\begin{equation}
N^{f_{X+Y}(\alpha)} \propto N^{f_X(\alpha)}P(\alpha_X<\alpha_Y|\alpha_X=\alpha) + N^{f_Y(\alpha)}P(\alpha_Y<\alpha_X|\alpha_Y=\alpha).
\end{equation}

Without loss of generality, let's assume that $\alpha_0(X)<\alpha_0(Y)$, i.e. that typically $Y$ is negligible. Then there are three cases: $\alpha<\alpha_0(X)$, $\alpha_0(X) < \alpha < \alpha_0(Y)$, and $\alpha>\alpha_0(Y)$. Let's consider them one by one.

The case $\alpha<\alpha_0(X)$ is simple: since $N^{-\alpha}$ is bigger than the typical values of both $X$ and $Y$,

\begin{equation}\label{eq:sum_of_tailes}
N^{f_{X+Y}(\alpha)} \propto N^{f_X(\alpha)}\cdot 1 + N^{f_Y(\alpha)}\cdot 1 \implies f_{X+Y}(\alpha) = max\{f_X(\alpha),f_Y(\alpha)\};
\end{equation}

here, both $P(\alpha_X<\alpha_Y|\alpha_X=\alpha)$ and $P(\alpha_Y<\alpha_X|\alpha_Y=\alpha)$ are of the order of unity as the typical values $\alpha_0(X)$, $\alpha_0(Y)$ are in the interval $\alpha_{X,Y}>\alpha$.

The case $\alpha>\alpha_0(Y)>\alpha_0(X)$ is not much different: since now the typical values of $X$ and $Y$ are both bigger than $N^{-\alpha}$,

\begin{equation}\label{eq:sum_of_tailes-1}
N^{f_{X+Y}(\alpha)} \propto N^{f_X(\alpha)}N^{f_Y(\alpha)-1} + N^{f_Y(\alpha)}N^{f_X(\alpha)-1} \implies f_{X+Y}(\alpha) = f_X(\alpha)+f_Y(\alpha)-1;
\end{equation}

here, the corresponding conditional probabilities are denoted[It is the consequence of the assumed convexity of the functions $f_{X,Y}(\alpha)$. Otherwise, we would have to write the probabilities as $\max_{\omega\geq\alpha}N^{f_{X,Y}(\omega)-1}$.] by $N^{f_{X,Y}(\alpha)-1}$.
And, as both $f_X(\alpha)$ and $f_Y(\alpha)$ 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, $\alpha_0(X) < \alpha < \alpha_0(Y)$, the probability $P(\alpha_X<\alpha_Y|\alpha_X=\alpha)$ is of the order of one, while the probability $P(\alpha_Y<\alpha_X|\alpha_Y=\alpha)$ is of the order of $N^{f_X(\alpha)-1}$, hence,

\begin{equation}
N^{f_{X+Y}(\alpha)} \propto N^{f_X(\alpha)} + N^{f_Y(\alpha)}N^{f_X(\alpha)-1} \implies f_{X+Y}(\alpha) = f_X(\alpha).
\end{equation}

In the last equality we used the fact that $f_Y(\alpha)-1<0$ at $\alpha < \alpha_0(Y)$.
A pictorial representation of this result is shown in Fig. <ref>.

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

§.§ 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, $\sum_j w_j=1$, of other probability density functions:

\begin{equation}
p_X(x) = \sum_i w_i p_{X_i}(x).
\end{equation}

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 \rightarrow p(i)di$ and $p_{X_i}(x) \rightarrow p(x|i)$, we arrive at the probably more familiar expression

\begin{equation}
p_X(x) = \int p(x|i)p(i)\rmd i.
\end{equation}

Now, the transformation of the PDF to the SFD in the saddle-point approximation gives the relation

\begin{equation}
f_X(\alpha) = \max\limits_i\{f_i(\alpha)+f(i)-1\},
\end{equation}

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 (<ref>) for two discrete $i=X,Y$ values with the same weights.
It was a particular case of this more general mix rule [The name `mix rule' highlights the fact that $X$ can be seen as a random variable realizing in itself different ensembles with corresponding probabilities, i.e., a mixture of ensembles.].

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^{\omega}$. 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.

This 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.

§.§ Product of two independent random variables

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.

To calculate the SFD of the product of two i.r.v.s $X\sim N^{-\xi}$ and $Y\sim N^{-\eta}$, let's consider the corresponding convolution in terms of $\xi$ and $\eta$:

\begin{equation}
N^{f_{XY}(\alpha)-1} \propto \int_{-\infty}^{\infty} \rmd \xi \rmd \eta \d(\alpha - \xi - \eta) N^{f_X(\xi)-1} N^{f_Y(\eta)-1} = \int_{-\infty}^{\infty} \rmd \xi N^{f_X(\xi) + f_Y(\alpha - \xi) - 2}.
\end{equation}

Calculating this integral in the saddle-point approximation, we arrive at the expression suspiciously similar to the one from the mix rule:

\begin{equation}
f_{XY}(\alpha) = \max\limits_\xi\{f_X(\xi) + f_Y(\alpha - \xi) - 1\}.
\end{equation}

And this is not a surprise: the product can be understood purely in terms of the previously introduced mix rule (Sec. <ref>). 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^{-\alpha}$ can be now considered as a superposition of different horizontal shifts $\alpha$ (multiplication by a constant) with different vertical shifts $f(\alpha)-1$ (the constant's weights), restoring the same mix rule, see Fig. <ref>.

§.§ Extensive sums of i.i.d random variables

In Sec. <ref>, 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 doesn't 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's consider a non-negative random variable $X$ defined by its spectrum of fractal dimensions $f_X(\alpha)$ and a random variable $S$ defined as a sum of $N^\beta$ independent copies of the random variable $X$:

\begin{equation}
s = \sum_{i=1}^{N^\beta} x_i.
\end{equation}

Our task now is to calculate $f_S(\alpha)$.

We start by focusing on such values of $x=N^{-\alpha}$ that appear in any typical sample of size $N^\beta$.
Specifically, let's consider $\alpha$ such that its typical number of realizations $N^\beta p_{\alpha(X)}(\alpha) \propto N^{f_X(\alpha)-1+\beta}$ is large enough, i.e., $\alpha \in \Omega = \{\alpha|f_X(\alpha)-1+\beta\geq0\}$.

Choosing $\alpha_i\in\Omega$, one can estimate a probability to find exactly $N^{g_i-1+\beta}\delta\alpha$ terms of the order of $N^{-\alpha_i}$ in a particular realization of the sum as $\exp\{-\left(N^{g_i}-N^{f_X(\alpha_i)}\right)^2/2N^{f_X(\alpha_i)+1-\beta}\delta\alpha\}$ which is double-exponentially small in $g_i$, provided $g_i \neq f_X(\alpha_i)$.
Thus, neglecting the double-exponentially small contributions, one can write a typical value of the sum as [20, 41]

\begin{equation}\label{eq:typ_extensive_sum}
s_{typ} \propto \int_\Omega \rmd \alpha N^{-\alpha} N^{f_X(\alpha)-1+\beta} \propto N^{-\min_\Omega\{\alpha-f_X(\alpha)+1-\beta\}},
\end{equation}

meaning that $\alpha_0(S) = \min_\Omega\{\alpha-f_X(\alpha)+1-\beta\}$. In addition, one can deduce that $f_S(\alpha>\alpha_0(S))=-\infty$, which is just a consequence of the `zeros-suppression' effect mentioned in Sec. <ref> in its extreme case.

At the same time, the values of $\alpha$ corresponding to $f_X(\alpha)+\beta<1$ with $\alpha<\alpha_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. <ref>): a probability density for such an event to contribute equals to $N^\beta p_{\alpha(X)}(\alpha)\sim N^{f_X(\alpha_i)+\beta-1}\ll 1$, resulting in the following graphical rules for obtaining the SFD of an extensive sum:

* Draw $f_X(\alpha) + \beta$.

* Draw a unit-slope line having exactly one common point with $f_X(\alpha) + \beta$ in the area where $f_X(\alpha) + \beta \geq 1$.

* A point $\alpha_0$ where this unit-slope line crosses the horizontal level $f=1$ is our new typical value, according to (<ref>). There is nothing to the right of this point due to the zeros-suppression effect.

* To the left of this point $f_S(\alpha)$ just equals $f_X(\alpha)+\beta$.

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

As one can see from this construction, all tails with a slope greater than unity eventually die as $\beta\rightarrow\infty$, 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 $\propto s^{-2}$, which is exactly what they do,
see Sec. <ref> for more details. For such distributions, the unit-slope line touches $f_X(\alpha)+\beta$ exactly at $f=1$ for any $\beta>0$, thus, $f_S(\alpha)$ 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 <ref>.

§ GAUSSIAN ROSENZWEIG-PORTER MODEL

The Gaussian Rosenzweig-Porter ensemble is essentially just an ensemble of $N\times N$ Gaussian orthogonal (or unitary) random matrices with broken rotational symmetry:

\begin{equation}\label{eq:ham_RP}
H_{RP} = H_0 + V, \quad [H_{0}]_{ij} = \e_i\d_{ij},\,\e_i\in\mathcal{N}(0,1), \quad V = N^{-\gamma/2}H_{GOE/GUE},
\end{equation}

where the elements of $H_{GOE/GUE}$ 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 [42, 37]. The idea of the cavity method is to use an exact expression relating diagonal elements of an $N\times N$ Green's function $G(z)=(z-H)^{-1}$ and an $(N-1)\times(N-1)$ reduced Green's function $G^{(i)}(z)=(z-H^{(i)})^{-1}$, where $H$ and $H^{(i)}$ differ by a single site $i$ ($H^{(i)}$ has a small “cavity"):

\begin{equation}\label{eq:exact_cavity}
G_{ii}(z) = \left(z-\e_i-\sum_{j,k\neq i} V_{ij}G_{jk}^{(i)}(z)V_{ki}\right)^{-1};
\end{equation}

here, $z=\e-\rmi\eta$ is a complex-valued parameter with a small imaginary part $\eta$ to ensure the existence of $G(z)$. For $\eta>0$, one can define a local density of states $\nu_i(\e)$ as

\begin{eqnarray}\label{eq:enrg_avg}
\nu_i(\e)=\frac{1}{\pi}\Im G_{ii}(z) = \sum_{n=1}^N \frac{\eta/\pi}{(\e - E_n)^2 + \eta^2} |\psi_n(i)|^2.
\end{eqnarray}

Following [42, 43], we choose $\eta=N^\beta\delta_{\e}$ with $0<\beta<D_1$ and $\delta_{\e}$ 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 $\delta_{\e}\propto N^{-1}$.

The idea of the diagonal approximation is to say that, in a thermodynamic limit, the sum in the denominator of (<ref>) is dominated by the diagonal elements of the reduced Green's function [42, 44]:

\begin{equation}\label{eq:diagonal_cavity}
G_{ii}(\e)\xrightarrow[N\rightarrow\infty]{} \left(z-\e_i -
\sum_{j \neq i} |V_{ij}|^2 G^{(i)}_{jj}(z)
\right)^{-1}.
\end{equation}

The approximation is considered valid, providing the system is in a non-ergodic phase [See Appendix <ref> for more details on the diagonal cavity method applicability conditions.].
However, 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 <ref>), for simplicity, here we employ the latter one:

\begin{equation}\label{eq:rp_ldos_self-consistent}
\nu_i(\e) \sim \frac{\Gamma_i(\e)}{(\e-\e_i)^2 + \Gamma_i(\e)^2},
\quad
\Gamma_i(\e) = \sum_{j\neq i} |V_{ij}|^2\nu_j(\e).
\end{equation}

In this expression, we neglected the real part of the diagonal self-energy $\Sigma_i = \sum_{j\neq i} |V_{ij}|^2 G^{(i)}_{jj}(\e)$ compared to the on-site disorder amplitude [The real part of the self-energy appears in the expression for LDOS as an energy shift renormalizing the on-site energies. Hence, the shift is not important until the renormalization affects the model's bandwidth in the ergodic phase, i.e. until $\gamma\leq 1$. This can be seen just by comparing the bandwidth $N^{(1-\gamma)/2}$ of the hopping matrix $V$ to the on-site disorder amplitude, which is 1.] and omitted $\eta$ compared to the broadening $\Gamma_i(\e)$. 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 $\nu_i(\e)$ appears as a self-consistent problem [29, 42, 44, 41], and below, we show how to solve it using the SFD graphical algebra introduced earlier.

First, let's consider the broadening $\Gamma_i(\e)$.

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 its typical value, i.e., $f_{\Gamma_i}(\alpha>\alpha_0(\Gamma_i))=-\infty$. On the other hand, the terms entering the sum do not correspond to any fat-tailed distribution: $V_{ij}$ 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 the broadening typical value. As a result, its SFD should look like

\begin{equation}\label{eq:RP_broadening_SFD}
\Gamma_i : \quad
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\clip(\g-1,-0.5) rectangle (\g+1,1.1);
\draw [line width=0.5pt,dash pattern=on 3pt off 3pt,domain=\g-1:2*\g+1] plot(\x,1);
\draw [line width=0.5pt,domain=\g-1:\g+1] plot(\x,0);

\draw [line width=1pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:1] plot(\g,\y);
\filldraw [blue] (\g,1) circle (2pt);
\begin{small}
\draw (\g,-0.3) node {$c$};
\end{small}
\end{tikzpicture}.
\end{equation}

The only unknown parameter left is the value of this typical value $\Gamma_i\sim N^{-c}$, which we parameterized by $c$.

The site index $i$ is used in the previous paragraph in $\Gamma_i$ and $\nu_i$ to specify the random quantities, corresponding to this site. Due to the statistical homogeneity, assumed by the self-consistent cavity formulation of the problem and respected by the model under consideration, further, we omit this unnecessary (double) indexing and write just $\Gamma$, or $f_\Gamma(\alpha)$, or, similarly, $p_\nu(x)$, where possible. Therefore, the index-free symbols $\Gamma$, $\nu$, 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 $\nu%_i
$, substitute it to the definition of $\Gamma%_i
$, and find $c$ self-consistently. To perform the first step, one should ensure that all terms entering the expression for $\nu%_i
$ are independent – which is not true because $\Gamma%_i
$ enters the expression two times. However, since, from the SFDs' point of view, $\Gamma%_i
$ is just a constant, we can still perform this step without introducing any further complications. The calculation is depicted below:

\begin{eqnarray}
\label{eq:standard_distribution_sfd}
|\e-\e_i| :&
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\clip(-0.5,-0.5) rectangle (2,1.1);
\draw [line width=0.5pt,dash pattern=on 3pt off 3pt,domain=-1:3] plot(\x,1);
\draw [line width=0.5pt,domain=-1:3] plot(\x,0);

\draw [line width=1pt,variable=\y,color=orange,dash pattern=on 3pt off 3pt,domain=0:1] plot(0,\y);
\draw [line width=1.5pt,variable=\x,color=orange,domain=0:3] plot(\x,1-\x);
\filldraw [orange] (0,1) circle (2pt);
\begin{small}
\draw (0,-0.3) node {$0$};
\draw (1.25,0.5) node {$\propto -\alpha$};
\end{small}
\end{tikzpicture}
\quad &\leftarrow\text{an SFD of the standard distribution};
\\
(\e-\e_i)^2 :&
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\clip(-0.5,-0.5) rectangle (3,1.1);
\draw [line width=0.5pt,dash pattern=on 3pt off 3pt,domain=-1:3] plot(\x,1);
\draw [line width=0.5pt,domain=-1:3] plot(\x,0);

\draw [line width=1pt,variable=\y,color=orange,dash pattern=on 3pt off 3pt,domain=0:1] plot(0,\y);
\draw [line width=1.5pt,variable=\x,color=orange,domain=0:3] plot(\x,1-0.5*\x);
\filldraw [orange] (0,1) circle (2pt);
\begin{small}
\draw (0,-0.3) node {$0$};
\draw (2.25,0.5) node {$\propto -\alpha/2$};
\end{small}
\end{tikzpicture}
\quad &\leftarrow\text{stretching, see Sec.~\ref{sec:exponentiation}};
\\
(\e-\e_i)^2 + \Gamma_i^2:&
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\clip(-0.5,-0.5) rectangle (3,1.1);
\draw [line width=0.5pt,dash pattern=on 3pt off 3pt,domain=-1:3] plot(\x,1);
\draw [line width=0.5pt,domain=-1:3] plot(\x,0);

\draw [line width=1pt,variable=\y,color=orange,dash pattern=on 3pt off 3pt,domain=0:1] plot(0,\y);
\draw [line width=1.5pt,variable=\x,color=orange,domain=0:2*\g-2] plot(\x,1-0.5*\x);
\draw [line width=1pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:2-\g] plot(2*\g-2,\y);
\filldraw [orange] (0,1) circle (2pt);
\filldraw [blue] (2*\g-2,2-\g) circle (2pt);
\begin{small}
\draw (0,-0.3) node {$0$};
\draw (2.25,0.5) node {$\propto -\alpha/2$};
\draw (2*\g-2,-0.3) node {$2c$};
\end{small}
\end{tikzpicture}
\quad &\leftarrow\text{sum rule, see Sec.~\ref{sec:sum_rule}};
\\
\frac{1}{(\e-\e_i)^2 + \Gamma_i^2}:&
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\clip(-3,-0.5) rectangle (0.5,1.1);
\draw [line width=0.5pt,dash pattern=on 3pt off 3pt,domain=-3:3] plot(\x,1);
\draw [line width=0.5pt,domain=-3:3] plot(\x,0);

\draw [line width=1pt,variable=\y,color=orange,dash pattern=on 3pt off 3pt,domain=0:1] plot(0,\y);
\draw [line width=1.5pt,variable=\x,color=orange,domain=2-2*\g:0] plot(\x,1+0.5*\x);
\draw [line width=1pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:2-\g] plot(2-2*\g,\y);
\filldraw [orange] (0,1) circle (2pt);
\filldraw [blue] (2-2*\g,2-\g) circle (2pt);
\begin{small}
\draw (0,-0.3) node {$0$};
\draw (-2.25,0.5) node {$\propto \alpha/2$};
\draw (2-2*\g,-0.3) node {$-2c$};
\end{small}
\end{tikzpicture}
\quad &\leftarrow\text{reflection, see Sec.~\ref{sec:exponentiation}};
\\ \label{eq:RP_SFD_shape}
\frac{\Gamma_i}{(\e-\e_i)^2 + \Gamma_i^2}:&
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\clip(-3,-0.5) rectangle (0.5,1.1);
\draw [line width=0.5pt,dash pattern=on 3pt off 3pt,domain=-3:3] plot(\x,1);
\draw [line width=0.5pt,domain=-3:3] plot(\x,0);

\draw [line width=1pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:1] plot(0,\y);
\draw [line width=1.5pt,variable=\x,color=orange,domain=2-2*\g:0] plot(\x,1+0.5*\x);
\draw [line width=1pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:2-\g] plot(2-2*\g,\y);
\filldraw [blue] (0,1) circle (2pt);
\filldraw [blue] (2-2*\g,2-\g) circle (2pt);
\begin{small}
\draw (0,-0.3) node {$c$};
\draw (-2.25,0.5) node {$\propto \alpha/2$};
\draw (2-2*\g,-0.3) node {$-c$};
\end{small}
\end{tikzpicture}
\quad &\leftarrow\text{horizontal shift, see Sec.~\ref{sec:product}}.
\end{eqnarray}

As one can see from a comparison with the result obtained in [36], the shape of the SFD we have just found is correct. Now, let's write the self-consistency equation:

\begin{eqnarray}
|V_{ij}| :&
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\clip(-0.5,-0.5) rectangle (2,1.1);
\draw [line width=0.5pt,dash pattern=on 3pt off 3pt,domain=-1:3] plot(\x,1);
\draw [line width=0.5pt,domain=-1:3] plot(\x,0);

\draw [line width=1pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:1] plot(0,\y);
\draw [line width=1.5pt,variable=\x,color=blue,domain=0:3] plot(\x,1-\x);
\filldraw [blue] (0,1) circle (2pt);
\begin{small}
\draw (0,-0.3) node {$\gamma/2$};
\draw (1.25,0.5) node {$\propto -\alpha$};
\end{small}
\end{tikzpicture}
\quad &\leftarrow\text{an SFD of the normal distribution};
\\ \label{eq:rp_hopping_squared_sfd}
|V_{ij}|^2 :&
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\clip(-0.5,-0.5) rectangle (3,1.1);
\draw [line width=0.5pt,dash pattern=on 3pt off 3pt,domain=-1:3] plot(\x,1);
\draw [line width=0.5pt,domain=-1:3] plot(\x,0);

\draw [line width=1pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:1] plot(0,\y);
\draw [line width=1.5pt,variable=\x,color=blue,domain=0:3] plot(\x,1-0.5*\x);
\filldraw [blue] (0,1) circle (2pt);
\begin{small}
\draw (0,-0.3) node {$\gamma$};
\draw (2.25,0.5) node {$\propto -\alpha/2$};
\end{small}
\end{tikzpicture}
\quad &\leftarrow\text{stretching, see Sec.~\ref{sec:exponentiation}};
\\
|V_{ij}|^2\nu_j :&
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\clip(0.5,-0.5) rectangle (2*\g,1.1);
\draw [line width=0.5pt,dash pattern=on 3pt off 3pt,domain=-3:5] plot(\x,1);
\draw [line width=0.5pt,domain=-3:5] plot(\x,0);

\draw [line width=1pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:1] plot(2*\g-1,\y);
\draw [line width=1.5pt,variable=\x,color=orange,domain=1:2*\g-1] plot(\x,{1+0.5*(\x-2*\g+1)});
\draw [line width=1pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:2-\g] plot(1,\y);
\draw [line width=1.5pt,variable=\x,color=blue,domain=2*\g-1:5] plot(\x,{1-0.5*(\x-2*\g+1)});
\filldraw [blue] (2*\g-1,1) circle (2pt);
\filldraw [blue] (1,2-\g) circle (2pt);
\begin{small}
\draw (2*\g-1,-0.3) node {$\gamma+c$};
% \draw (-2.25,0.5) node {$\propto \alpha/2$};
\draw (1,-0.3) node {$\gamma-c$};
\end{small}
\end{tikzpicture}
\quad &\leftarrow\text{product rule, see Sec.~\ref{sec:product}};
\\ \label{eq:RP_self-consistensy}
\sum_{j\neq i}|V_{ij}|^2\nu_j :&
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\clip(-0.5,-0.5) rectangle (2*\g-0.25,2.15);
\draw [line width=0.5pt,dash pattern=on 3pt off 3pt,domain=-3:5] plot(\x,1);
\draw [line width=0.5pt,domain=-3:5] plot(\x,0);

\draw [line width=0.5pt,dash pattern=on 1.5pt off 1.5pt,domain=2*\g-1.4:2*\g-0.6] plot(\x,2);

\draw [line width=0.5pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:2] plot(2*\g-1,\y);
\draw [line width=0.5pt,variable=\x,color=orange,dash pattern=on 3pt off 3pt,domain=1:2*\g-1] plot(\x,{2+0.5*(\x-2*\g+1)});
% \draw [line width=0.5pt,variable=\y,color=orange,dash pattern=on 3pt off 3pt,domain=0:3-\g] plot(1,\y);
\draw [line width=0.5pt,variable=\x,color=blue,dash pattern=on 3pt off 3pt,domain=2*\g-1:5] plot(\x,{2-0.5*(\x-2*\g+1)});
\draw [line width=0.5pt,variable=\x,color=gray,dash pattern=on 3pt off 3pt,domain=0.5:1.3] plot(\x,{3-\g+(\x-1)});
\draw [line width=1pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:1] plot(\g-1,\y);
\filldraw [blue] (\g-1,1) circle (2pt);
\filldraw [blue] (2*\g-1,2) circle (2pt);
\filldraw [blue] (1,3-\g) circle (2pt);

\begin{small}
\draw (2*\g-1,-0.3) node {$\gamma+c$};
% \draw (-2.25,0.5) node {$\propto \alpha/2$};
% \draw (1,-0.3) node {$\gamma-c$};
\draw (\g-1,-0.3) node {$c^\ast=\gamma - 1$};%2\gamma-2-c$};
\draw (2*\g-0.5,2) node {$2$};
\end{small}
\end{tikzpicture}
\quad &\leftarrow\text{extensive sum, see Sec.~\ref{sec:one-sided_clt}}.
\end{eqnarray}

In the last equation, (<ref>), we demonstrated the elevation of the highest blue point with the small dashed line. Further we will use the same notation below. Comparing (<ref>) to (<ref>), one can find $c=\gamma-1$, in complete agreement with the previously known results. The ergodic and Anderson transitions then follow from the equations $c=0$ ($\Gamma%_i
\propto N^0$, thus, all $\nu_i$ are roughly of the same amplitude) and $c=1$ ($\Gamma%_i
\propto \delta_\e$ and, thus, the normalization of $\nu%_i
$ is contributed by a finite fraction of large components), giving, as expected,

\begin{equation}\label{eq:GRP_PT}
\gamma_{ET}=1, \quad \text{and} \quad \gamma_{AT}=2.
\end{equation}

While the result (<ref>) is not new, its graphical calculation already provides some insights. For example, the orange linear region with the slope $1/2$ between $\alpha=\pm 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 [Recently, an interest to non-homogeneous (Cantor-set-like) on-site energies distributions has arisen [41, 45].]. The only contribution from the off-diagonal elements to the final result (<ref>) 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_{\e_i}(x)$ is finite and non-zero as $x\to\e$.

While due to the introduced above restriction $\Gamma%_i
\gg \eta \gg \delta_\e$, we cannot apply our method in the localized phase, $\gamma>2$, let us see how the Anderson localization formally looks like from our graphical self-consistency equation's point of view. When $\gamma$ approaches $2$ from below, the blue points at $\alpha=2\gamma-2-c$ and $\alpha=\gamma-c$ from (<ref>) approach each other until they coincide for $\gamma=\gamma_{AT}=2$. For $\gamma>2$, the same picture looks like

\begin{equation}\label{eq:localized_broadening}
\sum_{j\neq i}|V_{ij}|^2\nu_j :
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\clip(0.5,-0.5) rectangle (2*\g,2.15);
\draw [line width=0.5pt,dash pattern=on 3pt off 3pt,domain=-3:5] plot(\x,1);
\draw [line width=0.5pt,domain=-3:5] plot(\x,0);

\draw [line width=0.5pt,dash pattern=on 1.5pt off 1.5pt,domain=2*\g-1.4:2*\g-0.6] plot(\x,2);

\draw [line width=0.5pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:2] plot(2*\g-1,\y);%\gamma+c

\draw [line width=0.5pt,variable=\x,color=orange,dash pattern=on 3pt off 3pt,domain=2*\g-3:2*\g-1] plot(\x,{2+0.5*(\x-2*\g+1)});
\draw [line width=1.5pt,variable=\x,color=orange,domain=1:2*\g-3] plot(\x,{2+0.5*(\x-2*\g+1)});

\draw [line width=0.5pt,variable=\x,color=blue,dash pattern=on 3pt off 3pt,domain=2*\g-1:5] plot(\x,{2-0.5*(\x-2*\g+1)});
% \draw [line width=0.5pt,variable=\x,color=gray,dash pattern=on 3pt off 3pt,domain=0.5:1.3] plot(\x,{3-\g+(\x-1)});
\draw [line width=1pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:1] plot(2*\g-3,\y);
\draw [line width=1pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:(3-\g)] plot(1,\y);

\filldraw [blue] (2*\g-3,1) circle (2pt);
\filldraw [blue] (2*\g-1,2) circle (2pt);
\filldraw [blue] (1,3-\g) circle (2pt);
\begin{small}
\draw (2*\g-1,-0.3) node {$\gamma+c$};
% \draw (-2.25,0.5) node {$\propto \alpha/2$};
% \draw (1,-0.3) node {$\gamma-c$};
\draw (2*\g-3,-0.3) node {$\gamma+c-2$};
\draw (1,-0.3) node {$\gamma-c$};

\draw (2*\g-0.5,2) node {$2$};
\end{small}
\end{tikzpicture},
\end{equation}

which already contradicts our initial guess (<ref>) for the level broadening. Nevertheless, one can try to use (<ref>) as a new guess, which results in a self-consistency equation $c=\gamma+c-2$. The unknown $c$ drops out of this equation, leaving us with the only point this construction can hold for, $\gamma=2$. And, while this attempt doesn't lead us to a correct solution, it hints at an important conclusion about the localized phase: the level broadening for $\gamma>2$ no longer has the
form (<ref>), originated from the collective contribution of many sum's terms $|V_{ij}|^2\nu_j$. Instead, as one may
assume from the previously known results [36], 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 [46],
a similar RP ensemble with the fat-tailed Lévy-distributed amplitudes was proposed [28, 29]. 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 (<ref>), with the uncorrelated diagonal on-site energies $\e_i$, distributed according to some narrow size-independent distribution like the normal or box distribution, and with the i.i.d. off-diagonal elements $V_{ij}$ of typical amplitudes $N^{-\gamma/2}$, distributed according to the following PDF with the  parameter

\begin{equation}\label{eq:student_pdf}
p(V_{ij}) \sim \frac{2\mu N^{-\mu\gamma/2}}{|V_{ij}|^{1+\mu}}\theta\lrp{|V_{ij}| - N^{-\gamma/2}}.
\end{equation}

Such a polynomial decay of the PDF tails makes the off-diagonal elements' distribution heavy-tailed for $\mu<2$ (variance is undefined) and fat-tailed for $\mu<1$ (mean of the absolute value is undefined). The SFDs of the distributions with such polynomial tails look like

\begin{equation}
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\clip(0.5*\g-1.75,-0.5) rectangle (0.5*\g+2,1.1);
\draw [line width=0.5pt,dash pattern=on 3pt off 3pt] plot(\x,1);
\draw [line width=0.5pt] plot(\x,0);

\draw [line width=1.5pt,color=blue,domain=0.5*\g:-10] plot(\x,{1+\m*(\x-0.5*\g)});
\draw [line width=1.5pt,color=blue,domain=0.5*\g:10] plot(\x,{1+(0.5*\g-\x)});
\draw [line width=1pt,color=blue,variable=\y,dash pattern=on 3pt off 3pt,domain=0:1] plot(0.5*\g,\y);
\filldraw [blue] (\g/2,1) circle (2pt);
\begin{small}
\draw (0.5*\g,-0.3) node {$\alpha_0$};
\draw (0.5*\g+1.25,0.5) node {$\propto -\alpha$};
\draw (0.5*\g-1,0.5) node {$\propto \mu\alpha$};
\end{small}
\end{tikzpicture}.
\label{eq:stable_sfd}
\end{equation}

These heavy- and fat-tail properties are precisely what led the authors of [29] 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_\infty$ 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. <ref>. As we are mostly interested in the (multi)fractal phase, which is expected [29] to exist for $1<\mu<2$ and $2/\mu<\gamma<2$, this is exactly the parameter range we consider.
Thus, we start by defining the fractal spectrum of $|V_{ij}|^2$ and trying to guess the SFD of the broadening $\Gamma%_i
$. Rescaling (<ref>) according to the exponentiation rule from Sec. <ref> and using $\alpha_0=\gamma/2$, we get

\begin{equation}\label{eq:levy-rp_hopping_squared_sfd}
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\clip(\g-2.5,-0.5) rectangle (\g+3,1.1);
\draw [line width=0.5pt,dash pattern=on 3pt off 3pt] plot(\x,1);
\draw [line width=0.5pt] plot(\x,0);

\draw [line width=1.5pt,color=blue,domain=\g:-10] plot(\x,{1+0.5*\m*(\x-\g)});
\draw [line width=1.5pt,color=blue,domain=\g:10] plot(\x,{1-0.5*(\x-\g)});
\draw [line width=1pt,color=blue,variable=\y,dash pattern=on 3pt off 3pt,domain=0:1] plot(\g,\y);
\filldraw [blue] (\g,1) circle (2pt);
\begin{small}
\draw (\g,-0.3) node {$\gamma$};
\draw (\g+2.25,0.5) node {$\propto -\alpha/2$};
\draw (\g-1.5,0.5) node {$\propto \mu\alpha/2$};
\end{small}
\end{tikzpicture}.
\end{equation}

A product of $|V_{ij}|^2$ and the LDOS $\nu_j$ is at least as heavy-tailed as $|V_{ij}|^2$. Taking also into account that $\nu%_j
$ is normalizable, and, i.e., bounded, we conclude that the extensive sum $\Gamma_i = \sum_{j\neq i} |V_{ij}|^2 \nu_j$ must have the SFD of the form

\begin{equation}
\Gamma%_i
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\clip(\g-2.5,-0.5) rectangle (\g+1,1.1);
\draw [line width=0.5pt,dash pattern=on 3pt off 3pt] plot(\x,1);
\draw [line width=0.5pt] plot(\x,0);

\draw [line width=1.5pt,color=blue,domain=\g:-10] plot(\x,{1+0.5*\m*(\x-\g)});
% \draw [line width=1.5pt,color=blue,domain=\g:10] plot(\x,{1-0.5*(\x-\g)});
\draw [line width=1pt,color=blue,variable=\y,dash pattern=on 3pt off 3pt,domain=0:1] plot(\g,\y);
\filldraw [blue] (\g,1) circle (2pt);
\begin{small}
\draw (\g,-0.3) node {$c$};
% \draw (\g+2.25,0.5) node {$\propto -\alpha/2$};
\draw (\g-1.5,0.5) node {$\propto \mu\alpha/2$};
\end{small}
\end{tikzpicture}.
\label{eq:levy_broadening_sfd}
\end{equation}

where the right-wing tails are, as usual, eaten by the extensive number of the elements in the sum. The only unknown here, like in the Gaussian RP case from Sec. <ref>, is the scaling $c$ of the typical value of $\Gamma%_i

Next, to obtain the SFD of $\nu%_i
$, we use the mix rule from Sec. <ref>. Indeed, taking for each fixed $\Gamma%_i
$-value the corresponding conditional distribution (<ref>) and the distribution of `weights', given by (<ref>), we get

\begin{equation}\label{eq:levy-rp_ldos_sfd}
\nu%_i
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\clip(-\G-0.5,-0.5) rectangle (\G+1.75,1.1);
\draw [line width=0.5pt,dash pattern=on 3pt off 3pt] plot(\x,1);
\draw [line width=0.5pt] plot(\x,0);

\draw [line width=1pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:1] plot(\G-1,\y);
\draw [line width=1pt,variable=\x,color=orange,domain=1-\G:\G-1] plot(\x,{1+0.5*(\x-\G+1)});
\draw [line width=1pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:2-\G] plot(1-\G,\y);
\filldraw [blue] (\G-1,1) circle (2pt);
\filldraw [blue] (1-\G,2-\G) circle (2pt);
\draw [line width=1.5pt,color=blue,domain=\G-1:10] plot(\x,{1-0.5*\m*(\x-1+\G)});
\begin{small}
\draw (1-\G-0.75,0.5) node {$\propto \alpha/2$};
\draw (1-\G,-0.3) node {$-c$};
\draw (\G-1+1,0.4) node {$\leftarrow 1-\mu c$};
\draw (\G-1+1.75,-0.3) node {$\propto -\mu\alpha/2$};
\draw (\G-1,-0.3) node {$c$};
\end{small}
\end{tikzpicture},
\end{equation}

where the new `zeros' part to the right of $c$ originates from the rare realizations of $\Gamma_i$ so large that $\Im[G_{ii}]\propto\Gamma_i^{-1}$. The label with the left arrow demonstrates that $f_{\nu%_i
}(c+0)=1-\mu c$. To finish the calculations, we write the self-consistency equation:

\begin{eqnarray}
\label{eq:Levy-RP_individual_term}
|V_{ij}|^2\nu_j :&
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\clip(\g-\G-1,-0.5) rectangle (\g+\G+2,1.1);
\draw [line width=0.5pt,dash pattern=on 3pt off 3pt,domain=-3:5] plot(\x,1);
\draw [line width=0.5pt,domain=-3:5] plot(\x,0);

\draw [line width=1pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:1] plot(\G-1+\g,\y);
\draw [line width=1.5pt,variable=\x,color=orange,domain=1-\G+\g:\G-1+\g] plot(\x,{1+0.5*(\x-(\G-1+\g))});
\draw [line width=1pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:2-\G] plot(1-\G+\g,\y);
\draw [line width=1.5pt,variable=\x,color=blue,domain=\G-1+\g:5] plot(\x,{1-0.5*(\x-(\G-1+\g))});
\draw [line width=1.5pt,color=blue,domain=-10:1-\G+\g] plot(\x,{2-\G+0.5*\m*(\x-(1-\G+\g))});
\filldraw [blue] (\G-1+\g,1) circle (2pt);
\filldraw [blue] (1-\G+\g,2-\G) circle (2pt);
\begin{small}
\draw (\G-1+\g+0.25,-0.3) node {$\gamma+c$};
\draw (\G-1+\g+2.25,0.5) node {$\propto -\alpha/2$};
\draw (1-\G+\g-0.25,-0.3) node {$\gamma-c$};
\draw (1-\G+\g-1.25,0.5) node {$\propto \mu\alpha/2$};
\end{small}
\end{tikzpicture}
\quad &\leftarrow\text{product rule, see Sec.~\ref{sec:product}};
\\ \label{eq:Levy-RP_self-consistensy}
\sum_{j\neq i}|V_{ij}|^2\nu_j :&
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\clip(\g-\G-2,-0.5) rectangle (\g+\G+1,2.15);
\draw [line width=0.5pt,dash pattern=on 3pt off 3pt] plot(\x,1);
\draw [line width=0.5pt] plot(\x,0);

\draw [line width=0.5pt,dash pattern=on 1.5pt off 1.5pt,domain=\G+\g-1.6:\G+\g-0.15] plot(\x,2);

\draw [line width=0.5pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:2] plot(\G-1+\g,\y);
\draw [line width=0.5pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:3-\G] plot(1-\G+\g,\y);
\draw [line width=0.5pt,variable=\x,color=orange,dash pattern=on 3pt off 3pt,domain=1-\G+\g:\G-1+\g] plot(\x,{2+0.5*(\x-(\G-1+\g))});
\draw [line width=0.5pt,color=blue,dash pattern=on 3pt off 3pt,domain=\G-1:1-\G+\g] plot(\x,{3-\G+0.5*\m*(\x-(1-\G+\g))});
\draw [line width=0.5pt,variable=\x,color=blue,dash pattern=on 3pt off 3pt,domain=\G-1+\g:5] plot(\x,{2-0.5*(\x-(\G-1+\g))});
\draw [line width=1pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:1] plot(\G-1,\y);
\draw [line width=1.5pt,color=blue,domain=-10:\G-1] plot(\x,{3-\G+0.5*\m*(\x-(1-\G+\g))});
\filldraw [blue] (\G-1,1) circle (2pt);
\filldraw [blue] (\G-1+\g,2) circle (2pt);
\filldraw [blue] (1-\G+\g,3-\G) circle (2pt);
\begin{small}
\draw (\G-1+\g+0.25,-0.3) node {$\gamma+c$};
\draw (\G-1+\g-2.75,0.5) node {$\propto \mu\alpha/2$};
\draw (1,-0.3) node {$\gamma-c$};
\draw (\G-1,-0.3) node {$c^\ast$};

\draw (\G+\g,2) node {$2$};
\end{small}
\end{tikzpicture}
\quad &\leftarrow\text{extensive sum, see Sec.~\ref{sec:one-sided_clt}}.
\end{eqnarray}

Here, $c^\ast = \gamma-c-(1-c)\cdot 2/\mu = c$, meaning that

\begin{equation}
\gamma_{eff} = c+1 = \frac{(2+\gamma)\mu-4}{2(\mu-1)}.
\end{equation}

As one can see from (<ref>), the SFD of the LDOS in the Lévy-RP model, similarly to the Gaussian RP case, corresponds to such $f(\alpha)$ that $f(\alpha<\alpha_1)=-\infty$ and, hence, it is just fractal, not multifractal. Moreover, calculating the fractal dimensions $D_1$ and $D_\infty$ using this SFD[To obtain the correct result, one should first normalize the local density of states to one (like if it would be a wave function) and only then use (<ref>), (<ref>), or Fig. <ref>.] and the self-consistent value of $\gamma_{eff}$, we get

\begin{equation}\label{eq:Levy-RP_D}
D_{1,\infty} = 1-c = 2-\gamma_{eff} = \frac{\mu(2-\gamma)}{2(\mu-1)},
\quad
2/\mu <\gamma<2,
\quad
\ ,
\end{equation}

which coincides with the results for $D_1$ from [29], but not with the ones for $D_\infty$.
The fractal phase with $0<D_1<1$, as in the Gaussian RP case from Sec. <ref>, gives a way to the ergodic one as soon as $c = 0$, providing

\begin{equation}
\gamma_{ET} = 2/\mu \ , \quad 1<\mu<2 \ ,
\end{equation}

while the localized phase appears when $c=1$, giving

\begin{equation}
\gamma_{AT} = 2 \ , \quad 1<\mu<2 \ .
\end{equation}

As one can see from these equations, the purpose of introducing $\gamma_{eff}$ was in preserving the analogy with the Gaussian RP model, as the fractal dimension takes a universal form $D_{q>1/2}=2-\gamma_{eff}$ as well as the corresponding phase diagram (<ref>).

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 $\mu>2$ and $\mu<1$. The former case of $\mu>2$ reproduces the results of the Gaussian RP: it follows from the fact that the hopping distribution ceases to be heavy-tailed, the slope $\mu/2$ from (<ref>) becomes larger than one, and, as a consequence, the extensive sum from (<ref>) produces the same expression for $c^\ast$ as in the Gaussian RP model, giving the ergodic and the Anderson localization transitions at $\gamma_{ET}=1$ and $\gamma_{AT}=2$ for any $\mu>2$. The latter case of $\mu<1$ is a bit more interesting: similarly to the situation with $\gamma>2$ from the end of Sec. <ref>, $c$ drops out of the corresponding self-consistency equation, giving the Anderson localization transition line as $\gamma_{AT}=2/\mu$. However, the support set dimensions (<ref>) 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. <ref>.

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

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 [41]. However, an attempt to verify our analytical prediction can be seen in Fig. <ref>: an extrapolation to infinite sizes [10, 36, 47] is shown by black dots, and our theoretical prediction is by the thick red line(s). The meaning of two red lines, dashed and solid, to the right of the typical value, $\alpha>\alpha_0$, is the following: the dashed line is plotted according to (<ref>), while the solid line shows the actual behavior of the SFD in this region supported by the extrapolation of our numerical results.

Spectrum of fractal dimensions for the local density of states in the Lévy-Rosenzweig-Porter model for $\gamma=1.6$ and $\mu=1.8$. For these parameters, $\gamma_{eff}=1.55$. A thick dashed red line shows the wrong prediction $f_{\nu%_i
}(c+0)=1-\mu c$ caused by us neglecting the real part of the self-energy. The correct prediction, $f_{\nu%_i
}(c+0)=2-\mu \gamma/2$, is shown by the solid red line and derived in Appendix <ref>.

The reason why the calculations above failed to capture the behavior at $\alpha>\alpha_0$ should originate from one of the approximations we made on the way. Based on the estimation from Appendix <ref>, we can conclude that the diagonal cavity approximation holds in the parameters range $2/\mu<\gamma<2$, which we defined at the beginning of Sec. <ref>, and its violation cannot be the reason for the inconsistency between (<ref>) and simulation in Fig. <ref>. The only other approximation we made following [29] was throwing away the real part of the
self-energy $\Sigma_i(\e-\rmi\eta)=\sum_{j\neq i} |V_{ij}|^2 G^{(i)}_{jj}(\e-\rmi\eta)$ compared to the diagonal disorder. This approximation seems reasonable until the typical value of $\Re\Sigma_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, providing the conditions above hold, the SFD of $|\e_i+\Re\Sigma_i|$ will be identical to that of $|\e_i|$.
Strictly speaking, to prove the above statement, one should generalize the notion of SFD to distributions supported on the entire real axis (see Appendix <ref> for details of this approach). However, here one can provide the simpler reasoning. Indeed, first, the shape of SFD for $|\e_i+\Re\Sigma_i|$ for the atypically large values of $\Re\Sigma_i\sim N^{-\alpha}$, $\alpha<\alpha_0(|\Re\Sigma_i|)$, trivially coincides with that of $|\e_i|+|\Re\Sigma_i|$. And second, the linear decay of the SFD with the  unit slope at $\alpha>\alpha_0(|\Re\Sigma_i|)$ follows from the finite value of the PDF for $|\e_i+\Re\Sigma_i|$ at zero, guaranteed, in turn, by the mutual independence of $\e_i$ and $\Sigma_i$.

In the case of the Lévy-RP, the heavy-tailed hopping distribution violates the second condition because it is more heavy-tailed than the on-site disorder distribution. This 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-valued random variables. This broad and difficult topic is covered extensively in Appendix <ref>, but the take-home message from that generalization is encouraging: the real part of the self-energy can only contribute to the untypically small values of LDOS, $\alpha>\alpha_0$. This fact can also be understood intuitively: while we increase the denominator of (<ref>) keeping $\Gamma_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 $|\e-\e_i| \sim O(1)$. Finally, as $\Re\Sigma_i$ starts to dominate when it is larger than $1$, it only contributes to the SFD part to the right of $\alpha_0$. Hence, we can continue neglecting it, providing we are only interested in $D_q$ with $q>0$, given by $\alpha\leq \alpha_0$.

§ LOG-NORMAL ROSENZWEIG-PORTER MODEL

Another candidate to host a genuinely multifractal phase, circulating in the literature, is the so-called log-normal Rosenzweig-Porter model [16, 20, 39].
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$:

\begin{equation}
p_{V_{i\neq j}}(v) \propto \frac{1}{|v|}\exp\left\{-\frac{\ln^2(|v|/v_{typ})}{2p\ln(v_{typ}^{-1})}\right\},\quad v_{typ}=N^{-\gamma/2}.
\end{equation}

Proceeding according to (<ref>), we find that its SFD is parabolic, $f_{|V_{i\neq j}|^2}(\alpha)=1-(\alpha-\gamma)^2/(4p\gamma)$, and hence truly multifractal, see Fig. <ref>. Indeed, if we want a truly multifractal wave function, why not to start from a truly multifractal hopping distribution?

A spectrum of fractal dimensions of intensities of hopping matrix elements in the log-normal RP model.

As we usually do at the beginning of the graphical calculations, let's define the range of parameters we consider and guess the first step's SFD. Defining the hopping SFD by the relation $f_{|V_{i\neq j}|^2}(\alpha)=1-(\alpha-\gamma)^2/(4a)$, let's start from considering $1<\gamma<2$ and small $a=p\gamma>0$. The `smallness' of $a$ will be defined later. But, considering that, for $a\rightarrow 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's start from the Gaussian RP LDOS SFD (<ref>) and multiply it by the squared hopping element, with the log-normal distribution:

\begin{equation}\label{eq:ln-rp_starting_sfd}
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\clip(\g-\c-\a-1,-0.5) rectangle (\g+\c+1.5,1.1);
\draw [line width=0.5pt,dash pattern=on 3pt off 3pt] plot(\x,1);
\draw [line width=0.5pt] plot(\x,0);

\draw [line width=1.5pt,color=blue,domain=\g-\c-\a-2:\g-\c-\a,samples=50] plot(\x,{1-\c-(\x-(\g-\c))^2/(4*\a)});
\draw [line width=1.5pt,variable=\x,color=orange,domain=\g-\c-\a:\g+\c-\a] plot(\x,{1-\a/4+0.5*(\x-(\g+\c-\a))});
\draw [line width=1.5pt,color=blue,domain=\g+\c-\a:\g+\c+2,samples=50] plot(\x,{1-(\x-(\g+\c))^2/(4*\a)});

\draw [line width=1pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:1-\c-\a/4] plot(\g-\c-\a,\y);
\draw [line width=1pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:1-\a/4] plot(\g+\c-\a,\y);
\draw [line width=1pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:1] plot(\g+\c,\y);

\filldraw [blue] (\g+\c,1) circle (2pt);
\filldraw [blue] (\g+\c-\a,1-\a/4) circle (2pt);
\filldraw [blue] (\g-\c-\a,1-\c-\a/4) circle (2pt);
\begin{small}
% \draw (\g-\c-\a-0.25,0.75) node {$\propto \alpha/2$};
\draw (\g-\c-\a,-0.3) node {$A$};
\draw (\g+\c-\a-0.05,-0.3) node {$B$};
\draw (\g+\c+0.05,-0.3) node {$C$};
\end{small}
\end{tikzpicture}.
\end{equation}

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 $C=\alpha_0(|V_{ij}|^2\nu_j)$ are, correspondingly, $\gamma-c-a$, $\gamma+c-a$, and $\gamma+c$. Next, we calculate the level broadening

\begin{equation}\label{eq:ln-rp_small-a_self-consistency}
\sum_{j\neq i}|V_{ij}|^2\nu_j :
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\clip(\g-\c-\a-1.5,-0.5) rectangle (\g+\c+1.5,2.15);
\draw [line width=0.5pt,dash pattern=on 3pt off 3pt] plot(\x,1);
\draw [line width=0.5pt] plot(\x,0);

\draw [line width=0.5pt,dash pattern=on 1.5pt off 1.5pt,domain=\c+\g:\c+\g+1.15] plot(\x,2);

\draw [line width=1.5pt,color=blue,domain=\g-\c-\a-2:\c,samples=50] plot(\x,{2-\c-(\x-(\g-\c))^2/(4*\a)});
\draw [line width=0.5pt,color=blue,domain=\c:\g-\c-\a,samples=50,dash pattern=on 3pt off 3pt] plot(\x,{2-\c-(\x-(\g-\c))^2/(4*\a)});
\draw [line width=0.5pt,variable=\x,color=orange,domain=\g-\c-\a:\g+\c-\a] plot(\x,{2-\a/4+0.5*(\x-(\g+\c-\a))});
\draw [line width=0.5pt,color=blue,domain=\g+\c-\a:\g+\c+2,samples=50,dash pattern=on 3pt off 3pt] plot(\x,{2-(\x-(\g+\c))^2/(4*\a)});

\draw [line width=0.5pt,color=gray,domain=\g-\c:\c-0.5,samples=50,dash pattern=on 3pt off 3pt] plot(\x,{2-\c-((\g-\c-2*\a)-(\g-\c))^2/(4*\a)+(\x-(\g-\c-2*\a))});

\filldraw [blue] (\c+\g,2) circle (2pt);
\filldraw [blue] (\c,1) circle (2pt);
\filldraw [blue] (\g-\c-\a,2-\c-\a/4) circle (2pt);
\filldraw [gray] (\g-\c-2*\a,{2-\c-((\g-\c-2*\a)-(\g-\c))^2/(4*\a)}) circle (2pt);

\draw [line width=0.5pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:2] plot(\c+\g,\y);
\draw [line width=1.5pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:1] plot(\c,\y);
\draw [line width=0.5pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:2-\c-\a/4] plot(\g-\c-\a,\y);
\draw [line width=1pt,variable=\y,color=gray,dash pattern=on 3pt off 3pt,domain=0:{2-\c-((\g-\c-2*\a)-(\g-\c))^2/(4*\a)}] plot(\g-\c-2*\a,\y);

\begin{small}
\draw (\c+\g,-0.3) node {$\gamma+c$};
\draw (\c,-0.3) node {$c^\ast$};
\draw (\g-\c-2*\a,-0.3) node {$X$};
\draw (\g-\c-\a,-0.3) node {$A$};

\draw (\c+\g+1.3,2) node {$2$};
\end{small}
\end{tikzpicture}.
\end{equation}

The gray dot located at $X=\gamma-c-2a$ marks the point on the SFD (<ref>) where its slope is equal to one. The tangent line with a unit slope, touching this point, crosses the level $f(\alpha)=1$ at the point $c^\ast=X-(1-c-a)$. Substituting $c^\ast\rightarrow c$,
we get

\begin{equation}\label{eq:c&gamma_at_p<1/2}
c=\gamma-a-1 \text{ and }  \gamma_{eff}=\gamma-a \ .
\end{equation}

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 $\alpha<\gamma_{eff}$. At the same time, for $\alpha>\gamma_{eff}$ it should show a very low probability of zeros. This is the only place where the parabolic input reveals itself. Formally speaking, this is already multifractality for $D_q$ with negative $q$, but from the physical application perspective and according to the definition we gave at the end of Sec. <ref>, we consider it as a trivial fractal phase.

From Eq. (<ref>), one can find the part of the phase diagram, Fig. <ref>. 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

\begin{equation}\label{eq:gamma_ET_small_a}
\gamma_{ET} = a+1 \quad\Leftrightarrow \quad \gamma_{ET} = 1/(1-p), \quad p<1/2 \ .
\end{equation}

In the r.h.s. we have used the standard substitution $a = p\gamma$.

Let's now consider the question of the parameter's $a$ smallness. The geometric construction pictured in (<ref>) holds only until $X>c^\ast$. This implies the restriction $\gamma<2$. Given that, we can now draw a part of our model's phase diagram, Fig. <ref>.

A part of the phase diagram for the LN-RP LDOS. The question marks signify the parameter range which we haven't yet described.

Next, let's try to move to higher values of $\gamma$. To do that, let's start with some $\gamma$ in the fractal phase, $1+a<\gamma<2$, and gradually increase it until we reach the unexplored territory. When the gray dot from Eq. (<ref>) crosses the horizontal dashed line $f(\alpha)=1$, we can no longer determine $c^\ast$ as a crossing point between the unit-slope tangent line and the level $f(\alpha)=1$. Instead, we have to use the following construction:

\begin{equation}\label{eq:ln-rp_large-gamma_self-consistency}
\sum_{j\neq i}|V_{ij}|^2\nu_j :
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\clip(\g-\c-\a-2,-0.5) rectangle (\g+\c+0.5,2.15);
\draw [line width=0.5pt,dash pattern=on 3pt off 3pt] plot(\x,1);
\draw [line width=0.5pt] plot(\x,0);

\draw [line width=0.5pt,dash pattern=on 1.5pt off 1.5pt,domain=\c+\g-1:\c+\g+0.25] plot(\x,2);

\draw [line width=1.5pt,color=blue,domain=\g-\c-\a-2:\c,samples=50] plot(\x,{2-\c-(\x-(\g-\c))^2/(4*\a)});
\draw [line width=0.5pt,color=blue,domain=\c:\g-\c-\a,samples=50,dash pattern=on 3pt off 3pt] plot(\x,{2-\c-(\x-(\g-\c))^2/(4*\a)});
\draw [line width=0.5pt,variable=\x,color=orange,domain=\g-\c-\a:\g+\c-\a,dash pattern=on 3pt off 3pt] plot(\x,{2-\a/4+0.5*(\x-(\g+\c-\a))});
\draw [line width=0.5pt,color=blue,domain=\g+\c-\a:\g+\c+2,samples=50,dash pattern=on 3pt off 3pt] plot(\x,{2-(\x-(\g+\c))^2/(4*\a)});

\filldraw [blue] (\c+\g,2) circle (2pt);
\filldraw [blue] (\c,1) circle (2pt);

\draw [line width=0.5pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:2] plot(\c+\g,\y);
\draw [line width=1.5pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:1] plot(\c,\y);

\begin{small}
\draw (\c+\g,-0.3) node {$\gamma+c$};
\draw (\c,-0.3) node {$c^\ast$};

\draw (\c+\g-1.15,2) node {$2$};
\end{small}
\end{tikzpicture}
\end{equation}

From this we find that $c^\ast=\gamma-c-2\sqrt{a(1-c)}$, and the self-consistent solution for $c$ is now

\begin{equation}\label{eq:c_LNRP}
c = \frac{\gamma-a-\sqrt{a(4+a-2\gamma)}}{2} \ .
\end{equation}

Since the part of the SFD for the level broadening to the left of $c$ doesn't 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 $\alpha=c$.
Note that after the substitution of $a=p\gamma$ Eq. (<ref>) agrees well with (49) in [39] and (C3) in [29] for $\tau^*\equiv \alpha(\gamma,p)\equiv c$.

The square root in the self-consistent expression (<ref>) for $c$ immediately tells us that the result is valid only until $\gamma\leq 2+a/2$. From the geometrical point of view, the line $\gamma=2+a/2$ corresponds to the situation when the orange segment of (<ref>) touches the level $f(\alpha)=0$. Notice a similarity with the case $\gamma> 2$, discussed at the end of Sec. <ref> in the context of the Gaussian RP model: the change in geometry happening for $\gamma>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. <ref>, 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_{\nu%_i
}(\alpha=-1)=0$ and signaling the Anderson transition.
Thus, the expression for the Anderson transition from the fractal phase is given by

\begin{equation}\label{eq:gamma_AT_small_a}
\gamma_{AT}=2+\frac{a}2 \quad\Leftrightarrow\quad \gamma_{AT} = \frac2{1-p/2} \ ,
\end{equation}

where the latter expression is straightforwardly derived from $\gamma = 2+a/2$ after the substitution $a=p\gamma$.

Especially intriguing is that, like in the Lévy-RP model at $\gamma=2/\mu$ and $\mu<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's find the line corresponding to the ergodicity transition at $\gamma>2$. To do that, we solve the equation $c=0$ and get

\begin{equation}\label{eq:gamma_ET_large_a}
\gamma_{ET}=2\sqrt{a}
\quad\Leftrightarrow\quad
\gamma_{ET} = 4p \ ;
\end{equation}

the latter form, again, is the known expression obtained from the former one by the substitution $a=p\gamma$.
As a result, summarizing Eqs. (<ref>), (<ref>), and (<ref>), we obtain the full phase diagram for the LN-RP model, see Fig. <ref>, confirming the previous results from [16, 20, 29, 39].

A full LN-RP LDOS phase diagram. The dark blue line shows where the expression $\gamma=2\sqrt{a}$ is applicable, and the red point marks the tricritical point.

It has a tricritical point analogous to the Lévy-RP case.
According to [16], the case of RRG corresponds to the bisector $a=\gamma$ ($p=1$) in this phase diagram, which crosses the above tricritical point $a=\gamma=4$.
And again, the diagram may only host multifractal states on the border of the localized phase. The rest of the diagram contains only localized, fractal, or ergodic phases.

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, doesn't 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 [17, 18], the absence of multifractality in the LDOS SFD doesn't 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 $\psi_n(i)$ instead of the local density of states $\nu_i(\e)$. And, while we cannot calculate the corresponding SFD $f_{|\psi|^2}(\alpha)$ directly, we can infer restrictions for this function implied by the shape of $f_\nu(\alpha)$. Surprisingly, this analysis can be performed for any random Hamiltonian, and the result of this section goes beyond the Rosenzweig-Porter family.

By definition (<ref>), the local density of states $\nu_i(\e)$ is proportional to the average of the squared wave functions' amplitudes $|\psi_n(i)|^2$ over the energy window $\eta$ around the energy $\e$. This energy window for the averaging is controlled by the Lorentzian function meaning that its tails decrease as $(\e-E_n)^{-2}$.
First, for simplicity, let's consider a box kernel instead of the Lorentzian one and define the box LDOS $\tilde{\nu}_i(\e)$ as

\begin{equation}\label{eq:nu(E)_box}
% \tilde{\nu}_i(\e) = \delta_{\e}^{-1}\mean{|\psi_n(i)|^2}_{\e\pm\eta}
% = \delta_{\e}^{-1}\frac{\sum_{n\atop |E_n-\e|<\eta} |\psi_n(i)|^2}{\sum_{n\atop |E_n-\e|<\eta}1}.
\tilde{\nu}_i(\e) = \delta_{\e}^{-1}\mean{|\psi_n(i)|^2}_{\e\pm\eta}
= \delta_{\e}^{-1}\frac{\sum_{n=1}^N |\psi_n(i)|^2\theta(\eta - |E_n-\e|)}{\sum_{n=1}^N \theta(\eta - |E_n-\e|)}.
\end{equation}

Later, we return back to the standard Lorenzian LDOS.

From now on, let's fix a specific site $i$ of our system. In the RP family's models, all sites are statistically equivalent, but, in general, they don't have to be. Each realization of a random Hamiltonian $H$ will then produce a single realization of $\tilde{\nu}_i(\e)=\tilde{\nu}_i(\e;H)$ and of the order of $N^\beta=\eta/\delta_{\e}\gg 1$ realizations of $|\psi_n(i)|^2=|\psi_n(i;H)|^2$ contributing to the value of $\tilde{\nu}_i(\e;H)$, Eq. (<ref>). Our first task is to relate the distribution of $|\psi_n(i)^2|$ at $|\e-E_n|<\eta$ to the distribution of $\tilde{\nu}_i(\e)$.
For this, we consider $0<\beta\ll 1$ in order to assume that the wave functions in this small energy window are statistically equivalent. [If they are not, the result still holds but has a different physical meaning. It then relates the distribution of LDOS to the marginal distribution of all eigenstate coefficients from the considered energy window.]

The values of $|\psi_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. <ref> to this case of Eq. (<ref>). Having said that, let's consider the whole ensemble of $H$ and a set $\Tilde{\Omega}_\alpha$ consisting of all $|\psi_n(i;H)|^2$ contributing to $\tilde{\nu}_i(\e;H)$ such that $\tilde{\nu}_i(\e;H)=N^{-\alpha}$:

\begin{eqnarray}\label{eq:Omega_set}
\Tilde{\Omega}_\alpha = \left\{|\psi_n(i;H)|^2 \,\Big|\, |E_n-\e|<\eta \, ,\, \tilde{\nu}_i(\e;H)=N^{-\alpha}\right\}.
\end{eqnarray}

The unconditional distribution $p_{|\psi|^2}(x)$ of $|\psi_n(i)|^2$ from our energy window can be obtained from the conditional distribution $p_{|\psi|^2}(x|x\in\Tilde{\Omega}_\alpha)$ of $|\psi_n(i)|^2\in\Tilde{\Omega}_\alpha$ by a probability chain rule as

\begin{equation}
p_{|\psi|^2}(x) = \int p_{|\psi|^2}(x|x\in\Tilde{\Omega}_\alpha) p_{\tilde{\nu}}(N^{-\alpha}) \rmd N^{-\alpha}.
\end{equation}

Transposing from the probability density functions to the corresponding spectra of fractal dimensions, we recover the mix rule from Sec. <ref>:

\begin{equation}
f_{|\psi|^2}(\alpha) = \max\limits_\xi\left\{f_{|\psi|^2}(\alpha|N^{-\alpha}\in\Tilde{\Omega}_\xi)+f_{\tilde{\nu}}(\xi)-1\right\}.
\end{equation}

This formula directly relates the SFD of the eigenvalues $f_{|\psi|^2}(\alpha)$ to the SFD of the box LDOS $f_{\tilde{\nu}}(\xi)$.

Now, let's consider $f_{|\psi|^2}(\alpha|N^{-\alpha}\in\Tilde{\Omega}_\xi)$ and infer what it may look like. First, from the definition of $\Tilde{\Omega}_\xi$, Eq. (<ref>), we know that $|\psi|^2$ from our conditional distribution $p_{|\psi|^2}(x|x\in\Tilde{\Omega}_\xi)$ cannot exceed $N^{-\xi}\delta_{\e}=N^{-\xi-1}$, see Eq. (<ref>).

Indeed, otherwise, such $|\psi|^2$ would give rise to $\tilde\nu>N^{-\xi}$ at least for some realizations of $H$ related to $\Tilde{\Omega}_\xi$. Second, we know that the point $\alpha=\xi-\ln_N\d_\e=\xi+1$, $f_{|\psi|^2}(\alpha)=1$ belongs to $f_{|\psi|^2}(\alpha|N^{-\alpha}\in\Tilde{\Omega}_\xi)$. Otherwise, we would get $\tilde\nu<N^{-\xi}$ for some $H$ contributing to our conditional distribution. Hence, the conditional SFD is constrained to have the following form:

\begin{equation}\label{eq:conditional_SFD}
|\psi_n(i)|^2 \in \Tilde{\Omega}_\xi :
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\clip(-0.5,-0.5) rectangle (2,1.1);
\draw [line width=0.5pt,dash pattern=on 3pt off 3pt,domain=-1:3] plot(\x,1);
\draw [line width=0.5pt,domain=-1:3] plot(\x,0);

\draw [line width=1.5pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:1] plot(0,\y);
\draw [line width=0.75pt,variable=\x,color=pink,domain=0:3] plot(\x,1-\x);
\draw [line width=0.75pt,variable=\x,color=orange,domain=0:3,samples=150] plot(\x,{1-0.8*\x*sin(120*\x^2)^2});
\draw [line width=0.75pt,variable=\x,color=green,domain=0:3] plot(\x,{1-0.1*\x^2});
\filldraw [blue] (0,1) circle (2pt);
\begin{small}
\draw (0,-0.3) node {$1+\xi$};
\end{small}
\end{tikzpicture}.
\end{equation}

Here, the blue part of the graph is the necessary part: $f_{|\psi|^2}(\alpha|N^{-\alpha}\in\Tilde{\Omega}_\xi)$ has to be equal to $1$ at $\alpha=1+\xi$ and to $-\infty$ at $\alpha<1+\xi$. In the interval $\alpha>1+\xi$, our conditional SFD can have any possible shape allowed by our definition of SFD, $f_{|\psi|^2}(\alpha|N^{-\alpha}\in\Tilde{\Omega}_\xi)\leq 1$. Several examples of that are depicted in (<ref>) 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_{|\psi|^2}(\alpha)$ can differ from the corresponding SFD $f_{\tilde{\nu}}(\alpha)$ of the box LDOS, Eq. (<ref>), only in the region $\alpha>\alpha_0(\tilde{\nu})$. To the left of this point, these two SFDs must coincide [Here we have considered the convex $f_{\tilde{\nu}}(\alpha)$. For non-convex ones, the difference may appear to the right, $\alpha>\alpha_*$, of each maximum $f_{\tilde{\nu}}(\alpha_*) = f_*$, but with the deviated values $f_{|\psi|^2}(\alpha)\leq f_*$].

Finally, let's 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 $\e\pm\eta$. This possibility lifts the restriction for the point $\{1+\xi,1\}$ to belong to $f_{|\psi|^2}(\alpha|N^{-\alpha}\in\Tilde{\Omega}_\xi)$ for all $\xi$ except $\xi=\alpha_1(\nu)-1\equiv D_1(\nu)-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:

\begin{equation}
f_{|\psi|^2}(D_1) = f_{\nu}(D_1-1), \quad\text{and}\quad
f_{|\psi|^2}(\alpha) \leq f_{\nu}(\alpha-1), \quad \alpha<\alpha_0(\nu).
\end{equation}

For the log-normal and Lévy Rosenzweig-Porter models, it means that $D_q(|\psi|^2)=D_q(\nu)$ for $q\geq 1/2$. This conclusion follows from the fact that, for both of these models, $f_\nu(\alpha<D_1(\nu))=-\infty$, while the derivative of $f_\nu(\alpha)$ at $\alpha=D_1(\nu)+0$ is equal to $1/2$.
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's 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's consider a random Hamiltonian $H=H_{RP} + A_{ER}$, where $H_{RP}$ is a Gaussian RP model Hamiltonian from (<ref>), 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

\begin{equation}\label{eq:er-rp_hopping_squared_sfd}
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\clip(\g-2.5,-0.5) rectangle (\g+3,1.1);
\draw [line width=0.5pt,dash pattern=on 3pt off 3pt] plot(\x,1);
\draw [line width=0.5pt] plot(\x,0);

\draw [line width=1.5pt,color=blue,domain=\g:10] plot(\x,{1-0.5*(\x-\g)});
\draw [line width=1pt,color=blue,variable=\y,dash pattern=on 3pt off 3pt,domain=0:1] plot(\g,\y);
\filldraw [blue] (\g,1) circle (2pt);
\filldraw [blue] (0,0) circle (2pt);
\begin{small}
\draw (0,-0.3) node {$0$};
\draw (\g,-0.3) node {$\gamma$};
\draw (\g+2.25,0.5) node {$\propto -\alpha/2$};
\end{small}
\end{tikzpicture};
\end{equation}

hereafter, let's assume $\gamma>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^{-\gamma/2}$, see, e.g., [48] for the case of correlated hopping of this kind.
Performing the calculations of the SFD for $\nu%_i
$, one can straightforwardly show that for $\gamma>1+c$ any SFD curve, respecting Mirlin-Fyodorov symmetry [49], $f_\nu(\alpha)=\alpha+f_\nu(-\alpha)$, which is finite only on the support $|\alpha|<c$, i.e., $f_\nu(|\alpha|>c) =-\infty$, satisfies the self-consistent cavity equation for such hopping.
For example, one can take the parabolic one $f_\nu(\alpha)=1-(\alpha-c)^2/4c$ for $|\alpha|<c$:

\begin{equation}\label{eq:multifractal_ldos_sfd}
\nu%_i
\begin{tikzpicture}[baseline={(current bounding box.center)}]
\clip(-\c-1,-0.5) rectangle (\c+1,1.1);
\draw [line width=0.5pt,dash pattern=on 3pt off 3pt] plot(\x,1);
\draw [line width=0.5pt] plot(\x,0);

\draw [line width=1.5pt,color=blue,domain=-\c:\c] plot(\x,{1-(\x-\c)^2/(4*\c)});
\draw [line width=1pt,color=blue,variable=\y,dash pattern=on 3pt off 3pt,domain=0:1] plot(\c,\y);
\draw [line width=1pt,color=blue,variable=\y,dash pattern=on 3pt off 3pt,domain=0:1-\c] plot(-\c,\y);
\filldraw [blue] (\c,1) circle (2pt);
\filldraw [blue] (-\c,1-\c) circle (2pt);
\begin{small}
\draw (-\c,-0.3) node {$-c$};
\draw (\c,-0.3) node {$c$};
\end{small}
\end{tikzpicture}.
\end{equation}

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., [15] for the RRG. Analogously, the different lattice dimensionalities of the standard Anderson model drastically change the localization diagram [50, 51].
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's 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. <ref>: we assume the solution found, trace back its features to the input distributions, and conclude if such a self-consistent solution can actually exist or not.

So, for concreteness, let's start with the shape of the LDOS SFD given, e.g., by (<ref>). As will be shortly seen, this particular choice doesn't affect the argument. Being a part of the iteration procedure, this shape ought to originate from mixing together different Gaussian-RP-like Lorenzian shapes of LDOS (<ref>) corresponding to different fixed values of $\Gamma$: please see the orange dashed line in  (<ref>) for $\Gamma\sim N^{-\xi}$ and recall, e.g., how we obtained (<ref>). Schematically, this inheritance can be illustrated by

\begin{equation}\label{eq:multifractal_ldos_sfd_mix}
\nu%_i
\begin{tikzpicture}[baseline={(current bounding box.center)}, scale = 2]
\clip(-\c-0.25,-0.5) rectangle (5*\c+0.25,1.1);

\draw [line width=0.5pt,dash pattern=on 3pt off 3pt, domain=-\c-0.25:\c+0.25] plot(\x,1);
\draw [line width=0.5pt,domain=-\c-0.25:\c+0.25] plot(\x,0);

\draw [line width=1.5pt,color=magenta,domain=-\c:0] plot(\x,{1-(\x-\c)^2/(4*\c)});
\draw [line width=1.5pt,color=blue,domain=0:\c] plot(\x,{1-(\x-\c)^2/(4*\c)});
\draw [line width=0.5pt,color=black,variable=\y,dash pattern=on 3pt off 3pt,domain=0:1] plot(0,\y);

\draw [line width=1pt,color=blue,variable=\y,dash pattern=on 3pt off 3pt,domain=0:1] plot(\c,\y);
\draw [line width=1pt,color=magenta,variable=\y,dash pattern=on 3pt off 3pt,domain=0:1-\c] plot(-\c,\y);
\filldraw [blue] (\c,1) circle (1pt);
\filldraw [magenta] (-\c,1-\c) circle (1pt);

\draw [line width=1pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:{1-(\c/2-\c)^2/(4*\c)}] plot(\c/2,\y);
\draw [line width=1pt,variable=\x,color=orange,dash pattern=on 3pt off 3pt, domain=-\c/2:\c/2] plot(\x,{1-(\c/2-\c)^2/(4*\c)+0.5*(\x-\c/2)});
\draw [line width=1pt,variable=\y,color=magenta,dash pattern=on 3pt off 3pt,domain=0:{1-(-\c/2-\c)^2/(4*\c)}] plot(-\c/2,\y);
\filldraw [blue] (\c/2,{1-(\c/2-\c)^2/(4*\c)}) circle (1pt);
\filldraw [magenta] (-\c/2,{1-(-\c/2-\c)^2/(4*\c)}) circle (1pt);

\begin{small}
\draw (-\c,-0.15) node {$-c$};
\draw (\c,-0.15) node {$c$};
\draw (0,-0.15) node {$0$};

\draw (-\c/2,-0.15) node {$-\xi$};
\draw (\c/2,-0.15) node {$\xi$};
\end{small}

\draw [line width=0.5pt,dash pattern=on 3pt off 3pt, domain=3*\c-0.25:5*\c+0.25] plot(\x,1);
\draw [line width=0.5pt,domain=3*\c-0.25:5*\c+0.25] plot(\x,0);

\draw [line width=0.5pt,dash pattern=on 3pt off 3pt, domain=\c/2:{\c/2+4*\c}] plot(\x,{1-(\c/2-\c)^2/(4*\c)});

\draw [line width=1.5pt,color=blue,domain=4*\c:5*\c] plot(\x,{1-(\x-\c-4*\c)^2/(4*\c)});
\draw [line width=0.5pt,color=black,variable=\y,dash pattern=on 3pt off 3pt,domain=0:1] plot(4*\c,\y);
\draw [line width=1pt,color=blue,variable=\y,dash pattern=on 3pt off 3pt,domain=0:1] plot(5*\c,\y);

\filldraw [blue] (4*\c+\c/2,{1-(\c/2-\c)^2/(4*\c)}) circle (1pt);
\filldraw [blue] (5*\c,1) circle (1pt);
\draw [line width=1pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain=0:{1-(\c/2-\c)^2/(4*\c)}] plot(4*\c+\c/2,\y);

\draw [line width=0.75pt,color=green,domain=3*\c:4*\c] plot(\x,{1-(0.5*(\x-4*\c)-\c)^2/(4*\c)});
\draw [line width=0.75pt,color=gray,domain=3.5*\c:4*\c] plot(\x,{1-(2*(\x-4*\c)-\c)^2/(4*\c)});
\draw [line width=0.75pt,color=pink,domain=3*\c:4*\c,samples=150] plot(\x,{(1-(\x-4*\c-\c)^2/(4*\c))*cos(500*(\x-4*\c))^2});

\begin{small}
\draw (3*\c-0.75*\c,0.35) node {$\Gamma:$};

\draw (5*\c,-0.15) node {$c$};
\draw (4*\c,-0.15) node {$0$};

\draw (4*\c+\c/2,-0.15) node {$\xi$};
\end{small}

\end{tikzpicture};
\end{equation}

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_\nu(\alpha)$ is controlled by the same blue part of $f_\Gamma(\alpha)$, the thin differently-colored curves of the  SFD for the broadening demonstrate the unimportance of $f_\Gamma(\alpha<0)$ as it can only affect $f_\nu(\alpha>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_\Gamma(\alpha)$ 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 $\propto \alpha/2$ of the Poisson distribution we have already used to.

In its turn, the broadening SFD is obtained from the SFD of $|V_{ij}|^2\nu_j$ by the extensive summation according to Sec. <ref>. Because of the zeros suppression effect, its tail, $f_\Gamma(\alpha<\alpha_0(\Gamma))$, may only originate from the tail of $|V_{ij}|^2\nu_j$ such that, for $\alpha<c$, $f_\Gamma(\alpha) = f_{|V|^2\nu}(\alpha) + 1$. In our case, it must look like

\begin{equation}\label{eq:mf_sfd_self-energy_summand}
\begin{tikzpicture}[baseline={(current bounding box.center)}, scale = 2]
\clip(4*\c-0.5,0.5) rectangle (5*\c+0.5,1.5);

\draw [line width=0.5pt,domain=3*\c-0.25:5*\c+0.5] plot(\x,1);

\draw [line width=1.5pt,color=blue,domain=4*\c:5*\c] plot(\x,{1-(\x-\c-4*\c)^2/(4*\c)});
\draw [line width=0.5pt,color=black,variable=\y,dash pattern=on 3pt off 3pt,domain=0:2] plot(4*\c,\y);
\draw [line width=0.5pt,color=black,variable=\y,dash pattern=on 3pt off 3pt,domain=0:2] plot(5*\c,\y);
\draw [line width=1pt,color=blue,variable=\y,domain=0.95:1.05] plot(5*\c,\y);

\filldraw [blue] (4*\c+\c/2,{1-(\c/2-\c)^2/(4*\c)}) circle (1pt);
\draw [line width=1pt,variable=\y,color=blue,dash pattern=on 3pt off 3pt,domain={1-(\c/2-\c)^2/(4*\c)}:1] plot(4*\c+\c/2,\y);

\draw [line width=0.75pt,color=green,domain=3*\c:4*\c] plot(\x,{1-(0.5*(\x-4*\c)-\c)^2/(4*\c)});
\draw [line width=0.75pt,color=gray,domain=3.5*\c:4*\c] plot(\x,{1-(2*(\x-4*\c)-\c)^2/(4*\c)});
\draw [line width=0.75pt,color=pink,domain=3*\c:4*\c,samples=150] plot(\x,{(1-(\x-4*\c-\c)^2/(4*\c))*cos(500*(\x-4*\c))^2});

\draw [line width=0.75pt,color=green,domain=5*\c:6*\c] plot(\x,{1+(\x-\c-4*\c)^2/(4*\c)});
\draw [line width=0.75pt,color=gray,domain=5*\c:6*\c] plot(\x,{1+2*(\x-\c-4*\c)^2/(4*\c)});
\draw [line width=0.75pt,color=gray,domain=5*\c:6*\c] plot(\x,{1+3*(\x-\c-4*\c)^2/(4*\c)});

\begin{small}
\draw (3*\c-0.75*\c,0.35) node {$\Gamma:$};

\draw (5*\c + 0.15,1+0.15) node {$c$};
\draw (4*\c - 0.15,1+0.15) node {$0$};

\draw (4*\c+\c/2,1+0.15) node {$\xi$};
\end{small}

\end{tikzpicture}\ ;
\end{equation}

here, as usual, the solid horizontal axis marks the level $f(\alpha)=0$. What happens to the right of $\alpha=c$ is, again, unimportant, as long as $c$ stays the typical value of the extensive sum $\sum_{j\neq i} |V_{ij}|^2\nu_j$.
Then, what the hopping distribution of such a model should be?

To answer this question, we need to consider the product rule from Sec. <ref>. 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,

\begin{equation}
\alpha_q(|V|^2\nu) = \alpha_q(|V|^2) + \alpha_q(\nu),
\quad
f_{|V|^2\nu}(\alpha_q(|V|^2\nu)) =
f_{|V|^2}(\alpha_q(|V|^2)) + f_{\nu}(\alpha_q(\nu)) - 1.
\end{equation}

From (<ref>) and (<ref>), we know that, for any point from the blue region, $\alpha_q(|V|^2\nu)=\alpha_q(\nu)$, and $f_{|V|^2\nu}(\alpha_q(|V|^2\nu)) = f_{\nu}(\alpha_q(\nu)) - 1$, leaving us with the requirement for $f_{|V|^2}(\alpha)$ 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_\nu(\alpha)$ 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.[We, of course, do not consider the cases with the on-site energy distributions described by PDFs that are non-analytic on some finite interval. While such cases may lead to multifractal LDOS SFD, we don't know any physical model that such exotic on-site energies' distributions can describe.] Finally, we use the results of Sec. <ref> 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\geq 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(\alpha)$ 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(\alpha)$ we have found is beyond the point with the tangent slope $1/2$, $\alpha>\alpha_{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.
Statistically non-homogeneous distributions of the on-site disorder (like in [41]) 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 [52] and on-site disorder correlations [53, 54].

The latter, though, does not exclude non-trivial and anomalously slow dynamics in Rosenzweig-Porter models, see, e.g., [28, 39], 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 [39], 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 [39]) 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\to\infty$ (which has been done in this paper) but also may allow calculating sub-leading orders of the finite-size scaling for $N\gg 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.

§ ACKNOWLEDGEMENTS

We are grateful to V. E. Kravtsov for fruitful discussions and to him together with B. L. Altshuler and L. B. Ioffe for the works on the related topics.

Funding information
I. M. K. acknowledges the support by
the Russian Science Foundation, Grant No. 21-12-00409.

[1]
M. L. Mehta,
Random Matrices and the Statistical Theory of Energy Levels,
https://doi.org/10.1016/C2013-0-12505-6 (1967).

[2]
F. Evers and A. D. Mirlin,
Anderson transitions,
Rev. Mod. Phys 80, 1355 (2008),

[3]
V. N. Smelyanskiy, K. Kechedzhi, S. Boixo, S. V. Isakov, H. Neven and
B. Altshuler,
Nonergodic delocalized states for efficient population transfer
within a narrow band of the energy landscape,
Phys. Rev. X 10, 011017 (2020),

[4]
K. Kechedzhi, V. N. Smelyanskiy, J. R. McClean, V. S. Denchev, M. Mohseni,
S. V. Isakov, S. Boixo, B. L. Altshuler and H. Neven,
Efficient population transfer via non-ergodic extended states
in quantum spin glass (2018), 1807.04792.

[5]
J. Z. Imbrie, V. Ros and A. Scardicchio,
Local integrals of motion in many-body localized systems,
Annalen der Physik 529(7), 1600278 (2017).

[6]
G. De Tomasi, I. M. Khaymovich, F. Pollmann and S. Warzel,
Rare thermal bubbles at the many-body localization transition
from the fock space point of view,
Phys. Rev. B 104, 024202 (2021),

[7]
A. Morningstar, L. Colmenarez, V. Khemani, D. J. Luitz and D. A. Huse,
Avalanches and many-body resonances in many-body localized
Phys. Rev. B 105, 174205 (2022),

[8]
N. Macé, F. Alet and N. Laflorencie,
Multifractal scalings across the many-body localization
Phys. Rev. Lett. 123, 180601 (2019),

[9]
G. Biroli, A. C. Ribeiro-Teixeira and M. Tarzia,
Difference between level statistics, ergodicity and
localization transitions on the Bethe lattice (2012), 1211.7334.

[10]
A. De Luca, B. L. Altshuler, V. E. Kravtsov and A. Scardicchio,
Anderson localization on the bethe lattice: Nonergodicity of
extended states,
Phys. Rev. Lett. 113, 046806 (2014),

[11]
B. L. Altshuler, E. Cuevas, L. B. Ioffe and V. E. Kravtsov,
Nonergodic phases in strongly disordered random regular
Phys. Rev. Lett. 117, 156601 (2016),

[12]
B. L. Altshuler, L. B. Ioffe and V. E. Kravtsov,
Multifractal states in self-consistent theory of localization:
analytical solution,
arXiv:1610.00758  (2016).

[13]
V. E. Kravtsov, B. L. Altshuler and L. B. Ioffe,
Non-ergodic delocalized phase in Anderson model on Bethe
lattice and regular graph,
Annals of Physics 389, 148  (2018),

[14]
I. García-Mata, O. Giraud, B. Georgeot, J. Martin, R. Dubertrand and
G. Lemarié,
Scaling theory of the Anderson transition in random graphs:
Ergodicity and universality,
Phys. Rev. Lett. 118, 166801 (2017),

[15]
G. Parisi, S. Pascazio, F. Pietracaprina, V. Ros and A. Scardicchio,
Anderson transition on the Bethe lattice: an approach with
real energies,
Journal of Physics A: Mathematical and Theoretical 53(1),
014003 (2019),

[16]
V. Kravtsov, I. Khaymovich, B. Altshuler and L. Ioffe,
Localization transition on the random regular graph as an
unstable tricritical point in a log-normal rosenzweig-porter random matrix
arXiv preprint arXiv:2002.02979  (2020).

[17]
K. S. Tikhonov and A. D. Mirlin,
Fractality of wave functions on a cayley tree: Difference
between tree and locally treelike graph without boundary,
Physical Review B 94(18), 184203 (2016).

[18]
M. Sonner, K. S. Tikhonov and A. D. Mirlin,
Multifractality of wave functions on a Cayley tree: From root
to leaves,
Phys. Rev. B 96, 214204 (2017),

[19]
I. García-Mata, J. Martin, R. Dubertrand, O. Giraud, B. Georgeot and
G. Lemarié,
Two critical localization lengths in the Anderson transition
on random graphs,
Phys. Rev. Research 2, 012020 (2020),

[20]
I. M. Khaymovich, V. E. Kravtsov, B. L. Altshuler and L. B. Ioffe,
Fragile extended phases in the log-normal rosenzweig-porter
Phys. Rev. Research 2, 043346 (2020),

[21]
K. S. Tikhonov, A. D. Mirlin and M. A. Skvortsov,
Anderson localization and ergodicity on random regular graphs,
Phys. Rev. B 94, 220203(R) (2016),

[22]
G. Biroli and M. Tarzia,
Delocalization and ergodicity of the Anderson model on
Bethe lattices (2018), 1810.07545.

[23]
K. S. Tikhonov and A. D. Mirlin,
Statistics of eigenstates near the localization transition on
random regular graphs,
Phys. Rev. B 99, 024202 (2019),

[24]
K. S. Tikhonov and A. D. Mirlin,
Eigenstate correlations around many-body localization
Phys. Rev. B 103, 064204 (2021),

[25]
K. S. Tikhonov and A. D. Mirlin,
From Anderson localization on random regular graphs to
many-body localization,
Annals of Physics p. 168525 (2021),

[26]
C. Vanoni, B. L. Altshuler, V. E. Kravtsov and A. Scardicchio,
Renormalization group analysis of the Anderson model on
random regular graphs,
<http://arxiv.org/abs/2306.14965> (2023),

[27]
P. von Soosten and S. Warzel,
Non-ergodic delocalization in the Rosenzweig–Porter
Letters in Mathematical Physics pp. 1–18 (2018),

[28]
C. Monthus,
Multifractality of eigenstates in the delocalized non-ergodic
phase of some random matrix models: Wigner–weisskopf approach,
Journal of Physics A: Mathematical and Theoretical 50(29),
295101 (2017).

[29]
G. Biroli and M. Tarzia,
Lévy-rosenzweig-porter random matrix ensemble,
Phys. Rev. B 103, 104205 (2021),

[30]
N. Rosenzweig and C. E. Porter,
"Repulsion of energy levels" in complex atomic spectra,
Phys. Rev. B 120, 1698 (1960),

[31]
A. Pandey,
Brownian-motion model of discrete spectra,
Chaos Solitons Fractals 5, 1275 (1995),

[32]
E. Brézin and S. Hikami,
Correlations of nearby levels induced by a random potential,
Nucl. Phys. B 479, 697 (1996),

[33]
T. Guhr,
Transitions toward quantum chaos: With supersymmetry from
Poisson to Gauss,
Ann. Phys. 250, 145 (1996),

[34]
A. Altland, M. Janssen and B. Shapiro,
Perturbation theory for the Rosenzweig-Porter matrix
Phys. Rev. E 56, 1471 (1997),

[35]
H. Kunz and B. Shapiro,
Transition from Poisson to Gaussian unitary statistics: The
two-point correlation function,
Phys. Rev. E 58, 400 (1998),

[36]
V. Kravtsov, I. Khaymovich, E. Cuevas and M. Amini,
A random matrix model with localization and ergodic
New Journal of Physics 17(12), 122002 (2015).

[37]
E. Bogomolny and M. Sieber,
Eigenfunction distribution for the Rosenzweig-Porter
Phys. Rev. E 98, 032139 (2018),

[38]
D. Venturelli, L. F. Cugliandolo, G. Schehr and M. Tarzia,
Replica approach to the generalized Rosenzweig-Porter
SciPost Phys. 14, 110 (2023),

[39]
I. M. Khaymovich and V. E. Kravtsov,
Dynamical phases in a “multifractal” Rosenzweig-Porter
SciPost Phys. 11, 045 (2021),

[40]
A. De Luca, A. Scardicchio, V. E. Kravtsov and B. L. Altshuler,
Support set of random wave-functions on the Bethe lattice
(2014), 1401.0019.

[41]
M. Sarkar, R. Ghosh and I. Khaymovich,
Tuning the phase diagram of a rosenzweig-porter model with
fractal disorder,
arXiv preprint arXiv:2305.13370  (2023).

[42]
D. Facoetti, P. Vivo and G. Biroli,
From non-ergodic eigenvectors to local resolvent statistics and
back: A random matrix perspective,
Europhys. Lett. 115, 47003 (2016),

[43]
A. G. Kutlin and I. M. Khaymovich,
Emergent fractal phase in energy stratified random models,
SciPost Phys. 11, 101 (2021),

[44]
G. De Tomasi and I. M. Khaymovich,
Non-hermitian rosenzweig-porter random-matrix ensemble:
Obstruction to the fractal phase,
Phys. Rev. B 106, 094204 (2022),

[45]
B. Altshuler and V. Kravtsov,
Random cantor sets and mini-bands in local spectrum of quantum
Annals of Physics p. 169300 (2023).

[46]
M. Tarzia,
Many-body localization transition in hilbert space,
Phys. Rev. B 102, 014208 (2020),

[47]
P. A. Nosov, I. M. Khaymovich and V. E. Kravtsov,
Correlation-induced localization,
Physical Review B 99(10), 104203 (2019),

[48]
G. L. Celardo, R. Kaiser and F. Borgonovi,
Shielding and localization in the presence of long-range
Phys. Rev. B 94, 144206 (2016),

[49]
A. D. Mirlin, Y. V. Fyodorov, A. Mildenberger and F. Evers,
Exact relations between multifractal exponents at the anderson
Physical review letters 97(4), 046803 (2006).

[50]
P. W. Anderson,
Absence of diffusion in certain random lattices,
Phys. Rev. 109, 1492 (1958),

[51]
E. Abrahams, P. W. Anderson, D. C. Licciardello and T. V. Ramakrishnan,
Scaling Theory of Localization: Absence of Quantum Diffusion
in Two Dimensions,
Phys. Rev. Lett. 42, 673 (1979),

[52]
W. Tang and I. M. Khaymovich,
Non-ergodic delocalized phase with poisson level statistics,
Quantum 6, 733 (2022).

[53]
F. Liu, S. Ghosh and Y. D. Chong,
Localization and adiabatic pumping in a generalized
aubry-andré-harper model,
Phys. Rev. B 91, 014108 (2015),

[54]
J. Wang, X.-J. Liu, G. Xianlong and H. Hu,
Phase diagram of a non-abelian aubry-andré-harper model with
$p$-wave superfluidity,
Phys. Rev. B 93, 104504 (2016),

[55]
Z. Burda, J. Jurkiewicz, M. A. Nowak, G. Papp and I. Zahed,
Free random lévy and wigner-lévy matrices,
Physical Review E 75(5), 051126 (2007).

[56]
P. Cizeau and J.-P. Bouchaud,
Theory of lévy matrices,
Physical Review E 50(3), 1810 (1994).

[57]
F. L. Metz, I. Neri and D. Bollé,
Localization transition in symmetric random matrices,
Physical Review E 82(3), 031135 (2010).

[58]
E. Tarquini, G. Biroli and M. Tarzia,
Level statistics and localization transitions of lévy
Physical review letters 116(1), 010601 (2016).

§ EXTENSIVE SUM OF NON-NEGATIVE I.I.D. R.V.S

Consider a non-negative random variable $X$ defined by a fractal spectrum $f_X(\alpha)$ and a random variable $S$ defined as a sum of $N^\beta$ independent copies of the random variable $X$:

\begin{equation}
s = \sum_{i=1}^{N^\beta} x_i.
\end{equation}

Our goal here is to calculate $f_S(\alpha)$.

We start by introducing a size-independent coarse-graining of the horizontal axis in $\alpha$. Thus, we won't make any difference between $x_p$ and $x_q$ if both $-\log_N(x_p)$ and $-\log_N(x_q)$ lie between $\alpha_i$ and $\alpha_{i+1}=\alpha_i+\D \alpha$. After that, we take a single realization of $S$ or, equivalently, $N^\beta$ realizations of $X$, and count how many $X$'s in this particular sample correspond to the specific $\alpha_i$:

\begin{equation}
n_i = \#(x_i | \alpha_i<-\log_N(x_i)<\alpha_i+\D\alpha).
\end{equation}

Finally, we approximate the value of $S$ using the above empirical counts $n_i$ and the corresponding empirical fractal spectrum $g_i$ defined via $n_i = N^{g_i}\D\alpha$:

\begin{equation}\label{eq:sum_redefinition}
s \propto \sum_i n_i N^{-\alpha_i} \propto \sum_i N^{g_i-\alpha_i}\D\alpha \propto N^{\max\limits_i(g_i-\alpha_i)}\D\alpha.
\end{equation}

The result is closely related to Laplace's method of approximate integration: the point corresponding to the maximum of $g_i-\alpha_i$ lies close to the point where the line with unit slope touches the histogram $g_i$, see Fig. <ref>. We will use this last expression in all the following constructions as it allows a very natural approach to the probability distribution of $S$. Namely, we can define a probability that $s$ is less than $N^{-\xi}$ as a probability that all $g_i$ lie below the line $\alpha-\xi$, see Fig. <ref>. So now, we should calculate the likelihood of the empirical counts $n_i$ to be described by a vector $\bm{n}$, and then sum up the probabilities of all realizations respecting the condition above:

\begin{equation}\label{eq:exact_CLT_CDF}
\sum_{0\leq n_i<N^{\alpha_i-\xi}\D\alpha}
P\left(\bm{n} \Big| \sum_in_i = N^\beta\right).
\end{equation}

A given spectrum of fractal dimensions $f(\alpha)$ (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^{-\xi}$ follows from that the blue histogram lies below the dashed orange line with unit slope $\alpha-\xi$. The plotted configuration illustrates the case when $\alpha-\xi>f_X(\alpha)-1+\beta$ for all $\alpha>\xi$. The opposite case with $\alpha-\xi< f_X(\alpha)-1+\beta$ for all $\alpha>\xi$ 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  suppression of zeros, $f_S(\alpha>\alpha_0)=-\infty$.

In principle, this approach could work, but in practice, it probably doesn't because of the total count conservation condition $\sum_i n_i = N^\beta$. Indeed, given that the probability to have $x \propto N^{-\alpha_i}$ is $p_i \propto N^{f_X(\alpha_i)-1}\D\alpha$, the probability $P(\bm{n}|\sum_in_i = N^\beta)$ itself is easy to write down:

\begin{equation}\label{eq:fixed_Sigma_PDF}
P\left(\bm{n} \Big| \sum_in_i = N^\beta\right) = N^\beta! \prod_{i} \frac{p_i^{n_i}}{n_i!}.
\end{equation}

However, if we want to sum it over different $\bm{n}$, we would have a hard time doing that because of the factor $N^\beta!$ and inability to sum up all factors $p_i^{n_i}/n_i!$ independently. That is why, following the ideas of classical statistical mechanics, we now introduce a variable “number of particles" $\Sigma=\Sigma_i n_i$ and a narrow distribution over different values of $\Sigma$ peaked at $\Sigma=N^{\beta}$. Choosing this distribution to be the Poisson distribution

\begin{equation}\label{eq:poisson_weight}
P(\Sigma;\beta) = \frac{(N^\beta)^\Sigma \rme^{-N^\beta}}{\Sigma!},
\end{equation}

we cancel the nasty factorial and leave ourselves with the much easier expression to operate with:

\begin{equation}\label{eq:variable_Sigma_PDF}
P(\bm{n};\beta)=P(\Sigma;\beta)P\left(\bm{n} \Big| \sum_in_i = \Sigma\right) =  \frac{(N^\beta)^\Sigma \rme^{-N^\beta}}{\Sigma!}\Sigma! \prod_{i} \frac{p_i^{n_i}}{n_i!} = \rme^{-N^\beta} \prod_{i} \frac{\nu_i^{n_i}}{n_i!};
\end{equation}

here, we introduced $\nu_i=p_i N^\beta$ as the expected filling of the $i$'th bin. For sufficiently large $N^\beta$, the Poisson weight (<ref>) becomes asymptotically close to a Gaussian one with mean $N^{\beta}$ and the standard deviation $N^{\beta/2}$, meaning that the probabilities of having $\bm{n}$ such that $\Sigma_i n_i=N^{\beta'}$ with $\beta'\neq\beta$ are double-exponentially suppressed in the thermodynamic limit. Hence, being interested only in the large-$N^\beta$ scaling behavior of the distribution of $S$, we can use (<ref>) instead of (<ref>) and write $P(s<N^{-\xi})$ from (<ref>) in a large-$N^\beta$ limit as

\begin{equation}\label{eq:variable_count_CLT_CDF}
\begin{split}
P(s<N^{-\xi}) \propto \sum_{n_i<M_i(\xi)}
P(\bm{n};\beta) &= \rme^{-N^\beta} \prod_{i|M_i(\xi)>0} \sum_{n_i=0}^{M_i(\xi)}\frac{\nu_i^{n_i}}{n_i!}
\\
&= \rme^{-N^\beta} \prod_{i|M_i(\xi)>0} \rme^{\nu_i}\frac{\Gamma(1+M_i(\xi),\nu_i)}{\Gamma(1+M_i(\xi))}
\\
&= \rme^{-N^\beta \overline{p}(\xi)} \prod_{i|M_i(\xi)>0} \frac{\Gamma(1+M_i(\xi),\nu_i)}{\Gamma(1+M_i(\xi))},
\end{split}
\end{equation}

where $M_i(\xi) = N^{\alpha_i - \xi}\D\alpha$ is the largest possible count in the bin $i$, obeying the condition $s<N^{-\xi}$,
$\Gamma(1+M_i(\xi),\nu_i)$ is the upper incomplete Gamma functions, and

\begin{equation}\label{eq:prohibited_events_probabilities}
\overline{p}(\xi) = 1 - \sum_{i|M_i(\xi)>0}p_i
\end{equation}

is a probability of all "prohibited" for the given $\xi$ 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 $\alpha-\xi>f_X(\alpha)-1+\beta$ for all $\alpha>\xi$ or, in other words, when $M_i(\xi) \gg \nu_i$ for all $i|M_i(\xi)>0$, see, again, Fig. <ref>. In this limit, we can write

\begin{equation}
\begin{split}
\prod_{i|M_i(\xi)>0} \frac{\Gamma(1+M_i(\xi),\nu_i)}{\Gamma(1+M_i(\xi))}
&= \prod_{i|M_i(\xi)>0} \left(1-\frac{\gamma(1+M_i(\xi),\nu_i)}{\Gamma(1+M_i(\xi))}\right)
\\
&\sim \exp\left(-\sum_{i|M_i(\xi)>0}\frac{\gamma(1+M_i(\xi),\nu_i)}{\Gamma(1+M_i(\xi))}\right),
\end{split}
\end{equation}

where $\gamma(1+M_i(\xi),\nu_i)$ is the lower incomplete Gamma function. In the given limit, all terms in the sum are exponentially small:

\begin{equation}
\frac{\gamma(1+M_i(\xi),\nu_i)}{\Gamma(1+M_i(\xi))} \sim \frac{\rme^{M_i(\xi)-\nu_i}}{\sqrt{2\pi M_i(\xi)}}\left(\frac{\nu_i}{M_i(\xi)}\right)^{1+M_i(\xi)}.
\end{equation}

Moreover, since $\D\alpha$ doesn't scale with $N$ and the number of terms stays the same as $N\rightarrow\infty$, and because of this increasingly small $M_i(\xi)^{-M_i(\xi)}$ factor, each subsequent term is exponentially smaller than the previous one. In other words,

\begin{equation}
\sum_{i|M_i(\xi)>0}\frac{\gamma(1+M_i(\xi),\nu_i)}{\Gamma(1+M_i(\xi))} \sim \frac{\gamma(2,\nu(1))}{\Gamma(2)}\sim \frac{\nu(1)^2}{2},
\end{equation}

where $\nu(1)=N^{f_X(\alpha_i)-1+\beta}\D\alpha$ is the expected filling at the point where $M_i(\xi)=1$, or, in other words, at $\alpha_i=\xi-\log_N(\D\alpha)\rightarrow \xi$. Since in this case $f_X(\xi)-1+\beta<0$, $\nu(1)\ll 1$, and

\begin{equation}
\prod_{i|M_i(\xi)>0} \frac{\Gamma(1+M_i(\xi),\nu_i)}{\Gamma(1+M_i(\xi))}
\sim 1-N^{2(f_X(\xi)-1+\beta)}/2.
\end{equation}

Taking into account also the factor $\rme^{-N^\beta \overline{p}(\xi)}$, we get

\begin{equation}
\begin{split}
N^{f_S(\xi)-1}\propto \frac{\rmd}{\rmd x}P(s<N^{-\xi})
&\propto \frac{\rmd}{\rmd x} \rme^{-N^\beta \overline{p}(\xi)} (1-N^{2(f(\xi)-1+\beta)}/2)
\\
&\propto N^{f_X(\xi)-1+\beta}+N^{2(f(\xi)-1+\beta)} \propto N^{f_X(\xi)-1+\beta},
\end{split}
\end{equation}

and, thus, when $\alpha-\xi>f_X(\alpha)-1+\beta$ for all $\alpha>\xi$,

\begin{equation}
f_S(\alpha) = f_X(\alpha) + \beta.
\end{equation}

Now, consider the second case, i.e., when the line $\alpha-\xi$ has one or more intersections with $f_X(\alpha)-1+\beta$ in the region $\alpha>\xi$. Because in this case, there are such $i$ at $M_i(\xi)>0$ that $\nu_i \gg M_i(\xi)$, the probability $P(s<N^{-\xi})$ becomes exponentially small. As written at the end of the caption to Fig. <ref>, this can be understood in simple entropic terms. Alternatively, it can be extracted from (<ref>) and (<ref>), using the Gamma function asymptotic, similarly to how we have done it for the intersection-free case.
Moreover, as $\xi$ and $f_X(\xi)$ enter all the expressions as the exponents of $N$, this smallness is, in fact, double exponential, meaning that $P(s<N^{-\xi})\propto \exp(N^{-\xi})$. Thus, without any involved calculations, we can say that, in this case, $f_S(\alpha)\rightarrow-\infty$ as $N\rightarrow\infty$. 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 <ref>.

§ APPLICABILITY OF THE DIAGONAL CAVITY APPROXIMATION TO THE LÉVY-RP MODEL

The diagonal cavity approximation consists in substituting the exact self-energy $\sum_{j,k\neq i}V_{ij}G_{jk}^{(i)}V_{ki}$ with the diagonal cavity self-energy $\sum_{j\neq i}G_{jj}^{(i)}|V_{ij}|^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., [42, 44]. 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 $\braket{i}{n}\braket{n}{j}$ of the order of $N^{-D_1}$ is of the order of $N^{-2(1-D_1)}\ll 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's proceed to more rigorous reasoning.

To obtain the (asymptotic) equation (<ref>), all we need to do is to neglect, for some reason, the off-diagonal terms in the exact block matrix inversion formula (<ref>). This simplification was used in numerous papers, including [29, 55, 56, 42, 57, 58]. However, it isn't so simple to justify the simplification in each particular case. Following [56] and [55], one may write the analogous to (<ref>) expression for the off-diagonal Green's function elements

\begin{eqnarray}
G_{ij}=-G_{ii}\sum_{k\neq i} V_{ik} G_{kj}^{(i)}
\end{eqnarray}

and try to estimate the off-diagonal contribution based on the generalized central limit theorem. Assume the off-diagonal contribution is negligible: for $V_{ij}$ following a distribution with the PDF $p(x)$ decreasing like $\sim x^{-1-\mu}$ and a typical value equal to $N^{-\gamma/2}$, we get

\begin{eqnarray}
|G_{ii}^{(i)}|_{typ} \sim \frac{1}{1 + N^{-\gamma} N^{2/\mu} |G_{ii}^{(i)}|_{typ}},
\\
|G_{ij}^{(i)}|_{typ} \sim N^{-\gamma/2}|G_{ii}^{(i)}|_{typ}^2,
\\
\text{\ and}\quad V_{ij}G_{jj}^{(i)} + \sum_{k\neq i,j} V_{ik}G_{kj}^{(i)} \sim N^{-\gamma/2}|G_{ii}^{(i)}|_{typ} + N^{-\gamma/2}N^{1/\mu}|G_{ij}^{(i)}|_{typ}.
\end{eqnarray}

Comparing the terms in the r.h.s. of the last equation, we get that the assumption is true for $\gamma \mu > 2$ and false otherwise.[We base our conclusion on the assumption that the tail of the PDF of the off-diagonal contribution decays not slower than the diagonal contribution's one. Given the graphical algebra discussed above and the explicit exact expressions for the quantities involved, this assumption's validity is out of the question.] Still, the diagonal approximation is sometimes used even for $\gamma \mu < 2$, as, e.g., in [29]. We don't know any analytical justification for that.

§ GRAPHICAL ALGEBRA OF (SFD-)SYMMETRICALLY DISTRIBUTED R.V.S

Before trying to generalize the rules from Sec. <ref> to the case of complex random variables, let's 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 $\alpha\in\mathbb{R}$ only describes an absolute value of the random variable $|x|=N^{-\alpha}$. 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^{\zeta}(x)$, $\zeta=\pm 1$, equals zero for $\zeta x < 0$ correspondingly, and assign conventional SFDs $f^{\pm}(\alpha)$ to each of the one-sided probability density functions $p^{\pm}(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. Here, by SFD-symmetric, we mean that $f^+(\alpha)=f^-(\alpha)=f(\alpha)$ and $w^+\propto w^- \propto N^0$, which is a weaker condition than $p(x)=p(-x)$.

After restricting ourselves to this case, we only need to adjust the meaning of $f(\alpha)$, 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. <ref>, all odd moments in terms of the symmetric SFD are now ill-defined. The exponentiation rule from Sec. <ref> 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. <ref> and the ensemble-mixture rule from Sec. <ref> stay exactly the same, and only the sum rule from Sec. <ref> requires a substantial modification [The central limit theorem from Sec. <ref> also needs to be modified. Moreover, for strictly symmetric, not only SFD-symmetric, distributions, the modification is sufficiently different from the original. However, since these modifications will not be of any use for us, we leave them as exercises for the reader.]. This 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}(\alpha)$, it is enough to know the probability density of the absolute value of the sum, which can be expressed as

\begin{equation}\label{eq:SFD-symm_sum_decomposition}
p_{|X+Y|}(s) = (w_X^+ w_Y^+ + w_X^- w_Y^-)p_{|X|+|Y|}(s) + (w_X^+ w_Y^- + w_X^- w_Y^+)p_{|X|-|Y|}(s).
% p_{|X+Y|}(x) = \frac{1}{2}p_{|X|+|Y|}(x) + \frac{1}{2}p_{|X|-|Y|}(x).
\end{equation}

The SFD corresponding to $p_{|X|+|Y|}(s)$ can be calculated using the original sum rule from Sec. <ref>, the SFD corresponding to the weighted sum of different PDFs is given by the mix rule from Sec. <ref>, so the only thing missing here is a rule for the subtraction of non-negative random variables. Let's fill this gap. To do that, we write

\begin{equation}
p_{|X|-|Y|}(s) = w_{|X|-|Y|}^+ p_{|X|-|Y|}^+(s) + w_{|X|-|Y|}^- p_{|X|-|Y|}^-(s),
%= w_{|X|-|Y|}^+ p_{|X|-|Y|}^+(s) + w_{|X|-|Y|}^- p_{|Y|-|X|}^+(-s),
\end{equation}

\begin{equation}
% p_{|X|-|Y|}^+(s) = \theta(s)\int_s^\infty \rmd\chi p_{|X|}(\chi)p_{|Y|}(\chi-s),
% \quad\text{and}\quad
% p_{|Y|-|X|}^+(-s) = \theta(-s)\int_{-s}^\infty \rmd\chi p_{|Y|}(\chi)p_{|X|}(\chi+s).
p_{|X|-|Y|}^\pm(s) = \theta(\pm s)\int_{\max\{0,s\}}^\infty \rmd\chi p_{|X|}(\chi)p_{|Y|}(\chi-s).
\end{equation}

Then, we substitute $|s|=N^{-\alpha}$ and $\chi=N^{-\xi}$ and go from the PDFs to the SFDs:

\begin{equation}\label{eq:subtract_rule}
p_{|X|-|Y|}^\pm(|s|=N^{-\alpha}) \propto \int_{-\infty}^\alpha \rmd\xi N^{-\xi} N^{f_{X}(\xi)+\xi-1} N^{f_{Y}(\xi)+\xi-1} \propto N^{\max_{\xi\leq\alpha}\{f_{X}(\xi) + f_{Y}(\xi) + \xi\} - 2}.
\end{equation}

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 (<ref>) and assuming, as earlier, that $\alpha_0(X)<\alpha_0(Y)$ and, for simplicity, that both $f_X(\alpha)$ and $f_Y(\alpha)$ are convex, we get

\begin{equation}
\propto
\begin{cases}
%N^{f_X(\alpha)+f_Y(\alpha)+\alpha-2} +
N^{\max\{f_X(\alpha),f_Y(\alpha)\}+\alpha-1} & \alpha < \alpha_0(X), \alpha_0(Y)
\\
N^{\max\limits_{\xi\leq\alpha}\{f_X(\xi)+f_Y(\xi)+\xi\}-2} + N^{f_X(\alpha)+\alpha-1} & \alpha_0(X) < \alpha < \alpha_0(Y)
\\
%+ N^{f_X(\alpha)+f_Y(\alpha)+\alpha-2}
& \alpha > \alpha_0(X), \alpha_0(Y)
\end{cases}
\label{eq:symmetric_sum}
\end{equation}

using the mix rule from Sec. <ref> with equal weights $\propto 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}(\alpha)=(1-\alpha)\tilde{\theta}(\alpha)$, where $\tilde{\theta}(\alpha)$ is $-\infty$ for negative $\alpha$ and $1$ for positive ones, we get back another normally distributed one as we should:

\begin{equation}\label{eq:symmetric_gaussian_sfd}
f_{X+Y}(\alpha) = \begin{cases}
\max\limits_{\xi\leq\alpha}\{2f(\xi)+\xi\}-\alpha-1 = 1-\alpha, & \alpha\geq 0
\\
-\infty, &\alpha<0
\end{cases}.
\end{equation}

Here, we used (<ref>), given as $p_{X+Y}(|s|=N^{-\alpha})\equiv N^{f_{X+Y}(\alpha)+\alpha-1}$.

§ 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 don't 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's see what the operations are:

* Sum of independent real and complex random variables. This operation is necessary to compute the denominator of the r.h.s. of (<ref>) and add together the real on-site energy and the complex self-energy.

* 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.

* 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.

* 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.

* 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.