Dynamical phases in a"multifractal"Rosenzweig-Porter model

We consider the static and dynamic phases in a Rosenzweig-Porter (RP) random matrix ensemble with the tailed distribution of off-diagonal matrix elements of the form of the large-deviation ansatz. We present a general theory of survival probability in such a random-matrix model and show that the {\it averaged} survival probability may decay with time as the simple exponent, as the stretch-exponent and as a power-law or slower. Correspondingly, we identify the exponential, the stretch-exponential and the frozen-dynamics phases. As an example, we consider the mapping of the Anderson model on Random Regular Graph (RRG) onto the"multifractal"RP model and find exact values of the stretch-exponent $\kappa$ depending on box-distributed disorder in the thermodynamic limit. As another example we consider the logarithmically-normal RP (LN-RP) random matrix ensemble and find analytically its phase diagram and the exponent $\kappa$. In addition, our theory allows to compute the shift of apparent phase transition lines at a finite system size and show that in the case of RP associated with RRG and LN-RP with the same symmetry of distribution function of hopping, a finite-size multifractal"phase"emerges near the tricritical point which is also the point of localization transition.


Introduction
Random matrix theory (RMT) has tremendous range of applications in physics spanning from the spectra of heavy nuclei to physics of glassy matter. Specifically, it was a basis for mesoscopic physics with a very successful applications to the description of ergodic phases in disordered condensed matter systems both on the single-particle level, as in diffusive metals, and in the many-body problem of thermalization in generic isolated quantum systems. In the latter problem, the paradigm of RMT has been further extended to the eigenstate thermalization hypothesis (ETH) [1][2][3] claiming the thermalization in quantum many-body systems isolated from a thermal bath. Indeed, ETH uses the fact [1] that close-in-energy eigenstates of a chaotic Hamiltonian hybridize with essentially random phases under a small perturbation and thus can be represented by eigenvectors of the Gaussian random matrix ensembles (GRME) or by the so-called random pure states. In both descriptions the states have the (nearly) i.i.d. Gaussian-distributed wave function coefficients and they are used not only for the description of the equipartition over degrees of freedom, or ETH, but also provide the ground for thermalization based on the scrambling [4] and entanglement [5] of these degrees of freedom. In all these fields, RMT provides an excellent statistical description of the problems highlighting the universal properties shared by all kinds of chaotic systems.
The phenomenon of Anderson localization was the simplest example of the failure of the classical RMT of Wigner and Dyson (WD RMT) on the single-particle level. In order to describe Anderson localization that happens in a certain (laboratory) frame one has to introduce the basis preference in RMT. The simplest modification is the Rosenzweig-Porter RMT [6] which is nothing but the WD RMT with the special, stronger fluctuating diagonal matrix elements setting the basis preference.
During the recent decade the phenomenon of the many-body localization (MBL) was discovered [7,8] and quickly developed into a whole field of the thermalization breakdown in generic isolated quantum many-body systems under the sufficiently strong quenched disorder [9]. At such conditions ETH is strongly violated raising a question about a model which is as simple as WD RMT but capable of describing the generic features of the MBL and pre-MBL states.
One may start by the Fock/Hilbert space formulation of the MBL problem for models with local interactions in the coordinate basis. A crucial simplification on this way was suggested in the seminal paper [7] where it was assumed that the Fock/Hilbert space of such systems form a hierarchical graph like the Cayley tree (CT) or the Random Regular Graph (RRG) [10] and the MBL state is the Anderson localized state in the Hilbert space. This boosted the field of Anderson localization on hierarchical graphs with the local tree-like structure suggested [11][12][13][14] as a toy model of MBL in the Hilbert space of such many-body models.
However, in a number of recent works [15][16][17][18][19] such a picture has been put in doubt. Namely, due to the locality of the Hamiltonian in the coordinate (often one-dimensional) space, absence of transport in the coordinate space does not mean localization in the Hilbert space. A simple model of the small fraction 0 < f < 1 of frozen spins [15] in a 1D s-spin chain of length L with nearest neighbor interaction shows that the wave function support set in the Hilbert space has a volume (2s + 1) (1−f ) L = N (1−f ) , where N = (2s + 1) L is the total volume of the Hilbert space, i.e. it has a Hausdorff dimension 0 < D = 1 − f < 1 and thus is a fractal. At the same time frozen spins block the spin-excitation transport completely. A little more elaborate model that involves the localized one-particle states with finite localization radius and the Slatterdeterminant many-body states shows that the inverse-participation ratio of such many-body states corresponds to the fractal dimension 0 < D q < 1 even without interaction [18]. The role of interaction in this case is played by the entanglement and Pauli Principle encoded in the structure of Slatter determinants. Small repulsive interaction does the same job for vanishing localization radius. Thus multifractality is intrinsic to many-body interacting systems. Such states are neither localized nor fully ergodic. They occupy a sub-extensive part of the whole accessible Hilbert space and form an eigenstate multifractality in the Hilbert space in a whole MBL phase.
The above mentioned Anderson models on the hierarchical graphs have been suggested [11][12][13] as the toy models that support such a robust non-ergodic extended (NEE) phase, but the existence of such a phase in the thermodynamic limit has been disputed during recent years [13,14,[20][21][22][23][24].
In contrast, in the random-matrix theory such systems with the NEE phase have been recently found. It appears that the simplest basis non-invariant Rosenzweig-Porter ensemble [6] supports the NEE phase if the variance of the i.i.d. diagonal matrix elements is parametrically larger H 2 nn = N γ |H n =m | 2 than that of the off-diagonal ones drawn from the GRME, where N is the matrix size. It was discovered in Ref. [25] with the further rigorous mathematical proof in Ref. [26] that the eigenvectors of the RP model with γ in the interval 1 < γ < 2 have a fractal support set with the Hausdorff dimension D = 2 − γ. This example has been intensively studied [27][28][29][30][31][32] over the past few years.
It is conceptually important that the non-trivial NEE phase in the RP model arises only if the statistics of matrix elements depends on the matrix size in a proper way described above. We will show in this paper that the models with short-range, size-independent couplings, such as models with the nearest-neighbor coupling on RRG, can be mapped onto the infiniterange RP models with the distribution of matrix elements depending on the matrix size. This dependence is the price for a relative ease to treat the infinite-range matrix models. The first demonstration of such type of mapping was done by Efetov [33] when he reduced the d-dimensional Anderson model with a sparse matrix Hamiltonian to the infinite-range GRME for frequencies ω smaller than the Thouless energy by means of the supersymmetric non-linear sigma-model. In this case rescaling of all the matrix elements by any N -depending factor does not change the properties of GRME. However, to make the spectral bandwidth finite E BW = 1 one has to make N -dependent the Gaussian statistics of matrix elements in GRME.
Reducing the Anderson models on the hierarchical graphs to a certain kind of the RP models is the next non-trivial step. However, the Gaussian RP (GRP) model considered in Ref. [25] is rather oversimplified to be a good candidate for this goal.
Because of the infinite connectivity and self-averaging of the matrix element fluctuations the exact solution for this model can be written in the standard Breit-Wigner form [27,30,31] with the certain normalization constant C and the level broadening Γ given by the escape rate from the Fermi Golden rule The self-averaging of the above sum in the thermodynamic limit (confirmed by the cavity method calculations [27,31]) immediately supports the Breit-Wigner approximation and the diffusive transport [29,32] in such a model. In more realistic local many-body systems, the corresponding matrix in the reference Hilbert space basis (e.g. in the "computational" basis of strings, or "configurations" of up/down spins in the spin-1/2 spin chain) is sparse and the effective connection of far-away configurations (strings) is determined by the presence of a long series of quantum transitions and anomalously strong resonances. These effective long-range hopping amplitudes are in general correlated and fat-tailed distributed [17,18,34].
The main idea [17,35] of a random matrix representation of the many-body problem in the Hilbert space is to replace the sparse matrix of the many-body Hamiltonian by the full matrix of the RP type with off-diagonal elements representing the effective connections. This is equivalent to reduction of a quantum disordered problem on a hierarchical graph to that on a complete graph with a proper statistics of matrix elements. This mapping is essentially a reduction of an Anderson localization problem on an infinite-dimensional lattice to that on a zero-dimensional one, in which the underlying hierarchical graph structure disappears but remains encoded in the statistics of matrix elements.
An important simplification comes from the "infinite-temperature" condition. Indeed, the statistics of highly-excited states is supposed to be the same for a broad band of such states which occupy a certain part of the Hilbert space. In this fraction of the Hilbert space there is an approximate statistical homogeneity over the corresponding configurations, very much like the statistical homogeneity over sites on the Random Regular Graph. This homogeneity is a crucial point for the application of the RP-type random matrix theory with statistical homogeneity over matrix indices to realistic many-body systems. Apparently, a reduction to RP random matrix theory is not an exact mapping. Some information is inevitably lost after the "folding" of the hierarchical Hilbert space into a "zero-dimensional" RP model. For instance it is no longer possible to define the diffusion over configurations in the Hilbert space but it is possible to study the time-dependence of survival probability in a given configuration which contains the same information as the diffusion in the (Hilbert) space. It is also possible to study the Hilbert space multifractality via distribution function of the many-body wavefunction coefficients.
The main difference between the Gaussian RP model and those effective RP models that arise in the localization problem on hierarchical graphs [35] and in realistic many-body systems [17] is that the distribution of off-diagonal matrix elements in the latter is heavily tailed. One consequence of that is that the wave functions in the NEE phase of such RP models are genuinely multifractal and not mono-fractal, as in Gaussian RP model [25].
In this paper we concentrate on another important consequence of the broad distribution of off-diagonal matrix elements in such RP models. Namely, we show that the survival probability in a broad region of parameters shows the stretch-exponential behavior rather than the simple exponential one present in the Gaussian RP model [32]. Since the diffusion on a hierarchical graph (like in a well-known example of a Cayley tree) results in the simple exponential decay of survival probability, the stretch-exponential behavior should be associated with the subdiffusion. This observation is important because even on the extended side of the MBL transition in some many-body systems the diffusive character of transport is not completely confirmed by several numerical calculations [36][37][38][39][40][41][42][43][44]. This leaves a room for the subdiffusive dynamics, probably by forming the corresponding phase with the anomalously slow (glassy) transport properties due to some bottleneck originating from the Griffiths-type physics (see reviews [45,46] and references therein).
In this paper we demonstrate that in some important and quite wide class of RP models with a heavily-tailed, "multifractal" distribution of hopping matrix elements such phases with a stretch-exponential dynamics of survival probability are generically present. We derive analytically the lines of phase transitions between the exponential (E) and stretch-exponential (SE) dynamics, as well as the expression for the stretch-exponent κ for a generic "multifractal" RP model.

The model
We consider a RP model with a "large-deviation", or "multifractal" distribution F(H mn ) of i.i.d. off-diagonal (hopping) matrix elements of the form: with the typical coupling H typ ∼ N −γ/2 , γ > 0 which is polynomially small in the dimension N of the Hilbert space of the system. Below we show (see also Ref. [35]) that it is exactly the form that emerges in the effective RP model for the Anderson localization model on RRG. In a generic many-body problem, the non-trivial N -dependence in Eq. (3) can be understood from the following observation: when one forms the effective long-range hopping in the Hilbert space the transition amplitude is typically determined by the product of individual local hops and thus decays exponentially H mn ∼ e −Armn with the Hamming distance r mn between configurations m and n, while the number of available configurations n to hop from a certain configuration m grows exponentially with r, resulting in F(H mn ) ∼ e Brmn . As the most frequent inter-configuration distances r max are of the order of the Hilbert space diameter, r r max ∼ ln N , this leads to H mn ∼ N −g/2 and F(H mn ) ∼ N B . Finally, taking into account large deviations of highly fluctuating matrix elements H mn from the typical ones, one has in general to write F (g) instead of F = B. The properties of the function F (g) are determined by the underlying local structure of the sparse-matrix Hamiltonian in the Hilbert space and may be subject to certain conditions and symmetries.
Note that the diagonal matrix elements ε n = H nn in the MF-RP model are identical to onsite energies in the underlying short-range model. In what follows we consider the box-shaped distribution of on-site energies: In contrast to the distribution of the off-diagonal matrix elements, Eq. (3), the distribution Eq. (4) is independent of the matrix size N and is not tailed. Like in the previously considered Log-normal (LN) [35,47] and Levy [30,48] extensions of the RP model (which are the particular cases of F (g) being quadratic or linear function of g), the generic MF-RP has a rich phase diagram with the robust genuinely-multifractal phase (see Fig. 1). This phase diagram can be easily understood in terms of the Breit-Wigner formula (1) (see, e.g., [31,47,49]) where the broadening Γ is determined self-consistently via the cavity method (see, e.g., [48]).

Dynamical phases
We focus on the dynamical properties of the generic MF-RP model in order to reveal anomalously slow transport and possible subdiffusion -diffusion phase transition. Our analytical approach is based on the Wigner-Weisskopf approximation [30] which allows us to demonstrate the breakdown of the Fermi Golden rule approximation (2) and explicitly calculate the survival/return probability R(t) for a quantum system to return back to the initial configuration. The latter is shown to decay either exponentially (corresponding to the Fermi Golden rule and the diffusive transport) or stretch-exponentially corresponding to the sub-diffusion (similar to those in RRG [23,24]). This stretch-exponential (SE) decay with the exponent 0 < κ < 1 persisting in the extensive time interval 1 (E 0 t) κ ln N , is an emergent property of a generic MF-RP model. It may be present both in a non-ergodic (multifractal) and in a weakly-ergodic phases 1 The only difference in the SE decay in these two phases is in the scaling with N of the characteristic energy (or time) scale E 0 : in the multifractal phase this scale E 0 /E BW (measured in units of the total spectral band-width E BW ) decays to zero in the thermodynamic limit, while in the weakly-ergodic phase this ratio E 0 /E BW ∼ O(1) remains finite.
The spreading of the stretch-exponential decay of R(t) over the extensive time interval allows one to define (in the thermodynamic limit N → ∞) a set of dynamical phases characterized by the exponent κ. In a general case we identify the following phases: (i) the "frozen-dynamics phase" (F), κ = 0, 1 Here we call wave-function "weakly ergodic" if it occupies a finite fraction of the total Hilbert space. Such states play an important role in several recent works [16,23,24,[49][50][51][52][53]. Frozen-dynamics (F), Stretch-Exponential (SE) and Exponential (E). Albeit RP ensemble corresponding to RRG is not exactly log-normal, it has the same phase diagram as LN-RP with p = 1 and γ = 2λ typ / ln K ≥ 1. The corresponding path in the parameter space respecting the symmetry Eq. (6) is shown by a red dashed line. As the Lyapunov exponent on a Cayley tree λ typ → (1/2) ln K in the clean limit, the part of the phase diagram γ < 1 is not relevant for RRG. Thus in the thermodynamic limit on RRG there is only the localized and the weakly-ergodic phase with a stretch-exponential dynamics. In the inset: finite-size phases in the vicinity of the tri-critical point. Like on RRG, a finite-size multifractality (FSMF) emerges near the localization transition. Arrows show the direction and "speed" of evolution as the system size N increases.
Note that the F phase does not exclude a polynomial in time relaxation of the initial configuration and thus it is not necessarily coinciding with the localized phase.
Furthermore, we derive explicit conditions for each of the phase to occur in terms of the function F (g). For the particular example of the LN-RP model, F (g) = −(g − γ) 2 /(4pγ), we show that the transitions between the above dynamical phases do not coincide with the ones established in Ref. [47] on the basis of analysis of ergodicity and its violation. We show that both SE/E and F/SE transitions may occur in the multifractal phase as well as in the weaklyergodic one (see Fig. 1) and thus the above classification of dynamical phases reflects a property of the eigenfunction statistics which is complementary to ergodicity and multifractality.
An important particular case is the LN-RP model with the symmetry p = δg 2 /(2 g ) = 1 which is a particular case of a general symmetry This symmetry is also present in MF-RP model associated with RRG (see Sec. 7). It corresponds to the trajectory shown by a red vertical dashed line on the phase diagram of Fig. 1. This trajectory passes through the tri-critical point γ = 4, p = 1 on the (γ, p) plane where the localized, multifractal and weakly-ergodic phases merge together [47]. The dynamical analysis of this paper shows that in this very point the merging of the stretch-exponential and the frozen-dynamics phases also occurs.

Finite-size multifractality
Our analytical approach based on the Wigner-Weisskopf approximation [30] allows us to find the finite-size deformations of the phase diagram. In the inset of Fig. 1 we show that a "finitesize multifractal phase" (FSMF) appears in the vicinity of the tricritical point just filling in the emerging gap between the localized and weakly-ergodic phases. The width of the gap in terms of the diagonal disorder strength W is given by: where c ∼ 1, W c is the localization transition point in the thermodynamic limit and W AT and W ET are the finite-size apparent transition points for the localization and the ergodic transitions, respectively. This allows to define the characteristic length: such that for 1 N N MF (at fixed disorder strength W < W c ) one observes approximately a multifractal eigenfunction statistics, while for N N MF it becomes weakly-ergodic. This finding finally resolves the long-lasting dispute about the existence of the non-ergodic extended phase on RRG and gives an expression for the crossover system size N MF that is much larger than the correlation volume at the localization transition ln N c ∝ |1 − W/W c | −1/2 [21,47,[54][55][56]. This new length arises because of the proximity of the tri-critical point (shown by a red point in Fig. 1) to the true multifractal phase emerging when the symmetry, Eq. (6), is infinitesimally broken. The corresponding range of disorder strength (7) (with the parameters appropriate for K = 2 RRG) for the currently available system sizes N = 10 3 . . . 10 5 shown in Fig. 2 is in good agreement with the finite-size multifractality observed in the exact diagonalization earlier [11][12][13]57].
Finally, in order to confirm the applicability of our mapping and analytically derived results for the stretch exponent κ, we studied in great detail the particular cases of LN-RP and MF-RP associated with RRG. In the first case we compute analytically the phase transition lines and the dependence of the stretch exponent κ on the parameters γ and p. In the second case, we find the function F (g) from the exact theory of localization on a Cayley tree of Ref. [58] and on its basis we compute the dependence the stretch-exponent κ on the disorder strength W . In both cases we found a good agreement with the exact diagonalization numerics whenever convergence to the thermodynamic limit may be reached at available system sizes or the extrapolation of data to N → ∞ limit can be performed under controllable conditions. The paper is organized as follows. In Sec. 2 we remind the criteria for the localized, multifractal and ergodic static phases in the models with the full random matrix Hamiltonians based on the principles discussed in Refs. [47,49] and apply these principles to the MF-RP model with a generic function F (g). Next, in Sec. 3 we present the main steps of the derivation of the Wigner-Weisskopf approximation for the return probability R(t) and apply it in Sec. 4 to calculate the mean R(t) in the case of MF-RP. In the latter section we will show the general nature of the stretch-exponential decay of return probability as an emerging phenomenon and derive in Sec. 5 the transition lines between the frozen-dynamics, the stretch-exponential, and the exponential dynamical phases. Section 6 is devoted to consideration of the overlap correlation function and its relation to the stable distribution emerging beyond a certain highfrequency cutoff.
In the rest of the paper we apply the developed above general theory of MF-RP to a particular example of MF-RP that stems from the Anderson localization model on RRG (Secs. 7.1, 7.2) and inherits the symmetry Eq. (6) from the corresponding Cayley tree. More precisely, we relate the properties of the function F (g) to the eigenvalue ε β of the linearized transfer-matrix operator for the Cayley tree and to the corresponding Lyapunov exponents (Secs. 7.3, 7.4). In the end of the paper we present the results of extensive numerical simulations for RRG and the corresponding LN-RP model with a symmetry parameter p = 1 (Sec. 8).

Localization, ergodic and fully-weakly ergodic transitions
Before studying the dynamical phases of the model (3), (4) we briefly review its static phases shown in Fig. 1.
In Refs. [47,49] we formulated the criteria for the phase transitions seen in the eigenfunction statistics of full random matrices. We identified localized, multifractal, weakly-ergodic and fully-ergodic phases and the transitions between them: Anderson localization transition (AT), ergodic transition (ET) between the multifractal and ergodic phases and the fully-weakly ergodic (FWE) transition between the the fully-and the weakly-ergodic phases. Remarkably, all the corresponding criteria can be written in a universal way using the properly defined averages of the off-diagonal matrix elements H nm and comparing them with the mean level spacing δ = E BW /N .

• Localization transition:
• Ergodic transition where ... E BW denotes the upper cutoff at |H nm | ∼ E BW , and C is the normalization constant The quantities |H nm |, its upper cutoff and δ are measured in units of the total spectral bandwidth E BW .
In the saddle-point approximation the integrals (9) and (10) are dominated by g = max g 1/2 , 0 and g = max (g 1 , 0), respectively, while the normalization integral (12) is given by g = g 0 . Here we denote as g β the solution of equation The corresponding phase diagram for the LN-RP in (γ, p) plane (first obtained for this model in Ref. [47]) is shown in Fig. 1 where the blue, green, and red solid lines show the corresponding transition points for the Anderson ergodic and fully-weakly ergodic transitions 3 Wigner-Weisskopf approximation.
In order to complement the above static phase diagram by the dynamical phases, we turn in this section to the main goal of this paper: the calculation of return probability in a generic MF-RP model exhibiting an emergent stretch-exponential decay. We start by the description of the Wigner-Weisskopf (WW) approximation for the quantum dynamics [59] following the rout of Ref. [30].
In this approach the time-dependent Schroedinger equation: with the initial population only on site n is solved approximately as follows. As the first step we introduce the "hopping representation": where ε i = H ii and the r.h.s. contains only off-diagonal matrix elements: Now in the first order of perturbation theory: while the equation for Φ (n) (n, t): is solved in the Markovian approximation Φ (n) (n, τ ) → Φ (n) (n, t) justified by the stronglyoscillating integrand: Thus the return probability is given by: where we approximate: with δ(N ) being the mean level spacing. The crucial point of approximation Eqs. (23), (24) is that the dependence on time of a random decay rate Γ(t) depends non-trivially on the interval of summation n t ∼ 1/t of the random square of the matrix element |H nj | 2 , see Fig. 3. If all matrix elements |H jn | are typical, the sum in Eq. (24) is proportional to |H typ | 2 n t , that is Γ(t) = const. Thus we obtain the exponential decrease of return/survival probability with time. However, if the distribution of |H jn | is broad, then large matrix elements occur less frequently as time goes on and n t decreases. This makes the sum decaying faster than 1/t and Γ(t) decreases with time (see Fig. 3).
In the next section we present the detailed derivation of this fact and show that no matter what is the function F (g), the mere multifractal character of the distribution Eq. (3) makes the average return probability stretch-exponential in a broad interval of time. At the same  24), is shrinking. This leads to reduction of the interval of values of |H nm | at the tails of the distribution (shown by progressively shrinking yellow shaded areas) that one can "catch" having n t "trial events". As the decay rate Γ(t) is determined by rare large |H nm |, decreasing the maximal value of the set of n t random matrix elements |H nm | (shown by black dashed lines) results in decreasing Γ(t) with time, i.e. the behavior of R av (t) slower than exponential. For a distribution with no tails (e.g. when the purple parabola is very narrow) the typical maximal value of |H nm | is time-independent, and the dynamics of R av (t) is exponential.
time, the typical value of R(t) = e −tΓ(t) remains exponential, since it is determined by the typical values of |H jn |.
Concluding this section we would like to discuss the applicability of the Wigner-Weisskopf approximation. The main approximation here is the first order of perturbation theory in Eq. (20). This approximation works well for the full random matrix Hamiltonian (coordination number N → ∞) and fails for the sparse hopping matrix like in the initial Hamiltonian of the Anderson model on RRG (finite coordination number K + 1). The situation here is similar to the perturbative derivation of multifractality in the Gaussian RP model [25] which appeared to be exact [26,27]. That is why one should first derive the RP model equivalent to the short-range sparse problem and then apply the Wigner-Weisskopf approximation to this RP model. For the latter model WW approximation works as long as n t 1, e.g. as long as

Averaging of return probability.
Given that all matrix elements |H nj | 2 in Eq. (24) are independent random variables, the standard way of computing the distribution function G(Γ) of −(1/t) ln R(t) = Γ(t) is to compute first its characteristic function Θ(v) = e ivΓ G(Γ) dΓ: Then the average return/survival probability is while the typical one is The typical return/survival probability R typ is pure exponential, as n t ∝ t −1 . The situation is less trivial with the popular numerical measure of anomalous transport, namely the averaged survival/return probability R av (t). Let us consider this measure in detail. First of all we represent the average in Eq. (26) in terms of the distribution function Eq. (3): The integrand in this equation is vanishingly small at g < 2τ , while at g > 2τ one can write: Thus we obtain: Here we consider the case when: 2τ < g 0 .
As we will see below this time interval is relevant for all dynamical phases. Eq. (30) implies that (i) the first integral in Eq. (29) is dominated by the saddle-point g = g 0 and unity with a polynomial correction of the order of N F (2τ ) 1, while (ii) in the second integral the main contribution comes either from the vicinity of the lower cutoff g = 2τ , when the saddle point g = g 1 < 2τ is outside of the region of integration, or from the point g = g 1 itself.
The latter case, g 1 > 2τ , recovers the standard diffusion, R av (t) = e −Γavt , with the Fermi Golden rule decay rate The polynomial correction N F (2τ )−2τ to Γ av is subleading in this case due to the property (13) of g 1 : max In the more interesting case when the saddle point g = g 1 < 2τ is outside of the region of integration in (29), the return probability is governed by the boundary term At this point it is clear that the dynamics of R av (t) separates in three intervals, see Fig. 4 (i) At small times (shown by light blue) the evolution is perturbative as Eq. (32) shows very fast evolution, but it will be immediately terminated at F (2τ ) − τ + ∆ ∼ ln ln N/ ln N by the finite-size saturation of return probability of the order of − ln R av (t → ∞) ∼ ln N (see below), while (iii) The range of non-perturbative dynamics of our interest (shown by pink) is located in the vicinity δτ 1 of the point given by the solution of the equation Here the crucial step in understanding of the emergent stretch-exponential behavior is the expansion in δτ = τ − τ * 1 of F (2τ ) in Eq. (32) in the vicinity of τ = τ * (upper panel of where The expression (35) defines the only characteristic energy scale E 0 separating the three dynamical regions discussed above: where τ * is the solution to Eq. (34). In ergodic phases E 0 is of the order of the total bandwidth E 0 ∼ E BW = N δ(N ), while in the multifractal phase E 0 has a physical meaning of a width of the mini-band in the local spectrum. In this case it decreases polynomially with increasing N , i.e. τ * > 0.
In order to find the scaling exponent τ * one should use the knowledge about the scaling of the total bandwidth E BW ∼ N 1−∆ . It is well-known [25,47,53] that in the non-ergodic (e.g. localized and multifractal) phase E BW = W/2 ∼ N 0 , so that ∆ = 1. However, in the weaklyor fully-ergodic phases ∆ < 1. Indeed, using the fact that in this case E 0 ∼ E BW one obtains for ergodic phases ∆ = 1 + τ * , where −1 < τ * < 0. Therefore: Now, Eq. (34) can be written in a general form: For us relevant is the smaller root of Eq. (39) with 2τ * < g 0 where the dynamics is not yet saturated.
Eqs. (36), (37), (39) uniquely determine the stretch exponent 0 < κ < 1 in the survival/return probability: This is the main result of the paper. To find the true region of its validity we note that in a finite system there is a saturation of R av (t) = α |ψ α (n)| 4 ∼ N −D 2 at large times expressed in terms of the complete set of eigenstates ψ α (n) of the MF-RP problem. This saturation is beyond the Wigner-Weisskopf approximation and this reduces its validity range from t 1/δ(N ) to It is important to note that this domain of validity corresponds to 0 < δτ κ −1 ln ln N/ ln N 1 (upper panel of Fig. 4). Thus despite the width of the interval of E 0 t diverges as N → ∞ (lower panel of Fig. 4), the width of δτ shrinks to zero, hence the expansion of F (g) around g = 2τ becomes exact in the thermodynamic limit. This is the consequence of the "multifractal" character of the distribution Eq. (3). It is generic to any distribution of this type, no matter what the function F (g) ∼ N 0 of g ∼ N 0 . Thus we have shown that the stretch-exponential behavior of the survival/return probability is an emergent property of a large "multifractal" RP random matrices when the thermodynamic limit N → ∞ is taken before E 0 t → ∞.

Exponential/stretch-exponential and frozen-dynamics transitions
In order to obtain the stretch-exponential behavior Eqs. (40) it is crucially important that the lower cut-off 2τ ≈ 2τ * in the integrals in Eq. (29) is larger than the saddle point g 1 in the second integral (but smaller than the saddle point g 0 in the first one). If this is not the case then the survival probability R av (t) is pure exponential, κ = 1 with the decay rate (31) given by the Fermi Golden rule. Thus the transition between the exponential and the stretch-exponential dynamics happens at: Eq. (42) is general for any MF-RP model with the function F (g) obeying only the normalization condition F (g 0 ) = 0. It can be rewritten in a physically transparent form similar to (9), (10), (11). For this we note that the LN-RP distribution with the parabolic function F (g) obeys a special symmetry: This symmetry includes the symmetry, Eq. (6), as a particular case, but it is wider than that which allows only one value of p = δg 2 /(2 g ) = 1. However, it does not cover all the possible functions F (g).
Using the normalization F (g 0 ) = 0, the symmetry (43), and the definition H typ = N −g 0 /2 , we reduce the condition, Eq. (42), to the one covering any LN-RP and RRG We believe that this condition is the basic one for all the physically-motivated applications. The physics behind it is that at H typ δ(N ) all the states are typically hybridized and the Fermi Golden Rule holds leading to the pure exponential dynamics. In contrast, at H typ δ(N ) only atypically large couplings with |H mn | > δ(N ) can hybridize the corresponding states, all other couplings may be safely neglected. They constitute an effective sparse matrix of couplings, and so the stretch-exponential behavior of R av (t) emerges exactly at the point where the effective matrix of couplings becomes sparse (see [60]). Now let us consider the condition that the stretch-exponent κ turns to zero. At this point the dynamics of survival probability is at most some power law of time. Such an ultra-slow dynamics is essentially a glassy behavior which we will refer to as the frozen-dynamics phase.
Its physical meaning is related to the mobility edge. Indeed, according to Eq. (52) the return probability saturates at large t when different states close in the energy do not overlap in an observation point of an infinite system. Obviously, this happens when the states are localized. However, it may happen also if all extended states which are close in the energy do not populate the observation point. The reason is that there is a remote in energy state beyond the mobility edge that is localized exactly at this point. Then by orthogonality of states all other (extended) states must have a population hole at this site [61]. This is a somewhat exotic situation, as the localized state must be localized on one single site in the limit N → ∞, otherwise the orthogonality condition may be met by matching of phases of wave function (or its signs). However, this is known to happen in the RP models [32].
According to Eq. (36) the stretch-exponent κ turns to zero when: Using Eq. (34) one can express the condition for the transition to the frozen-dynamics phase in the form very similar to the Anderson transition (9) with the only difference of the absence of the upper cutoff at |H mn | ∼ O(N 0 ) which is present in (9). This immediately means that F/SE-transition coincides with the Anderson transition as soon as 2τ * = g 1/2 ≥ 0, i.e. when the cutoff is not important. On the other hand, τ * > 0 implies the multifractal phase. We conclude, therefore, that in a generic MF-RP model the F/SE transition is not possible inside the multifractal phase. Indeed, it coincides with the Anderson localization transition when the multifractal phase is present. However, there is no prohibition for the F/SE transition to be inside the weakly-ergodic phase, and, indeed, it happens there (see Fig. 1). The fact that the criteria of all the phase transitions can be presented in the common form (9), (10), (11), (44), and (46) allows to establish a general sequence of transitions as disorder increases. By writing the corresponding general inequalities for averages of |H nm | and |H nm | 2 , we conclude that the SE/E transition (44) has to appear at smaller disorder compared to the F/SE transition (46), while both of them should be located between the Anderson (9) and the fully-weakly ergodic transition (11), irrespectively of the location of ergodic transition.
For this model the function F (g) is parabolic: and one can easily find The smaller of the two solutions of Eq. (39) corresponding to this function is: Then solving the equation 2τ * = g 1 one obtains for the E/SE transition two different expressions depending on whether the transition happens inside the weakly-ergodic or inside the multifractal phases: Analogously, from 2τ * = g 1/2 one obtains the line of transition to the frozen-dynamics phase: Eqs. (50), (51) should be complemented by the lines of the Anderson localization (14), the ergodic transition (15), and the transition between the weakly-ergodic and the fully ergodic phases (16) first obtained for this model in Ref. [47]. All these transition lines are shown in Fig. 1. Please notice a complexity of this diagram which possesses four multi-critical points and seven different phases, including three types of weakly-ergodic phases and two types of multifractal phases. Notice also that the phase with frozen dynamics appears inside the weakly-ergodic phase, and the phase with the stretch-exponential dynamics appears both inside the weakly-ergodic and inside the multifractal phase. This implies that the classification by dynamical phases reflects different, complementary features of eigenfunction statistics compared to the classification by the extent of ergodicity violation. While the localized, multifractal, weakly-and fully-ergodic phases are classified according to the fractal dimension of the support set of a single eigenfunction and the fraction of the occupied sites in it, the classification by dynamical phases accounts for the correlation between eigenfunctions of different energies.

Overlap correlation function and Levy stable distribution
The survival/return probability R av (t) is the Fourier-transform of a correlation function where ψ α (r) is the eigenfunction corresponding to eigenenergy E α at a site r. To show this one expands the wave function ψ (n) (r, t) = α a α ψ α (r) e −iEαt in the eigenbasis. If the coefficients are chosen as a α = ψ α (n), the completeness imposes the correct initial conditions ψ (n) (r, t = 0) = δ r,n . Then one immediately obtains In order to eliminate the effect of level repulsion, we normalize C(ω) by the two-level correlation function Q(ω) = N −1 α,β δ(ω − E α + E β ) and define the eigenfunction overlap correlation function: The functions C(ω) and Q(ω) both have a correlation hole at ω δ due to level repulsion, while for ω δ the function Q(ω) = δ −1 is a constant. Therefore, for exact diagonalization at a finite N the function K(ω) is preferred over C(ω), as it allows to separate the effects of eigenfunction statistics from the trivial level repulsion. This poses an advantage, as it extends the region of the stretch-exponential behavior to larger times compared to direct computation of survival probability in the time domain.
Taking the inverse Fourier transform of Eq. (53) one can find K(ω) at ω δ:

K(ω)
The corresponding parts of the curves for K(ω) and µ(ω) (blue solid lines) are the part of the symmetric about ω = 0 Levy α-stable distribution with α = κ and its log-log derivative, respectively (their continuations are shown by the blue dashed lines) . The parts of the curve K(ω) and µ(ω) on the right of the minimum ω min ∼ E 0 (orange solid lines) are K(ω) ∼ ω −1 (ln(1/ω)) −c , which was analytically derived for RRG in Ref. [21], and its log-log derivative (their continuations are shown by the orange dashed lines). This part of K(ω) is relevant for the Fourier-transforming at shorter times E −1 In order to get a more detailed information we also consider the log-log derivative of K(ω): It gives the "running power" of the locally power-law function K(ω).
Taking into account a well-known characteristic function of the Levy α-stable distribution: we conclude from Eqs. (55),(40) that the stretch-exponential R av (t) corresponds to the symmetric Levy α-stable distribution: with α = κ, the skewness parameter β = 0, the shift parameter µ = 0 (the distribution is symmetric about ω = 0) and the scale parameter b = E 0 . This result is valid for ω corresponding to the region of validity of SE, Eq. (41): For E BW ω E 0 the main contribution to the Fourier-transform of R av (t) comes from t ∼ ω −1 : This is consistent with the analytical result of Ref. [21] for RRG: The parts of K(ω) corresponding to the Levy α-stable distribution and Eqs. (60),(61) are shown in Fig. 5 by the blue and the orange solid lines, respectively. A superficial look at K(ω) could give a wrong impression that the entire falling part of K(ω) is described by K(ω) ∝ |ω| −1 , possibly with some logarithmic corrections [13,21,25]. To see that this is not so, one has to plot the log-log derivative of K(ω), Eq. (56), shown on the lower panel of Fig. 5. On this plot the part corresponding to the log-log derivative of the Levy α-stable distribution (blue solid line) matches with the log-log derivative of Eq. (61) (orange solid line) near the minimum of the log-log derivative µ(ω) of K(ω). Either of the two parts of µ(ω) do not have a minimum and their continuations (shown by the dashed lines of the corresponding color) tend asymptotically to −(1 + κ) and −1, respectively. However, at the matching point a minimum in µ(ω) emerges with the minimal value of µ(ω) being |µ min | ≤ 1. Thus the minimum in µ(ω) is a natural marker of termination of the emergent Levy α-stable distribution in K(ω).
The numerical results for µ(ω) shown in Appendix A (see Fig. 9) demonstrate that the minimum in µ(ω) indeed exists in the weakly-ergodic, stretch-exponential phase of LN-RP model with p = 1, as well as on RRG. However, in the multifractal, stretch-exponential phase of LN-RP model at p = 1/2 there is, instead, a shoulder. This signals again of the two different regimes in K(ω), the one that is described by the Levy α-stable distribution and another one, which, however, is different from Eq. (61).

Correspondence between Anderson model on RRG and MF-RP
In the rest of the paper we apply the general idea of mapping to MF-RP model to a particular example of the Anderson tight-binding model on a random regular graph (RRG). We argue that a RP random matrix model with a properly chosen "multifractal" distribution of i.i.d. off-diagonal (hopping) matrix elements, is in the same universality class as RRG. The corresponding function F (g) in the distribution, Eq.(3), is determined by the eigenvalue ε β of the linearized transfer-matrix operator for the corresponding Cayley tree (CT) with the branching number K = m − 1, where m is the coordination number of RRG. It obeys a special symmetry that stems from the basic Abou-Chacra-Thouless-Anderson (ACTA) symmetry ε β = ε 1−β [58]. It is this symmetry which ensures that the Anderson transition point in the RP ensemble equivalent to RRG is in fact the tri-critical point of a more generic RP model where the localized, multifractal and ergodic phases meet together [47]. This is the reason for existence of a transient finite-size multifractal phase on RRG [11,12] with ergodicity of states restored in the thermodynamic limit N → ∞ [14,21]. However, the extended states on RRG are not fully ergodic like in the classic Wigner-Dyson random matrices. We show that the survival probability R av (t) in the RP ensemble corresponding to RRG is stretch-exponential in the entire extended phase [47], like it was earlier conjectured in Ref. [23]. This statement rests on the proposed equivalence between the RP and RRG and on the ACTA symmetry and is valid for any distribution of on-site energy with no fat tails and for any branching number K of RRG.
We show that effective RP ensemble which is equivalent to RRG is of the form Eq. (3). However, it must also respect the ACTA symmetry (6): We will show below that this symmetry enforces that the F/SE transition on RRG coincides with the Anderson localization transition and the SE/E transition coincides with the point where the Lyapunov exponent λ typ = −(1/2)∂ β β | β=0 on the corresponding Cayley tree takes its minimal value λ typ = (1/2) ln K characterizing the fully-ergodic phase [62]. For the conventional Anderson model on a CT with one orbital per site this value is reached only in the limit of vanishing disorder, so that the entire weakly-ergodic phase on RRG is characterized by a stretch exponential dynamics 2 . Furthermore, we will show that SE phase on RRG takes exactly the place of MF phase on a corresponding Cayley tree. Thus when passing from the finite CT to the RRG (loosely speaking "by connecting the leaves"), the multifractal phase is transformed [20,63] into the weakly-ergodic phase with stretch-exponential dynamics. The latter is thus the remnant of the finite-size multifractality discovered in Refs. [11,12].
We would like to emphasize that the form of the distribution F(x = H nm ) corresponding to RRG is strongly disorder-dependent and is not the log-normal [47] in general. Thus the details of the phase diagram do depend on the precise form of the dependence of β on β and the function F (g) that follows from it. However, the topology of the phase diagram is generic and is well represented by the "symmetric" log-normal RP model of Fig. 1 with p = 1.

Mapping of tight-binding Anderson model on RRG onto RP model with long-range hopping
Before justifying the mapping of the Anderson localization model with nearest neighbor hopping on RRG to the Rosenzweig-Porter random matrix model with infinite-range hopping we review the basic facts about RRG. The random K-regular graph (RRG) of N sites is a graph in which any site is connected to K + 1 other sites in a random way. The hopping amplitude V = 1 between nearest neighbor sites is fixed for all links. The Anderson localization model on RRG adds a random potential ε n fluctuating independently at any site n with the site-independent distribution P(ε) characterized by the variance ε 2 n ∼ W 2 . RRG is known to have locally a Cayley-tree structure with the branching number K. However, in contrast to a finite Cayley tree which is a hierarchical graph growing from a certain root onwards, any point of RRG can be considered as a root of a local Cayley tree. This is because RRG has loops (of which overwhelming majority are long loops with the length of the order of the diameter of the graph) which connect one local tree with the other thus making all points of the graph statistically equivalent.
The local tree structure and the predominance of long loops on RRG lead to the exponential growth of the distribution p N (r) ∼ K r /N of distances r between pairs of vertices on this graph. This growth lasts nearly up to the maximal distance on RRG (the diameter d ≈ ln N/ ln K), followed by an abrupt drop to zero. The most probable distance, r = r 0 , differs from the graph diameter d only by few units and for large graphs of N vertices they are approximately equal r 0 d (see the left panel of Fig. 6).
Moreover, due to the exponential growth of p N (r) with r, in the thermodynamic limit N → ∞ a finite fraction of pairs of sites on RRG is at the most probable distance, p N →∞ (r 0 ) → These "halo" sites are connected by multiple loops on RRG. The sites with small distances r d can be considered as a "medium" through which an effective hopping between the "halo" sites occurs in the high-order of perturbation theory in the nearest neighbor hopping on RRG. These sites form an almost loop-less, tree structure which allows to compute the distribution function of i.i.d. effective large-distance hopping matrix elements between the "halo" sites. We conjecture that the Rosenzweig-Porter random matrix model obtained by completing the random matrix of effective hopping over the "halo" sites to full matrix (complete graph) captures all the essential large-distance f > 0. This "condensation of large distances" is the crucial point for the mapping of the Anderson model on RRG onto a RP model. Indeed, let us take any vertex ("site") n on RRG and consider the set of vertices m(n) at distances r 0 − a < |n − m(n)| < r 0 (a r 0 ) around the most probable distance r 0 from it. We call them the "halo" sites. The other sites j at distances |n − j| < r 0 − a will be referred to as the "local tree" sites. Suppose that the local tree structure of RRG with the root at the site n holds strictly for all the sites j (i.e. the loops are absent). Then there is one single path from the site n to any site j. Thus one can compute the Green's function G n,j as a product of the (real) single-site Green's functions G l,l of a Cayley tree along the path from n to j. We assume that this property holds up to the distance |n − j| = r 0 − a (with a proper choice of a), so that the transmission amplitude from the site n to the "halo sites" m(n) can be computed within this "tree approximation". In this way one finds an effective hopping matrix element H nm = G nm from n to m(n). The matrix element H nm is a strongly fluctuating quantity of random sign. However, the probability distribution F(y = ln |H nm |) is identical for all n, m(n) because of the nearly fixed distance r nm ≈ r 0 ≈ d to any of the "halo" sites. It is also independent of the initial site n, since all the sites on RRG are statistically equivalent. So, we obtain a random matrix ensemble with the i.i.d. distribution of the offdiagonal matrix elements H nm = G nm fully determined by the local tree structure of RRG. However, the two "halo" sites m 1 (n 1 ) and m 2 (n 2 ) are not necessarily at the distance ∼ d from each other. This means that the random matrix H nm is not a full matrix representing a complete graph. Rather, it has a structure of a large-distance matrix on RRG (see middle panel of Fig. 6) with the matrix elements H nm between the "halo" sites shown in dark blue. The complete sub-graphs ("cliques") corresponding to this matrix have maximal sizes ∼ ln N N . Notwithstanding this we conjecture that the scaling with N of statistics of eigenfunctions and the stretch-exponential dynamics in so defined random matrix ensemble is identical to that in the Rosenzweig-Porter random matrix ensemble, in which all the off-diagonal matrix elements are i.i.d. with the distribution F(y = ln |H nm |). The reason for that is that even in typically full random matrix drawn from the RP ensemble these properties are determined by atypically large 3 hopping matrix elements (shown by a cyan strip in right panel of Fig. 6). Since the relevant large matrix elements in a RP matrix form anyway a sparse matrix, the sparseness of the matrix H nm = G nm does not matter. What matters, instead, is the broad, large-deviation character of the distribution F(y = ln |H nm |), Eq. (3) (see right panel of Fig. 6).
We would like to note that the replacement of the short-range hopping matrix on RRG by the long-range hopping matrix in the equivalent RP ensemble is not an exact transformation. Rather, it has the meaning similar to the replacement of the sparse matrix of hopping in the three-dimensional Anderson model (3D AM) by the full matrix of classical Wigner-Dyson (WD) ensemble describing well the level-and eigenfunction statistics in 3D AM at smallenergies/large times in the extended state. However, it does not apply to the "high-energy" properties like the spectral bandwidth E BW . The latter is independent of N in 3D AM and ∝ √ N in the Wigner-Dyson semi-circle. The correspondence between the two models is established after the level spacing is measured in units of the mean-level spacing δ = E BW /N and on the scale of energies smaller than the Thouless energy. Likewise, a proposed correspondence between the the Anderson model on RRG and the RP random matrix model requires that all the energies are measured in units of the spectral bandwidth and are much smaller than this scale.
Furthermore, in the case of long-range, semi-classical impurities in 3D metals one should distinguish between the Thouless energy (E T h ) below which WD level statistics is valid and the Ehrenfest energy (E Ehr ) above which perturbation theory works. In the interval of times 1/E Ehr < t < 1/E T h (or in the corresponding interval of energies) the dynamics is diffusive. In the case of RRG the role of the Ehrenfest energy is played by E 0 ∼ E BW , while the role of the Thouless energy is played by In the interval of times E −1 Ehr t E −1 T h the survival probability R av (t), Eq. (5), is stretchexponential (see Eqs. (40), (41)).
As the diffusion on the Cayley tree results in the exponentially decreasing survival probability, the stretch-exponential one may be considered as a sub-diffusion.

"Multifractal" distribution of hopping
Next, we consider generic properties of the distribution function F r (y = ln |G nj |) of the twopoint Green's functions G nj on an infinite Cayley tree.
The crucial point here is that in the absence of loops the two-point cavity Green's function G nj,l→k (short-hand notation G nj ) can be expressed in terms of the product of one-point cavity Green's functions G p→p (short-hand notation G p ) along a (unique) path of length r = |n − j| that connects the points n and j: In order to proceed further we use the result of Ref. [13] where it is shown that the moments I β = G 2β are expressed in terms of the eigenvalue β of the linearized transfer-matrix operator (analytically continued from the segment β = [0, 1] to entire complex plane β ∈ C), first introduced in the seminal work of Abou-Chacra, Thouless and Anderson (ACTA) [58]: Then the Mellin transform yields: where the integration is performed over the Bromwich β = c + iz contour which goes parallel to the imaginary axis ( Im β ∈ [−∞, +∞]) on the positive side of the real one (c > 0). The eigenvalue β obeys a celebrated ACTA symmetry [13,58]: and the identities: As a result [13], in the limit of long path r 1 the distribution F r (y = ln |G nj |) of the product Eq. (64) has a special symmetry: which coincides with Eq. (62). At large r one can compute the Mellin transform in the saddle-point approximation: where β * is found from the stationarity condition: Eq. (71) implies that β * is a function of y/r. It follows from Eq. (70) that where F (g) is some independent of r (in the limit of large r) function (see Fig. 7) which precise form depends on β (for the dependence of β on disorder for K = 2 Cayley tree see Appendix B) and which is obtained from the following equations: The distribution F(y = ln |H nm |) for the RP model equivalent to RRG is obtained from Eq. (72) by setting r = d ≈ ln N/ ln K: Eqs. (70), (71) and (73) define the distribution of y = ln |H nm | of the form Eq. (3). This is a broad distribution which characteristic width grows with ln N and its typical value decreases as a negative power of N . This form is known as the large deviation, or multifractal ansatz. It appears in many different problems of statistical mechanics (see e.g. Ref. [64] and references therein) and is a non-trivial generalization of Central Limit Theorem when the logarithm of the fluctuating quantity is a sum of many terms with special correlations between them. An example of such a distribution is shown on the right panel of Fig Let us now establish some properties of the function F (g) and some particular values of g β that follow from Eqs. (62), (69). Differentiating Eq. (62) w.r.t. g one obtains: Substituting g = g 1 in this relation we easily obtain: If the function F (g) is monotonous and given that F (g 0 ) = 0 one obtains: Now, substituting g 1 for g in Eq. (62) and using Eq. (79) we obtain: Here we also used the normalization condition of the PDF Eqs. (3),(75) that requires: Now let us plug g 1/2 into Eq. (77) and take into account the definition F (g 1/2 ) = 1/2. Then: Now, again, assuming the monotonic behavior of F (g) one obtains: Next, from Eq. (83),Eq. (73) and Eq. (76) we obtain: Another useful relationship one can obtain from Eqs. (74),(13): In particular:

The Lyapunov exponents
In this section we express the fundamental quantities, the Lyapunov exponents, on an infinite Cayley tree in terms of the function F (g) and the eigenvalue β . The Lyapunov exponents are defined as: and where G(r) = G nj at |n − j| = r. The typical Lyapunov exponent, Eq. (87) is given by the typical value y typ /r = (g 0 /2) ln K of F r (y): The "average" Lyapunov exponent, Eq. (88) is expressed in terms of F (0) as follows: To arrive at this result we used the saddle point approximation in computing the average |G(r)| over the distribution F r (y), Eq. (72), and Eqs. (83),(84). Eq. (90) allows to establish the limits of variation of F (0) in the extended phase on the Cayley tree. Indeed, the localization transition point in an infinite tree is determined by [58]: Thus the lower limit of F (0) in the extended phase is F (0) = −1. On the other hand, in the case of single orbital per site in the clean limit both Lyapunov exponents tend to (1/2) ln K (see, e.g. [13]). This corresponds to the upper limit F (0) = −1/2. Note that for a granular tree with n → ∞ orbitals per site (described by a sigma-model) the limit F (0) = −1/2 is reached at a finite value of inter-granula conductance. This happens at a point of a transition to an ergodic phase [20] in which F (0) remains equal to −1/2. We conclude therefore that for any model on a tree in the extended phase:

Finite-size multifractality in RRG
Now we are in the position to show the presence of the finite-size multifractality in the MF-RP model with ACTA symmetry (62). Indeed, computing in the saddle-point approximation the averages (9), (10), (11), and (12) with the distribution Eq. (3) and using the properties of the function F (g) established in Sec. 7.3 we obtain the following criteria: • Localization transition: • Ergodic transition: • FWE transition: One can see from Eqs. (93),(94) that in the thermodynamic limit N → ∞ the criteria for the localization and the ergodic transitions become identical, i.e. the transitions merge together at the tricritical point (see Fig. 1). The absence of a gap between them implies that the multifractal phase collapses in this limit. We would like to emphasize that this conclusion is heavily based on the ACTA symmetry, Eq. (62), since it is this symmetry that ensures g 1/2 = 0 and g 1 < 0. For a distribution, Eq. (3), that does not respect this symmetry (e.g. for the logarithmically-normal one, Eq. (47), with p = 1) the true multifractal phase does exist if F (g 1/2 ) > −1 but F (0) < −1 (e.g. for p < 1 in Fig. 1).
However, at any finite ln N the finite-size multifractality appears even if the ACTA symmetry is respected: Note that the leading double-logarithmic finite-size correction is present for the ergodic transition and is absent for the localization one. This happens because the integral that determines the average in Eq. (9) is dominated by the saddle-point g = g 1/2 = 0, while the saddle point g = g 1 < 0 in the average in Eq. (10) is outside the region of integration. This average is dominated by the lower cut-off g = 0 where F (g) = 1. Therefore, the measure of the saddlepoint integration ∼ (ln N ) −1/2 in Eq. (9) is canceled out by the corresponding factor in the normalization constant C, Eq. (12), while the measure ∼ (ln N ) −1 of the end-point integration in Eq. (10) is not. One can see from Eq. (96) and Fig. 2 that the leading finite-size correction (−1/2) ln ln N/ ln N makes W AT higher than W ET . The next correction appears to be also of the same sign at K = 2. As the result a finite-size multifractality appears inside a gap between W AT and W ET which is shown in the inset of Fig. 1. At conventional sizes N = 32000 accessible for exact diagonalization the finite-size multifractality may extend itself as far as to W = 14 down from the infinite-size localization transition W c = 18.2 at K = 2 (see Fig. 2). The emergence of finite-size multifractality for system sizes smaller than the one given by Eq. (96) is of principal importance. For very large ln N one can expand F (0) = ln 1/2 / ln K in W around the infinite-size tri-critical point W c and define the critical length L MF = ln N MF which is associated with finite-size multifractality, (8): For K = 2 the coefficient k ≈ 0.6, c ≈ 1.9.
Note that this length is principally different from the length L c ∼ ln N c ∼ (W c − W ) −ν with ν = 1/2 which emerged in the sigma-model description of the granular Cayley tree and RRG in Refs. [54,55,65] and in many other works thereafter (e.g. in the recent work [66]). The critical length L c with the exponent ν = 1/2, is indeed, present in the finite-size scaling (FSS) in localization problem on RRG. It is also present in FSS in an RP model with the ACTA symmetry, e.g. in the logarithmically-normal RP ensemble with p = 1 [47]. The critical length L MF that controls finite-size multifractality was earlier suggested [56] but it was believed that it is the same length as L c . However, recently a numerical evidence appeared that the additional critical length with the exponent ν = 1 is also present in problem considered [13,67,68].
Eq. (8) that follows from Eq. (97) finally solves this problem and provides a physical meaning for the additional critical length L MF that is much larger than L c .
8 Results of exact diagonalization on RRG and on LN-RP ensemble.
In order to confirm the correspondence between RRG and MF-RP developed in the previous section, here we present the results of exact diagonalization for R av (t) both in the logarithmically-normal RP ensemble, Eq. decreases with increasing N approximately as κ(N ) = κ ∞ + a/ ln N . This extrapolation gives κ ∞ = 0.23 which is in a reasonably good agreement with the 'theoretical value' κ theor = 0.193, see Appendix B. However, the convergence to the thermodynamic limit is almost reached for RRG at W = 8 at system sizes N = 32768, 65536 (see Fig. 9 and Fig. 10 in Appendix A).
In this case the stretch-exponent κ(N = 65536) ≈ 0.50, which should be compared with the theoretical value κ theor ≈ 0.473 (see Appendix A and B). In Appendix A we also present the results for RRG at W = 6 where the convergence to N → ∞ limit is good already at the system size N = 32768. In this case we obtain κ(N = 32768) = 0.64 vs. the theoretical value κ theor = 0.689. These results convincingly demonstrate that the mapping of the Anderson model on RRG onto the multifractal RP ensemble is quantitatively correct.

Conclusions and Discussion
In summary, in this paper we introduce a generic set of dense Rosenzweig-Porter randommatrix models with fat-tailed distributions of the off-diagonal elements. The consideration of such models is motivated by the fact that similar models emerge in the mapping of the sparse matrix Hamiltonians of disordered interacting many-body systems [17] and of the Anderson localization model on hierarchical graphs [35] onto the dense random matrix models. The above mapping is characterized by the large deviation ("multifractal") character of the distribution of the off-diagonal matrix elements with the corresponding dependence of the typical value and the width of the distribution on the system size.
We focus on the dynamical properties of such random-matrix models and identify analytically three distinctly different dynamical regimes: (ii) Sub-diffusion, and (iii) "Frozen" dynamics.
Since the above mentioned mapping of models with short-range coupling to the infinite-range RP ones invalidates the notion of distance, in order to describe the dynamics we use, instead of the wave packet spreading, the probability R(t) of a particle initially prepared in a certain graph node (representing a configuration in the Hilbert space) to return back to it.
In the single-particle models on d-dimensional lattices, such return probability decays as the inverse (sub)diffusion volume occupied by a wave packet at a time t. Instead, the Hilbert space of the many-body systems and the Anderson localization model on the hierarchical graphs correspond to the infinite dimensionality and show the exponential growth of the number of nodes with the Hamming distance. This immediately implies the exponential decay of the return probability in the diffusive phase R(t) ∼ e −Dt . The corresponding sub-diffusive dynamics should then be given by the stretch-exponential decay We have shown that the sub-diffusive, stretch-exponential dynamics of R(t) emerges quite generically in the ergodic phase of the "multifractal" RP model in the proximity to the Anderson localization transition. In the toy model of a quantum particle on disordered Random Regular Graph, the sub-diffusive regime spreads throughout the entire ergodic phase down to vanishing disorder. The mechanism of sub-diffusion is related to the rare large couplings which contribution to the decay rate decreases with time. This is essentially a kind of a quantum Griffiths effect which leads to the failure of the Fermi Golden Rule and as such it goes in line with the explanation [36,45] of slow dynamics in certain quantum many-body systems [69,70]. As a next step in this direction, one can consider the mapping of a true many-body system onto the Rosenzweig-Porter model of the corresponding symmetry, similarly to an oversimplified approach used in Ref. [17] to map several many-body problems onto the Gaussian Rosenzweig-Porter ensemble. This challenging problem consists of two parts. One should first take into account the essential graph structure of the Hilbert space of a realistic manybody system, eliminating system-specific details (e.g. the weak links as in the percolation problem [71,72]). Second, the disordered potential which is uncorrelated on different sites in the coordinate space in the interacting systems, imposes correlations and constraints in the corresponding disorder distribution in the Hilbert space [34]. The effects of the correlated resonances and the corresponding avalanche structure of delocalization transition in such systems may drastically affect the dynamical phase diagram.
Another relevant direction is to focus on the (sub)diffusive dynamics in the coordinate space of the interacting systems instead of the Hilbert space. In order to do this in terms of the local observable like R(t) one should consider the probability of a system to return not to a certain Hilbert space configuration but to a fixed state of a certain local operator (e.g., the occupation number on one site). This fixed state corresponds to an entangled manifold of quantum configurations and should be therefore described by a density matrix. The corresponding generalization of the Wigner-Weisskopf approximation for such a case is the subject of the further investigations.
A Numerics on LN-RP model and on Anderson localization model on RRG.
In this Appendix we present the results of the exact diagonalization of 30% to 100% of states in the LN-RP model defined by Eq. (47) and compare these results with the corresponding results for the Anderson model on RRG with K = 2.
We start by studying the functions K(ω) and µ(ω), Eqs. (54) and (56), which definitions and principle properties were discussed in Section 6. In Fig. 9 we plot the results of the exact diagonalization for the functions K(ω) and µ(ω) for LN-RP with p = 1 respecting the symmetry, Eq. (62) and for RRG with W = 12 and W = 8 at different system sizes up to N = 131072. One can see that there are two different regions of ω: the "low-energy" region, where K(ω) and µ(ω) exhibit similar universal behavior for RRG and the "symmetric" LN-RP model, and the "high-energy" region where the behavior is system-specific and independent of the system size. The energy scale separating the high-and the low-energy parts is the energy scale of establishing of the universal "RP regime". It does not depend on N but depends on disorder strength W . In some vague sense it is similar to the inverse mean free path −1 in the problem of diffusion in dirty metals. For distances r < the motion of particle is ballistic while for r > a universal diffusive regime is formed. For describing the localization effects one has to take into account interaction between diffusion propagators (of particle-hole and particle-particle type) which is conveniently done in the framework of the nonlinear sigmamodel. For such a coarse-grained model, the length is the low-distance cutoff. Analogously, the "coarse-graining" in our problem is done by mapping of the full Anderson model on RRG onto an effective Rosenzweig-Porter model.
Let us focus at the "low-energy" part of K(ω) and µ(ω). One can see that at small sizes the curves "move" to the left with increasing the system size N signaling on the N -dependence of the characteristic energy scale E 0 ∼ N −a , (a > 0). In general, such a behavior is a signature of multifractality, including a finite-size multifractality.
In the case of a true MF phase in LN-RP model with p = 0.5 without the symmetry, Eq. (62), shown on the lower left panel of Fig. 9, the "speed" of this motion approaches a finite limiting value, which is demonstrated in the lower inset of that panel. At large enough N the position of the inflection point ln ω inf in µ(ω) moves to the left with a "speed" d ln ω inf /d ln N that is extrapolated to N → ∞ linearly in 1/ ln N to be τ * = −0.4. This value is reasonably close to the predicted value τ * = −0.35 (see Eq.(100) of Appendix C).
However, in the case of the finite-size multifractality in the weakly-ergodic phase, as the system size N increases, the "speed" of this "motion" slows down to zero (see upper-right panel of Fig. 9) and may even reverse the sign (upper-left panel of Fig. 9). On RRG for small enough disorder (e.g. for W = 8 in upper-right panel of Fig. 9) the "motion" stops at accessible system sizes showing convergence to some limiting curve in the thermodynamic limit N → ∞. For LN-RP in the weakly-ergodic phase (see upper left panel of Fig. 9) the sign of a(N ) changes at some N and becomes positive for large enough N . This is in agreement with E 0 ∼ E BW → N −τ * with τ * < 0 (see Eq.(100) in Appendix C) in LN-RP model and E 0 ∼ E BW ∼ N 0 on RRG. Now let us consider the shape of the curves µ(ω) in the universal "low-energy" regime. Again, there is a qualitative difference between µ(ω) in the weakly-ergodic phase shown in the right and upper panels of Fig. 9 and that in the multifractal phase shown on the lower left panel of that figure. In weakly ergodic phase all the curves for µ(ω) on the main plot show a minimum, while in the multifractal phase there is only a bending point on the corresponding level µ ≈ −1.
A similar minimum in µ(ω) can be extracted from the population dynamics numerics presented in Fig.10 of Ref. [13] (see Fig. 11).
In order to get the correct behavior of survival/return probability R av (t) we Fouriertransformed K(ω) in the entire "low-energy" region of ω (see Fig. 9) and checked the results by the direct calculation of R av (t) in the time domain (see lower panel of Fig. 8). The results are shown in the upper panels of Fig. 8 in the main text and in Fig. 10. All those plots are very similar, no matter in the weakly-ergodic or in the multifractal phase. This implies that in the non-symmetric case of LN-RP when both multifractal and weakly-ergodic phases exist, the stretch-exponential decay of survival/return probability is not sensitive (apart from changing the characteristic scale E 0 ) to the ergodic transition (e.g., when crossing the green line on the phase diagram of Fig. 1 ). This is an argument in favor of the common origin and hidden similarity of both these phases.
B Theoretical values of the stretch-exponent κ for K = 2 RRG.
The "theoretical" value of the stretch-exponent κ theor on RRG can be computed with any prescribed accuracy using the exact theory of localization on a Cayley tree of Ref. [58]. A reformulated variant of this theory and details of calculations of β will be published elsewhere [73]. Here we note that in these calculations we solved a non-linear equation for the modified effective distribution of on-site disorder prior to solve the eigenvalue problem for the linearized transfer-matrix equation with this modified distribution. The results for the eigenvalue β and the function F (g) in the "multifractal" distribution, Eq. (3), of hopping matrix elements is shown in Fig. 12 for several values of disorder. The "theoretical" value of the stretch-exponent κ theor is found from Eqs. (36) and (39) and is determined by the derivative F (g) at a point where F (g) = −1, as shown on Fig. 12. We obtained the following values of κ theor : The full curve κ theor vs W for RRG with K = 2 and box distribution of disorder is shown in Fig. 13. The inset on this figure demonstrates that the stretch-exponential phase persists up to vanishing disorder, like the multifractal phase on the finite Cayley tree with one orbital per site [13,22]. This implies that there is an intimate relation between the two phases, with the SE phase on RRG being a "remnant" of the MF phase on the finite Cayley tree.
C Theoretical values of exponents κ, τ * and ∆ for generic LN-RP model As it follows from Eq.(100) the characteristic energy scale E 0 ∼ N −τ * decreases with increasing N in the multifractal phases (in this case its physical meaning is the width of a mini-band in the local spectrum) and increases with N in the ergodic phases (in this case it is of the order of the total spectral bend-width E BW ). The band-width E BW ∼ N 1−∆ is blowing up with increasing N in the ergodic phases, while it is constant in the non-ergodic ones (see Eq.(101)).  47), with the ACTA symmetry at p = 1 and γ = 3.6 for different matrix sizes N = 512 − 65536. The running power µ(ω) interpolates between zero at low frequencies and µ ≈ −1 at the minimum. The dependence on the system size shows a reentrant behavior with the direction of evolutions shown by arrows. It reflects the reduction of the characteristic energy scale E 0 at small system sizes due to the finite-size multifractality and the growth of E 0 ∼ E BW with N at large system sizes due to N -dependence of the total spectral bandwidth E BW (for RRG E BW ≈ W/2 is size-independent). (Upper-Right panel): µ(ω) in the Anderson model on RRG with W = 8. Lower inset: the evolution of the inflection point on the main plot with increasing the system size. The main plot and this evolution show convergence to the N = ∞ limit. The "low-energy" part ln ω < 2 is similar to that for LN-RP model on the upper-left panel. In the limit N → ∞ we expect both curves to be identical in their low-energy parts if time is measured in units of the corresponding E −1 BW . (Lower-Left panel) µ(ω) for LN-RP model (p = 0.5, γ = 2.5) without the ACTA symmetry in its multifractal phase for N = 1024 − 131072. The curve for µ(ω) do not show a minimum as a function of ln ω but rather a bending point at the level of −1. Lower inset: the evolution with N of the inflection point on the main plot. It approaches the power-law E 0 ∼ N −0.4 in the limit of very large N . The characteristic scale E 0 decreases in the multifractal phase due to the reduction of the width of a mini-band. (Lower-Right panel): µ(ω) for RRG at W = 12. The "low energy" part ln ω −4 is similar to that in 'symmetric' LN-RP on the upper-left panel, though the convergence to the thermodynamic limit N → ∞ is much slower and does not show re-entrance at system sizes studied. The "high-energy" part at ln ω −4 is system specific and N -independent. The low-energy part ln ω < −4 was Fourier-transformed to obtain the stretch-exponential behavior on the upper-right panel of Fig. 8. (Upper insets) in all panels show K(ω) for the same parameters as the main panels. The color code (from purple to black) of the curves for system sizes N = 512, 1024, 2048, 4096, 8192, 16384, 32768, 65536, 131072 is the same on all the plots and given by the legend in the upper left panel.  N = 131072). It was obtained by the numerical Fourier-transform of the "low-energy" part ln ω < −1.5 of K(ω) in the lower-left panel of Fig. 9. The stretch-exponential behavior R av (t) ∼ exp[−(tE 0 ) κ ] is qualitatively similar to the one obtained in the weakly-ergodic phase of LN-RP model (upper-left panel of Fig. 8), despite the qualitative difference in behavior of µ(ω) for ω above the onset of a shoulder (see left panels of Fig. 9). This demonstrates that only ω below the onset of the shoulder is relevant for the Fourier transform in the time-interval of the stretch-exponential behavior, where the Levy α-stable distribution emerges in K(ω).