The Many-Body localization transition in the Hilbert space
Marco Tarzia
cond-mat.stat-mech cond-mat.dis-nn







LPTMC, CNRS-UMR 7600, Sorbonne Université, 4 Pl. Jussieu, F-75005 Paris, France

Institut  Universitaire  de  France,  1  rue  Descartes,  75231  Paris  Cedex  05,  France


In this paper we propose a new perspective to analyze the many-body localization (MBL) transition when recast in terms of a single-particle tight-binding model in the space of many-body configurations.
We compute the distribution of tunneling rates between many-body states separated by an extensive number of spin flips 
at the leading order in perturbation theory starting from the insulator, and determine the scaling of their typical amplitude with the number of accessible states 
in the Hilbert space. 
By using an analogy with the Rosenzweig-Porter random matrix ensemble, we propose an ergodicity breaking criterion for the MBL transition based on the Fermi Golden Rule. According to this criterion, in the MBL phase 
many resonances are formed at large distance 
from an 
infinite temperature initial state, but they are not enough for the quantum dynamics to decorrelate from it in a finite time.
This implies that, differently from Anderson localized states, in the insulating phase many-body eigenstates are multifractal in the Hilbert space, as they occupy a large but subexponential part of the total volume, in agreement with recent numerical results, perturbative calculations, and intuitive arguments. Possible limitations and implications of our interpretation are discussed in the conclusions.




The Many-Body localization transition in the Hilbert space
    M. Tarzia
        ^1 Department of Statistical Sciences, University College London, London WC1E 6BT, UK

    ^2 School of Mathematical Sciences, Queen Mary University of London, London E1 4NS, UK

^3 The Alan Turing Institute, The British Library, London NW1 2DB, UK 


===================================================================================================================================================================================================================================================================










§ INTRODUCTION


Quantum systems of interacting particles subject to sufficiently strong disorder will fail to come to thermal equilibrium when they are not coupled to an external bath even though prepared with extensive amounts of energy above their ground states. 
This  phenomenon, commonly referred to as Many-Body Localization (MBL),
was originally predicted by Anderson <cit.>, but firmly established only during the last 15 years, after the famous breakthrough of  <cit.>, 
and corresponds to a novel dynamical out-of-equilibrium quantum phase transition due to the interplay of disorder, interactions, and quantum fluctuations <cit.>. 
Its existence has received support from perturbative <cit.>, numerical <cit.>, and experimental studies <cit.>, as well as rigorous mathematical approaches <cit.>.
These investigations have shown that the main feature of the MBL phase is a robust effective integrability <cit.>: an extensive set of quasilocal integrals of motion emerges, providing an intuitive explanation of the breakdown of thermalization, 
and producing several unusual and remarkable consequences, such as the absence of dc transport <cit.>, the violation of the eigenstate thermalization hypothesis <cit.> along with common concepts of equilibrium statistical mechanics <cit.>, and the area-law entanglement of eigenstates <cit.>.

In the latest years these remarkable phenomena have attracted a huge interest (see <cit.> for recent reviews), predominantly from the fact that MBL can protect quantum correlations from decoherence even at finite energy density and for arbitrarily long times. Yet, despite an impressively wide amount of work and several significant progress, many important problems remain open, especially concerning the critical properties of the transition <cit.>, the existence of MBL in higher dimensions d>1 <cit.>, and the anomalous diffusion and out-of-equilibrium relaxation observed in the “bad metal” regime preceding the MBL <cit.>.
In this context, simplified effective models might naturally play an important role to sharpen these questions and provide a playground to explore the nature of MBL and improve our understanding of it. 

In this respect, a paradigmatic route which gives a very intuitive picture of MBL is obtained by recasting the many-body quantum dynamics in terms of a single-particle tight-binding problem in the Hilbert space (HS) <cit.>. Within this mapping many-body configurations are seen as site orbitals on a given graph (with strongly correlated diagonal disorder) and the interactions play the role of an effective hopping connecting them. 
Although the local structure and topology of the graph depend on the specific form of the many-body Hamiltonian and on the choice of the basis, the HS is generically a very high dimensional disorder lattice. 
 It is therefore tempting to argue that single-particle Anderson localization (AL) on the Bethe lattice <cit.> can be thought as a pictorial representation of MBL, as put forward in the seminal work of <cit.>, and later further investigated in Refs. <cit.>. On the one hand, several intriguing observations support this analogy: The critical point of the MBL transition is expected to be in the localized phase, as for the Anderson model in d →∞ <cit.>; Recent phenomenological renormalization group (RG) approaches <cit.> predict KT-like flows for the MBL transition characterized by two relevant localization lengths (the typical and the average one) which are expected to exhibit the same critical behavior as in AL on the Bethe lattice <cit.>; Moreover,  finite-size effects close to the Anderson transition on the Bethe lattice reveal a non-monotonicity <cit.> which is also characteristic of 
 the MBL transition.
 
 
 On the other hand, however, 
there is also a major 
difference regarding the spatial extension of many-body wavefunctions 
in the insulating regime. 

 The statistics of eigenstates is generally characterized by their fractal dimensions D_q, defined through the asymptotic scaling behavior of the moments 
 Υ_q = ⟨∑_i | ψ_α (i) |^2 q⟩ with the size of the accessible volume V:  Υ_q ≃ V^D_q(1 - q) [q=2  recovers  the  usual  inverse  participation (IPR)].
For a perfectly delocalized ergodic states similar to plane waves D_q = 1.  Conversely, if a state is localized on a finite volume, one gets D_q = 0, as for single-particle AL. In an intermediate  situation,  wavefunctions  are  extended  but nonergodic, with 0<D_q<1.
In contrasts  with  the  well established  case  of  AL,  where  the  spatial  extension  of single-particle orbitals is known to display genuine multifractality only at criticality (in any dimension d≥2 <cit.>), recent numerical results <cit.>, as well as perturbative calculations <cit.>, strongly indicate that the many-body wavefunctions are multifractal in the whole MBL regime. 
This result can be easily rationalized considering for instance a quantum spin chain at strong disorder: Deep in the insulating regime most of the spins are strongly polarized 

due to the random fields (paraller or antiparallel to it depending whether one is looking at the ground state or excited states); Yet a small but finite fraction ρ 
of them 
remains “active” on the sites where the random field is smaller than the typical tunneling rate 
for spin flips <cit.>. 
For a chain of length n, many-body wavefunctions then typically occupy 2^n ρ configurations in the spin configuration basis, which yields D_2 ≈ρ. 
	Note that at strong disorder the fraction of active spins is proportional to Γ / h <cit.>, where Γ is the transition rate for spin flips and h is the width of the disorder distribution (see also Sec. <ref>), and thus is expected to vanish only in the limit of infinite disorder.

Another way of rephrasing the same concept is as follows: In absence of interactions, single-particles orbitals are Anderson localized over a localization length ξ_ loc. Once the interactions are turned on, 
 at strong enough disorder many-body eigenstates are expected to be weak modifications of the Slater determinant of single-particles orbitals. 
 In other words, they must be eigenstates of  the LIOMS <cit.>, which 
 are essentially linear combinations of  local operators 
 over a finite length of order ξ_ loc. Hence, 
 
 in the Fock space  of the occupation number of single-particle localized orbitals  they roughly occupy a volume of the order order ξ_ loc^n, yielding roughly D_2 ≈log_2 ξ_ loc <cit.>.

These simple arguments indicate that many-body wavefunctions cannot be Anderson localized on a finite volume (except at infinite disorder) and are thus generically nonergodic multifractal in configuration space (with fractal dimensions 0< D_q  < 1), although for a rather simple reason.
Of course, the localization properties of many-body eigenstates depends crucially on the choice of the basis. However, this is true also for single-particle Anderson localization (e.g., fully delocalized eigenvectors in real space are fully localized in momentum space).  
	The results of Ref. <cit.> clerly indicate  that the multifractal nature of many-body eigenstates  is the same for two 
	relevant choices of the basis that diagonalize the many-body Hamiltonian in specific limits in which the system is completely localized, and are thus used as a starting point for the l-bits construction or efficient numerical simulations of 
	MBL.

In this paper we put forward a novel perspective to analyze the 
MBL transition in the HS, providing a clear explanation of the difference between AL and MBL. We focus on 

the non-local propagator 
G_0,nx which connect a T=∞ randomly chosen initial state | 0 ⟩ (e.g., in the middle of the many-body spectrum) with the configurations at large extensive distance nx from it (separated by nx spin flips, with 0<x<1). These matrix elements encode the probability that a system being in the state | 0 ⟩  at t=0 is found in a state which differs from it by nx spin flips after infinite time.
We evaluate these amplitudes using the Forward Scattering Approximation (FSA) <cit.>, i.e., at the leading order in perturbation theory starting from the insulator, and analyze the asymptotic scaling behavior of their typical value  with the number N_nx of accessible configurations at distance nx from | 0 ⟩ and in the same energy shell, obtaining |  G_0,nx^ FSA |_ typ∝( N_nx)^-γ/2. The exponent γ  increases with the disorder strength. Based on an analogy with the so-called Rosenzweig-Porter (RP) random matrix ensemble <cit.>, and following the ideas of <cit.>, we  put forward  a criterion for ergodicity breaking based on the Fermi Golden Rule (FGR): 
If γ < 1 the escape rate of the initial state | 0 ⟩ 
is much larger than the spread of energy level due to disorder and the system is in the fully ergodic phase; Conversely, if γ >1 the initial state only hybridize with a sub-extensive fraction  
V^1 - γ of the total configurations, 

thereby producing nonergodic multifractal wavefunctions which only occupy V^2 - γ configurations close in energy <cit.>.
According to this interpretation, although in the MBL phase many resonances are formed at large distance in the configuration space, they are not enough to ensure ergodicity and to allow  the quantum dynamics to decorrelate from the initial condition in a finite time.

From the point of view of single particle hopping in the HS, MBL is thus reminiscent of the transition from ergodic to multifractal states of the  RP random matrix ensemble <cit.> (and its generalizations <cit.>), 
and not to the transition to AL, which instead occurs at γ = 2, corresponding to the requirement that the number of resonances found from | 0 ⟩ stays finite in the thermodynamic limit <cit.>.
In fact we find that γ tends to 2 from below in the limit of infinite disorder (which can be treated analytically within the FSA), implying that the many-body wavefunctions become truly Anderson localized on a finite volume of the HS only when the density of “active spins” vanishes and the bare localization length ξ_ loc≲ 1 (i.e., when the LIOMS τ_i^z ≃ S_i^z), in agreement with the intuitive arguments given above.

We apply this approach to three different 1d models commonly used in the context of MBL <cit.>, showing that 

the ergodicity breaking criterion γ_c=1 yields an estimation of the critical disorder h_c in strikingly good agreement with the one obtained from the most recent numerical studies of systems of approximately the same size of the ones considered here. This observation supports the robustness of our conclusions and the validity of the ergodicity breaking criterion based on the FGR.
We also backup the perturbative analysis by inspecting the signature of nonergodic multifractal eigenstates by probing the non standard scaling limit of the spectral statistics in exact diagonalizations of the many-body Hamiltonians of small sizes <cit.>. 


All in all, our interpretation fully supports the picture recently proposed in Ref. <cit.> where MBL is seen as a fragmentation of the HS (see also Refs. <cit.> where similar ideas were promoted to explore the analogy between MBL  and a percolation transition in the configuration space).

The paper is organized as follows. In the next section we define the model and its HS representation as a single-particle tight-binding problem. In Sec. <ref> we present the results obtained within the FSA for the scaling of the matrix elements between distant states in the HS. In Sec. <ref> we examine the analogy with the RP model and discuss the ergodicity breaking criterion based on the FGR in the light of this analogy. In Sec. <ref> we recall the results obtained for  the Anderson model on the Bethe lattice within the FSA and analyze the differences between AL and MBL. In Sec. <ref> we present a strong disorder approximation 
for γ. In Sec. <ref> we study the signature of the presence of multifractal states in the anomalous scaling limit of the local spectral statistics obtained from exact diagonalizations of small systems. Finally, in Sec. <ref> we discuss  
the limitations of our interpretation and  in Sec. <ref> we describe some possible implications and perspectives for future investigations.
In App. <ref> we provide more details and supplemental information related to several points discussed in the main text.






§ THE MODELS AND HILBERT SPACE REPRESENTATION
 
We perform our analysis for three paradigmatic 1d models for MBL, 
namely the random-field Heisenberg XXZ spin chain, which has been used as a prototype for the MBL transition <cit.>, the  “Imbrie” model, for which the existence of the MBL transition has been proven rigorously <cit.>, 
and a model of interacting (spinless) fermions in a quasiperiodic (QP) potential <cit.>, similar to the one actually realized in cold atom experiments <cit.>. In the main text we will mostly focus on the disorder XXZ spin chain, although our results are valid for all the tree models, as shown in App. <ref>.
The Hamiltonian of the random-field XXZ S=1/2 chain is:

    H_ XXZ = ∑_i=1^n ( S_i^x S_i+1^x + S_i^y S_i+1^y +Δ S_i^z S_i+1^z + h_i S_i^z )   ,
 
with periodic boundary conditions and h_i independent and identically distributed uniformly in the interval [-h,h].
This model has been intensively studied <cit.> and its phase diagram is known for Δ=1, where a MBL transition takes place at a critical disorder within the interval h_c ∈ [3.7,4.5] in the middle of the many-body spectrum, E_0 ≈ - n /(4 n - 4), in the zero magnetization sector ∑_i=1^n S_i^z = 0, and for n ≲ 24 <cit.>.

By choosing as a basis of the HS the simultaneous eigenstates of the operators S_i^z, the Hamiltonian (<ref>) can be recast as a single-particle Anderson problem of the form

    H = ∑_a=1^ V E_a | a ⟩⟨ a | + t ∑_⟨ a, b ⟩( | a ⟩⟨ b | + | b ⟩⟨ a |  )   ,
 
where site orbitals represents many-body states in the spin configurations basis, | a ⟩ = | ↑↓↑⋯⟩, the constant hopping rate t=1/2 allows tunneling between states
a and b connected by the last term of (<ref>) which produces spin flips of two neighboring spins of opposite sign. The sums run over all V = nn/2 many-body configurations with zero magnetization. 
 The diagonal part of (<ref>)  yields the on-site random energies E_a = ⟨ a | ∑_i=1^n ( Δ S_i^z S_i^z+1 + h_i S_i^z)  | a ⟩, which are strongly correlated <cit.> (the V random energies are linear combination of n iid random numbers only).
Note that the spin configuration basis diagonalize H_ XXZ in the infinite disorder limit h →∞, where all many-body eigenfunctions are completely localized on single sites | a ⟩, and is thus suitable to study the stability of the insulating phase. The connectivity of the state | a ⟩ is equal to the number of domain walls 
(↑↓ or ↓↑) in the spin configuration, and ranges from 2 (for the n configurations with n/2 consecutive +1/2 spins and n/2 consecutive -1/2 spins) to n (for the two Neel states), with average value ⟨ z ⟩≈ (n+1)/2. Hence the network is sparse and high-dimensional, however, differently from the sparse random lattices, it is a deterministic graph with many regular local motifs and loops of all sizes (see Fig. <ref>).






§ THE FORWARD SCATTERING APPROXIMATION
 
A simple and powerful route to study and understand single-particle AL is by studying the convergence of perturbation theory starting from the insulator via the so-called the locator expansion <cit.>.
As shown by Anderson, in the insulating phase of the Anderson model on a d-dimensional lattice resonances do not proliferate at large distance in space and the locator expansion converges, implying that hopping only hybridize degrees of freedom within a finite volume of size ξ_ loc.
The FSA  consists in retaining only the lowest order terms in the locator expansion, which amounts in summing only over the amplitudes of the shortest paths connecting two points.
Within this approximation the non-local 
propagator (at energy E) between two points a and b at a given distance 
reads:

    G_a,b^ FSA (E) = ∑_∈_(a,b)∏_c ∈t/E - E_c  ,

	where ∈_(a,b) denotes a path among all Forward Scattering Paths connecting a and b.
Since at strong disorder the amplitude of each path decreases exponentially with its length, this approximation is expected to become more and more accurate as the disorder is increased.


In this paper we use the FSA to estimate the propagators G_0,nx for the XXZ disordered chain (<ref>) when recast as a single-particle tight binding problem (<ref>) in the spin configuration basis. We pick an infinite temperature many-body state | 0 ⟩ = {^0S_i^z 
} 
[with energy E_0 ≈ -n/(4n-4)] and determine the probability distribution of the matrix elements with configurations at distance nx from it (which differ by nx spins) at the lowest order in the “hopping” (i.e., the off-diagonal part of the many-body Hamiltonian in the spin configuration basis).

Of course one might argue that the FSA is a crude approximation for the true propagator. 
	However, on the one hand the FSA (i.e., calculating the Green's function by retaining only the lowest order in the off-diagonal terms) has been already successfully applied several times in the context of MBL,  yielding reasonably accurate  estimations of the boundaries of the insulating phase <cit.>, and providing an approximate strategy to construct the LIOMS <cit.>.
	Furthermore, since the rigorous results of <cit.> ensure that the perturbative expansion is in fact convergent in the MBL regime, taking only the leading order terms should not provide a too unreasonable starting point at least at strong enough disorder, as demonstrated by the recent quantitative analysis of Ref. <cit.>. However, since one of the main problem of the FSA is that the non-renormalized perturbative expansion has poles for any value of the energy within the support of the probability distribution of the random energies, in the following we will only focus on the typical value of the propagator, neglecting the effect of large matrix elements in the tails of the distribution. The limitations of this approach will be discussed further in Sec. <ref>.
	
The sum in (<ref>) over the _(0,nx) can be efficiently computed exactly using the transfer matrix technique described in <cit.>. In principle one should estimate the ergodic transition by requiring that any x-sector become ergodic. However for convenience in the following we only focus on the states {^0S_i^z 
 }  at zero overlap q = (1/n) ∑_i ^0S_i^z 
 · S_i^z=0 from the initial one (i.e., when half of the spins have flipped), as schematically depicted in Fig. <ref>, which is the largest sector on the HS. 
For the XXZ chain the shortest path to achieve q=0 (ignoring “loopy” terms in which spins are flipped  twice that contribute at higher order in perturbation theory) corresponds to Hamming distance n/4 on the graph.
The total number of configurations at distance n/4 and from | 0 ⟩ and in the same energy shell is N_n/4≈n/2n/4^2 Ω(E_0) δ E, where Ω(E_0) = e^S(E_0)∝ 1/√(n) is the many-body density of states defined in terms of the microcanonical entropy S(E_0) in the middle of the spectrum.

In Fig. <ref> we plot (the log of) the typical value of the matrix elements G_0,n/4^ FSA (computed within the FSA) as a function of (the log of) N_n/4  varying the system size from n=8 to n=28 (the whole probability distributions are shown in Fig. <ref> of App. <ref>). Different curves correspond to different values of the disorder strength across the MBL transition <cit.>.
The plot clearly shows that 

    | G_0,n/4^ FSA|_ typ= μ (  N_n/4)^-γ/2

at large n, 
with an exponent γ which increases as the disorder is increased (and μ of order 1).
The value of γ obtained by the linear fitting of ⟨log |  G_0,n/4^ FSA | ⟩ vs log N_n/4 (dashed lines of Fig. <ref>) is plotted in
Fig. <ref> as a function of the disorder h, showing that 
the disorder strength at which γ becomes larger than one (h_c ≈ 4.05) happens to be strikingly close to 
the critical disorder of the MBL transition determined by the most recent numerical works (for chains of about the same length of the ones considered here), h_c ∈ [3.7,4.5] <cit.>. Moreover, we find that in the limit of infinite disorder γ tends to 2. This exact same behavior is found for all the three models considered, as shown in Figs. <ref> and <ref> of App. <ref>.
The most intuitive way to interpret and rationalize these observations is by using an analogy with a very simple random matrix model, the RP model <cit.> and its generalizations <cit.>, which we detail in the next section.






§ ANALOGY WITH THE RP MODEL
 
The Hamiltonian of the RP model is a N × N symmetric matrix of the form <cit.>:

    H =  A + μ M/N^γ/2  ,

where A = a_i δ_ij is a diagonal matrix with iid entries drawn from some probability distribution p(a) (whose specific form is unimportant), M belongs to the Gaussian Orthogonal Ensemble (GOE) with unit variance (e.g., m_ij are iid random Gaussian variables with ⟨ m_ij⟩ = 0 and ⟨ m_ij^2 ⟩ = 1), and μ is a constant of O(1). 
The first term can be thought as the diagonal quenched on-site disorder, while M plays the role of the hopping that might create resonances between two (Poisson distributed) levels close in energy if |a_i - a_j| ≲ |m_ij|. The phase diagram of the RP model contains three different phases:
For γ > 2, standard second-order perturbation theory shows that the GOE term is a small regular perturbation, 
the Hamiltonian is close to A, and eigenstates are completely localized; Conversely, for γ < 1 the first term is a small regular perturbation, hence the rotationally invariant term M  dominates, the eigenstates are uniformly distributed on the unitary sphere, and are fully delocalized;
The regime 1 < γ < 2 is instead special, and provides an example of an extended nonergodic phase <cit.>: The eigenstates are supported over a large number of sites—hence they are delocalized—but only over a fraction N^1 - γ of them (D_2 = 2 - γ) which tends to zero in the thermodynamic limit—i.e., they are multifractal.

The transition taking place at γ_ AL = 2 corresponds to the standard AL and occurs when the amplitude Γ_ AL = N ⟨ |  H_ij | ⟩ vanishes in the thermodynamic limit. The physical interpretation of this criterion is that localization occurs when the number of sites j in resonance with a given site i stays finite for N →∞.
The transition from multifractal to ergodic eigenfunctions at γ_c = 1, instead, occurs when the amplitude

Γ_ ergo = N ⟨ H_ij^2⟩

diverges. 
This sufficient criterion for ergodicity has been proposed in Refs. <cit.> based on the idea that, using the FGR, the width Γ_ ergo essentially quantifies the escape rate of a particle created in i
(note however that the perturbative estimation of Γ_ ergo is valid as long as one can neglect the contribution of the off-diagonal elements to the density of states, i.e., γ > 1).
For γ<1 the width Γ_ ergo is much larger than the spreading of energy levels due to the disorder and the system is fully delocalized. For 1 < γ < 2 instead, Γ_ ergo vanishes as N^1 - γ in the thermodynamic limit, implying that eigenstates only occupy N^2 - γ sites close in energy.

Thus, adopting this analogy 
and using Γ_ ergo→ 0 as the criterion for ergodicity breaking based on the FGR, the FSA estimation of the effective exponents γ for the scaling of matrix elements between states at large distance in the HS, 
 
Eq. (<ref>) and Fig. <ref>), suggests that the MBL phase is similar to the intermediate phase of the RP model <cit.> (and its generalizations <cit.>) where Γ_ AL→∞ but Γ_ ergo→ 0, corresponding to delocalized but multifractal eigenstates: 
The initial state |0 ⟩ only hybridize with a sub-extensive fraction  
[  V e^S(E_0)]^1 - γ of the total configurations, 

thereby producing nonergodic multifractal wavefunctions which only occupy [ Ve^S(E_0)]^2 - γ states close in energy (D_2 ∼ 2 - γ) <cit.>.
In other words, the initial state |0 ⟩ is in resonances with many other states at large distance in the configuration space, but their number is not enough for the quantum dynamics to decorrelate in a finite time.
According to this interpretation, MBL in the HS does not correspond to AL, which instead would occur for γ >  2, when the number of resonances found from | 0 ⟩ stays finite in the thermodynamic limit <cit.>.





§ AL ON THE BETHE LATTICE WITHIN THE FSA
 
In order to elucidate the difference between the criterion Γ_ ergo→ 0 proposed above to detect MBL in the HS and standard AL in the limit of infinite dimensions, it is instructive to recall the paradigmatic case of the Anderson tight-binding model on the Bethe lattice <cit.>, described by the Hamiltonian (<ref>) with E_a iid in the interval [-W/2,W/2]. 
As for the MBL case, we determine the probability distribution of the propagator G_0,n^ FSA between a point 0 and one particular point at distance n from it within the FSA <cit.>.
Let us consider a Bethe lattice of branching number k (total connectivity k+1). Without loss of generality the site 0 can be considered as the root of the tree, and its energy can be set to E_0 = 0, in the middle of the band <cit.>.
Differently from the many-body problem, at the lowest order in perturbation theory there is only one path connecting two sites at distance n on the tree. According to Eq. (<ref>) one thus has that G_0,n^ FSA = ∏_a=1^n (t/E_a). One then immediately finds that on each single path the typical value of the matrix elements decays exponentially as ⟨log | G_0,n^ FSA| ⟩ = - n /ξ_ typ, over a typical length ξ_ typ^-1 = - log(2 e t/W). Since the number of sites at distance n from 0 (and in the same energy window) grows as N_n ≈ k^n ρ(0) δ E, one obtains that | G_0,n^ FSA|_ typ∝ ( N_n)^-γ/2, with an effective exponent γ = 2/(ξ_ typlog k). Adopting the ergodicity breaking criterion γ_c = 1 <cit.>, we get:

    W_ ergo^ FSA = 2 e t √(k)

(W_ ergo^ FSA/t ≈ 7.7 for k=2). However, such typical decay of the amplitude has nothing to do with Anderson localization, which is instead determined by the requirement that a particle created in 0 can escape on at least one of the k^n paths. The localization transition is thus obtained from the decay rate of the maximum amplitude among an exponential number of paths, which for large n is determined by the power law tails of the distribution of G_0,n^ FSA. This calculation yields the familiar result <cit.>

    W_ AL^ FSA≈ 2 e t k log k

(W_ AL^ FSA / t ≈ 29.1 for k=2, providing un upper bound for the true critical value W_ AL≈ 18.2 <cit.>), with a diverging localization length at the transition (with an exponent ν = 1), 
ξ_ loc = 1/[2 log (W/W_ AL)] ≫ξ_ typ. 

We argue that the threshold W_ ergo^ FSA signals the transition from ergodic to multifractal wavefunctions, which is a genuine phase transition (W_ ergo≈ 6.65 for k=2 <cit.>) only on (loop-less) Cayley trees <cit.> (and is related to the freezing glass transition of directed polymers in random media) and becomes a smooth crossover on the so-called random-regular graphs without boundaries (and loops whose typical size scales as the logarithm of the total number of sites of the graph), where full ergodicity is restored by loops larger than a characteristic correlation length which diverges at W_ AL <cit.>.

Note that in the limit of large connectivity, which is the relevant one for the many-body problem, W_ ergo and W_ AL have very different scaling with k, as √(k) and k log k respectively.
For a disordered XXZ chain of n interacting spins, the average connectivity of the HS grows as n/2, the effective width of the disorder is of order h √(n/12), and t=1/2.
Using the Bethe lattice estimation for W_ ergo^ FSA, Eq. (<ref>), one obtains a critical disorder h_c ≈ e √(6)≈ 6.7.
The scaling of W_ AL with the connectivity, Eq. (<ref>), would instead prohibit AL (note however the strong correlations of the potential between neighboring sites <cit.>).

Assuming that in the Hilbert space the typical value of the propagator decays exponentially with the distance over a typical length scale, | G_0,nx |_ typ∝ e^-n x / ξ_ typ, since the number of nodes in the Hilbert space at Hamming distance n/4 from |0⟩ is N_n/4 = n/2n/4^2 e^S(E_0)∼ 2^n/n^3/2, the amplitude 
	Γ_ ergo can be expressed as Γ_ ergo = | G_0,n/4 |_ typ^2   N_n/4∼ e^- n / (2 ξ_ typ) + n log (2) + O(log n). Hence, full delocalization and ergodicity occur when ξ_ typ≥ξ_ typ^c =  1/(2 log 2). This is in a certain sense the analogue of  the sufficient condition of delocalization obtained in <cit.> for the Anderson model on the BL, which correspond to the requirement that the exponential decay of  typical  correlations does not compensate anymore for the exponential proliferation of sites at large distance. Surprisingly enough, the  existence of such universal value of  ξ_ typ^c at the MBL transition predicted by   this simple argument is in perfect agreement with the recent results of the phenomenological RG approach of Ref. <cit.>.

This simple example clearly illustrates the differences that arise between AL on the Bethe lattice and MBL in the HS (at least within the FSA).
As mentioned above, random on-site effective energies of the many-body problem are not independent variables and are strongly correlated <cit.>. Moreover, 

the  number of forward-scattering paths connecting two many-body configurations grows factorially with the length of the paths, while on the Bethe lattice two points are connected by a unique path at the lowest order in t/W. However our analysis indicates that probably the most important difference consists in the fact that while AL on the Bethe lattice occurs when the number of resonances  found at large distance from a given site are finite in the thermodynamic limit (i.e., single-particle eigenstates occupy a finite volume in the thermodynamic limit), 
MBL instead takes place when the number of resonances found at distance nx from a given many-body state is still large but not enough to ensure e ergodicity and to allow the quantum dynamics to decorrelate from |0⟩ in a finite time (i.e., many-body eigenstates are extended but non-ergodic and occupy a sub-extensive portion [ Ve^S(E_0)]^D_2 of the  accessible volume in the HS). 
Hence, while single-particle AL on the BL is governed by the tails of the distribution P( G_0,n), MBL in the HS is governed by its bulk properties. 
Of course, focusing only on the typical value of the matrix elements, as done in Sec. <ref> is a drastic, and possibly wrong, assumption, since it neglects the effect of strong rare resonances that are known to play a very important role in MBL <cit.> (and probably increases the estimate for the critical disorder h_c <cit.>). We will come back to this issue in Sec. <ref> and in the concluding section, arguing that it might be corrected by mapping the MBL problem in the HS onto suitable generalizations of the RP model with power-law distributed off-diagonal elements <cit.>, and possibly refining the FSA computation by adding the self-energy corrections <cit.> in the denominators of (<ref>) and/or higher order terms of the locator expansion <cit.> in order to describe more accurately the tails of the distribution P( G_0,nx).



§ STRONG DISORDER APPROXIMATION
 
In this section we discuss how the effective exponent γ describing the scaling of the typical tunneling rates between two many-body states separated by nx spin flips with the number of configurations N_n x at distance nx from a given configuration in the HS, Eq. (<ref>), can be estimated analytically in the limit of strong disorder within the FSA.

A first very naive estimation 
might be obtained by recalling the intuitive argument given in the introduction for the origin of the multifractality of the many-body eigenstates in the HS, due to the presence of a finite density of active spins on the sites where the random fields is smaller than the energy required for spin flips <cit.>. For the disordered XXZ spin chain the energy needed to flip two neighboring spins with opposite sign, e.g. i and i+1, is Δ E = ± (h_i - h_i-1) + m Δ /2, with m=0,+1,-1 depending on how many domain walls 
have been created (annihilated) in the process. In the limit h ≫Δ the second term can be neglected, and Δ E is then a random variable of zero mean and variance 2 h^2/3. The density of active pairs of spins is thus proportional to the probability that this random variable is smaller than t=1/2, ρ∼√(3)/(2 √(2) h)

Since many-body states have typically n/2 domain walls, the volume occupied by many-body wavefunctions in this limit is 2^n/2ρ. Assuming, by analogy with the RP model, that D_2 ∼ρ/2 ∼ 2 - γ <cit.>, one obtains that 

    γ∼ 2 - √(3)/4 √(2) h  .

However this expression gives a very poor approximation of the numerical results plotted in Fig. <ref>, and overestimates γ by a large amount due to the fact that resonances and hybridization beyond the nearest neighboring spins are completely neglected.

A slightly more refined calculation can be performed as described below.
Let us consider a random initial state | 0 ⟩ = | ↑↓↑⋯⟩ with energy E_0 in the middle of the many-body spectrum.

In order to reach a state at zero overlap from it 
at the lowest order in the hopping we have to flip n/4 pair of spins in different locations of the chain.
Consider one particular path joining | 0 ⟩ with a given state at distance n/4 from it. The on-site random energies on the states visited along the path evolve as E_0 → E_0 + Δ E_1 → E_0 + Δ E_1 + Δ E_2 →… E_0 + ∑_i =1^n/4Δ E_i (recall that for ℓ large ∑_i =1^ℓΔ E_i is a Gaussian random variable of zero mean and variance 2 ℓ h^2 /3). The contribution to the sum (<ref>) coming from this specific path is then:

    1/2 Δ E_1×1/2(Δ E_1 + Δ E_2)×⋯×1/2 ∑_i =1^n/4Δ E_i  ,

which is a random variable.
All the (n/4)! permutations of the sequence by which spins are flipped yield a different path contributing to the same matrix elements between the same two states. 
In the following we assume that the contributions of different paths are uncorrelated and that the variance of the product (<ref>) on a single path is finite. Both assumptions are of course wrong. Yet, since we are only interested in the typical scaling of the matrix elements, they might give a reasonable approximation for the effective exponent γ.
The total number of configurations in the HS at distance n/4 from  | 0 ⟩ and in the same energy shell is N_n/4∝n/2n/4^2 Ω(E_0), with Ω(E_0) = e^S(E_0)∝ 1/√(n) being the many-body density of states in the middle of the many-body spectrum (and S(E_0) is the microcanonical entropy).
An approximate estimation of the exponent γ defined in Eq. (<ref>) is then given by:

    γn→∞≈1/2logn/4! - ⟨log(2 Δ E_1)⟩  - ⋯ - ⟨log(2 ∑_i=1^n/4Δ E_i)⟩/2 logn/2n/4 - 1/2log n  .


The prediction of Eq. (<ref>) is plotted in Fig. <ref>, showing a good agreement with the numerical results obtained by linear fitting of the data of Fig. <ref> at large n in the whole range [h_c,∞) (similar results are found also for the Imbrie model, Fig. <ref> of App. <ref>).





§ STATISTICS OF THE LOCAL DENSITY OF STATES
 
In this section we backup the perturbative results obtained within the FSA by probing the non-standard scaling limit of the local density of states (LDoS) which gives direct access to the nonergodic features of wavefunction statistics and provides an independent estimation of the fractal spectral dimension D_1.

The LDoS computed on a particular “site orbital” |a ⟩ of the HS is:

    Im G_aa (E_0 + i η) = ∑_β=1^ V | ψ_β (a) |^2 η/(E_β - E_0)^2 + η^2  ,

where |a ⟩ = | ↑↓↑⋯⟩, E_0 =  -n/(4n-4) is the energy in the middle of the many-body spectrum, V = nn/2 is the dimension of the HS, E_β are the eigenvalues of the Hamiltonian (<ref>), ψ_β (a) are the eigenfunctions' amplitudes expressed in the spin configuration basis, and η is an additional imaginary regulator.
The average value of this quantity gives the many-body density of states (1/π) ⟨ Tr Im G (E_0) ⟩ = e^S(E_0). 
In contrast, the typical value of Im G is controlled by the matrix element that couples a given site to the resonance sites at energy E_0.
Let us imagine a situation in which the sum in (<ref>) contains only [ V e^S(E_0)]^D_1 peaks of significant weight. Then the LDoS becomes smooth only if the broadening η exceeds the typical spacing between the peaks which, in the middle of the many-body spectrum, is typically of the order [ V /√(n)]^-D_1.
This implies that the typical value of the LDoS should exhibit the characteristic localized behavior (i.e., ⟨log Im G⟩∝η) up to a characteristic crossover scale η_⋆ much larger than the mean level spacing δ≈√(n)/ V.

In Fig. <ref> we plot the logarithm of the typical value of the LDoS, ⟨log Im G⟩, as a function of the imaginary regulator measured in units of the mean level spacing, η/δ, for several system sizes (n from 10 to 16) and for four values of the disorder strength h across the MBL transition. These plots are obtained by inverting exactly the many-body Hamiltonian (<ref>) in presence of the imaginary regulator, and averaging over several (about 2^40-2n) independent realizations of the disorder. The curves clearly show the existence of the crossover scale η_⋆ such that

    {[      Im G_ typ∝η        for η≪η_⋆; Im G_ typ∝η^ζ(h)        for η≫η_⋆ ] .

with an exponent ζ(h)>0 (except at very small h) which depends on the disorder (see inset of Fig. <ref>) but not on the system size. [Note that for completely AL eigenstates one should observe instead Im G_ typ∝η up to η of O(1).] Concretely, we have measured η_⋆ by performing linear fits of ⟨log Im G⟩ at small and large η with slope 1 and ζ(h) respectively, and determining where the two straight lines cross (orange squares of Fig. <ref>).
The crossover scale η_⋆ obtained by applying this procedure is plotted in Fig. <ref> as a function of the mean level spacing δ for several values of h.
In the MBL phase, h ≥ h_c, log (η_⋆ / δ) increases linearly by decreasing logδ (i.e., increasing n), consistently with the presence of multifractal eigenfunctions which only  occupy a subextensive part of the HS.
By fitting η_⋆ / δ∝δ^D_1 - 1 one obtains a 
measure of the fractal dimension D_1 which, by analogy with the RP model, gives a rough estimation of the effective exponent γ = 2 - D_1 <cit.> (blue squares in Fig. <ref>).
Instead on the metallic side, h < h_c, one observes a deviation from the straight line for the largest system sizes, signaling the recovery of a fully ergodic behavior.
In particular at small enough disorder (h=2.6, top-left panel of Fig. <ref>) we find that η_⋆≈δ, implying that D_1 = 1, as expected for fully ergodic eigenstates. Similar results are found also for the Imbrie model, Fig. <ref> of App. <ref>.
 
Note that the anomalous scaling limit of the spectral statistics in the multifractal regime 1 < γ < 2 of the RP model 
has been analyzed in full details <cit.>, and turns out to be slightly different from the one observed for the MBL system and shown in Fig. <ref>. In particular for the RP model one expects a region η_⋆ <η < η_ th where Im G_ typ is independent of η due to the presence of mini-bands of eigenfunctions close in energy and occupying N^D_1 sites. The Thoules energy η_ th is the width of these mini-bands and is simply obtained by multiplying the number of levels within a mini-band times the mean level spacing, i.e., η_ th/δ = N^D_1. Hence the typical value of the LDoS, when plotted as a function of η / δ should exhibit a flat part for N^D_1-1≪η/δ≪ N^D_1 which becomes broader and broader as the system size is increased.
The absence of such flat region in Fig. <ref> indicates that, differently from the RP model, for the many-body Hamiltonian (<ref>) the mini-band in the LDoS are multifractal <cit.>.  





§  LIMITATIONS OF THE FSA AND THE EXAMPLE OF TWO-DIMENSIONAL SYSTEMS
 

The main problem of the approach put forward in this work comes from the fact that the non-renormalized perturbative expansion has poles for any value of the energy within the support of the probability distribution of the random energies. In fact the non-renormalized perturbative expansion of the resolvent is always (i.e. with probability 1 in the thermodynamic limit) divergent even in the localized phase, due to local resonances, i.e. the sites in the expansion (<ref>)  where |E-E_c|<t, whose presence is inevitable in the thermodynamic limit however strong the disorder might be. Yet physically this is not a problem for localization. In fact the exact poles  of the Green's functions should be found at the eigenenergies of the Hamiltonian and not at the random on-site energies. As shown by Anderson <cit.>, this issue should be solved by re-summing the closed paths in the series expansion through the self-energy corrections. Since this re-summation cannot be done exactly on a generic lattice, in this work we have chosen to retain only the leading order terms of the perturbative expansion and to focus only on the scaling of the typical value of the propagator, which is only weakly affected by the presence of the poles.
Nonetheless, by doing so we are possibly overlooking  the effect of rare large amplitudes of the tunneling rates in the tails of the distribution (see e.g. Fig. <ref> of App. <ref>), 
whereas  rare delocalizing process (also called “thermal inclusions”) are known to play a  crucial role in the context of MBL.
In fact, in the latest years a phenomenological description for the many-body delocalization was proposed <cit.>, which relies on the “avalanche” instability, i.e., proliferation of an initial effectively thermal seed which grows until it swallows the whole system for h<h_c <cit.>.
The  avalanche picture predicts for instance that thermalization avalanches should destabilize the MBL phase in any dimension larger than 1 <cit.>.


On the one hand, the ergodicity breaking criterion obtained from the typical decay of the matrix elements, N_n x | G_0,nx|_ typ^2 → 0, 
yields a critical disorder which is in strikingly good agreement with the most recent numerical results obtained from exact diagonalizations of chains of about the same length of the ones considered here <cit.>.
On the other hand, however,  
recent results obtained in by applying time-dependent variational principle to matrix product states that allow to study 1d chains of a length up to n=100 indicate a substantial increase of the estimate for the critical disorder (h_c ≳ 5.5) that separates the ergodic and many-body localized regimes <cit.>. Such enhancement of ergodicity in large systems is likely to be due to  the  existence  of rare 
delocalizing processes <cit.>—possibly involving degrees of freedom distant in real space—that are not typically present in small systems and that are not detectable by typical amplitudes.



In order to have a concrete benchmark example of the possible limitations of our approach, in this section we apply the ideas discussed above  to a two-dimensional setting, where thermalization avalanches are expected to destabilize the insulating phase in the thermodynamic limit. More specifically, 
we consider hard-core bosons with nearest-neighbor interactions on a quasi-1d square lattice of length L and width s and total number of sites n=Ls <cit.>. The Hamiltonian is given by:

    H_  HCB = ∑_⟨ i, j ⟩[ - J/2( b_i^† b_j + h.c.) + U n_i n_j ] - ∑_i ϵ_i n_i   ,

where b^†_i creates a boson at site i, n_i  = b^†_i b_i (the occupation of each site is restricted to n_i = 0,1), and 
the summation over ⟨ i, j ⟩ couples neighboring sites on the strip. The on-site potentials ϵ_i are iid random variables taken from a uniform distribution on [-h,h ]. We set J=1 and U=1 in the half-filling sector, ∑_i n_i = n/2 .

For s=1 this model is equivalent via a Jordan-Wigner transformation to the Heisenberg XXZ random-field spin chain (<ref>) considered above, while upon increasing the width s of the strip we move towards a 2d geometry. Below we repeat the same analysis described in Sec. <ref> for this model with s=2 and s=4. 
We choose as a basis in the HS the tensor product of the simultaneous eigenstates of the number operators n_i, 
| a ⟩ = | ∙∘∙⋯⟩, 
such that for J=0 the many-body eigenstates are perfectly localized.
We pick an infinite temperature many-body state | 0 ⟩ at random with energy close to the middle of the many-body
spectrum, and we evaluated the effective tunneling rates between such state and states at  Hamming distance n/4 in the HS, i.e., such that half of the bosons have moved compared to the initial configuration. The computation of the propagator  is again performed at the level of the FSA, i.e., retaining only the ∼ (n/4)! leading-order terms in the perturbative expansion (i.e., only the shortest paths in the HS).
By comparing the scaling of the typical value of the matrix elements | G^ FSA_0,n/4|_ typ with the number of accessible nodes of the HS at Hamming distance n/4 from | 0 ⟩ and in the same energy shell, N_n/4 = n/2n/4^2 e^S(E_0), we obtain the effective exponent γ plotted in Fig. <ref> for s=1 (same curve as in Fig. <ref>), s=2 and s=4. 

Applying the  ergodicity breaking criterion N_n/4 | G_0,n/4|_ typ^2 → 0 inspired by the RP model <cit.> and its generalizations <cit.>, we obtain that the critical disorder h_c where γ becomes larger than one does systematically (possibly exponentially) increase with the width of the strip s, as predicted by the avalanche approach. Yet, such increase is much weaker than the one recently obtained in <cit.>  by applying time-dependent variational principle to very large quasi-1d strips, and is also much weaker than the analytical prediction of the effect of  avalanches in 2d, h_c ∝ 2^s <cit.>. 

It is also instructive to study how the whole probability distribution of the propagator is modified upon increasing s (at the same disorder strength in the MBL regime, h=6, and for the same total number of sites n=Ls =24), as plotted in Fig. <ref> of App. <ref>. One clearly observes a strong enhancement of the tails of the distribution when s is increased, corresponding to rare large tunneling amplitudes, which, however, have only a moderate effect on the typical value.

All in all, this analysis indicates that the  approach put forward in this work is  able to capture some some mild signature of  the avalanche instability. Yet, focusing only on the typical value of the propagator evaluated at the leading order of the perturbative expansion does not allow to recover the full correct quantitative behavior, especially in situations in which the thermal inclusions are expected to have a strong effect (i.e., very large systems and/or d>1). 



§ DISCUSSION AND CONCLUSIONS
 
In this paper we have proposed  a novel perspective to analyze the properties of the MBL transition in the HS <cit.>.
We evaluated the tunneling rates G_0,nx^ FSA between two many-body states at extensive distance nx (when a finite fraction x of the spin are flipped) at the lowest order in perturbation theory starting from the insulator, and compared their (typical) amplitude to the number of accessible configurations N_n x at distance nx from a given initial state (and in the same energy shell). 
We have shown that in the MBL phase, although typically many resonances are formed, they are not enough to allow the quantum dynamics to decorrelate from the initial condition in a finite time.
Concretely, we have put forward a 
criterion for ergodicity breaking 
based on the ideas of Refs. <cit.> and on the FGR. This criterion 
is much weaker than requiring AL of many-body eigenstates in the HS and
suggests that the MBL transition takes place when the amplitude Γ_ ergo =  N_n x | G_0,nx|_ typ^2 → 0 <cit.>. 
This implies that many-body eigenfunctions in the HS are delocalized but multifractal, and typically only occupy a subxtensive portion [  V e^S(E_0)]^D_2 of the total accessible configurations, with 0<D_2<1 and D_2 → 0 only in the limit of infinite disorder. 
MBL in the HS is thus reminiscent of the transition from ergodic to multifractal states of the RP matrix ensemble <cit.> (and its generalizations <cit.>), 
and not to AL in the limit of infinite dimension <cit.>, which instead occurs when the number of resonances found from a given configuration stays finite in the thermodynamic limit <cit.> (in fact we find that many-body eigenstates become AL only in the limit of infinite disorder, as expected from intuitive arguments <cit.>).

Our interpretation fully supports the picture recently proposed in Ref. <cit.> where MBL is seen as a fragmentation of the HS, as well as similar ideas promoted to explore the analogy between MBL  and a percolation transition in the configuration space <cit.>.
It is also in agreement with the most recent numerical results on the wavefunctions' statistics <cit.>, with perturbative calculations <cit.>, and with intuitive arguments that strongly indicate that the many-body eigenstates are multifractal in the whole MBL regime.  

The approach presented in this paper has several advantages. On the one hand, it provides a clear view of the MBL transition in the HS, which is conceptually of prime interest and gives a transparent explanation of the difference between AL and MBL; On the other hand, it yields a simple and quantitatively predictive tool to estimate the critical disorder and the properties of the eigenstates in the MBL phase, since the transfer matrix algorithm <cit.>  used to determine the amplitude of the propagator (<ref>) is computationally much easier than exact diagonalizations and allows to investigate larger system sizes.

In the following we discuss several possible  limitations and implications of our approach, as well as some perspectives for future work.

(1)As already discussed in details in Sec. <ref>, in this work we only focused on the asymptotic scaling behavior of the typical value of the amplitude of the propagator evaluated at the lowest order of the perturbative expansion. 
In this way we are clearly overlooking the effect of strong rare resonances in the tails of the distributions of the propagators (see Fig. <ref>). By doing so we find that the ergodicity breaking criterion built on the FGR, γ_c = 1, yields a critical disorder which is in strikingly good agreement with the most recent numerical results obtained from exact diagonalizations of chains of about the same length of the ones considered here <cit.>.
However, 

our approach is not able to capture the enhancement of ergodicity observed  in very large chains <cit.> and in two-dimensional systems <cit.> which is likely to be due to  the  existence  of rare 
delocalizing processes <cit.>—possibly involving degrees of freedom distant in real space—that are not typically present in small systems and that are not detectable by typical amplitudes.
It would be therefore helpful to go beyond the FSA either including higher order terms in the perturbative expansion (see e.g. <cit.> for a recent attempt in this direction), or developing  some approximate treatment for the  self-energy corrections in the denominators of (<ref>), as, for instance, recently proposed in Refs. <cit.>. This might  allow one to  describe more accurately the tails of the distribution P( G_0,nx) and to  take into account the effect of rare resonances which are known to play a very important role in the context of MBL <cit.>.

(2) A tightly related issue is that the RP model is certainly oversimplified: The mini-band in the local spectrum is not multifractal, the spectrum of fractal dimension is degenerate, and strong resonances are absent (see above).
In fact, we have already noticed in Sec. <ref> that the typical value of the LDoS of the many-body Hamiltonian, Fig. <ref>, behaves differently from the one of the RP model <cit.>. 
In this respect it would be useful to study suitable generalizations of the RP ensemble with broadly distributed off-diagonal elements (see, e.g., <cit.>) that might provide a better effective description for 
the MBL transition in the HS. 
It would be also desirable to complete the present computation by studying the asymptotic scaling of the propagator in all the x-sectors, which might  allow one to obtain a more precise estimation of the effective exponent γ.

(3) Another related problem is that the ergodicity breaking criterion used here, together with the mapping onto the RP model, seem to predict that the spectral dimension D_2 ∼ 2 - γ is continuous at the MBL transition (i.e., γ→ 1 for h → h_c^+), while recent numerical results <cit.> as well as theoretical arguments <cit.> indicate a discontinuous jump of D_q at h_c. 
It would be interesting to understand whether going beyond the FSA by including higher order corrections to Eq. (<ref>) and/or considering generalizations of the RP ensemble with broadly distributed off-diagonal matrix elements  as effective descriptions of the MBL transition in the HS can lead to a scenario in which the fractal dimensions exhibit a discontinuous jump at the transition. 

(4) Another important aspect concerns the implications of the interpretation proposed here on the unusual properties of the bad metal phase preceding MBL <cit.>.
In some recent works it was in fact suggested that the subdiffusive transport and the anomalously slow out-of-equilibrium relaxation observed in numerical simulations and experiments might be explained in terms of the apparent nonergodic features of many-body wavefunctions in the HS <cit.>, while the interpretation proposed here and recent numerical results <cit.> indicates that eigenstates of (<ref>) become fully ergodic on the metallic side of the transition. This issue might be explained in terms of strong finite-size effects. 
The scaling analysis of <cit.> indicates indeed that for h<h_c the finite-size effects controlling the asymptotic scaling behavior of the inverse participation ratios Υ_q are dominated by a nonergodic volume which diverges very fast as the transition is approached (see also <cit.>). 
As a result, many-body wavefunctions might behave as if they were multifractal in a very broad range of system sizes even before h_c, especially in the region in which γ is close to 1 and the spectral band-width associated to the off-diagonal tunneling rates is of the same order of the spreading of the energy levels due to the disorder, thereby producing anomalous diffusion and slow power-law relaxation of physical correlations on a very large time-window spanning many decades.



(5) This problem is in fact tightly related 
to the critical properties of the MBL transition.
It is well known that the FSA yields the mean-field exponent ν = 1 for the divergence of the localization length at the Anderson transition irrespectively of the dimension and of the structure of the underlying graph <cit.>.
The same exponent ν=1 governs the transition from ergodic to multifractal eigenstates of the RP model taking place at γ_c = 1 <cit.>, while in the critical exponent for the generalized RP ensemble with log-normal distributed off-diagonal elements was recently found to vary between 1 and 1/2 depending on the parameters <cit.>. However these critical behaviors are not compatible with the  
 the recent phenomenological RG studies for the MBL transition <cit.> which instead predict a KT-like scenario with an exponential divergence of the localization length. 



(6) A promising direction for future research would be to exploit the simplicity of our approach to address important questions such as 
the stability of MBL with respect to rare thermal inclusions of weak disorder that occur naturally inside an insulator and that may trigger a thermalization avalanche <cit.>. 
It would be interesting for instance to insert by hand large but finite ergodic bubbles of weak disorder in the 1d XXZ spin chains and investigate the signature of quantum avalanches <cit.> by studying the effect of these bubbles on the scaling behavior of the amplitude of the propagators.

(7) Finally, we would like to comment on the possible implications of our results on the recent debate on quantum chaos vs MBL <cit.>.
In fact a recent paper <cit.> has claimed that MBL is not a phase of matter, but rather a finite-size regime that yields to ergodic behavior in the thermodynamic limit. This conclusion was reached on the basis of a finite-size-scaling analysis of small 1d spin models using diagnostics from quantum chaos that probe the statistics of level spacing only, such as the structure form factor (SSF) S(τ) and the average level spacing ratio ⟨ r ⟩. 
In the light of the new interpretation proposed here, it is instructive to recall the known results for the statistics of eigenvalues of the RP model <cit.>:
In the intermediate regime 1 < γ <2 of delocalized but nonergodic wavefunctions, although the average DoS asymptotically converges to the distribution of the  diagonal energies p(a) (and not to the Wigner semicircle), the nearest-neighbor level statistics is described by the Wigner-Dyson statistics <cit.>. In particular the (unfolded) spectral form factor was shown to be universal, i.e. independent of the specific form of p(a), and to converge to the Wigner-Dyson form for γ < 2, and to Poisson only for γ > 2 <cit.>. Similarly, ⟨ r ⟩ approach the GOE universal value ⟨ r ⟩≈ 0.53 in the thermodynamic limit for 1 < γ < 2 and the Poisson one ⟨ r ⟩≈ 0.38 for γ > 2. In fact the N^D_1 states close in energy inside each mini-band exhibit level repulsion and GOE-like correlations. The crossover from GEO-like behavior to Poisson statistics occurs thus on the scale of the Thouless energy η_ th = N^D_1×δ = N^1-γ, which is vanishingly small for N →∞ but is still much larger than the mean level spacing.

As discussed in Sec. <ref>, the behavior of the spectral statistics of the many-body problem exhibits several important differences with respect to the RP model. The fact that the crossover scale η_⋆ below which one observes  the characteristic localized behavior (see fig. <ref>)  is much larger than the mean level spacing indicates that  consecutive energy levels are not hybridized by the off-diagonal perturbation, whose effect only sets in  on an larger energy scale ≫δ, thereby implying that the level statistics on the scale of the mean level spacing should be of Poisson type.

Moreover, differently from the RP model, for the many-body Hamiltonian (<ref>) the mini-band in the LDoS are multifractal.
Yet,  if one takes the mapping to the RP model seriously beyond the qualitative level, one might be tempted to argue that any observable related to the level statistics on the scale of the mean level spacing only might be uninformative on the existence of a MBL transition in the thermodynamic limit. In order to illustrate this idea, we have measured the spectral form factor (SFF) defined as <cit.>:

    S(τ) = 1/Z⟨|∑_β = 1^ V g(Ẽ_β) e^- i Ẽ_β t|^2 ⟩  ,

where Ẽ_β are the unfolded eigenvalues of (<ref>), such that ⟨Ẽ_β + 1 - Ẽ_β⟩ = 1, g(Ẽ_β)=exp[- (Ẽ_β - Ẽ_0)^2/(2 λ^2 σ^2_Ẽ)] is a Gaussian filter (λ = 0.2 is a dimensionless parameter that controls the effective fraction of eigenstates included in the S(τ), Ẽ_0 and σ^2_Ẽ are the average energy and the variance of the unfolded eigenvalues, respectively, for a given disorder realization), and the normalization Z = ⟨∑_β| g(Ẽ_β) |^2 ⟩ is such that S(τ) → 1 for τ→∞. The numerical results for S(τ) 
computed from exact diagonalizations of small (n=10, 12, and 14) XXZ chains (<ref>)  and for h=6, deep into the MBL phase, are shown in Fig. <ref>.
Inspired by the analogy with the RP random matrix ensemble discussed above, we benchmark these results with the SFF obtained for the RP model (<ref>) with parameters chosen in such a way to mimic as closely as possible the interacting one: We set N = nn/2 (equal to the size of the HS of the interacting model), and consider a Gaussian distribution of diagonal energies p(a) of variance ∝ n and equal to the variance of the random energies E_i in the spin configuration basis of the XXZ spin chain; We set γ=1.35 (which is approximately the value of the exponent found numerically at h=6, see Fig. <ref>), and μ = 0.1 (which is the value of the prefactor in Eq. (<ref>) obtained by fitting the numerical data of Fig. <ref> at h=6).
The SSFs computed for the RP model (dashed lines in Fig. <ref>) turn out to be very similar to the ones of the interacting model at the same value of n, although we know rigorously that they should approach the GOE result 2 τ - τlog(1 + 2 τ) (black line) for N →∞. However, as shown by Fig. <ref>, the convergence to the GOE asymptotic result is very slow and finite-size effect are very big at finite N due to the fact that the Thouless energy is still too close to the  mean level spacing.
All in all, this analysis suggests that one should take extreme caution when using diagnostics for the MBL transition based on the statistics of energy levels on the scale of the mean level spacing only, as they might slowly drift towards a GOE-like behavior in the thermodynamic limit even deep inside the MBL phase <cit.>.
This might also explain why the critical disorder estimated from the crossing of the curves of the average level spacing ratio ⟨ r ⟩ at different system sizes drifts to larger h when increasing n <cit.>.


	I would like to thank  G. Biroli, D. Facoetti, I. V. Gornyi,  I. Khaymovich, G. Lemarié, D. J. Luitz, A. D. Mirlin, M. Schiró, and V. Ros for many enlightening and helpful discussions.






§ RESULTS FOR THE IMBRIE MODEL AND FOR INTERACTING FERMIONS IN A QP POTENTIAL
 

In this appendix we provide more details and supplemental information related to several points discussed in the main text. In particular
we present the results obtained for two other models for the MBL transition described below.






 §.§ The models


The  “Imbrie” model is defined by the following Hamiltonian:

    H_I = ∑_i=1^n ( Δ_i σ_i^z σ_i+1^z + h_i σ_i^z + Γ_i σ_i^x )   .

We follow Ref. <cit.> and set Γ_i=1, and Δ_i and h_i uniformly distributed in Δ_i ∈ [0.8,1.2] and h_i ∈ [-h,h].
The existence of the MBL transition has been proven rigorously for this model under the minimal assumption of absence of level attraction <cit.>.
The numerical results of <cit.> seem to indicate that the critical disorder strength should be in the interval h_c ∈ [3.75,4.5].

We also studied a one-dimensional model of spinless fermions on a QP lattice <cit.>:

    H_ QP = ∑_i=1^n [ t ( c_i^† c_i+1 +  h.c.) + Δ n_i n_i+1 + V_i n_i  ]   ,

where V_i is a QP potential of the form:

    V_i = h cos (2 πω i + ϕ)   ,


and ϕ is a random phase.
(Note that H_ QP exactly maps to the XXZ spin chain in a QP magnetic field under a Jordan-Wigner transformation.)
This model is very similar to the one realized in cold atom experiments <cit.>.
As in Ref. <cit.>, we set t=1, the irrational number ω to be the inverse of the golden mean ω = (√(5) - 1)/2, Δ = 1, and only consider the half-filling sector ∑_i n_i = n/2. For this choice of the parameters previous studies <cit.> have established the presence of a MBL transition with a critical strength of the QP potential in the interval h_c ∈ [3.5,4.5].
For both models we consider periodic boundary conditions.



 §.§ Scaling of the matrix elements within the FSA

In the following we discuss the results obtained for these two models by applying the same analysis described in Sec. <ref> of the main text. Concretely, we compute the probability distributions of the propagators
G_0,nx^ FSA when the many-body systems are recast as single-particle tight binding problems (<ref>). 



By choosing the spin configuration basis, the HS of (<ref>) is a
n-dimensional hypercube of V = 2^n sites (the total magnetization is not conserved by H_I).
Each configuration | a ⟩ = | ↑↓↑⋯⟩ of n spins corresponds to a corner  
of the  hypercube  by  considering σ_i^z = ±1 as  the  top/bottom  face  of  the  cube's n-th  dimension. 
The random part of the Hamiltonian is by definition diagonal on this basis, and gives correlated 
random energies on each site orbital of the hypercube, E_a = ⟨ a | ∑_i=1^n [ Δ_i σ_i^z σ_i^z+1 + h_i σ_i^z] | a ⟩.
The interacting part of H_I acts as single spin flips on the configurations {σ_i^z }, and plays the role the hopping rates connecting “neighboring” sites in the configuration space (with t=1).
At Γ = 0 the many-body eigenstates of (<ref>) are simply product states of the form |σ_1^z ⟩⊗|σ_2^z ⟩⊗⋯⊗|σ_n^z ⟩, and the system is fully localized in this basis.

Similarly, for the QP model we choose as a basis the tensor product of the simultaneous eigenstates of the number operators n_i, 
| a ⟩ = | ∙∘∙⋯⟩, 
such that for t=0 the many-body eigenstates are perfectly localized. The HS of (<ref>) is then represented by the same graph as for the XXZ random-field Heisenberg model considered in the main text, and its size is V = nn/2. The diagonal part of (<ref>)  yields the on-site quasi-random energies E_a = ⟨ a | ∑_i=1^n [ Δ n_i n_i+1 + V_i n_i ] | a ⟩, while the interacting part 
allows tunneling between “neighboring” configurations with hopping rate t=1.
The connectivity of the state | a ⟩ is equal to the number of pairs ∙∘ or ∘∙ 
in the state, and ranges from 2 to n with average value ⟨ z ⟩≈ (n+1)/2.

For both model we pick an infinite temperature many-body state | 0 ⟩  with energy E_0 in the middle of the spectrum, 
and determine the probability distribution of the matrix elements with configurations at distance nx from it  at the lowest order in the hopping.
As for the XXZ model presented in the main text, instead of considering all x-sectors separately, for simplicity we only focus on the states at zero overlap 

from the initial one, i.e. when half of the spins have flipped or half of the particles have moved respectively. (Specifically the overlap is defined as q = (1/n) ∑_i ^0σ_i^z ·σ_i^z and q = (1/n) ∑_i (n_i^0 - 1/2) (n_i - 1/2) for the two models, where the random initial state is denoted as | 0 ⟩ = {^0σ_i^z  } and | 0 ⟩ = { n_i^0  } respectively.) 
The shortest path to achieve q=0 corresponds to Hamming distance n/2 on the hypercube for the Imbrie model and to Hamming distance n/4 on the graph for the QP model.
The total number of configurations at q=0 from | 0 ⟩ and in the same energy shell is 
N_n/2≈nn/2 e^S(0)δ E for the Imbrie model and N_n/4≈n/2n/4^2 e^S(E_0)δ E for the QP model, where S is the microcanonical entropy in the middle of the spectrum, defined as the logarithm of the number of states at that energy.



In Fig. <ref> we plot (the log of) the typical value of the propagator (computed within the FSA) as a function of (the log of) N_q=0 for the two models varying the system size. Different curves correspond to different values of the disorder strength across the MBL transition <cit.>.
The plots clearly shows that Eq. (<ref>) holds. 
The value of the effective exponent γ is obtained by the linear fitting of the numerical results at large N_q=0 (dashed lines of Fig. <ref>), and is plotted in
Fig. <ref>. This figure shows that for both models
the critical disorder determined by the most recent numerical works is consistent with the value of the disorder such that γ becomes larger than one, while γ tends to 2 in the limit of infinite disorder. This is exactly the same behavior found for the XXZ spin chain and discussed in the main text, Figs. <ref> and <ref>, supporting the robustness of our conclusions and the validity of the criterion  built on the FGR, Γ_ ergo =  N_q=0 | G_q=0|_ typ^2 → 0, for the MBL transition.

Note that for the Imbrie model one can repeat the strong disorder approximation discussed in Sec. <ref>, and straightforwardly obtain the following expression:

    γn→∞≈1/2logn/2! - ⟨log(2 h_1)⟩  - ⋯ - ⟨log(∑_i=1^n/2 2 h_i)⟩/lognn/2 - 1/2log n  ,

which is in good agreement with the numerical results (red dashed curve in the top panel of Fig. <ref>).
The Bethe lattice estimation of the ergodicity breaking transition based on the exponential decay of the matrix element on a single branch, Eq. (<ref>), yield instead h_c = 2 √(3) e ≈ 9.4.



In Fig. <ref> we show the full probability distributions of the amplitude of the tunneling rates for the three models considered in this work and for several value of the disorder strength across the MBL transition.
The distributions have power-law tails P(| G_q=0^ FSA|) ∼ | G_q=0^ FSA|^-(1+μ) as expected for any value of E in Eq. (<ref>) within the support of the probability distribution of the random energies E_c. The exponent μ varies from one model to another, but does not depend (or depends very weakly) on the disorder strength and on the system size.

Finally, in Fig. <ref> we  plot the probability distribution of the propagator P(| G_0,n/4^ FSA|)  for the quasi-1d model of hard-core bosons with nearest neighbor interactions described by the Hamiltonian (<ref>), obtained when the width of the strip from s=1 to s=4. The curves correspond to disorder strength h=6, inside the MBL regime,  and for the same total number of sites n=Ls =24). One clearly observes a strong enhancement of the tails of the distribution when s is increased, corresponding to rare large tunneling amplitudes, accompanied by a moderate increase of the typical value.





 §.§ Signature of  eigenstates' multifractality  in the spectral statistics of the Imbrie model
 
As done in Sec. <ref> for the XXZ random-field spin chain, one can investigate the signatures of the multifractality of the eigestates from the unusual scaling limit of the spectral statistics of the Imbrie model.
In the top panels of Fig. <ref> we plot the logarithm of the typical value of the LDoS, ⟨log Im G⟩, as a function of the imaginary regulator measured in units of the mean level spacing, η/δ, for four system sizes (n from 8 to 14) and for two values of the disorder strength h across the MBL transition. These plots are obtained by inverting exactly the many-body Hamiltonian (<ref>) in presence of the imaginary regulator, and averaging over several (about 2^36-2n) independent realizations of the disorder. The curves show the existence of the crossover scale η_⋆ separating the behavior of Im G at small and large η, as described by Eq. (<ref>). We have extracted such crossover scale from the data by applying the procedure described in Sec. <ref> of the main text. In the bottom-left panel of Fig. <ref>
 we plot log(η_⋆ / δ)  as a function of δ for several values of h ≳ h_c.
We observe that log (η_⋆ / δ) increases linearly by decreasing logδ (i.e., increasing n), consistently with the presence of multifractal eigenfunctions which only  occupy a subextensive part of the HS.
By fitting η_⋆ / δ∝δ^D_1 - 1 one obtains an independent estimation of the fractal dimension D_1 and, by analogy with the RP model, of the effective exponent γ = 2 - D_1 <cit.> (blue squares in Fig. <ref>).
For the Imbrie model, however, such numerical estimation of γ 
does not agree well with the numerical results of the FSA. This discrepancy is possibly due to strong finite-size effects, since simple intuitive arguments suggest that one should find γ∼ 2 - c/h at strong disorder.

99

anderson P. W. Anderson, Phys. Rev. 109, 1492 (1958).

BAA D. M. Basko, I. L. Aleiner, and B. L. Altshuler, Annals of Physics 321, 1126 (2006).

Gornyi I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Phys. Rev. Lett. 95, 206603 (2005).

reviewMBL E. Altman and R. Vosk, Annu. Rev. Condens. Matter Phys. 6, 383 (2015).

reviewMBL2
R. Nandkishore and D. A. Huse, Annu. Rev. Condens. Matter Phys. 6, 15 (2015).

reviewMBL3 D. A. Abanin and Z. Papić,  Annalen  der  Physik 529, 1700169 (2017).

reviewMBL4
	F. Alet and N. Laflorencie, Comptes Rendus Physique 19, 498 (2018).
	
reviewMBL5  D.  A.  Abanin,  E.  Altman,  I.  Bloch, and  M. Serbyn, Rev. Mod. Phys. 91, 021001 (2019).

Huse V. Oganesyan and D. A. Huse, Phys. Rev. B 75, 155111 (2007). 

pal A. Pal and D. A. Huse, Phys. Rev. B 82, 174411 (2010).

alet D. J. Luitz, N. Laflorencie, and F. Alet, Phys. Rev. B 93, 060201(R) (2016).

experiments1 M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and 
I. Bloch, Science 349, 842 (2015).

experiments2 P. Bordia, H. P. Lüschen, S. S. Hodgman,  M. Schreiber, I. Bloch,  and U. Schneider, 
Phys. Rev. Lett. 116, 140401 (2016).

experiments3 J.-Y. Choi, S. Hild,   J.  Zeiher,   P. Schauss, A. Rubio-Abadal, T. Yefsah, V. Khemani,  D.  A.  Huse,  I.  Bloch,    and  C. Gross, Science 352, 1547 (2016).

experiments4 J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, Nat. Phys. 12, 907 (2016). 

experiments5 G. Kucsko, S. Choi, J. Choi, P. C. Maurer, H. Sumiya, S. Onoda, J. Isoya, F. Jelezko, E. Demler, N. Y. Yao, and M. D. Lukin, arXiv:1609.08216. 

experiments6 K. Xu et  al., Phys.  Rev.  Lett. 120, 050507 (2018).

LIOMS J. Z. Imbrie, Phys. Rev. Lett. 117 027201 (2016).

LIOMSb
	J. Z. Imbrie, V. Ros, and A. Scardicchio, Annalen der Physik 1600278 (2017).

LIOMS1 M. Serbyn, Z. Papić, and D. A. Abanin, Phys. Rev. Lett. 111, 127201 (2013); Phys. Rev. B 90, 174302 (2014).

LIOMS2 V. Ros, M. Müller, and A. Scardicchio, Nuclear Physics B 891, 420 (2015).
	
ETH M. Srednicki, Phys. Rev. E 50, 888 (1994); M. Rigol, V. Dunjko, and M. Olshanii, Nature 481, 224 (2012).

violation I. L. Aleiner, B. L. Altshuler, and G. V. Shlyapnikov, Nat. Phys. 6, 900 (2010); D. A. Huse, R. Nandkishore, V.Oganesyan, A. Pal, and S. L. Sondhi, Phys. Rev. B 88, 014206 (2013).

entanglement J.  H.  Bardarson,  F.  Pollmann,  and  J.  E.  Moore,  Phys. Rev. Lett. 109, 017202 (2012).

entanglement1 R. Vosk and E. Altman, Phys.  Rev.  Lett. 110,  067204,  (2013).  

critical  R. Vosk,  D. A. Huse,  and E. Altman,  Phys. Rev. X 5, 031032 (2015); A. C. Potter, R. Vasseur, and S. A. Parameswaran, Phys.Rev. X 5, 031033 (2015);  
P.  T.  Dumitrescu,  R.  Vasseur  and  A.  C.  Potter,  Phys. Rev. Lett. 119, 110604 (2017).

thiery T. Thiery, F. Huveneers, M. Müller, and W. De Roeck, Phys. Rev. Lett. 121, 140601 (2018); T. Thiery, M. Müller, W. De Roeck, arXiv:1711.09880

KT A. Goremykina, R. Vasseur, and M. Serbyn, Phys. Rev. Lett. 122, 040601 (2019); A. Morningstar and D. A. Huse, Phys. Rev. B 99, 224205 (2019);  

KT1 P. T. Dumitrescu, A. Goremykina, S. A. Parameswaran, M.  Serbyn,  and  R.  Vasseur,  Phys.  Rev.  B 99,  094205 (2019).

mirlin K. S. Tikhonov, A. D. Mirlin, M. A. Skvortsov,
Phys. Rev. B 94, 220203(R) (2016); K. S. Tikhonov and A. D. Mirlin, Phys. Rev. B 99, 024202 (2019).

avalanches W. De Roeck and F. Huveneers, Phys. Rev. B 95, 155129 (2017);  
D.  J.  Luitz,  F.  Huveneers,  and  W.  De  Roeck, Phys. Rev. Lett. 119, 150602 (2017).

gopala S. Gopalakrishnan and David A. Huse, 
Phys. Rev. B 99, 134305 (2019).

doggen E. V. H. Doggen, I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, arXiv:2002.07635

bad_metal1  D. J. Luitz and Y. Bar Lev, Ann. Phys. (NY) 529, 1600350 (2017).

bad_metal2 K. Agarwal, E. Altman, E. Demler, S. Gopalakrishnan, D. A. Huse, and M. Knap, Ann. Phys. (NY) 529, 1600326 (2017).

dot B. L. Altshuler, Y. Gefen, A. Kamenev, L. S. Levitov, Phys. Rev. Lett. 78, 2803 (1997).

spectral_diffusion I. V. Gornyi, A. D. Mirlin, D. G. Polyakov, and A. L. Burin, Annalen der Physik 529, 1600360 (2017).

abou R Abou-Chacra, P. W. Anderson, and D. J. Thouless, 
J. Phys. C. 6, 1734 (1973).

Jacquod Ph. Jacquod and D. L. Shepelyansky, Phys. Rev. Lett. 79, 1837 (1997).

scardicchio_bethe_mbl A. De Luca and A. Scardicchio, Europhys. Lett. 101, 37003 (2013).

dinamica G. Biroli and M. Tarzia, Phys. Rev. B 96, 201114(R) (2017).

logan1 D. E. Logan and S. Welsh, Phys. Rev. B 99, 045131 (2019). 

SUSY_Bethe A.  D.  Mirlin  and  Y.  V.  Fyodorov,  Nucl.  Phys.  B 366, 507 (1991); A. D. Mirlin and Y. V. Fyodorov, Phys. Rev. B 56, 13393 (1997).

large_d  E. Tarquini, G. Biroli, and M. Tarzia, Phys. Rev. B 95, 094204(2017).

lemarie I. García-Mata, J. Martin, R. Dubertrand, O. Giraud, B. Georgeot, and G. Lemarié,
Phys. Rev. Research 2, 012020 (2020).

RRGmulti Note however that the Anderson model on the RRG is known to display a very peculiar form of strong multifractaliy in the whole localized phase <cit.> for the wavefunction moments Υ_q with q<q_⋆ = 0.5, while Υ_q = 0 for q>q_⋆. This corresponds to the multiftactal behavior of very the small amplitudes of the exponentially localized eigenstates. 

mace N. Macé, F. Alet, and N. Laflorencie, Phys. Rev. Lett. 123, 180601 (2019).

laflorencie F. Pietracaprina and N Laflorencie, arXiv:1906.05709

evers F. Evers and A. D. Mirlin, Rev. Mod. Phys. 80, 1355 (2008).

resonances I. V. Gornyi, A. D. Mirlin, D. G. Polyakov, and A. L. Burin, Annalen der Physik 529, 1600360 (2017);
K. S. Tikhonov and A. D. Mirlin, Phys. Rev. B 97, 214205 (2018).



entanglement_D2 This behavior is closely related to the volume-law scaling of the long-time saturation value of the entanglement entropy for the case when the initial state is a basis state <cit.>.

pietracaprina F. Pietracaprina, V. Ros, and A. Scardicchio, Phys. Rev. B 93, 054201 (2016).

kravtsov V. E. Kravtsov, I. M . Khaymovich, E. Cuevas, M. Amini,
New Journal of Physics 17, 122002 (2015).

bogomolny E. Bogomolny and M. Sieber, Phys. Rev. E 98, 042116 (2018).

nosov P. A. Nosov, I. M. Khaymovich, and V. E. Kravtsov, Phys. Rev. B 99, 104203 (2019).

kravtsov1 V. E. Kravtsov, I. M. Khaymovich, B. L. Altshuler, L. B. Ioffe, arXiv:2002.02979;  I. M. Khaymovich, V. E. Kravtsov, B. L. Altshuler, and L. B. Ioffe, arXiv:2006.04827

caveat Note that in fact the accessible volume  at a given extensive energy E_0 is VΩ(E_0), where Ω(E_0) = e^S(E_0) is the many-body density of states, given by the exponential of the microcanonical entropy, S(E_0).

facoetti D. Facoetti, P. Vivo, and G. Biroli, EPL 115, 47003 (2016).

chalker S. Roy, D. E. Logan, and J. T. Chalker, Phys. Rev. B
 bf 99, 220201(R) (2019); S. Roy, J. T. Chalker, and D. E. Logan, Phys. Rev. B
99, 104206 (2019).

znidaric M. Žnidarič, T. Prosen, and P. Prelovšek, Phys. Rev. B
77, 064426 (2008).

serbyn M. Serbyn, Z. Papić, and D. A. Abanin,
Phys. Rev. X 5, 041047 (2015).

devakul T. Devakul  and R.  R.  P.  Singh,  Phys.  Rev.  Lett. 115,  187201 (2015).

doggen_sub E. V. H. Doggen, F. Schindler, K. S. Tikhonov, A. D. Mirlin, T. Neupert, D. G. Polyakov, I. V. Gornyi, Phys. Rev. B 98, 174202 (2018).

abanin  D. A. Abanin, J. H. Bardarson, G. De Tomasi, S. Gopalakrishnan, V. Khemani, S. A. Parameswaran, F. Pollmann, A. C. Potter, M. Serbyn, and R. Vasseur, arXiv:1911.04501

huseQP S. Iyer, V. Oganesyan, G. Refael, and D. A. Huse, Phys. Rev. B, 87 (2013) 134202.

roscilde T. Roscilde,  P. Naldesi, and E. Ercolessi, SciPost Phys. 1 010 (2016).

barlev Y. Bar Lev, D. M. Kennes, C. Klöckner, D. R. Reichman, and C. Karrasch, EPL 119, 37003 (2017).

logan S. Roy and D. E. Logan, arXiv:1911.12370

baldwin C. L. Baldwin, C. R. Laumann, A. Pal,  and A. Scardicchio, Phys. Rev. B 93, 024202 (2016); C. L. Baldwin, C. R. Laumann, A. Pal,  and A. Scardicchio, Phys. Rev. Lett.118, 127201 (2017); C.  L.  Baldwin  and  C.  R.  Laumann,  Phys.  Rev.  B 97, 224201 (2018).

colmenarez L. Colmenarez, P.A. McClarty, M. Haque, and D.J. Luitz,
SciPost Physics 7 (2019).

rosenzweig N. Rosenzweig and C. E. Porter, Phys. Rev. 120, 1698 (1960).

disclaimer Note that here we are implicitly generalizing the results obtained at q=0 (i.e., Hamming distance n/4 on the HS for the XXZ model) to all other x-sector, assuming that the same effective exponent γ describes the scaling of the matrix elements at all distances, while one should instead compute the exponent γ in each sector separately, and verify that ergodicity is broken for all x.

disclaimer1 Note that we are implicitly assuming that the loops of the Bethe lattice (if any) are much larger than n. This means that the total number of sites must be such that N ≫ k^n.

tikhonov_critical K. S. Tikhonov and A. D. Mirlin, Physical Review B 99, 214202 (2019).

mirlinCT K. S. Tikhonov and A. D. Mirlin,
Phys. Rev. B 94, 184203 (2016); M. Sonner, K. S. Tikhonov, A. D. Mirlin,
Phys. Rev. B 96, 214204 (2017).

Ioffe B. L. Altshuler, E. Cuevas, L. B. Ioffe, V. E. Kravtsov, Phys. Rev. Lett. 117, 156601 (2016); B. L. Altshuler, L. B. Ioffe, V. E. Kravtsov, arXiv:1610.00758;
V. E. Kravtsov, B. L. Altshuler, L. B. Ioffe, Annals of Physics 389, 148 (2018).

MG C. Monthus and T. Garel, J. Phys. A: Math. Theor. 44, 145001 (2011).

DPRM G. Biroli and M. Tarzia, arXiv:2003.09629

gabriel I. Garcia-Mata, O. Giraud, B. Georgeot, J. Martin, R. Dubertrand, G. Lemarié,
        Phys. Rev. Lett. 118, 166801 (2017).

Bethe G. Biroli and M. Tarzia, arXiv:1810.07545

aizenman M. Aizenman and S. Warzel, 
Europhys. Lett. 96, 37004 (2011);
M. Aizenman and S. Warzel, 
Phys. Rev. Lett. 106, 136804 (2011).

QREM The only exception to this is provided by the quantum version of the Derrida's Random Energy Model, studied in L. Faoro, M. V. Feigel'man,  and L. Ioffe, Annals of Physics 409, 167916 (2019) and in Refs. <cit.>.

yukalov D. Cohen, V. I. Yukalov, and K. Ziegler, Phys. Rev. A 93, 042101 (2016).

BarLevnonergo D. J. Luitz, I. M. Khaymovich, and Y. Bar Lev, arXiv:1909.06380

pino M. Pino, J. Tabanera, P. Serna,
J.Phys. A: Math. and Theor. 52, 475101 (2019).

altman  I.-D. Potirniche, S. Banerjee, and E. Altman,  Phys. Rev. B 99, 205149 (2019).

prosen J. Šuntajs, J. Bonča, T. Prosen, and L. Vidmar, arXiv:1905.06345 


