Out of equilibrium Phase Diagram of the Quantum Random Energy Model
Giulio Biroli, Davide Facoetti, Marco Schiró, Marco Tarzia, Pierpaolo Vivo
OOut of equilibrium Phase Diagram of the Quantum Random Energy Model
Giulio Biroli, Davide Facoetti,
1, 2
Marco Schir´o, Marco Tarzia,
4, 5 and Pierpaolo Vivo Laboratoire de Physique de l’Ecole normale sup´erieure, ENS, Universit´e PSL,CNRS, Sorbonne Universit´e, Universit´e de Paris, F-75005 Paris, France Department of Mathematics, King’s College London, Strand, London, WC2R 2LS, United Kingdom JEIP, USR 3573 CNRS, Coll`ege de France, PSL Research University, F-75321 Paris, France LPTMC, CNRS-UMR 7600, Sorbonne Universit´e, 4 Pl. Jussieu, F-75005 Paris, France Institut Universitaire de France, 1 rue Descartes, 75231 Paris Cedex 05, France
In this paper we study the out-of-equilibrium phase diagram of the quantum version of Derrida’sRandom Energy Model, which is the simplest model of mean-field spin glasses. We interpret itscorresponding quantum dynamics in Fock space as a one-particle problem in very high dimensionto which we apply different theoretical methods tailored for high-dimensional lattices: the Forward-Scattering Approximation, a mapping to the Rosenzweig-Porter model, and the cavity method.Our results indicate the existence of two transition lines and three distinct dynamical phases: acompletely many-body localized phase at low energy, a fully ergodic phase at high energy, anda multifractal “bad metal” phase at intermediate energy. In the latter, eigenfunctions occupy adiverging volume, yet an exponentially vanishing fraction of the total Hilbert space. We discuss thelimitations of our approximations and the relationship with previous studies.
I. INTRODUCTION
As discovered over 10 years ago by the seminal work ofBasko, Aleiner, and Altshuler, isolated disordered inter-acting many-body systems can show absence of transportand thermalization even at finite energy density if the dis-order is strong enough. This is known as Many-Body Lo-calization (MBL) and is a purely quantum phenomenonwhich occurs due to Anderson localization in the Fockspace as the result of the interplay of disorder, quantumfluctuations, and interactions, and gives rise to a com-pletely new mechanism for ergodicity breaking, that pro-duces a robust dynamical phase of matter which is stablewithin a range of interaction and other Hamiltonian pa-rameters. This remarkable phenomenon has attractedconsiderable interest recently—see Refs. [4–8] for recentreviews—as it implies that the long-time properties ofMBL systems cannot be described by the conventionalensembles of quantum statistical mechanics: They canremember, forever and locally, information about theirinitial conditions.Although significant and exciting progress has beenmade in understanding these phenomena in recent years,both in theory and experiment, there still remainmany open issues. A first set of open questions is aboutthe nature (the universality class) of the MBL phase tran-sition between the thermal and localized phases as therandomness is increased. This transition is an eigen-state phase transition, marked by a sharp change inproperties of the many-body wave-functions and thus inthe dynamics of the system. However, the behavior ofmany-body eigenstates in the Hilbert space is not firmlyestablished: For instance, it is still debated whetherthere is only one phase transition, or could there pos-sibly be some sort of intermediate phase that is nei-ther fully localized nor fully thermal, where eigenstatesare delocalized but non-ergodic, called the “badmetal” phase. Investigations of the MBL transition so far have mostly been numerical studies based on ExactDiagonalization (ED) of relatively small one-dimensionalsystems, and how to do a proper finite-size scal-ing analysis of these numerical data remains unclear.
Also in experiments it is challenging to access the verylong time, and possibly also long length scales, on whichthe critical behavior develops. Another frontier of thefield is directed towards understanding the existenceof MBL in higher dimensions, and its relationshipwith other form of ergodicity breaking such as quantumglassiness.
In this context, exactly solvable, mean-field-like, andsimplified toy models might naturally play an impor-tant role in making some progress, at least partially,in these directions and in improving our understand-ing of MBL. Also developing new techniques and toolsto tackle analytically or semi-analytically the transitionand its properties might provide an important step for-ward to shed new lights on some of the problems men-tioned above. With this in mind, in this paper we studythe out-of-equilibrium phase diagram of the quantumversion of Derrida’s Random Energy Model, which isthe simplest toy model of mean-field spin glasses. Thequantum model’s equilibrium phase diagram has beenstudied before and a glassy phase was identified atlow temperature and small transverse magnetic field.The MBL transition of the QREM was also investigatedpreviously. In Refs. [25] and [26] the presenceof a mobility edge separating ergodic eigenstates frommany-body localized ones was established, based on EDsof small systems and on a perturbative expansion built onthe Forward-Scattering Approximation (FSA). A laterstudy identified three distinct dynamical phases, referredto as trapped, tunneling and excited phases in the con-text of quantum optimization problems. The interpre-tation of those phases, and in particular of the phaseright above the MBL localized phase, has been recentlyput into question. In particular, in Ref. [28] the MBL a r X i v : . [ c ond - m a t . d i s - nn ] S e p phase has been identified with an hyperglass where dy-namics is practically absent while the entire phase above T MBL with a dynamical glass phase (or bad metal) char-acterized by non-ergodic extended (NEE) eigenstates. InRef. [34] the authors derive an estimate of the transi-tion to NEE eigenstates in agreement with Ref. [28], andargue that the NEE phase is layered in an alternatingsequence of two distinct subphases. The dynamical pop-ulation transfer protocol on the QREM was further ana-lyzed in Ref. [33], yielding a numerical estimation of thedynamical phase diagram and of the fractal dimensionsof the eigenstates in the NEE regime.In this work we revisit the problem of the out of equi-librium phase diagram of the QREM using two comple-mentary techniques, the first based on the FSA and ona mapping onto a paradigmatic random matrix model,the Rosenzweig-Porter (RP) model, and the secondbased on a generalization of the self-consistent theory oflocalization (hereafter called the “cavity approach”)designed to take into account the local structure of theHilbert space of the QREM. In agreement with Ref. [28]we find a multifractal “bad metal” phase in a broadrange of intermediate energies, where eigenfunctions aredelocalized but non-ergodic, and out-of-equilibrium re-laxation to thermal equilibrium is expected to be veryslow (exponential in the system size). We also ob-tain a second transition into a fully delocalized ergodicphase at higher temperatures.The paper is organized as following: In section II weintroduce the QREM model and recall basic propertiesof its equilibrium phase diagram. Section III describesthe mapping to a single particle tight-binding Andersonproblem in Hilbert space. Section IV provides qualitativearguments for the phase diagram within the FSA. Sec-tion V contains the cavity approach and the numericalresults found with this method. In Section VI we summa-rize the results found with our approximations and dis-cuss their relationship with previous results. Section VIIput forward an interpretation of the results based on afamily of auxiliary Anderson “toy” models. Finally, inSec. VIII we provide concluding remarks and perspectivefor future studies. Several technical details are reportedin the Appendices A-C.
II. THE MODEL
The Quantum Random Energy Model (QREM) for N spin-1 /
2s is defined by the following Hamiltonian: H QREM = E ( { ˆ σ za } ) − Γ (cid:88) a ˆ σ xa , (1)where Γ is the transverse field, and E ( { ˆ σ za } ) is a ran-dom operator diagonal in the { ˆ σ za } basis, which takes 2 N different values for the 2 N configurations of the N spinsin the z -basis, identically and independently distributed according to: P ( E ) = e − E /N √ πN . (2)Such natural choice of the scaling of the random many-body energies insures that they are with high probabil-ity contained in the interval [ − N √ log 2 , + N √ log 2] inthe thermodynamic limit. Throughout the paper we willdenote by ε = E/N the intensive energy per spin cor-responding to the extensive energy E . A concrete im-plementation of the E ( { ˆ σ za } ) is given by the p → ∞ limit of the fully-connected p -spin model: E ( { ˆ σ za } ) =lim p →∞ (cid:80) a ,...,a p J a ,...,a p ˆ σ za · · · ˆ σ za p , where the J a ,...,a p are i.i.d. Gaussian random variables. In consequence,the diagonal part of the Hamiltonian corresponds to amean-field spin-glass model, and exhibits all the fun-damental features of the so-called “random first-ordertheory”, with a 1-step replica symmetry breaking glasstransition. The equilibrium properties (in the canonical ensem-ble) of the QREM are well established.
At low trans-verse field, it displays the same transition as the classicalmodel between the paramagnetic and the glass phase atthe Kauzmann temperature T K = 1 / (2 √ log 2). All ther-modynamic quantities are identical to the Γ = 0 classicalREM. For T < T K the system freezes in its ground stateat ε GS = −√ log 2. The cassical model has also a dy-namical ergodicity-breaking transition below which thetime to reach thermal equilibrium is exponentially largein the system size . However, differently from the p -spinmodels with finite p , for which the dynamical transitiontemperature T d is finite, in the p → ∞ limit T d → ∞ , dueto the fact that the random energies and spin configura-tions are totally uncorrelated and flipping a single spincan change the energy by a large amount. Within thesemiclassical approximation of Ref. [30] T d stays infiniteeven when quantum fluctuations are turned on (Γ > > Γ c ( T ) the system is astandard quantum paramagnet, and the REM term in theHamiltonian does not influence the equilibrium physics ofthis phase. The first-order transition between these tworegions takes place at Γ c ( T ), which is equal to √ log 2 ≈ .
833 for T → √ / ≈ .
707 for T → ∞ . Inthis paper we will only focus on the small-Γ region of thephase diagram (i.e., Γ < √ / III. MAPPING TO ANDERSONLOCALIZATION ON THE HYPERCUBE
The QREM defined by Eq. (1) can be viewed as thesimplest many-body model that displays Anderson local-ization in its Hilbert’s space: If one chooses as a basisthe tensor product of the simultaneous eigenstates of theoperators σ za , the Hilbert space of the many-body Hamil-tonian is a N -dimensional hypercube of V = 2 N sites.One can map a configuration of N spins to a corner ofthe N -dimensional hypercube by considering σ za = ± a -th dimension. Therandom part of the Hamiltonian is by definition diagonalon this basis, and gives uncorrelated random energies oneach site orbital of the hypercube: At Γ = 0 the many-body eigenstates of Eq. (1) are simply product states ofthe form | σ z (cid:105) ⊗ | σ z (cid:105) ⊗ · · · ⊗ | σ zN (cid:105) , and the system is fullylocalized. The interacting part of the Hamiltonian actsas single spin flips on the configurations { σ za } , and playsthe role the hopping rates connecting “neighboring” sitesin the configuration space. The many-body quantum dy-namics is then recast as a single-particle non-interactingtight-binding Anderson model for spinless electrons in adisordered potential living on the 2 N corners of an hy-percube in N dimensions, with the spin configurationsbeing “lattice sites”, | i (cid:105) ≡ |{ σ za }(cid:105) , and the transversefield playing the role of the hopping amplitude betweenneighboring sites: H = − Γ (cid:88) (cid:104) i,j (cid:105) ( | i (cid:105)(cid:104) j | + | j (cid:105)(cid:104) j | ) + V (cid:88) i =1 E i | i (cid:105)(cid:104) i | , (3)where (cid:104) i, j (cid:105) denotes nearest-neighbors on the hypercube, V = 2 N is the total number of sites, and Γ is the hop-ping kinetic energy scale. This mapping is exact , in thesense that the Hamiltonians (1) and (3) have the same eigenvalues and the eigenvectors (when the simultaneouseigenstates of the operators σ zi is chosen as a basis). How-ever, for a generic interacting many-body Hamiltonian infinite dimensions the random energies defined on neigh-boring corners of the hypercube are strongly correlated,as they correspond to many-body configurations whichonly differ by a single spin flip, while for the QREM E i are i.i.d. random variables distributed according toEq. (2), since its distinguishing feature is precisely theabsence of such correlations. IV. ESTIMATE OF THE OUT OFEQUILIBRIUM PHASE DIAGRAM WITHIN THEFORWARD-SCATTERING APPROXIMATION
As discussed in Section III, the QREM can be mappedto an Anderson model on a hypercube with V = 2 N sites,labeled by σ z configurations. The typical tunneling ratebetween two configurations depends on their energy andon the Hamming distance x between them ( N x is theminimum number of spin flips which separate the twoconfigurations). Since the energy levels are independent,the typical number of configurations of energy | E | = N ε and at distance x from a given configuration is N ε ( x ) = (cid:18) NN x (cid:19) e − Nε √ πN . (4)Here we estimate the matrix elements M ( ε, x ) betweenthese two configurations by perturbation theory in Γ, us-ing the FSA. This consists in assuming that the ma-trix element between two configurations at distance x is given by the product of the matrix elements obtainedalong the N x spin flips that connect the two configura-tions, ignoring “loopy” contributions in which spins areflipped twice since they contribute at higher order in per-turbation theory. Almost all states have energy O ( √ N ),while E ≈ O ( N ), therefore we take the energy differencesappearing in the denominators appearing in the pertur-bation theory to be E . Since there are ( N x )! such con-tributions, corresponding to the (
N x )! to connect twoconfigurations
N x spin flips away, the resulting matrixelement reads M ( ε, x ) ≈ (cid:18) Γ N ε (cid:19) Nx ( N x )! . (5)Based on this analogy with the Anderson localizationtransition one can estimate the point at which Many-Body delocalization on states with the same intensive en-ergy takes place using the criterion N ε ( x ) |M ( ε, x ) | → Analogously, the ergodicity breaking transition occurswhen the transition rate obtained from the Fermi GoldenRule N ε ( x ) |M ( ε, x ) | → This means that al-though a given state is in resonances with many otherstates of energy ε and at distance x from it, the num-ber of those resonances is not enough for the quantumdynamics to decorrelate in a finite time from the initialcondition.This can be illustrated using an effective Rosenzweig-Porter random matrix, with off-diagonal matrixelements roughly scaling as M ( ε, x ) and connecting anumber of levels which equals the number of states atfixed energy density ε and Hamming distance x . Thisleads to a N ε ( x ) × N ε ( x ) RP random matrix. For thesake of clarity, here we briefly recall the definition of theRP model, whose Hamiltonian is a matrix of size V × V given by the sum of two terms, H RP = E + µ V γ/ G , (6)where E ij = E i δ ij is diagonal with i.i.d. entries (the dis-tribution does not matter as long as it has finite vari-ance), µ is a constant of order one (whose value is unim-portant), and G is a GUE (or GOE) matrix with unitvariance. The latter mimics an ergodic systems (e.g. theclean lattice), while E represents the on-site disorder.Within the analogy sketched above, the effective γ forthe RP model associated to a given set of configurationswith ( ε, x ) is γ = − M ( ε, x ) / log N ε ( x ).The parameter γ acts in the RP model as a proxy of thedisorder strength: at large γ >
2, the GUE contribu-tion is suppressed, and the systems is localized; at small γ <
1, the system is ergodic, while the regime 1 < γ < V − γ sites close in energy. Thecriteria for these two transitions are exactly the ones weintroduced and used before for the QREM.Going back to the QREM, we therefore expect thatthe transition from Anderson localization to nonergodicextended states for the QREM then happens when atleast one x -sector becomes delocalized, i.e.max x N ( ε, x ) M ( ε, x ) ≈ max x e Nf ( x,ε,γ ) (cid:38) , (7)where f ( x, ε, Γ) = x log (cid:18) Γ eε (cid:19) − (1 − x ) log(1 − x ) − ε . (8)If Γ < ε , f is always negative. Otherwise, it has a non-negative maximum at x ∗ = 1 − ε/ Γ, which determinesthe mobility edge Γ
MBL ( ε ) through the implicit equation f ( x ∗ , ε, Γ MBL ) = ε Γ MBL − ε + log Γ MBL e ε = 0 . (9)This is the same result obtained through a similar argu-ment in Ref. [25].Similarly, full ergodicity is recovered by requiring thatat least one x -sector becomes ergodic,max x N ( ε, x ) |M ( ε, x ) | ≈ max x e Nf ( x,ε,γ ) (cid:38) , (10)with f ( x, ε, Γ) = x log x − (1 − x ) log(1 − x ) − ε +2 x log (cid:18) Γ eε (cid:19) . (11)If Γ < ε , f is always negative. Otherwise, it has anon-negative maximum at x ∗ = 1 / (cid:112) − ε/ Γ) ),which gives a different implicit equation for the ergodictransition Γ ETH ( ε ).Expanding the solutions to the two implicit equa-tions around ε = 0, we find Γ MBL ≈ ε , , andΓ ETH ≈ ε/
2. The analogy with the RP model thereforeindicates that the QREM also undergoes two separatelocalization and ergodicity transitions. The estimates forthe transition lines obtained in this way are shown inFigs. 1 and 2.Qualitatively, these results are in agreement with therecent analysis of Refs. [28] and [34], which are based ona different (and more thorough) strategy to estimate theoff-diagonal tunneling rates between different spin config-urations in the Hilbert space, and are also in agreementwith the numerical estimations of Ref. [33]. However,while Eq. (11) predicts that the transition line from thefully ergodic to the dynamical glass (bad metal) regimeis at finite energy density ε ETH , within the approach ofRefs. [28] and [34] the transition line is squeezed to zeroenergy density, i.e., infinite temperature. We will comeback to this point in Sec. VI.
V. CAVITY APPROACHA. Cluster approximation on the hypercube
In large spatial dimensions the neighbors of a givensite are organized in a hierarchical way (i.e., the fraction 0 . . . . − . . . ε Figure 1. Localization (blue) and ergodicity (orange) tran-sition lines, obtained from the mapping to the RP model,Eqs. (8,11). The limits of the y axis coincide with the edgesof the many-body spectrum ( | ε | < (cid:112) log(2)). The red dashedlines correspond to ε = ± Γ (see Refs. [25] and [26]). . . . . . T Figure 2. Ergodicity and localization transition lines (Fig. 1),transposed on the canonical phase diagram, T = 1 / ε . of short loops is suppressed) and their number grows ex-ponentially with the distance. These are distinctive fea-tures of tree-like structures. In fact, it was argued origi-nally in [2] that the (non-interacting) Anderson modelon the Bethe lattice, first introduced and studied inRef. [38], can be thought as a toy model for MBL (seealso Refs. [3, 47–50] for a similar analysis and Ref. [51]for a quantitative investigation of these ideas).Since on tree-like structures the model (3) allows, inprinciple, for an exact solution, which yield the diago-nal elements of the resolvent matrix, assuming that forlarge enough N the hypercube is well approximated bya Bethe lattice provides a very simple way to investigateanalytically, although approximately, the spectral prop-erties of the eigenvectors of the QREM. The simplest approximation consists in takingRandom-Regular Graphs (RRGs) of V = 2 N sites andfixed connectivity N as the underlying lattices mimick-ing the Hilbert space, i.e., random lattices which havelocally a tree-like structure but have loops whose typicallength scales as ln V ∝ N and no boundary, and which arestatistically translationally invariant. In practice, thiscorresponds to shuffling the position of the sites and/orrewiring the connections of the hypercube in a randomway, keeping the total number of sites and the local con-nectivity of the lattice fixed.There is however a potentially important issue relatedto the Bethe approximation, which we explain below.The spectrum of the kinetic term of (3) is given by theDensity of States (DoS) of the adjacency matrix of the N -dimensional hypercube, which coincides with the dis-tribution of the eigenvalues of the second term of theHamiltonian (1), i.e., a simple paramagnet, with ener-gies contained in the interval E ∈ [ − N Γ , N Γ]: ρ HC I ( E ) = Ω( E )Γ2 N +1 , (12)where the term Γ2 N − in the denominator is a normaliza-tion factor that insures that (cid:82) + N Γ − N Γ ρ HC I ( E ) d E = 1, andΩ( E ) is the number of spin configurations at energy E :Ω( E ) = (cid:18) N ( N + E/ Γ) / (cid:19) ∼ Ω e Ns ( ε/ Γ) . (13)Ω is a normalization factor and s ( ε/ Γ) is the entropyper spin at large N (apart from logarithmic corrections)for a polarization m = ε/ Γ = (cid:104) σ x (cid:105) of the spins in the x direction: s ( m ) = log(2) − m m ) − − m − m ) . (14)Approximating the hypercube as a tree-like latticeamounts to replacing Eq. (12) with the spectrum of theadjacency matrix of a Bethe lattice of connectivity N which, for large enough N tends asymptotically to a semi-circle law: ρ RRG I ( E ) ≈ √ N − E π Γ N , (15)with support in the interval E ∈ [ − √ N , √ N ]. Forenergies of order O ( √ N ), where the vast majority ofthe eigenvalues is, ρ RRG I provides in fact a reasonablygood estimation for the true DoS, ρ HC I (see fig. 10 ofAppendix A for a quantitative comparison). Yet, for ex-tensive energies of O ( N ) Eq. (15) completely neglects theexponentially small fraction of eigenvalues in the tails ofthe DoS, corresponding to strongly polarized spins in the x direction. Since at finite energy density ε and N suffi-ciently large, the energy N ε will fall outside the edges ofthe semicircle (15) which describes the spectrum of thedelocalizing kinetic term within the Bethe approxima-tion, the Anderson localization of the eigenfunctions ofthe Hamiltonian (3) will occur in the far Lifshits tails ofthe DoS, and this might affect its properties. In otherwords, the system might appear as localized within theBethe approximation due to the fact that some of the matrix elements associated to the kinetic term are arti-ficially suppressed by approximating the hypercube as aRRG.In order to overcome, at least partially, these limita-tions in this paper we put forward a cluster expansion,which takes into account, at least locally, the specificstructure of the hypercube up to a certain distance (inparticular, including all the shortest loops of length four,six, eight, etc.) and improves systematically the sim-plest, single-site, Bethe approximation. In practice, weconsider clusters of s = 2 n neighboring corners on thehypercube (corresponding to spin configurations whichdiffer by few spin flips only), and obtain self-consistentrecursion equations for the s × s elements of the localresolvent matrix on each cluster by assuming that theclusters are on a tree-like structure.The standard single-site Bethe approximation corre-sponds to n = 0, while the cases n = 2 is schematicallyrepresented in fig. 3 ( n = N corresponds, of course, tothe exact solution of the problem). We will take n = 2throughout.For n = 2 a plaquette of four corners corresponds tofour spin configurations such as, e.g., | (cid:105) = | ↑ , ↑ , σ z , . . . , σ zN (cid:105)| (cid:105) = | ↑ , ↓ , σ z , . . . , σ zN (cid:105)| (cid:105) = | ↓ , ↓ , σ z , . . . , σ zN (cid:105)| (cid:105) = | ↓ , ↑ , σ z , . . . , σ zN (cid:105) , for any configuration { σ za } a =3 ,...,N . The N − N − σ z , . . . , σ zN one by one. Two neighbor-ing plaquettes of four sites are connected by four edges.The Hilbert space will be then approximated as a RRGof 2 N / N −
2. Moredetails are given in Appendix A.One can show that within the cluster approximationthe support of the spectrum of the kinetic term of theHamiltonian (3) becomes indeed broader and broaderas n is increased, E max ≈ Γ(2 √ N − n + n ) (see Ap-pendix A).One can easily obtain, at least formally, the self-consistent recursion relations for the elements of the re-solvent matrix (or Green’s functions) of the Hamilto-nian (3), defined as G ( z ) = ( H − z I ) − , at any order n ofthe cluster expansion. The key objects are the so-called cavity resolvent matrices, G p → q ( z ) = ( H p ↔ q − z I ) − , i.e.,the resolvent matrices of modified Hamiltonians H p ↔ q where all the 2 n edges between the sites of the cluster p and the sites of the cluster q have been removed (graydashed lines of fig. 3). Let us assume, as explained aboveand as sketched in fig. 3, that at large enough N theclusters occupy the vertices of a tree-like structure (atleast locally), and let us imagine to take a given cluster p and its neighbors { q , . . . , q N − n } . If one removes suchcluster from the graph, then the clusters { q , . . . , q N − n } are (quasi-)uncorrelated, since the lattice would break in N − n semi-infinite (quasi-)disconnected branches (ne- q q q N- . . . pq N- Figure 3. Schematic representation of the recursion stepwhich yields the self-consistent equations for the s × s ele-ments of the (cavity) resolvent matrix for clusters of 4 sites.All loops up to length 12 of the hypercube are treated exactly. glecting the large loops of length of order N ). One thenobtains (e.g., by Gaussian integration) the following it-eration relations for the elements of the cavity resolventmatrix on the s sites of the cluster: (cid:2) G − p → q k ( z ) (cid:3) uv = H uv − zδ uv − Γ (cid:88) q l ∈ ∂p/q k [ G q l → p ( z )] uv , (16)where z = N ε + iη , η is an infinitesimal imaginary regula-tor which regularize the pole-like singularities in the righthand sides (see below), ε is the intensive energy densityaround which one chooses to study the spectral proper-ties, H uv are the matrix elements of the Hamiltonian (3)between the sites u and v belonging to the cluster p , and ∂p/q denotes the set of all N − n neighbors of the cluster p except q . The indices u, v = 1 , . . . , s identifying thesites belonging to each clusters are chosen as in fig. 3.(Note that for each cluster with N − n neighbors one candefine N − n cavity Green’s functions and N − n recursionrelations of this kind.)After finding the solution of Eqs. (16), one can finallyobtain the resolvent matrix of the original problem on agiven cluster p as a function of the cavity Green’s func-tions on the neighboring clusters: (cid:2) G − p ( z ) (cid:3) uv = H uv − zδ uv − Γ (cid:88) q k ∈ ∂p [ G q k → p ( z )] uv . (17)For n = 0 these equations simply give back the standardrecursion relations for the Anderson model on the Bethelattice (with connectivity N and Gaussian iid randomenergies). For n = 1 and n = 2 simple analytic ex-pressions of the inverse of the local s × s resolvent matri-ces are known, which allows one to write simple recursionequations for its s diagonal elements and its s ( s − / n ≥
3, however, the local inversion involved in Eqs. (16)and (17) must be done numerically at each iteration step.There are essentially two ways, that we detail in Ap-pendix B, to solve the recursion equations for the Green’sfunction and obtain information on the spectral statis-tics at finite N . The most accurate strategy, to whichwe will refer to as “Cluster Belief Propagation” (C-BP)algorithm (see Ref. [54] for a detailed explanation of thisapproach for the usual tight-binding Anderson model onthe Bethe lattice), is to solve directly Eqs. (16) and (17)on random realizations of the hypercube of 2 N sites (i.e., N spins) (see Appendix B for details). However, thanksto the fact that random energies E i of the QREM are un-correlated, in order to access larger system sizes one canadopt another strategy, to which we will refer hereafteras “Cluster Population Dynamics” (C-PD) algorithm, which consists in interpreting the recursion relations forthe Green’s functions as equations for their probabilitydistributions once the average over the disorder is taken.In fact, since G p → q and G p are random matrices, onecan assume that averaging over the on-site random en-ergies leads to functional equations on their probabilitydistribution Q ( G ) and Q ( G ) (see Appendix B for details).Fig. 11 of App. B shows that the cluster approximationprovides a quite accurate approximation of local observ-ables, such as the distribution of the LDoS, as comparedto exact diagonalization for small systems. B. Spectral statistics and the η → + limit The statistics of the diagonal elements of the resolventgives—in the η → + limit, see below—the spectral prop-erties of H . In particular, the probability distribution ofthe Local Density of States (LDoS) at energy E = N ε isgiven by: ρ i ( ε ) = (cid:88) α |(cid:104) i | α (cid:105)| δ ( N ε − E α ) = lim η → + Im G i ( N ε + iη ) π , (18)from which the average Density of States (DoS)is simply obtained as ρ ( ε ) = (1 / V ) (cid:80) i ρ i ( ε ) =lim η → + / ( V π )Tr Im G ( N ε + iη ). (We have defined G i =[ G p ] uu with i = 2 n p + u , p = 1 , . . . , N − n and u =1 , . . . , s = 2 n ). Similarly, the spectral representation ofthe Inverse Participation Ratio of the eigenstates | α (cid:105) ofenergy E α close to N ε can be obtained as:Υ ( ε ) = (cid:88) i |(cid:104) i | α (cid:105)| δ ( N ε − E α )= lim η → + η (cid:80) V i =1 |G i ( N ε + iη ) | (cid:80) V i =1 Im G i ( N ε + iη ) . (19)Another useful observables is the typical DoS: ρ typ ( ε ) = lim η → + exp (cid:16) V − (cid:80) V i =1 ln Im G i ( N ε + iη ) (cid:17) V − (cid:80) V i =1 Im G i ( N ε + iη ) . (20)However, at this point we encounter another difficultywhich is due to the very unusual (and simultaneous) scal-ing with N of the parameters of the Hamiltonian (3).In fact, the dependence of the random energies and theconnectivity of the lattice on N produces a density ofstates that strongly concentrate around zero energy den-sity, as naturally expected for many-body systems: Atsmall Γ one has that ρ ( E ) ≈ P ( E ) = e − E /N / √ πN ,while at large Γ one expects that ρ ( E ) ≈ Ω( E ) / (Γ2 N +1 ),see Eq. (12). Thus, in both cases the mean level spacing δ ( E ) = 1 / ( V ρ ( E )) is well-defined locally, but dependsstrongly (i.e., exponentially) on the local energy density.In particular, for small Γ one has that δ ( ε ) ≈ √ πN e N ( ε − log(2)) . (21)In order for Eqs. (18)-(20) to be well defined, the limit η → + should be taken in such a way that the imaginaryregulator goes to zero on the same scale as the mean levelspacing. Hence, studying the asymptotic behavior ofthe model at large N implies varying simultaneously thefollowing parameters:- The total number of sites diverges exponentially V = 2 N ;- The connectivity of the lattice grows as N ;- The variance of the random on-site energies growsas (cid:112) N/
2, according to Eq. (2);- The energy at which we study the system grows as
N ε , with ε of O (1);- The imaginary regulator vanishes exponentially as δ ( ε ) = 1 / ( V ρ ( ε )).Thus, the N → ∞ limit of the model (3) is quite differ-ent from the usual thermodynamic limit of the standard(non-interacting) Anderson model (where one just takesthe limit V → ∞ and η → ε , and N of the probabilitydistribution of the LDoS, Eq. (18), from which one cancompute several spectral quantities of interest, such as,e.g., the IPR (19) and the typical value of the LDoS (20).In particular, in order to assess the ergodicity of the wave-functions, it is custom to introduce the fractal dimensionsdefined through the asymptotic behavior of these two lat-ter quantities: Υ ( ε ) ∼ (cid:104) ˜ V ( ε ) (cid:105) − D ( ε ) ,ρ typ ( ε ) ∼ (cid:104) ˜ V ( ε ) (cid:105) D ( ε ) − . (22)This definition takes into account the actual volume ofthe portion of the phase space accessible at finite en-ergy density ε , ˜ V ( ε ) = 2 N ρ ( ε ), since at finite energy den-sity ε , ergodic eigenstates are uniformly spread over the
10 20 30 40 50 N D ε=0.08ε=0.1ε=0.12ε=0.14ε =0.2 ε =0.24 ε =0.26 ε =0.32
10 20 30 40 50 N D Figure 4. D (left) and D (right) versus N for Γ = 0 . n = 2), and shadow symbols correspond to the result ofthe C-PD approach (with n = 2). Two independent ED data-sets, obtained with two different algorithms, are shown, togive the idea of the typical size of the errorbars. The verticaldashed lines represent the limits of the range of applicabilityof ED ( N ≤
14) and the C-BP approach ( N ≤ D and D at large N in the bad metal phase. (Note that thegeneralized fractal dimensions can become larger than 1 forsome intermediate values of N and for ε small enough, due tologarithmic corrections to the many-body DoS.) hyper-surface at constant energy. In the small-Γ part ofthe phase diagram, as a first approximation one has that ρ ( E ) ≈ P ( E ). C. Numerical results for the fractal dimensions
In fig. 4, where we plot D (right panel) and D (leftpanel) as a function of N for Γ = 0 .
2, and several valuesof the energy density ε , mostly on the delocalized sideof the MBL transition and close to the mobility edge( ε (cid:46) Γ). D and D are obtained as numerical deriva-tives of ρ typ and Υ with respect to log ˜ V . The figureshows three data-sets, corresponding to the results ob-tained from EDs ( N ≤
14, filled symbols), the C-BPapproach ( N ≤ n = 2, open symbols), and the C-PDapproximation ( N ≤ n = 2, shadow symbols).First of all, these plots demonstrate that the C-BP andC-PD approximations are in reasonably good agreementwith the exact results for all values of ε , at least in therange of system sizes accessible via EDs (see also Fig. 11of App. B for a detailed comparison of the full probabilitydistribution of the LDoS).At moderate energy density both fractal dimensionsshow a clear non-monotonic dependence: D and D firstrapidly increase with N , and then start to decrease slowlyafter going through a maximum. At small enough energydensity, ε (cid:46) .
14, both D and D reach a finite plateaustrictly smaller than one at large N (horizontal dottedlines). This behavior corresponds to genuine multifrac-tal eigenstates, as recently predicted in Ref. [28], and isfound in a broad range of energy density. The lower theenergy, the higher are the plateau values reached by D and D , i.e., the system gets closer and closer to full er-godicity as the energy density is decreased. At largerenergies, instead, above the mobility edge ε MBL , D and D decay to zero in the large N limit.These results support the existence of two distinct non-ergodic regions of the phase diagram: a delocalized mul-tifractal phase (0 < D , D <
1) at intermediate energydensity, where eigenstates occupy a volume that diverges,yet is exponentially smaller than the total Hilbert space,and a Anderson localized phase ( D , D → N -indepentent volume onthe hypercube. We have repeated the same analysis forΓ = 0 .
1, finding similar results.Note, however, that the cavity approach does not allowone to determine sharply the phase boundaries betweenthe three phases because the numerical results are onlyavailable for systems of moderate size, N (cid:46)
50, and theasymptotic values of the fractal dimensions cannot befirmly established, especially in the vicinity of the tran-sition line between the bad metal and the MBL phases(see, e.g., the data for ε = 0 . . ε MBL ∈ (0 . , . D. Level statistics
A natural question that arises concerns the statisticsof the energy levels in the multifractal phase. In fact,in analogy with the RP model it is reasonable to ex-pect that in the mixed phase the level statistics shouldbe described by the GOE ensemble on the scale of themean level spacing, while it might crossover to a differ-ent, possibly non-universal behavior, on a larger energyscale ( ∝ ˜ V D − ) which goes to zero exponentially with N but stays much larger than δ . This scenario is alsosupported by general arguments based on the conver-gence of the Dyson Brownian Motion to its stationaryGOE distribution.
In order to check this idea wehave analyzed the behavior of the level compressibility χ ( ε ; ω ) for the number of energy levels inside the inter-val [ N ε − ω/ , N ε + ω/ which display different scal-ing behaviors for the ergodic, localized, and multifractalstates. The number of energy levels inside an en-ergy interval of width ω (and centered around N ε ) isdefined as L ( ε ; ω ) = (cid:82) ε + ω/ ε − ω/ (cid:80) N α =1 δ ( E (cid:48) − E α ) d E (cid:48) . Thelevel compressibility is then the ratio between the vari- N N χ / N ε=0.08ε=0.1ε=0.12ε=0.14ε=0.26 Figure 5. Level compressibility (obtained from the C-BPapproach with n = 2) on the scale of the mean level spacingrescaled by the GOE asymptotic behavior, 2 N χ ( ε ; 2 cδ ) /N , asa function of N for Γ = 0 . ε . The blackdashed line corresponds to χ = 1 (Poisson statistics). ance of L ( ε ; ω ) and its average: χ ( ε ; ω ) = ( L ( ε ; ω )) − L ( ε ; ω ) L ( ε ; ω ) , where · · · denotes the average over the disorder. In thediffusive regime of the standard ergodic metallic phase,described by the Wigner-Dyson statistics, energy lev-els strongly repel each other, and the variance scalesas ( L ( ε ; ω )) − L ( ε ; ω ) ∝ ln L ( ε ; ω ). The level com-pressibility thus vanishes as χ ( ε, ω ) ∝ N ln 2 / N for large N . Conversely, in the localized phase energy levels arethrown as random points on a line and are described bya Poisson distribution. Hence ( L ( ε ; ω )) − L ( ε ; ω ) = L ( ε ; ω ) and χ ( ε ; E ) → N → ∞ . Finally, for non-ergodic multifractal states the variance of the numberof energy levels inside an interval should scale linearlywith the average, at least in simplest scenarios, and χ ( ε ; ω ) is expected to converge to a (system-dependent)constant between 0 and 1 in the large N limit. In thefollowing for simplicity we will only focus on the behav-ior of the level compressibility when the energy interval ω is taken of the order of the mean level spacings. Inparticular we will set ω = 2 η = 2 cδ , where η is given inEq. (B1) and c = 64. As shown in Refs. [54] and [58] a simple spectral rep-resentation of L ( ε ; ω ) can be achieved in the frameworkof the C-BP approach, in terms of the resolvent matricesdefined on the clusters of the hypercube and of the cav-ity resolvent matrices defined on the edges between theclusters: L ( ε ; ω ) = 1 π lim η → + (cid:26) N − n (cid:88) p =1 (cid:2) Ψ p ( z + ) − Ψ p ( z − ) (cid:3) + (cid:88) (cid:104) p,q (cid:105) (cid:2) ϕ p ↔ q ( z + ) − ϕ p ↔ q ( z − ) (cid:3)(cid:27) , (23)where z ± = N ε ± ω/ iη , the angle Ψ p ( z ) is definedas the phase of det G p ( z ), det G p ( z ) = | det G p ( z ) | e i Ψ p ( z ) ,and the angle ϕ p ↔ q ( z ) is defined as the phase of det ( I s − Γ G q → p ( z ) G p → q ( z )) (we have chosen here to put thebranch-cut in the complex plane along the negative realaxis).In order to analyze the scaling properties of the levelcompressibility we then just need to compute the averageof L ( ε ; ω ) and its fluctuations over many independent re-alizations of the random energies of the hypercube. Thescaling behavior of χ when ω is taken on the scale of themean level spacing, ω = 2 cδ , are shown in fig. 5, wherewe plot the compressibility (divided by the GOE asymp-totic) versus N for Γ = 0 . N χ/N has a non-monotonic behavior roughlyon the same scale as D and D , and seems to approacha finite value at large N which grows as ε is increasedtowards the localized phase. For ε (cid:38) .
26, in the MBLphase, χ tends instead to 1 at large N , as expected forPoisson statistics. This implies that in the whole multi-fractal region, on energy scales proportional to δ the levelcompressibility goes to zero in the thermodynamic limitas in the diffusive regime of the standard metallic phase. VI. DISCUSSION
Above we have presented two complementary approxi-mate strategies to determine the out of equilibrium phasediagram of the QREM. The first approach is based on theFSA and on the mapping to the RP model, while thesecond approach is a generalization of the self-consistenttheory of AL adapted to take into account (at leastpartially) the local structure of the Hilbert space of theQREM. As discussed more in detail in App. B the latterpossibly provides a quite accurate approximation of localobservables, such as the distribution of the LDoS, whilethe former is expected to yield a better estimation of cor-relations, since it is able to capture the fact that thereis a factorial number of paths connecting two points atlarge distance in the configuration space.In this section we discuss and compare more in detailsthe two approaches for what concerns the behavior of thefractal dimensions D and D as a function of the energydensity ε , that we plot in fig. 6 for Γ = 0 . N limit.Within the analogy between the QREM and the RPmodel discussed in Sec. IV, the fractal dimension is ex-pected to be equal to one in the ergodic phase, | ε | < F r a c t a l d i m en s i on s ε D (cavity, Fig. 4)D (cavity, Fig. 4)D=2- γ [Ref. [26]]D [Ref. [34]] Figure 6. Fractal dimensions D (green squares) and D (yellow circles) obtained from the cavity approximation asa function of ε for Γ = 0 .
2, determined by estimating theheight of the plateaus of fig. 4 at large N . The yellow shadedregion indicates the energy interval within which the numer-ical results of fig. 4 suggest that the MBL transition occurs,0 . (cid:46) ε MBL (cid:46) .
19. The vertical dashed blue lines show theposition of ε ETH ≈ . ε MBL ≈ . D obtained in Refs. [28] (orange line) and [34](magenta line), which predict that D → ε → ε ETH → ε ETH , and to zero in the MBL phase, | ε | > ε MBL , and isconjectured to decrease from 1 to 0 in the intermediateNEE regime (the multifractal spectrum should be ob-tained as an “average” of the effective fractal dimensionsover all x -sectors).We also plot the estimations of Refs. [28] (orange line)and [34] where the effective spectral dimension D is ob-tained using similar (although probably more accurate)methods to evaluate the amplitude of the tunneling rates (cid:104){ σ za }| Γ σ xb |{ σ za } (cid:48) (cid:105) between two distant many-body config-urations. The symbols are the cavity-cluster predictions,extracted from the largest size available when reasonablyconverged to a plateau (see fig. 4). The shaded area in-dicates the energy interval within which the numericalresults of fig. 4 suggest that the MBL transition shouldtake place. Due to the limited range of system sizes ac-cessible via the cavity approach, we are not able to con-clude whether the fractal dimension would continuouslygo to zero at ε MBL or rather exhibit a finite jump at thetransition.All approaches agree in indicating the existence ofthree different phases of the QREM: a fully ergodicregime at low energy density, | ε | < ε ETH , a NEE (orbad metal, or dynamical glassy) one at intermediate en-ergy density, ε ETH < | ε | < ε MBL , where the time toreach thermal equilibrium is exponentially large in thesystem size and eigenvectors are extended but multifrac-tal, corresponding to 0 < D , <
1, and a fully localizedone, with Anderson localized eigenstates, at high energy,0 | ε | > ε MBL where D , → ε ETH →
0. The cavity approx-imation indicates that if ε ETH is finite, it is significantlysmaller than the estimation of Eqs. (10) and (11). As weare going to see in the next section, an argument based ona simplified solution of the cavity equations, equivalentto an auxiliary Anderson model in unconventional ther-modynamic limit, also seem to suggest that ε ETH mightindeed squeeze to zero energy in the thermodynamic limit(as (cid:112) log
N/N in the large N limit). The other difference is that within the cavity approachthe fractal dimensions D and D might possibly exhibitan abrupt jump from a finite value smaller then one tozero at the transition between the nonergodic delocal-ized and fully localized regime, while the mapping to theRP model indicates that if one identifies D = 2 − γ , the fractal dimensions should vanish continuously at theMBL mobility edge, as also found in Refs. [28] and [34].These are still open questions for future investigations. VII. AUXILIARY ANDERSON MODELS ANDUNCONVENTIONAL THERMODYNAMICLIMIT
In this section we further simplify the quantum cav-ity analysis and introduce a family of auxiliary “toy”Anderson tight-binding models on a Bethe lattice withconnectivity N (cid:29) V of the system istreated as an independent parameter from N and takenequal to infinity from the start. The basic idea behindthis approach is that in the original Anderson model onthe Hypercube (see section III) the scaling of the numberof sites and number of neighbors with respect to the size N of the original QREM is remarkably different, the for-mer being exponential V = 2 N while the latter is linear k = N . As such one could hope to get some insight onthe solution of the quantum cavity equations by takingthe volume to infinity first, possibly providing an esti-mation for the transition line between the fully ergodicETH phase and the multifractal bad metal one.Concretely, we consider a hybrid version of themodel (3) where the total number of sites of the latticeis sent to infinity from the start keeping N fixed. Hence,for any given choice of Γ and ε , this leads to a fam-ily of tight-binding Anderson model parametrized by theconnectivity N , with random on-site disorder of variance (cid:112) N/
2, given by Eq. (2). The advantage of this procedureis that now for any choice of Γ, ε , and N , the imaginaryregulator can be taken as infinitesimally small (since themean level spacing vanishes in the thermodynamic limit)and we can study whether the system is in the localizedor in the extended phase with standard techniques. We start by determining the mobility edge ε loc ( N ) = E loc ( N ) /N of the auxiliary models by computing theLyapunov exponent which describes the evolution of theimaginary part ∆ of the self-energy Σ, once the iterationrelations (16) have been linearized. At a given order n of the cluster expansion, the (cavity) self-energy on acluster p of s = 2 n sites (in absence of the 2 n edges withone of the neighboring clusters q ) is a s × s matrix definedas: Σ p → q = S p → q + i ∆ p → q = H p − z I s − G − p → q , where H p is the Hamiltonian (3) acting on the sites of thecluster, and I s is the s × s identity matrix. In the localizedphase its imaginary part vanishes for η → + . Hence, inthe thermodynamic limit and close to the localizationtransition, one can take the limit η → + from the startand linearize the recursive equations (16) with respect to∆:[ S p → q k ] uv = Γ (cid:88) q l ∈ ∂p/q k [ H q l − N ε I s − S q l → p ] − uv , [∆ p → q k ] uv = Γ (cid:88) q l ∈ ∂p/q k s (cid:88) w,y =1 [ H q l − N ε I s − S q l → p ] − uw × [∆ q l → p ] wy [ H q l − N ε I s − S q l → p ] − yv . (24)Since S p → q and ∆ p → q are random matrices, these equa-tions naturally lead to functional self-consistent equa-tions on their probability distribution (see also App. C),which can be solved with arbitrary numerical precisionusing a population dynamics algorithm for each valueof Γ, ε and N .The Lyapunov exponent Λ describes the exponentialgrowth or the exponential decrease of the imaginary partof the diagonal elements of the self-energy with the num-ber of recursion steps r as: ∆ typ ∝ e Λ r . The most accu-rate way to compute the value of Λ is provided by the “in-flationary” algorithm put forward in Ref. [64]. The ideais to include an additional step to the recursion Eqs. (24)where all the ∆’s of the population are multiplied by afactor e − Λ r so to keep the typical imaginary part fixedand small: ∆ typ = θ (with θ (cid:28) Q ( S, ∆) is reached in this recursiveprocedure, Λ r → Λ. In practice one has to consider thelimit θ → M → ∞ , where M is the number of el-ements of the population. In fig. 7 we report the resultsof an accurate numerical computation of the Lyapunovexponent for Γ = 0 . ε and of the connectivity N , performedwith populations’ size M ranging from 2 to 2 , andwith θ from 10 − to 10 − (and for n = 2).One observes that the critical energy “density” ε loc ( N )at which the Lyapunov exponent vanishes slowly butcontinuously decreases as N is increased. In Fig. 8 weplot the dependence of ε loc ( N ) as a function of N forthree different values of Γ: We find a similar behaviorfor all values of the transverse field, at least in the re-1 ε Λ N=32N=64N=96N=128N=160N=192N=224N=256N=320N=384N=512N=768N=1024
Figure 7. Lyapunov exponent Λ as a function of the energy ε = E/N of the hybrid Anderson tight-binding models when V is sent to infinity at fixed N , for several values of N rangingfrom 32 to 1024 and for Γ = 0 .
2. The results are obtainedwithin the cluster expansion with n = 2. The vertical graydotted line corresponds to the prediction of the FSA for thelocalization threshold ε loc = Γ (see App. C).
100 1000 N | d Λ / d ε | l o c
100 1000 N ε l o c Γ =0.1 Γ =0.2 Γ =0.4 Figure 8. Main panel: Mobility edge ε loc ( N ) (such thatΛ( ε loc ) = 0) of the family of models (3) when the thermody-namic limit is taken from the start, as a function of N forthree different values of Γ (and for n = 2). The continuouscurves correspond to the analytical prediction of Eq. (25),with the fitting parameters c and c smoothly varying withΓ as c ≈ − . − .
27, and − .
62, and c ≈ .
25, 4 . . .
1, 0 .
2, and 0 . ε loc = Γ. Inset: Slope of the Lyapunov exponent computedat ε loc as a function of N for the same values of Γ as in themain panel. The black dashed line is a pwer-law fit of thedata as | dΛ / d ε | ε loc (cid:39) A N γ , with γ ≈ . gion of the phase diagram where the physical proper-ties of the QREM are dominated by the random term,Γ < √ / In particular, we observe that ε loc ( N )seems to vanish in the large N limit as (cid:112) log N/N . Con-sistently, we find that the slope | dΛ / d ε | ε loc around ε loc grows roughly as √ N (inset). This behavior can be understood in terms of the ana-lytic computation, carried out in App. C in full details,of the largest eigenvalue of the integral operator associ-ated to the self-consistent equation for Q ( S, ∆) in thelarge-connectivity limit and for n = 0 (i.e., in the stan-dard single-site Bethe approximation, when the underly-ing lattice is taken as a RRG of connectivity N ), whichyield: ε loc = log (cid:104)(cid:112) N Γ /π (cid:0) log (cid:0) N/ Γ (cid:1) + c (cid:1)(cid:105) + c N / , (25)where c is a real constant of O (1) which only depend onΓ, and c can be expressed in terms of the solution µ (cid:63) ofthe self-consistent equations (C9) for the real part of theself-energy. The predictions of this equations are plottedin Fig. 8 on top of the numerical points, showing a goodagreement with the numerical results.It is interesting to compare the asymptotic behavior atlarge N of the localization threshold ε loc found here forthe family of auxiliary models with the large connectiv-ity limit of the standard Anderson model on the Bethelattice. In fact, for Bethe lattices of connectivity k + 1 and on-site random energies taken from a box dis-tribution of width W , in the large k (and large disorder)limit the localization transition (at zero energy) takesplace at disorder W L given by:4 ρ log( W L ) = 1 k , (26)where ρ is the density of state in the middle of thespectrum which, at strong disorder, is just given by ρ (0) (cid:39) /W . In order to translate this relation to ourcase, assuming that that ρ ( E ) ≈ P ( E ), one finds thatthe exponential dependence on N of the DoS is exactlycanceled for E loc = N ε loc given by Eq. (25): ρ ( ε loc ) ≈ e − c N Γ log( N/ Γ ) . (We have neglected the constant c for simplicity.) Thevariance of the random energies of the QREM scales as σ E = N/
2, which leads us to identify the effective disor-der as W ≈ √ N . Thus, from Eq. (26) one gets:4 e − c log √ NN Γ log( N/ Γ ) ≈ N , which is satisfied for c = log(2 / Γ).Going back now to the QREM and to its Hilbert spaceformulation (3) defined on finite boolean hypercubes of V = 2 N sites, we argue that the mobility edges of theauxiliary models provide asymptotically in the large N limit an estimation for the transition line between thefully ergodic ETH phase and the delocalized but multi-fractal one. The argument goes as follows: A naturalway to construct many-body multifractal eigenstates in2the Hilbert space is to assume that the wave-functions areexponentially localized on the N -dimensional hypercubewith a finite, N -independent, localization length ξ . Thenumber of sites occupied by such states scales with N roughly as ξ N , which corresponds to a fractal dimension D roughly equal to log ξ . If the volume of the system(2 N ) is sent to infinity concomitantly with the spatialdimension ( N ), these wave-functions will show a mul-tifractal behavior in the thermodynamic limit, as theytend to occupy an exponentially diverging volume with N , yet exponentially smaller than the total volume. If,instead, the thermodynamic limit is taken first keepingthe dimensionality N fixed, these states will appear asgenuinely Anderson-localized, as they tend to occupy afinite, O (1), volume.On the other hand, the FSA analysis of the linearizedrecursion relations of the auxiliary models on the Bethelattice, which consists in neglecting the real part of theself-energy in the denominators of Eqs. (24), simply pre-dicts (again in the simplest n = 0 setting) that ε FSAloc = Γ,irrespectively of N (see App. C for a detailed calculation).This is the same result for the many-body mobility edgeof the QREM at the lowest order in Γ . We argue thatthe localization threshold predicted by the FSA does notdepend on whether the V → ∞ limit is taken beforethe N → ∞ one or not. In fact, the FSA only keepsthe leading-order contribution to the wave-function am-plitude at each site, and determine the convergence ofthe perturbative expansion by counting the relative num-ber of resonances found at a given distance, comparedto the total number of sites accessible at such distance.In this sense, this approximation captures the transitionfrom the Anderson-localized regime (where the perturba-tive expansion is convergent, the eigenstates are weakly-dressed single configurations of spins and occupy a finitevolume on the hypercube) to a delocalized regime (wherethe perturbative expansion does not converge, and res-onances can be found at arbitrary large distances) irre-spectively of the multifractal nature of the eigenstates.These arguments thus suggest that, while ε FSAloc gives arough estimate of the mobility edge between the MBLphase and the NEE phase, ε loc given in Eq. (25) pro-vides an estimation of the transition between the fullyergodic ETH phase and the delocalized non-ergodic one: ε MBL ≈ Γ ,ε ETH ≈ (cid:34) log (cid:112) N Γ /π + O (log(log N )) N (cid:35) / . (27)According to this interpretation, fully ergodic eigenstatesof the QREM are only found in a narrow energy windowaround | ε | < ε ETH , which concentrate around zero inthe thermodynamic limit (in agreement with Refs. [28]and [ ). Yet, the fraction of ergodic eigenfunctions atlarge N is approximately given by (cid:82) Nε ETH − Nε ETH ρ ( E ) d E ≈ − C/ √ N , with ρ ( E ) ≈ e − E /N / √ πN , and C being aconstant of order 1 which depend on Γ. Hence, although ε ETH →
0, due to the scaling of the many-body energieswith N , only a fraction of order 1 / √ N of the 2 N eigen-states are delocalized but non-ergodic. Furthermore, afinite value of the mobility edge, ε MBL ≈ Γ, only corre-sponds to an exponentially small fraction of Anderson-localized eigenstates in the tails of the DoS.
VIII. CONCLUSIONS AND PERSPECTIVES
In this work we have revisited the dynamical phasediagram of the QREM, using a complementary set ofapproaches, the FSA coupled to a mapping to the RPmodel and the self-consistent theory of localization extended to include the local Hilbert space structure ofthe QREM. While the FSA is expected to yield a betterestimation of correlations and large distance physics thelatter provides a quite accurate approximation of localobservables. These approaches provide a qualitativelysimilar scenario for the phase diagram in the energy-transverse field plane, namely the existence of three dy-namical phases: a fully ergodic delocalized one, an inter-mediate non-ergodic extended regime with multifractalbehavior, and an Anderson localized one. For what con-cerns the quantitative features of the phase diagram andthe properties of the three phases there remain howevermany open questions. The FSA and our RP mappingseem to suggest that ergodic delocalized states exist inan entire region around zero energy density (see also ),while the analysis of the quantum cavity equations sug-gest that if such a region exists it is much narrower inenergy, in agreement with Refs. [28] and [34]. On theother hand, an approximate analytic solution of the cav-ity equations, corresponding to an auxiliary Andersonmodel on a Bethe lattice where the connectivity is sentto infinity after the thermodynamic limit is taken, mim-icking the exponential scaling of the number of sites ofthe hypercube, points toward a threshold energy for fulldelocalization squeezing to zero in the thermodynamiclimit. It is worth stressing upon concluding that certain fea-tures of the QREM make it very peculiar, and producesome specific and unique features compared to othergeneric interacting disordered models such as 1 d disor-dered spin chains or mean-field spin glass models. Firstamong all the absence of correlations between the many-body energies E i ’s and the spin configurations { σ za } . Theother unique feature of the QREM is the fact that inthe frozen glassy phase, T < T K , the Edwards-Andersonorder parameter is equal to one, implying that essen-tially no spin can be flipped with respect to the initialstate. This property implies that in the MBL phase of theQREM many-body wavefunctions are genuinely Ander-son localized and occupy a finite volume in the Hilbertspace, while for generic interacting models one expectsthat the volume occupied by many-body eigenstates inthe configuration space is subexponentially large due tothe presence of a finite fraction of active spins. This3makes interesting and worth pursuing the investigationusing techniques developed here of other disordered meanfield models.
ACKNOWLEDGMENTS
We warmly thank L. Ioffe for helpful discussions.This work was partially supported by the grant fromthe Simons Foundation ( ◦ q q q N- pq N- . . .
12 12
Figure 9. Schematic representation of the recursion stepwhich yields the self-consistent equations for the 2 × n = 1). Appendix A: Recursion relations for the matrixelements of the resolvent within the clusterapproximation
For n = 1 the clusters are simply made by two sites(corresponding to two spin configurations which differby a single spin flip) connected by an edge (Fig. 9). Thecavity resolvent matrix on such cluster can then just beparametrized by three complex numbers: G p → q = (cid:18) g p → q g p → q g p → q g p → q (cid:19) . The matrix elements of the Hamiltonian (3) on the sitesof the cluster are: H c = (cid:18) − E − Γ − Γ − E (cid:19) . Thus, Eq. (16) becomes: g p → q k det G p → q k = − E − z − Γ (cid:88) q l ∈ ∂p/q k g q l → p ,g p → q k det G p → q k = − E − z − Γ (cid:88) q l ∈ ∂p/q k g q l → p ,g p → q k det G p → q k = − Γ + Γ (cid:88) q l ∈ ∂p/q k g q l → p , where det G p → q k = g p → q k g p → q k − ( g p → q k ) . FromEq. (17), one can then write down the equations for the4elements of the resolvent matrix on the cluster: g p det G p = − E − z − Γ (cid:88) q k ∈ ∂p g q k → p , g p det G p = − E − z − Γ (cid:88) q k ∈ ∂p g q k → p , g p det G p = − Γ + Γ (cid:88) q k ∈ ∂p g q k → p , where g p and g p are the diagonal elements on sites i and j of the cluster p of the resolvent matrix, g p is the off-diagonal element, and det G p = g p g p − ( g p ) When g p → q k and g p → q k are set to zero, these equations give backthe standard (cavity) recursion equations for the single-site Anderson model on the Bethe lattice. Moreover,since the off-diagonal elements are proportional to Γ, inthe limit of small transverse field one might expand theseequations in powers of Γ to obtain the systematic cor-rections to the zeroth-order equations due to the smallloops: g p → q k , = − E − z − Γ (cid:88) q l ∈ ∂p/q k g q l → p , ,g p → q k , = − E − z − Γ (cid:88) q l ∈ ∂p/q k g q l → p , ,δg p → q k , = − Γ g p → q k , g p → q k , + O (Γ ) ,δg p → q k , = Γ (cid:0) g p → q k , (cid:1) g p → q k , + O (Γ ) ,δg p → q k , = Γ (cid:0) g p → q k , (cid:1) g p → q k , + O (Γ ) , where g p → q k , and g p → q k , are the diagonal elements of theresolvent at the 0-th order of the cluster expansion (i.e.,within the standard single-site Bethe approximation),and δg p → q k , , δg p → q k , , and δg p → q k , are the corrections for n = 1 up to the lowest order in Γ.One can proceed in a similar way for n = 2 and obtainclosed equations for the ten independent elements of thecavity resolvent matrices on each cluster in terms of theelements of the cavity resolvent matrices on the neigh-boring clusters. However the equations are much longerand we do not write them here explicitly. It is just worthto mention that in this case the off-diagonal elements ofthe resolvent matrix between pairs of sites of the clusterthat are connected by an edge on the hypercube (e.g., g p → q k , g p → q k , g p → q k , and g p → q k , using the notation offig. 3) are proportional to Γ. Conversely, the Green’sfunctions between pairs of sites that are not connectedby an edge on the hypercube (e.g., g p → q k and g p → q k ) are(as expected) proportional to Γ .
1. Spectrum of the kinetic term
It is instructive to study how the spectrum of the ki-netic term of the Hamiltonian is modified by the cluster ε ρ I ε ρ I N=16N=32N=64
Figure 10. Comparison between the exact DoS of the de-localizing interacting part of the QREM [i.e., Γ times theadjacency matrix of the N dimensional hypercube, Eq. (12),dotted red, blue, and green curves] and the spectra of the ki-netic term of Eq. (3) within the cluster approximation for n = 0 [i.e., the standard single-site Bethe approximation,which yields (Γ times) the adjacency matrix of a RRG ofconnectivity N , Eq. (15), dashed orange, violet, and darkgreen curves], n = 1 (dashed-dotted violet, and dark greencurves), and n = 2 (continuous violet, and dark green curves)for Γ = 0 . N . At each order of the clus-ter expansion the edge of the spectrum is shifted by Γ to theright (to the leading order in N ). The imaginary regulator η is set to be η = cδ , Eq. (B1), with δ = 1 / (2 N ρ ( E )). expansion. To this aim we consider the pure limit (in ab-sence of disorder) of the equations above for n = 1. In thepure case ( E i = E j = 0) the hypercube is translationallyinvariant. One can thus look for a uniform solution ofthe equations in the form: gg − h = − z − ( N − g ,gg − h = − Γ + ( N − h , where g = g p → q k = g p → q k and h = g p → q k for all p, q .One can then introduce the variables g + = g + h and g − = g − h in terms of which the equations above become:1 g + = − z + Γ − ( N − g + , g − = − z − Γ − ( N − g − , which coincide with the equations that one obtains for thestandard single-site Anderson model on the Bethe latticein the uniform limit, with energies shifted of ± Γ. TheDoS is thus modified accordingly. In particular the edgesof the spectrum are shifted as well by +Γ on the rightedge and by − Γ on the left edge. One can indeed define g = g p = g p and h = g p for all p , and introduce thevariables g ± = g ± h , which verify the following equations:1 g ± = − z ± Γ − ( N − g ± , ρ ( n =1) I =Im( g + + g − ) / (2 π ).For n = 2 a similar treatment of the equations yields:1 g ± = − z ± − ( N − g ± , g = − z − ( N − g , where g ± = g + f ± h and g = g − f , where g = g p → q k u , h = g p → q k = g p → q k = g p → q k = g p → q k , and f = g p → q k = g p → q k for all u, p, q . Similarly one defines g = g pu , h = g p = g p = g p = g p , and f = g p = g p for all u, p , andintroduces the variables g ± = g + f ± h and g = g − f ,which verify the following equations:1 g ± = − z ± − ( N − g ± , g = − z − ( N − g , in terms of which the DoS reads: ρ ( n =2) I = Im (cid:20) g + + g − + 2 g π (cid:21) . A comparison between the exact spectrum of the kineticterm of the Hamiltonian (3), i.e., the adjacency matrixof the N -dimensional hypercube (12), the DoS result-ing from the cluster approximation for n = 0, i.e., theadjacency matrix of a RRG of connectivity N , Eq. (15), n = 1, and n = 2 is shown in fig. 10 for Γ = 0 . N . The imaginary regulator is set to the valueused to solve the recursion equations within the C-BPapproximation in presence of the disordered many-bodyon site energies, Eq. (B1). At the order n of the clusterexpansion the (right) edge of the spectrum is shifted by+ n Γ ( − n Γ) to the right (left) to the leading order in N . Appendix B: Solution of the recursion equations andcomparison with exact diagonalization
As explained in the main text, the cluster approxi-mation allows us to derive a system of closed equations,Eqs. (16) and (17), for the diagonal elements of the re-solvent matrix of (3). A first, and crucial, question thatwe want to address here is to what extent this approxi-mation provides a good qualitative and quantitative de-scription of the spectral statistics of the QREM. To thisaim in this section we consider samples of moderate sizeand compare the probability distributions of the LDoScomputed from the numerical solution of the recursionrelations (16) and (17) with those obtained from EDs ofthe QREM.There are essentially two ways, that we detail below,to solve the recursion equations for the Green’s functionand obtain information on the spectral statistics at finite N .
1) C-BP on the hypercube.
The most accurate strat-egy, to which we will refer to as “Cluster Belief Propa-gation” (C-BP) algorithm (see Ref. [54] for a detailedexplanation of this approach for the usual tight-bindingAnderson model on the Bethe lattice), is to solve directlyEqs. (16) and (17) on random realizations of the hyper-cube of 2 N sites (i.e., N spins). In practice, one proceedsas follows: • One first generates a random instance of the hy-percube drawing the 2 N on-site energies from thedistribution (2); • One finds a partition of the hypercube in 2 N − n clus-ters of 2 n sites each (note that the choice of thepartition is not unique); • Then one finds the fixed point of Eqs. (16), whichconstitute a system of ( s + 1)( N − n )2 N − coupledequation for the s ( s + 1) / s +1)( N − n )2 N − ); • Using Eqs. (17) one obtains the s ( s + 1) / • One then repeat this procedure several times to av-erage over different realizations of the on-site dis-order (and over different choices of the cluster par-titioning).As discussed above (see also [54]), in order for the recur-sive equations to converge to the physical fixed point, thebroadening η must be larger than the mean level spacing.Hence, in order to implement the η → + limit correctly,for any given choice of the parameters Γ, ε and N , theimaginary regulator is self-consistently set to be a con-stant of order 1 times the mean level spacing: η = c N ρ η ( ε ) = cπ (cid:80) V i =1 Im G i ( N ε + iη ) . (B1)As shown in [54], for large enough system sizes and for c large enough, ρ η converges to its asymptotic value ob-tained in the limit V → ∞ and η → + . We will take c = 64 throughout. Within the C-BP approach onecan study hypercubes of sizes up to N = 26, which areconsiderably larger than the ones accessible via the mostefficient ED algorithms.
2) C-PD algorithm on the RRG.
In order to accesseven larger system sizes, one can adopt another strategy,to which we will refer hereafter as “Cluster PopulationDynamics” (C-PD) algorithm, which consists in inter-preting the recursion relations for the Green’s functionsas equations for their probability distributions once theaverage over the disorder is taken. In fact, since G p → q G p are random matrices, one can assume that av-eraging over the on-site random energies leads to func-tional equations on their probability distribution Q ( G )and Q ( G ). From Eq. (16) we naturally get: Q ( G ) = (cid:90) n (cid:89) u =1 d P ( E u ) N − n − (cid:89) q =1 d Q ( G q ) × δ (cid:32) G − + H c + z I s + Γ N − n − (cid:88) q =1 G q (cid:33) , (B2)where P ( E ) is given by Eq. (2), H c is the Hamiltonian (3)on the sites of the cluster (see, e.g., Appendix A for theexplicit expression of H c for n = 1), which contains the s diagonal random energies E , . . . , E s , and I s is the s × s identity matrix. (The notation d Q ( G q ) is just a short-cut for the integration over the s ( s + 1) / G q .) Once the fixed point of this equation isobtained, using Eq. (17) one can find an equation for theprobability distribution of the elements of the resolvent: Q ( G ) = (cid:90) n (cid:89) u =1 d P ( E u ) N − n (cid:89) i = q d Q ( G q ) × δ (cid:32) G − + H c + z I s + Γ N − n (cid:88) q =1 G q (cid:33) . (B3)As before, z = N ε + iη and the imaginary regulator isself-consistently set to be c times the mean level spacing: η = 2 − N cπ/ (cid:104) Im G(cid:105) , where the average is performed overthe distribution Q ( G ). This set of functional equationscan be solved numerically with an arbitrary degree ofprecision using a population dynamics algorithm .Hereafter we will show results obtained using populationsof M fields going from M = 2 to M = 2 (and for n = 2). Note that, differently from the C-BP approach,within the C-PD approximation the specific structure ofthe hypercube is completely lost for distances larger thanthe size of the clusters (apart from the local connectivityof each cluster equal to N − n ).On the other hand, from ED we can easily obtainthe matrix elements G i for a given instance of the QREMin terms of the eigenvalues E α and the eigenvectors | α (cid:105) of (1) as: G i ( N ε + iη ) = N (cid:88) α =1 |(cid:104) α | i (cid:105)| E α − N ε + iη ( E α − N ε ) + η . (B4)For each choice of the parameters Γ, N , and ε , the imag-inary regulator is set to the same value as the one usedto solve the self-consistent recursion relations, Eq. (B1).In fig. 11 we focus on the probability distribution ofthe imaginary part of the Green’s functions, and plot Q (log Im G ) for Γ = 0 . .
2, for N ranging from10 to 15, and for two values of ε which are supposed to beon the ergodic side of the MBL transition and close to the -8 -4 0 4 log ImG P ( l og I m G ) -8 -4 0 4 log ImG P ( l og I m G ) -8 -4 0 4 log ImG P ( l og I m G ) -8 -4 0 4 log ImG P ( l og I m G ) Γ=0.2ε=0.2 Γ=0.2ε=0.32ε=0.16Γ=0.1Γ=0.1ε=0.07
Figure 11. Probability distributions P (log Im G ) obtainedfrom ED (filled symbols), C-BP ( n = 2, continuous curves),and C-PD ( n = 2, dashed curves) for N = 8 (orange), 10(red), 11 (violet), 12 (maroon), 13 (blue), 14 (turquoise), and15 (green), and for Γ = 0 . ε = 0 . . ε = 0 .
32 (top right panel), Γ = 0 . ε = 0 . . ε = 0 .
16 (bottom rightpanel). mobility edge respectively.
In all cases, we observea good agreement between the probability distributionsfound from EDs and Eq. (B4), and their C-BP counter-part, found from the numerical solution of the recursionrelations (16) and (17) on the hypercube for n = 2. Inthe figure we also plot the distributions Q (log Im G ) ob-tained from the C-PD algorithm, Eqs. (B2) and (B3),which presents very small deviations from the C-BP re-sults only in the very far tails of the distributions at smallIm G , and only visible for some values of Γ, ε , and N . Appendix C: Analytic computation of thelocalization threshold of the auxiliary Andersonmodels in the large- N limit In this appendix we discuss the analytical computationof the localization threshold(s) of the family of auxiliaryAnderson tight-binding models described by the Hamil-tonian (3), when the thermodynamic limit,
V → ∞ , istaken from the start, while keeping N fixed. For simplic-ity, we will only consider the simplest setting n = 0, i.e.,the standard single-site Bethe approximation in whichthe hypercube is approximated by a tree-like structureof connectivity N .
1. Probability distribution of the real part of theself-energy
The fisrt step is to realize that the recursion relationsfor the real part of the self-energies in the linearizedregime is independent on the imaginary part, and can7be solved as explained below. It is useful to introducethe variables X i → j (i.e., the real part of the diagonal el-ements of the resolvent matrix in the linearized regime)defined as: X i → j = − E i + N ε + S i → j = G Ri → j . (C1)In terms of these variables at the 0-th order of the clusterexpansion Eqs. (24) become: S i → j = Γ (cid:88) j (cid:48) ∈ ∂i/j X j (cid:48) → i , ∆ i → j = Γ (cid:88) j (cid:48) ∈ ∂i/j X j (cid:48) → i ∆ j (cid:48) → i . (C2)Hence, the probability distribution of the real part ofthe self-energy R S ( S ) can be obtained in terms of theprobability distribution R X ( X ): R S ( S ) = (cid:90) N − (cid:89) i =1 d X i R X ( X i ) δ (cid:32) S − Γ (cid:88) i X i (cid:33) = (cid:90) d k π e ikS (cid:34)(cid:89) i d x i d k i π e i ( k i − Γ k ) x i ˆ R X ( k i ) (cid:35) , where ˆ R X ( k ) is the characteristic function of R X ( X ).Assuming that at small k it behaves as the characteristicfunction of a Cauchy distribution,ˆ R X ( k ) (cid:39) − A | k | − ikµ , (C3)implies that in the large N limit also R S ( S ) is given bya Cauchy distribution:ˆ R S ( k ) = (cid:104) ˆ R X (Γ k ) (cid:105) N − (cid:39) e − ( N − A Γ | k |− i ( N − µk ,R S ( S ) = A S π (cid:104) ( S − µ S ) + A S (cid:105) , (C4)with A S = N Γ A , and µ S = N Γ µ . (C5)(Throughout we will consider the large N limit N − ≈ N .) On the other hand, from Eq. (C1) we have that: R X ( X ) = (cid:90) d E d S P ( E ) R S ( S ) δ (cid:18) X + 1 E + N ε + S (cid:19) = 1 | X | (cid:90) d E P ( E ) R S (cid:18) − E − N ε − X (cid:19) . (C6)After some simple algebra we obtain that: R S (cid:18) − E − N ε − X (cid:19) = cX π (cid:104) ( X − X ) + c (cid:105) , (C7) with c = A S ( E + N ε + µ S ) + A S ,X = − E + N ε + µ S ( E + N ε + µ S ) + A S . As a result, from the second line of Eq. (C6) and fromthe relations above, we get: R X ( X ) = (cid:90) d E P ( E ) cπ (cid:104) ( X − X ) + c (cid:105) . (C8)We can now finally compute self-consistently the char-acteristic function of R X ( X ) by expanding the equationabove up to first order in k :ˆ R X ( k ) = (cid:90) d E P ( E ) e − c | k |− ikX (cid:39) − (cid:90) d E P ( E ) ( c | k | + ikX ) . From Eqs. (C3), (C5), and the last equation we can thusobtain two self-consistent relations for the coefficients A and µ : A (cid:63) = (cid:90) d E P ( E i ) N Γ A (cid:63) ( E + N ε + N Γ µ (cid:63) ) + ( N Γ A (cid:63) ) ,µ (cid:63) = − (cid:90) d E P ( E ) E + N ε + N Γ µ (cid:63) ( E + N ε + N Γ µ (cid:63) ) + ( N Γ A (cid:63) ) . (C9)These equations can be easily solved numerically. Infig. 12 we show the solutions A (cid:63) and (minus) µ (cid:63) ofEqs. (C9) for Γ = 0 . ε = 0 . ε = 0 . N . While A (cid:63) de-cay very fast (exponentially) with N , µ (cid:63) decreases muchslowlyer, roughly as 1 /N (blue dashed lines). The verticalthick gray dashed lines indicate the localization thresholdwhere the Lyapunov exponent vanishes for these partic-ular values of Γ and ε (see fig. 7).In fig. 13 we plot the probability distribution of the realpart of the self-energy, R S ( S ), for Γ = 0 . ε = 0 .
2, andseveral values of N across the localization threshold. Wefocus on the negative real axis since the peak of the distri-bution is located in S = µ S which turns out to be nega-tive. Filled symbols correspond to the numerical solutionfound using the C-PD algorithm (for n = 0) of the lin-earized recursion equations for the self-energy (24), whilethe continuous lines correspond to the analytic predic-tion (C4), with A S and µ S given by Eqs. (C5) and (C9).The agreement between the numerical results and theanalytic solution is excellent, and it improves for large N .
2. Computation of the Lyapunov exponent
Once the probability distribution of the real part ofthe self-energy has been obtained, we can focus on the8
32 64 128 256 N A * , - µ ∗ A * - µ ∗
16 64 256 1024 4096 N A , µ Figure 12. Solutions A (cid:63) (red) and (minus) µ (cid:63) (blue) ofEqs. (C9) for Γ = 0 . ε = 0 . ε = 0 . N . The blue dashed lines cor-respond to µ (cid:63) ∼ − /N . The vertical thick gray dashed linesindicate the localization threshold where the Lyapunov ex-ponent vanishes for these particular values of Γ and ε (seefig. 7). -S R S (- S ) N=32N=64N=96N=128N=160
Figure 13. Probability distributions R S ( − S ) for Γ = 0 . ε = 0 .
2, and several values of N across the localization tran-sition of the auxiliary models. Filled symbols are obtainedas the numerical solution of the linearized recursion rela-tions for the self energy with the C-PD algorithm for n = 0,while continuous lines correspond to the analytic predictionof Eqs. (C4), (C5), and (C9). integral equation for the joint distributions of the real and the imaginary part (for n = 0): Q ( S, ∆) = (cid:90) N − (cid:89) i =1 [d E i P ( E ) d S i d∆ i Q ( S i , ∆ i )] × δ (cid:32) S + Γ (cid:88) i E i + N ε + S i (cid:33) × δ (cid:32) ∆ − Γ (cid:88) i ∆ i ( E i + N ε + S i ) (cid:33) . (C10)We replace the δ -functions by their integral representa-tion in the Fourier space and also write Q ( S i , ∆ i ) as theinverse Fourier transform of ˆ Q ( S i , k i ) with respect tothe second argument, defined as:ˆ Q ( S, k ) = (cid:90) + ∞−∞ d∆ e − ik ∆ Q ( S, ∆) , (C11)yielding:ˆ Q ( k , k ) = (cid:90) (cid:89) i =1 (cid:20) d E i P ( E i ) d S i d∆ i d k i π ˆ Q ( S i , k i ) × e ik Γ / ( E i + Nε + S i ) e ∆ i [ k i − k Γ / ( E i + Nε + S i ) ] (cid:21) . We can now perform the integration over d∆ i , whichgives 2 πδ ( k i − k Γ / ( E i + N ε + S i ) ), and then integrateover k i :ˆ Q ( k , k ) = (cid:20) (cid:90) d E P ( E ) d S ˆ Q (cid:18) S, k Γ ( E + N ε + S ) (cid:19) × e ik Γ / ( E + Nε + S ) (cid:21) N − . (C12)Similarly to Ref. [38], we assume that in the localizedphase the following asymptotic form of Q ( S, ∆) holds forlarge enough ∆: Q ( S, ∆) (cid:39) A ( S )∆ β for ∆ → ∞ , ˆ Q ( S, k ) (cid:39) ˆ Q ( S, − α | k | β A ( S ) for k → , (C13)where ˆ Q ( S,
0) is by definition the marginal of Q ( S, ∆)once we integrate over ∆, i.e., Q ( S, ∆) = R S ( S ).Plugging the asymptotic form (C13) into both sides ofEq. (C12) we obtain:ˆ R S ( k ) − α | k | β ˆ A ( k ) (cid:39) (cid:20) (cid:90) d E P ( E ) d S (cid:18) R S ( S ) − α (cid:12)(cid:12)(cid:12)(cid:12) k Γ ( E + N ε + S ) (cid:12)(cid:12)(cid:12)(cid:12) β A ( S ) (cid:19) e ik Γ / ( E + Nε + S ) (cid:21) N − (C14)We can now expand the r.h.s of the equation above inpowers of k up to the order | k | β and define: I = (cid:90) d E P ( E ) d S R S ( S ) e ik Γ / ( E + Nε + S ) ,I = (cid:90) d E P ( E ) d S Γ β | E + N ε + S | β A ( S ) e ik Γ / ( E + Nε + S ) . I − α | k | β I ) N − (cid:39) I N − [1 − α ( N − | k | β I /I ] = I N − − α ( N − | k | β I N − I . From Eq. (24) we have that by defini-tion I N − = ˆ R S . Hence, in the large- N limit, I N − (cid:39) I N − = ˆ R S ( k ), we get:ˆ A ( k ) (cid:39) N Γ β ˆ R S ( k ) (cid:90) d E P ( E ) d S × A ( S ) e ik Γ / ( E + Nε + S ) | E + N ε + S | β . Changing variable to w = E + N ε + S , replacing A ( S )by the inverse Fourier transform of ˆ A ( k ) and integratingover d E we obtain:ˆ A ( k ) (cid:39) N Γ β ˆ R S ( k ) (cid:90) d w d k π e − Nk − ikNε × ˆ A ( k ) e ikw + ik Γ /w | w | β . (C15)For a given choice of the parameters Γ, ε , and N , the lo-calization threshold at the 0-th order of the cluster expan-sion is thus given by the value of the energy ε loc (Γ , ε, N )such that the largest eigenvalue λ β of the integral op-erator defined by the equation above, becomes equal toone. As first noticed in Ref. [38] (see also Refs. [65], [66],and [70]) the kernel of the integral operator is symmetricaround β = 1 / β → − β ,which implies that λ = 1 if and only if β = 1 /
2. Sincethis is the value of interest for the transition, hereafterwe will focus on the case β = 1 / w can then be performed in terms of modified Besselfunctions: (cid:90) d w e ikw + ik Γ /w | w | = − πY (cid:16) (cid:112) kk (cid:17) . Following Ref. [65] we now assume that in the large con-nectivity limit the eigenvector of the integral operatordefined by Eq. (C15) for β = 1 / R ( k ). Assuming that the localization transi-tion occurs on such energy scales (apart from logarithmiccorrections) ε = ˜ ε √ N , (C16)with ˜ ε of O (1), the equation for the mobility edge be-comes:1 = N Γ (cid:90) d k π e − Nk − N Γ A (cid:63) | k |− ik ( √ N ˜ ε + N Γ µ (cid:63) ) × (cid:104) − πY (cid:16) (cid:112) kk (cid:17)(cid:105) , where we have used Eqs. (C5). A (cid:63) is exponentially smallin N and can be neglected, while N µ (cid:63) is of order 1 atthe transition (see fig. 12), and gives a correction of order1 / √ N to ˜ ε as ˜ ε (cid:48) = ˜ ε +Γ √ N µ (cid:63) . Since the integral over k is cut-off on a scale 1 / √ N we can then expand the Besselfunction keeping only the leading logarithmic divergenceat small k , Y ( x ) ≈ k ) /π . In the N → ∞ limit wecan then change variable to ˜ k = √ N k , yielding:1 = √ N Γ √ π e − (˜ ε (cid:48) ) (log N − O (1)) . Putting everything together we finally obtain the equa-tion for the mobility edge (25) given in the main text.
3. The Forward-Scattering Approximation
The FSA consists in neglecting the real part of theself-energy in the denominators of Eqs. (24). FromEq. (C6), one can then compute the probability distri-bution R X ( X ) for n = 0: R X ( X ) = 1 | X | e − N ( X − Nε ) √ πN , which does not verify exactly the asymptotic form (C3)for its behavior at large X in presence of the real part S . This implies that, differently from what happensfor the large connectivity limit of the usual Andersontight-binding model on the Bethe lattice, the distribution R S ( S ) found within the FSA does not coincides exactlywith the distribution obtained in Sec. C 1 in presence ofthe real parts in the denominators.Following the steps of the calculation detailed abovefor the largest eigenvalue of the linearized recursion rela-tions, one can compute the distribution of the imaginarypart of the self-energy as: Q (∆) = (cid:90) N − (cid:89) i =1 [d E i P ( E ) ∆ i Q (∆ i )] × δ (cid:32) ∆ − Γ (cid:88) i ∆ i ( E i + N ε ) (cid:33) , which gives the following self-consistent equation for itsFourier transform:ˆ Q ( k ) = (cid:20)(cid:90) d E P ( E ) ˆ Q (cid:18) k Γ ( E + N ε ) (cid:19)(cid:21) N − . Assuming, as before, the asymptotic form ˆ Q ( k ) ≈ − α | k | β , and expanding the last equation at large N , onegets the equation for the localization threshold within theFSA: 1 ≈ N Γ β (cid:90) d E P ( E ) 1 | E + N ε | β . (C17)The symmetry λ ( β ) = λ (1 − β ) is now lost. Hence thetransition point is not achieved at β = 1 /
2, but ratherat a given point β (cid:63) ∈ [0 , /
2] which depend on the other0parameters of the auxiliary models. Eq. (C17) can beeasily solved numerically for any choice of Γ, ε , N , and β , and gives the localization threshold ε FSAloc = Γ, with β (cid:63) → / N → ∞ . Indeed, since the random energy E is typically of order √ N , if ε is of order 1, one canexpand the denominator in powers of E/ ( N ε ) yielding:1 ≈ N − β (cid:18) Γ ε (cid:19) β (cid:18) β (2 β − N ε + . . . (cid:19) , which, in the N → ∞ limit and β → /
2, gives ε FSAloc = Γ. D. Basko, I. Aleiner, and B. Altshuler, Annals of Physics , 1126 (2006), cond-mat/0506617. B. L. Altshuler, Y. Gefen, A. Kamenev, and L. S. Levitov,Phys. Rev. Lett. , 2803 (1997), cond-mat/9609132. I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Phys. Rev.Lett. , 206603 (2005), cond-mat/0506617. E. Altman and R. Vosk, Annual Review of Condensed Mat-ter Physics , 383 (2015), arXiv:1408.2834 [cond-mat.dis-nn]. R. Nandkishore and D. A. Huse, Annual Review of Con-densed Matter Physics , 15 (2015), arXiv:1404.0686[cond-mat.stat-mech]. D. A. Abanin and Z. Papi´c, Annalen der Physik ,1700169 (2017), arXiv:1705.09103 [cond-mat.dis-nn]. F. Alet and N. Laflorencie, Comptes Rendus Physique ,498 (2018), quantum simulation / Simulation quantique,arXiv:1711.03145 [cond-mat.str-el]. D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev.Mod. Phys. , 021001 (2019), arXiv:1804.11065 [cond-mat.dis-nn]. M. Schreiber, S. S. Hodgman, P. Bordia, H. P. L¨uschen,M. H. Fischer, R. Vosk, E. Altman, U. Schneider, andI. Bloch, Science , 842 (2015), arXiv:1501.05661 [cond-mat.quant-gas]. P. Bordia, H. P. L¨uschen, S. S. Hodgman, M. Schreiber,I. Bloch, and U. Schneider, Phys. Rev. Lett. , 140401(2016), arXiv:1509.00478 [cond-mat.quant-gas]. J.-y. Choi, S. Hild, J. Zeiher, P. Schauß, A. Rubio-Abadal, T. Yefsah, V. Khemani, D. A. Huse, I. Bloch,and C. Gross, Science , 1547 (2016), arXiv:1604.04178[cond-mat.quant-gas]. D. J. Luitz, N. Laflorencie, and F. Alet, Phys. Rev. B ,081103 (2015), arXiv:1411.0660 [cond-mat.dis-nn]. D. J. Luitz and Y. Bar Lev, Annalen der Physik ,1600350 (2017), arXiv:1610.08993 [cond-mat.dis-nn]. N. Mac´e, F. Alet, and N. Laflorencie, Phys. Rev. Lett. ,180601 (2019), arXiv:1812.10283 [cond-mat.dis-nn]. M. Pino, L. B. Ioffe, and B. L. Altshuler, Proceedingsof the National Academy of Sciences , 536 (2016),arXiv:1501.03853 [cond-mat.dis-nn]. M. Pino, V. E. Kravtsov, B. L. Altshuler, and L. B. Ioffe,Phys. Rev. B , 214205 (2017), arXiv:1704.07393 [cond-mat.dis-nn]. E. J. Torres-Herrera and L. F. Santos, Annalen der Physik , 1600284 (2017), arXiv:1610.02035 [cond-mat.dis-nn]. V. Oganesyan and D. A. Huse, Phys. Rev. B , 155111(2007), cond-mat/0610854. A. Pal and D. A. Huse, Phys. Rev. B , 174411 (2010),arXiv:1010.1992 [cond-mat.dis-nn]. J. ˇSuntajs, J. Bonˇca, T. Prosen, and L. Vidmar, arXiv:1905.06345 [cond-mat.str-el] (2019). D. A. Abanin, J. H. Bardarson, G. D. Tomasi, S. Gopalakr-ishnan, V. Khemani, S. A. Parameswaran, F. Pollmann,A. C. Potter, M. Serbyn, and R. Vasseur, arXiv:1911.04501[cond-mat.str-el] (2019). W. De Roeck and F. Huveneers, Phys. Rev. B , 155129(2017), arXiv:1608.01815 [cond-mat.dis-nn]. D. J. Luitz, F. Huveneers, and W. De Roeck, Phys. Rev.Lett. , 150602 (2017), arXiv:1705.10807 [cond-mat.dis-nn]. P. Ponte, C. Laumann, D. A. Huse, and A. Chandran, Phil.Trans. R. Soc. A , 20160428 (2017), arXiv:1707.00004[cond-mat.dis-nn]. C. R. Laumann, A. Pal, and A. Scardicchio, Phys. Rev.Lett. , 200405 (2014), arXiv:1404.2276 [cond-mat.stat-mech]. C. L. Baldwin, C. R. Laumann, A. Pal, and A. Scardicchio,Phys. Rev. B , 024202 (2016), arXiv:1509.08926 [cond-mat.stat-mech]. C. L. Baldwin, C. R. Laumann, A. Pal, and A. Scardicchio,Phys. Rev. Lett. , 127201 (2017), arXiv:1611.02296[cond-mat]. L. Faoro, M. V. Feigel’man, and L. Ioffe, Annals of Physics , 167916 (2019), arXiv:1812.06016 [cond-mat.dis-nn]. B. Derrida, Phys. Rev. Lett. , 79 (1980). Y. Y. Goldschmidt, Phys. Rev. B , 4858 (1990). T. J¨org, F. Krzakala, J. Kurchan, and A. C. Maggs, Phys.Rev. Lett. , 147204 (2008), arXiv:0806.4144 [quant-ph]. C. L. Baldwin and C. R. Laumann, Phys. Rev. B ,224201 (2018), arXiv:1803.02410 [cond-mat.dis-nn]. T. Parolini and G. Mossi, arXiv:2007.00315 [cond-mat.dis-nn] (2020). V. Smelyanskiy, K. Kechedzhi, S. Boixo, H. Neven, andB. Altshuler, arXiv:1907.01609 [cond-mat.dis-nn] (2020). P. W. Anderson, Phys. Rev. , 1492 (1958). V. E. Kravtsov, I. M. Khaymovich, E. Cuevas, andM. Amini, New Journal of Physics , 122002 (2015),arXiv:1508.01714 [cond-mat.dis-nn]. D. Facoetti, P. Vivo, and G. Biroli, EPL (Europhysics Let-ters) , 47003 (2016), arXiv:1607.05942 [cond-mat.dis-nn]. R. Abou-Chacra, D. J. Thouless, and P. W. Anderson,Journal of Physics C: Solid State Physics , 1734 (1973). D. E. Logan and S. Welsh, Phys. Rev. B , 045131 (2019),arXiv:1806.01688 [cond-mat.str-el]. T. R. Kirkpatrick, D. Thirumalai, and P. G. Wolynes,Phys. Rev. A , 1045 (1989). M. Baity-Jesi, G. Biroli, and C. Cammarota, Journalof Statistical Mechanics: Theory and Experiment ,013301 (2018), arXiv:1708.03268 [cond-mat.dis-nn]. F. Pietracaprina, V. Ros, and A. Scardicchio, Phys. Rev.B , 054201 (2016), arXiv:1508.05097 [cond-mat.dis-nn]. E. Bogomolny and M. Sieber, Phys. Rev. E , 042116(2018), arXiv:1805.06675 [math-ph]. P. A. Nosov, I. M. Khaymovich, and V. E. Kravtsov, Phys.Rev. B , 104203 (2019), arXiv:1810.01492 [cond-mat.dis-nn]. V. E. Kravtsov, I. M. Khaymovich, B. L. Altshuler, andL. B. Ioffe, arXiv:2002.02979 [cond-mat.dis-nn] (2020). G. De Tomasi, M. Amini, S. Bera, I. M. Khay-movich, and V. E. Kravtsov, SciPost Phys. , 14 (2019),arXiv:1805.06472 [cond-mat.dis-nn]. P. Jacquod and D. L. Shepelyansky, Phys. Rev. Lett. ,1837 (1997), cond-mat/9706040. D. E. Logan and P. G. Wolynes, Phys. Rev. B , 4135(1987). D. E. Logan and P. G. Wolynes, The Jour-nal of Chemical Physics , 4994 (1990),https://doi.org/10.1063/1.458637. R. Bigwood, M. Gruebele, D. M. Leitner, and P. G.Wolynes, Proceedings of the National Academy of Sciences , 5960 (1998). A. De Luca and A. Scardicchio, EPL (Europhysics Letters) , 37003 (2013), arXiv:1206.2342 [cond-mat.str-el]. The properties of random-regular graphs have been exten-sively studied. For a review see ? . G. Biroli, G. Semerjian, and M. Tarzia, Progressof Theoretical Physics Supplement , 187 (2010),arXiv:1005.0342 [cond-mat.dis-nn]. G. Biroli and M. Tarzia, arXiv:1810.07545 [cond-mat.dis-nn] (2018). M. M´ezard and G. Parisi, The European Physical Jour-nal B - Condensed Matter and Complex Systems , 217(2001), cond-mat/0009418. L. Erd˝os, A. Knowles, H.-T. Yau, and J. Yin, Electron. J.Probab. , 58 pp. (2013), arXiv:1212.0164 [math.PR]. M. L. Mehta,
Random matrices , 3rd ed. (Academic Press,2004). F. L. Metz and I. P. Castillo, Phys. Rev. B , 064202(2017), arXiv:1703.10623 [cond-mat.dis-nn]. B. Altshuler, I. K. Zharekeshev, S. Kotochigova, andB. Shklovskii, Zh. Eksp. Teor. Fiz , 343 (1988). J. T. Chalker, V. E. Kravtsov, and I. V. Lerner, Journalof Experimental and Theoretical Physics Letters , 386(1996), cond-mat/9609039. E. Bogomolny and O. Giraud, Phys. Rev. Lett. ,044101 (2011), arXiv:1011.3686 [nlin.CD]. A. D. Mirlin, Physics Reports , 259 (2000), cond-mat/9907126. We have checked that varying c from 16 to 128 do notmodify the results. B. L. Altshuler, E. Cuevas, L. B. Ioffe, and V. E. Kravtsov,Phys. Rev. Lett. , 156601 (2016), arXiv:1605.02295[cond-mat.dis-nn]. V. Bapst, Journal of Mathematical Physics , 092101(2014), arXiv:1303.4908 [math.PR]. G. Parisi, S. Pascazio, F. Pietracaprina, V. Ros, andA. Scardicchio, Journal of Physics A: Mathematical andTheoretical , 014003 (2019), arXiv:1812.03531 [cond-mat.dis-nn]. M. Tarzia, Phys. Rev. B , 014208 (2020),arXiv:2003.11847 [cond-mat.stat-mech]. We have checked that varying c from 16 to 128 do notmodify the results. In fact, the diagonal elements of the resolvent can be ob-tained by matrix inversion, which is slightly faster thanED. B. Altshuler and V. Prigodin, Sov. Phys. JETP68