Many-body localization transition in large quantum spin chains: The mobility edge
MMany-body localization transition in large quantum spin chains: The mobility edge
Titas Chanda, ∗ Piotr Sierant, † and Jakub Zakrzewski
1, 2, ‡ Instytut Fizyki Teoretycznej, Uniwersytet Jagiello´nski, (cid:32)Lojasiewicza 11, 30-348 Krak´ow, Poland Mark Kac Complex Systems Research Center, Uniwersytet Jagiello´nski, Krak´ow, Poland. (Dated: June 5, 2020)Thermalization of random-field Heisenberg spin chain is probed by time evolution of densitycorrelation functions. Studying the impacts of average energies of initial product states on dynamicsof the system, we provide arguments in favor of the existence of a mobility edge in the large system-size limit.
I. INTRODUCTION
Many-body localization (MBL) [1, 2] is a robust mech-anism that prevents reaching of thermal equilibrium byquantum many-body systems [3–5]. The phenomenon,originating from an interplay of interactions and disor-der [6–8], has been studied numerically in various mod-els: spin chains [9–12] that map onto spinless fermionicchains, spinful fermions [13–16] or bosons [17–19]. De-spite those efforts, a complete understanding of the tran-sition between ergodic and MBL phases is still lack-ing. While the recent works [20–24] suggest a Kosterlitz-Thouless scaling at the MBL transition, it became clearthat the exact diagonalization studies are subject tostrong finite size effects [25–27] that prevent one fromreaching unambiguous conclusions about the thermody-namic limit [28, 29].Alternatively, time evolution of large [30] (or even in-finite [31]) disordered many-body systems can be sim-ulated with tensor network algorithms. Reaching largetime scales, necessary to assess thermalization proper-ties [32] is challenging especially in the vicinity of tran-sition to ergodic phase. Nevertheless, such an approachallows to obtain estimates for critical disorder strengthfor large system sizes [30, 32], in quasiperiodic systems[33], or even beyond one spatial dimension [34, 35]. Anadvantage of such an approach is that it directly mimicsexperimental observations of MBL [36–43].Exact diagonalization study [12] of the Heisenbergchain showed that there exists, for a fixed disorderstrength W , a threshold energy such that eigenstatesabove (below) that energy are extended (localized). Thismeans that the system possesses the so-called many-bodymobility edge . The question of precise location of mobil-ity edges, even in a well understood case of single par-ticle Anderson localization [44, 45] is not fully resolved[46, 47]. Expectedly, the situation is even more com-plex for interacting systems. The existence of many-bodymobility edges was questioned in the work [48] which ar-gues that local fluctuations in presence of a mobility edgewould induce a global delocalization. That would leave ∗ [email protected] † [email protected] ‡ [email protected] FIG. 1. Phase diagram of random-field Heisenberg spin chain,disorder strength W on horizontal axis, rescaled energy (cid:15) onvertical axis. Background shows color-coded value of aver-age gap ratio r for system size L = 16. Solid lines show theposition of boundary between ergodic and MBL phases ob-tained in study of decay of density correlations in systemsof size L = 20 , , , , , L → ∞ . out only the two possibilities: either the states at all en-ergies are delocalized or localized. The mobility edgesfound in earlier works [12, 49, 50], were argued in [48] tobe indistinguishable from finite-size effects.The aim of our work is to study many-body mobilityedges at much larger system sizes than those availableto exact diagonalization in an attempt to verify conclu-sions of [48]. To that end, we employ Chebyshev poly-nomial expansion of the evolution operator [51–53] andthe time-dependent variational principle (TDVP) appliedto matrix product states (MPS) [54–57] to simulate timedynamics of random-field Heisenberg spin chain. Our ap-proach is, in spirit, similar to that of [58] (and used forbosons in [59]). However, instead of considering an injec-tion of controllable amount of energy into ground stateof the system, we consider time evolution of initial prod-uct states with specified average energies, exactly similarto what was done recently in spin quantum simulator[60]. Probing time decay of density correlation functionsallows us to estimate the critical disorder strength as a a r X i v : . [ c ond - m a t . d i s - nn ] J un function of energy of the initial state. Studying systemsof size up to L = 100, we perform a finite size scaling ofour results which provides arguments in favor of existenceof mobility edge even in large systems.The paper is organized as follows: Sec. II briefly intro-duces the disordered model as well numerical techniquesconsidered in this study. We present our main resultsregarding the estimation of many-body mobility edge viathe time evolution of large systems in Sec. III . Finally,we conclude in Sec. IV. II. THE MODEL AND METHOD
We consider 1D random-field Heisenberg (
XXZ ) spinchain with the Hamiltonian given by H = J L − (cid:88) i =1 (cid:0) S xi S xi +1 + S yi S yi +1 + S zi S zi +1 (cid:1) + L (cid:88) i =1 h i S zi (1)where S αi , α = x, y, z, are spin-1/2 matrices, J = 1 isfixed to be the unit of energy, and h i ∈ [ − W, W ] areindependent, uniformly distributed random variables. Inthis work, we consider open boundary conditions in theHamiltonian (1). The random-field Heisenberg spin chainhas been widely studied in the MBL context [12, 31, 61–67], which has made it the de facto standard model ofMBL studies.The transition between ergodic and MBL phases is re-flected in change of statistical properties of energy levelsof the system. A common approach is to consider thegap ratio r i = min { E i +2 − E i +1 ,E i +1 − E i } max { E i +2 − E i +1 ,E i +1 − E i } , where E i are theenergy eigenvalues of the system. Averaging the gap ra-tio over part of the spectrum of the system and overdisorder realizations, one obtains an average gap ratio r , which differentiates between level statistics of ergodicsystem [10, 68], well described by Gaussian orthogonalensemble of random matrices, for which r ≈ .
53 and be-tween Poissonian statistics of eigenvalues in MBL phase(for which r ≈ . r i are averaged over only a certainnumber of eigenvalues with energies close to a rescaled en-ergy (cid:15) = ( E − E min ) / ( E max − E min ), where E min ( E max )it the energy of the ground (highest excited) state. Sucha calculation of average gap ratio (supported with resultsfor other probes of localization) for random-field Heisen-berg spin chain reveals that the ergodic region has shapeof a characteristic lobe on the phase diagram in variablesof the rescaled energy (cid:15) and disorder strength W [12] .The average gap ratio, obtained in exact diagonalizationof random field Heisenberg spin chain of size L = 16, isplotted as a function of (cid:15) and W in the background ofFig. 1.To probe the transition between ergodic and MBLphases with time evolution, we propose the following pro- FIG. 2. Disorder averaged standard deviation ∆ (cid:15) ψ of rescaledenergy of the initial states as a function of disorder strength W for three exemplary rescaled energies (cid:15) . Left: system size L = 20, Right: system size: L = 50. tocol. We consider an initial state | ψ (cid:105) = | σ , . . . , σ L (cid:105) ,where σ i = ↑ , ↓ are chosen randomly with constraintthat the average rescaled energy (cid:15) ψ = ( (cid:104) ψ | H | ψ (cid:105) − E min ) / ( E max − E min ) of this state lies withing the range[ (cid:15) − δ(cid:15), (cid:15) + δ(cid:15) ] corresponding to a given rescaled energy (cid:15) , where δ(cid:15) is a small tolerance (we take δ(cid:15) = 0 . (cid:15) ψ for L (cid:54)
26 we find E max , E min withthe standard Lanczos algorithm [76]. For larger sys-tem sizes, E min and E max are calculated using densitymatrix renormalization group (DMRG) algorithm [77–81](see Appendix A).Subsequently, we calculate time evolved state | ψ ( t ) (cid:105) = e − iHt | ψ (cid:105) with the standard Chebyshev expansion of theevolution operator [53] for L (cid:54)
26. For larger systemsizes, we use the recently developed TDVP algorithm [54–57]. Technically, we follow [32, 82] and employ a hybridof two-site and one-site versions of TDVP [57, 83] (seeAppendix A).Our quantity of interest is the density correlation func-tion C ( t ) = D L − l (cid:88) i =1+ l (cid:104) ψ ( t ) | S zi | ψ ( t ) (cid:105)(cid:104) ψ | S zi | ψ (cid:105) , (2)where the constant D assures that C (0) = 1 and l > l = 2). The standard deviation of therescaled energy∆ (cid:15) ψ = (cid:16) (cid:104) ψ | (( H − E min ) / ( E max − E min ) − (cid:15) ψ ) | ψ (cid:105) (cid:17) / (3)is smaller than 0 . (cid:15) can be well probedby time evolution of the state | ψ (cid:105) and reflected, in par-ticular, by the density correlation function C ( t ). FIG. 3. Quench dynamics in disordered
XXZ spin chain.Density correlation function C ( t ) for system size L = 20and various disorder strengths W = 2 . , ... (cid:15) = 0 .
5, power-law fits C ( t ) ∝ t − β for t ∈ [100 , III. QUENCH DYNAMICSA. Disorder strength dependence
Fig. 3 shows the density correlation functions C ( t ) ob-tained for the random-field Heisenberg spin chain of afixed size L = 20, for rescaled energy (cid:15) = 0 . W = 2 . C ( t ) t →∞ → W = 5, a non-zero sta-tionary value of the correlation function C ( t ) t →∞ → c > imbalance [36], quantity analo-gous to the density correlation function – for quantitativecomparison of the two quantities see Appendix B.At large times ( t > C ( t ) ∝ t − β .Griffiths rare regions are one possible explanation of thisbehavior [86]. However, it was shown experimentally andnumerically that time dynamics in quasiperiodic poten-tials, where Griffiths regions are necessarily absent, haveanalogous features [64, 87, 88]. Regardless of the origin ofthe power law decay of the correlation function, the dis-order strength dependence of the exponent β can be usedto locate the onset of ergodicity breaking in the system.The exponent β governing the decay of the densitycorrelation function is shown in Fig. 4. Let us first con-centrate on the results in the middle of the spectrum( (cid:15) = 0 . W , the exponent decreases exponentially with W with agood approximation β ∝ e − W/ Ω . The large number of FIG. 4. The exponent β , obtained in fitting the density corre-lation function C ( t ) with an algebraic decay a t − β in interval t ∈ [100 , W .The errorbars represent the 1 σ errors of the fitting obtainedfrom statistical resampling of disorder realizations. The sys-tem size is L = 20, results for various rescaled energies (cid:15) ofthe initial state are shown. The dashed line shows the cut-offexponent β . disorder realizations (10000) used in calculation of C ( t )allows us to see that even at the large disorder strength W = 5 the exponent β = 4 . · − is non-vanishing. Ifthe power-law decay C ( t ) = a t − β prevailed for t → ∞ ,the density correlation function would vanish in the long-time limit and the system would be ergodic. This, how-ever, does not happen for L = 20, as after the so-calledHeisenberg time t H discreteness of spectrum manifestsitself in saturation of the correlation functions [89–92],so that one would observe C ( t ) t →∞ → c > W = 5and L = 20. The Heisenberg time t H increases exponen-tially with the system size L . This illustrates a difficultyin locating the MBL transition using time dynamics oflarge systems on time scales of few hundred J − accessi-ble to tensor network methods (or to current experimentswith e.g., ultra-cold atoms): one cannot predict whethera slow decay of correlation functions governed by an ex-ponent β (cid:28) t ∈ [100 , C ( t ) t →∞ → β must be vanishing within er-ror bars to be compatible with saturation of correlationfunctions in the long time limit. The drawback of thiscriterion is that the error bar of β depends on the numberof disorder realizations used in calculation of the corre-lation function. Therefore, we introduce a cut-off β :disorder strength W C ( L ) for which β = β is regarded asdisorder strength for transition to MBL phase at systemsize L . Exact diagonalization results show that: i) col-lapse of data for L (cid:54)
22 gives a critical disorder strength W C ≈ .
7; ii) the similar values W C ≈ . W C ≈ . W C = 3 . L = 20 [29]. The obtained results for FIG. 5. Time evolution of density correlation function C ( t )for rescaled energy of initial state (cid:15) = 0 . (cid:15) = 0 .
8) and dis-order strength W = 3 . W = 2 .
7) in panel right (left). Thesystem size L varies from 22 to 100. The dashed lines denotepower-law fits C ( t ) = a t − β in the t ∈ [100 , β at L = 20 and (cid:15) = 0 . β = 0 .
014 is consistent with the above estimates forthe critical disorder strength obtained from exact diago-nalizations. Consequently, throughout this work, we use β = 0 .
014 as a threshold value which separates ergodicand MBL regimes for all system sizes and energies of theinitial state we consider.The values of β presented in Fig. 4 show that the in-crease of disorder strength W slows down the dynamicsmore severely for rescaled energies of initial state differentthan (cid:15) = 0 .
5. Notably, the exponent β decreases expo-nentially with W : β ∝ e − W/ Ω (where Ω is a constant) ina wide regime of disorder strengths. This resembles thescaling of Thouless time t T h ∝ e − W/W observed in exactdiagonalization data in [25]. B. System-size dependence
Density correlation function C ( t ) for larger systemsizes are shown for two exemplary pairs of disorderstrength W and initial rescaled energy (cid:15) in Fig. 5. Thedecay of C ( t ) at large times is well fitted by an algebraicdependence C ( t ) ∝ t − β . The exponents β obtained inthe fitting of power-law decay to C ( t ) are shown for twoexemplary values of the rescaled energy (cid:15) of the initialstates in Fig. 6. For a given disorder strength W , weobserve a clear increase of β with increasing system size.Interestingly, the shift is, to a good approximation, uni-form for all disorder strengths so that the exponentialdecrease β ∝ e − W/ Ω (at sufficiently large W ) is observedfor all considered system sizes. Let us mention here thatwe consider 400 realizations of disorder for L = 34 , L = 100 for each values of (cid:15) and W. For small system sizes ( L = 20 , ,
26) we considerbetween 10000 and 500 disorder realizations.We obtain estimates for disorder strength W C ( L ) fortransition to MBL phase by finding the crossings of β ( W )curve for given system size L with the β = β line. Re-sults of this procedure are shown in Fig. 7. The disor-der strength W C ( L ) depend, within the estimated error FIG. 6. The exponent β obtained in fitting the density corre-lation function C ( t ) with an algebraic decay a t − β in interval t ∈ [100 , (cid:15) = 0 . (cid:15) = 0 .
5) ofthe initial state shown in the left (right) panel. Data shownfor system sizes L = 20 , , , , σ errors of β obtained in resampling of disorder realizations.The dashed lines show the cut-off exponent β .FIG. 7. Disorder strength W C ( L ) for which decay of correla-tion function is governed by power-law with β = β plotted asfunction of 1 /L where L is the system size. Results shown forvarious rescaled energies (cid:15) of initial state. Available data arefitted with linear functions W C (1 /L ) = a/L + W C ( ∞ ) whichallow extrapolation to L → ∞ . bars, linearly on the inverse of the system size L withclear growth of W C ( L ) as the system size increases. Onone hand, this trend allows us, by means of a linear fit W C (1 /L ) = A/L + W C ( ∞ ), to extrapolate the results to L → ∞ and to obtain the estimate of critical disorderstrength W C ( ∞ ) for transition to MBL phase.On the other hand, we observe that the slopes A aresimilar for all of the considered rescaled energies of theinitial state. Thus, the shape of the boundary betweenergodic and MBL regimes observed for L = 20 does notchange considerably when the system size is increased.This is visible in Fig. 1. The points for various systemsizes L are precisely the values of W C ( L ) obtained fromthe condition β = β . The characteristic shape of thelobe does not change when the system sizes increasesfrom L = 20 to L = 100 and is preserved even afterthe extrapolation to L → ∞ . Therefore, there exists acertain range of disorder strengths such that the statesfor (cid:15) < (cid:15) LME are localized, states for (cid:15)
LME < (cid:15) < (cid:15)
UME are extended and states for (cid:15) > (cid:15)
UME are again localized.Thus, our results indicate that the system indeed pos-sesses a many-body mobility edge in the thermodynamiclimit.
IV. DISCUSSION AND OUTLOOK
Chebyshev polynomial expansion of the time evolutionoperator and the TDVP method applied to MPS allowedus to study the problem of energy dependence of thetransition between ergodic and MBL phases in large dis-ordered quantum spin chains. Introducing a cut-off valueof exponent β of power-law decay in time of density cor-relation function, we were able to probe the transitionfor different rescaled energies of the initial state. Forsmall system-sizes (e.g., L = 20), our approach gives re-sults consistent with exact diagonalization. Importantly,our method allows to consider much larger system sizes( L = 100) for which it predicts an existence of a mobilityedge.The disorder strength W C ( L ) is a lower bound on thetransition to MBL phase: the residual decay of densitycorrelations with exponent β is insufficient to restore theuniform density profile for system size L = 20, but it ispossible that it leads to an eventual decay of correlationfunction for larger system sizes.The protocol we considered is, in principle, experimen-tally realizable. Our results can be verified experimen-tally if the setup of [60] was scaled to larger system sizes.Many-body mobility edge arises also in disordered Bose-Hubbard models [93]. It can be probed by a quenchprotocol analogous to the one considered in this work.Since the bosonic models allow for occupations in eachsite larger than unity, density wave-like states that areeasier to obtain in an experiment with ultra-cold atomscan be use to probe the many-body mobility edge [59, 93]. ACKNOWLEDGMENTS
Support of the Polish National Science Cen-tre via grants Unisono 2017/25/Z/ST2/03029(T.C.) (under QTFLAG Quantera collaboration),Opus 2015/19/B/ST2/01028 (P.S.), and Opus2019/35/B/ST2/00034 (J.Z.) is acknowledged. P.S.thanks Polish National Science Centre for an additionalsupport via Etiuda programme 2018/28/T/ST2/00401.The partial support by PL-Grid Infrastructure isalso acknowledged. The MPS-based techniqueshave been implemented using ITensor library v2( https://itensor.org ). FIG. 8. The exponent β C obtained in fitting the density corre-lation function C ( t ) with an algebraic decay a t − β in interval t ∈ [100 , (cid:15) = 0 . β Imb obtained from the fitting of I ( t )for N´eel state in the same time interval. Error bars dictate1 σ errors in β obtained from statistical resampling procedure.Left: system size L = 26, Right: system size L = 50. Appendix A: Brief descriptions of the MPScalculations
For large system sizes, i.e.,
L >
26, we use matrixproduct states (MPS) ansatz [80, 81] to represent thequantum state of the disordered Heisenberg chain. Theenergies of the ground state ( E min ) and highest excitedstate ( E max ) are obtained using standard density matrixrenormalization group method (DMRG) [77–81] . E min can simply be found from the ground-state search of H ofEq. (1) using DMRG, while that of − H gives us − E max .In these calculations, we keep truncation error due tosingular value decomposition (SVD) below 10 − , i.e., wediscard states corresponding to smallest singular values λ n that satisfy (cid:80) n ∈ discarded λ n (cid:80) n ∈ all λ n < − .For time evolution large systems, i.e., L >
26, weemploy MPS based time-dependent variational principle(TDVP) method [54–57]. Specifically, we use ‘hybrid’scheme of TDVP [32, 57, 82, 83] with step-size δt = 0 . χ max . One-site variant of TDVPis employed in the subsequent dynamics. In our calcula-tions, we fix χ max to be 384, which gives reliable resultsup to t = 500. For a more detailed analysis of the conver-gence of TDVP in the disordered Heisenberg chain, see[32]. Appendix B: Imbalance vs density correlationfunction
The N´eel state | ψ N´eel (cid:105) = | ↑ , ↓ , ↑ , . . . (cid:105) (B1)was used in previous works [30, 32] to locate the transi-tion to MBL phase. It is worth to mention that densitywave states, analogous to the N´eel state were used in ex-periment with ultra-cold fermions [36] that demonstratedsignatures of MBL transition.Imbalance I ( t ) is the common name for the densitycorrelation function (2) calculated for the N´eel state.The imbalance behaves similarly to the density corre-lation function: initially it decreases and oscillates, sub-sequently it decays algebraically in time. Fig. 8 showscomparison of exponents β C and β Imb that govern thepower-law decay of C ( t ) correlation function for rescaledenergy of initial state (cid:15) = 0 . β Imb areconsistently bigger than β C . This fact may be tracedback to the number of states coupled by the Hamiltonian(1) to state | ψ N´eel (cid:105) which is higher than the average num-ber of states coupled to a typical Fock state with (cid:15) = 0 . β C < β Imb justifies the cut-off beta value β Imb = 0 .
02 for decay of imbalance in N´eel state assumedin [32]. [1] I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Phys.Rev. Lett. , 206603 (2005).[2] D. Basko, I. Aleiner, and B. Altschuler, Ann. Phys. (NY) , 1126 (2006).[3] J. M. Deutsch, Phys. Rev. A , 2046 (1991).[4] M. Srednicki, Phys. Rev. E , 888 (1994).[5] M. Rigol, V. Dunjko, and M. Olshanii, Nature , 854EP (2008).[6] R. Nandkishore and D. A. Huse, Ann. Rev. Cond. Mat.Phys. , 15 (2015).[7] F. Alet and N. Laflorencie, Comptes Rendus Physique , 498 (2018).[8] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev.Mod. Phys. , 021001 (2019).[9] L. F. Santos, G. Rigolin, and C. O. Escobar, Phys. Rev.A , 042304 (2004).[10] V. Oganesyan and D. A. Huse, Phys. Rev. B , 155111(2007).[11] A. Pal and D. A. Huse, Phys. Rev. B , 174411 (2010).[12] D. J. Luitz, N. Laflorencie, and F. Alet, Phys. Rev. B , 081103 (2015).[13] R. Mondaini and M. Rigol, Phys. Rev. A , 041601(2015).[14] P. Prelovˇsek, O. S. Bariˇsi´c, and M. ˇZnidariˇc, Phys. Rev.B , 241104 (2016).[15] J. Zakrzewski and D. Delande, Phys. Rev. B , 014203(2018).[16] M. Kozarzewski, P. Prelovˇsek, and M. Mierzejewski,Phys. Rev. Lett. , 246602 (2018).[17] P. Sierant, D. Delande, and J. Zakrzewski, Phys. Rev.A , 021601 (2017).[18] T. Orell, A. A. Michailidis, M. Serbyn, and M. Silveri,Phys. Rev. B , 134504 (2019).[19] M. Hopjan and F. Heidrich-Meisner, “Many-body local-ization from a one-particle perspective in the disordered1d bose-hubbard model,” arXiv:1912.09443.[20] A. Goremykina, R. Vasseur, and M. Serbyn, Phys. Rev.Lett. , 040601 (2019).[21] A. Morningstar and D. A. Huse, Phys. Rev. B , 224205(2019).[22] P. T. Dumitrescu, A. Goremykina, S. A. Parameswaran,M. Serbyn, and R. Vasseur, Phys. Rev. B , 094205(2019).[23] N. Laflorencie, G. Lemari´e, and N. Mac´e, “Chain break-ing and Kosterlitz-Thouless scaling at the many-body lo-calization transition,” arXiv:2004.02861.[24] J. ˇSuntajs, J. Bonˇca, T. Prosen, and L. Vidmar, “Ergod-icity breaking transition in finite disordered spin chains,” (2020), arXiv:2004.01719.[25] J. ˇSuntajs, J. Bonˇca, T. Prosen, and L. Vidmar,arXiv:1905.06345.[26] P. Sierant, D. Delande, and J. Zakrzewski, Phys. Rev.Lett. , 186601 (2020).[27] D. A. Abanin, J. H. Bardarson, G. D. Tomasi,S. Gopalakrishnan, V. Khemani, S. A. Parameswaran,F. Pollmann, A. C. Potter, M. Serbyn, and R. Vasseur,arXiv:1911.04501.[28] R. K. Panda, A. Scardicchio, M. Schulz, S. R. Taylor,and M. ˇZnidariˇc, EPL (Europhysics Letters) , 67003(2020).[29] P. Sierant, M. Lewenstein, and J. Zakrzewski, “Polyno-mially filtered exact diagonalization approach to many-body localization,” arXiv:2005.09534.[30] E. V. H. Doggen, F. Schindler, K. S. Tikhonov, A. D.Mirlin, T. Neupert, D. G. Polyakov, and I. V. Gornyi,Phys. Rev. B , 174202 (2018).[31] T. Enss, F. Andraschko, and J. Sirker, Phys. Rev. B ,045121 (2017).[32] T. Chanda, P. Sierant, and J. Zakrzewski, Phys. Rev. B , 035148 (2020).[33] E. V. H. Doggen and A. D. Mirlin, Phys. Rev. B ,104203 (2019).[34] C. Hubig and J. I. Cirac, SciPost Phys. , 31 (2019).[35] E. V. H. Doggen, I. V. Gornyi, A. D. Mirlin, and D. G.Polyakov, “Slow many-body delocalization beyond onedimension,” arXiv:2002.07635.[36] 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).[37] P. Bordia, H. L¨uschen, U. Schneider, M. Knap, andI. Bloch, Nature Physics , 460 EP (2017), article.[38] T. Kohlert, S. Scherg, X. Li, H. P. L¨uschen,S. Das Sarma, I. Bloch, and M. Aidelsburger, Phys. Rev.Lett. , 170403 (2019).[39] 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).[40] M. Rispoli, A. Lukin, R. Schittko, S. Kim, M. E. Tai,J. L´eonard, and M. Greiner, Nature , 385 (2019).[41] 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).[42] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess,P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, NaturePhysics , 907 (2016).[43] 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).[44] P. W. Anderson, Phys. Rev. , 1492 (1958).[45] F. Evers and A. D. Mirlin, Rev. Mod. Phys. , 1355(2008).[46] D. Delande and G. Orso, Phys. Rev. Lett. , 060601(2014).[47] M. Pasek, G. Orso, and D. Delande, Phys. Rev. Lett. , 170403 (2017).[48] W. De Roeck, F. Huveneers, M. M¨uller, and M. Schiulaz,Phys. Rev. B , 014203 (2016).[49] J. A. Kj¨all, J. H. Bardarson, and F. Pollmann, Phys.Rev. Lett. , 107204 (2014).[50] I. Mondragon-Shem, A. Pal, T. L. Hughes, and C. R.Laumann, Phys. Rev. B , 064203 (2015).[51] H. Tal-Ezer and R. Kosloff, J. Chem. Phys. , 3967(1984).[52] C. Leforestier, R. Bisseling, C. Cerjan, M. Feit, R. Fries-ner, A. Guldberg, A. Hammerich, G. Jolicard, W. Kar-rlein, H.-D. Meyer, N. Lipkin, O. Roncero, andR. Kosloff, J. Comput. Phys. , 59 (1991).[53] H. Fehske and R. Schneider, Computational many-particle physics (Springer, Germany, 2008).[54] J. Haegeman, J. I. Cirac, T. J. Osborne, I. Piˇzorn, H. Ver-schelde, and F. Verstraete, Phys. Rev. Lett. , 070601(2011).[55] T. Koffel, M. Lewenstein, and L. Tagliacozzo, Phys. Rev.Lett. , 267203 (2012).[56] J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken,and F. Verstraete, Phys. Rev. B , 165116 (2016).[57] S. Paeckel, T. K¨ohler, A. Swoboda, S. R. Manmana,U. Schollw¨ock, and C. Hubig, Annals of Physics ,167998 (2019).[58] P. Naldesi, E. Ercolessi, and T. Roscilde, SciPost Phys. , 010 (2016).[59] R. Yao and J. Zakrzewski, “Many-body localization inBose-Hubbard model: evidence for the mobility edge,”arXiv:2002.00381.[60] Q. Guo, C. Cheng, Z.-H. Sun, Z. Song, H. Li, Z. Wang,W. Ren, H. Dong, D. Zheng, Y.-R. Zhang, R. Mondaini,H. Fan, and H. Wang, “Observation of energy resolvedmany-body localization,” arXiv:1912.02818.[61] T. C. Berkelbach and D. R. Reichman, Phys. Rev. B ,224429 (2010).[62] K. Agarwal, S. Gopalakrishnan, M. Knap, M. M¨uller,and E. Demler, Phys. Rev. Lett. , 160401 (2015).[63] S. Bera, H. Schomerus, F. Heidrich-Meisner, and J. H.Bardarson, Phys. Rev. Lett. , 046603 (2015).[64] S. Bera, G. De Tomasi, F. Weiner, and F. Evers, Phys.Rev. Lett. , 196801 (2017).[65] L. Herviou, S. Bera, and J. H. Bardarson, Phys. Rev. B , 134205 (2019).[66] L. Colmenarez, P. A. McClarty, M. Haque, and D. J.Luitz, arXiv:1906.10701.[67] P. Sierant and J. Zakrzewski, Phys. Rev. B , 104201(2020).[68] Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux,Phys. Rev. Lett. , 084101 (2013).[69] M. Serbyn, Z. Papi´c, and D. A. Abanin, Phys. Rev. Lett. , 127201 (2013).[70] D. A. Huse, R. Nandkishore, and V. Oganesyan, Phys.Rev. B , 174202 (2014).[71] V. Ros, M. Mueller, and A. Scardicchio, Nuclear PhysicsB , 420 (2015).[72] J. Z. Imbrie, Phys. Rev. Lett. , 027201 (2016).[73] T. B. Wahl, A. Pal, and S. H. Simon, Phys. Rev. X ,021018 (2017).[74] M. Mierzejewski, M. Kozarzewski, and P. Prelovˇsek,Phys. Rev. B , 064204 (2018).[75] S. J. Thomson and M. Schir´o, Phys. Rev. B , 060201(2018).[76] C. Lanczos, Journal of Research of the National Bureauof Standards (1950).[77] S. R. White, Phys. Rev. Lett. , 2863 (1992).[78] S. R. White, Phys. Rev. B , 10345 (1993).[79] U. Schollw¨ock, Rev. Mod. Phys. , 259 (2005).[80] U. Schollw¨ock, Annals of Physics , 96 (2011).[81] R. Or´us, Annals of Physics , 117 (2014).[82] T. Chanda, J. Zakrzewski, M. Lewenstein, and L. Tagli-acozzo, Phys. Rev. Lett. , 180602 (2020).[83] S. Goto and I. Danshita, Phys. Rev. B , 054307 (2019).[84] D. J. Luitz, N. Laflorencie, and F. Alet, Phys. Rev. B , 060201 (2016).[85] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol,Advances in Physics , 239 (2016).[86] K. Agarwal, E. Altman, E. Demler, S. Gopalakrishnan,D. A. Huse, and M. Knap, Annalen der Physik , (2016).[87] H. P. L¨uschen, P. Bordia, S. Scherg, F. Alet, E. Altman,U. Schneider, and I. Bloch, Phys. Rev. Lett. , 260401(2017).[88] F. Weiner, F. Evers, and S. Bera, Phys. Rev. B ,104204 (2019).[89] E. J. Torres-Herrera and L. F. Santos, Phys. Rev. B ,014208 (2015).[90] E. J. Torres-Herrera and L. F. Santos, PhilosophicalTransactions of the Royal Society A: Mathematical,Physical and Engineering Sciences , 20160434 (2017).[91] E. J. Torres-Herrera, A. M. Garc´ıa-Garc´ıa, and L. F.Santos, Phys. Rev. B , 060303 (2018).[92] M. Schiulaz, E. J. Torres-Herrera, and L. F. Santos,Phys. Rev. B , 174313 (2019).[93] P. Sierant and J. Zakrzewski, New Journal of Physics20