Emergent fractal phase in energy stratified random models

We study the effects of partial correlations in kinetic hopping terms of long-range disordered random matrix models on their localization properties. We consider a set of models interpolating between fully-localized Richardson's model and the celebrated Rosenzweig-Porter model (with implemented translation-invariant symmetry). In order to do this, we propose the energy-stratified spectral structure of the hopping term allowing one to decrease the range of correlations gradually. We show both analytically and numerically that any deviation from the completely correlated case leads to the emergent non-ergodic delocalization in the system unlike the predictions of localization of cooperative shielding. In order to describe the models with correlated kinetic terms, we develop the generalization of the Dyson Brownian motion and cavity approaches basing on stochastic matrix process with independent rank-one matrix increments and examine its applicability to the above set of models.


Introduction
Fractality and multifractality are intriguing and rich phenomena quite widely spread in real world, including coastline profiles, turbulence or even heartbeat and financial dynamics. In quantum systems, multifractal wave functions sharing the properties of metallic (ergodic) and insulating phases open the gap for the ergodicity breaking without a complete localization. First found at the Anderson localization transition [1,2] in single-particle disordered systems, such non-ergodic extended states have been recently claimed to describe wave-function properties of a whole many-body localized (MBL) phase in the Hilbert space of strongly interacting disordered quantum systems [3][4][5][6][7][8]. The complexity of the above systems makes the numerical simulations of them to be difficult. Indeed, the description of the MBL transition is a highly demanding problem which shows controversial numerical results not only on the critical disorder value [9][10][11], but on the existence of the MBL phase itself [11][12][13].
Therefore it is of particular concern and high demand to model essential attributes of such systems in less complicated single-particle or random matrix settings. Here, in order to simulate the multifractal eigenstate structure of the MBL phase one has to realize the random matrix models with the robust multifractal phase(s) in them.
Recently the interest to such long-range models and robust multifractality in them has been stirred up by their relevance and an immediate applicability of their multifractal paradigm to the description of black-hole physics based on the Sachdev-Ye-Kitaev model with diagonal disorder [8,29] and to the speeding up the algorithms of quantum annealing [30][31][32] and teaching of weak learners in artificial neural networks [33].
Nevertheless, the disordered long-range models undergoing Anderson localization transition are not just purely idealized mathematical objects living in the brains of scientists. By contrast, being suggested in 1980s [2,[34][35][36], they have been realized in the variety of physical systems spreading from trapped ions [37,38], ultra cold polar molecules [39], magnetic [40,41], and Rydberg [42,43] atoms to nitrogen vacancies in diamond [44], photonic crystals [45], nuclear spins [46], and Frenkel excitations [47]. All these realizations confirm that the presence of the long-range hopping matrix elements indeed may delocalize the system and restore the localization transition even in low-dimensional systems where the short-range models are known to be completely localized [48].
In addition, the disordered many-body systems considered in the Hilbert space [6] and their counterparts on the hierarchical graph structures [?,22,23,49] have been recently mapped (within some approximations) to the above Rosenzweig-Porter-like models. In these mappings the all-to-all hops which amplitudes depend on the Hilbert space dimension N emerge naturally from the short-range models in the hierarchical (Hilbert) space. Indeed, the all-to-all coupling is obtained in these models from the consideration of long series of quantum transitions to (at least) exponentially growing number of available configurations at the certain hopping distance, while the size-dependence of the hopping amplitudes is related to the exponentially decaying propagators [6]. Proliferating this process up to the Hilbert space diameter, one immediately restores the complete graph with the N -dependent amplitudes [49].
However, with all the above mentioned success in applications [6,8,22,23,[29][30][31][32][33]49], the RP model has some caveats. Indeed, it is formed by a complete graph with the independent identically distributed (i.i.d.) hops, while in the models mapped to the short-range manybody disordered models the hopping terms in the Hilbert space are strongly and non-trivially correlated due to the dominance of far-away resonances in the corresponding long series of quantum transitions [?].
Taking account of such correlations brings another surprise into play. Indeed, in the case of completely correlated hopping amplitudes (see, e.g., disordered Richardson's model [50,51] and Burin-Maksimov model [52,53]) even the long-range character of hopping terms cannot delocalize the majority of the states. These effects called in the literature the cooperative shielding [54,55] and the correlation-induced localization [21,[56][57][58][59][60] are based on the observation that in the disorder-free versions of such systems due to the correlated nature of the kinetic long-range terms there is measure zero of high-energy states that take the most spectral weight of the hopping and effectively screen the bulk states from the off-diagonal matrix elements. The coexistence of few high-energy states with the nearly degenerate bulk states forms a kind of energy stratification, when measure zero of states are separated from each other and from the rest of the spectrum. In the disordered setting, the above screening effectively returns the majority of the eigenstates to the short-range disordered situation with the localization described by the previously known results [48]. The high energy of a few stratified spectral edge states prevent them from being localized.
Here, a question immediately arises: to what extent such a correlation-induced localization is fragile to the reduction of the correlations? On one hand, the answer to this question has been partially given in [56,57,60,61]: (i) the localization properties of the systems with the correlations respecting only the translation-invariant (TI) symmetry are nearly indistinguishable from the ones in the uncorrelated systems [57]; (ii) the addition of a tiny fraction of delocalized uncorrelated hopping to completely correlated models breaks down the localization in them [61], while (iii) the correlation-induced localization survives under the addition of finite staggered potential [56] or anisotropy-induced quasi-disorder [60] terms. But on the other hand, the partial correlations (not just mixture of fully-correlated and fully-uncorrelated terms) were out of scope of these papers.
Therefore in this paper we would like to combine the above mentioned energy stratification of measure zero of levels in the disorder-free case, common for both the correlation-induced localization and cooperative shielding, with the partial correlations in the hopping terms. As a remarkable example we consider the full set of models interpolating between the fully correlated Richardson's model (RM) [50,51] and the Rosenzweig-Porter with translation-invariant correlations (TIRP) [57].
Our choice of the TIRP model which has the same wave-function statistics as the conventional uncorrelated RP model is motivated by the simple (plain wave) eigenbasis of the kinetic term allowing one straightforwardly interpolate between the above models keeping this hopping term basis fixed. On the other hand, although the wave-function phase diagram of TIRP is the same as the one of RP, there is no rigorous method to describe this model. Indeed, both the Dyson Brownian motion (DBM) [62] and the cavity method [63] working for fully uncorrelated models like RP [15,24] are not applicable to the set of models of our interests. Therefore in this paper we develop another method based on the spectral decomposition of the kinetic term and the Sherman-Morrison formula for the Green's functions, generalizing the above methods, and benchmark this powerful tool with respect to the numerical results.

Model
We consider an ensemble of random N × N Hamiltonians H = H 0 + V with the on-site disorder given by the matrix H 0 , diagonal in the coordinate basis of sites |i , having i.i.d. random entries ε i with zero mean and unit variance on the diagonal and perturbed by a hopping matrix V . This hopping matrix V is assumed to be infinitely long-range with a certain eigenbasis {|g k }. Further for simplicity we consider {|g k } to be plain waves 1 with an integer wave number 0 ≤ g k < N to the state's index k correspondence chosen randomly for each realization, i.e.
In order to describe the RM and TIRP (as well as the whole class of models interpolating between them) we consider the hopping matrix V to have the rank R = N β , i.e., the number of non-zero real eigenvalues v k = 0, see Figs. 1(d-f). Then RM is characterized by the rank R = 1 matrix V (β = 0) when all the hopping elements are completely correlated to each other (see Fig. 1(a)): e.g., for g 1 = 0, i.e. i|g 1 = N −1/2 = const while TIRP model corresponds to R = N (β = 1) with N i.i.d. Gaussian random numbers v k providing as well N i.i.d. Gaussian random TI-coefficients Fig. 1(c) 1 All the results are not sensitive to this choice: the only important assumption is that the coefficients i|g k ∼ 1/N 1/2 are of the same order and, thus, ergodic in the coordinate basis.

Richardson ( = 0)
Translation-invariant Rosenzweig-Porter ( = 1) Intermediate models with stratified levels The phase diagram of both above models is determined by the scaling of V i,j ∼ N −γ/2 with a certain parameter γ taken from the consideration of RP model [14]. The standard locator expansion converges at V i,j ∼ N −1 , γ = 2, providing the localization of all eigenstates at γ > 2. However the difference in correlations (2) and (3) makes the phase diagrams of RM and TIRP completely different at γ < 2, see the top and bottom lines in Fig. 5.
In the above parametrization when V ij ∼ N −γ/2 for a general value of R = N β the i.i.d. Gaussian random variables v k have zero mean and the fixed variance given by Therefore, like in RM with β = 0, Fig. 1(d), in a more general case with 0 < β < 2 − γ one realizes the energy stratification in the hopping term V where a zero fraction R/N ∼ N −(1−β) of states has an enhanced energies v k ∼ N 1−(γ+β)/2 , Fig. 1(e), while the rest of the states are degenerate in the disorder-free setting H 0 = 0. The presence of R N non-zero v k provides the correlation between TI-hopping terms Fig. 1(b). Indeed, where . . . {v k } stands for the average over the hopping eigenvalues {v k } and the function f (x) for each fixed choice of {g k } ↔ {k} correspondence is given by the sum jumping from f (0) = 1 to f (x = 0) ∼ R −1/2 . Of course for different sets of {g k } the correlations (5) given by f (x = 0) have different signs. 2 A set of such Hamiltonians with different values of β smoothly connects RM (β = 0) with the TIRP model (β = 1), i.e. interpolates between a maximally structured translationinvariant model to a TI model with a maximally unstructured hopping while preserving an amplitude of the off-diagonal entries. This is of course not the only way to go between these two models but definitely the interesting one since by a similar interpolation one can in principle, rank by rank, reach any model, not only TIRP one.
In this work, we focus on the calculation of the fractal dimension D(β) defined via the inverse participation ratio of the eigenstates |ψ E in the spectral bulk for the model (1) with 0 < β < 1 as 3 From the previous studies (see, e.g., [57]) the fractal dimensions for the limiting cases of β = 0 and β = 1 are known. Indeed, as for β = 0 and γ < 2 the eigenstates are critically localized, one should have D(0) = 0, while for β = 1 and γ < 2, In the intermediate regime, 0 < β < 1, following the paradigm of the cooperative shielding [54] one expects to have zero fractal dimension D(β) = 0 as soon as there are energy stratified eigenstates with the extensive level spacing 2 Formally, on average over the sets of as in the TIRP model. Moreover, this works at any β including the Richardson's model (β = 0) due to the randomness of the basis vectors |g k . However, as we will see below the presence of the higher-order correlations like V∆iV∆jV * ∆k V * ∆i+∆j−∆k interpolates the phase diagram between TIRP and Richardson's models. The washing out of the two-point correlations in this case are not important for the further consideration. 3 Since TIRP doesn't possesses multifractality [57], we do not consider Dq for different q. Another definition of this quantity which we actually use will be given via self-energy in Sec. 5.4. i.e., at least for β < (2 − γ)/3. These N β levels are detached from the bulk of the spectrum, while the typical states in the spectral bulk are shielded from the long-ranged hopping by these high-energy modes and stay virtually unchanged after an orthogonalization to the detached ones. However, surprisingly, this guess fails. Indeed, as we show further analytically and numerically, Fig. 2, the fractal dimension is given by the simple formula and shows the immediate emergence of eigenstate fractality of mid-spectrum states as soon as the number N β of energy stratified states becomes extensive, β > 0 (unlike the case of the Burin-Maksimov model [56][57][58]). This is the main result of our paper. Such a fractality emergence shows the fragile nature of the correlation-induced localization to partial reduction of correlations like in the case of mixture [61] of completely correlated and completely uncorrelated hopping terms. Before going to the technical calculations, let us give qualitative arguments in favor of the above result. As one can see from the Fig. 2, the model has two regimes which differ by the ratio of the bandwidths of diagonal σ(H 0 ) ∼ N 0 and hopping σ(V ) ∼ N 1−γ/2−β/2 matrix terms.
1. In the case of σ(H 0 ) σ(V ), β > 2 − γ, there is no energy-stratified states and the matrix V can be considered as a perturbation in the sense of the Fermi's golden rule. The corresponding fractal dimension can be estimated similarly to the usual RP case as D = 2 − γ. Though, one cannot take this estimate on faith without a reservation in a generic model as the hopping matrix V is correlated.
2. Unlike this, as soon as σ(H 0 ) σ(V ), β < 2−γ, almost all N β finite-energy eigenvectors of V have extensive energies v k N 0 and form a high-energy subspace S = span({|g k }) of states which cannot be significantly hybridized by the presence of a "small" matrix elements of H 0 . However, this H 0 matrix does hybridize the rest N − N β degenerate eigenstates of V and push them forward towards the localization. Though, since these hybridized eigenstates should be orthogonal to the high-energy subspace S and also orthonormal between each other, their typical support set is of the order of dim(S) = N β giving the desired linear dependence D = β 4 .
In order to derive Eq. (11) mathematically, one has to generalize the known methods (such as the cavity method of the Dyson Brownian motion) to the case of the correlated hopping terms (5). This we will do in the next sections.

Local resolvent as observable of interest
Before going to a newly developed method, we would like to describe an observable of our interest which can be used to extract the information about the eigenstate fractal dimensions (8). Like in the paper for the RP model [15], we would like to focus on the statistics of the local resolvent G(z) = (z − H) −1 with z = ε − iη. Following [15], in order to probe the non-ergodic properties of the system we consider η = N α δ ∼ 1/N 1−α , α > 0, to be much larger than the mean level spacing δ determined by the diagonal disorder 5 . In this case we are able to detect an arbitrary small non-zero fractal dimension by choosing a finite positive constant α < D(β). This works when the global resolvent Tr [G(z)] /N (the imaginary part of which is just the global density of states ρ ∼ (N δ) −1 ) is a featureless function converging to an order-one number in the thermodynamic limit N → ∞. 6 For the opposite case, please see in Appendix B.
The local density of states (LDOS) given by the imaginary part of the resolvent's diagonal element G ii (z) has the form where ψ n (i) = i|ψ n is the eigenstate's amplitude on site i; the sum represents a weighted average of |ψ n (i)| 2 over the energy window (ε − η, ε + η) and can be linked to the spatial structure of the eigenstates in this window. Now, consider an ensemble of Hamiltonians H with different hopping matrices V but fixed H 0 . If we assume the close in energy eigenstates from different realizations to live on almost the same fractal, see the bottom panel of Fig. 3, we can extract the information about the broadening Γ of the Lorenzian, determining the fractal  13), show more or less the same Lorenz-like profile with the width Γ provided η Γ and the close-in-energy eigenstates live on close fractals. dimension (8), from the resolvent averaged over the hopping disorder (please see Eqs. (30,31) for more details) This last assumption looks reasonable when the energy landscape of H 0 plays a dominant role in forming of the eigenstates profile, e.g., at least in case with all hopping elements of roughly the same amplitude. Hence, by choosing G(z) as our observable, we are restricting ourselves to such models.

Previously known methods
In this section we consider a couple of previously known methods and explain why they are not applicable for the models of our interest.
The cavity method is based on the block-matrix inversion formula for the diagonal part of the local resolvent written in terms of the resolvent G (i) jk (z) of the Hamiltonian with the ith row and column being removed (see Fig. 4(a)). Both for the tree-like graphs [64,65] and RP-like models [15,24] in the thermodynamic limit N → ∞ it is known that the contribution of off-diagonal elements G This method works well for the dense matrices with the uncorrelated entries V ij [15,19]. Indeed, in this case G (i) jk (z) is independent from V ij and V ki for any j, k and the sum in r.h.s. can be tackled using the central limit theorem. For N → ∞ the local resolvent averaged over V ij takes the form However, in the case of the correlated V ij (5) it is hard to proceed beyond Eq. (15). Another method used for the description of the RP model is the Dyson Brownian motion approach [15]. It is based on the stochastic process of adding random Gaussian matrices dV (t) to the special diagonal H 0 (see Fig. 4 This method allows to calculate the instant eigenvectors and eigenvalues of the problem from the stochastic differential equations and even write a closed equation for the local resolvent G ii (z) averaged over the uncorrelated Gaussian off-diagonal elements V ij [15] Like the cavity method, the Dyson Brownian motion approach is limited by the statistical independence of the additions dV ij (t) to the matrix elements both at different time instants and different coordinates (i, j).
5 The rank-one increment method

The main idea
In order to overcome the problem with the correlated character of the hopping matrix entries V ij we consider the matrix V as a sum of uncorrelated rank-one matrices. For example, we can use a natural eigenbasis decomposition (1b). Then we get a series of Hamiltonians with H 0 = H 0 (V 0 = 0) and H R = H. Exploiting the so-called Sherman-Morrison formula, one obtains a corresponding series of resolvents We use this exact recursion to obtain a differential equation for the resolvent G(z; r) averaged over the off-diagonal disorder. Note that in the case of uncorrelated hopping entries this equation resembles the aforementioned Dyson Brownian motion approach result. The graphical representation of the method is given by Fig. 4(c). Unlike two previous methods, this one does not suffer from the correlated nature of the added matrices and uses the natural decomposition of energy stratified models. Apart from an excellent quantitative description of the model's characteristics, the resulting differential equations give deep insights into how exactly the fractal phase emerges with lifting more and more hopping eigenvalues from zero.

Self-averaging
A key point which allows us to proceed with the above method is a statement about self-averaging of a quantity g|G|g which enters a denominator of (20b). 7 Indeed, if this statement fails, the r.h.s. of (20a) becomes highly nonlinear and extremely difficult to tackle analytically while in case of aforementioned self-averaging the matrix S depends only linearly on |g g| and doesn't depend on G at all.
To get an idea when and why this self-averaging condition holds, let us write the considered matrix element explicitly in the eigenbasis {|n } of H r−1 : Due to the imaginary part η of z, the main contribution to this sum goes from ∼ η/δ ∼ N α states in a vicinity of the energy ε (assuming that we are in the spectral bulk). Then if the corresponding components of |g are "ergodic", i.e., all | n|g | 2 are of the same order 1/N , then it is clear that the whole sum must only slightly fluctuate around some definite value of the order of unity. The above "local ergodicity" condition is likely to hold until the corresponding fractal dimension of |n (with |E n − ε| η) in the basis |i is less than one while for |g the fractal dimension is assumed to be equal to one. As a result, we come to a conclusion that this self-averaging holds at least as long as we work in the spectral bulk and in one of the non-ergodic phases. In other words, we claim that in the aforementioned cases (21) can be replaced by In the last equality we drop the real part of this quantity out due to its unimportance to the localization effects, and ρ = 1/(N δ) is a global density of states. For more details on the self-averaging derivation (complimented by the numerical evidences) please look into the Appendix A. Substituting (22) into (20), we obtain the self-averaged version of the Sherman-Morrison formula This formula is somehow similar both to the cavity method on the tree (15) (but non-linear in G r−1 (z)) and to the Dyson Brownian Motion approach (18). In next sections we consider separately two distinguished cases of (i) β < 1 with the measure zero of non-zero energy levels in the hopping term, and (ii) β = 1 with the finite fraction of such levels.
This separation is dictated by the following fact. The averaging over the hopping matrix V , Eq. (13), contains not only the averaging over its i.i.d. eigenvalues v r (4), but also over eigenvectors |g r . According to Eq. (23) this averaging can be done step by step over v r and g r with increasing index 1 ≤ r ≤ N β . However, at each rth step the averaging over |g r should imply the orthogonality of this vector to all the previous ones |g k<r . For β < 1 the measure of such vectors is zero among the overall basis |g m and therefore their effect on the averaging is negligible in the thermodynamic limit. This leads us to the simple results from Sec. 5.3 analogous, but not equivalent to the ones in the uncorrelated RP model [15].
Unlike this for β = 1 the orthogonality of the vectors |g r play a significant role and might make the results different in the case β = 1 with respect to the one given by a limit β → 1. 8 Therefore in Sec. 5.5 we develop the method combining the Sherman-Morrison formula (23) with the Dyson Brownian motion.

Cavity-like method for an intermediate regime β < 1.
Using the arguments given in the previous section for the case of β < 1 we average Eq. (23) over independent v r and g r with fixed r and obtain the recursive equation similar to the one of the cavity method on the tree structure (15) Here the integral in σ is given by the averaging over the probability distribution p(v) of the hopping eigenvalues v, while the averaged projector |g r g r | gives the unity matrix multiplied by 1/N 9 . In order to proceed with the averaging over all other g k and v k with k < r, one can replace the squared resolvent by the derivative G(z) 2 ≡ −∂ z G(z). After this trick the equation becomes linear in G(z) and can be straightforwardly averaged over g k and v k with k < r Finally, assuming that σG r−1 (z) is sufficiently small, we arrive at the desired differential equation for the mean resolvent: From this differential equation it is apparent that G(z; r) = F (z − σr) is a function of the propagating wave argument z − σr. With the initial condition G(z; 0) = (z − H 0 ) −1 , we get the resolvent averaged over the hopping in the form 8 For example, consider a model with V = cI = N k=1 c|g k g k |, i.e. with V being a non-trivial hopping term at any β < 1, but collapsing to just a finite energy shift at β = 1. It is clear that in this case D(β = 1) is zero while the β → 1 limit of D(β) is one. This example illustrates a highly correlated case when the aforementioned orthogonality really matters. 9 One may argue that for r > 1 the projector |gr gr| cannot be averaged to identity matrix due to the orthogonality conditions gp|gq = δpq. However, for β < 1, a number of the orthogonality conditions is of measure zero, and we will neglect their effect here and after. So, as the title of the current subsection implies, this result holds only for β < 1, but we will see below in Sec. 5.5 how it can be generalized to a full-rank case.
From the result one can see that the assumption of the smallness of σG r−1 (z) is valid until σ 1. The recursive equation (25) has much in common with the cavity method (15) as well as Eq. (26) is similar to the one of the Dyson Brownian motion (18). Even the result (27) for β → 1 is formally the same as the one of the cavity method (16). However unlike these methods both equations (25), (26) are derived for the correlated matrix elements V ij for which both cavity and DBM methods fail. In any case we nickname this method "cavity-like" due to the similarity between recursive equations.

A fractal dimension D
As it is apparent from the results of the previous section, it is only the parameter σ which matters for the determination of the self-energy Σ for the models (1) with β < 1 and independently distributed eigenvalues v k , and The real part of σ results in the energy shift. For the symmetric distribution p(v) it vanishes, while in general for the mid-spectrum energies ε, it can be neglected provided Re Σ(β) = Re σN β 1 is less than the bandwidth. Assuming this is the case, let's focus on the imaginary part which then determines the localization properties of the wave functions. The imaginary part, in turn, can be estimated up to the unimportant prefactor as where the first term in r.h.s. corresponds to the small v 2 1/ρ 2 , while the second gives the cutoff for the large values v 2 1/ρ 2 ; δ = 1/(ρN ) is the mean level spacing in the spectral bulk. Now, recalling the argument about the average site population following (12) and taking the imaginary part of (27), we get Together with (29) and the standard definition of the fractal dimension D via the support set of the wave function the Eq. (30) finally gives Substituting particular values (4), v 2 = N 2−γ−β , corresponding to our model (1) at γ < 2, we get restoring the dependence (11) already shown in Fig. 2. The overall phase diagram in the plane (β, γ) is given by Fig. 5. It shows how the Richardson's model (noted as RM) looses its localization properties at γ < 2 as soon as one adds an extensive number N β>0 of energy stratified hopping terms. Indeed, for any fixed 0 < β < 1 the localized phase at γ < 2 is replaced by the fractal one with the maximal fractal dimension given by D = β, see Fig. 5. This maximal fractal dimension can be physically understood as follows. The N − N β low-energy eigenstates of the translation-invariant model (without the diagonal disorder H 0 ) are degenerate and therefore should transform to the delta-function localized ones in the coordinate basis as soon as one adds the diagonal disorder H 0 to the problem. However this delta-function localization is limited by the presence of N β energy-stratified states. It is straightforward to show that these N β states are localized in the momentum {|g k } basis due to their extensive energies, while all the rest states should be orthogonal to them. Considering also that the hybridized low-energy states should be orthonormal to each other, we get that such a state typically occupies at least N β sites. This finally puts the minimal support set N β of the mid-spectrum states in our models and provides the necessary understanding of the result for the fractal dimension D = β instead of delta-peak localization with D = 0, see Fig. 1(e). The above result works for γ < 2 − β as soon as the high-energy states are energetically stratified from the bulk of the spectrum, Fig. 1(a-b).
On the other hand, in the interval of 2 − β < γ < 2 the fractal dimension is given by D = 2−γ as in the TIRP model, see Fig. 5. This corresponds to the bandwidth v 2 ∼ N 2−γ−β of the high-energy states to be of the order of the diagonal disorder or smaller. As a result, in such a case these N β levels are not anymore energy-stratified and the whole picture goes back to the TIRP case which practically shows the same result as the Fermi's golden rule consideration.
The above result in the limit β → 1 is also consistent with the previously known results for β = 1 despite the fact that this derivation is formally applicable only for β < 1. In the next subsection we make sure that this is not a coincidence.

Dyson-Brownian-motion-like stochastic process for TIRP: β = 1
Now, consider the full-ranked hopping matrix V of size N × N . As it was mentioned in the end of Sec. 5.2, the averaging over |g r in the form presented in Eq. (24) doesn't work for r ∼ N because during the averaging we have to exclude the subspace spanned by |g k with k < r, and such effects cannot be neglected anymore. It provides some correlations between the averaging on the current step and all the previous ones and makes the analysis much more involved.
In order to proceed, in the full-ranked version of our model (1) we use the analogy between the recursive representation (19) and the one of the Dyson Brownian motion (17). Indeed, like in the Dyson Brownian motion, one can consider a stochastic random matrix process over auxiliary time, 0 ≤ t ≤ 1, defined as a ∆t → 0 limit of the recursion (20a) where the random eigenvalue increment ∆v ∼ √ N ∆t and the random wavevector g = 1, N are chosen independently on each step. It is important to emphasize that here, unlike the cavity-like method case, Sec. 5.3, the parameter ∆v is not given by eigenvalue of the matrix V , but it should be found in such a way to sample the desired matrix ensemble at t = 1. This, in particular, means that each individual eigenvector |g k appears extensive amount of times ∼ (N ∆t) −1 1 (but not once) during the random process. In order to emphasize this we everywhere use N ∆t 1 as a small parameter. Hence, since this approach doesn't force the eigenvector |g added on some k'th step to be orthogonal to all previously added ones, the aforementioned orthogonality problem is not a problem anymore, and we can use the analogue of the averaged Sherman-Morrison formula (23) with v k replaced by ∆v k for each individual step of the stochastic matrix process.
In the following, we focus on the β = 1 TIRP model which is characterized by the Gaussian distribution of the eigenvalues (4). This immediately simplifies our stochastic matrix process to a process with the Gaussian eigenvalue increments ∆v = v √ N ∆t, v ∼ N (0, N −γ ). Now, using the self-averaged version (23) of the Sherman-Morrison formula and expanding it analogously to the Dyson Brownian motion in powers of N ∆t up to the first power, one obtains With all the above simplifications given by the random process, one averages the latter equation in the same way as (24) and derives the analogous equation From this equation one obtains Recalling that the desired full-ranked model arises from our stochastic process at t = 1 we finally get Σ(β = 1) = iπρN 1−γ which gives similarly to the previous section D = 2 − γ as expected [57]. One may notice the similar paradigm of this section as in the Dyson Brownian motion, therefore we use this name of "DBM-like" for this method. However, like in Sec. 5.3, the equivalence between Eq. (36) and (18) and formally the same results should not confuse a reader as they are derived for the TIRP model with translation-invariant correlations of the matrix elements V ij = V i−j for which the DBM method does not work.

Conclusions
To sum up, in this paper we consider the set of disordered random matrix models with translation-invariant correlations in the hopping term. This set of models interpolates between the completely correlated case of the Richardson's model and the translation-invariant version of the Rosenzweig-Porter (TIRP) model. The former model is Bethe ansatz integrable and known to have all (except one) eigenstates to be localized even beyond the convergence of the locator expansion, while TIRP model shows the same fractal wave-function statistics as the usual uncorrelated RP one. The interpolation between these models is given by a number N β , 0 ≤ β ≤ 1, of finite-energy levels in the spectrum of the translation-invariant hopping in the combination with the energy-degenerate rest eigenstates of the disorder-free translationinvariant model. Such kind of the energy stratification is known to play an important role in cooperative shielding and correlation-induced localization.
However quite surprisingly in contradiction to the cooperative shielding arguments, we have found that as soon as the number of energy-stratified levels starts to grow extensively with the system size N , the mid-spectrum states become delocalized beyond the locator expansion, but fractal.
In order to describe the above set of models and calculate the fractal dimension of the bulk spectral states, we develop the generalizations of the cavity and the Dyson Brownian motion methods, based on the natural spectral decomposition of the correlated hopping matrix complemented by the Sherman-Morrison formula for the local resolvent.
Our methods uncover the whole phase diagram of the above-mentioned set of models and shows that for all 0 < β < 1 the localization beyond the locator expansion, which is present in the Richardson's model, is replaced by the emergent fractal phase in the same parameter range.
The numerical simulations check the validity of the above methods and confirm the analytical results for the fractal dimensions. The first of these methods has a wider application than just a Gaussian distribution of the eigenvalues that will be considered in more detail in further publications.

Outlook and discussions
As a next step in the development of the approach suggested in this paper, one can consider the cases of either or both eigenvalues v k and eigenvectors g k of the hopping matrix to be correlated between each other. This direction has immediate applications in the description of short-range graph models and their mapping to the random-matrix ensembles (see, e.g. [49]).
Another promising direction is to consider the fat-tailed distribution of eigenvalues v k in the full-rank case, similarly to uncorrelated models [22][23][24]49] and uncover the possible nature of the emergent multifractality for such models. This is the subject of further investigations, but here we would like to notice that the DBMlike method developed in Sec. 5.5 cannot be immediately applied for the fat-tailed distributed v k . Indeed, as soon as the PDF of ∆v is not narrow and the fluctuations are comparable with the bandwidth, ∆v 2 − ∆v 2 ρ −2 the entire self-averaging argument breaks down. An origin of this fact lies in a behaviour of g k |G|g k on later steps of the random process when the state |g k was already used earlier. This quantity can be computed exactly using (20a) and looks like where v k (t) is the corresponding hopping matrix eigenvalue on step t. Since v k (t) is not zero for the later stages of the process, it can break the convenient self-averaging property of g k |G r+1 |g k and, if the fluctuations of ∆v k (and, hence, v k (t + ∆t)) are large, we cannot neglect this effect. As a result, S(v, g, G) no longer can be naively averaged to a multiple of unity independent of G which leads to significant complications and lies beyond the scope of the present paper.

A Self-averaging
Consider the matrix element g r |G r−1 (z)|g r in the denominator of the fraction from (20b). We claim that this quantity self-averages almost always for a range of parameters we are interested in. This claim is based on the assumption that the added on the r'th step hopping eigenstate |g r is generic with respect to the eigenstates |n r−1 of the Hamiltonian H r−1 from the previous step, i.e., that its fractal dimension in the basis of the Hamiltonian's eigenstates is one. And, if |g r is ergodic in the basis of H 0 than this is indeed the case at least until |n r−1 are not ergodic in that basis, i.e. until D(r) < 1. The only obvious thing which can break this logic is the orthogonality condition g r |g p = δ rp between |g r and |g p when |g p was added on the previous steps of the recursion (20a) thereby making |n and |g r correlated. As we know from the analysis of the Richardson model, for a very large eigenvalue v p , the corresponding detached eigenstate will be almost collinear to |g p and, hence, indeed orthogonal to all other |g r with r = p. However, as long as this happens only for detached levels, there is no reason to think about such effects if we talk about bulk states.
So now, when we agreed to work in the spectral bulk and in a non-ergodic phase, we can approximate n|g r as taken from the normal distribution and write for g r |G r−1 (z)|g r and some fixed G r−1 (z) that where x n are independently distributed standard normal random variables which, together with the factor 1/N , approximate the uniform over the unitary group distribution of |g in the large-N limit [66]. The distribution of τ (z) can be estimated by computing separately the characteristic functions of its real and imaginary parts. For example, Then, taking δ(ε) η ρ(ε)/ρ (ε), where δ(ε) = 1/N ρ(ε) is the energy-dependent mean level spacing and ρ(ε) is the density of states, we can approximate the sum in the exponent by the integral This integral can be approximately calculated for a fixed N in two limiting cases: in the case ξ N η when the logarithm can be expanded in the Taylor series and in the case ξ N E 2 BW /η where E BW is the bandwidth of the Hamiltonian. The desired asymptotics then take the form where is a constant independent of ξ. Now, recalling from Sec. 5.1 that the proper choice of η is η(ε) = δ(ε)N α with α > 0, we see that the region where the characteristic function Q Im (ξ) oscillates with a constant amplitude grows with the system size N , and beyond this region it decays polynomially with the power of N/2. All this points to the fact that the distribution of Im τ (z) becomes narrower with the increasing system size and concentrates around the mean point Im τ (z) = πρ(ε). The distribution function of Re τ (z) behaves in a similar way concentrating around π∂ ε C(ε, η), and this finally justifies our self-averaging approximation giving until we work in the spectral bulk and beyond the ergodic phase. The second and the fourth rows of Figure 6 demonstrate the self-averaging of the corresponding distributions. It is also worth mentioning that, in the limit η δ(ε), one cannot transform the sum, Eq. (40), into the integral, Eq. (41), since particular values of E n now matter. Instead of a single frequency 2πρ(ε), Q Im (ξ) has now a set of equally important realization-dependent frequencies and, hence, no self-averaging, see the first and the third rows of Figure 6.
Finally, since Re τ (z) = π∂ ε C(ε, η) for the bulk energies is of the order of the bandwidth or even parametrically smaller than that, we drop this real part and write just g r |G r−1 (ε − iη(ε))|g r ∼ iπρ(ε) which will simplify the calculations without qualitatively affecting the results or losing the generality.

B A self-consistent approach
In the case of not smooth density of states, one cannot drop the energy dependence from the self-averaging matrix element g r |G r−1 (z)|g r . In such a case one should consider the self-consistent formulation of the problem. In case of β < 1 it can be written as where σ now depends both on the energy z and the rank r. In case of Dyson Brownian motion-like approach for β = 1 it is ideologically similar and differs only by a definition of σ and "time" r. The self-consistency condition can be extracted from this system and treated independently: after taking a trace of the first equation and defining τ (z, r) = Tr G(z, r) /N we get a closed non-linear system    ∂ r τ (z, r) + σ(τ )∂ z τ (z, r) = 0 and the equation for G(z, r) then gives G(z, r) = G 0 (z − σ(τ )r) with the self-consistently determined σ(τ (z, r)). In turn, the system (46) can be easily solved using the method of characteristics which for the initial conditions τ (z, 0) = τ 0 (z) gives the solution in the form τ (z, r) = τ 0 (z − σ(τ )r) = τ 0 z − rσ τ 0 z − rσ τ 0 (. . .) .

B.1 Order-independence property
In general, there are R! different ways to obtain the same G(z) using the exact recursion (20a); the ways differ by the order we choose to lift the hopping eigenvalues v k from zero. It is reasonable now to check if our self-consistent equation (45) for the mean resolvent shares the same property. To do this, let's define a set of probability density functions p(v; k) such that Substituting these k-dependent functions to σ instead of p(v), we emulate different orders of v k in (20a) and arrive to the spectral factor σ explicitly depending on the rank r. Now, the equations for the characteristics of (46) take the form dτ dr = 0, dz dr = 1 N vp(v; r)dv 1 − vτ , (49) which leads to z(r) = z 0 + 1 N vdv 1 − vτ 0 (z 0 ) r 0 p(v; k)dk = z 0 + Σ(τ 0 (z 0 ); r) (50) and τ (z, r) = τ 0 (z − Σ(τ ; r)).
And, since the integral over k in (50) behaves very much like the sum from (48), the selfenergy Σ(τ ; R) (and, hence, the quantity τ (z, R)) is indeed order-independent. Finally, since the mean resolvent G(z, r) from our equations also has the form G(z, r) = G 0 (z − Σ(τ ; r)), its value for r = R is order-independent too as it should be.