Characterizing many-body localization via exact disorder-averaged quantum noise
CCharacterizing many-body localizationvia exact disorder-averaged quantum noise
Michael Sonner, ∗ Alessio Lerose, ∗ and Dmitry A. Abanin Department of Theoretical Physics, University of Geneva,Quai Ernest-Ansermet 24, 1205 Geneva, Switzerland (Dated: 18th December 2020)Many-body localized (MBL) phases of disordered quantum many-particle systems have a numberof unique properties, including failure to act as a thermal bath and protection of quantum coherence.Studying MBL is complicated by the effects of rare ergodic regions, necessitating large system sizesand averaging over many disorder configurations. Here, building on the Feynman-Vernon theoryof quantum baths, we characterize the quantum noise that a disordered spin system exerts on itsparts via an influence matrix (IM). In this approach, disorder averaging is implemented exactly, andthe thermodynamic-limit IM obeys a self-consistency equation. Viewed as a wavefunction in thespace of trajectories of an individual spin, the IM exhibits slow scaling of temporal entanglementin the MBL phase. This enables efficient matrix-product-states computations to obtain temporalcorrelations, providing a benchmark for quantum simulations of non-equilibrium matter. The IMquantum noise formulation provides an alternative starting point for novel rigorous studies of MBL.
Introduction —
Many-body localization (MBL)in strongly disordered interacting quantum systemsrepresents one of the rare known examples of a genu-inely non-ergodic phase of quantum matter [1–3].The phenomenology of the MBL phase [4–7], its per-sistence in periodically driven systems [8–10], andnew phases of matter enabled by MBL [11] are be-ing intensely investigated.Unlike in conventional phase transitions, thermalensembles bear no signatures of the MBL trans-ition. Localization effects instead manifest them-selves in the properties of individual highly-excitedenergy eigenstates, as well as in the coherent, far-from-equilibrium dynamics, which differ drasticallyfrom those in the thermalizing phase. In particu-lar, MBL eigenstates exhibit boundary-law entan-glement scaling [4, 12], similar to ground states ofgapped systems [13]. Local memory of the initialstate is preserved, owing to the emergence of localintegrals of motion [4–7] – a signature that has beenwidely used to diagnose MBL in quantum quenchexperiments with ultracold atoms [14–16], trappedions [17], and superconducting qubits [18].The simultaneous presence of disorder and inter-actions, combined with a necessity to describe highlyexcited eigenstates, poses a major challenge for thetheoretical description of MBL. Rigorous perturbat-ive approaches are difficult, as they require treat-ing disorder probabilistically [7] and incorporatingthe effects of rare regions [19]. Exact diagonaliza-tion (ED) studies [20–23] give access to excited ei-genstates, but are limited in system size. Tensor-network approaches to MBL, relying on low real-space entanglement of MBL states, have been usedto approximately construct eigenstates [24–27] and simulate non-equilibrium dynamics [28] beyond thereach of ED. While such investigations gave insightsinto MBL physics, their power is limited by finite-size and/or -time effects, and by the practical costof sampling disorder configurations.In this Article, we introduce a method for de-scribing disorder-averaged dynamical properties ofmany-body systems, and apply it to a periodicallydriven (Floquet) spin model of MBL. We charac-terize the system’s dynamics by its influence matrix (IM) [29], inspired by the Feynman-Vernon theory ofa quantum particle interacting with a bath of har-monic oscillators [30, 31]. We use the IM to char-acterize how a quantum many-body system affectsthe time evolution of its local subsystems, i.e., howit acts as a “bath” on itself. An IM contains full in-formation on temporal correlations of local operat-ors, and can be viewed as a kind of generating func-tional of the system’s self-induced quantum noise.This formalism is advantageous for MBL, as it nat-urally accommodates the thermodynamic limit and exact disorder averaging.In one-dimensional homogeneous systems, the IMcan be found from a linear self-consistency equa-tion [29, 32]. While individual realizations of dis-ordered systems have IMs that depend on position,averaging over a translationally invariant disorderdistribution leads to a translationally invariant IM.Disorder-averaging, however, leads to terms whichare non-local in time in the self-consistency equation.Below we argue that, in spite of this nonlocality, thedisorder-averaged IM of a MBL system is charac-terized by low temporal entanglement . This opensthe door to efficient matrix-product states (MPS)methods for computing disorder-averaged dynam- a r X i v : . [ c ond - m a t . d i s - nn ] D ec ≡≡≡ e ihσ α δ αβγδ α β δ γ α β δ γ e iJσ α σ β cos g δ σ γ , σ δ + i sin g δ σ γ ,− σ δ a)b) c) d)e) p ℐ L , p ℐ R , p = } ̂ F space } ρ ̂ T p +2 time ⋯ + p p − p ⋯ + p } ̂ F * t . . . . S . . . . g t .
98 0 .
99 1 . ̂𝒯 ⟨ℐ⟩ = ⟨ℐ⟩ Figure 1. a) Circuit representation of a local temporal correlation function in the Keldysh path integral represent-ation, Eq. (2), with circuit elements for model (1) defined in panel b). The two layers reflect forward and backwardpropagation of the system. Panel a) pictures Eq. (3). The blue shaded tensor represents the transfer matrix, definedin Eq. (4). Panel c) illustrates Eq. (5) (self-consistency equation). The green shaded tensor corresponds to non-localin time “interactions” arising from exact disorder averaging. d) Scaling of temporal entanglement of the IM for J = g = 0 . , . , . , . , . (DMRG) and J = g = 0 . , . (ED), bottom to top. e) Expectation value of thetransfer matrix applied to the MPS with bond dimension χ = 128 obtained from DMRG. ical properties. As a first application of the method,we compute the dynamical correlation functions ofa MBL system up to long times.The IM undergoes a drastic change when the sys-tem transitions from the ergodic to the MBL phase;MBL IMs are characterized by persistent quantum-interference effects, which express the fact that MBLsystems are not efficient thermal baths. This can beviewed as the emergence of temporal long-range or-der in the statistical ensemble of local trajectoriesgoverned by the IM. Below, we will use an MPS ap-proach to demonstrate this phenomenon [33]. Model —
For concreteness, in this Article we willfocus on the disordered kicked Ising chain (KIC),which provides a Floquet model of MBL [34, 35].Time evolution of this system is governed by re-peated applications of the Floquet operator, ˆ F = exp (cid:18) i X j g ˆ σ xj (cid:19) exp (cid:18) i X j J ˆ σ zj ˆ σ zj +1 + h j ˆ σ zj (cid:19) , (1)where ˆ σ αj , α = x, y, z , are Pauli matrices acting onsite j ∈ Z of a linear chain. The phases h j are in- dependently drawn from a uniform distribution in [ − π, π ) . The Hamiltonian version of this model,obtained by substituting J, g, h j τ J, τ g, τ h j inEq. (1) and taking the continuous-time limit τ → ,is similar to the model where MBL was rigorously es-tablished in the regime g (cid:28) [7]. MBL behavior isknown to persist in the same regime for finite drivingperiod τ > [8, 9, 34, 36, 37]. Setting J = g in Eq.(1), ED studies indicate that the MBL phase extendsto | g | < g ∗ ≈ . [35]. For weaker disorder strength | J | = | g | > g ∗ the disordered KIC is ergodic. Inparticular, at the self dual points | J | = | g | = π/ signatures of chaotic behavior in spectral correla-tions have been obtained [38]. We note that otherkicked Floquet models of MBL have been investig-ated [8, 37, 39]. Disorder-averaged influence matrix —
The influ-ence matrix encodes the full set of temporal cor-relations of local operators [29, 30]. To illustratethis, we consider the dynamical structure factor h O p ( t ) O p (0) i = Tr ( ˆ F − t O p ˆ F t O p ρ ) of a local observ-able O p = ⊗ · · · ⊗ ⊗ O ⊗ ⊗ · · · ⊗ acting onspin p , using Keldysh path integral representation,graphically illustrated in Fig. 1(a,b), h O p ( t ) O p (0) i = X { σ τj } , { ¯ σ τj } [ O p ] { ¯ σ t } , { σ t } (cid:2) O p ρ (cid:3) { σ } , { ¯ σ } Y j t − Y τ =0 W σ τ +1 j σ τj W ∗ ¯ σ τ +1 j ¯ σ τj e iJ ( σ τj σ τj +1 − ¯ σ τj ¯ σ τj +1 ) + ih j ( σ τj − ¯ σ τj ) (2)where W σ σ = h σ | e ig ˆ σ x | σ i . We say that configura-tions σ τj associated with the operator ˆ F t are on the forward time path and those ¯ σ τj associated with ˆ F − t on the backward path. The summation is over allspin trajectories { σ τj = ± } , { ¯ σ τj = ± } extending from time τ = 0 to τ = t . Assuming that the initialdensity matrix ρ = ⊗ j ρ j is a product operator, weformally perform the summation over trajectories ofall spins on the left { σ τj
p } , { ¯ σ τj>p } of spin p in Eq. (2), obtaining h O p ( t ) O p (0) i = X { σ τ } , { ¯ σ τ } I L,p [ { σ τ , ¯ σ τ } ] (cid:18) [ O ] ¯ σ t σ t t − Y τ =0 W σ τ +1 σ τ W ∗ ¯ σ τ +1 ¯ σ τ e ih p ( σ τ − ¯ σ τ ) (cid:2) Oρ p (cid:3) σ ¯ σ (cid:19) I R,p [ { σ τ , ¯ σ τ } ] (3)where we have denoted the result of the summationover spins on the left (right) as the left (right) in-fluence matrix I L ( R ) ,p acting on the forward andbackward trajectories { σ τ } , { ¯ σ τ } , ≤ τ ≤ t − ofspin p only. These objects, graphically representedin Fig. 1(a), capture the influence the rest of thesystem has on the dynamics of spin p . Clearly, ex-pression (3) can be straightforwardly generalized toarbitrary time-ordered temporal correlations of local operators h O ( n ) p ( t n ) · · · O (1) p ( t ) i .In a chain (or, more generally, in loop-free geomet-ries), an influence matrix can be recursively com-puted from influence matrices of smaller subsys-tems, as illustrated from in Fig. 1(a,c). For a fi-nite system of N spins, this can be concisely ex-pressed by introducing transfer matrices ˆ T j actingalong the space direction, i.e., I L,p = (cid:16)Q p − j =1 ˆ T j (cid:17) I and I R,p = (cid:16)Q p +1 j = N ˆ T j (cid:17) I , with [ T j ] { σ, ¯ σ } , { s, ¯ s } = (cid:18) t − Y τ =0 e iJ ( σ τ s τ − ¯ σ τ ¯ s τ ) (cid:19) e ih j P t − τ =0 ( s τ − ¯ s τ ) (cid:18) δ s t ¯ s t t − Y τ =0 W s τ +1 s τ W ∗ ¯ s τ +1 ¯ s τ (cid:2) ρ j (cid:3) s ¯ s (cid:19) , (4)and open boundary condition I [ { s, ¯ s } ] ≡ .We are interested in disorder-averaged temporalcorrelations such as hh O p ( t ) O p (0) ii . Since the ran-dom field h j is uncorrelated and equally distributedat each lattice site, the averaged transfer matrices h ˆ T j i = ˆ T are translationally invariant, provided theinitial density matrices ρ j are the same. Unitarity oftime evolution for a periodic chain of arbitrary size N gives Tr( ˆ T N ) ≡ hh ii = 1 ; hence, ˆ T has a singlenon-vanishing and non-degenerate eigenvalue equalto . The “bulk” disorder-averaged influence matrix h I i = ˆ T ‘ I , with ‘ > t , is thus real [40], and canbe identified with the single eigenvector of ˆ T , i.e.,by the characteristic self-consistency equation ˆ T h I i = h I i . (5)Due to the absence of correlations in the initial stateand the strictly linear light-cone effect in this Flo-quet model, finite-size effects are present only at adistance ‘ ≤ t from the boundaries [41].Disorder averaging in model (1) can be carriedout explicitly [38]: since disorder is random in spacebut uniform in time, averaging introduces non-local in time “interactions”. Indeed, integrating over h j in Eq. (4) transforms the corresponding phase terminto a constraint δ (cid:0) P t − τ =0 ( s τ − ¯ s τ ) (cid:1) ≡ δ M,M , whichcouples the configurations of the considered spinbetween all times. Such a term cancels interferencebetween forward and backward trajectories with dif-ferent total magnetization M = M , and is respons-ible for the development of long-range temporal cor-relations in the MBL phase. MPS approach and temporal entanglement —
Togain an insight into the structure of the IM in theMBL phase, we approximate the solution of the self-consistency equation (5) by an MPS ansatz with amaximal bond dimension χ . The reliability of thisapproximation depends on the amount of “bipartiteentanglement” in the IM interpreted as a many-body“wavefunction” in the t -dimensional space of for-ward/backward spin trajectories. This temporal en-tanglement (TE) may be considered as a quantifierof the system’s dynamics computational complexity.Previously, we have shown [29] that the max-imal (half-chain) von Neumann entanglement en-tropy S ( t/ of the infinite-temperature folded IM —obtained by considering σ τ , ¯ σ τ as a single 4-dit, cf.Fig. 1(a-c) — is exactly zero when | J | = | g | = π/ for any h j [42] thanks to the fact that the systemacts as a perfectly dephasing (PD) Markovian bath.However, detuning away from these PD points, theentanglement S ( t/ acquires “volume-law” scalingwith t , albeit with a prefactor that vanishes as thosepoints are approached. This scaling is expected tobe a generic feature of thermalizing phases, and gen-erally prevents MPS description from being efficientat long times, except near “special” points [29].In contrast, we find that MPS methods remain ef-ficient in the MBL phase up to long times, thanksto the slow scaling of TE. To implement MPS al-gorithms, we work in the folded picture [29, 32, 43].The disorder-averaged transfer matrix ˆ T is repres-ented as a matrix product operator (MPO), whichconsists of the diagonal two-site matrices W , theglobal projection operator δ M,M originating from av-eraging over the random phase h , and the factor-ized operator dependent on the interaction strength J [see Eq. (4)]. The first and last operators canbe expressed as MPOs with bond dimension and , respectively. The projection operator can be ex-pressed as a MPO with maximal bond dimension t [44]. Thus, the maximal bond dimension of ˆ T is t .To find the IM one can iteratively apply the trans-fer matrix to the boundary vector I . This ap-proach was used to obtain the ED results, but itdoes not yield a good approximation of the optimalMPS at a fixed bond dimension, mainly due to therelatively large bond dimension of the MPO. To mit-igate this problem we refine the MPS afterwardsby using the density matrix renormalization group(DMRG) algorithm [45]. We estimated the qualityof the MPS representation using several metrics [44],including the proximity of the eigenvalue obtainedfrom DMRG to the exact value (see Fig. 1(e)).We have used a combination of the MPS methodand ED to analyze the infinite-temperature IM’sTE across the MBL transition in model (1), seeFig. 1(d). In the ergodic phase, TE increases fastwith t , similarly to generic non-disordered thermalsystems [29], which restricts us to the ED approachand thus limited time. However, this behaviordrastically changes in the MBL phase, where TE ex-hibits an initial rise followed by a crossover to a veryslow growth. Interestingly, the crossover occurs atlonger times at smaller values of g . We remark thatthe TE patterns in the MBL phase and near PDpoints are qualitatively different: in [44], we showthat TE of unfolded IM remains low in the formercase, and is high in the latter case. MBL and temporal long-range order —
TheIM contains information about the full disorder-averaged quantum noise spectrum of a system, allow-ing for a computation of time-dependent correlationfunctions. Here we consider the infinite temperaturecorrelator of the local magnetization, hh ˆ σ zp ( t )ˆ σ zp (0) ii ,whose long-time behavior provides a direct probeof ergodicity breakdown in the MBL phase, widelyused in experiments [14]; in this case, this correlatordoes not decay to zero as t → ∞ , indicating reman-ent magnetization. From a theory standpoint, thetime-averaged remanent magnetization is of key im-portance, since it reflects the emergence of LIOMs,providing a dynamical order parameter of MBL [46].Analysis of this quantity [47] is challenging due toinevitable presence of rare resonances, which haveto be treated non-perturbatively [7]. In contrast,disorder-averaged IM contains the contribution ofall resonances that are effective up to time t . Us-ing formula (3) and the MPS representation of theIM obtained by DMRG, we calculate the disorder-averaged correlator hh ˆ σ z ( t )ˆ σ z (0) ii [Fig. 2-(a)]. Weobserve that at strong disorder the magnetizationsaturates to a finite value, signalling MBL, while inthe critical region it continues decaying at accessibletime scales.Next, we inquire into the difference of the IMbetween the MBL and the ergodic phase. The gen-eral structure of the IM of a thermalizing bath hasbeen studied for an ensemble of non-interacting har-monic oscillators [30, 31], and more recently forquantum spin chains [29]. In this case the IMstrongly suppresses “quantum” trajectories where σ τ = ¯ σ τ , behaving similarly to a source of clas-sical noise. This causes the system to damp localquantum-interference effects and dephase its spins,thus erasing local memory after a finite correlationtime. An MBL system, in contrast, does not actas an efficient bath upon itself, producing quantumnoise that does not fully erase local memory of initialstates. We find that here the quantum trajectoriesare only weakly suppressed, reflecting the key roleof persistent interference processes.The onset of remanent magnetization can belinked to the appearance of temporal long-rangeorder in the IM. To that end, we write hh Z ii ≡ lim t →∞ t P tτ =1 hh ˆ σ zj ( τ )ˆ σ zj (0) ii as follows: hh Z ii = lim t →∞ t t X τ,τ =0 X { σ s } , { ¯ σ s } σ τ σ τ × δ (cid:16) X s ( σ s − ¯ σ s ) (cid:17) W σ s +1 ,σ s W ∗ ¯ σ s +1 , ¯ σ s h I i [ { σ, ¯ σ } ] , (6) t . . . . . hh ˆ σ z ( ) ˆ σ z ( t ) ii − . − . . . . m p ( m ) − . − . . . . m . . . . p ( m ) − . − . . . . m p ( m ) Figure 2. (a) Remanent magnetization hh ˆ σ z (0)ˆ σ z ( t ) ii vs time for different disorder strength equally spaced along theline J = g = 0 . . . . . (top to bottom curves). (b-d) Probability density p ( m ) of time-averaged magnetizationsectors in the ensemble of local spin trajectories defined by the IM: in the MBL phase ( J = g = 0 . ) (b), in thetransition region ( J = g = 0 . )(c), and in the ergodic phase ( J = g = 0 . ) (d). [cf. Eq. (3)] where we took into account that left andright IM are equal to h I i . Next, we represent the r.-h.s. of Eq. (6) as a sum over sets of trajectories withfixed magnetization M = P s σ s = P s ¯ σ s [48]. Thatenables us to take the sum over τ, τ , which, withina given magnetization sector, yields M . Thus, hh Z ii = lim t →∞ t X M = − t ( M/t ) P ( M ) = Z − dm m p ( m ) , (7)where P ( M ) is a positive “weight” representing thesum over trajectories in the sector with magnetiza-tion M . Unitarity dictates that P M P ( M ) = 1 , andtherefore P ( M ) may be viewed as a probability. Inthe last equation, we have switched to the magnet-ization density m = M/t in time and rescaled theprobability density p ( m ) = t P ( M ) .Expression (7) gives a necessary and sufficient cri-terion for MBL: ergodicity is broken, and local in-tegrals of motion exist, if the probability distribu-tion p ( m ) for the time-averaged magnetization ofindividual-spin trajectories has a finite width in theinfinite-time limit. In the ergodic phase, the widthof p ( m ) shrinks around its average m = 0 as /t ,satisfying central limit theorem scaling. We con-firmed this by an exact computation at the self-dual points, which gives a binomial distribution P ( M ) = 2 − t (cid:0) ( t + M ) / t (cid:1) . However, as the MBL phaseis approached, p ( m ) develops two symmetric peaksat finite values ± m ∗ . This highlights that the MBLphase dynamically breaks the Z -symmetry of thedisorder-averaged chain: Selecting “up” or “down”boundary conditions for the trajectory (i.e., initialand final state of spin p ) produces a finite bias inthe average magnetization towards the positive ornegative side, respectively.In Fig. 2(b-d) we report the results for the distri-bution p ( m ) obtained with MPS for the MBL phase,and by ED at the transition and in the ergodic phase.It is apparent that p ( m ) follows the above expect- ations, developing increasingly sharp peaks close to m = ± in the MBL phase as the observation timewindow t is enlarged [panel (b)]; in contrast, in theergodic phase [panel (d)], the distribution is single-peaked at m = 0 and narrows as t is increased.In the critical region between them [panel (c)], apeak at m = 0 is observed which slowly developsupon increasing t , reflecting the tendency to restor-ing thermal behavior at long times. Summary and outlook —
We have characterizeda many-body, disordered system via its influencematrix, which fully describes its properties as aquantum bath. This approach allows for exact dis-order averaging, therefore incorporating effects ofrare regions on MBL. The slow increase of temporalentanglement in the MBL phase reflects that the sys-tem fails to act as a thermal bath on itself, and opensthe door to efficient tensor-network approaches.We have implemented a proof-of-principle MPS al-gorithm and used it to extract time-dependent cor-relation functions that provide a benchmark for cur-rent quantum simulation experiments.Looking forward, the slow entanglement scaling ofIM in the MBL phase paves the way to constructingvariational, and possibly exact, solutions for the self-consistency equation (5) in the limit t → ∞ . Such asolution, and its breakdown at weaker disorder dueto proliferating resonances, will shed new light onMBL and the MBL-thermal transition. Finally, fur-ther applications of the IM approach may includeother form of ergodicity breaking such as quantumscars, time crystals as well as circuits that combineunitary evolution with measurements or dissipation. Acknowledgments —
We thank S. Choi, S. Gar-ratt and L. Piroli for insightful discussions. A.L.acknowledges valuable discussions with G. Giudiciand F. M. Surace on tensor-network implementa-tions. The MPS computations in this work wereperformed using TeNPy [49]. This work is suppor-ted by the Swiss National Science Foundation. ∗ These two authors contributed equally to this work[1] R. Nandkishore and D. A. Huse, Annual Review ofCondensed Matter Physics , 15 (2015).[2] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn,Rev. Mod. Phys. , 021001 (2019).[3] E. Altman, Nature Physics , 979 (2018).[4] M. Serbyn, Z. Papić, and D. A. Abanin, Phys. Rev.Lett. , 127201 (2013).[5] D. A. Huse, R. Nandkishore, and V. Oganesyan,Phys. Rev. B , 174202 (2014).[6] V. Ros, M. Müller, and A. Scardicchio, NuclearPhysics B , 420 (2015).[7] J. Z. Imbrie, Journal of Statistical Physics , 998(2016).[8] P. Ponte, Z. Papić, F. Huveneers, and D. A. Abanin,Phys. Rev. Lett. , 140401 (2015).[9] A. Lazarides, A. Das, and R. Moessner, Phys. Rev.Lett. , 030402 (2015).[10] P. Bordia, H. Lüschen, U. Schneider, M. Knap, andI. Bloch, Nature Physics , 460 (2017).[11] R. Moessner and S. L. Sondhi, Nature Physics ,424 (2017).[12] B. Bauer and C. Nayak, Journal of StatisticalMechanics: Theory and Experiment , P09005(2013).[13] J. Eisert, M. Cramer, and M. B. Plenio, Rev. Mod.Phys. , 277 (2010).[14] M. Schreiber, S. S. Hodgman, P. Bordia, H. P.Lüschen, M. H. Fischer, R. Vosk, E. Altman,U. Schneider, and I. Bloch, Science , 842 (2015).[15] 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).[16] A. Lukin, M. Rispoli, R. Schittko, M. E. Tai, A. M.Kaufman, S. Choi, V. Khemani, J. Léonard, andM. Greiner, Science , 256 (2019).[17] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W.Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Mon-roe, Nat. Phys. , 907 (2016).[18] B. Chiaro, C. Neill, A. Bohrdt, M. Filippone,F. Arute, K. Arya, R. Babbush, D. Bacon,J. Bardin, R. Barends, S. Boixo, D. Buell, B. Bur-kett, Y. Chen, Z. Chen, R. Collins, A. Dun-sworth, E. Farhi, A. Fowler, B. Foxen, C. Gidney,M. Giustina, M. Harrigan, T. Huang, S. Isakov,E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi, J. Kelly,P. Klimov, A. Korotkov, F. Kostritsa, D. Landhuis,E. Lucero, J. McClean, X. Mi, A. Megrant,M. Mohseni, J. Mutus, M. McEwen, O. Naaman,M. Neeley, M. Niu, A. Petukhov, C. Quintana,N. Rubin, D. Sank, K. Satzinger, A. Vainsencher,T. White, Z. Yao, P. Yeh, A. Zalcman, V. Smely-anskiy, H. Neven, S. Gopalakrishnan, D. Abanin,M. Knap, J. Martinis, and P. Roushan, “Directmeasurement of non-local interactions in the many-body localized phase,” (2020), arXiv:1910.06024[cond-mat.dis-nn].[19] W. De Roeck and F. m. c. Huveneers, Phys. Rev. B , 155129 (2017). [20] A. Pal and D. A. Huse, Phys. Rev. B , 174411(2010).[21] D. J. Luitz, N. Laflorencie, and F. Alet, Phys. Rev.B , 081103 (2015).[22] X. Yu, D. J. Luitz, and B. K. Clark, Phys. Rev. B , 184202 (2016).[23] M. Serbyn, Z. Papić, and D. A. Abanin, Phys. Rev.B , 104201 (2017).[24] A. Chandran, J. Carrasquilla, I. H. Kim, D. A.Abanin, and G. Vidal, Phys. Rev. B , 024201(2015).[25] M. Serbyn, A. A. Michailidis, D. A. Abanin, andZ. Papić, Phys. Rev. Lett. , 160601 (2016).[26] F. Pollmann, V. Khemani, J. I. Cirac, and S. L.Sondhi, Phys. Rev. B , 041116 (2016).[27] T. B. Wahl, A. Pal, and S. H. Simon, Phys. Rev.X , 021018 (2017).[28] M. C. Bañuls, N. Y. Yao, S. Choi, M. D. Lukin, andJ. I. Cirac, Phys. Rev. B , 174201 (2017).[29] A. Lerose, M. Sonner, and D. A. Abanin, “Influencematrix approach to many-body floquet dynamics,”(2020), arXiv:2009.10105 [cond-mat.str-el].[30] R. Feynman and F. Vernon, Annals of Physics ,118 (1963).[31] A. J. Leggett, S. Chakravarty, A. T. Dorsey,M. P. A. Fisher, A. Garg, and W. Zwerger, Rev.Mod. Phys. , 1 (1987).[32] M. C. Bañuls, M. B. Hastings, F. Verstraete, andJ. I. Cirac, Phys. Rev. Lett. , 240603 (2009).[33] We note that for problems with discrete disorderpotential, an alternative tensor-network approach todisorder averaging exists [50].[34] L. Zhang, V. Khemani, and D. A. Huse, Phys. Rev.B , 224202 (2016).[35] M. Sonner, M. Serbyn, Z. Papic, and D. A. Abanin,“Thouless energy in a floquet model of many-bodylocalization,” (2020), to appear.[36] D. A. Abanin, W. D. Roeck, and F. Huveneers,Annals of Physics , 1 (2016).[37] C. Sünderhauf, D. Pérez-García, D. A. Huse,N. Schuch, and J. I. Cirac, Phys. Rev. B , 134204(2018).[38] B. Bertini, P. Kos, and T. c. v. Prosen, Phys. Rev.Lett. , 264101 (2018).[39] A. Chan, A. De Luca, and J. T. Chalker, Phys.Rev. Lett. , 060601 (2018).[40] By virtue of the Z symmetry of the model.[41] Mathematically, this means that the null subspaceof ˆ T fractures into nilpotent Jordan blocks of size ≤ t .[42] L. Piroli, B. Bertini, J. I. Cirac, and T. c. v. Prosen,Phys. Rev. B , 094304 (2020).[43] A. Müller-Hermes, J. I. Cirac, and M. C. Banuls,New Journal of Physics , 075003 (2012).[44] “Supplemental online material,”.[45] U. Schollwöck, Annals of physics , 96 (2011).[46] A. Chandran, I. H. Kim, G. Vidal, and D. A.Abanin, Phys. Rev. B , 085425 (2015).[47] V. Ros and M. Müller, Phys. Rev. Lett. , 237202(2017).[48] We remind that the delta-function there originates from disorder-averaging, which effectively cancelsout interference terms between sectors with differentmagnetization.[49] J. Hauschild and F. Pollmann, SciPost Phys. Lect. Notes , 5 (2018), code available from https://github.com/tenpy/tenpy , arXiv:1805.00055.[50] B. Paredes, F. Verstraete, and J. I. Cirac, Phys.Rev. Lett. , 140501 (2005). upplemental MaterialCharacterizing many-body localization via exact disorder-averaged quantum noise Michael Sonner, ∗ Alessio Lerose, ∗ and Dmitry A. Abanin Department of Theoretical Physics, University of Geneva,Quai Ernest-Ansermet 24, 1205 Geneva, Switzerland (Dated: December 18, 2020)
In this Supplemental Material, we provide additional details on the numerical calculations used in the main text.
EXACT DIAGONALIZATION
As a benchmark for MPS calculations and to explore the ergodic side of the MBL transition, we employ an exactdiagonalization (ED) code. To this end, we interpret the influence matrix (IM) I [ { σ, ¯ σ } ] as a “wavefunction” in the2 t -dimensional vector space spanned by the basis {| σ = ± , . . . , σ t − = ± , ¯ σ t − = ± , . . . , ¯ σ = ±i} . Accordingly, weexpress the transfer matrix [ T ] { σ, ¯ σ } , { s, ¯ s } in Eq. (5) of the main text [obtained by taking the average of Eq. (4) of themain text over the random phase h j ] as a product of three operators,ˆ T = ˆ V ˆ P ˆ W , (S1a)[ V ] { σ, ¯ σ } , { s, ¯ s } = t − Y τ =0 exp ( iJσ τ s τ ) t − Y τ =0 exp ( − iJ ¯ σ τ ¯ s τ ) , (S1b)[ P ] { σ, ¯ σ } , { s, ¯ s } = t − Y τ =0 δ σ τ s τ t − Y τ =0 δ ¯ σ τ ¯ s τ (cid:16) δ P t − τ =0 s τ , P t − τ =0 ¯ s τ (cid:17) , (S1c)[ W ] { σ, ¯ σ } , { s, ¯ s } = t − Y τ =0 δ σ τ s τ t − Y τ =0 δ ¯ σ τ ¯ s τ δ s t − ¯ s t − t − Y τ =0 W s τ +1 s τ W ∗ ¯ s τ +1 ¯ s τ ρ s ¯ s ! , (S1d)where W σ σ = h σ | e ig ˆ σ x | σ i ≡ cos g δ σ ,σ + i sin g δ σ , − σ . This representation is illustrated in Fig. S1-(a,b). We note thatthe last kick has been erased compared to Fig. 1-(c) of the main text, exploiting its unitarity; a further simplificationcan be done to the first longitudinal field when ρ = / s = ¯ s in ˆ P .The combination of these two operations decreases the total dimensionality of the input and output vector space bya factor 4.The matrices ˆ W and ˆ P are diagonal in the computational basis, which we can interpret as the ˆ σ z product basisfor our chain of 2 t spins-1 /
2. Within this interpretation, the operator ˆ V is diagonal in the ˆ σ x product basis. Weexploit this fact for an efficient implementation of the applications of ˆ T via the Fast Walsh Hadamard transformation(FWHT) S1 , which can be used to convert between these two bases with O ( t t − ) basic operations. Note that theFWHT is an involution. The algorithm to apply ˆ T to an arbitrary vector thus reads:1. multiply its components by the diagonal components of ˆ P and ˆ W ;2. apply FWHT;3. multiply the vector’s components by the (now diagonal) components of ˆ V ;4. apply FWTH again.To find the thermodynamic-limit IM we simply use this procedure t times starting from the boundary vector I [ { s, ¯ s } ] ≡
1. For the infinite-temperature initial ensemble, t/ T , thus allowing us to push our ED results up to t = 12 with modest resources. MATRIX PRODUCT STATES
For the matrix product state (MPS) method, we choose to work in the folded picture where each spin on theforward trajectory is paired up with its corresponding spin on the backward trajectory such that they form a 4-dimensional local space S2 . This gives an open chain geometry, thus avoiding complications arising from periodic a r X i v : . [ c ond - m a t . d i s - nn ] D ec a) c) ≡ δ ∑ τ s τ ,∑ τ ¯ s τ ∏ τ δ s τ , σ τ δ ¯ s τ ,¯ σ τ ≡≡ α β δ γ e iJσ α σ β cos g δ σ γ , σ δ + i sin g δ σ γ ,− σ δ b) s ¯ s s t −1 ¯ s t −1 σ ¯ σ σ t −1 ¯ σ t −1 } ̂𝒫 } ̂𝒱 } ̂𝒲 ⇝ ̂𝒫 MPO ≡ δ sσ δ ¯ s ¯ σ δ B + s −¯ s , B ′ s ¯ s σ ¯ σ BB ′ FIG. S1. Panels a), b): Graphical representation of Eqs. (S1). Panel c): Graphical illustration of the MPO representation ofˆ P in Eq. (S4) boundary conditions of the closed Keldysh contour, present in the unfolded representation. We label and order thebasis vectors of this local space by S = ( s, ¯ s ) = ( ↑ , ↑ ) , ( ↓ , ↓ ) , ( ↑ , ↓ ) , ( ↓ , ↑ ). It is straightforward to write the factors ˆ V and ˆ W in the transfer matrix decomposition in Eq. (S1) as matrix product operators (MPOs) of bond dimensions 1and 4, respectively:[ V ] { S τ } , { Σ τ } = t − Y τ =0 e − iJ e iJ e iJ e − iJ e − iJ e iJ e iJ e − iJ S τ , Σ τ (S2)[ W ] { S τ } , { Σ τ } = X { A τ } A t t − Y τ =0 cos g sin g i sin g cos g − i sin g cos g sin g cos g − i sin g cos g i sin g cos gi sin g cos g − i sin g cos g cos g sin g − i sin g cos g i sin g cos g sin g cos g A τ +1 ,A τ δ A τ ,S τ δ S τ , Σ τ ρ A (S3)In the second equation, the labels { A τ = ( a τ , ¯ a τ ) } of virtual bonds run over the bases of four-dimensional ancillaryspaces isomorphic to the local folded spin spaces. The projection operator ˆ P which arises from disorder averagingcan be represented as a MPO with maximum bond dimension t . Here the virtual index B T on each bond ( T, T + 1)represents the difference in magnetization P Tτ =0 s τ − ¯ s τ between the portions of forward and backward paths to theleft of this bond. In other words, the local matrices P S τ , Σ τ B τ ,B τ +1 = δ S τ , Σ τ δ B τ + s τ − ¯ s τ ,B τ +1 (S4)composing the MPO, read the incoming virtual index B τ and add the local magnetization difference s τ − ¯ s τ to it toproduce the outgoing virtual index B τ +1 . The value of the virtual index B τ can thus remain unchanged, increase byone or decrease by one when proceeding to B τ +1 , making the local bond dimension χ τ = 2 τ + 1. Since ˆ P globallyprojects onto the zero magnetization sector P τ ( s τ − ¯ s τ ) = 0, we only need to carry virtual indices which are consistentwith this sector. Thus the maximal bond dimension is t at the central bond(s) of the chain. This MPO is representedin Fig. S1-(c).Due to the large bond dimension it is not possible to iteratively apply the transfer matrix MPO for ˆ T to an MPS atonce and compress the result. Instead the MPS needs to be compressed during the application of the MPO, using thezip-up method S3 . We can improve the iterative MPS by feeding it into a two-site DMRG code where we target theeigenvector with largest absolute eigenvalue. As this method is variational, it avoids compounding errors in contrastto the iterative method S4 .Different metrics of the quality of the MPS are displayed in Fig. S2. The role of energy in conventional DMRG istaken by the eigenvalue of the transfer matrix MPO in the last DMRG step. Due to the pseudo-projection propertyof ˆ T (see the main text), the exact eigenvalue is 1. We can also inspect the discrepancy of the IM entry on classical . . . . g T . . . . .
000 0 . . . . g T . . . . .
00 0 .
05 0 .
10 0 . g . . . . . . T . . . . . . FIG. S2. Different metrics to assess the quality of the MPS representation of the IM, as a function of the model parameter g = J and of the evolution time T . Left panel: DMRG eigenvalue of the transfer matrix; the exact value is 1 (see the maintext). Middle panel: average IM elements of classical trajectories { σ = ¯ σ } ; the exact value of each of them is 1 (see Ref. S2).Right panel: expectation value of the identity matrix hh ( t ) (0) ii [cf. Fig. 1-(a) of the main text] computed using the DMRGIM. For all the three metrics, depletion from the exact result is consistently observed for large T and large J = g .
10 20 30 t . . . . . . hh σ z ( ) σ z ( t ) ii MPS, χ = 128MPS, χ = 64ED FIG. S3. Comparison between the evolution of the dynamical structure factor hh ˆ σ z ( t )ˆ σ z (0) ii computed with ED and MPSmethods, for J = g = 0 . trajectories { σ = ¯ σ } from the exact value 1 (see Ref. S2). Last, we can compute the full path-integral of Fig. 1-(a) ofthe main text, with all observables of spin p set to identity, using the DMRG IMs. The exact value of this quantityis 1. We observe consistent deviations of the three metrics from their exact values for intermediate coupling strength(i.e., towards the MBL transition) and large times beyond the reach of ED, the last metric being the most sensitive.Furthermore, Fig. S3 shows the agreement between converged MPS data and ED data for the dynamical structurefactor hh ˆ σ z ( t )ˆ σ z (0) ii .Lastly, we elucidate the differences in entanglement patterns between the MBL phase and the vicinity of PDpoints. Since both regions show low temporal entanglement (albeit scaling differently with the evolution time), itis legitimate to question the nature of this apparent similarity. The latter, however, disappears upon unfolding theIM wavefunction, as illustrated in Fig. S4. Here, we have the possibility to consider a bipartition separating forwardand backward spins. The entanglement entropy S f/b ( t ) associated with such a bipartition is only low in the MBLphase, whereas it is large in the ergodic phase. At the self-dual point, in particular, we have S f/b ( t ) = ( t + 1) log 2,because the IM wavefunction is a product state of t + 1 maximally entangled Bell pairs between each spin on theforward branch and its equal-time partner on the backward branch. This occurrence is interpreted as follows. In theMBL phase, there exist as much correlations between spins on the same time branch as between spins on oppositetime branches. Conversely, ergodicity produces strong correlations between forward and backward trajectories: amanifestation of the suppression of quantum interference. ∗ These two authors contributed equally to this work[S1] T. L. Lezama, S. Bera, and J. H. Bardarson, Physical Review B , 161106 (2019). . . . . . g S foldedunfolded FIG. S4. Entanglement pattern across the MBL transition. Data are obtained from ED computations of the IMs for t = 12.The two curves show the bipartite entanglement entropy of the IM wavefunction corresponding to two structurally differentbipartitions. The “folded” curve shows S ( t/
2) of the folded IM, as in the main text; this quantity is low both near thefully decoupled line J = 0 and near the self-dual point J = g = π/
4. The “unfolded” curve S f/b ( t ) corresponds instead tobipartitioning the chain into forward and backward spins; this quantity is only low in the MBL phase, but is very large in theergodic phase.[S2] A. Lerose, M. Sonner, and D. A. Abanin, “Influence matrix approach to many-body floquet dynamics,” (2020),arXiv:2009.10105 [cond-mat.str-el].[S3] E. Stoudenmire and S. R. White, New Journal of Physics , 055026 (2010).[S4] U. Schollw¨ock, Annals of physics326