Perturbative instability towards delocalization at phase transitions between MBL phases
PPerturbative instability towards delocalizationat phase transitions between MBL phases
Sanjay Moudgalya, David A. Huse, and Vedika Khemani Department of Physics, Princeton University, Princeton, NJ 08544, USA Department of Physics, Stanford University, Stanford, CA 94305, USA (Dated: August 21, 2020)We examine the stability of marginally Anderson localized phase transitions between localizedphases to the addition of many-body interactions, focusing in particular on the spin-glass to param-agnet transition in a disordered transverse field Ising model in one dimension. We find evidence fora perturbative instability of localization at finite energy densities once interactions are added, i.e. evidence for the relevance of interactions - in a renormalization group sense - to the non-interactingcritical point governed by infinite randomness scaling. We introduce a novel diagnostic, the “suscep-tibility of entanglement”, which allows us to perturbatively probe the effect of adding interactionson the entanglement properties of eigenstates, and helps us elucidate the resonant processes that cancause thermalization. The susceptibility serves as a much more sensitive probe, and its divergencecan detect the perturbative beginnings of an incipient instability even in regimes and system sizesfor which conventional diagnostics point towards localization. We expect this new measure to be ofindependent interest for analyzing the stability of localization in a variety of different settings.
I. INTRODUCTION
Many-body localization (MBL) was born in investiga-tions of the stability of Anderson localization – the phe-nomenon that strong enough disorder exponentially lo-calizes non-interacting wavefunctions – to the additionof interactions [1, 2]. This stability was demonstrated toall orders in perturbation theory, following early precur-sors [3–5]. While the non -perturbative stability of MBLremains an open question in various settings [6–9], it hasbeen proven that a stable MBL phase can exist in lo-cal, strongly-disordered one dimensional spin chains [10].More generally, understanding the stability of phenom-ena to small changes, such as the introduction of interac-tions, is a central enterprise in theoretical physics. In thiswork, we add to this important literature by examiningthe effect of interactions on marginally Anderson local-ized critical points between Anderson localized phases.While phases and phase transitions are traditionallystudied in the framework of equilibrium statistical me-chanics, recent work has shown that there is a rich no-tion of phase structure even within the out-of-equilibriumMBL phase [11, 12]. Different MBL phases representdistinct types of novel dynamical phenomena that maybe completely invisible to, or forbidden by, equilibriumthermodynamics — a paradigmatic example being the re-cently discovered Floquet MBL time-crystal phase [13–15]. Localized phases can be understood as eigenstatephases characterized by distinct patterns of long-rangeorder (LRO), both symmetry-breaking and topological,in individual highly-excited MBL eigenstates [11, 12, 16–19]; a phase transition between different localized phasesrequires singular changes in the eigenspectrum proper-ties. Indeed, the passage from localization to thermal-ization is itself a dynamical phase transition involving asingular change in the entanglement properties of highlyexcited many-body eigenstates. While the nature of the MBL-to-thermal phase transition has been a subject ofintense study [20–33], transitions between different MBLphases have recently received considerably less attentionand are the focus of this work.While our conclusions are quite general, for specificity,the majority of our analysis will be presented for a disor-dered transverse-field Ising model (TFIM) in one dimen-sion. This model exhibits a localized, symmetry-broken“spin-glass” (SG) phase with LRO, and a localized para-magnetic (PM) phase with no order [11, 12, 22, 34, 35].In the absence of interactions, the phase transition be-tween the SG and PM phases is governed by an infinite-randomness fixed point that is studied using the strongdisorder renormalization group (SDRG) [12, 35]. Whileall single-particle (SP) eigenstates are exponentially lo-calized in either phase, the critical point (CP) is only marginally localized. The SP eigenstates correspondingto SP energies
E → stretched exponentiallylocalized at the CP, and both the SP density of statesand localization length diverge in this limit [35].Once weak interactions are added, localization remainsstable deep in the PM and SG phases, on account of theusual arguments for the stability of Anderson localiza-tion with strong enough disorder. The stability of local-ization near the CP, however, requires careful consider-ation. The CP exhibits “marginal” Anderson localiza-tion due to the presence of (weakly) extended states inthe SP spectrum, which could aid in the formation oflong-range “resonances” that make the CP more suscep-tible to thermalization [6, 36]. We note that the instabil-ity towards thermalization may be visible in a perturba-tive treatment of the interactions [5, 36], or it may havea subtler non-perturbative origin, in which case it willonly be detectable at asymptotically large system sizesand times [6]. Previous SDRG studies of the interactingIsing model explicitly treat interactions as irrelevant –even at finite many-body energy densities – and so donot consider the resonances that could destabilize the a r X i v : . [ c ond - m a t . d i s - nn ] A ug CP [12, 37]. The stability of this transition i.e. the rel-evance of interactions to the non-interacting CP in RGlanguage, has been a major outstanding question in theliterature and is the subject of this work. We find evi-dence in favor of a perturbative instability of the CP tothe addition to arbitrarily weak interactions.This paper is organized as follows. In Sec. II, we studya self-dual TFIM using various diagnostics accessible tomany-body exact diagonalization (ED). We find thateven for the small system sizes accessible to ED, the CPthermalizes at a tiny interaction strength, λ c ∼ / λ c , together with the incipient thermalization alreadynoticeable at small sizes, strongly suggests a perturba-tive instability of the CP. We study this further usinga novel measure introduced in Sec. III, the susceptibilityof entanglement χ S , which serves as a sensitive probeof the effect of interactions on the entanglement proper-ties of many-body eigenstates. We derive a perturbativeexpansion for χ S in terms of the non-interacting eigen-states (Sec. IV), and use this to elucidate the low-orderprocesses that could destabilize localization at the CP(Sec. V). We expect χ S to be of independent interest inexamining the perturbative effect of interactions in dif-ferent contexts. II. MODEL AND DIAGNOSTICS
We study a disordered statistically self-dual TFIM ina one dimensional spin-1/2 system of length L with openboundary conditions: H = (cid:80) i J i σ xi σ xi +1 + (cid:80) i h i σ zi + λ (cid:80) i ( h ∗ σ zi σ zi +1 + J ∗ σ xi σ xi +2 ) ≡ H + λV. (1)The couplings { J i } ∈ [0 , J ∗ ] and fields { h i } ∈ [0 , h ∗ ] aredrawn independently and randomly from uniform distri-butions. When λ = 0, the model maps on to Andersonlocalized non-interacting Majorana fermions via a JordanWigner transformation, see App. A for a review. We ap-propriately scale terms in the interaction by the strengthsof the non-interacting couplings J ∗ and h ∗ , which al-lows us to study a dimensionless parameter λ that setsthe interaction strength while preserving statistical self-duality. As J ∗ and h ∗ are tuned to sweep across thephase diagram, we pick max[ J ∗ , h ∗ ] = 1, which simplysets an overall scale in the Hamiltonian (1) because of themanner in which the interaction strength is scaled. Wenote that localization is considered truly unstable onlyif thermalization happens for any strength or configura-tion of the disordered couplings, and we indeed observethe same qualitative behavior for the analysis presentedin the subsequent sections even for “stronger” power-law − . − . . . . Jh = log J ∗ − log h ∗ . . . . . . λ MBL spin-glassMBL paramagnet Thermal h r ih s ihGi FIG. 1. Finite size phase boundaries of the disordered self-dual Ising model (1) in the vicinity of the CP, numericallyobtained using three diagnostics discussed in Sec. II. If thereis a direct MBL PM-to-SG transition for a non-zero inter-action strength λ >
0, it will be at the self-dual J ∗ = h ∗ line (red). However, the most sensitive diagnostic, (cid:104)G(cid:105) , al-ready shows the onset of thermalization at a small value ofthe dimensional interaction strength, λ c ∼
1% at the self-dual J ∗ = h ∗ line. These phase boundaries are estimated using ex-act diagonalization on systems of size L ≤
14. We also findthat the boundaries strongly drift towards towards smaller λ as L is increased, suggesting an instability of the CP for anynon-zero λ in the infinite size limit. The white shaded boxrepresents this possibility i.e. that the thermal phase extendsall the way to λ c = 0, which is suggested by our data andanalysis but is inaccessible to finite-size numerics. disorder. .The Hamiltonian in Eq. (1) has Z Ising symmetry P = (cid:81) i σ zi , [ H, P ] = 0. Deep in the SG phase, J ∗ (cid:29) h ∗ , all many-body (MB) eigenstates look like pairs of Isingsymmetric superposition “cat” states with P = ± x axis: | n ± (cid:105) ≈ √ ( | ←→→ · · · ←(cid:105) ± | →←← · · · →(cid:105) ). Domainwall excitations between oppositely oriented spins are lo-calized, and the system shows long-range spin glass or-der: connected correlation functions of σ x evaluated in | n ± (cid:105) are non-zero for arbitrarily distant correlators, butwith random sign [11, 22]. By contrast, without disor-der, eigenstates at any finite temperature have a non-zero density of delocalized domain walls which destroy However, power-law disorder distributions with large exponentsare plagued with strong finite size effects since they have along tail to very small values; this leads to large separations ofscale between local couplings, effectively “cutting” the chain intosmaller pieces in finite-size numerics. Nevertheless, while morechallenging to demonstrate numerically, we expect the qualita-tive arguments and resonant processes identified in this work toalso hold for power law disorder distributions that have beenconsidered in various works [38, 39] − − λ . . h r i (a) L = 8L = 10L = 12L = 14PoissonGOE L − h G i (b) λ = 0.05 λ = 0.03 λ = 0.02 λ = 0.01 λ = 0.005 λ = 0.0025 FIG. 2. Diagnostics of thermalization with increasing inter-actions for ∆ Jh = 0. (a) Average (cid:104) r (cid:105) plotted against λ . Notethe strong drift of the finite-size crossings towards smaller λ with increasing system size. (b) (cid:104)G(cid:105) as a function of L for dif-ferent λ . Averages are over the middle half of all eigenstatesin 3000 to 300 independent disorder samples for sizes rangingfrom L = 8 to L = 14. LRO in one dimension. In the opposite limit in the PM, h ∗ (cid:29) J ∗ , the eigenstates resemble random product statesin the z basis. Far enough from the CP, the localized non-interacting SG and PM phases are stable to weak inter-actions, thermalizing only at a strong enough interactionstrengths λ = λ c > Jh ≡ log J ∗ − log h ∗ = 0 due to the Kramers-Wannier duality which maps the PM and SG phases toeach other. This transition takes place across the entire MB spectrum, and displays the same infinite random-ness scaling at all energy densities. It has been arguedthat interactions are irrelevant to the ground-state tran-sition (where there are no MB resonances) [35]; we areinstead concerned about the fate of the CP to the addi-tion of interactions at high energy densities. We studythe phase boundary to thermalization for different ∆ Jh : λ c (∆ Jh ), focusing on ∆ Jh = 0. The question is whether λ c (∆ Jh = 0) = 0, so that infinitesimal interactions ther-malize the non-interacting CP in the limit of an infinitesystem. We now turn to a numerical exploration of thephase diagram of this model as a function of ∆ Jh and λ .We start by probing the level repulsion in the eigenen-ergies { E n } within one Z symmetry sector P = 1, viathe ratio r = min( δ n ,δ n +1 )max( δ n ,δ n +1 ) where δ n = E n +1 − E n [40].In Fig. 2a, we plot (cid:104) r (cid:105) vs. λ for ∆ Jh = 0. Theaverage is taken over independent disorder realizationsand the middle half of the many-body eigenstates cen-tered around the energy density corresponding to infinite-temperature. The ratio changes from the localized Pois-son value, (cid:104) r (cid:105) ∼ = 0 .
39, to the thermal value, (cid:104) r (cid:105) ∼ = 0 .
53, as λ is increased [41]. The finite-size crossover to thermal-ization takes place at the interaction strength λ c ( L , L )where the curves for L , cross. Note that these crossingsare happening at rather small value of the dimensionlessinteraction, λ c (12 , ∼ = 0 .
05 even for the small sizes un-der study. While the sizes are too small to estimate the − − Jh h S i (a) L = 8L = 10L = 12L = 14 − − λ . . . . h s i (b) L = 8L = 10L = 12L = 14
FIG. 3. Entanglement diagnostics (a) Mean value of the en-tanglement entropy across a horizontal cut on the phase dia-gram of Fig. 1 for λ = 0 .
5. The entanglement entropy clearlychanges from 0 in the PM phase to a volume-law in the ther-mal phase to log 2 in the spin-glass phase. (b) Mean value ofthe entanglement entropy density s = S/S max across a ver-tical cut on the phase diagram of Fig. 1 for ∆ Jh = 0. Theentropy density decreases to zero in the critical phase andincreases towards 1 in the thermal phase. Averages are overthe middle half of all eigenstates in 3000 to 300 independentdisorder samples for sizes ranging from L = 8 to L = 14. asymptotic critical interaction strength λ c at large L , thelocations of the finite-size crossings strongly drift towardssmaller λ c with increasing L , consistent with the possi-bility that infinitesimal interactions destabilize the CP.The finite-size λ c ( L , L ) from the crossings between thetwo largest system sizes for different ∆ Jh are shown inFig. 1 (blue).A separate diagnostic for the MBL transition is G =log (cid:16) |(cid:104) n +1 | (cid:98) O | n (cid:105)| ( E n +1 − E n ) (cid:17) , which probes the ratio of matrix ele-ments to energy gaps for a local operator (cid:98) O in nearbyenergy eigenstates [42]. The MB energy spacings are ex-ponentially small in L ; (cid:98) O strongly mixes nearby eigen-states in the thermal phase resulting in (cid:104)G(cid:105) ∼ L , whilethe mixing is exponentially suppressed in the localizedphase resulting in (cid:104)G(cid:105) ∼ − L . Fig. 2b plots (cid:104)G(cid:105) aver-aged over the middle half of the spectrum for various λ and ∆ Jh = 0 with (cid:98) O = σ zL/ . λ c ( L ) is diagnosed bythe change in the sign of the slope of (cid:104)G(cid:105) ( L ). Repeatingthis for different ∆ Jh gives the finite-size λ c estimatesshown in Fig. 1 (black). Note that this diagnostic is moresensitive to thermalization than (cid:104) r (cid:105) and, among the con-sidered diagnostics, gives the smallest λ c (∆ Jh ) finite-sizeestimates, λ c ( L = 14) (cid:39) .
01 at ∆ Jh = 0. We emphasizethat λ is dimensionless, so it is quite striking that inter-actions only 1 / λ c ( L ) and the strong finite-sizedrifts towards smaller λ c ( L ) as L is increased, combinedwith the typical tendency of thermalization to dominatewith increasing system size [27, 43, 44], strongly suggestthat interactions are a relevant perturbation to the non-interacting Ising CP at finite energy densities.Finally, the von Neumann entanglement entropy (EE)of eigenstates is another widely-used diagnostic to probethe MBL transition [20–22]. The EE displays a “volumelaw” scaling in the thermal phase, an area-law scalingin the MBL phase, and a log scaling in the (putative)critical phase [45, 46]. Deep in the MBL PM and SGphases, since the eigenstates are product and cat statesrespectively, the eigenstates EEs are approximately 0 andlog 2, respectively. The existence of three phases in thissystem is clearly detected by the average half-chain eigen-state EEs plotted in Fig. 3a along a path where ∆ Jh istuned from -4 to 4 for a fixed value of (relatively large)interaction strength λ = 0 .
5. Next, to probe the criti-cal behavior, we plot in in Fig. 3b the average entropydensity,
S/S max , as a function of λ at ∆ Jh = 0, where S max = ( L/
2) log 2 is the maximum EE for a subsystem ofsize L/
2. This quantity decreases to 0 with increasing L in the localized/critical phase, while it saturates towards1 in the thermal phase. The finite size estimates for λ c are again obtained from where these curves cross for thelargest pair of L ’s, and are shown in Fig. 1 (green). Thismeasure gives λ c (12 , ∼ = 0 .
04, and hence detects ther-malization more sensitively than (cid:104) r (cid:105) but less sensitivelythan (cid:104)G(cid:105) . III. SUSCEPTIBILITY OF ENTANGLEMENT
The incipient thermalization at small λ motivates usto introduce a new diagnostic, the susceptibility of en-tanglement χ S ( E, L ), which perturbatively probes theeffect of adding weak interactions on the entanglementof eigenstates at energy density
E/L — specifically, theperturbative relevance of adding interactions may mani-fest as a divergence in χ S ( E, L ).To define χ S ( E, L ), we expand the half-chain EE, S ( n ) ( λ, L ), of the n -th many-body eigenstate | n (cid:105) at en-ergy density E/L near interaction strength λ = 0: S ( n ) ( λ, L ) = S ( n )0 ( L ) + λS ( n )1 ( L ) + λ S ( n )2 ( L ) + · · · ,χ S ( E, L ) ≡ (cid:104) S ( n )2 ( L ) (cid:105) typ( E ) , (2)where (cid:104)·(cid:105) typ( E ) denotes the typical (median) value oversamples and eigenstates within small energy densitywindows of δ = 0 .
05 centered around
E/L . While S ( n ) ( λ, L ) can be chosen to be any measure of entan-glement of the eigenstate | n (cid:105) , in this work we choose S ( n ) to be the second Renyi entropy of | n (cid:105) , i.e. S ( n ) ≡− log(Tr[(Tr B | n (cid:105)(cid:104) n | ) )]), where B is taken to be halfthe system, and we observe the same qualitative behav-ior for the von Neumann entropy. This choice of entropyenables a simple analytical expression of S ( n )2 in termsof SP eigenstates (App. B), allowing us to better eluci-date the different contributions to χ S ( E ). It is known Note that S ( n ) j in the expansion of Eq. (2) should not be confusedwith the j th or n th Renyi entropy. − . − . − . . . . . E/L χ S ( E , L ) L = 12L = 14 − − . . . L = 8L = 10
FIG. 4. χ S ( E, L ), obtained using Eqs. (2), (3) as a function ofenergy density at ∆ Jh = 0 (CP, main plot), ∆ Jh = − Jh = +1 (SG, right inset) with (cid:15) = 10 − , ob-tained by taking the median over disorder samples and smallenergy density windows of δ = 0 .
05 across the spectrum. No-tice the different scale of the axes of the insets. The strongincreasing trend with L in the middle of the spectrum at theCP points towards a perturbative instability. Errorbars de-note a 2-percentile interval around the median. that S ( n )0 ( L ) ∼ log( L ) for the marginally localized non-interacting CP at ∆ Jh = 0 [45, 46].A key aim of using this susceptibility as a diagnostic isto “subtract out” the non-interacting contribution to theEE, which overwhelms the value of the EE at the CP atsmall λ at small sizes, for example in Fig. 3b. If interac-tions are perturbatively irrelevant, we expect the higherorder corrections to S ( n )0 ( L ) to be small, and they shouldnot grow with L faster than S ( n )0 ( L ). On the other hand,if interactions are relevant, this may be manifested in astrong growth of S ( n )2 ( L ) with increasing L . Importantly,the behavior of S ( n )2 ( L ) detects a susceptibility towardsthermalization even when the total second order correc-tion — obtained by multiplying S ( n )2 ( L ) by λ / λ and L . Hencethis measure serves as a much more sensitive test of anincipient instability.To start, we can estimate S ( n )1 ( L ) and S ( n )2 ( L ) numer-ically using many-body ED and finite-difference deriva-tives: S ( n )1 ( L ; (cid:15) ) = S ( n ) ( (cid:15),L ) − S ( n ) ( − (cid:15),L )2 (cid:15) S ( n )2 ( L ; (cid:15) ) = S ( n ) ( (cid:15),L )+ S ( n ) ( − (cid:15),L ) − S (0 ,L ) (cid:15) , (3)where S ( n ) ( (cid:15), L ) is the second Renyi entropy of the n thMB eigenstate of the Hamiltonian in Eq. (1) with in-teraction strength λ = (cid:15) . In doing the finite differencecalculations for a given disorder realization, all valuesof the local fields and Ising couplings are kept the samewith only the interaction strength changing. The expres-sions above will agree with the “true” perturbative cor-rections in the RHS of Eq. (2) in the limit (cid:15) →
0. Instead,when the interaction strength (cid:15) is larger than the typicalmany-body level spacing, we expect (avoided) level cross-ings in the many-body spectrum of the finite system. Inthis regime, the interacting MB eigenstate | n ( λ ) (cid:105) can nolonger be interpreted as being perturbatively related tothe non-interacting MB eigenstate | n ( λ = 0) (cid:105) . Hence, toprobe the regime where non-degenerate perturbation the-ory is valid, we use (cid:15) ≤ − , which we numerically findis smaller than the typical many-body level spacing forall the system sizes L ≤ (cid:15) = 10 − which isthree orders of magnitude smaller than the smallest es-timate of λ c ≈ .
01 at these sizes, we find that χ S ( E, L )shows a strong increasing trend with L for finite energydensities at the CP (∆ Jh = 0), strongly indicative of aperturbative instability of localization (Fig. 4). In con-trast, far from the CP in the SG/PM phases, the typi-cal value of χ S ( E, L ) is two orders of magnitude smallerand there is no strong trend with L (shown in the in-sets of Fig. 4), consistent with the stability of localiza-tion in these regimes. Interestingly, we continue to seethis strong growth of χ S ( E, L ) at the CP even for muchlarger values of (cid:15) that are in the non-perturbative regime(but still well below the finite-size λ c estimates).Finally, we note that we have defined χ S ( E, L ) in termsof the second order correction to the entropy rather thanthe first. This is because we find that S ( n )1 ( L ) appearswith a random fluctuating sign, and its typical value isthus small for large system sizes, while S ( n )2 ( L ) typicallygives a strong positive correction. Even upon taking anabsolute value, the typical value of | S ( n )1 ( L ) | is muchsmaller than S ( n )2 ( L ) with a much weaker growth withsystem size. This is also because, S ( n )1 ( L ) admits an ex-pansion as a sum of contributions with fluctuating signsthat give a suppressed correction, as we show in App. B. IV. PERTURBATIVE EXPRESSION OF χ S ( E, L ) Given that a strong growth of χ S ( E, L ) with systemsize is already visible in the regime of non-degenerate per-turbation theory, we now try to understand its behaviorusing a perturbative expansion of the second R´enyi en-tropy S ( n ) ( λ, L ) derived using second-order perturbationtheory on the non-interacting model H (App. B). Themany-body eigenstates and eigenvalues of H (obtainedby ‘filling’ SP orbitals in the fermionic language) are de-noted | ψ n (cid:105) and E n respectively. The expression for S ( n )2 obtained by perturbing | ψ n (cid:105) to second order is of theform (see Eq. (B12)) S ( n )2 = (cid:88) k (cid:54) = n c kn g kn + (cid:88) k (cid:54) = l (cid:54) = n d kln g kn g ln + (cid:88) k (cid:54) = n e kn α kn , (4)where the functional dependence of the quantities on L is implicit. The sums in Eq. (4) run over all many-body eigenstates n and k of H , c kn , d kln , and e kn are O (1)numbers that are related to the properties of the many-body wavefunctions (see Eqs. (B14)-(B16) for their def-initions), and g kn and α kn involve ratios of matrix ele-ments of the interaction and energy denominators thatcan lead to large values of S ( n )2 (see Eq. (B4) for theirdefinitions). Since the full expressions for S ( n )1 / are quitecomplicated, we are not aware of an efficient algorithmto compute all terms in the perturbative correction and exact numerics are still limited to small system sizes ac-cessible to many-body ED. Hence, in the rest of this work,we use numerical observations to perform various approx-imations that allow us to speculate on the behavior of χ S ( E, L ) at larger system sizes than those accessible toexact numerics.First, note that the g kn ’s and α kn ’s in Eq. (4) typicallyappear with fluctuating signs. If one neglects correla-tions between different terms, we expect the second andthird sums to not result in large contributions to S ( n )2 .Indeed, we numerically observe that the typical value of S ( n )2 receives its dominant contribution from the first sum(involving g kn ) when n is a highly excited many-bodyeigenstate, thereby allowing us to approximate χ S ( E, L ) = (cid:104) S ( n )2 ( L ) (cid:105) typ(E) ∼ (cid:104) (cid:88) k (cid:54) = n c kn g kn (cid:105) typ(E) . (5)This approximation is not valid when restricted to theground state, for which we find that the different termsin the expansion conspire and cancel to give a small S ( n )2 without a strong system size dependence, consistent witharguments for the perturbative stability of the groundstate at and away from the CP [35].Next, note that g kn is defined as (see Eq. (B4)) g kn ≡ V kn E kn ≡ (cid:10) ψ k (cid:12)(cid:12) V (cid:12)(cid:12) ψ n (cid:11) E n − E k , (6)where V kn is the matrix element of V between the MBeigenstates of H , | ψ n/k (cid:105) , with MB energies E n/k . Analo-gous to the diagnostic (cid:104)G(cid:105) studied earlier, g kn probes theability of V to generate resonances between the eigen-states of H , and can systematically grow or shrink with L . When written in terms of fermions, V is quarticand fermion parity preserving (see Eqs. (A7) and (A13)).Hence V kn vanishes unless (cid:12)(cid:12) ψ k (cid:11) and (cid:12)(cid:12) ψ n (cid:11) differ in the oc-cupation of two or four single-particle orbitals, and we re-fer to these as “two-orbital” and “four-orbital” processesrespectively. This allows for an efficient polynomial in L computation of all g kn for a given n (App. C).Finally, we observe that c kn is an O (1) number (de-fined in Eq. (B14)) that is bounded as − ≤ c kn ≤ c kn . How-ever, note that: (i) c kn is an O (1) number, so it cannotbe the source of any divergence in χ S ( E, L ) with systemsize and (ii) we have found that the primary role of c kn isto enforce the condition that the entanglement only re-ceives a substantial contribution from processes that tog-gle the occupation of SP orbitals with weight on oppositeor both sides of the entanglement cut. The latter state-ment has been numerically verified, and we observe that c kn is strongly bimodal, peaked at 0 and 4 (Fig. 5a), andthese peaks respectively correspond to whether or not theorbitals involved straddle the entanglement cut. This isconsistent with the fact that resonances between orbitalslocalized on the same side of the entanglement cut can-not contribute to any extra entanglement. As a result,despite the fact that c kn can be negative, the sum in thefirst term of Eq. (4) is typically strongly positive. In thenext section, we examine the behavior of g kn to identifythe processes that could lead to a growth of χ S ( E, L ). V. RESONANT PROCESSES IN THENON-INTERACTING LIMIT
We begin by noting that the distributions of g kn forboth two-orbital and four-orbital processes are broad ona log-scale, so that the behavior of the typical value ofthe sum in Eq. (5) is captured by the behavior of thetypical value of the maximum c kn g kn : χ S ( E, L ) ∼ (cid:104) max k (cid:54) = n (cid:0) c kn g kn (cid:1) (cid:105) typ( E ) . (7)for highly excited n . We have numerically verified thatthis approximation captures the dominant contributionto χ S ( E, L ) and tracks its dependence on L for the sys-tem sizes accessible to ED for which we can explicitlycompute c kn . Using this approximation for χ S ( E, L ), wenow elucidate the non-interacting processes that couldlead to a perturbative instability of the CP for excitedstates at larger sizes than those accessible to many-bodyED. Eq. (7) implies that, for a highly-excited eigenstate | n (cid:105) , the susceptibility of entanglement is dominated bythe most resonant two/four orbital process relativeto | ψ n (cid:105) , with orbitals straddling the entanglement cut.Perturbative stability is determined by the scaling of χ S ( E, L ) with system size L , so we will be concernedwith the scaling of the typical most resonant processesthat straddle the cut. We find that these diverge with L for highly excited states at the CP and, for the sizeswe have been able to study, we find that the four-orbitalprocesses are dominant over the two-orbital ones and cap-ture almost the entire growth of χ S ( E, L ). Note that thissensitivity to the typical most resonant processes is qual-itatively distinct from measures such as (cid:104)G(cid:105) , which areinstead analogous to the typical g nk which is exponen-tially decaying with L at and away from the CP.We now identify possible sources of the instability ofthe non-interacting Ising CP at ∆ Jh = 0. As reviewedin App. D, the SP eigenstates at the CP have a di-verging localization length as the SP energy E → stretched exponentially local-ized: | ψ α ( x ) | ∼ e − √ | ( x − x α ) | /ξ , where ψ α ( x ) denotes the α -th SP wavefunction. Likewise, the density of statesdiverges near
E →
0, so that the typical energy spacingbetween these stretched exponentially localized orbitals also scales as a stretched exponential: δ E ∼ e −√ sL . Bycontrast, away from E = 0 at the CP (or away from theCP, ∆ Jh (cid:54) = 0, at any E ), the SP orbitals are exponen-tially localized, with the energy spacings only fall off aslaws in L . The presence of extended orbitals near zeroenergy can mediate two and four orbital resonances, aswe explain next.We consider first the two-orbital processes which al-low for a more transparent explanation, even though thedominant contribution comes from four-orbital terms.The relevant two-orbital processes are those in whichtwo stretched exponentially localized states, p and q ,with E p,q (cid:39) pq = |E p ± E q | ∼ e − √ R pq / ∆ ,have their occupations toggled between | ψ n (cid:105) and | ψ k (cid:105) .Since V is a sum of local operators, the matrix element V kn ∼ e − √ R pq / Ξ is itself stretched exponentially decay-ing. The localization centers between these states aretypically separated by R pq ∼ L . Thus, g kn = V kn / ∆ pq can be divergent depending on the relative sizes of ∆and Ξ. As we discuss in App. D, the distributions of thestretched-exponential forms of the end-end correlationsand energy gaps for the ground state at the CP were de-rived in Ref. [47], and the scales of those were shown tobe comparable. We similarly expect a comparable scalingfor Ξ and ∆ for the matrix elements of local operatorsand gaps involving extended SP orbitals near E = 0, sothat there is a finite probability for a resonance with adivergent g kn . Fig 5b samples the energy distributionsof the SP orbitals involved in strongly resonant processes( g kn > E = 0. By contrast, in typ-ical processes involving exponentially localized orbitals,the matrix element V kn typically decays exponentiallywith R pq , while the energy difference only decays as apower law so that there is typically no resonance.Note that we are able to go to much larger sizes of L ∼
60 because of the polynomial in L scaling of thetime for computing g kn . However, the random sam-pling of resonant states is not the same as sampling (cid:104) max k (cid:54) = n g kn (cid:105) typ( E ) , which is more computationally in-tensive. However, the random sampling still qualitativelyillustrates the types of processes that can lead to reso-nances at the CP. Note also that the two body processesconsidered here are first-order in the interaction strength,and they probe a qualitatively distinct process from thesecond-order processes analysed in Ref. [36], which alsostudied the stability of marginal Anderson localized sys-tems. Indeed, the Ising CP is stable to the processesconsidered in Ref. [36] (App. E). More precisely, the orbitals at SP energy E are stretched ex-ponentially localized on a length scale ζ ( E ), with a crossover toexponential tails on longer length scales, and ζ ( E ) → L as E → − − − − c kn − − − P ( c k n ) (a) L = 8L = 10L = 12L = 14 . . . . . . . . E α P ( E α ) (b) α = pα = q . . . . . . . . E α P ( E α ) (c) L = 20L = 40L = 60 L = 20L = 40L = 60 α = pα = qα = rα = s FIG. 5. (a) Bimodal distribution of c kn for randomly sampled infinite temperature states (b-c) Distribution of SP energies oforbitals involved in strongly resonant ( g kn > E p < E q ) (c) four orbitals processes ( E p < E q < E r < E s ). Data is obtained by sampling 100 random processes for each of 30-60 many-body eigenstates for 80-250 disorderrealizations for various system sizes, and choosing the strongly resonant processes ( g kn > Next, we consider four-orbital processes involving or-bitals p, q, r, s with energies E p,q,r,s . The energy denom-inator thus reads ∆ pqrs = |E p ± E q ± E r ± E s | , wherethe signs depend on the occupations of the four orbitalsin (cid:12)(cid:12) ψ n (cid:11) , as we explain in App. C. As shown in Fig 5c,typical resonant four-orbital processes at these sizes in-volve 1-2 extended orbitals near SP energy E = 0, and2-3 exponentially localized orbitals with energies awayfrom zero. The wavefunctions of the localized orbitalsoverlap with the stretched exponential wavefunctions inorder to obtain a substantial matrix element mediated bythe delocalized orbital(s). As the system size is increasedand the SP spectrum “fills in”, there are more delocal-ized orbitals near E ≈
0, and the resonant processes mayinvolve more extended orbitals. At the sizes amenableto our analysis, we find that four-orbital processes aredominant over two-orbital ones at the CP. Away fromthe CP, all SP orbitals are exponentially localized andfour-orbital resonances are highly suppresed.As a cautionary remark, we note that while we havefocused on identifying processes at the CP that can giverise to resonant g kn ’s, this alone is not enough to arguefor a diverging χ S ( E, L ), and the c kn ’s in Eq. (7) do playa significant role. It is important that the most reso-nant processes at the CP typically involve one or moreextended orbitals with weight on both sides of the en-tanglement cut, and hence also give a positive c kn andcontribute to S ( n )2 via Eq. (7). In contrast, we find thatthe typical most resonant processes deep in the PM/SGphases, obtained via flipping a few rare resonant SP or-bitals, have a vanishingly small c kn ≈ S ( n )2 . If, instead, we consider theRHS of Eq. (7) which includes c kn then, in the PM/SGphases, the resonances that contribute typically involvelocalized orbitals within an O (1) distance of the entangle-ment cut and there is no strong system size dependence, consistent with the perturbative stability of the PM/SGphases to interactions. Likewise, even though the pres-ence of extended SP orbitals at the CP may, in principle,also destabilize the ground state of the CP, we do notsee any signs of this in practice because (i) the energydenominators in g kn when n is the ground state alwaysinvolve the sum of SP eigenenergies (as opposed to sumsand differences in excited states) and are less likely to beresonant and (ii) when n is the ground state, the approx-imation in Eq. (5) is not a good one and one needs toconsider the full expression Eq. (4) on account of variouscorrelations and cancellations between these terms.In summary, it required a conspiracy of many factorsto conclude that the typical most resonant two/four bodyprocess relative to | ψ n (cid:105) — which involves changing theoccupation of one or more extended SP orbitals at the CP— is able to dominantly capture the behavior of S ( n )2 forexcited states at the CP. We also note that while Eq. (5)only considers a single MB state k relative to n , we arenot suggesting that a single pair of resonant MB states issufficient to thermalize the system at large sizes. Rather,the resonance between the many-body eigenstates ( n, k )and the corresponding strong growth of S ( n )2 ( L ) with L is only capturing the perturbative beginnings of an incip-ient instability towards thermalization. We remind thereader that, at the sizes amenable to ED, the product λ S ( n )2 ( L ) / S ( n )0 ( L )is still very small in absolute terms. However, the strongtrend of growth of S ( n )2 ( L ) with L strongly suggests thatthe nascent signature in χ S ( E, L ) may lead to a ther-malizing cascade across the entire MB energy spectrumas higher order processes and larger sizes are considered.
VI. CONCLUDING REMARKS
We have presented a study of the stability of themarginally Anderson localized spin-glass to paramag-net critical point in a disordered transverse-field Isingmodel. Within many-body ED, we obtain a finite-size estimate for the critical (dimensionless) interactionstrength, λ c ( L ) ∼ λ c ( L )even for modest system sizes L ≤
14, coupled with a drifttowards smaller λ c with increasing L , point to a pertur-bative instability of the CP to the addition of interactions i.e. λ c = 0 in the asymptotic infinite size limit.We introduced a new measure, the “susceptibility ofentanglement” χ S ( E, L ), which perturbatively probesthe effect of adding interactions on the entanglement ofnon-interacting eigenstates. This serves as a much moresensitive probe of an incipient instability, and shows astrong divergence with L at the CP even when estimatedto leading order in an arbitrarily weak interaction. Us-ing a perturbative expansion for χ S ( E, L ), we relatedthe susceptibility to the ratio of matrix elements andenergy differences in the non-interacting problem, andidentified that resonances mediated by extended single-particle states at the CP are the leading order processescontributing to the growth of χ S ( E, L ) with system size L . At these sizes and interaction strengths, the absolute(in magnitude) correction to the non-interacting entan-glement is still very small; but the strong divergence of χ S ( L ) with L points to the perturbative beginnings of anincipient instability that could lead to full thermalizationacross the MB energy spectrum as higher-order processesand larger sizes are considered.An important point is that a divergence in χ S ( E, L )may be caused even if the addition of interactions leadsto a discontinuous change in the critical eigenstate en-tanglement entropy, for example from S ∼ c log L in thenon-interacting limit to S ∼ c log L in the interactingproblem. While such a change does not correspond tothermalization, it still points to a relevance of interac-tions at the non-interacting infinite-randomness criticalpoint. This (weaker) effect seems unlikely for highly-excited eigenstates at the CP, which already show signa-tures of thermalization for small λ in finite size ED stud-ies. However, it would be interesting to examine whetherinteractions might be relevant in this weaker sense for theground state phase transition. Indeed, revisiting variousstrong randomness RG treatments of the Ising transi-tion [12, 37] - which explicitly ignore the possibility thatthe interactions are relevant - with these considerationsin mind is an important direction for future work.While this paper has focused on a perturbative insta-bility of the CP to the addition of interactions, we notethat the system may also be subject to non-perturbativeinstabilities on account of the diverging SP localizationlength near E = 0 [6]. These effects, if present, would besubdominant to the perturbative processes, and wouldonly be visible at asymptotically larger sizes than thoseconsidered here. However, these may play a dominant role once the system is perturbed slightly away fromcriticality, in which case the non-interacting localizationlength ξ ( E ) remains finite but may get very large as E →
0. Due to the energy dependence of the localizationlength, this case is not directly covered by the argumentsin Ref. [6] which assume a uniform ξ across the SP spec-trum and predict an instability once ξ exceeds an O (1)threshold. Instead, this requires a more nuanced analy-sis along the lines of the recent study in Ref. [48], whichallowed for a distribution in ξ ( E ).We reiterate that an instability of the marginally An-derson localized CP corresponds to (asymptotic) ther-malization in the presence of interactions, independentof the details or strength of the disorder configurations.However, stronger disorder, for instance generated bypower law distributions, is much harder to analyze nu-merically because of strong finite-size effects. The power-law disorder generates tails to very small couplings, whicheffectively decouples the system into smaller pieces andcan make the system look more localized than it is. Wehave repeated our analysis for power-law disorder andbroadly found qualitatively similar behavior. However,even for (not too strongly disordered) power-law distri-butions with exponent 2, we find that the many-bodyenergy spacings become small enough that a perturba-tive estimate of S ( L ; (cid:15) ) for small (cid:15) runs into machineprecision issues. Hence prior numerical studies of MBL-to-MBL phase transitions that use power-law distribu-tions with large exponents [38, 39] should be interpretedwith caution.Finally, while we have focused on the PM to SGphase transition in a disordered TFIM, our considera-tions are expected to apply more broadly to various pu-tative MBL-to-MBL phase transitions with a marginallyAnderson localized non-interacting counterpart. Someof these transitions, for example those between the “ π spin-glass” (or discrete time crystal) phase and variousparamagnetic phases in a driven Ising model [13] are inthe infinite randomness universality class as the disor-dered TFIM [49] and will be subject to a similar insta-bility once interactions are added. Separately, Ref. [39]considered a “particle-hole” symmetric disordered XXZchain; they found that the addition of interactions tothe non-interacting (critical) disordered XX chain couldproduce a localized spin glass phase with spontaneouslybroken particle-hole symmetry. In the fermion language,the non-interacting limit of their model corresponds totwo decoupled and critical Majorana chains with an addi-tional reflection symmetry, which enabled the possibilityof localization via symmetry breaking on adding interac-tions — a possibility that is absent in the critical randomTFIM which corresponds to a single critical Majoranachain. The disordered XXZ model represents a particle-hole symmetric slice through a broader class of disorderedXYZ spin chains with localized spin-glass phases pointingalong the x/y/z directions [38]. If we consider the phasediagram in this broader parameter space, our considera-tions are expected to apply to transitions between glassyphases ordered in different directions. This model wasstudied via “spectrum-bifurcation” RG in Ref. [38], andthis RG scheme again did not include the possibility ofinstabilities at the critical lines. Ref. [38] also presentedan ED analysis, and found that they needed power-lawdisorder with a large exponent (greater than 4) to pre-vent thermalization at the critical lines — but this regimeis not trustworthy for the small sizes amenable to theirED analysis, as discussed earlier.In all, our work adds to the growing body of work onnon-equilibrium quantum criticality, and addresses long-standing open questions about the nature and stabilityof MBL-to-MBL phase transitions. We expect χ S ( E, L )may be of independent interest for studying the effect ofinteractions in a variety of settings, for example to probethe existence of many-body mobility edges, another ma-jor outstanding question in the literature. While we hadto rely on various approximations to study χ S ( E, L ), itwould also be interesting to see if more exact methodscould be developed to study this quantity at larger sizes,given its perturbative nature. Separately, different tech-niques such as matrix product state based methods [50–53], or numerical linked cluster expansions [43], or ma-chine learning techniques [44, 54–56] may prove useful formore large scale analyses.
Note Added—
During the completion of this work, webecame aware of complementary work on the presenceof intervening thermal phases between MBL transitionswhich will appear in the same arXiv posting [57].
Acknowledgements—
We thank Subroto Mukerjee for aninitial collaboration, and Trithep Devakul, Chris Lau-mann, Siddharth Parameswaran, Shivaji Sondhi and Ro-main Vasseur for helpful discussions. This work was sup-ported in part by the DARPA DRINQS program.
Appendix A: Jordan-Wigner Transformations ofEq. (1)
In this appendix, we perform a Jordan-Wigner trans-formation on the Hamiltonian Eq. (1) with open bound-ary conditions. We split the Hamiltonian of Eq. (1) as H = H + λV (A1), where H is the disordered transverse field Ising model H = L (cid:88) j =1 h j σ zj + L − (cid:88) j =1 J j σ xj σ xj +1 , (A2)and V is the Kramers-Wannier self-dual interaction thatreads V = L − (cid:88) j =1 ( h ∗ σ zj σ zj +1 + J ∗ σ xj σ xj +2 ) . (A3) H has a Z parity symmetry, defined by the operator P = L (cid:89) j =1 σ zj , [ P, H ] = 0 , [ P, V ] = 0 . (A4)To perform the Jordan-Wigner transformation tofermions with creation and annhilation operators c j and c † j , we apply the transformations: σ + j = ( − (cid:80) k 1) + J j ( c † j c j +1 + c † j c † j +1 + h.c.) (cid:17) , (A6)and the interaction V maps onto V = L − (cid:80) j =1 ( h ∗ (2 n j − 1) (2 n j +1 − J ∗ ( − n j +1 ( c † j c j +2 + c † j c † j +2 + h.c.)) . (A7)The expression for V can be simplified using the prop-erty ( − n j = 1 − n j . Furthermore, under the Jordan-Wigner transformation of Eq. (A5), the parity operator P of Eq. (A4) maps onto P = ( − L × ( − L (cid:80) j =1 n j (A8)Since a superconducting Hamiltonian is better ex-pressed in terms of Majorana fermions, we apply a secondset of substitutions c † j = 12 ( χ j − + iχ j ) c j = 12 ( χ j − − iχ j ) (A9)where the χ j − ’s and χ j ’s are Majorana fermions thatobey the commutation relations { χ a , χ b } = 2 δ ab . (A10)In terms of the Majorana fermions, the fermion numberand parity operators read n j = 12 (1 − iχ j − χ j ) , P = ( − i ) L L (cid:89) l =1 χ l , (A11)and H and V simplify to H = − i L − (cid:80) j =1 ( h j χ j − χ j + J j χ j +1 χ j ) (A12) V = − L − (cid:80) j =1 ( h ∗ χ j − χ j χ j +1 χ j +2 + J ∗ χ j χ j +1 χ j +2 χ j +3 ) (A13)0At the critical point (when h ∗ = J ∗ ≡ C ), V can bewritten more elegantly as V = C L − (cid:88) j =1 χ j χ j +1 χ j +2 χ j +3 , (A14) Appendix B: Perturbative expansion of the secondR´enyi entropy In this appendix, we compute the expression for thesusceptibility of entanglement entropy in terms of ma-trix elements of the non-interacting problem. For sim-plicity, we use choose the second R´enyi entropy to be theentanglement entropy. That is, S ( λ, L ) ≡ − log (cid:0) Tr ρ (cid:1) (B1)is the entanglement entropy of a particular eigenstate atan interaction strength λ and system size L . We expand S ( λ, L ) in powers of λ as shown in Eq. (2). S ( λ, L ) = S ( L ) + λS ( L ) + λ S ( L ) + . . . (B2)where S ( L ) is the entanglement entropy in the non-interacting limit and S ( L ) and S ( L ) are the first and second derivatives of S ( λ, L ). To compute S ( L ) and S ( L ) perturbatively, we start with a non-interactingmany-body wavefunction (cid:12)(cid:12) ψ n (cid:11) , and write the perturbedwavefunction | ψ n (cid:105) up to O (cid:0) λ (cid:1) as | ψ n (cid:105) = (cid:32) − λ (cid:80) k (cid:54) = n | g kn | (cid:33) (cid:12)(cid:12) ψ n (cid:11) + λ (cid:80) k (cid:54) = n g kn (cid:12)(cid:12) ψ k (cid:11) + λ (cid:80) k (cid:54) = n α kn (cid:12)(cid:12) ψ k (cid:11) + O (cid:0) λ (cid:1) , (B3)where { (cid:12)(cid:12) ψ k (cid:11) } is the set of unperturbed non-interactingmany-body wavefunctions and g kn ≡ (cid:104) ψ k | V | ψ n (cid:105) E n − E k ,α kn ≡ (cid:80) l (cid:54) = n (cid:18) (cid:104) ψ k | V | ψ l (cid:105)(cid:104) ψ l | V | ψ n (cid:105) ( E n − E l )( E n − E k ) (cid:19) − (cid:104) ψ n | V | ψ n (cid:105)(cid:104) ψ k | V | ψ n (cid:105) ( E n − E k ) , (B4)where V is the interaction, E n and E k are the many-body energies of (cid:12)(cid:12) ψ n (cid:11) and (cid:12)(cid:12) ψ k (cid:11) respectively. Note thatin the cases we are working with, the Hamiltonian H andinteraction V are time-reversal symmetric, and thus g kn and α kn are real numbers, which simplifies the followinganalysis. Using Eq. (B3) we expand the density matrixof | ψ n (cid:105) as | ψ n (cid:105) (cid:104) ψ n | = (cid:32) − λ (cid:80) k (cid:54) = n g kn (cid:33) (cid:12)(cid:12) ψ n (cid:11) (cid:10) ψ n (cid:12)(cid:12) + λ (cid:80) k (cid:54) = n (cid:2) g kn (cid:12)(cid:12) ψ k (cid:11) (cid:10) ψ k (cid:12)(cid:12)(cid:3) + λ (cid:80) k (cid:54) = l (cid:54) = n (cid:2) g kn g ln (cid:0)(cid:12)(cid:12) ψ k (cid:11) (cid:10) ψ l (cid:12)(cid:12) + (cid:12)(cid:12) ψ l (cid:11) (cid:104) ψ k | (cid:1)(cid:3) + (cid:32) − λ (cid:80) k (cid:54) = n g kn (cid:33) λ (cid:80) k (cid:54) = n (cid:2) g kn (cid:0)(cid:12)(cid:12) ψ k (cid:11) (cid:10) ψ n (cid:12)(cid:12) + (cid:12)(cid:12) ψ n (cid:11) (cid:10) ψ k (cid:12)(cid:12)(cid:1)(cid:3) + λ (cid:80) k (cid:54) = n α kn (cid:0)(cid:12)(cid:12) ψ k (cid:11) (cid:10) ψ n (cid:12)(cid:12) + (cid:12)(cid:12) ψ n (cid:11) (cid:10) ψ k (cid:12)(cid:12)(cid:1) (B5)Defining the matrices ρ red ≡ Tr B ( | ψ n (cid:105) (cid:104) ψ n | ) , ρ ll ≡ Tr B (cid:0)(cid:12)(cid:12) ψ l (cid:11) (cid:10) ψ l (cid:12)(cid:12)(cid:1) ρ lm ≡ Tr B (cid:0)(cid:12)(cid:12) ψ l (cid:11) (cid:10) ψ m (cid:12)(cid:12) + (cid:12)(cid:12) ψ m (cid:11) (cid:10) ψ l (cid:12)(cid:12)(cid:1) , (B6)where Tr B represents the trace over a subsystem, using Eq. (B5) we obtain ρ red = (1 − λ (cid:88) k (cid:54) = n g kn ) ρ nn + λ (cid:88) k (cid:54) = n (cid:2) g kn ρ kk + 2 α kn ρ kn (cid:3) + λ (cid:88) k (cid:54) = k (cid:48) (cid:54) = n [ g kn g k (cid:48) n ρ kk (cid:48) ] + 2 λ (cid:88) k (cid:54) = n [ g kn ρ kn ] + O ( λ ) , (B7)and consequently ρ = (cid:34)(cid:32) − λ (cid:80) k (cid:54) = n g kn (cid:33) ρ nn (cid:35) + 4 λ (cid:80) k (cid:54) = n g kn ρ kn ρ nn + 2 λ (cid:32) (cid:80) k (cid:54) = n (cid:2) g kn ρ kk ρ nn + 2 α kn ρ kn ρ nn (cid:3) + (cid:80) k (cid:54) = k (cid:48) (cid:54) = n g kn g k (cid:48) n ρ kk (cid:48) ρ nn (cid:33) + 4 λ (cid:34) (cid:80) k (cid:54) = n g kn ρ kn (cid:35) + O (cid:0) λ (cid:1) . (B8)1Restricting to O (cid:0) λ (cid:1) , taking a trace, and rearranging terms, we obtainTr (cid:0) ρ (cid:1) = Tr( ρ nn ) + 4 λ (cid:80) k (cid:54) = n [ g kn Tr ( ρ kn ρ nn )] + 2 λ (cid:32) (cid:80) k (cid:54) = k (cid:48) (cid:54) = n [ g kn g k (cid:48) n Tr ( ρ kk (cid:48) ρ nn + 2 ρ kn ρ k (cid:48) n )]+ (cid:80) k (cid:54) = n (cid:2) g kn (cid:0) Tr (cid:0) ρ kk ρ nn + 2 ρ kn − ρ nn (cid:1)(cid:1) + 2 α kn ρ kn ρ nn (cid:3)(cid:33) + O (cid:0) λ (cid:1) . (B9)Thus, using Eq. (B1) an expression for S ( λ, L ) that reads S ( λ, L ) = − log Tr (cid:0) ρ nn (cid:1) − λ (cid:80) k (cid:104) g kn Tr( ρ kn ρ nn )Tr( ρ nn ) (cid:105) + 8 λ (cid:32) (cid:80) k (cid:54) = n (cid:104) g kn Tr( ρ kn ρ nn )Tr( ρ nn ) (cid:105)(cid:33) − λ (cid:32) (cid:80) k (cid:54) = n (cid:20) g kn Tr ( ρ kk ρ nn +2 ρ kn − ρ nn ) Tr( ρ nn ) (cid:21) +2 (cid:80) k (cid:54) = n (cid:104) α kn Tr( ρ kn ρ nn )Tr( ρ nn ) (cid:105) + (cid:80) k (cid:54) = l (cid:54) = n (cid:104) g kn g ln Tr ( ρ kl ρ nn +2 ρ kn ρ ln )Tr( ρ nn ) (cid:105)(cid:33) ≡ S ( L ) + λS ( L ) + λ S ( L ) + O (cid:0) λ (cid:1) . (B10)Consequently, S ( L ) = − (cid:80) k (cid:54) = n (cid:104) g kn Tr( ρ kn ρ nn )Tr( ρ nn ) (cid:105) , (B11) S ( L ) = − (cid:80) k (cid:20) g kn Tr ( ρ kk ρ nn +2 ρ kn − ρ nn ) Tr( ρ nn ) (cid:21) + 16 (cid:32) (cid:80) k (cid:54) = n (cid:104) g kn Tr( ρ kn ρ nn )Tr( ρ nn ) (cid:105)(cid:33) − (cid:80) k (cid:54) = l (cid:54) = n (cid:104) g kn g ln Tr( ρ kl ρ nn +2 ρ kn ρ ln )Tr( ρ nn ) (cid:105) + 4 (cid:80) k (cid:54) = n (cid:104) α kn Tr( ρ kn ρ nn )Tr( ρ nn ) (cid:105) (B12) S ( L ) in Eq. (B12) can further be written as S ( L ) = (cid:80) k (cid:54) = n g kn (cid:18) − ρ kk ρ nn )Tr( ρ nn ) − ( ρ kn ) Tr( ρ nn ) + 16 Tr( ρ kn ρ nn ) Tr( ρ nn ) (cid:19) + (cid:80) k (cid:54) = l (cid:54) = n g kn g ln (cid:16) 16 Tr( ρ kn ρ nn )Tr( ρ ln ρ nn )Tr( ρ nn ) − ρ kl ρ nn )Tr( ρ nn ) + ρ kn ρ ln )Tr( ρ nn ) (cid:17) + 4 (cid:80) k (cid:54) = n α kn Tr( ρ kn ρ nn )Tr( ρ nn ) ≡ (cid:80) k (cid:54) = n c kn g kn + (cid:80) k (cid:54) = l (cid:54) = n d kln g kn g ln + (cid:80) k (cid:54) = n e kn α kn , (B13)where c kn , d kln , and e kn , where they are all bounded O (1) quantities defined as c kn ≡ (cid:18) − Tr( ρ kk ρ nn )Tr( ρ nn ) − ( ρ kn ) Tr( ρ nn ) + ρ kn ρ nn ) Tr( ρ nn ) (cid:19) (B14) d kln = 4 (cid:16) ρ kn ρ nn )Tr( ρ ln ρ nn )Tr( ρ nn ) − Tr( ρ kl ρ nn )Tr( ρ nn ) + ρ kn ρ ln )Tr( ρ nn ) (cid:17) (B15) e kn = Tr( ρ kn ρ nn )Tr( ρ nn ) . (B16) Appendix C: Computation of matrix elements andenergy gaps in the non-interacting model In this section we exemplify the computation of ma-trix elements of the interaction and energy differencesbetween the eigenstates of the non-interacting Hamilto-nian of Eq. (A2) (Eq. (A12) in the Majorana fermionlanguage). We now recall the construction of many-bodyeigenstates of the non-interacting Hamiltonian. H inEq. (A12) can be written as a 2 L × L Hermitian matrix in the basis of the Majorana fermions as H = χ † M χ, (C1)where χ is ( χ χ . . . χ L ) T , a vector of Majoranafermions. M is Hermitian as well as anti-symmetric, thusits eigenvalues occur in pairs of real (+ E , −E ), where weassume E ≥ 0. If the eigendecomposition of M reads M = Q Λ Q † , (C2)where Λ is a diagonal, we obtain H = χ † Q Λ Q † χ. (C3)2The creation and annihilation operators of the single-particle eigenstates of H { b n } and { b † n } are then encodedin the vectors b = ( b b · · · b L b † b † · · · b † L ) T and b † ,where b = Q † χ, b † = χ † Q. (C4)Thus the eigenstates corresponding to single-particle en-ergies + E α and −E α are b † α | (cid:105) and b α | (cid:105) , where | (cid:105) is theFock vacuum that satisfies c j | (cid:105) ∀ j . Using Eq. (C4),with an appropriate labelling of the indices of Q , theMajorana fermions can be written in terms of b α ’s and b † α ’s as χ i = L (cid:88) α =1 (cid:0) Q i,α b † α + Q i, − α b α (cid:1) ≡ (cid:88) α (cid:54) =0 Q i,α b tα , (C5)where for convenience we introduced the notation b tα ≡ (cid:26) b †| α | if α > b | α | if α < , (C6)In Eq. (C5), since χ i ’s are real, the Q i,α ’s satisfy Q ∗ i, − α = Q i,α . (C7)The interaction of Eq. (A14) at the critical point canthen be written in terms of the b tα ’s as V = C (cid:88) i,α,β,γ,δ Q i,α Q i +1 ,β Q i +2 ,γ Q i +3 ,δ b tα b tβ b tγ b tδ ≡ (cid:88) α,β,γ,δ V ( α,β,γ,δ ) b tα b tβ b tγ b tδ , (C8)where we have defined V ( α,β,γ,δ ) ≡ C (cid:88) i Q i,α Q i +1 ,β Q i +2 ,γ Q i +3 ,δ . (C9)A similar expression can be obtained for the interactionaway from the critical point. Using Eq. (C8), we want toobtain matrix elements between many-body eigenstates (cid:12)(cid:12) ψ n (cid:11) and (cid:12)(cid:12) ψ k (cid:11) . That is, we want to compute V kn ≡ (cid:10) ψ k (cid:12)(cid:12) V (cid:12)(cid:12) ψ n (cid:11) = (cid:88) α,β,γ,δ V ( α,β,γ,δ ) (cid:10) ψ k (cid:12)(cid:12) b tα b tβ b tγ b tδ (cid:12)(cid:12) ψ n (cid:11) . (C10) In Eq. (C10), it is clear that all the matrix elements (cid:10) ψ k (cid:12)(cid:12) b tα b tβ b tγ b tδ (cid:12)(cid:12) ψ n (cid:11) vanish unless (cid:12)(cid:12) ψ n (cid:11) and (cid:12)(cid:12) ψ k (cid:11) differ inthe occupation of four or two of the single-particle or-bitals, the computation of which we illustrate separately.We first introduce the notations and conventions usedin the following subsections. For any tuple A =( x , x , · · · , x n ), we introduce − A ≡ ( − x n , − x n − , · · · , − x ) . (C11)We further introduce products → (cid:89) x ∈ A f ( x ) = f ( x ) f ( x ) · · · f ( x n ) , ← (cid:89) x ∈ A f ( x ) = f ( x n ) f ( x n − ) · · · f ( x ) . (C12)We also introduce a tuple concatenation operator ◦ that acts on tuples A = ( x , x , · · · , x n A ) and B =( y , y , · · · , y n B ) as C = A ◦ B = ( x , x , · · · , x n A , y , y , · · · , y n B ) . (C13)We will also be using the usual set operation / (differ-ence) for tuples instead. Furthermore, for any n -tuple A , for any permutation σ that belongs to the permuta-tion group S n , we denote the corresponding permutationof A as σ ( A ). Finally, we define tuples { Λ n } with ele-ments ordered in ascending order as the tuple containingthe indices of the single-particle orbitals that are occu-pied in the many-body states { (cid:12)(cid:12) ψ n (cid:11) } . Consequently, theexpression for the many-body state (cid:12)(cid:12) ψ n (cid:11) reads (cid:12)(cid:12) ψ n (cid:11) = → (cid:89) α ∈ Λ n b † α | θ (cid:105) , (C14)where | θ (cid:105) is the Bogoliubov vacuum defined by b α | θ (cid:105) = 0 ∀ α. (C15) 1. Four-orbital processes We first consider the matrix element of the interaction between two many-body eigenstates (cid:12)(cid:12) ψ n (cid:11) and (cid:12)(cid:12) ψ k (cid:11) whichdiffer in the occupation of four of the single-particle orbitals, say the orbitals O = ( p, q, r, s ) where 0 < p < q < r < s .3Here, the matrix element can be non-zero only if (cid:12)(cid:12) ψ n (cid:11) and (cid:12)(cid:12) ψ k (cid:11) are of the forms (cid:12)(cid:12) ψ n (cid:11) = → (cid:89) α ∈ Λ n b † α | θ (cid:105) = ( − (cid:80) α ∈O nα − (cid:80) β =1 n β (cid:124) (cid:123)(cid:122) (cid:125) η n → (cid:89) α ∈O n b † α | ψ (cid:105) , (cid:12)(cid:12) ψ k (cid:11) = → (cid:89) α ∈ Λ k b † α | θ (cid:105) = ( − (cid:80) α ∈O kα − (cid:80) β =1 n β (cid:124) (cid:123)(cid:122) (cid:125) η k → (cid:89) α ∈O k b † α | ψ (cid:105) , (C16)where O n , O k ⊆ O are disjoint tuples such that O n = O / O k , and | ψ (cid:105) is a many-body eigenstate in which the orbitals p , q , r and s are unoccupied, i.e. | ψ (cid:105) = → (cid:89) α ∈ Λ n / O n b † α | θ (cid:105) , (C17)and n β is the occupation number of the single-particle orbital β in | ψ (cid:105) . The matrix element of Eq. (C10) then reads V kn = η k η n (cid:88) σ ∈ S V σ ( O k ◦−O n ) (cid:104) ψ | ← (cid:89) α ∈O k b α → (cid:89) α ∈ σ ( O k ◦−O n ) b tα → (cid:89) α ∈O n b † α | ψ (cid:105) = ( − (cid:80) α ∈O kα − (cid:80) β =1 n β + (cid:80) α ∈O nα − (cid:80) β =1 n β (cid:88) σ ∈ S (cid:0) sgn ( σ ) V σ ( O k ◦−O n ) (cid:1) (cid:104) ψ | ← (cid:89) α ∈O k b α → (cid:89) α ∈O k b † α ← (cid:89) α ∈O n b α → (cid:89) α ∈O n b † α | ψ (cid:105) = ( − (cid:80) α ∈O α − (cid:80) β =1 n β (cid:88) σ ∈ S (cid:0) sgn ( σ ) V σ ( O k ◦−O n ) (cid:1) = ( − q − (cid:80) β = p n β + s − (cid:80) β = r n β (cid:88) σ ∈ S (cid:0) sgn ( σ ) V σ ( O k ◦−O n ) (cid:1) . (C18)The energy differences between E n and E k then read E kn ≡ E n − E k = (cid:88) α ∈O n E α − (cid:88) α ∈O k E α . (C19)An important observation is that the magnitudes of V kn and E kn in Eqs. (C18) and (C19) do not depend on theoccupations of any of the single-particle orbitals apart from the ones involved. 2. Two-orbital processes We now discuss the case where (cid:12)(cid:12) ψ n (cid:11) and (cid:12)(cid:12) ψ k (cid:11) differ in the occupation of two single-particle orbitals, O = ( p, q )where 0 < p < q . Here too, the matrix element of the interaction can be non-zero only if (cid:12)(cid:12) ψ k (cid:11) and (cid:12)(cid:12) ψ n (cid:11) are of theforms shown in Eq. (C16), where O = ( p, q ). The two-orbital matrix element of the interaction operator reads V kn = η k η n (cid:88) γ (cid:88) σ ∈ S V σ ( O k ◦ ( γ, − γ ) ◦−O n ) (cid:104) ψ | ← (cid:89) α ∈O k b α → (cid:89) α ∈ σ ( O k ◦ ( γ, − γ ) ◦−O n ) b tα → (cid:89) α ∈O n b † α | ψ (cid:105) = η k η n (cid:34)(cid:88) γ (cid:88) σ ∈ S (cid:0) sgn ( σ ) V σ ( O k ◦ ( γ, − γ ) ◦−O n ) (cid:1) (cid:104) ψ | ← (cid:89) α ∈O k b α → (cid:89) α ∈O k b † α (cid:0) b † γ b γ (cid:1) ← (cid:89) α ∈O n b α → (cid:89) α ∈O n b † α | ψ (cid:105) Θ ( σ ( O k ◦ ( γ, − γ ) ◦ −O n ) , γ, − γ )+ (cid:88) γ (cid:88) σ ∈ S (cid:0) sgn ( σ ) V σ ( O k ◦ ( − γ,γ ) ◦−O n ) (cid:1) (cid:104) ψ | ← (cid:89) α ∈O k b α → (cid:89) α ∈O k b † α (cid:0) b γ b † γ (cid:1) ← (cid:89) α ∈O n b α → (cid:89) α ∈O n b † α | ψ (cid:105) Θ ( σ ( O k ◦ ( − γ, γ ) ◦ −O n ) , − γ, γ ) (cid:35) , (C20)where we have defined a Θ symbol for a tuple A and elements a , b of the tuple asΘ ( A, a, b ) = (cid:26) a appears to the left of b in A a appears to the right of b in A . (C21)4Thus, the matrix element of Eq. (C20) can be written as V kn = ( − (cid:80) α ∈O kα − (cid:80) β =1 n β + (cid:80) α ∈O nα − (cid:80) β =1 n β (cid:34)(cid:88) γ (cid:88) σ ∈ S (cid:0) sgn ( σ ) V σ ( O k ◦ ( γ, − γ ) ◦−O n ) δ n γ , Θ ( σ ( O k ◦ ( γ, − γ ) ◦ −O n ) , γ, − γ ) (cid:1) + (cid:88) γ (cid:88) σ ∈ S (cid:0) sgn ( σ ) V σ ( O k ◦ ( − γ,γ ) ◦−O n ) δ n γ , Θ ( σ ( O k ◦ ( − γ, γ ) ◦ −O n ) , − γ, γ ) (cid:1)(cid:35) = ( − (cid:80) α ∈O α − (cid:80) β =1 n β (cid:34)(cid:88) γ (cid:88) σ ∈ S (cid:0) sgn ( σ ) V σ ( O k ◦ ( γ, − γ ) ◦−O n ) δ n γ , Θ ( σ ( O k ◦ ( γ, − γ ) ◦ −O n ) , γ, − γ ) (cid:1) − (cid:88) γ (cid:88) σ ∈ S (cid:0) sgn ( σ ) V σ ( O k ◦ ( γ, − γ ) ◦−O n ) δ n γ , Θ ( σ ( O k ◦ ( γ, − γ ) ◦ −O n ) , − γ, γ ) (cid:1)(cid:35) = ( − q − (cid:80) β = p n β (cid:88) γ ( − n γ (cid:88) σ ∈ S (cid:0) sgn ( σ ) V σ ( O k ◦ ( γ, − γ ) ◦−O n ) Θ (cid:0) σ ( O k ◦ ( γ, − γ ) ◦ −O n ) , ( − n γ γ, ( − n γ +1 γ (cid:1)(cid:1) , (C22)where we have used the facts that (cid:104) ψ | ← (cid:89) α ∈O k b α → (cid:89) α ∈O k b † α (cid:0) b † γ b γ (cid:1) ← (cid:89) α ∈O n b α → (cid:89) α ∈O n b † α | ψ (cid:105) = δ n γ , , (cid:104) ψ | ← (cid:89) α ∈O k b α → (cid:89) α ∈O k b † α (cid:0) b γ b † γ (cid:1) ← (cid:89) α ∈O n b α → (cid:89) α ∈O n b † α | ψ (cid:105) = δ n γ , . (C23)Meanwhile, the energy difference between E n and E k reads E kn ≡ E n − E k = (cid:88) α ∈O n E α − (cid:88) α ∈O k E α . (C24)Unlike four-orbital processes, the magnitude of V kn in Eq. (C22) does depend on the occupations of the orbitals otherthan the ones directly involved in the process. Appendix D: Review of the non-interacting model In this appendix, we review the properties of the non-interacting ( λ = 0) limit of the Hamiltonian Eq. (1).In this limit, the Hamiltonian is the well-known dis-ordered transverse field Ising model of Eq. (A2) whichhas been widely studied in literature in various con-texts [35, 47, 58–63]. For concreteness, we assume uni-form distributions in the couplings { J i } ∈ [0 , J ∗ ] andfields { h i } ∈ [0 , h ∗ ] in Eq. (A2) and open boundaryconditions. Since H is statistically Kramers-Wannierself-dual, its eigenstates undergo a phase transition be-tween a “spin-glass” phase to the paramagnetic phase atlog J ∗ = log h ∗ . As shown in Eq. (A12) of App. A, theHamiltonian of Eq. (A2) can be written as a quadraticHamiltonian of 2 L Majorana fermions. The propertiesof H can thus be understood using the single particleeigenstates of Eq. (A12).The model of Eq. (A2) has been solved using theStrong Disorder Renormalization Group (SDRG), a real-space renormalization group [35, 64]. This RG procedureproceeds by diagonalizing the strongest on-site/nearest-neighbor term in the Hamiltonian, projecting onto itslow-energy subspace and decimating the site/bond asso-ciated with that term. Each step of the RG moves lowerin energy and changes the distributions of the couplings { J i } and the fields { h i } [64]. The nature of the RG fixedpoint is different at the critical point (where { h i } and { J i } are chosen from identical distributions) and awayfrom the critical point [35].At the critical point, typical end-end correlation func-tions (cid:104) σ x σ xL (cid:105) GS in the many-body ground state scale stretched exponentially , and so does the typical energygap ∆ E GS above the many-body ground state [35, 47,59, 64]. That is, for a system size of L with open bound-ary conditions, (cid:104)(cid:104) σ x σ xL (cid:105) GS (cid:105) typ ∼ e − √ L Ξ , (cid:104) ∆ E GS (cid:105) typ ∼ e − √ L ∆ , (D1)where the distributions of Ξ and ∆ are derived inRef. [47]. In terms of the single-particle eigenstates ofEq. (A12), we find that a few of the single-particle eigen-states close to E = 0 are “extended”, or more preciselystretched exponentially localized at the critical point,consistent with the scaling of the typical correlations inEq. (D1). That is, we find that the wavefunction of the α -th single-particle eigenstate (one with energy E α ) hasthe form | ψ α ( x ) | ∼ exp (cid:18) − (cid:113) | x − R α | ξ ext (cid:19) if E α ≈ (cid:16) − | x − R α | ξ loc (cid:17) otherwise , (D2)5 i − − − | Q i | (a) i − − − | Q i | (b) − − − − E . . . ρ ( E ) (c) − − − E . . . h I P R i (d) − − − E h ξ i (e) L h ξ e x t i (f) FIG. 6. Various properties of the non-interacting model at the critical point ∆ Jh = 0 or log J ∗ = log h ∗ (a-b) Spatial profileof a typical single-particle eigenstate with energy (a) E ≈ E ∼ O (1) for L = 5000. Figures show the weights | Q i | of theeigenstates on the i ’th Majorana fermion χ i . (c) Single-particle density of states (d) Typical Inverse Participation Ratio (IPR)of the single-particle eigenstates. A lower IPR indicates lesser localization. (e) Typical second moment of the single-particleeigenstates. (f) Growth of the typical second moment of the single-particle eigenstates closest to E = 0 with system size L . where R α is the localization center, ξ ext and ξ loc arethe second moments of the wavefunctions. Examples ofsingle-particle wavefunctions at the critical point close toand away from E = 0 are shown in Figs. 6a-b. Indeed,the typical IPR’s (resp. second moments) of the single-particle orbitals appear to decrease (resp. increase), asshown in Fig. 6d (resp. Fig. 6e). However, the secondmoments of the stretched exponentially localized orbitalssaturate to a constant for large L . The growth of thesecond moment of the single-particle orbital closest to E = 0 with system size is shown in Fig. 6f. Furthermore,as shown in Fig. 6c the single-particle density of statesat the critical point diverges as E → 0, and scales as ρ ( E ) ∼ dkd E ∼ L dLd E ∼ − E (log E ) , (D3)where k denotes momentum. These properties can alsobe derived using the fact that the (positive part of the)single-particle spectrum of Eq. (A12) at the critical pointis identical to that of the well-studied one-dimensionalfermion random hopping model [34, 65–67].Meanwhile, away from the critical point in the Hamil-tonian Eq. (A2), the typical correlation functions andenergy gaps in the ground state [35, 59] (cid:104)(cid:104) σ x σ xL (cid:105) GS (cid:105) typ ∼ exp (cid:18) − LM (cid:19) , (cid:104) ∆ E GS (cid:105) typ ∼ L δ . (D4)Furthermore, the single-particle spectrum away from thecritical point shows a uniform density of states, and thesingle-particle eigenstates are all exponentially localized,with the form | ψ α ( x ) | ∼ exp (cid:18) − | x − R α | ξ loc (cid:19) , (D5)where R α is the localization center and ξ loc is the secondmoment of the wavefunction. Appendix E: Nandkishore-Potter delocalizationmechanism In this appendix, we briefly comment on a delocaliza-tion mechanism due to resonances mediated by extendedstates in the single-particle spectrum, as exemplified byNandkishore and Potter (NP) in Ref. [36]. NP considerednon-interacting fermion models where the single-particleenergy eigenstates are exponentially localized with a lo-calization length (i.e. the second moment) that scaleswith the single-particle energy E as ξ ( E ) ∼ E − ν for some ν > E → 0. In summary,they show that single-particle orbitals localized far awayfrom each other hybridize at second order in perturbationtheory via a hopping process mediated by the extendedstates at E = 0, showing that delocalization is inevitableif νd > d is the dimension of the systemand the single-particle density of states as E → ρ ( E ) ∼ E Υ .To apply the NP argument to the HamiltonianEq. (A2), note that the (positive part of the) single-particle spectrum of the non-interacting model Eq. (A2)is identical to that of a non-interacting model of spinlessfermions with random hopping strengths [66]. This isevident when is written in terms of Majorana fermions,shown in Eq. (A12). The localization length of the single-particle eigenstates in these models thus scales as − log E ,as can be derived from Eq. (D3) using the Thouless the-orem that relates the single-particle density of states tothe localization length in one-dimensional random hop-ping systems [66, 68]. Thus, since as E → ρ ( E ) ∼ − E (log E ) , ξ ( E ) ∼ − log E (E1)one might naively conclude that Υ = − ν = 0for the present model, which is marginal according to6NP. However, a closer look shows that this model is notdelocalized due to NP.We first reproduce the intuitive version of the NP con-dition [36] with a slight generalization. Using Eq. (E1),the number of orbitals N ( R ) within a distance R that agiven localized orbital can be connected in second-orderperturbation theory via the mediation of extended or-bitals close to E → N ( R ) ∼ (cid:90) E max ( R )0 d E ρ ( E ) ξ ( E ) (E2)where E max ( R ) is the maximum energy for which the localization length ξ ( E ) ≥ R , and the factor of ξ ( E ) inEq. (E2) should be interpreted as the number of orbitalsof energy E that connect to a given localized orbital. Notethat E max ( R ) ∼ exp( − R ) using Eq. (E1), Thus, usingEqs. (E1) and (E2), we obtain N ( R ) ∼ (cid:90) e − R d EE (log E ) ∼ R . (E3)Thus, the number of orbitals that a given orbital connectsto does not increase with increasing R , and a delocaliza-tion at the critical point in the present model is differentfrom the Nandkishore-Potter mechanism. [1] R. Nandkishore and D. A. Huse, “Many-body localiza-tion and thermalization in quantum statistical mechan-ics,” Annu. Rev. Condens. Matter Phys. , 15–38 (2015).[2] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, “Col-loquium: Many-body localization, thermalization, andentanglement,” Rev. Mod. Phys. , 021001 (2019).[3] Philip W Anderson, “Absence of diffusion in certain ran-dom lattices,” Physical Review , 1492 (1958).[4] L. Fleishman and P. W. Anderson, “Interactions and theanderson transition,” Phys. Rev. B , 2366–2377 (1980).[5] DM Basko, IL Aleiner, and BL Altshuler, “Metal–insulator transition in a weakly interacting many-electronsystem with localized single-particle states,” Annals ofphysics , 1126–1205 (2006).[6] Wojciech De Roeck and Fran¸cois Huveneers, “Stabil-ity and instability towards delocalization in many-bodylocalization systems,” Physical Review B , 155129(2017).[7] David J. Luitz, Fran ¸cois Huveneers, and WojciechDe Roeck, “How a small quantum bath can thermal-ize long localized chains,” Phys. Rev. Lett. , 150602(2017).[8] Wojciech De Roeck, Francois Huveneers, Markus M¨uller,and Mauro Schiulaz, “Absence of many-body mobilityedges,” Physical Review B , 014203 (2016).[9] Ionut-Dragos Potirniche, Sumilan Banerjee, and EhudAltman, “Exploration of the stability of many-body lo-calization in d > , 205149 (2019).[10] John Z Imbrie, “On many-body localization for quantumspin chains,” Journal of Statistical Physics , 998–1048(2016).[11] David A Huse, Rahul Nandkishore, Vadim Oganesyan,Arijeet Pal, and SL Sondhi, “Localization-protectedquantum order,” Physical Review B , 014206 (2013).[12] David Pekker, Gil Refael, Ehud Altman, Eugene Demler,and Vadim Oganesyan, “Hilbert-glass transition: Newuniversality of temperature-tuned many-body dynami-cal quantum criticality,” Physical Review X , 011052(2014).[13] Vedika Khemani, Achilleas Lazarides, Roderich Moess-ner, and Shivaji L Sondhi, “Phase structure of drivenquantum systems,” Physical Review Letters , 250401(2016).[14] Dominic V. Else, Bela Bauer, and Chetan Nayak, “Flo-quet time crystals,” Phys. Rev. Lett. , 090402 (2016). [15] Curt W von Keyserlingk, Vedika Khemani, and Shiv-aji L Sondhi, “Absolute stability and spatiotemporallong-range order in floquet systems,” Physical Review B , 085112 (2016).[16] Bela Bauer and Chetan Nayak, “Area laws in a many-body localized state and its implications for topologicalorder,” Journal of Statistical Mechanics: Theory and Ex-periment , P09005 (2013).[17] Anushya Chandran, Vedika Khemani, CR Laumann,and SL Sondhi, “Many-body localization and symmetry-protected topological order,” Physical Review B ,144201 (2014).[18] Yasaman Bahri, Ronen Vosk, Ehud Altman, and AshvinVishwanath, “Localization and topology protected quan-tum coherence at the edge of hot matter,” Nature com-munications (2015).[19] SA Parameswaran, Andrew C Potter, and RomainVasseur, “Eigenstate phase transitions and the emer-gence of universal dynamics in highly excited states,”Annalen der Physik (2017).[20] Arijeet Pal and David A Huse, “Many-body localizationphase transition,” Physical Review B , 174411 (2010).[21] David J Luitz, Nicolas Laflorencie, and Fabien Alet,“Many-body localization edge in the random-field heisen-berg chain,” Physical Review B , 081103 (2015).[22] Jonas A Kj¨all, Jens H Bardarson, and Frank Pollmann,“Many-body localization in a disordered quantum isingchain,” Physical Review Letters , 107204 (2014).[23] Ronen Vosk, David A Huse, and Ehud Altman, “The-ory of the many-body localization transition in one-dimensional systems,” Physical Review X , 031032(2015).[24] Andrew C. Potter, Romain Vasseur, and S. A.Parameswaran, “Universal properties of many-body de-localization transitions,” Phys. Rev. X , 031033 (2015).[25] Thimoth´ee Thiery, Fran ¸cois Huveneers, Markus M¨uller,and Wojciech De Roeck, “Many-body delocalization asa quantum avalanche,” Phys. Rev. Lett. , 140601(2018).[26] Xiongjie Yu, David J. Luitz, and Bryan K. Clark, “Bi-modal entanglement entropy distribution in the many-body localization transition,” Phys. Rev. B , 184202(2016).[27] Vedika Khemani, Say-Peng Lim, DN Sheng, andDavid A Huse, “Critical properties of the many-body localization transition,” Physical Review X , 021013(2017).[28] Anna Goremykina, Romain Vasseur, and Maksym Ser-byn, “Analytically solvable renormalization group forthe many-body localization transition,” Phys. Rev. Lett. , 040601 (2019).[29] Alan Morningstar and David A. Huse, “Renormalization-group study of the many-body localization transition inone dimension,” Phys. Rev. B , 224205 (2019).[30] J. ˇSuntajs, J. Bonˇca, T. Prosen, and L. Vidmar, “Quan-tum chaos challenges many-body localization,” arXive-prints , arXiv:1905.06345 (2019), arXiv:1905.06345[cond-mat.str-el].[31] D. A. Abanin, J. H. Bardarson, G. De Tomasi,S. Gopalakrishnan, V. Khemani, S. A. Parameswaran,F. Pollmann, A. C. Potter, M. Serbyn, and R. Vasseur,“Distinguishing localization from chaos: challenges infinite-size systems,” arXiv e-prints , arXiv:1911.04501(2019), arXiv:1911.04501 [cond-mat.str-el].[32] R. K. Panda, A. Scardicchio, M. Schulz, S. R. Taylor,and M. ˇZnidariˇc, “Can we study the many-body localisa-tion transition?” EPL (Europhysics Letters) , 67003(2020).[33] Piotr Sierant, Dominique Delande, and Jakub Za-krzewski, “Thouless time analysis of anderson and many-body localization transitions,” Phys. Rev. Lett. ,186601 (2020).[34] Daniel S Fisher, “Random antiferromagnetic quantumspin chains,” Physical Review B , 3799 (1994).[35] Daniel S Fisher, “Critical behavior of random transverse-field ising spin chains,” Physical Review B , 6411(1995).[36] Rahul Nandkishore and Andrew C Potter, “Marginalanderson localization and many-body delocalization,”Physical Review B , 195115 (2014).[37] Ronen Vosk and Ehud Altman, “Dynamical quantumphase transitions in random spin chains,” Physical Re-view Letters , 217204 (2014).[38] Kevin Slagle, Yi-Zhuang You, and Cenke Xu, “Dis-ordered xyz spin chain simulations using the spectrumbifurcation renormalization group,” Phys. Rev. B ,014205 (2016).[39] Romain Vasseur, Aaron J. Friedman, S. A.Parameswaran, and Andrew C. Potter, “Particle-hole symmetry, many-body localization, and topologicaledge modes,” Phys. Rev. B , 134207 (2016).[40] Vadim Oganesyan and David A Huse, “Localization ofinteracting fermions at high temperature,” Physical Re-view B , 155111 (2007).[41] YY Atas, E Bogomolny, O Giraud, and G Roux, “Distri-bution of the ratio of consecutive level spacings in randommatrix ensembles,” Physical Review Letters , 084101(2013).[42] Maksym Serbyn, Z Papi´c, and Dmitry A Abanin, “Cri-terion for many-body localization-delocalization phasetransition,” Physical Review X , 041047 (2015).[43] Trithep Devakul and Rajiv RP Singh, “Early breakdownof area-law entanglement at the many-body delocaliza-tion transition,” Physical Review Letters , 187201(2015).[44] Elmer VH Doggen, Frank Schindler, Konstantin STikhonov, Alexander D Mirlin, Titus Neupert, Dmitry GPolyakov, and Igor V Gornyi, “Many-body (de) lo-calization in large quantum chains,” arXiv preprint arXiv:1807.05051 (2018).[45] Gil Refael and Joel E Moore, “Entanglement entropy ofrandom quantum critical points in one dimension,” Phys-ical Review Letters , 260602 (2004).[46] Trithep Devakul, Satya N Majumdar, and David A Huse,“Probability distribution of the entanglement across acut at an infinite-randomness fixed point,” Physical Re-view B , 104204 (2017).[47] Daniel S Fisher and AP Young, “Distributions of gapsand end-to-end correlations in random transverse-fieldising spin chains,” Physical Review B , 9131 (1998).[48] Philip J. D. Crowley and Anushya Chandran, “Avalancheinduced co-existing localised and thermal regions indisordered chains,” arXiv e-prints , arXiv:1910.10812(2019), arXiv:1910.10812 [cond-mat.dis-nn].[49] William Berdanier, Michael Kolodrubetz, S. A.Parameswaran, and Romain Vasseur, “Floquetquantum criticality,” Proceedings of the NationalAcademy of Sciences , 247204(2016).[51] Frank Pollmann, Vedika Khemani, J Ignacio Cirac, andSL Sondhi, “Efficient variational diagonalization of fullymany-body localized hamiltonians,” Physical Review B , 041116 (2016).[52] Xiongjie Yu, David Pekker, and Bryan K Clark, “Find-ing matrix product state representations of highly excitedeigenstates of many-body localized hamiltonians,” Phys-ical Review Letters , 017201 (2017).[53] Trithep Devakul, Vedika Khemani, Frank Pollmann,David A Huse, and SL Sondhi, “Obtaining highly ex-cited eigenstates of the localized xx chain via dmrg-x,”Phil. Trans. R. Soc. A , 20160431 (2017).[54] Frank Schindler, Nicolas Regnault, and Titus Neupert,“Probing many-body localization with neural networks,”Physical Review B , 245134 (2017).[55] Jordan Venderley, Vedika Khemani, and Eun-Ah Kim,“Machine learning out-of-equilibrium phases of matter,”Physical Review Letters , 257204 (2018).[56] Evert PL van Nieuwenburg, Yuval Baum, and GilRefael, “From bloch oscillations to many body local-ization in clean interacting systems,” arXiv preprintarXiv:1808.00471 (2018).[57] R. Sahay, F. Machado, B. Ye, C. R. Laumann, and N. Y.Yao, to appear.[58] Barry M McCoy and Tai Tsun Wu, “Theory of a two-dimensional ising model with random impurities. i. ther-modynamics,” Physical Review , 631 (1968).[59] AP Young and H Rieger, “Numerical study of the randomtransverse-field ising spin chain,” Physical Review B ,8486 (1996).[60] Ferenc Igl´oi and Heiko Rieger, “Random transverse isingspin chain and random walks,” Physical Review B ,11404 (1998).[61] Olexei Motrunich, Siun-Chuon Mau, David A Huse, andDaniel S Fisher, “Infinite-randomness quantum ising crit-ical fixed points,” Physical Review B , 1160 (2000).[62] Ross H McKenzie, “Exact results for quantum phasetransitions in random xy spin chains,” Physical ReviewLetters , 4804 (1996). [63] Wade DeGottardi, Diptiman Sen, and Smitha Vishvesh-wara, “Majorana fermions in superconducting 1d sys-tems having periodic, quasiperiodic, and disordered po-tentials,” Physical Review Letters , 146404 (2013).[64] Gil Refael and Ehud Altman, “Strong disorder renormal-ization group primer and the superfluid–insulator transi-tion,” Comptes Rendus Physique , 725–739 (2013).[65] TP Eggarter and R Riedinger, “Singular behavior oftight-binding chains with off-diagonal disorder,” Physi-cal Review B , 569 (1978).[66] Chenggang Zhou and RN Bhatt, “One-dimensional chain with random long-range hopping,” Physical Review B ,045101 (2003).[67] Akshay Krishna and R. N. Bhatt, “Beyond univer-sal behavior in the one-dimensional chain with ran-dom nearest-neighbor hopping,” Physical Review B (2020), 10.1103/physrevb.101.224203.[68] D J Thouless, “A relation between the density of statesand range of localization for one dimensional random sys-tems,” Journal of Physics C: Solid State Physics5