The manybody localized phase of the quantum random energy model
Abstract
The random energy model (REM) provides a solvable meanfield description of the equilibrium spin glass transition. Its quantum sibling (the QREM), obtained by adding a transverse field to the REM, has similar properties and shows a spin glass phase for sufficiently small transverse field and temperature. In a recent work, some of us have shown that the QREM further exhibits a manybody localization  delocalization (MBLD) transition when viewed as a closed quantum system, evolving according to the quantum dynamics. This phase encloses the familiar equilibrium spinglass phase. In this paper we study in detail the MBLD transition within the forwardscattering approximation and replica techniques. The predictions for the transition line are in good agreement with the exact diagonalization numerics. We also observe that the structure of the eigenstates at the MBLD critical point changes continuously with the energy density, raising the possibility of a family of critical theories for the MBLD transition.
Contents
I Introduction
Experiments on fewatom systems show striking contradictions with the extension of classical laws to arbitrarily small distances and energies. Quantum mechanics has proven to be a successful theory for the dynamics of these microscopic systems. Further developments make clear that quantum effects are also important for macroscopic systems, although typically at low temperature or high density.
On the other hand, at sufficiently high temperature or low density, the semiclassical principle appears to state that quantum and classical dynamics yield the same results for most measurable quantities. It is then interesting to examine the cases in which this reconciliation does not occur. Disordered systems make a particularly interesting playground. Noninteracting particles in a disordered potential can exhibit Anderson localization Anderson (1958): the gas’ diffusion coefficient(s) are completely suppressed, even in a regime where the classical dynamics remains diffusive. This difference is so spectacular that in low dimensions the dynamics at all energy scales supported by the system is localized. Here the classical description is never, not even qualitatively, accurate.
Whether Anderson localization survives the introduction of interactions between particles or not has been a subject of debate since the beginning of the field Mott (1969); Fleishman and Anderson (1980). Yet in the last decade, the work generated by the seminal paper of Basko, Aleiner and Altshuler Basko et al. (2006) on interacting, disordered systems has shed considerable light on the issue Oganesyan and Huse (2007); Žnidarič et al. (2008); Pal and Huse (2010).
It is now clear that a sufficiently small interaction is for many purposes irrelevant (one exception being entanglement Bardarson et al. (2012); Nanduri et al. (2014)). The system behaves as if the singleparticle occupation numbers are “perturbatively dressed” into the interacting phase Huse et al. (2014); Serbyn et al. (2013); Imbrie (2014); Ros et al. (2015); Chandran et al. (2015). Excitations in the “manybody localized” (MBL) regime are localized, transport is suppressed, and the dynamics is not thermalizing. For some spin chains, researchers have even been able to pinpoint a transition between an MBL region and an ergodic region Pal and Huse (2010); De Luca and Scardicchio (2013); Kjäll et al. (2014); Luitz et al. (2015); MondragonShem et al. (2015); Goold et al. (2015). The properties of this transition are not well understood. A few elementary constraints have been imposed Grover (2014); Chandran et al. (2015) and similarities with infinite randomness fixed points have been found in onedimensional studies Pal and Huse (2010); Vosk and Altman (2013); Potter et al. (2015). Like the Anderson transition, this is not a usual thermodynamic phase transition but rather a dynamical phase transition. There is no local order parameter in terms of which to write a LandauGinzburg free energy density. The transition itself is the breakdown of the hypothesis under which one derives the statistical description from the underlying microscopic equations of motion. Consequently, this transition can be observed even at infinite temperature.
An even more recent line of research studies how ergodicity breaking in the quantum dynamics compares to that of more canonical (classical and quantum) glassy phases. Since MBL is easily observed in spin chains with quenched disorder (and also a phase more akin to configurational glasses has been conjectured Schiulaz and Müller (2013); Pino et al. (2015); Yao et al. (2014); Papic et al. (2015)), it is natural to look for MBL in spin glasses.
With this in mind, some of us have recently looked Laumann et al. (2014) at a quantum version of Derrida’s random energy model (REM) Derrida (1981), which is a simplified model of meanfield spin glass. The quantum model’s equilibrium phase diagram has been studied before Goldschmidt (1990) and a glassy phase was identified at low temperature and small transverse magnetic field ^{1}^{1}1A similar phase diagram is found in more realistic, still meanfield quantum spin glasses Laumann et al. (2008). This should be generic for a large family of models, including combinatorial optimization problems in quantum annealers Laumann et al. (2015)..
In Laumann et al. (2014) it was found that the quantum dynamics is ergodic for high temperatures and large transverse field, but ergodicity breaks down upon lowering the temperature and the transverse field. The ergodicitybroken phase, which we identify as the MBL phase, encompasses the glassy region. Therefore the quantum dynamics becomes nonergodic before the glassy phase sets in, disentangling the concept of ergodicity breaking from that of replica symmetry breaking in these quantum models. A qualitatively similar observation was recently made by studying RokhsarKivelsontype wavefunctions derived from the REM Chen et al. (2015). However, on further thought, what is really surprising is that there is an ergodic phase for the simple transversefield quantum dynamics, since the classical Monte Carlo dynamics of the REM is never ergodic. This observation is relevant for the science of quantum annealers, for which the MBL phase could be a significant stumbling stone Altshuler et al. (2010); Laumann et al. (2015).
In this paper we analyze in much more detail the MBL phase of the QREM and the transition between the ergodic and MBL phases. We also study in detail the application of the forwardscattering approximation (FSA) in the MBL phase. We find that the localized phase is consistently distinct from the ergodic phase in its levelspacing statistics, observables, and eigenstate structure. Naive perturbation theory cannot accurately characterize the MBL phase, but by carefully handling neardegeneracies in the FSA, we quantitatively describe both the localized eigenstates and the phase boundary within perturbation theory. We accomplish this using a combination of simple approximations, numerics and a replica treatment of the forwardscattering wavefunctions.
The QREM is the first model in which a manybody mobility edge was clearly observed in the numerics accompanied by an analytical prediction within the forwardscattering approximation. It is an ideal test bed to discuss the properties of the mobility edge. One of the things that we observe numerically is that the critical statistics of the eigenvalues changes continuously along the mobility edge. It is worth mentioning the analogy with the critical properties of meanfield spin glasses Sompolinsky and Zippelius (1982); Parisi and Rizzo (2013), which also change continuously with lowering the temperature.
As a final remark, we notice that on lowering the transverse field the mobility edge shifts to higher temperatures. Furthermore, a mobility edge opens up even for an infinitesimal transverse field, in contrast to the situation in onedimensional systems, where the MBL phase at infinite temperature is stable upto a finite value of the interaction. This is most probably a special feature of the infinite dimensionality of the QREM.
Ii Phenomena
The quantum randomenergy model (QREM) for spins is defined by
(1) 
Here, is the transverse field, and is a random operator diagonal in the basis, with the diagonal entries identically and independently distributed according to
(2) 
With this normalization the spectrum of is with high probability contained in . We note that throughout the paper, capitalized represents extensive energy and the corresponding energy density.
ii.1 Equilibrium Phase Diagram
Goldschmidt Goldschmidt (1990) determined the canonical phase diagram of the QREM by using the SuzukiTrotter expansion and replica trick. The results we give here are taken from his paper, in which detailed derivations are found as well. See also Jörg et al. Jörg et al. (2008) for a simple perturbative derivation of some parts of the phase diagram.
The QREM has three equilibrium phases:

The REM paramagnet. The free energy density is
(3) This is equal to the free energy density of the classical REM at temperature above
(4) All thermodynamic quantities are identical to the zerofield REM, in particular the energy density and the entropy density .

The quantum paramagnet. The free energy density is
(5) Note that this is the free energy density for noninteracting spins in a field . The REM term in the Hamiltonian does not influence the equilibrium physics of this phase. Thermodynamic properties such as the energy, entropy, and magnetization density are given by the standard formulae.

The spinglass. The free energy density is simply
(6) identical to the classical REM. The system is frozen in its ground state at , where there are states and thus .
The transition lines between these phases are located where the free energies are equal. The boundary between the REM paramagnet and the spinglass is as in the classical model: the system is paramagnetic for and frozen in the spinglass phase for ( given by Eq. (4)). The system undergoes a firstorder transition to the quantum paramagnet when increases past , where is defined for by and for by . The blue dashed lines of Fig. 1b mark these boundaries. In particular, note that at is the field strength at which the ground state energy of the quantum paramagnet matches the ground state of the classical REM.
The essence of the canonical phase diagram is that the system makes a firstorder transition between “REM” physics and “quantum paramagnet” physics. REM physics results from an otherwise structureless system freezing into an intensive number of configurations at nonzero temperature. Quantum paramagnet physics is that of noninteracting spins in a magnetic field. The QREM does not compromise between these two regimes. It exhibits only one or the other, at least thermodynamically.
ii.2 Dynamical Phase Diagram
When treated as an isolated system, the QREM has two dynamical phases:

The ergodic phase. Eigenstates satisfy the Eigenstate Thermalization Hypothesis (ETH) Srednicki (1994); Rigol et al. (2012): expectation values of local observables agree with the microcanonical ensemble, which, here, is paramagnetic. Thus for all spins . Fluctuations in are large but decay exponentially in time. In addition, these eigenstates are delocalized in the following sense: one can map a configuration of spin’s to a corner of an dimensional hypercube by considering as the top (bottom) face of the cube’s ’th dimension. The QREM Hamiltonian is then an Anderson model on the corners of this hypercube, with the spin configurations being “lattice sites”. Ergodic eigenstates are delocalized over these sites. The probability overlap between two ergodic states decays exponentially with their energy difference, as observed in the Anderson model Evers and Mirlin (2008).

The manybodylocalized phase. Eigenstates are weaklydressed single configurations of spins. Thus . Fluctuations within an eigenstate are small (and decrease with system size), but fluctuations between eigenstates over energy and realization of disorder are order. These eigenstates are Andersonlocalized on the hypercube. The probability overlap between a pair of eigenstates decays as , where is the energy difference of the two states. This is welldescribed by firstorder perturbation theory.
A major focus of the present work is to quantitatively locate and characterize the boundary between these phases in the  plane ( is eigenstate energy density). We define order parameters with welldefined localized and ergodic limits for each characteristic described above. See Figs. 2 and 3 for examples and Sec. III for details. Since the states at are statistically equivalent to those at , the phase boundary must be symmetric about . We focus on the negative portion, denoted . Numerical evidence and perturbation theory each give constraints on . We show these, overlaid on our conjectured phase diagram, in Fig. 1.
Our numerical results indicate a transition between two phases and provide a lower bound for . The dark green markers and shading in Fig. 1 indicate where order parameters computed via exact diagonalization (Sec. III) transition between the two limits. Similarly, the light green markers and shading are the corresponding regions from a numerical perturbation series in (Sec. V.3). The transitions sharpen as increases, and a finitesize scaling analysis indicates that they become sharp in the thermodynamic limit. The energy density around which the transitions sharpen is consistent within numerical error amongst all order parameters. Thus we identify this energy density, once properly extrapolated to infinite , as . It is difficult to make the extrapolation from accessible system sizes, and we observe a slight finitesize drift in the transition towards smaller as increases. For this reason, the marked points in Fig. 1 are actually lower bounds for .
Expanding the eigenstates perturbatively in sets upper bounds. A rough analytical treatment of the perturbation series sets a clear upper bound for at (Sec. V.1). Proceeding more carefully but still analytically (Sec. V.2), we obtain the green dashed curve in Fig. 1. This curve is still only an upper bound on , but it behaves as at small .
lies between the green markers and the green dashed curve in Fig. 1. In addition, we conjecture that the bound as is tight, and conjecture that
(7) 
for finite as well. This is the black dashed line in Fig. 1. Jörg et al. observed Jörg et al. (2008) that the QREM ground state undergoes a quantum phase transition between the REM ground state and at ( given by Eq. (4)). Our conjecture for implies that the dynamical phase boundary is consistent with this result. It is consistent with all of our numerical and analytical bounds as well.
The critical eigenstates, i.e., those with , have order parameter values intermediate between the localized and ergodic limits. We are able to fit the dimensionless order parameters to a finitesize scaling ansatz ( is the order parameter). The critical exponent , independent of . The critical amplitudes , though, do depend on . This could be due to finitesize effects, but if not, the critical eigenstates are described by a line of fixed points.
Finally, it is interesting to compare the QREM’s manybody mobility edge to that of a transition between ensembles of random matrix theory (RMT). In particular, the random symmetric matrix , where is a diagonal random matrix whose elements are i.i.d. random variables of order 1 and is a Gaussian random matrix, has been studied extensively as a model for the crossover from Poisson to GOE statistics. The extent to which exhibits Poisson versus GOE statistics is determined by the ratio of ’s mean level spacing, denoted , to the typical offdiagonal matrix elements. If has eigenvalues distributed uniformly over a fixed interval, for matrix dimension and thus the scale on which ’s level statistics transition from Poisson to GOE is .
We, however, are interested in when ’s eigenvalues are distributed according to Eq. (2). The mean level spacing is welldefined locally, but it depends on the local energy density . In particular,
(8) 
By setting , we see a smooth interpolation from Poisson to GOE statistics at the energy density . Yet for this choice of scaling for , all other energy densities flow exponentially fast towards one of the limiting statistics as increases. The phase boundary is horizontal in the plane. This is in marked contrast to the QREM, where we find a nontrivial . The locality of the transverse field in Eq. (1) thus plays an important role in the physics of this transition.
Iii Exact Diagonalization
The most direct way to study the QREM is through exact diagonalization. We generate between and realizations of the Hamiltonian in Eq. (1) for system sizes ranging from to and ranging from to . We obtain the entire spectrum and complete set of eigenstates through full exact diagonalization. The eigenstates are then binned according to energy density so as to study how average properties depend on energy.
We use three types of averages, often simultaneously. Quantummechanical averages within single eigenstates are denoted by angular brackets , averages among eigenstates of a single sample within an energy density window are denoted by a bar , and averages between realizations of disorder are denoted by square brackets . Unless noted otherwise, every quantity that we describe in this section is averaged both within an energy density window and over disorder.
iii.1 Spectral Statistics
The most straightforward numerical approach to the localizationdelocalization transition is to study the statistical properties of the spectra. Each spectrum can be split into two regions: WignerDyson statistics dominate the distribution of energy levels for close to , whereas Poisson statistics dominate for far from . From a randommatrixtheory perspective, WignerDyson statistics govern the spectra of matrices in which every element is an independent random variable (strictly speaking, the WignerDyson distribution gives the probability of level spacings in GOEdistributed matrices). Poisson statistics, on the other hand, describe the energy gaps when the eigenvalues themselves are independent and uniformly distributed. As is commonly done, we identify WignerDyson statistics with delocalized eigenstates and Poisson statistics with localized eigenstates. We find that the transition between these two is sharp in the thermodynamic limit. Thus we identify the large region as a localized phase and the small region as a delocalized phase.
To be more quantitative, we use three measures of the energy gap distribution. First, we calculate the level spacing ratio Oganesyan and Huse (2007). We define
(9) 
and consider the average, . represents the degree of level repulsion in the system: there is less variation in the separation between energy levels, and thus a larger , when those levels repel each other. The Poisson distribution, which corresponds to independent energy levels, has . Compare this to the WignerDyson distribution, for which level repulsion is significant and . The values of that we obtain from (1) lie between these two limits and tell us the strength of level repulsion at the target energies. We in turn interpret this as the degree of delocalization.
In addition to , we determine the cumulative distribution function of level spacings . We characterize the CDF by two quantities used in the literature Zharekeshev and Kramer (1995): , i.e, the fraction of spacings that are less than (roughly) half the mean, and , i.e, the fraction that are greater than twice the mean. Both these values are much smaller in the WignerDyson distribution than in the Poisson distribution, since level repulsion suppresses the appearance both of packed and isolated eigenvalues. Thus and , like , quantify the degree of level repulsion, by examining the frequency of small gaps and by examining the frequency of large gaps. We present each as a deviation from WignerDyson statistics relative to that of the Poisson distribution, i.e., we report
(10) 
and
(11) 
where the superscript refers to WignerDyson statistics and the superscript refers to Poisson statistics.
Fig. 4 shows , , and at various as a function of energy density (at , a representative field strength), successfully collapsed using the scaling form . From these we obtain and , which are shown as functions of in Fig. 5a,b. Except for at , all three statistics predict a common and independent . Fig. 5c shows the corresponding level statistic values at , which we obtain from the collapsed curves. The statistics become more WignerDysonlike as increases. It is in this sense that the character of the transition depends on .
iii.2 Local Eigenstate Observables
The magnetization of single spins in each phase clearly distinguishes the two. Since a delocalized eigenstate has and a localized eigenstate has , we consider . Fig .6 shows the values as a function of . We see clear finitesize flow towards in the delocalized phase and towards in the localized phase. The crossover region’s location is consistent with the spectral statistics.
The limit in the localized phase violates the ETH. We further demonstrate that the localized eigenstates fail to thermalize by studying how fluctuates from state to state within a sample. More precisely, we define
(12) 
Fig. 7 shows the distribution of , for all within an energy density window, as a function of . In the localized phase, the distributions have weight at . The total weight at is as much as at , signifying that is as likely to switch sign from one eigenstate to the next as it is to not. In the delocalized phase, the entire weight is centered around . This is further evidence that each individual delocalized eigenstate is thermal.
iii.3 Connected Autocorrelations
The operators evolve over time as a result of the transverse field, so correlations in time are important characteristics of the two phases. We quantify this by studying the connected autocorrelation function
(13) 
This quantity is computed exactly from the full exact diagonalization, up to in steps of . See Fig. 3 for the results. is sufficiently close to 0 at all times in the localized phase that it is hard to extract a meaningful decay. This is consistent with the localized eigenstates being weakly dressed single configurations of spins with magnetization .
In the delocalized phase, is initially close to 1 and decays exponentially to 0. The decay becomes slower as nears the transition. We fit the longtime behavior of to a straight line and extract the slope. From this we study how the decay time depends on . See Fig. 8. At , the decay time saturates exponentially with to a finite value (not shown). As , increases monotonically. Very close to , is diverging with system size. We obtain good scaling collapse with the form . This form captures ’s linear dependence on near the transition. We again find that and agree with the spectral statistics.
iii.4 Eigenstate Structure
As mentioned in Section II.2, configurations of spin1/2s map to the corners of an dimensional hypercube. The basis is then the “coordinate” basis, and we consider diagnostics from Anderson localization Evers and Mirlin (2008). These have already been used in the MBL context Buccheri et al. (2011); De Luca and Scardicchio (2013), particularly the inverse participation ratio (IPR), defined as
(14) 
The IPR quantifies how many sites of the hypercube is distributed over. corresponds to concentrated on a single site and corresponds to with equal weight on all sites. Fig. 9 shows , a measure of ’s typical order of magnitude. flows towards 0 with system size in the localized phase, signifying that the eigenstates are concentrated onto single sites. decreases linearly with in the delocalized phase, signifying that the eigenstates extend over Hamming distances of order . As with , there is a welldefined crossover region consistent with the spectral statistics.
We also compute the eigenstate correlation function , defined as
(15) 
where the sum over and is over eigenstates. This quantity measures the overlap between eigenstates’ probabilities as a function of their location in the spectrum and energy separation. is known as the twoparticle spectral function in the context of finitedimensional lattices Chalker and Daniell (1988). There, decays exponentially with for and both in the delocalized phase. The scale on which decays is the Thouless energy ( is the diffusion coefficient and is the linear system size). Fig. 10 shows the same behavior in the QREM’s delocalized phase. The characteristic scale is independent of , which is reasonable since our fullyconnected model has an effective “linear size” of 1. also decays exponentially with , for the same reason that does. In the localized phase, decays as . We justify this in Section IV, where we apply firstorder perturbation theory to this model. The results in Fig. 11 agree remarkably well with Eq. (23) below.
Iv Naive Perturbation Theory
The QREM eigenstates differ drastically across a welldefined boundary. The spectral statistics and eigenstate properties characterize and distinguish the two phases. However, the numerical results by themselves are somewhat opaque and limited to finite sizes. We get a more transparent picture of the MBLD transition by considering the perturbative structure. Section V systematically investigates the perturbativein expansion of the eigenstates. Here we study the firstorder corrections heuristically and compare to the exact diagonalization results for the localized phase.
To first order in , each QREM eigenstate has probability amplitude 1 on its initial site of the hypercube (with unperturbed energy ) and probability amplitude on the ’th adjacent corner. All other sites have probability amplitude 0. For the IPR and of such an eigenstate, we obtain:
(16) 
(17) 
We now study the distributions of Eq. (16) and Eq. (17) over samples, treating as a fixed parameter of order . With probability 1 (in the thermodynamic limit), the random variables are all and we expand in powers of . This is not strictly allowed because, although its typical values are small, the moments of diverge. For now, we assume that the typical values of Eq. (16) and Eq. (17) will be satisfactory estimates and so assume that all .
To determine the disorder averages, we need only keep the first term that is even in all . To determine the disorder variances, we in addition require the first term odd in . Carrying this out,
(18) 
(19) 
where , . All have variance , so all have variance and has variance . Thus
(20) 
(21) 
In addition, note that at this order.
We also estimate in this manner. Suppose that two eigenstates have energies and . In order for their wavefunctions to overlap at first order in , their zerothorder eigenstates must be on adjacent corners of the hypercube. This occurs in a fraction of the samples. When the eigenstates do overlap at first order,
(22) 
Thus
(23) 
These are simple estimates, and we can generalize them to higher orders. However, most turn out to be incompatible with our exact diagonalization results, even in how they scale with and . The one exception is , for which Eq. (23) is very accurate in the welllocalized region. We expect that the overall discrepancy is due to our limited range of system sizes. In particular, we have expanded in , which is in typical samples. Yet even for our largest system (), . Similarly, the exponentially small number of “atypical” samples is not that small. Eq. (23) successfully describes because there is no need to control the random variables at first order, since by definition must be equal to . When we remove complications arising from the randomness in , the leadingorder estimates become highly accurate, but sadly we cannot do this for any but .
With our limited systems sizes we must turn to more sophisticated approximation schemes to understand the localizationdelocalization transition. We have two goals in doing so: to quantitatively predict the critical energy density in a procedure more efficient than exact diagonalization, and to understand how the structure of localized eigenstates changes as approaches . The ForwardScattering Approximation accomplishes both.
V Forward Scattering Analysis
The ForwardScattering Approximation (FSA) treats the transverse field perturbatively and keeps only the leadingorder contribution to the wavefunction amplitude at each site (see Huse et al. (2015); Pietracaprina et al. (2015) for recent applications of the method to the Anderson problem). We approach this approximation from the hypercube perspective. The term in the Hamiltonian is a random onsite potential and the term allows hopping between adjacent sites. Without loss of generality, our unperturbed state has all . Denote it by , with unperturbed energy . The perturbed state is denoted with energy . From timeindependent perturbation theory,
(24) 
where . The leadingorder contributions to are the “directed” paths from to , i.e.,
(25) 
where the sum is over the directed paths from the initial site to the site and is the unperturbed energy of the th site along a given path . We treat as a tunable parameter. All are random variables distributed according to Eq. (2).
Eq. (25) is the forwardscattering approximation to the eigenstates of the REM in a transverse field. We use it to investigate signatures of delocalization. As long as the perturbation theory is valid, i.e., all , localization persists. Strictly speaking, the appearance of order1 , which we refer to as “resonances”, does not mean that the exact eigenstates are delocalized, only that Eq. (25) cannot describe them. Regardless, we take the energy density at which resonances proliferate in the perturbed wavefunction as a strong lower bound for .
Although the FSA is an approximation, as it stands, Eq. (25) must still be evaluated numerically. We do so below, but first consider further approximations. In Subsection V.1 we make quick estimates as to where the perturbation series breaks down. We refine these estimates in Subsection V.2. In Subsection V.3 we evaluate Eq. (25) numerically. Finally, in Sec. V.4 we briefly consider the replica treatment of the FSA, which reproduces many of the results obtained more directly in the previous sections.
v.1 Rough Estimates
From Eq. (2), most configurations have energies of order . Thus most terms in Eq. (25) are to leading order in , and a rough estimate comes from replacing all terms with . Then
(26) 
As a function of , first exceeds at , if at all. We estimate by setting :
(27) 
We have written the estimate above as a bound on the position of the transition curve because it clearly underestimates the possibility of smalldenominator resonances. Even though most REM energies are of order , an exponentiallylarge number of states are still within any finite window near . Such states contribute much more than to Eq. (25). Nevertheless, this estimate agrees with ED in predicting that the transition lies at finite energy density/temperature. It also suggests that as . Localization persists up to higher temperatures as the transverse field weakens.
This analysis is clearly invalid at , where the weights are large for most sites. In this regime, the probability distributions for path contributions (the terms in (25)) are longtailed. Thus, we expect the largestweighted path to dominate the sum over them. The “greedy algorithm” estimates this largestweighted path as follows: most sites have REM energies of order , and the smallest energy among sites will be . The smallest energy of sites neighboring the initial site is , as the initial site has neighbors. The corresponding site has a neighbor at distance 2 with energy , for whom the smallest energy of neighboring sites at distance 3 is , etc., such that the largest path is typically
(28) 
At large ,
(29) 
The wavefunction amplitude decays exponentially with distance for , where
(30) 
Any nonzero independent field causes the perturbation series to break down. The crossover region goes as , and once exceeds , the breakdown occurs immediately, i.e., at .
v.2 Single Resonance Approximation
Consider the probability that a site at distance produces , given that the closer sites all have energies of order . This checks the consistency of the treatment above, and we expect single sites to produce resonances in this way because the distribution of is longtailed. From Eq. (26), all sites at a distance have amplitude
(31) 
Then the site at distance has amplitude
(32) 
The probability of is
(33) 
Since the are independent, the for all sites at distance are (under the current assumptions) independent, and thus the probability of no resonances at distance is
(34) 
where . In the large limit, this is with
(35) 
and . There will not be resonances in the thermodynamic limit and localization will persist if for all . is maximized at with a value of
(36) 
The zero of Eq. (36) is at . We find that
(37) 
The phase diagram in Fig. 1 displays this estimate as a function of . We get a much larger estimate than that in Eq. (27), yet is still proportional to as .
In addition to an estimate for the critical energy density, this shows how suddenly and drastically the perturbation theory breaks down. Suppose we had defined a resonant site to be one whose amplitude exceeds . Then the probability of no resonances at distance would be . The large asymptotics only change if grows/decays faster than . Thus the resonances that appear at are exponentially large in system size, and all sites have exponentially small amplitudes at .
We similarly study the infinitetemperature case, where the greedy algorithm estimated . As described above, we cannot consider single sites being resonant but we can consider single paths being resonant. A path to a site at distance has amplitude
(38) 
To simplify notation, we write this as
(39) 
where
(40) 
(41) 
For , is distributed as
(42) 
Thus , i.e., , occurs with probability
(43) 
We use this estimate to determine the probability that none of the paths to sites at distance give resonances. The paths are not independent, but we assume that we can treat them independently in estimating whether or . Then the probability of no resonances at distance is
(44) 
Consider the exponent , and in particular the ratio
(45) 
The righthand side of Eq. (45) is decreasing as a function of . When it drops below 1, the exponents in Eq. (44)) decrease as increases, i.e., resonances become less likely as we move farther into the hypercube. In order to declare that resonances are unlikely throughout the entire hypercube, we require that . This is only possible when is such that
(46) 
The critical is, to leading order,
(47) 
decays to 0 faster than the greedy estimate in Eq. (30) by a factor of , but otherwise Eq. (47) confirms the qualitative description. Perturbation theory inevitably breaks down in the thermodynamic, hightemperature limit at any nonzero . We see the resonances immediately, i.e., at an order in the expansion. This contrasts with the finitetemperature case, where the perturbative description persists until finite and the resonances appear after terms in the series.
We have a description of when and how QREM eigenstates become nonperturbative in , but we have made some uncontrolled approximations. Order wavefunction amplitudes need not be due solely to single resonant sites or paths. This then implies that the amplitudes are not independent of each other. Below, we study the statistics of the wavefunctions in Eq. (25) numerically, without any additional approximations. We check not only the estimate for the critical energy density but also whether the nature of the eigenstates agrees with the description above.
v.3 Numerical Treatment
To quantify resonance proliferation, we generate sets of REM energies and evaluate Eq. (25). First, we directly compute what we attempted to estimate in Sec. V.2: the probability of a sample at fixed and containing at least one resonance. In addition, we calculate the sampleaveraged IPR and singlespin magnetization, for comparison with the exact diagonalization. These statistics each provide an independent measure of when resonances show up and how significant the resonances are. See Fig. 12 for scaled results, and refer back to Section III for definitions and notation.
To begin, we count the fraction of samples that contain a resonance at each value of and . There is a crossover from at large to at small (Fig. 12a). This crossover sharpens as increases, consistent with a zeroone law in the thermodynamic limit. We estimate the critical by finitesize scaling with the form . The curves collapse well, and is consistent with the exact diagonalization results. , however, is significantly larger than ED predicts. The curve in Fig. 12a has .
We next compute the disorderaveraged IPR (Fig. 12b). For , increases towards as increases. This is consistent with localization and with the exact diagonalization results. It also confirms that only a negligible fraction of samples in this phase contain resonances. For , however, is roughly independent of both and , and stays at . This signifies resonances, and presumably delocalization, although the signature is very different from that of the exact eigenstates. The probability distributions for the wavefunction amplitudes are very longtailed. Thus in a given sample, a few amplitudes will be much larger than the rest. In the localized phase, these largest amplitudes still do not compare to those of the initial site. As approaches , those large amplitudes approach and decreases. Once passes , the initial site no longer contributes to the IPR and is set by the largest noninitial amplitudes. Increasing (and ) even further will make the wavefunction amplitudes larger, but it won’t change that only a few sites dominate the IPR. Hence is independent of both and once the perturbation theory breaks down.
Keep in mind that the FSA’s prediction for in the delocalized phase has no connection to the actual IPRs in the QREM’s delocalized phase. The FSA is a perturbative expansion, and all we learn from the initial site no longer contributing to is that we cannot trust any results from the FSA in this region. Fig. 13 plots the FSA’s disorderaveraged IPR against the disorderaveraged IPR from exact diagonalization. The two curves agree for all system sizes when . The deviations become significant at , though, and the curves bear no relation to each other for . The FSA successfully describes the structure of eigenstates in the localized phase but not in the delocalized phase. Yet even though the FSA IPR is not accurate for much of the spectrum, a scaling collapse gives critical parameters that agree with those of the resonantsample fraction .
The magnetization of a single spin behaves analogously (Fig. 12c). Recall that we arbitrarily take the unperturbed state to have all , and sites closer to the initial site have more spins pointing down. Thus is roughly constant at for large , and it starts to increase significantly once passes . As before, the FSA and ED predictions for agree in the localized phase but not in the delocalized phase. The curves collapse onto each other quite well, with scaling parameters that again agree with the others from the FSA.
v.4 Replica Treatment
The statistical properties of the forward scattering wavefunctions can also be studied using the replica method, which provides complementary understanding to the other numerical and analytical approaches Mézard et al. (1987); Mézard and Montanari (2009). In this approach, we view the amplitude as the partition sum of a directed random polymer living on the hypercube with the longtailed random weights . As these weights do not have any finite moments, we expect the directed random polymer to condense onto a small number of large weight paths Kardar (1994). This also justifies our neglect of the sign of the weights as their destructive interference is unimportant in this regime. The replica approach is especially useful as it naturally regulates the divergence of these moments.
Within the forward scattering approximation, the wavefunction amplitude at distance is
(48) 
where is the random weight of the piece of path going through site on path .
The typical value of is provided by averaging its logarithm, , which may be calculated using the formal replica trick
(49) 
Thus, we need to calculate
(50) 
where gives the number of times that the paths cross site :
(51) 
Replica symmetric ansatz—
The replica symmetric ansatz consists of assuming that each of the paths contributes independently to the ’th moment of . That is, for each of the sites visited by the paths and otherwise. Thus,
(52)  
(53) 
Thus, as , we find
(54) 
This is illdefined if , which is clearly true for the weights arising in the QREM. However, for , the weight in the tail of is parametrically suppressed by . If we simply replace by its typical value , we find the ‘typical weight RS’ result,
(55) 
where . This indicates that the typical amplitude decays exponentially with at all points inside the hypercube for .
This clearly overestimates the critical for delocalization as it neglects the possibility of small denominators and the concomittant atypical resonances. Nonetheless, even at this level it shows that delocalization should take place for infinitesimal at . This estimate agrees precisely with the rough estimate made directly in Sec. V.1.
Replica symmetry breaking ansatz—
In the 1RSB ansatz, the dominant configurations contributing to consist of tightly bound groups of paths each. Thus,