Local integrals of motion and the quasiperiodic many-body localization transition
Hansveer Singh, Brayden Ware, Romain Vasseur, Sarang Gopalakrishnan
LLocal integrals of motion and the quasiperiodic many-body localization transition
Hansveer Singh , Brayden Ware , Romain Vasseur , and Sarang Gopalakrishnan , Department of Physics, University of Massachusetts, Amherst, MA 01003, USA Department of Physics, The Pennsylvania State University, University Park, PA 16820, USA College of Staten Island, Staten Island, NY 10314, USA
We study the many body localization (MBL) transition for interacting fermions subject toquasiperiodic potentials by constructing the local integrals of motion (LIOMs) in the MBL phaseas time-averaged local operators. We study numerically how these time-averaged operators evolveacross the MBL transition. We find that the norm of such time-averaged operators drops discontinu-ously to zero across the transition; as we discuss, this implies that LIOMs abruptly become unstableat some critical localization length of order unity. We analyze the LIOMs using hydrodynamicprojections and isolating the part of the operator that is associated with interactions. Equippedwith this data we perform a finite-size scaling analysis of the quasiperiodic MBL transition. Ourresults suggest that the quasiperiodic MBL transition occurs at considerably stronger quasiperiodicmodulations, and has a larger correlation-length critical exponent, than previous studies had found.
Intuition suggests that isolated many-body systemsinitialized out of equilibrium should “thermalize” undertheir intrinsic unitary dynamics, in the sense of approach-ing a state in which local observables and correlationfunctions exhibit equilibrium behavior [1–4]. Since An-derson’s work [5], it has been understood that thermal-ization is not fully generic: systems subject to strongquenched randomness can instead exhibit a many-bodylocalized (MBL) phase, in which thermalization fails, andthe system instead retains a local memory of its ini-tial conditions to arbitrarily late times [4, 6–11]. Thismemory is due to the existence in the MBL phase of(quasi-)local operators that are exact integrals of mo-tion, called LIOMs or l-bits [11–18]. The existence ofthe MBL phase and of LIOMs has been established un-der minimal assumptions in one-dimensional random spinchains [11]; experimental evidence for the MBL phase ex-ists in many different settings [19–32] (but see Refs. [33–39]). Rare regions—i.e., regions of a sample in whichthe disorder is anomalously weak or strong—play a cen-tral part in our understanding of MBL, determining thenature of the MBL transition [40–49], response on bothsides of the transition [50–61], and even the stability ofthe MBL phase in higher dimensions [54, 62]. However,many experiments on the MBL phase treat systems sub-ject to quasiperiodic (QP) rather than random poten-tials [19, 23, 63–68]. Noninteracting 1d QP systems havea localized phase [69], which appears to be perturbativelystable in the presence of interactions [63]. However, rareregions are absent in QP systems, so it seems that theMBL transition—and the response near it—must differqualitatively from the random transition [64]. The nu-merical evidence on the QP-MBL transition is mixed,with some studies casting doubt on whether a transitionexists at all [70], while others find a breakdown of diffu-sion [65, 71], potentially even in a regime where single-particle states are delocalized [72, 73].Most work on QP-MBL systems has worked in theSchr¨odinger picture, considering the properties of typical individual eigenstates across the transition. The responseof typical eigenstates in the MBL phase to a probe willinvolve both the external QP potential and self-generated configurational randomness, from the random patternof occupation of localized orbitals (which exert randomHartree shifts on one another). Thus from the eigenstateperspective there is no clear distinction between randomand QP MBL systems; since the transition is really aninstability of the MBL phase, one might be led to con-clude that the transition should also be the same. In thepresent work, we instead take the
Heisenberg perspec-tive and focus on the properties of the LIOMs as opera-tors [74]. In the QP-MBL phase, the pattern of LIOMsis quasiperiodic, with LIOMs approximately repeating atregular distances that are rational approximants to theQP pattern; there are no rare regions with anomalous LI-OMs. Thus from this operator perspective the QP andrandom MBL phases differ, and one would also expectthe transition at which LIOMs cease to exist to differ, ifrare regions are indeed important. (Whether this transi-tion coincides with the transition into the thermal phaseis an issue we revisit below.)In the present work we explicitly construct LIOMs inQP systems by time-averaging local operators, as firstproposed in Ref. [15]. We perform the infinite-time aver-age explicitly, via full diagonalization (we also exploretensor-network methods [75]). We analyze these LI-OMs by computing the fraction of the operator normthat comes from n -fermion terms in the expansion O = A ij c † i c j + B ijkl c † i c † j c k c l + . . . , using tensor-network meth-ods to efficiently extract these quantities [75]. These n -fermion weights give us a handle on the specifically many-body effects that occur at the transition: unlike transportand entanglement, it is not contaminated by the single-particle critical point, which lies somewhat near the ap-parent many-body transition. We find that the n -fermionweights and the norm of the LIOMs give us new ways ofanalyzing the transition, pointing to a transition thatoccurs at larger values of the QP potential, with differ- a r X i v : . [ c ond - m a t . d i s - nn ] J a n - - - - - - | l og Q |
4. Inset shows the evolution of N sub with system size, indicating an instability that sets inas the system size is increased, even for h = 3 . h c , ν . The yellow dot marks the minimum of the cost function; by eye, the best collapse is for slightly differentparameters (green dot) [these collapses are shown in panels (e) and (f)]. The transition seen in level statistics is marked by thered dot; our data are clearly inconsistent with these values. For all figures, averaging was performed over 200 phases equallyspaced in the interval [0 , π ). The color scheme for different system sizes is shared for figures (b), (c), (e) and (f). ent critical exponents, than previously expected. Thistransition has notable similarities to the random case: inparticular, the LIOMs slightly on the MBL side of thetransition are tightly localized. Thus, as in the randomcase, it seems that the QP-MBL transition is an instabil-ity of the localized phase, which sets in at some criticalvalue of the localization length. The microscopic originof this instability remains unclear. Model.–
We consider the following model, H = L − (cid:88) i =1 σ xi σ xi +1 + σ yi σ yi +1 + V σ zi σ zi +1 + L (cid:88) i =1 h i σ zi , (1)where σ x,y,z denote the Pauli matrices and h i = h cos(2 π/ϕ ( i − L/
2) + φ ) where ϕ = √ and φ isa phase offset that we tune to translate our window.When V = 0, this is the noninteracting Aubry-Andr´emodel which has localized eigenstates for h > h <
2. At nonzero V , finitesize exact-diagonalization studies [63–65] of the averageeigenstate entanglement and level statistics ratio havefound an MBL transition at h c ≈ ν ∼
1. (However, studies on longer spin-chainsusing the time-dependent variational principle have seen a larger critical point, consistent with ours [68].) In thisletter we set V = 1 / O , which we choose to be σ zL/ . The time average of O is given by¯ O ≡ lim T →∞ T (cid:90) T dt O ( t ) = (cid:88) E (cid:104) E | O | E (cid:105)| E (cid:105)(cid:104) E | , (2)where | E (cid:105) are eigenstates of H . In the MBL phase, weexpect ¯ O to be an approximately local operator with ex-ponential tails, i.e., we expect there is some operator O n with support on n sites such that (cid:107) O − O n (cid:107) ≤ e − n/ξ where ξ is a characteristic localization length. In theergodic phase, the time average instead produces a non-local integral of motion, predominately the projection of O onto conserved charges. We construct ¯ O by full exactdiagonalization, using Eq. (2). We explore finite-time av-erages, performed using exact diagonalization as well asmatrix-product methods, in [75]. Fermion weights. — We analyze the LIOMs by expand-ing them in a basis of n -fermion operators. These arerelated to the Pauli operators by a Jordan-Wigner trans-formation, and evidently form a complete basis:¯ O = (cid:88) α c α w α w α · · · w α L L , (3)where { w i , w j } = 2 δ ij are Majorana fermions. In whatfollows we will focus on two quantities: the Frobeniusnorm of the operator, N ≡
Tr( ¯ O ), and its two-bodyweight f , defined by f = 1 N (cid:88) | α | =2 | c α | , N = (cid:88) α | c α | . (4)We can also define four- and six-body weights accord-ingly. As f is 1 for quadratic fermion operators, 1 − f = f + f + . . . measures the many-body content of the op-erator. (These weights can be efficiently computed us-ing matrix-product operator methods; we outline thesemethods and present results on the weights f n with n > N and f probe complementary as-pects of the time evolution: N addresses how much ofthe initial-state information survives in the time aver-age, while f addresses what fraction of this informationis encoded in “simple” (i.e., few-body) operators. Hydrodynamic modes .— Fig. 1(a) shows the averageand distribution of f . As we might expect, f ap-proaches 1 and we also find that N is of order unityin the MBL phase, since in this phase the LIOMs areapproximately single-site occupation numbers. Less ob-viously, f also approaches 1 deep in the thermal phase(although N , not shown, goes to zero with system size).One can understand this as follows. The operator O = Z i has some overlap with the total conserved charge, I ≡ (cid:80) i Z i , which is conserved (and is a two-fermion opera-tor), and also with the Hamiltonian I = H (which con-tains two- and four-point operators). More generally, ina system of size L , there are 2 L nonlocal conserved oper-ators, i.e., projectors onto eigenstates, while the operatorHilbert space is 4 L -dimensional. Since the operator O atlate times under chaotic dynamics is essentially random,its projection onto the conserved eigenstates would be ex-ponentially small in L if it were not for local conservationlaws. Neglecting these exponentially small components,one can write ¯ O in the thermal phase as its projectiononto hydrodynamic modes using the (super)projector P = (cid:88) l,k =1 , | I k (cid:105)(cid:105) C − kl (cid:104)(cid:104) I l | , (5)where I = (cid:80) i Z i and I = H the conserved charges,acting on a Hilbert space H , are now viewed as stateson the doubled Hilbert space H ⊗ H , and C kl = (cid:104)(cid:104) I k | I l (cid:105)(cid:105) the susceptibility matrix with (cid:104)(cid:104) A | B (cid:105)(cid:105) ≡ − L tr( A † B ).Since H and Q are both composed of two- and four-body operators, f remains of order unity throughoutthe thermal phase.Since the hydrodynamic modes exist on both sides ofthe transition and the projection of an operator onto these modes is a property of the t = 0 operator thatis insensitive to critical properties, we subtract the pro-jection onto this hydrodynamic subspace and define the“subtracted operator”¯ O sub ≡ ¯ O − P ( O ) (cid:107) O − P ( O ) (cid:107) . (6)The denominator in Eq. (6) corrects for the fact that thehydrodynamic projection of O smoothly increases withincreasing disorder, since the Hamiltonian is dominatedby single-site potential terms. (Empirically, we find thatnot fixing the normalization of ¯ O sub leads to spuriousfinite-size drifts in small-system numerics.) The norm N sub is defined as for the full operator. We also introducethe subtracted two body component, ˜ f , defined as˜ f , sub ≡ (cid:80) i 1. InFig. 1(d) we have plotted the figure of merit [specifically,the log quality factor | log Q | extracted from the Pythonpackage pyfssa [79]] for attempted data collapses withvarious possible combinations of h c and ν : our numericaldata for the transition in the LIOMs are evidently incon-sistent with previous predictions for the critical point. Ifanything, our scaling collapses show weak drift to largervalues of h c , suggesting that the transition might occureven deeper in the apparent MBL phase than our esti-mates above.A clue as to why our results look so different from pre-vious studies can be gleaned from the inset to Fig. 1(b).For h = 3 , . N sub remains largefor all the system sizes we study. However, it decreaseswith system size in a way that is accelerating at larger L . This pattern is not consistent with the expectedfinite-size effects in the MBL phase, which should scaleas e − L/ξ , where ξ is the correlation length, and shouldtherefore flatten out at larger sizes. Rather, these resultssupport a scenario in which the system seems localizedon short scales, but then (beyond some critical scale ξ )realizes it is unstable. The scale ξ diverges as h is in-creased toward h c , and can be regarded as a correlationlength. Discussion .— In this work we studied the propertiesof LIOMs across the quasiperiodic many-body localiza-tion transition, constructing them as time-averaged localoperators. In the thermal phase, these time-averaged op-erators are just projectors onto global conserved quanti-ties like the total energy and charge; once these hydro-dynamic parts are subtracted out, the remainder of theoperator vanishes rapidly. In the MBL phase, instead, a time-averaged local operator retains a finite norm, sinceit has non-hydrodynamic projections onto the LIOMs.There is a transition at which these LIOMs cease to existand the norm of the subtracted time-averaged operatorvanishes. This apparent transition has a critical pointand critical exponents that are inconsistent with the ap-parent transition in observables such as eigenstate entan-glement and nearest-neighbor level statistics. Notably,the apparent correlation-length exponent ν ∈ (2 , . ν ≥ ν = 1, saturating the Luck bound).Indeed, we should emphasize that the results we havepresented do not constitute strong evidence for the exis-tence of an MBL phase at all, and are in principle con-sistent with a transition that occurs at h c = ∞ ; how-ever, the MBL phase is perturbatively stable for suffi-ciently large h , and no nonperturbative instabilities haveyet been identified, so we take the point of view thatthere is a transition in the window where we see one.A counterintuitive implication of our results is that if ν ≥ 2, the QP-MBL critical point is perturbatively sta-ble against weak randomness by the Harris criterion [80].(We note that a similar result was found using a real-space RG scheme in Ref. [66].) This could suggest, ei-ther that there is a critical value of randomness requiredto change the universality class of the QP-MBL transi-tion, or that both the QP-MBL critical point and somepart of the QP-MBL phase undergo a nonperturbativeinstability for infinitesimal randomness.How can we reconcile our observations with the resultson level statistics and entanglement? One possibility isthat there are two separate transitions with distinct crit-ical properties, one at which the level statistics changesits character and another at which LIOMs cease to exist.This could happen, for example, if there were an inter-mediate phase with a many-body mobility edge [81, 82].However, it is unclear whether such many-body mobilityedges can exist [82], and even if they do, the transitionin entanglement should occur once the entire spectrumis localized. Thus it is not clear how this scenario couldapply to our case. A second possibility is that the LI-OMs we study here have weaker finite-size effects thanthe level statistics, because they are less affected by state-to-state fluctuations that are large in finite systems [64].(All known finite-size effects favor the MBL phase, so ahigher h c value is more plausible, assuming there is asingle transition.)Our results shed light on the nature of this transitionat which LIOMs cease to exist, which we tentatively iden-tify with the MBL transition. In particular, we find thatLIOMs even slightly on the MBL side are mostly fermionbilinears with large norm; thus they overlap strongly withmicroscopic spins. The QP transition, like the randomone, appears to be an instability of the MBL phase thatsets in at some critical localization length as one increasesthe system size. In random systems, such an instabil-ity is thought to be seeded by rare regions that are lo-cally thermal. Although rare regions do not exist, strictlyspeaking, in QP systems, one might still expect the in-stability to occur first in some parts of the sample. Onemight expect LIOMs to be unusually delocalized in sam-ples that contain these parts. However, we do not seemuch evidence of enhanced heterogeneity at the transi-tion (Fig. 1). It therefore seems that the instability weare seeing is due to the proliferation of many-body res-onances in typical regions of the sample. The origin ofthese resonances remains to be identified.Our work, like all ED studies, is inherently limitedto small system sizes. An important question for futurework is whether one can construct LIOMs for much largersystems. We attempted to do so by time evolving lo-cal operators via time-evolving block decimation (TEBD)applied to matrix-product operators [75], and averagingover finite time windows. Unfortunately, to get a goodapproximation to the LIOMs away from the deeply lo-calized limit, one must average over such long time win-dows (comparable to the Heisenberg time) that TEBD isimpractical [75], because the bond dimension needed todescribe the operator grows intractably large. Whetherother forms of explicit time-averaging, e.g., based onKrylov-space methods [59, 83], can provide access tolarger systems and sufficiently long times is an interestingquestion for future work.The authors thank U. Agrawal, P. Dumitrescu, D.Huse, V. Khemani and V. Oganesyan for helpful dis-cussions. S.G. acknowledges support from NSF DMR-1653271. R.V. acknowledges support from the US De-partment of Energy, Office of Science, Basic Energy Sci-ences, under Early Career Award No. DE-SC0019168and the Alfred P. Sloan Foundation through a Sloan Re-search Fellowship. [1] J. M. Deutsch, Phys. Rev. A , 2046 (1991).[2] M. Srednicki, Phys. Rev. E , 888 (1994).[3] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii,Phys. Rev. Lett. , 050405 (2007).[4] R. Nandkishore and D. A. Huse, Annual Review of Con-densed Matter Physics , 15 (2015).[5] P. W. Anderson, Phys. Rev. , 1492 (1958).[6] D. M. Basko, I. L. Aleiner, and B. L. Altshuler, Annalsof Physics , 1126 (2006).[7] I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Phys.Rev. Lett. , 206603 (2005).[8] D. M. Basko, I. L. Aleiner, and B. L. Altshuler, cond-mat/0602510.[9] V. Oganesyan and D. A. Huse, Phys. Rev. B , 155111(2007).[10] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev.Mod. Phys. , 021001 (2019).[11] J. Z. Imbrie, Journal of Statistical Physics , 998 (2016).[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] I. H. Kim, A. Chandran, and D. A. Abanin, ArXiv e-prints (2014), arXiv:1412.3073 [cond-mat.dis-nn].[15] A. Chandran, I. H. Kim, G. Vidal, and D. A. Abanin,Phys. Rev. B , 085425 (2015).[16] V. Ros, M. M¨uller, and A. Scardicchio, Nuclear PhysicsB , 420 (2015).[17] J. Z. Imbrie, V. Ros, and A. Scardicchio, Annalen derPhysik , 1600278 (2017).[18] S. Gopalakrishnan and S. Parameswaran, Physics Re-ports (2020).[19] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. L¨uschen,M. H. Fischer, R. Vosk, E. Altman, U. Schneider, andI. Bloch, Science , 842 (2015).[20] S. S. Kondov, W. R. McGehee, W. Xu, and B. DeMarco,Phys. Rev. Lett. , 083002 (2015).[21] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess,P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, NaturePhysics , 907 (2016).[22] J.-y. Choi, S. Hild, J. Zeiher, P. Schauß, A. Rubio-Abadal, T. Yefsah, V. Khemani, D. A. Huse, I. Bloch,and C. Gross, Science , 1547 (2016).[23] P. Bordia, H. L¨uschen, S. Scherg, S. Gopalakrishnan,M. Knap, U. Schneider, and I. Bloch, Physical ReviewX , 041047 (2017).[24] P. Bordia, H. L¨uschen, U. Schneider, M. Knap, andI. Bloch, Nature Physics , 460 (2017).[25] P. Roushan, C. Neill, J. Tangpanitanon, V. Bastidas,A. Megrant, R. Barends, Y. Chen, Z. Chen, B. Chiaro,A. Dunsworth, et al. , Science , 1175 (2017).[26] K. Xu, J.-J. Chen, Y. Zeng, Y.-R. Zhang, C. Song,W. Liu, Q. Guo, P. Zhang, D. Xu, H. Deng, K. Huang,H. Wang, X. Zhu, D. Zheng, and H. Fan, Phys. Rev.Lett. , 050507 (2018).[27] A. Lukin, M. Rispoli, R. Schittko, M. E. Tai, A. M. Kauf-man, S. Choi, V. Khemani, J. L´eonard, and M. Greiner,Science , 256 (2019).[28] B. Chiaro, C. Neill, A. Bohrdt, M. Filippone, F. Arute,K. Arya, R. Babbush, D. Bacon, J. Bardin, R. Barends, et al. , Preprint at https://arxiv. org/abs/1910.06024(2019).[29] M. Rispoli, A. Lukin, R. Schittko, S. Kim, M. E. Tai,J. L´eonard, and M. Greiner, Nature , 385 (2019).[30] J. L´eonard, M. Rispoli, A. Lukin, R. Schittko, S. Kim,J. Kwan, D. Sels, E. Demler, and M. Greiner, arXivpreprint arXiv:2012.15270 (2020).[31] I. Tamir, T. Levinson, F. Gorniaczyk, A. Doron, J. Lieb,and D. Shahar, Phys. Rev. B , 035135 (2019).[32] T. Nguyen, N. Andrejevic, H. C. Po, Y. Tsurimaki, N. C.Drucker, A. Alatas, E. E. Alp, B. M. Leu, A. Cunsolo,Y. Q. Cai, et al. , arXiv preprint arXiv:2008.02257 (2020).[33] J. ˇSuntajs, J. Bonˇca, T. Prosen, and L. Vidmar, arXive-prints , arXiv:1905.06345 (2019), arXiv:1905.06345[cond-mat.str-el].[34] J. ˇSuntajs, J. Bonˇca, T. Prosen, and L. Vidmar, arXivpreprint arXiv:2004.01719 (2020).[35] D. Sels and A. Polkovnikov, arXiv preprintarXiv:2009.04501 (2020).[36] T. LeBlond, D. Sels, A. Polkovnikov, and M. Rigol,arXiv preprint arXiv:2012.07849 (2020). [37] D. Abanin, J. Bardarson, G. De Tomasi, S. Gopalakr-ishnan, V. Khemani, S. Parameswaran, F. Pollmann,A. Potter, M. Serbyn, and R. Vasseur, arXiv preprintarXiv:1911.04501 (2019).[38] R. K. Panda, A. Scardicchio, M. Schulz, S. R. Taylor,and M. ˇZnidariˇc, EPL (Europhysics Letters) , 67003(2020).[39] P. J. Crowley and A. Chandran, arXiv preprintarXiv:2012.14393 (2020).[40] R. Vosk, D. A. Huse, and E. Altman, Phys. Rev. X ,031032 (2015).[41] A. C. Potter, R. Vasseur, and S. A. Parameswaran, Phys.Rev. X , 031033 (2015).[42] P. T. Dumitrescu, R. Vasseur, and A. C. Potter, Phys.Rev. Lett. , 110604 (2017).[43] M. Serbyn, Z. Papi´c, and D. A. Abanin, Phys. Rev. X , 041047 (2015).[44] T. Thiery, F. m. c. Huveneers, M. M¨uller, andW. De Roeck, Phys. Rev. Lett. , 140601 (2018).[45] A. Goremykina, R. Vasseur, and M. Serbyn, Phys. Rev.Lett. , 040601 (2019).[46] P. T. Dumitrescu, A. Goremykina, S. A. Parameswaran,M. Serbyn, and R. Vasseur, Phys. Rev. B , 094205(2019).[47] A. Morningstar and D. A. Huse, Phys. Rev. B , 224205(2019).[48] A. Morningstar, D. A. Huse, and J. Z. Imbrie, Phys.Rev. B , 125134 (2020).[49] V. Khemani, S. P. Lim, D. N. Sheng, and D. A. Huse,Phys. Rev. X , 021013 (2017).[50] K. Agarwal, S. Gopalakrishnan, M. Knap, M. M¨uller,and E. Demler, Phys. Rev. Lett. , 160401 (2015).[51] Y. Bar Lev, G. Cohen, and D. R. Reichman, Phys. Rev.Lett. , 100601 (2015).[52] S. Gopalakrishnan, M. M¨uller, V. Khemani, M. Knap,E. Demler, and D. A. Huse, Phys. Rev. B , 104202(2015).[53] S. Gopalakrishnan, K. Agarwal, E. A. Demler, D. A.Huse, and M. Knap, Phys. Rev. B , 134206 (2016).[54] W. De Roeck and F. m. c. Huveneers, Phys. Rev. B ,155129 (2017).[55] D. J. Luitz, F. Huveneers, and W. De Roeck, Phys. Rev.Lett. , 150602 (2017).[56] L. Zhang, B. Zhao, T. Devakul, and D. A. Huse, Phys.Rev. B , 224201 (2016).[57] P. J. D. Crowley and A. Chandran, Phys. Rev. Research , 033262 (2020).[58] K. Agarwal, E. Altman, E. Demler, S. Gopalakrishnan,D. A. Huse, and M. Knap, Annalen der Physik ,1600326 (2017).[59] D. J. Luitz, N. Laflorencie, and F. Alet, Phys. Rev. B , 060201 (2016).[60] M. ˇZnidariˇc, A. Scardicchio, and V. K. Varma, Physicalreview letters , 040601 (2016).[61] T. L. M. Lezama, S. Bera, and J. H. Bardarson, Phys.Rev. B , 161106 (2019).[62] S. Gopalakrishnan and D. A. Huse, Phys. Rev. B ,134305 (2019).[63] S. Iyer, V. Oganesyan, G. Refael, and D. A. Huse, Phys.Rev. B , 134202 (2013).[64] V. Khemani, D. N. Sheng, and D. A. Huse, Phys. Rev.Lett. , 075702 (2017).[65] F. Setiawan, D.-L. Deng, and J. H. Pixley, Phys. Rev.B , 104205 (2017).[66] S.-X. Zhang and H. Yao, Physical review letters ,206601 (2018).[67] S. A. Weidinger, S. Gopalakrishnan, and M. Knap, Phys.Rev. B , 224205 (2018).[68] E. V. H. Doggen and A. D. Mirlin, Phys. Rev. B ,104203 (2019).[69] S. Aubry and G. Andr´e, Ann. Israel Phys. Soc , 18(1980).[70] M. ˇZnidariˇc and M. Ljubotina, Proceedings of the Na-tional Academy of Sciences , 4595 (2018).[71] Y. B. Lev, D. M. Kennes, C. Kl¨ockner, D. R. Reichman,and C. Karrasch, EPL (Europhysics Letters) , 37003(2017).[72] X. Li, S. Ganeshan, J. H. Pixley, and S. Das Sarma,Phys. Rev. Lett. , 186601 (2015).[73] S. Xu, X. Li, Y.-T. Hsu, B. Swingle, and S. Das Sarma,Phys. Rev. Research , 032039 (2019).[74] U. Agrawal, S. Gopalakrishnan, and R. Vasseur, Naturecommunications , 1 (2020).[75] See Online Supplementary Information.[76] P. Ponte, C. Laumann, D. A. Huse, and A. Chandran,Phil. Trans. R. Soc. A , 20160428 (2017).[77] L. Herviou, S. Bera, and J. H. Bardarson, Phys. Rev. B , 134205 (2019).[78] N. Laflorencie, G. Lemari´e, and N. Mac´e, Phys. Rev.Research , 042033 (2020).[79] A. Sorge, “pyfssa 0.7.6,” (2015).[80] A. B. Harris, Journal of Physics C: Solid State Physics , 1671 (1974).[81] D. J. Luitz, N. Laflorencie, and F. Alet, Phys. Rev. B , 081103 (2015).[82] W. de Roeck, F. Huveneers, M. M¨uller, and M. Schiulaz,ArXiv e-prints (2015), arXiv:1506.01505 [cond-mat.dis-nn].[83] K. H´emery, F. Pollmann, and D. J. Luitz, Phys. Rev. B , 104303 (2019). upplemental Material: Local integrals of motion and the quasiperiodic many-bodylocalization transition Hansveer Singh , Brayden Ware , Romain Vasseur , and Sarang Gopalakrishnan , Department of Physics, University of Massachusetts, Amherst, MA 01003, USA Department of Physics, The Pennsylvania State University, University Park, PA 16820, USA College of Staten Island, Staten Island, NY 10314, USA (Dated: January 10, 2021) I. COMPUTATION OF THE k -BODY FERMION WEIGHT Here we provide a formula to compute the k -body fermion weight. To start we write an operator, O , in the basisof Majorana fermions O = X α c α w α w α · · · w α L L , (S1)where w i are hermitian operators which satisfty { w i , w j } = 2 δ ij . The k -body fermion weight is given by, f k = 1 N X | α | = k | c α | , N = X α | c α | . (S2)It is straightforward to show that, X | α | = k | c α | = 14 L X | α | = k | tr( O † w α w α · · · w α L L ) | , (S3)where | α | ≡ P i α i , and it is also straightforward to see that X α | c α | = 12 L tr( O † O ) . (S4)Plugging in the above relations into Eq. (S2), one finds f k = 12 L tr( O † O ) X | α | = k | tr( O † w α w α · · · w α L L ) | . (S5) II. CONSTRUCTING L-BITS WITH TEBD As discussed in the main text the l-bits are obtained through time averaging, i.e.¯ O = lim T →∞ T Z T dt O ( t ) . (S6)¯ O is expected to be approximately local (starting from a local operator) deep in the MBL phase, so ¯ O can beefficiently represented as an MPO [S1]. Our main text suggests that deep in the thermal phase ¯ O is dominated bythe hydrodynamic modes of our system and so it should be well described by an MPO since the only (dominant)hydrodynamic modes are H and M . Near the transition it is less clear how accurate an MPO description is of ¯ O .Denoting ξ as the characteristic length over which the l-bits have support (this is only well defined in the MBLphase) then intuition from the random case suggests that ξ should remain finite as one approaches the delocalizationtransition from the MBL phase (see Section VII). Motivated by this idea we attempted to construct the l-bits usingTEBD. A. Method Description Proceeding with the discussion of ¯ O ’s construction, we define the finite time average ¯ O ( T ) ≡ T R T dtO ( t ). Dis-cretizing time into intervals of ∆ t and denoting M ≡ T / ∆ t , then¯ O ( T ) ≈ M M X n =1 O ( n ∆ t ) + O ((∆ t ) ) . (S7)Using this approximation we have the following relation,¯ O ( T ) = 1 M (( M − 1) ¯ O ( T − ∆ t ) + O ( T )) . (S8)Keeping track of the finite time average, we can use this formula to construct ¯ O ( T ) by adding the two MPO repre-sentations of ¯ O ( T − ∆ t ) and O ( T ). If ¯ O ( T − ∆ t ) is described by a bond dimension ψ MPO and O ( T ) is describedby a bond dimension χ MPO then ¯ O ( T ) is described by a bond dimension χ + ψ MPO. Since the enlarged bonddimension, χ + ψ , may not be the optimal bond dimension for ¯ O ( T ) or it may be too large to handle numerically, onegenerally needs a compression scheme to lower the bond dimension.The standard method of compressing an MPO is done by applying an SVD to each tensor in the MPO and discardingthe appropriate number of singular values until either we are below some error threshold, (cid:15) , or above some maximalbond dimension, D . This method will give the optimal tensor of bond dimension D to use for a given site but thismay not be the best global representation of ¯ O at bond dimension D . An improvement to this method is to usean iterative compression scheme [S2] which takes an initial guess for the globally best MPO and then minimizes theFrobenious norm between the truncated average, ¯ O ( T ) approx , and the exact average, i.e. we iteratively minimize || ¯ O ( T ) − ¯ O ( T ) approx || where || A || = tr( A † A ).Let us summarize the above steps: we calculate O ( T ), then using ¯ O ( T − ∆ t ) we construct ¯ O ( T ) using Eq. (S8)and using the SVD compressed approximation of ¯ O ( T ) as our initial guess, we employ iterative compression to bring¯ O ( T ) down to a lower bond dimension.To perform the standard TEBD step we pick ∆ t = 0 . 05 and use a fourth order trotterization scheme of the timeevolution operator given by [S3, S4], e − iH ∆ t ≈ U ( τ ) U ( τ ) U ( τ ) U ( τ ) U ( τ ) + O ((∆ t ) ) , (S9)where U ( t ) = e − iH odd t/ e − iH even t e − iH odd t/ , (S10) τ = ∆ t − / , τ = ∆ t − τ and H even / odd are the terms in the Hamiltonian acting on even/odd bonds in the chain. Aswe perform time evolution we keep track of the operator norm square of ¯ O ( T ) and the two body weight of ¯ O ( T ). Whilethe former quantity is easily computable via standard matrix-product techniques, the fact that f can be efficientlycomputed is less trivial. Not only do we show how this can be done in Section IV for f but also for all f k where k iseven. B. Results Unfortunately we found two big computation walls to this method: convergence time of observables and bonddimension. In Fig. (S1) we demonstrate that already at L = 8 we see a long convergence time is required for f andwhen we are closer to the transition this time scale increases. Furthermore, we see that when we are deeper into theMBL phase we have fairly good convergence with bond dimension-seeing convergence occuring close to χ = 32. Onthe other hand when we are near the transition even for χ = 64 (25% of the Hilbert space dimension) we are unableto converge at long times which suggests that we cannot get near the transition if we were to increase system size. Inaddition we remark that the inaccuracy that develops in the description of ¯ O ( T ) is due to the truncation error from O ( T ) which is a well known issue in performing long time simulations with TEBD. We also note that it may be thecase that ¯ O still has an approximate MPO representation but it is not practical to construct it using TEBD. Thuswe find that using matrix-product methods is not a practical tool to study the MBL transition but could be used toexamine quantities deep within the MBL phase as previous studies have done [S1, S5]. We will however use MPOtechniques to compute the k -body weight even for system sizes accessible to exact diagonalization (see below). FIG. S1. The ( left ) and ( middle ) plots show the time evolution for f with L = 8, zero global phase shift, and h = 4 . h = 5 respectively. It is evident from these two plots that convergence with bond dimension is exceedingly better when we aredeeper inside the MBL phase rather than near the transition. The color scheme for both the ( left ) and ( middle ) are the same.( right ) The four body weight averaged over 200 phases uniformly spaced in the interval [0 , π ) and the error bars signify thestandard deviation. Clear strong finite size effects can be seen in f and based on the maximum of f the critical point seemsto be closer to h ≈ III. ED COMPUTATION OF f k While MPOs are not a viable route to study the full phase diagram, the MPO construction of f k turned out tobe a useful way to compute f k using ED. The standard brute force method would require using Eq. (S5) which isnot very efficient even though the w i are sparse matrices. An improvement to this would be to utilize the U (1)symmetry present in the problem and compute the overlap with normal ordered fermion operators c † i c j , c † i c † j c k c l .,etc. Another way we found that was not extremely expensive was to utilize Eq. (S16) (in the next section) and create O + by arranging all the blocks in a 2 L × L matrix and build H k by contracting the tensors. The cost of matrixmultiplication is O (4 L (2 k + 1) ) but in practice we found that diagonalization was still the dominant computation.In Fig. (S1) we present the average of f and its behavior as a function of quasiperiodic modulation. Contrary to˜ f , sub , this suffers from much larger finite size effects, which accounts for the smaller critical point value. Further finitesize effects do not weaken if we consider its subtracted counterpart (not shown), f , sub , which is defined in Section V. IV. EXTENDING THE k -BODY WEIGHT CALCULATION FOR MATRIX PRODUCT OPERATORS Here we show that the k -body weight can be efficiently computed for matrix product operators. We first write anoperator, O , in the pauli basis, O = X β b β σ β σ β · · · σ β L L . (S11)Using the Jordan-Wigner transformation, w j − = Y k A BC D E FG H I XY Z ZX Y Z XY Y XZ XY Z Z YXZ Y X FIG. S2. Diagrams which construct the weighted finite automata for the ( left ) two body hamiltonian and ( right ) the four bodyhamiltonian. Using this diagrammatic notation and Eq. (S19), one finds that only two diagrams correspond to the two bodyoperators: ←→ X j k − Y i = j +1 Z i Y k , Y j k − Y i = j +1 Z i Y k , X j k − Y i = j +1 Z i X k , Y j k − Y i = j +1 Z i X k ←→ Z j . This suggests the following interpretation: two vertical lines connected by a horizontal line contributes two fermionoperators and an ”X” contributes two fermion operators.One can now write down the diagrams for the four bodyoperators, ←→ two non-overlapping two body strings ←→ a single site Z j with a non-overlapping two body string ←→ two non-overlapping single site Z j ←→ a two body string broken by an insertion of a Z j . The last diagram tells us that a circle contributes two fermion operators while a vertical line with a horizontal lineattached contributes one fermion operator. Accounting for the possible spatial orderings of the diagrams and the factthat a vertical line can either be X j or Y j , one arrives at the diagram that constructs the weighted finite automatonshown in Fig. (S2).The diagrams can be used to construct the matrices as shown in Ref. [S6]. Denoting W σ j σ j [ j ]2 / the matrices in theMPO representation of H / , then W σ σ [1]2 = (cid:0) X Y Z (cid:1) , (S20) W σ j σ j [ j ]2 = X Y Z Z X Y Z Y X 00 0 0 0 , < j < L, (S21) W σ L σ L [ L ]2 = ZX + YX + Y , (S22)and W σ σ [1]4 = (cid:0) X Y Z (cid:1) , (S23) W σ j σ j [ j ]4 = X Y Z Z X Y Z Y X X Y Z 00 0 0 0 X Y Z 00 0 0 0 0 Z X Y Z Y X 00 0 0 0 0 0 0 0 , < j < L, (S24) W σ L σ L [ L ]4 = ZZX + YX + Y . (S25)One can see in Fig. (S2) that the four body diagram is essentially two two body diagrams (the first two bodydiagram is with states A − D and the second is formed with states E − I ) joined together by extra transitions betweenstates. This extends to arbitrary k body weight for k even and the matrices for the k body weight with k even havebond dimension 2 k + 1 and the initial and final MPO tensors are the same but with 0s padding to the right of the firsttensor and the top of the last tensor. The bulk MPO tensor is given by repeating the first four rows of eq. (S24) for2 k − × × V. COMPUTATION OF THE k -BODY WEIGHT FOR THE SUBTRACTED L-BIT Here we calculate the k -body weight of the subtracted l-bit. Given a Hamiltonian with conserved quantities, I k ,which are not necessarily linearly independent, then the (super)projector onto the subspace spanned by the I k is givenby P = X lk | I k ii C − kl hh I l | , (S26)where the I k , acting on a Hilbert space H , are now viewed as states on the Hilbert space H ⊗ H and C kl = hh I k | I l ii with hh A | B ii ≡ − L tr( A † B ). Using Eq. (S5), the k -body weight for the subtracted l-bit, denoted as f k, sub is given by, f k, sub = 12 L tr( ¯ Z ) X | α | = k | tr( ¯ Z sub w α w α · · · w α L L ) | , ¯ Z sub = ¯ Z k − P ( ¯ Z k ) . (S27)In this section we heavily abuse notation and will usually replace ¯ Z k with ¯ Z and Z k with Z . We first calculate thesubtracted norm, i.e. the denominator of Eq. (S27).tr( ¯ Z ) = tr( ¯ Z ) + tr( P ( ¯ Z ) ) − Z P ( ¯ Z ))= tr( ¯ Z ) + tr( P ( ¯ Z ) ) − P ( ¯ Z ) + Q ( ¯ Z )) P ( ¯ Z ))= tr( ¯ Z ) − y T C − y. Going from the second to the third line we wrote ¯ Z = P ( ¯ Z ) + Q ( ¯ Z ) where Q is the (super)projector onto theorthogonal complement of the hydrodynamic subspace and the fact that hh I k | ¯ Z ii = hh I k | Z ii . In the last line weintroduced the vector y k ≡ hh I k | Z ii . Using the fact that H and the magnetization, M ≡ P i Z i , are the only(dominant) hydrodynamic modes, we find C HM = L X k =1 h k , C HH = 2( L − 1) + V ( L − 1) + L X k =1 h k , C MM = L, (S28) y H = h k , y M = 1 . (S29)Before calculating the numerator we remark that since P ( ¯ Z ) only contains operators with only two or four bodyfermion operators (in our model) then the numerator of with ¯ Z sub is the same as ¯ Z . Thus we only need to calculatethe numerator for k = 2 and k = 4. A. Computation of the subtracted 2-body weight Straightforward computation shows that, X | α | =2 | tr( ¯ Z sub w α · · · w α L L ) | = X i 1) + L X k =1 h k ) , F MM = 4 L L. (S36)We now calculate the third term in Eq. (S30). Plugging in the expression for P ( ¯ Z ), we find X i The calculation is carried out in a similar fashion as the 2-body weight. X | α | =4 | tr( ¯ Z sub w α · · · w α L L ) | = X i The l-bits we used were constructed by performing time evolution on a local operator followed by a time average.These l-bits are different than the usual τ α i i which are orthonormal and only exist in the MBL phase. First, theyare not necessarily orthonormal but as was shown in Ref. [S5] the overlaps between l-bits exponentially decay as onegoes deep into the MBL phase for the random case. Second, our l-bits are conserved operators regardless of wherewe sit in the phase diagram. Here we numerically demonstrate that the l-bits are indeed linearly independent in thequasiperiodic case by computing the inverse participation ratio (IPR) and participation ratio (PR) of the singularvalues of the matrix ‘ jk ≡ hh ¯ Z j | ¯ Z k ii , i.e. IPR = P i s i and PR = 1 / IPR. In Fig. S3 we show the average IPR andaverage PR as well as the IPR and PR from the singular values of the averaged matrix h ‘ jk i φ . Recall the number ofnon-zero singular values directly translates to the rank of ‘ jk (or its averaged counterpart). This means deep in thethermal phase we expect the time averaged operators to be linearly dependent so the rank should be close to 1 and sothe IPR should be close to 1. In the MBL phase we expect the time averaged operators to be linearly independent andso all the singular values should be of the same order. Since we normalize the singular values to satisfy P Li =1 s i = 1then each s i should approximately scale as 1 /L and so the IPR should vanish in the MBL phase. As one can see inFig. S3 both of these limiting behaviors are exhibited in the IPR and in the PR (since it is just the reciprocal of theIPR). In addition, a crossing seems to occur for h ∈ (2 , 3) which is considerably lower than the estimate in the maintext but this is due to stronger finite size effects in the IPR than to our observables in the main text.0 - - - - - - - - - - - - - - - - - - FIG. S4. Operator locator density averaged over 176 phase realizations uniformly distributed between [0 , π ). x is the latticesite shifted by half the system size so that all system sizes have the same origin. One can see that the l-bit remains fairlylocalized even as we are approaching the vicnity of the critical point given by either the subtracted norm or subtracted twobody component. FIG. S5. [( left ) and ( right ) Average half entanglement entropy divided by the page value, S T , and the average level statisticsratio, r for the quasiperiodic case ( left ) and the random case ( right ). One can see that a crossing occurs around h ≈ h ≈ M = 0 sector removing 10% of the eigenstates from the extremes ofthe spectrum. - - - FIG. S6. ( left ) Distribution of f of ¯ Z L/ over 899 disorder realizations picked uniformly from [ − h, h ]. Solid lines indicate theaverage of f . One can see stronger finite size effects and broader distributions compared to the quasiperiodic case consistentwith previous studies. ( middle ) Here is the average subtracted norm squared of ¯ Z L/ which seems to exhibit a crossing near h ≈ . 5, but the number of realizations averaged over is too low to make any definitive prediction for h c . The inset shows theunsubtracted norm squared of Z L/ which does not exhibit a crossing. Overall there is qualitatively similar behavior betweenthe random case and the quasiperiodic case but quantitatively the subtracted norm squared is lower than the quasiperiodiccase. ( right ) Here we perform a collapse of the subtracted norm squared data for ν = 2. Although the low sample size doesnot make this estimate precise, we did find that ν = 1 does not gives a worse collapse than ν = 2. VII. SPATIAL STRUCTURE OF THE L-BITS Here we present a characterization of the spatial structure of our l-bits. To do so, we decompose our l-bit in thePauli basis given by Eq. (S11). One can then compute the left-right weight distribution given by, ρ ( z, y ) = 1 N X | β | = | z − y | | b β | , N = X β | b β | , (S45)which gives us the fraction of weight the l-bit has in Pauli strings which start at site z and end at site y . Finally, onecan define another quantity called the operator locator density defined as, L ( x ) = X z ≤ x ≤ y ρ ( z, y ) y − z + 1 , (S46)which gives a weighted average of Pauli strings which have support on site x . This quantity gives a picture of wherethe support of the l-bit is. As can be seen in Fig. (S4), at large h the l-bit only has support near the center of thechain since we constructed our l-bit from Z L/ , while at small h the l-bit has support across the whole chain since thel-bit becomes a random operator in the ergodic phase. As we go towards the transition 4 < h < . VIII. COMPARISON TO THE RANDOM CASE In this section we show how the random case deviates from the quasiperiodic case. We first show in Fig. (S5) theaverage half entanglement entropy, S , divided by the Page value for a random pure state, S T = 1 / L log 2 − r ≡ min(∆ n , ∆ n +1 ) / max(∆ n , ∆ n +1 ) where ∆ n = E n − E n +1 is the spacing betweenenergy eigenvalues for the M = 0 sector. Examining Fig. (S5), we observe a crossing at h ≈ h ≈ f and N sub . Similar to Refs. [S7, S8] we see broader distributionsin f compared to the quasiperiodic case as well as stronger finite size effects. While we only averaged over a smallsample size of 899 disorder realizations, one can see that a crossing seems to occur at higher disorder strengths thanpreviously seen before and with a critical exponent of ν ≈ [S1] D. Pekker and B. K. Clark, Phys. Rev. B , 035116 (2017).[S2] U. Schollw¨ock, Annals of Physics , 96 (2011).[S3] M. Suzuki, Journal of Mathematical Physics , 400 (1991), https://doi.org/10.1063/1.529425.[S4] M. Suzuki, Progress of Theoretical Physics , 1454 (1976).[S5] A. Chandran, I. H. Kim, G. Vidal, and D. A. Abanin, Phys. Rev. B , 085425 (2015).[S6] G. M. Crosswhite and D. Bacon, Phys. Rev. A , 012356 (2008).[S7] V. Khemani, D. N. Sheng, and D. A. Huse, Phys. Rev. Lett. , 075702 (2017).[S8] F. Setiawan, D.-L. Deng, and J. H. Pixley, Phys. Rev. B96