Non-ergodic delocalized phase with Poisson level statistics

Motivated by the many-body localization (MBL) phase in generic interacting disordered quantum systems, we develop a model simulating the same eigenstate structure like in MBL, but in the random-matrix setting. Demonstrating the absence of energy level repulsion (Poisson statistics), this model carries non-ergodic eigenstates, delocalized over the extensive number of configurations in the Hilbert space. On the above example, we formulate general conditions to a single-particle and random-matrix models in order to carry such states, based on the transparent generalization of the Anderson localization of single-particle states and multiple resonances.


Introduction
The phenomenon of quantum thermalization in generic isolated quantum systems is highly non-trivial. Indeed, the unitarity of quantum evolution keeps all the information about the initial state (unlike the dynamics in classical chaotic systems). Therefore the relaxation of isolated quantum systems to the equilibrium thermal state, independently of the initial configuration, gives rise to an obvious paradox. This paradox is partially resolved by the ergodicity hypothesis [19,27,55,57,66,67], claiming that as time goes, the quantum information diffusively scrambles among extensive number of highly non-local degrees of freedom and, thus, most of physical (local) observables acting on a part of the system thermalize, because the rest of the system plays a role of a bath for them.
However, as now accepted in the literature, quantum thermalization may generically fail as one applies a strong enough disorder to the system. This happens due to the emergence of local integrals of motion and leaves a way for a new phase of matter, called many-body localization (MBL) [2,6,30,53,54]. The presence of an extensive number of the local integrals of motion [33,65] in the MBL phase prevents the quantum information scrambling and also localizes the (particle) transport across the system. In addition, the MBL phase is characterized by the Poisson level statistics (absence of level repulsion), while many-body wave functions, localized in the real space, occupy a zero fraction, but extensive number of Hilbert-space configurations [23,[43][44][45]70] 1 . Such states bear the name of non-ergodic extended or multifractal wave functions in the Hilbert space.
The complexity of the above mentioned many-body problems prevents one from analytical description of them. Thus, it is of particular concern and high demand to model essential attributes of these phenomena in a universal manner. In this respect, the randommatrix setting provides a perfect and more tractable playground, which can be used as quantum simulations of many-body systems. In the quantum thermalizing phase, the standard Gaussian random matrices [46] with the Wigner-Dyson level repulsion do their job. However going to the non-ergodic (especially MBL) phase, we face an immediate problem as the overwhelming majority of the random-matrix ensembles shows multifractal properties in the Hilbert space only at the very point of the Anderson transition [28]. There are only few ensembles like the Rosenzweig-Porter model [36,61] or some related ones [12,24,34,35,37,39,40,[49][50][51] which provide an entire phase of non-ergodic delocalized eigenstates. However, all such models show the standard Wigner-Dyson level repulsion in the whole non-ergodic phase, like in Gaussian random matrix ensembles 2 .
The origin of the non-ergodic Hilbert-space delocalization of MBL wave functions is at least two-fold. First, the structure of the Hilbert space matters. The number of available 1 There are many other physical quantities characterizing the ergodicity breaking like fluctuations of local observables [27,66,67], entanglement entropy [1,5,31], out-of-time-ordered correlators [29,32] and operator spreading [1], showing deviations from their ergodic values (sometimes differently from other observables [17,21]). In this paper, instead, we focus on those of these quantities which have a clear meaning in the random-matrix setting, namely, eigenstate localization properties in the Hilbert space and eigenlevel statistics. 2 There are also models based on the quasiperiodic potential instead of uncorrelated on-site disorder which do show the deviations from the Wigner-Dyson level repulsion of the non-ergodic extended states [56,[62][63][64], but unfortunately they also do not show Poisson level statistics either.
configurations N in a generic quantum system scales exponentially with the system size L, i.e. N ∼ e cL , and even if each single particle is localized in the coordinate space, it can still spread over few lattice sites there corresponding to a finite localization length ξ. It is this finiteness of ξ > 0, which leads to the delocalization in the Hilbert space. Indeed, one can estimate the number of configurations occupied by the M -particle wave-function with a finite filling fraction ν = M/L as ξ M ∼ N D , with D ≃ (ν/c) ln ξ > 0 providing the above non-ergodic delocalization (see, e.g., [23,45]). Second, the correlated nature of the diagonal matrix elements in the Hilbert-space configurations (given by the random uncorrelated on-site disorder and few-particle diagonal interactions in the coordinate space) leads to a rather complicated structure of multiple resonances [13,14] (unlike the case of the Anderson localization), which, in turn, delocalizes the many-body wave function in the Hilbert space.
To sum up, due to the presence of the Hilbert space in addition to the coordinate space, MBL provides an example of wave functions with Poisson level statistics, which are (non-ergodically) extended in the Hilbert space. Such a combination is typically absent in the picture of the Anderson localization or in any random-matrix ensemble.
In this paper we provide such an example in random-matrix setting and unveil its main ingredients, summarized below 3 .
1. In order to realize the structure of multiple resonances similar to MBL, one has to make the ratio of the amplitude of the hopping term j to the disorder amplitude W scaling up j/W ∼ N γ/2 , 0 < γ < 1, with the matrix size N . 4 2. The Poisson level statistics can be realized if resonant hopping terms between sites i and j, satisfying the inequality |j ij | > |ε i − ε j |, form either an effectively block diagonal matrix with the extensive number of blocks N B ≫ 1, or are only few per row and only at the large, but sub-extensive distance 1 ≪ R ≪ N from the site i.
An immediate straightforward example satisfying the above conditions can be realized in a one-dimensional Anderson model by the scaling up ratio j/W ∼ N γ/2 . Indeed, as soon as the wave functions exponentially decay on the scale of the localization length ξ ∼ (j/W ) 2 ∼ N γ , which, in turn, scales sub-extensively with the system size, 0 < γ < 1, the eigenstates are non-ergodic and delocalized, while as soon as all such states are effectively split into a sub-extensive number N B ∼ N/ξ = N 1−γ of blocks of independent eigenvalues the corresponding level statistics should be Poisson. This simple example provides some intuition of how the MBL eigenstate structure can be mimicked, however, in somehow trivial manner. In the following, we will focus on a more realistic models, relevant to many-body systems and experimental realizations.
In terms of the many-body systems, the above ingredients 1, 2 find applications in physically relevant random-matrix models [34,35,37]. Indeed, as soon as one considers the Hilbert-space structure of a generic many-body system [23,48,68], one immediately sees that the main contribution to the many-body delocalization is given by the Hilbert-space configurations separated by the extensive (Hamming) distance ∼ L ∼ ln N . The statistics of the corresponding matrix elements, coupling such configurations, appears to have two main properties related to the above mentioned ingredients: First, the typical matrix elements should scale as a power of the Hilbert-space dimension N (like in the item 1), (1)- (4). The brightness of blue-colored diagonal elements (left) and vertices (right) corresponds to the on-site disorder ε n , while the color from red to yellow and light blue of the off-diagonal elements (left) and edges (right) shows the power-law decay ∼ N γ/2 ||m−n|| a of hopping from/to the topmost vertex. In both panels we plot a realization of the matrix for a = 0.5, γ = 1, and N = 17. and, second, the distribution of them should be fat-tailed (at least log-normal [34,37,48]), leading to rare resonances (like in the item 2).
From the experimental point of view, the random-matrix models provide an excellent description of the Anderson transition in d-dimensional dipolar systems. Indeed, since 1989 it is known that in the random-matrix setting one can realize an Anderson transition at any d-dimensional lattice by considering a generalized dipolar hopping term (like in [41,42,47]) with the exponent a of the power-law decay j m,m+R ∼ 1/R a . For large a > d all the states are known to be (power-law) Anderson localized, while for small enough a < d they become delocalized and ergodic, and this transition is barely affected by any finite amplitude W of on-site disorder.
However, the latter delocalization property at small powers a < d of power-law longrange systems might fail in the dipolar systems, where all the dipoles are not randomly oriented [41,42], but aligned, e.g., by an electric field [15,20,25,26,38,51]. In the literature this is called the Burin-Maksimov model by the names of the authors of [15] first suggested it. In this case, the localization scenario is different: for all a < 3d/2 at moderate disorder W (and for all W at a < d) there are high-energy states which are delocalized and ergodic, but for d ≤ 2 their fraction in thermodynamic limit, N → ∞, is zero among the entire spectrum [20] 5 . Unlike that measure zero of ergodic states, the spectral-bulk eigenfunctions are all power-law localized for any a > 0 with the effective power-law decay exponent being symmetric with respect to a = d [25,51] Already this very Burin-Maksimov (BM) model, Fig. 1, places the question of the presence of non-ergodic extended states with Poisson statistics. Indeed, the localized bulk states, showing Poisson statistics, are separated from the spectral edge ergodic ones, which have spectral degeneracies and extensive energies scaling with the matrix size. This hints that the states in the buffer between the bulk (E ∼ O(1)) and the edge (E ∼ N d−a ) spectral states may combine the properties of both: they can be non-ergodic, delocalized, and demonstrate Poisson statistics at the same time. Unfortunately, they form again just a zero fraction of all the states in the spectrum. In order to overcome the last drawback of BM model, in this paper we consider its modification, with the bulk majority of eigenstates being non-ergodic, delocalized, and showing Poisson statistics by utilizing the first of the above ingredients, namely, by making the hopping amplitude scaling as j/W ∼ N γ/2 with respect to the on-site disorder.
Namely, we consider the above Burin-Maksimov model determined by the N × N matrix, Fig. 1

(left)
integers 0 ≤ m, n < N , and N independently identically distributed random numbers ε n with zero mean and the following variance To satisfy the condition of multiple resonances we use the ingredient 1 and scale the hopping amplitude t with the matrix size N For simplicity of analytical consideration we consider periodic boundary conditions and determine the distance ||m − n|| as the minimal length on the one-dimensional (1d) circle, see

Known methods
In this section, we remind several methods to determine localization known in the literature and apply it to the above model (focusing mostly on γ = 0). We consider the methods written in the papers [41,42,51]. Those who is interested in the results and current analytical consideration, can go directly to Secs. 3 and 4.

Anderson's resonance counting
Here we count the number of sites m which are in resonance, |ε n − ε m | < j mn , with a certain nth one. As soon as this number of resonances is finite, the state n is Anderson localized. Similarly to [41], we count the resonances at each distance ||m − n|| ∼ R where in 1d the level spacing ε n − ε m between closest levels is given by δ R ≃ W/R, while j m,m+R = N γ/2 /R a : if j m,m+R > δ R , the number of resonances will scale as ∼ R.
As one can see from above expressions, the number of resonances is finite only if most of resonances appear at small distances R ≲ O(1), i.e., for a > 1 and γ ≤ 0, or don't appear at all, j m,m+R < δ R for 1 ≤ R < N , i.e., for γ ≤ −2.
For the other parameters, the number of resonances formally diverges with increasing N . However, unlike the random (uncorrelated) models [41,42,47], the resonance counting fails to describe the wave-function structure in this case. The origin of this failure follows from the fact that this divergence in the models with the correlated kinetic terms leads to the delocalization of a measure zero of eigenstates (or even one state in Richardson model) that take the most weight of this kinetic term due to their high energies (for the details please see [51]). Therefore in the next subsections we consider another method invented in [51] and called the matrix-inversion trick.

Momentum-space localization and matrix-inversion trick
The matrix (1) corresponds to the Hamiltonian H =ε +ĵ (5) written in the coordinate basis {|n⟩} via the diagonal disorder ε = n ε n |n⟩ ⟨n| (6) and the translation-invariant hopping term The latter can be diagonalized in a momentum basis as the translation invariance of the hopping term j mn = j ||m−n|| and periodic boundary conditions, Fig. 1(right), allow one to find the spectrum ofĵ via the Fourier transform with q = 2πk/N and integer 0 ≤ k < N . Here and further all the summations over n and m are taken over the entire interval 0 ≤ m, n < N . General results for the spectrum j q given, e.g., in [51] take the form Here A a = Γ 1−a sin πa 2 , B a = 2(1 − 2 3−a )ζ a−2 ≃ a/2, C π = 2(2 1−a − 1)ζ a , ζ a and Γ a are the zeta and Gamma functions. Here and further we restrict our consideration to the range 0 < a < 2 for simplicity where the maximal value of j q is given by and the minimum is negative and equal to C π 7 . In the momentum basis where j q plays a role of the diagonal potential, the disorderε transforms to a translation-invariant hopping term Further in this part, we consider the localization of the states in the momentum space due to the large diagonal elements j q of the high-energy spectral edge states followed by the method to understand why the presence of such states may localize the spectral bulk states in the coordinate basis.

Momentum-space localization
Similarly to the previous section 2.1, one can count the number of resonances in the momentum space, but taking into account the fact that the level spacing δj q = |j q −j q+π/N | is highly q-dependent It is given by in the bulk q ≃ 1, while the ones at the edges of the q-range take the form Thus, the minimal level spacing is always (in our range of interest, 0 < a < 2) given by q = π, while the maximal is δj 0 . Taking all this into account, one can find at which qs the states are localized in the momentum space (therefore mostly represented by the plane waves in the coordinate basis): As we are not interested in the localization in the momentum space of most (measure one of all) bulk states, further we focus on γ < 1 and 0 < a < 2. In this case only zero fraction q * , of eigenstates are localized in the momentum basis, i.e. nearly plane waves in the coordinate basis Moreover they have an extensive energy for a < 1 These factors together are necessary for a so-called correlation-induced localization [51], when a measure zero of high-energy states causes the localization of the bulk spectrum states beyond the convergence of the resonance counting.

Correlation-induced localization and matrix inversion trick
The presence of (measure zero of) high-energy states, delocalized in the coordinate basis, leads to the reduction of the effective hopping term for the bulk spectral states [16,26,49,51]. Indeed, on one hand, using the localization of the high-energy states |q| ≪ q * in the momentum space, one can find the exact eigenstates with the energy E q ≃ j q of the entire Hamiltonian (5) perturbatively Thus, the hopping term (7) can be rewritten using the above eigenstates aŝ where due to the divergence of eigenenergies (13) at a < 1 the first term gives the main contribution, while the term j eff , corresponding to the summation over the states, ergodic in the momentum space, |q| > q * , and to the perturbative terms from Eq. (19), is drastically reduced with respect to the initial one. In a sense, j eff corresponds to the effective hopping term in the Hilbert space sector, orthogonal with respect to the high-energy states, which are non-ergodic in the momentum space. On the other hand, as |ψ q ⟩ are the exact eigenstates of the entire Hamiltonian (5) they should be orthogonal, ⟨ψ q |ψ E ⟩ = 0, to all other eigenstates |ψ E ⟩. Therefore the effective hopping for |ψ E ⟩ is just j eff This already hints that these bulk eigenstates can be localized in the coordinate space even beyond localization perturbation theory of Sec. 2.1 due to the presence of correlations in the hopping term j mn . In order to estimate j eff , in [51] the authors invented a so-called matrix-inversion trick 8 Here we take the freedom to remind it to the reader. The main idea behind the matrixinversion trick is to get rid of the divergence in the spectrum j q (13) by inverting the hopping matrix. Before that, one should shiftĵ by an identity multiplied by a certain scalar E 0 ∼ O(1) in order to not only make former large-energy terms small, but also to avoid adding some large-energy terms to the inverted model due to the possible resonances: In such a way one can immediately show localization of the bulk eigenstates as the above effective hopping term for a < 1 decays as together with the eigenstate tails confirming the duality suggested in [25]. In the further sections we generalize this method to t ∼ N γ/2 .

Spatial renormalization group
In this section, we consider a renormalization group (RG) derived for the BM model in [15] and developed further in [38] beyond the strong-disorder limit. The main idea of the RG is as follows. First, one cuts off the tunneling beyond a certain scale R 0 and calculate the wave functions (R 0 modes) for this scale. Then one chooses a new cutoff R > R 0 and constructs new modes (R modes) a superposition of R 0 modes, taking into account the resonances between R 0 modes. The localization length increases from l 0 ≲ R 0 to l 1 ≲ R due to the presence of these resonances. Due to the smallness of the parameter j eff /W ≪ 1 only pairs of resonances are taken into account. The R modes |ψ (1) k ⟩ can be written via the initial site basis vectors |n⟩ as follows Thus, the hopping term rewritten in the new operators takes the form According to RG assumption the modes ψ (1) k (m) are localized within the interval |k − m| < l 1 at the length l 1 ≲ R for the deterministic potential j mn = t/||m−n|| a . Therefore one can neglect the difference between j mn and j kl (||i − j| − |k − l|| < |i − k| + |j − l| < 2l 1 ≲ R) due to the smoothness of the potential (see [38] for more rigorous consideration). As a result, Eq. (26) reads as In order to estimate the renormalized hopping term tℓ k ℓ * l /||k − l|| a one should consider the typical amplitude of ℓ k at a certain energy E as follows where the global density of states (DOS) is given by G m−n and G q stand for the Green's function in the real and momentum space satisfying the equation (E −Ĥ)Ĝ(E) =1, and the bar on top it stands for the disorder averaging over the diagonal elements ε n . Using the standard coherent potential approximation (see, e.g., [73]), for completeness considered in Appendix A, one can find both the averaged Green's function, diagonal in the momentum space, and its self-energy part Σ given by the following self-consistency equationsĜ where P (ε) is the distribution of the diagonal disorder ε n .
Note that the above results can be checked numerically by calculating the eigenstate spatial decay from the maximum k = n |ψ typ | 2 (R) = exp ⟨ln |ψ n (k = n + R) The calculation of the moment of the wave function ⟨|ψ| q ⟩(R) has a divergence over the energy difference ε n − ε k in the denominator and, thus, one can use the Breit-Wigner approximation in order to calculate it (see, e.g., [38]): and, thus, From this, one can easily calculate the spectrum of fractal dimension f (α) via the number of points at the distance R where the typical wave function scales as |ψ typ |(R) 2 ∼ N −α R . With this general description, in the further sections we focus on the Burin-Maksimov model [15] in 1d [25,51] generalized to the N -dependent prefactor of the hopping term.

Numerical results
Starting with the numerical simulations, we show the eigenstate spatial structure, including the short-range behavior and the long-range power-law decay (31). Then we explicitly determine the spectrum of fractal dimensions and show it for mid-spectrum states. In order to demonstrate the similarity to the MBL problem in the end of the section we also show the adjacent gap ratio statistics of eigenvalues. Figure 2 shows the long-distance power-law decay, Eq. (31), from the eigenstate maximum In order to emphasize the power-law nature of the decay both panels are shown in log-log scale. Like in the standard BM model [51] or in its Euclidean version [38] the exponent of this power-law decay is symmetric with respect to the critical point a = d = 1. However, unlike the above models one can see that the short-range behavior of the wave functions depends on the system size.
In order to clarify the system-size dependence of the short-range eigenstate decay, in Fig. 3 we show the zoom of this distance range for several power-law exponents a = 0.5 and 1.5 and hopping scaling γ = 0.25 and 0.5. From the log-linear plots one can immediately recognize an exponential decay of eigenstates, similar to the one in the shortrange Anderson models, |ψ E (R)| ∼ e −R/ξ . In addition, one can see that the effective localization length ξ of this exponential decay grows with the system size scaling as a power-law ξ ∼ N γ .  This scaling is extracted from the exponential fits of the data in the main panels and shown in the corresponding insets.
Next, we consider the spectrum of fractal dimensions, which is defined as the power f (α) of the scaling of the distribution P This function f (α) characterizes the fractal dimensions D q via the wave-function moments averaged over disorder and via the corresponding Legendre transform Among many properties of f (α) (like the normalization conditions of the probability distribution, and of the wave function) we would like to focus on the so-called multifractal symmetry [28] f which relates wave-function tails (large α > 1) with peaks (small α < 1). This symmetry is known to work for the extended states, fractal or multifractal with some (at least critical) level repulsion. However as we show in Fig. 4, the spectrum of fractal dimension in the bulk spectral states of the generalized BM model, (1) -(4), violates this symmetry by combining an approximate plateau f (α) ≃ γ at α ≳ γ corresponding to the exponential decay shown in Fig. 3 with the linear ramp f (α) = (α − γ)/(2a eff ) governed by the powerlaw decay (35). The violation of the symmetry (39) implicitly confirms the Poisson level statistics.
In order to demonstrate the Poisson eigenlevel statistics explicitly, in Fig. 5 we show the adjacent gap ratio statistics introduced in [3,53] as Here we will consider this r-statistics averaged over disorder realizations, r = ⟨r n ⟩. The absence of the level repulsion, i.e. Poisson statistics, corresponds to r = r P = 2 ln 2 − 1 ≃    The data averaged over 400 disorder realizations is extrapolated from the system sizes N = 2 12 -2 14 to infinity using the standard extrapolation techniques (see, e.g., [51]).
0.3863, while the Wigner-Dyson level repulsion [46] leads to r = r GOE ≃ 0.5307 for the orthogonal symmetry to which we stick here 9 . Figure 5 shows r-statistics, averaged over disorder, versus the energy. As the spectrum in BM model is highly inhomogeneous (13), in addition we show the corresponding fraction of states at each energy, given by the integration of DOS and called a cumulative distribution function. One sees the Poisson level statistics, r = r P , in the extensive interval of energies corresponding to most of the states, except a small fraction (which is of the order of 20% for the considered system size) of the topmost energy states, corresponding to the momentum-space localization, (17) 10 .
In order to uncover global spectral correlations, we have also considered a so-called power spectrum [8,9,18,58,[58][59][60], see Fig. 6. Unlike r-statistics, this measure, similar to the level rigidity, is a global spectral characteristics, which can distinguish not only Wigner-Dyson and Poisson level statistics, but also can pinpoint the energy scale of the crossover between them [8,9,18]. In the simplest incarnation [8], which does not need any unfolding, the power spectrum P S k is the intensity λ 2 k of the eigenvalues of a singular value decomposition of the rectangular matrix E, containing N R realizations, E k,n , 1 ≤ k ≤ N r , of the spectrum, with N eigenvalues in each, 1 ≤ n ≤ N : Here unitary matrices U and V are N r × N r and N × N , respectively, and we consider N r ≤ N and order eigenvalues λ k in the non-increasing order λ k ≥ λ k+1 . Note that the singular values at k ≪ N r capture the global trends of the spectra, while the ones at k ∼ N r represent the local fluctuations (similar to the r-statistics).
In such a definition, the Wigner-Dyson level statistics corresponds to the 1/f -"noise" behavior of the power spectrum [71,72] (dotted black line in Fig. 6, f = 2πk/N r , while the pure Poisson level statistics gives P S ∼ 1/f 2 (dashed black line).
From the Fig. 6, one can see that for all considered 0 < γ < 1 the the high-k behavior of the power spectrum coinsides indeed to the Poisson one in agreement with the r-statistics, while the low-k P S-values show the tendency of even faster decay. The finite-size powerlaw exponent (extracted from a power-law fit, green line) is significantly larger than 1 and flows towards 2 with increasing system size. All this confirms the Poisson level statistics in the bulk of the spectrum, while the low-k contribution is related to the degeneracy of the high-energy momentum-space localized states.

Analytical consideration
In this section, we provide the analytical consideration of the generalized BM model, (1) -(4), based on the renormalization group analysis from Sec. 2.3 supplemented by the coherent potential approximation, Appendix A.
First, we consider the long-range asymptotic power-law decay of the eigenstates (31). For this purpose, we take into account the Lorenzian form (30) of the imaginary part of the Green's function in the momentum basis where for simplicity we consider the box distribution of the diagonal disorder ε n and focus on the limit in the spectral bulk, whereḠ mm (E) ≪ 1/W , i.e. the level broadening, Γ = Im Σ, being the imaginary part of the self-energy Σ, is given by Eq. (85), Γ = πg(E)W 2 /12. In this case, one can calculate the effective hopping matrix element as follows with the DOS given by (29): leading to In the latter equality of (44) we substitute the expression for the spectrum j q , Eq. (13), to Eq. (42) and used the solution q E of the equation j q E = E, while the latter equality of Eq. (45) should be considered as the definition of the scaling parameter b.
In addition, the above substitution determines the crossover length scale R * , separating the specific short-range behavior from the universal power-law decay. For a < 1 this crossover length scale takes the form and at large distances R ≫ R * , the effective power-law decay switches to In the opposite case of a > 1 at R ≫ R * only the prefactor may change with respect to the bare hopping, but not the power-law decay Next, in order to understand the short-range structure of eigenstates, one should compare the maximal J * and typical j typ q hopping, Eq. (13), as well as their level spacing with the broadening Γ. Here and further we will focus on the scaling with N in order to distinguish ergodic, fractal, and localized states. Those who is not interested in the overall analysis of the cases, but would like to focus on the bulk-spectrum delocalization with the Poisson statistics, can directly go to the item 3: 1. In the case of large broadening Γ ≫ J * > |j q |, corresponding to γ+2b < min (0, 2a − 2), one obtains in the bulk the Anderson localization with g(E) ∼ 1/W , Γ ∼ W , b = 0, and the non-rescaled result at γ < min (0, 2a − 2) ≤ 0 similarly to the usual BM-model at a > 1.

2.
In the case when the broadening is in between the maximal and typical j q , j typ q ≪ Γ ≪ J * , corresponding to 2a − 2 < γ < −2b ≤ 0 (i.e., a < 1) the Green's function has two contributions The first one is reminiscent to the box DOS, with being the fraction of q with j q < Γ, while the second term gives the high-energy tail of the DOS. Here 1/N ≫ q * ≫ 1 is given by thus, f → 1 in the thermodynamic limit. Consider separately the bulk and the edge spectral states.
In the bulk of the spectrum, |E| < Γ, the second term in (50) is not resonant j q > Γ > E and is given by Assuming that b < 1 one will obtain the usual box term to be dominant giving again Γ ∼ W , with b = 0. The effective hopping in this case At the edge of the spectrum |E| ≫ W the main contribution is given by the corresponding pole E = j q at q = q E and gives the result g(E) ≃ |dq/dj q | q=q E ≪ 1/W and Γ ∼ W 2 g(E). This pole results in the hopping term with the resonance at the extensive distance R ∼ 1/q E ≫ 1. The value at the resonance is very sensitive to the position of the energy E between the adjacent levels j 1/R , with R closest to 1/q E . Thus the variable varies in the range As soon as |E| ≫ W , ∆E ≫ Γ and the average resonance hopping takes the form In addition, as ∆E ≫ Γ such hopping term is resonant for the only site from the initial site, which is located at the distance R ∼ 1/q E . As a result, such a resonant consideration is equivalent to the effective 1d Anderson model with the exponential localization and the localization length given by the resonant hopping j res as |E| ≫ W and, thus, q E ≪ q * = N − |γ| 2(1−a) . Note that here the normalization condition of the wave function is taken into account by the prefactor 1/ξ.
As soon as localization length ξ scales up with the system size N , the corresponding eigenstates at the spectral edge are delocalized. These states become ergodic, when the localization length becomes of the order of the system size, ξ ∼ N 1 , i.e at 11 The states in the middle are delocalized, but non-ergodic. This works as well for the conventional Burin-Maksimov model, γ = 0, however, as mentioned in the introduction, such non-ergodic delocalized states form a measure zero of all states.
3. In the most interesting case when the broadening is much smaller than the typical hopping, but much larger than its level spacing, δj typ q ≪ Γ ≪ j typ q , corresponding to 0 < γ + 2b < 2, the fractal properties of the edge states come to the bulk of the spectrum. Indeed, on one hand, in Eq. (50) the fraction f → 0 goes to zero in the thermodynamic limit and, thus, the density of states in the spectral bulk is given by the pole contributions, g(E) ≃ |dq/dj q | q=q E ≪ 1/W and Γ ∼ W 2 g(E).
On the other hand, the broadening makes this spectrum continuous and the effective hopping reads as has not only the resonance hopping like (58) at R ≃ 1/q E , but also the small powerlaw tail at R ≫ 1/q E from j 1/R ∼ tR 1−a (47) for a < 1 and (48) for a > 1 Both these power-law terms are small compared to the diagonal disorder W and, thus, perturbative and cannot lead to delocalization.
Eventually in the bulk of the spectrum E ∼ t, q E ∼ N 0 ≪ 1 and, thus, Γ ∼ W 2 /t, corresponding to b = −γ/2. This limits the corresponding case to the range However, the presence of resonant hops j res ∼ t at R ≃ 1/q E ∼ O(1) should delocalize the wave function similarly to the edge spectral states in the previous item. Indeed, the resonances form the short-range Anderson model with the random hopping terms on top of the weak rescaled power-law tails. Thus, we should have the exponentially decaying wave function (59) with the N -dependent localization length (60) This localization length exceeds the system size for γ > 1 and leads to the localization at γ < 0 in agreement with the considered range of γ and the properties of the bulk spectral states in the previous items.
At the top part of the spectrum similarly to the previous section there is the ergodicnonergodic crossover in energy giving the same effective momentum (61). The bottom of the spectrum where j q ∼ t(π −q) 2 is also a bit special, but we skip its consideration here for brevity.
4. In the last case when Γ ≪ δj typ q most of the states are just localized in the momentum space.
The above analytical results for the case 3, 0 < γ < 1, can be summarized as follows: both for a < 1 and a > 1 the effective hopping term (62) after renormalization combines the resonances at j 1/R ≃ E leading to the exponential decay (59) of the bulk eigenstates with the localization length (36) growing with the system size N . At larger distances, R ≳ ξ ∼ N γ , the hopping (62) becomes power-law decaying and perturbative giving for a < 1 Eq. (63) and for a > 1 Eq. (64) similarly to the case of the conventional BM model (35). This agrees very well with the numerical simulations of the typical wave-function decay, Eq. (31), shown in Figs. 2 and 3. The above exponential (59) and the power-law (62) decay result in the piecewise linear f (α), Eq. (34). The exponential decay with the prefactor 1/ξ corresponds to the plateau at r ≪ ξ, where |ψ E (r)| ∼ 1/ξ and to a sharp cutoff at r > ξ, corresponding to α < γ. The power-law decay leads to the linear increase This result is in fair agreement with Fig. 4. The absence of the level repulsion, shown in Fig. 5 is confirmed by the analytical consideration of the resonances. Indeed, as effectively the model including only resonance couplings is equivalent to the 1d Anderson model with the hopping on the distance ∼ 1/q E ≫ 1, it effectively separates the spectrum into N q E ≪ N blocks of the size 1/q E and therefore leads to the Poisson level statistics.

Conclusion and outlook
In this paper, we develop the framework in the random-matrix setting, which mimics the properties of the eigenstates in the many-body localized phase, namely, the non-ergodic delocalization of the wave functions in the Hilbert space combined with the Poisson level statistics. For this, we formulate the ingredients needed for the realization of the above wave-function structure in random-matrix ensembles and apply them to a concrete model, relevant for the experimental observation in trapped ions, ultra cold and Rydberg atoms.
Namely, we consider a so-called Burin-Maksimov model with the long-range correlated hopping term which does not lead to the delocalization of the bulk spectral states, due to the renormalization of the power-law tails of the hopping, similar somehow to the structure of the local integrals of motion emerging in the many-body localized phase of disordered interacting models. We generalize this model to the case of the non-ergodic delocalized eigenstates in the bulk of the spectrum, by rescaling the hopping term with the matrix size N (mimicking the Hilbert space dimension). The developed analytical approach allows us to investigate the eigenstate statistics and spatial profile and to find the structure of resonances in the effective hopping. It is shown that the eigenstates both for shortand long-range models, first, decay exponentially with the localization length ξ, scaling up with N , and then go to the power-law decay. It is this localization length scale-up, ξ ∼ N γ , which makes the wave function non-ergodically delocalized for 0 < γ < 1 and forms effectively the extensive number of blocks of repulsing eigenvalues. The latter leads to the Poisson statistics. The numerical simulations explicitly confirm the above analytical investigations.
These results provide the way to mimic the wave-function structure in the Hilbert space of the many-body localized phase of matter with rather simple random matrix ensembles and open many possibilities to investigate further properties of the eigenstates using such probes as the subdiffusive wave-packet spreading in the Hilbert space (see, e.g., [7,10,11,22,52]), the logarithmic entanglement growth in the static and driven setting [5], the relation of the entanglement to non-ergodicity [17,21] or even out-of-time-ordered correlators [29]. Although the presence of the extensive number of ergodic eigenstates at the spectral edge move our model closer to the many-body mobility edge, present just below the MBL phase, the zero fraction of these states makes our model more suitable for the description of the very vicinity of the MBL transition and finite-temperature MBL for smaller disorder amplitudes.
The more deep understanding of global spectral properties would be also very helpful in order to understand how to distinguish the Poisson statistics of the localized and nonergodic delocalized wave functions. or after multiplying by ⟨k| whereḠ kk (E) = ⟨k|Ĝ(E) |k⟩. Note thatḠ kk (E) ≡Ḡ 0 (E) does not depend on k if the basis {|q⟩} momentum basis | ⟨n|q⟩ | 2 = 1/N , where N is the size of the system. Averaging the latter equation over a certain distribution P (ε) of ε = ε k and dividing by a common factor ⟨k|Ĝ(E) we get the self-consistency equation Decomposing 1 into P (ε)dε we get In the left equation we divided byḠ 0 (E), while in the right one we take into account (79). Eventually, substituting the right equation of (80) to the right one of (75), we obtain • For the Lorenzian distribution P (ε) = (W/π)[ε 2 + W 2 ] −1 we get from (79) leading to Σ = −iW for anyḠ 0 (E).
• On the other hand, the box distribution leads to The first result is compatible with the Fermi Golden rule as in this case Im Σ = π ε 2 n g(E), ε 2 n = W 2 /12.