Quasi-Many-Body Localization of Interacting Fermions with Long-Range Couplings
LLocalisation of Interacting Power-Law Random Banded Fermions
S. J. Thomson
1, 2, ∗ and M. Schir´o † Centre de Physique Th´eorique, CNRS, Institut Polytechnique de Paris, Route de Saclay, F-91128 Palaiseau, France Institut de Physique Th´eorique, Universit´e Paris-Saclay, CNRS, CEA, F-91191 Gif-sur-Yvette, France JEIP, USR 3573 CNRS, Coll´ege de France, PSL Research University,11 Place Marcelin Berthelot, 75321 Paris Cedex 05, France (Dated: May 1, 2020)We study a model of one-dimensional fermions with random long-range hoppings and interaction,both Gaussian-distributed with power-law widths which decay with different exponents, and on-site disorder. This model realizes an interacting many-body version of the well-studied power-lawbanded random matrix model. Using a truncated flow equation approach, we study static anddynamical properties as a function of hopping and interaction exponents. We show that, at largeon-site disorder and for short-range interactions, a transition from a delocalised phase to quasi many-body-localised (MBL) behavior exists upon decreasing the hopping range. This quasi-MBL phase ischaracterized by intriguing properties such as algebraically decaying l -bit interactions. Surprisinglywe find that a crossover survives, albeit broadened, upon increasing the range of the interactions. Introduction -
The combination of disorder and inter-actions in many-body quantum systems can lead to a richvariety of physics far from equilibrium, from many-bodylocalisation (MBL) [1–5] to quantum glasses [6–11]. MBLis by now rather well understood for one-dimensionalmodels with short-range interactions, where a set of mu-tually commuting local integrals of motion (LIOMs, or l -bits) can be identified [12–15], while its fate in higherdimensions or in the presence of long-range couplings isless well established.Theoretical investigations of localisation in long-rangesystems date back to Anderson’s original work [16]. Onewell-understood example is the non-interacting randomhopping problem, where the hopping terms decay as apower-law with exponent α , also known as Power-LawRandom Banded Matrix (PRBM) model. In this case,localisation is destroyed for α < d (where d is the spa-tial dimension) and the system is critical at α = d [17–25]. In the interacting many-body case, the the effect oflong-range couplings are certainly less well-understood.Most studies have focused on quantum spin models withpower-law decaying exchange couplings of random signs.Estimates based on the locator expansion and its break-down suggest an instability of the (many-body) localisedphase for slowly decaying transverse exchange with ex-ponent β < d [26–28], independently of the longitudi-nal exponent α which controls the degrees of freedominvolved in resonance formation [27, 29]. The robustnessand generality of those perturbative arguments howeverhas not been fully discussed. In particular, convergenceof the locator expansion provides at most a sufficient con-dition for localisation but does not usually guarantee de-localisation. Different scenarios have emerged recentlywhich are consistent with localised behavior even in pres-ence of slowly decaying power-law interactions, for whichthe locator expansion does not converge. Examples in-clude order-enabled localisation [30] cooperative shield-ing [31–33] or correlation-induced localisation in single particle problems [34, 35]. A better theoretical under-standing of disordered quantum many-body systems withpower-law couplings is therefore much needed. Exact di-agonalisation studies, which played a crucial role in un-derstanding short-ranged MBL, are limited to small sizesand suffer from strong finite size effects in long-rangemodels, therefore alternative methods are highly desir-able. This is particularly pressing given the experimentalrelevance of long-ranged interacting quantum simulators,from dipolar systems to trapped ions [36–39] which pro-vides evidence of MBL phenomenology.In this Letter we investigate a model of one-dimensional fermions with diagonal disorder, randomlong-range hopping and random long-range interactionswhich decay like power-law functions of distance, withexponents α and β respectively. Using a truncated flowequation approach able to describe both the short-rangedMBL phase [40] and the delocalization of non-interactingfermions with power-law hopping [41], we show thatfor rapidly decaying power laws ( α, β >
2) the sys-tem at large on-site disorder is in a quasi-MBL phase[70] characterized by algebraically decaying l -bit interac-tions [33, 42] that we explicitly construct. This quasi-MBL phase turns delocalised upon decreasing α . In-terestingly, this transition survives upon decreasing therange of the interactions β . The model -
We start from a Hamiltonian describing aone-dimensional chain of interacting fermions given by: H = X i h i n i + 12 X ij V ij n i n j + X ij J ij c † i c j (1)where the on-site disorder is drawn from a box distri-bution h i ∈ [0 , W ]. The couplings J ij and V ij arealso random and drawn from Gaussian distributions withstandard deviations which decay with distance as σ J = J / | i − j | α and σ V = V / | i − j | β respectively. We fix J = 0 . V = 0 . W = 5, such that the modelwith short-ranged hopping and interactions (respectively a r X i v : . [ c ond - m a t . d i s - nn ] A p r α = β = ∞ ) would be in the MBL phase, and vary thepower-law exponents α and β only.In absence of interactions, the model reduces to thePower-Law Random Banded Matrix (PRBM) model,previously studied in Refs. [17–25]. In this case, it isknown that long-range hoppings favour the delocalisedphase and give rise to a delocalization transition evenin d = 1, for α ≤
1. In the following we discuss whathappens to the PRBM model in presence of interactionsand study the interplay/competition between power-lawhoppings and power-law interactions. We first considerthe two effects separately, fixing α = ∞ and varying β and vice versa, while later we present a complete phasediagram in the ( α, β ) plane. Method -
We attack the problem using the flow equa-tion approach [43–53] which we have recently used tostudy MBL in the short ranged case [40] as well as thenon-interacting PRBM model [41]. The main idea is todiagonalise the Hamiltonian through a series of infinites-imal unitary transforms parametrised by a fictitious ‘flowtime’ l which runs from l = 0 (initial basis) to l → ∞ (diagonal basis). The Hamiltonian flow readsd H d l = [ η ( l ) , H ( l )] . (2)where η ( l ) is the generator of the flow and the initialcondition at l = 0 is given by the Hamiltonian in Eq. (1).While for quadratic problems the approach is exact, inpresence of interactions the flow generates higher-ordercouplings not present in the original microscopic model.We use here a truncated scheme, originally introduced inRef [40], that captures the essential physics of the MBLphase. We make an ansatz for the form of the runningHamiltonian H ( l ) = H ( l ) + V ( l ), with H ( l ) = X i h i ( l ) : c † i c i : + 12 X ij ∆ ij ( l ) : c † i c i c † j c j + V ( l ) = X ij J ij ( l ) : c † i c j : , (3)which allows us to close the hierarchy of flow equations,discarding all newly generated higher-order terms. InEq. (3) the : O : notation signifies normal-ordering [71].Given the above ansatz, here we use Wegner’s choice forthe generator and choose it to be the commutator ofthe diagonal and off-diagonal parts of the Hamiltonian, η ( l ) = [ H ( l ) , V ( l )]. The flow equations for the runningcouplings can be read off from d H / d l = [ η ( l ) , H ( l )] [54].This choice of generator, although not unique [41, 48, 55],guarantees [43, 56] that the off-diagonal terms vanish inthe l → ∞ limit; other choices of generator are also pos-sible [41, 48, 55]. We can quantify the accuracy of ourtruncation by computing quantities which are preservedby unitary transforms, such as traces of integer powers ofthe Hamiltonian. By computing how precisely they are preserved by this ansatz we can get a measure for theerror in this truncation scheme - see [54] for details. Decay of l -bit interactions and real-space support - In the l → ∞ limit, the off-diagonal terms will van-ish and we will obtain a diagonal Hamiltonian given by˜ H = P i ˜ h i n i + P ij ˜∆ ij n i n j . In all of the following, thetilde notation indicates quantities in the l → ∞ diago-nal l-bit basis. By computing the effective Hamiltonian,we can extract two different aspects of its long-range be-haviour. First, we extract the distance-dependence of thecoefficients ˜∆ ij . These coefficients, which decay expo-nentially in short-range systems [40, 57, 58], are stronglymodified by the existence of long-range couplings.We can also compute the real-space support of the l -bitoperators directly. Starting from a local density operator˜ n i defined in the diagonal l → ∞ basis with support onlyon a single site, we can transform it back into the physical(i.e. real space) basis by inverting the unitary transformused to diagonalise the Hamiltonian. This operator canbe transformed according to:d˜ n i ( l )d l = [ η ( l ) , ˜ n i ( l )] (4)here starting with l = ∞ and flowing backwards to l = 0.To parameterise the flow of this operator, we make thefollowing ansatz for the running number operator:˜ n i ( l ) = X j A ( i ) j ( l ) n j + X jk B ( i ) jk ( l ) c † j c k (5)Note that higher-order terms cannot be consistently in-cluded at this order of the truncation scheme. Thenormal-ordering procedure employed as part of this con-struction [54] does, however, allow us to take into accountthe leading effects of the interactions even at this order.In Fig. 1, we show these quantities in the case of power-law hopping and nearest-neighbor interactions (corre-sponding to β = ∞ ). The ∆ ij retain their exponentially-decaying nature at short distances, but acquire power-law tails at long range, with a decay exponent ζ ≈ α for α ≥
1. This follows immediately from the structure ofthe eigenstates of the PRBM problem, which are indeedexponentially localised at short distance with power-lawtails [19]. The real-space support of the l -bits also showpower-law tails characteristic of delocalisation, after aninitial exponential decay at short range. The precise dis-tance where the decay crosses from exponential to power-law depends on the exponent, as well as both the disorderand interaction strength. As α → ∞ , the real-space sup-port of the l -bits decays exponentially over a larger rangebefore the power-law tail appears, and the resulting l -bitsclosely match the nearest-neighbour case (black dashedline). In Fig. 2, we show the case of power-law inter-actions and nearest-neighbor hopping (corresponding to α = ∞ ). The ∆ ij retain their initial power-law distri-bution at all distances and at all stages during the flow FIG. 1: l -bit interactions (top) and real-space support (bot-tom) for power-law hopping and nearest-neighbor interactions( α ∈ [0 . , .
0] (top to bottom) in increments of 0 .
25, and β = ∞ ). a) The disorder-averaged (median) ˜∆ ij decay as apower-law at long distances (notice log-log scale, dashed lineis a power law guide to eye) and as an exponential at shortdistances (see inset, semi-log scale, for α ∈ [3 . , l -bits exhibit an exponential decay (most visible for large α )crossing over to an extended behavior with long power-lawtails. The dashed line is the ( α → ∞ , β → ∞ ) short-rangelimit. Chain size L = 128, disorder realisations N s = 256. procedure. Surprisingly, we find that the real-space sup-port of the l -bits is essentially unmodifed by the range ofthe interactions: they retain their exponentially decay-ing character even in the limit of β = 0, with only anextremely small extended ‘tail’ appearing following thestrong initial exponential decay. This may be an effect ofthe truncation in Eq. 3 suppressing degrees of freedomresponsible for delocalisation, or it may be that delocal-isation is only seen in (weak) higher-order contributionsto Eq. 5, corresponding to multipole processes. Dynamics of Imbalance and Phase Diagram -
We nowmove on to study the effect of power-law couplings on thequantum dynamics of the system. This can be accessedusing the FE approach [40, 54]. We set up an initialcharge density wave (CDW) state and see how it relaxesunder its own quantum dynamics. To monitor this, wedefine the imbalance as: I ( t ) = 2 L X i ( − i h n i ( t ) i (6)The long time behavior of the imbalance is often usedas a proxy for the MBL transition, since in a localisedphase any initial inhomogeneity persists at long time dueto enhanced memory of initial conditions while in a ther-mal, delocalised phase the imbalance is expected to decayto zero as a power law with a disorder-dependent expo-nent, vanishing at the transition [59, 60]. Using the time-dependent mean-field decoupling on the effective l -bitHamiltonian, the results for the relaxation dynamics of FIG. 2: l -bit interactions (top) and real-space support (bot-tom) for nearest-neighbor hopping and power-law interaction( α = ∞ , β ∈ [0 . , .
0] (top to bottom) in increments of 0 . ij retain their initial power-lawdistribution for all β , except at very short distance and large β (see inset, semi-log scale, for β ∈ [3 . , l -bits re-main exponentially localised in real space, with no almost nodependence on β . The dashed line is the same quantity for ashort ranged many-body localised model (( α → ∞ , β → ∞ )).Chain size L = 128, disorder realisations N s = 256. the imbalance are shown in Figure 3, for chains of length L = 64 in the cases of power-law hopping with nearest-neighbour interactions (panel a), and nearest-neighbourhopping with power-law interactions (panel b). In Fig.3a, we see that for α & α theimbalance continuously decrease toward zero, a behav-ior that is reminscent of the PRBM model. For α = 0,the decay of the imbalance is approximately exponential,while for α > β , i.e. makingthe range of interactions larger, has little to no effect onthe long-time imbalance and the system remains local-ized, with small values of β leading to the appearance ofa short plateau that vanishes at longer times.Having examined their effects separately, we now com-pute the imbalance in the presence of both long-rangedinteractions and long-range hopping, and obtain a qual-itative phase diagram shown in Fig. 4 where we showthe imbalance I ( t ) at a time t ∗ = 100 after the quenchas a function of α, β and super-impose lines at fixed im-balance as guide to the eye. In the upper-right corner,corresponding to fast decaying hopping and interactions( α, β ≥ β ≥ α the imbalance displaysa sharp crossover from localized to delocalized behavior,the interacting many-body analog of the PRBM previ-ously studied with FE techniques in Ref. [41]. We can FIG. 3: Relaxation of the imbalance following a quench from aCDW state with a) power-law hopping α ∈ [0 . , .
5] in incre-ments of 0 .
25 (bottom to top) and β = ∞ (nearest-neighbourinteractions) and b) power-law interactions β ∈ [0 . , .
5] inincrements of 0 .
25 (top to bottom) with α = ∞ (nearest-neighbour hopping). Decreasing α makes the long-time im-balance go to zero (as a power law in time for small α , see topinset) whereas changing β has almost no effect on the long-time dynamics of the imbalance which approaches a finiteplateau almost exponentially (see inset). Chain size L = 64,disorder realisations N s = 256. now ask what happens to those two phases as we increasethe range of the interaction, i.e. decreases β toward zero.The ergodic phase is expected to be robust to long-rangeinteractions, and indeed we see that the imbalance for α < β (see the almost vertical contour lines) . On theother hand, and quite surprisingly, we find the imbalanceto remain strongly unaffected by long range interactionseven for α &
1, consistently with the results of Figure 3for the α = ∞ case. However the lines at fixed imbalancebends toward right for small β , suggesting that the lowerright corner of the phase diagram is less localised thanthe upper right corner. Discussion -
Our results show that upon increasing therange of the hopping, a transition from delocalisation toquasi-MBL exists, both for short ranged interactions aswell as for β <
2, in a regime where perturbative ar-guments based on locator expansion would exclude it.While one could suspect of an artefact of our truncationscheme, we have performed extensive checks to validateour approach in this regime, including comparison withexact numerics for small system sizes and monitoring theflow invariant, a sensitive probe of the validity of ourscheme (see [54] for details). Still we cannot exclude thatcertain processes leading to delocalisation are not cap-tured by our scheme. This quasi-MBL phase could alsobe metastable for finite size and finite time: recent workssuggest that in the intermediate regime 1 < β <
FIG. 4: Qualitative phase diagram of model (1) as a functionof α (hopping exponent) and β (interaction exponent). Thecolour scale shows the imbalance I ( t ) at a time t ∗ = 100following a quench. The dotted lines show contours of theimbalance I ( t ∗ ) = 0 . , . , . , .
45: the solid white linesare guides to the eye. The system size is L = 64, with 50 ≤ N s ≤
128 disorder realisations, as required for convergence. increasing system size L (or equivalently, exhibit a size-dependent critical disorder W c ( L )) [26, 28, 61, 62]. Ourresults show [54] that the quasi-MBL phase shrinks asthe system size is increased. The Gaussian distributionof couplings could also play a role, as long-range randomsign couplings, commonly studied in quantum spin mod-els, exhibit enhanced delocalization [54]. Finally, it isworth noticing that in the α, β → α and β = 0 induces a finite rangehopping and interaction which could lead, in addition tothe random onsite field, to an increased localised behav-ior [64], at least for finite systems. Conclusion -
We have used the flow equation methodto study a model of one dimensional fermions with gaus-sian distributed, power-law decaying, hoppings and inter-actions and diagonal box disorder, an interacting manybody version of the celebrated PRBM model. For largediagonal disorder, compared to typical scales of inter-actions and hoppings, we have provided evidence of atransition from a delocalised ergodic phase to a quasiMBL phase upon increasing the exponent α controllingthe range of hopping. Such a crossover survives even forslowly decaying interactions, β <
2, although it appearsto become less sharp. This quasi-MBL phase, althoughpossibly metastable, appears to have intriguing proper-ties such as algebraically decaying l -bit interactions. Anopen question is the stability of such a phase to the prop-agation of ergodic bubbles: avalanche arguments sug-gest that MBL should disappear in any dimension forinteractions decaying slower than exponential [65] andour model and approach could provide insights into thislargely unexplored question, for example by studying thecoupling of this quasi-MBL phase to an ergodic bath [66]. Acknowledgements -
The computations were per-formed on the Coll´ege de France IPH computer cluster.We acknowledge use of the QuSpin exact diagonalisationlibrary for benchmarking our FE code [67, 68], and sup-port from the grant DynDisQ from DIM SIRTEQ. ∗ Electronic address: [email protected] † Electronic address: [email protected]; On Leavefrom: Institut de Physique Th´eorique, Universit´e ParisSaclay, CNRS, CEA, F-91191 Gif-sur-Yvette, France[1] D. M. Basko, I. L. Aleiner, and B. L. Altshuler, Annalsof physics , 1126 (2006).[2] E. Altman and R. Vosk, Annu. Rev. Condens. MatterPhys. , 383 (2015).[3] R. Nandkishore and D. A. Huse, Annu. Rev. Condens.Matter Phys. , 15 (2015).[4] F. Alet and N. Laflorencie, Comptes Rendus Physique , 498 (2018).[5] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev.Mod. Phys. , 021001 (2019).[6] M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S.Fisher, Phys. Rev. B , 546 (1989).[7] L. F. Cugliandolo and G. Lozano, Phys. Rev. Lett. ,4979 (1998).[8] G. Biroli and L. F. Cugliandolo, Phys. Rev. B , 014206(2001).[9] G. Biroli and O. Parcollet, Phys. Rev. B , 094414(2002).[10] L. Berthier and G. Biroli, Rev. Mod. Phys. , 587(2011).[11] S. J. Thomson, P. Urbani, and M. Schiro, arXiv e-printsarXiv:1904.03147 (2019), 1904.03147.[12] M. Serbyn, Z. Papi´c, and D. A. Abanin, Phys. Rev. Lett. , 127201 (2013).[13] D. A. Huse, R. Nandkishore, and V. Oganesyan, Phys.Rev. B , 174202 (2014).[14] J. Z. Imbrie, Journal of Statistical Physics , 998(2016), ISSN 1572-9613.[15] J. Z. Imbrie, V. Ros, and A. Scardicchio, Annalen derPhysik , 1600278 (2017).[16] P. W. Anderson, Physical review , 1492 (1958).[17] C. Yeung and Y. Oono, EPL (Europhysics Letters) ,1061 (1987).[18] L. S. Levitov, Physical review letters , 547 (1990).[19] A. D. Mirlin, Y. V. Fyodorov, F.-M. Dittes, J. Quezada,and T. H. Seligman, Physical Review E , 3221 (1996).[20] L. Levitov, Annalen der Physik , 697 (1999).[21] I. Varga and D. Braun, Physical Review B , R11859(2000).[22] A. D. Mirlin and F. Evers, Physical Review B , 7920(2000).[23] F. Evers and A. D. Mirlin, Physical Review Letters ,3690 (2000).[24] V. E. Kravtsov, O. Yevtushenko, and E. Cuevas, Journalof Physics A: Mathematical and General , 2021 (2006). [25] F. Evers and A. D. Mirlin, Reviews of Modern Physics , 1355 (2008).[26] A. L. Burin (2006), cond-mat/0611387.[27] N. Y. Yao, C. R. Laumann, S. Gopalakrishnan, M. Knap,M. Mueller, E. A. Demler, and M. D. Lukin, Physicalreview letters , 243002 (2014).[28] A. L. Burin, Physical Review B , 104428 (2015).[29] D. B. Gutman, I. V. Protopopov, A. L. Burin, I. V.Gornyi, R. A. Santos, and A. D. Mirlin, Physical ReviewB , 245427 (2016).[30] R. M. Nandkishore and S. L. Sondhi, Physical Review X , 041021 (2017).[31] L. F. Santos, F. Borgonovi, and G. L. Celardo, Physicalreview letters , 250402 (2016).[32] S. Roy and D. E. Logan, SciPost Phys. , 42 (2019).[33] X. Deng, G. Masella, G. Pupillo, and L. Santos (2019),1912.08131.[34] X. Deng, V. E. Kravtsov, G. V. Shlyapnikov, and L. San-tos, Phys. Rev. Lett. , 110602 (2018).[35] P. A. Nosov, I. M. Khaymovich, and V. E. Kravtsov,Phys. Rev. B , 104203 (2019).[36] G. A. ´Alvarez, D. Suter, and R. Kaiser, Science , 846(2015), ISSN 0036-8075.[37] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess,P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, NaturePhysics , 907 (2016).[38] J. Zeiher, J.-y. Choi, A. Rubio-Abadal, T. Pohl, R. vanBijnen, I. Bloch, and C. Gross, Phys. Rev. X , 041063(2017).[39] S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya,F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, et al.,Nature , 221 (2017).[40] S. J. Thomson and M. Schir´o, Physical Review B ,060201(R) (2018).[41] S. J. Thomson and M. Schir´o, Eur. Phy. J. B (2020).[42] G. De Tomasi, Phys. Rev. B , 054204 (2019).[43] S. Kehrein, The flow equation approach to many-particlesystems , vol. 217 (Springer, 2007).[44] M. Moeckel and S. Kehrein, Phys. Rev. Lett. , 175702(2008).[45] A. Hackl and S. Kehrein, Phys. Rev. B , 092303 (2008).[46] A. Hackl and S. Kehrein, Journal of Physics: CondensedMatter , 015601 (2009).[47] M. Eckstein, A. Hackl, S. Kehrein, M. Kollar,M. Moeckel, P. Werner, and F. Wolf, The EuropeanPhysical Journal Special Topics , 217 (2009).[48] C. Monthus, Journal of Physics A: Mathematical andTheoretical , 305002 (2016).[49] V. L. Quito, P. Titum, D. Pekker, and G. Refael, Phys.Rev. B , 104202 (2016).[50] D. Pekker, B. K. Clark, V. Oganesyan, and G. Refael,Physical Review Letters , 075701 (2017).[51] S. Savitz, C. Peng, and G. Refael, Phys. Rev. B ,094201 (2019).[52] X. You, D. Pekker, and B. K. Clark (2019),arXiv:1909.11097.[53] S. P. Kelly, R. Nandkishore, and J. Marino, NuclearPhysics B , 114886 (2020).[54] See Supplementary Material for more details on the flowequations and for comparison with exact diagonalization.[55] S. Savitz and G. Refael, Physical Review B , 115129(2017).[56] F. Wegner, Annalen der physik , 77 (1994). [57] L. Rademaker and M. Ortu˜no, Phys. Rev. Lett. ,010404 (2016).[58] L. Rademaker, M. Ortu˜no, and A. M. Somoza, Annalender Physik pp. 1600322–n/a (2017), ISSN 1521-3889,1600322.[59] D. J. Luitz, N. Laflorencie, and F. Alet, Phys. Rev. B , 060201(R) (2016).[60] G. Biroli and M. Tarzia, Phys. Rev. B , 201114(R)(2017).[61] K. S. Tikhonov and A. D. Mirlin, Physical Review B ,214205 (2018).[62] S. Gopalakrishnan and D. A. Huse, Physical Review B , 134305 (2019).[63] S. Sachdev and J. Ye, Phys. Rev. Lett. , 3339 (1993).[64] A. M. Garc´ıa-Garc´ıa and M. Tezuka, Physical Review B , 054202 (2019).[65] W. De Roeck and F. Huveneers, Phys. Rev. B , 155129(2017). [66] P. J. D. Crowley and A. Chandran, arXiv e-printsarXiv:1910.10812 (2019), 1910.10812.[67] P. Weinberg and M. Bukov, SciPost Phys. , 003 (2017).[68] P. Weinberg and M. Bukov, SciPost Phys. , 20 (2019).[69] F. Wegner, Journal of Physics A: Mathematical and Gen-eral , 8221 (2006).[70] We use the term “quasi-MBL” because we cannot ruleout the disappearance of localisation at longer times, orat larger system sizes, due to the effects of avalancheinstabilities - see Ref. [65] and the Discussion at the endof our manuscript for more details.[71] We adopt normal ordering using the : ˆ O : notation in or-der to i) ensure a consistent ordering of operators whencomputing commutation relations, and ii) efficiently re-sum contributions from higher-order terms to turn theflow equation method into a non-perturbative scheme -see Refs. [43, 69] and the Supplementary Material fordetails. upplementary Material to “Localisation of Interacting Power-Law Random BandedFermions” S. J. Thomson
1, 2, ∗ and M. Schir´o
2, 3, † Centre de Physique Th´eorique, CNRS, Institut Polytechnique de Paris, Route de Saclay, F-91128 Palaiseau, France Institut de Physique Th´eorique, Universit´e Paris-Saclay, CNRS, CEA, F-91191 Gif-sur-Yvette, France JEIP, USR 3573 CNRS, Coll´ege de France, PSL Research University,11 Place Marcelin Berthelot, 75321 Paris Cedex 05, France (Dated: May 1, 2020)
Normal-ordering
A key ingredient in the calculation is the adoption of a normal-ordering procedure [1–3], which allow us to con-sistently group together terms at each order of the Hamiltonian, and to incorporate corrections from higher-orderterms which are then discarded from our variational manifold. We will assume all contractions will be computed withrespect to a product state, and the relevant contractions will be denoted: { c † i , c j } = G ij + ˜ G ji = δ ij (1) G ij = h c † i c j i = δ ij h n i i (2)˜ G ji = h c j c † i i = δ ij − h c † i c j i = δ ij (1 − h n i i ) (3)To calculate the commutators of normal-ordered strings of operators, we need to use the following theorem [2]:: O ( A ) :: O ( A ) : =: exp X ij G ij ∂ ∂A j ∂A i O ( A ) O ( A ) : (4)which, for example, leads to the following commutation relation for pairs of fermion operators:[: c † α c β : , : c † γ c δ :] = ( G γβ + ˜ G βγ ) : c † α c δ : − ( G αδ + ˜ G δα ) : c † γ c β : +( G αδ ˜ G βγ − G ) γβ ˜ G δα ) (5)= δ βγ : c † α c δ : − δ αδ : c † γ c β : +( G αδ ˜ G βγ − G γβ ˜ G δα ) (6)which is just the regular commutator plus a constant. All necessary commutators can be computed from Eq. 4, thoughthe calculation is extremely tedious and will not be shown here: for further details, see Refs. [1–3]. In principle, oneshould define an l -dependent state and recompute the normal-ordering corrections at each flow timestep accordingly,however to capture the main physics it is sufficient to simply pick a target state and compute the corrections withrespect to that state [2]. In the main text, we compute the contractions with respect to an infinite-temperatureproduct state such that h n i i = 0 . ∀ i . This has the advantage that many of the normal-ordering corrections (e.g. thefinal terms in Eq. 6 above) vanish identically. Generator
The generator η ( l ) is scale-dependent and changes throughout the flow. Here we use Wegner’s choice for thegenerator such that it is given by the commutator of the diagonal and off-diagonal parts of the Hamiltonian, η ( l ) =[ H ( l ) , V ( l )]. This choice guarantees [2, 4] that the off-diagonal terms vanish in the l → ∞ limit; other choices ofgenerator with different convergence properties are also possible [5–7]. WIth this choice, the generator is given by: η = X ij F ij : c † i c j : + X ijk ζ kij : c † k c k c † i c j : (7)with F ij ≡ J ij [( h i − h j ) − ∆ ij ( h n i i − h n j i )] and ζ kij ≡ J ij (∆ ik − ∆ jk ), where the scale-dependence of the coefficientshas been suppressed for clarity, The interaction term will lead to the generation of new higher-order terms in theHamiltonian during the flow. In practice, the successive generation of these higher-order terms quickly renders thecalculation analytically intractable, so we make the truncation specified in the main text and discard all newlygenerated terms outside of this variational manifold. a r X i v : . [ c ond - m a t . d i s - nn ] A p r Flow Equations
The flow of the Hamiltonian coefficients can be read off from d H / d l = [ η ( l ) , H ( l )], following a lengthy calculation.Explicit expressions for the flow equations are as follows:d h i ( l )d l = 2 X j J ij ( h i − h j ) − X j J ij ∆ ij ( h n i i − h n j i ) + X jk J jk (∆ ik − ∆ ij )( h n k i − h n j i ) (8)d J ij ( l )d l = − J ij ( h i − h j ) − X k J ik J kj (2 h k − h i − h j ) + 2 J ij ∆ ij ( h i − h j )( h n i i − h n j i ) − J ij ∆ ij ( h n i i + h n j i − h n i ih n j i ) − X k J ij (∆ ik − ∆ jk ) h n k i (1 − h n k i )+ X k J ik J kj [(∆ ij − jk )( h n j i − h n k i ) + (∆ ij − ik )( h n i i − h n k i )] (9)d∆ ij ( l )d l = 2 X k = i,j (cid:2) J ik (∆ ij − ∆ kj ) + J jk (∆ ij − ∆ ik ) (cid:3) (10)We numerically integrate these equations until the off-diagonal elements have decayed to the required accuracy,typically using l max ≈ , discarding couplings which have reached zero below some cutoff (typically 10 − or less).Note that in the low-disorder, large-interaction limit where the system may become delocalised, the neglected higher-order terms may carry significant weight. In this case, the normal-ordering corrections due to these neglected higher-order terms can become large enough to overpower the dominant Wegner decay terms, i.e. the − J ij ( h i − h j ) term inEq. 9. In order to obtain physically reasonable flows in this regime, one must neglect the normal-ordering correctionswhenever they threaten to overwhelm the dominant Wegner rules. Retaining them will result in an unphysical anduncontrolled flow. The normal-ordering corrections should always be sub-leading to the main terms in the flowequations. Also note that in the same low-disorder, large-interaction regime, Eq. 10 can exhibit spurious divergenceswhich must be handled carefully in order to obtain physically reasonable results. We emphasise, however, that whenthese problems are encountered, the validity of the ansatz will likely have already broken down. Accuracy: Eigenvalue Comparison with Exact Diagonalisation
In order to benchmark the accuracy of our results, we compared the static properties (i.e. the eigenvalues) withExact Diagonalisation (ED) results obtained using the QuSpin package [8, 9]. We define the averaged relative erroras: δε = 1 N N X i | ε F Ei − ε EDi | ε EDi (11)where the sum runs over states in the many-body Hilbert space. We can compute this quantity, here restrictingourselves to the half-filled states in the centre of the spectrum, for a variety of power-law exponents α and β in orderto benchmark the accuracy of our results. The results are summarised in Fig. 1, where we show the average relativeerror across the phase diagram shown in Fig. 4 of the main text, here for a system size of L = 12 and with N s = 512disorder realisations. We also verified that the error decreases rapidly with increasing disorder strength, as expected,though we do not reproduce this data here.We do not attempt to compute level spacing distributions in this work - while it is possible to resolve exponentiallysmall level spacings using flow equation methods (see Ref. [7]), it requires integration to extremely large flow times,and in any case the truncation of the running Hamiltonian to only a polynomial number of couplings likely preventsus from resolving exponentially small features, as previously found in Ref [3] for a short-range MBL system. Accuracy: Invariants of the Flow
As with any other unitary transform, there are a variety of conserved quantities of the flow equation formalism.Specifically, traces of integer powers of the Hamiltonian I p = Tr[ H p ] are commonly known as ‘invariants of the flow’, FIG. 1: The logarithm of the disorder-averaged relative error δε plotted across the same parameter values as the phase diagramin Fig. 4 of the main text, averaged over N s = 512 disorder realisations. The error is largest in the case where all couplingsare both long-range, and decreased sharply when either or both exponents have a value greater than zero. Note that and are preserved by an exact implementation of the flow equation formalism. As we have seen, however, in order forthe calculation to remain tractable we must make an approximation for the running Hamiltonian of the system. Theneglect of any terms not contained within the ansatz Hamiltonian introduces an error: this error may be quantifiedby computing the invariants of the flow at the start and end of the procedure, and then computing the differencebetween them. This difference is zero if the unitary transform is exact, and non-zero if the truncation has introducedan error. Here we focus on the second invariant ( p = 2) and define the truncation error as: δI = | I ( l = 0) − I ( l = ∞ ) | ( I ( l = 0) + I ( l = ∞ )) (12)The main source of error in this scheme is the strength of the interactions, which contribute to the generation ofhigher-order terms not included in our variational manifold. In the present case, as the truncated higher-order termsscale approximately with integer powers of the interaction strength ∆ (cid:28)
1, the high order terms are typically smalland the accuracy very good. For example, in the results shown in Figs. 1 and 2 of the main text, δI ≤ .
015 at alltimes. However, in the limit of β →
0, there are a large number of interaction terms and the neglected terms canbegin to become significant. To get an idea of the accuracy of our results, we can compute this quantity across thephase diagram in Fig. 4 of the main text: the result is shown in Fig. 2. We find that the transform is almost perfectlyunitary across the entire phase diagram, with the main deviations away from unitarity occuring close to β = 0. Dynamics
We can compute the dynamics of any operator by transforming the operator into the basis which diagonalises theHamiltonian, time-evolving with respect to the diagonal Hamiltonian, and then flowing the operator back into thephysical basis. We demonstrate this with the number operator n i ( t ). We make the following ansatz for the flow ofthis operator: n i ( l, t = 0) = X j A ( i ) j ( l ) : c † j c j : + X jk B ( i ) jk ( l ) : c † j c k : (13)where the coefficients are explicit functions of both the fictitious flow time l and the real time t . FIG. 2: Behaviour of the flow invariant across the phase diagram, using the same data as Fig. 4 of the main text, with L = 64.The flow invariant is maximal for β = 0. Note that the colour scale shows the logarithm of δI : the deviation of the flowequation transform from perfect unitarity is less than one percent across the phase diagram. We remind that reader that eachof the 11 ×
11 points in this phase diagram is the result of 50 ≤ N s ≤
128 disorder realisations, as required for convergence.The reader may also note that for β = 0 . , .
25 and 2 .
0, we took additional data points at double the resolution along the α axis in order to ensure that our resolution was sufficient to resolve the main features in the phase diagram. The flow equations for this operator can be obtained by computing n i ( l ) = [ η ( l ) , n i ] and are given by:d A ij d l = − X k J jk ( h k − h j ) B kj , (14)d B jk d l = − J jk ( h k − h j )( A ik − A ij ) − X n [ J nj ( h n − h j ) B nk + J nk ( h n − h k ) B nj ] , (15)After transforming n i ( t = 0) into the diagonal basis, we can time-evolve it with respect to the diagonal Hamiltonian˜ H = P k ˜ h k ˜ n k + P jk ˜∆ jk ˜ n j ˜ n k by solving the Heisenberg equation of motion. As the Hamiltonian is still interact-ing, despite being diagonal, a closed form cannot be obtained. To close the equation, we choose a time-dependentdecoupling of the interaction term, and the time-evolved operator is given by:˜ n i ( l = ∞ , t ) = X j A ( i ) j ( l ) n j + X jk B ( i ) jk ( l )e iφ jk ( t ) c † j c k (16) φ jk ( t ) = Z t d t " (˜ h k − ˜ h j ) + X m ( ˜∆ km − ˜∆ jm ) h n m ( t ) i (17)where the expectation values are calculated self-consistently at each timestep. We then use the flow equations (Eqs.15) to transform the number operator back into the original basis, where it will take the form: n i ( l = 0 , t ) = X j A ( i ) j ( t ) n j + X jk B ( i ) jk ( t ) c † j c k (18)At this point, the expectation value of this operator may be computed with respect to the desired initial state. Inthe main text, our initial state is always a product state of the form | ... i . Note that this procedure also gives usa very natural way to understand operator spreading [10–14], as the real-space support of the operator is encoded inthe time-dependent coefficients A ( i ) j ( t ) and B ( i ) jk ( t ), which can be straightforwardly extracted [7]. Accuracy: Dynamical Comparison with Exact Diagonalisation
In order to verify the accuracy of our method, we benchmarked the dynamic results with exact diagonalisation(ED). For this, we again employed the QuSpin package [8, 9]. Sample results for the density dynamics on a singlesite are shown in Fig. 3 for a variety of values of α and β across the phase diagram. The agreement in all cases isexcellent, with flow equations differing only very slightly from the exact results. FIG. 3: The density dynamics on the central site of a chain of length L = 12 when quenched from a CDW initial state andaveraged over 512 disorder realisations, comparing ED (blue) with FE (orange). a) α = 0 . , β = 2 .
0, b) α = 2 . , β = 2 . α = 0 . , β = 0 .
5. d) α = 3 . , β = 0 .
5. In all cases, the results are close, but the FE method slightly overestimates thelocalisation. In the more strongly localised phases for α, β (cid:29)
1, the FE and ED results agree very closely. The insets show thedecay of fluctuations around their long-time mean value, with δn = σ ( h n L/ i − ˜ n ) and ˜ n = h n L/ ( t ) i t →∞ : note the power-lawdecay in the ED data which is not seen in the FE data, due to the mean-field decoupling employed for the dynamics. Despite this excellent agreement, it is interesting to note that the results from the flow equation method do notcapture the decay of fluctuations around their mean values (shown in the insets of Fig. 3). The reason for this is dueto the mean-field decoupling, which does not allow for the slow build-up of correlations that leads to the power-lawdecay of fluctuations (or to the logarithmic growth of entanglement entropy). Similar results are seen in the quantumFisher information (not shown), a proxy for the entanglement entropy, which does not display the expected slowincrease with time due to the nature of the mean-field decoupling used here in computing the dynamics.
FIG. 4: The same quantity as in Fig. 4 of the main text, here for system size L = 36 and averaged over N
00 disorderrealisations. The solid white line represents I ( t ∗ = 100) = 0 .
25 (half the maximum value) for the L = 36 system, and is arough indicator of the position of the transition, while the dashed white line is the same quantity for the L = 64 system shownin Fig. 4 of the main text. There is a clear drift of the boundary towards larger values of α as we increase the system size,however the main features are robust. Phase Diagram: Effect of System Size
To verify our conclusions, we have also computed the phase diagram for a chain of L = 36 sites averaged over N s = 100 disorder realisations, shown in Fig. 4. The phase boundary moves, as expected, but the general conclusionis the same. This demonstrates that the main features of the phase diagram presented in the main text are robust.The flow invariant remains below a maximum value of δI max = 0 .
012 at all points in this figure. This data suggeststhat, all other things being equal, there is a slow growth of the number of resonances as the system size is increased,consistent with the resonance counting arguments in the existing literature. Our results are an indication that evenfor large system sizes, localisation still persists over a large region of the phase diagram. Note however that thereversal of curvature seen in Fig. 4 of the main text for α > L = 36 systemis more localised in this region, with a larger imbalance. This is consistent with the idea that larger systems exhibitmore delocalising resonances, destabilising the localised phase. Random-Sign Disorder
Previous works on long-range couplings in spin chains have considered so-called ‘random sign disorder’, in whichthe couplings are fixed in magnitude but allowed to vary in sign, i.e. J ij = ± J / | i − j | α and V ij = ± V / | i − j | β wherethe signs are chosen randomly. These works have predicted the absence of a localised phase in the regime d ≤ β ≤ d ,whereas we find clear signs of localisation in this regime. While this could be a finite-size effect, or equivalently wemay simply be below the critical disorder threshold for this system size, we have nonetheless simulated this type ofdisorder as well in order to compare with our (zero mean) Gaussian-distributed random couplings. The results areshown in Fig. 5.Remarkably, we find that the case of Gaussian-distributed random couplings is indeed significantly more localisedthan the pure random-sign disorder, both quantiatively and qualitatively. This difference, while striking at firstsight, can be explained simply by the typical magnitude of the coupling terms being large (and, crucialy, non-zero)in the case of random-sign disorder, while the typical value is identically zero for the Gaussian-distributed disorderconsidered in the main text. FIG. 5: Various static properties of the fixed-point Hamiltonian with random-sign disorder, rather than Gaussian-distributeddisorder. All data here is taken for system sizes L = 64 with N s = 128 disorder realisations, and the colour schemes are thesame as in the main text. The left column shows data for long-range hopping, while the rigth column shows data for long-rangeinteractions. a) Fixed-point couplings ∆ ij in the case of power-law hopping and nearest-neighbour interactions β → ∞ , againwith α ∈ [0 ,
5] as in the main text. The black dashed lines are the same as in the main text. b) The same quantity plotted forthe case of power-law interactions (with α → ∞ and β ∈ [0 ,
5] as before). c) The real-space support of the l -bits in the caseof long-range hopping. The black dashed line is the same as in the main text (the α, β → ∞ limit with Gaussian distributeddisorder), while the blue dots show the α, β → ∞ limit of the random-sign disorder. d) The same quantity plotted for long-rangeinteractions, with α → ∞ and β ∈ [0 , α, β → ∞ limit of the random-sign disorder. ∗ Electronic address: [email protected] † Electronic address: [email protected][1] F. Wegner, Journal of Physics A: Mathematical and General , 8221 (2006).[2] S. Kehrein, The flow equation approach to many-particle systems , vol. 217 (Springer, 2007).[3] S. Thomson and M. Schir´o, Physical Review B , 060201 (2018).[4] F. Wegner, Annalen der physik , 77 (1994).[5] C. Monthus, Journal of Physics A: Mathematical and Theoretical , 305002 (2016).[6] S. Savitz and G. Refael, Physical Review B , 115129 (2017).[7] S. J. Thomson and M. Schir´o, Eur. Phy. J. B (2020).[8] P. Weinberg and M. Bukov, SciPost Phys. , 003 (2017).[9] P. Weinberg and M. Bukov, SciPost Phys. , 20 (2019).[10] M. Foss-Feig, Z.-X. Gong, C. W. Clark, and A. V. Gorshkov, Phys. Rev. Lett. , 157201 (2015).[11] A. Nahum, S. Vijay, and J. Haah, Phys. Rev. X , 021014 (2018).[12] V. Khemani, A. Vishwanath, and D. A. Huse, Phys. Rev. X , 031057 (2018).[13] S. Gopalakrishnan, D. A. Huse, V. Khemani, and R. Vasseur, Phys. Rev. B , 220303 (2018).[14] D. J. Luitz and Y. Bar Lev, Phys. Rev. A99