Almost-conserved operators in nearly many-body localized systems
Nicola Pancotti, Michael Knap, David A. Huse, J. Ignacio Cirac, Mari Carmen Bañuls
AAlmost-conserved operators in nearly many-body localized systems
Nicola Pancotti, Michael Knap, David A. Huse, J. Ignacio Cirac, and Mari Carmen Ba˜nuls Max-Planck-Institut f¨ur Quantenoptik, Hans-Kopfermann-Str. 1, 85748 Garching, Germany Department of Physics and Institute for Advanced Study,Technical University of Munich, 85748 Garching, Germany Physics Department, Princeton University, Princeton, New Jersey 08544, USA (Dated: April 12, 2018)We construct almost conserved local operators, that possess a minimal commutator with theHamiltonian of the system, near the many-body localization transition of a one-dimensional disor-dered spin chain. We collect statistics of these slow operators for different support sizes and disorderstrengths, both using exact diagonalization and tensor networks. Our results show that the scalingof the average of the smallest commutators with the support size is sensitive to Griffiths effects inthe thermal phase and the onset of many-body localization. Furthermore, we demonstrate that theprobability distributions of the commutators can be analyzed using extreme value theory and thattheir tails reveal the difference between diffusive and sub-diffusive dynamics in the thermal phase.
I. INTRODUCTION
Isolated interacting quantum systems can undergo anunconventional transition from an ergodic thermal phasein which the system approaches thermal equilibrium irre-spective of the initial conditions to a many-body localized(MBL) phase in which it does not.
There exists firmevidence, from numerical, rigorous mathematical, andexperimental work, that the MBL phase is real-ized in strongly disordered and interacting quantum sys-tems. Recently, the MBL transition between the ther-mal and the MBL phase has been started to be investi-gated using general arguments, mean-field theory, and renormalization-group schemes . Since the MBLtransition is inherently dynamical and can occur at arbi-trary energy density, it is beyond the scope of a conven-tional thermodynamic description.A possible strategy to explore MBL is to directly studyproperties of the Hamiltonian. In particular, a phe-nomenological description of the Hamiltonian in the MBLphase using an extensive set of local integrals of motion(LIOMs) has been proposed. This approach explains,for instance, the dephasing and entanglement dy-namics in the localized phase. Moreover, LIOMs havebeen constructed explicitly for some models using ana-lytical techniques, exact diagonalization , stochas-tic methods, and tensor networks. Upon reducingthe disorder, LIOMs become more and more extended inspace and eventually cease to exist at the MBL transi-tion.Here, we take another route and study the dynam-ics of a disordered and interacting Heisenberg spin chainby constructing almost conserved local operators with fi-nite support, which minimize the commutator with thesystem Hamiltonian. Deep in the MBL phase the slowoperators resemble LIOMs. However, our procedure canbe directly extended to the thermal phase as well. Theslow operators can potentially be used to study dynami-cal properties, because the value of the commutator givesa lower bound on the thermalization times of the corre- sponding operator. We compute the slowest operatorsfor different support sizes and disorder strengths, usingboth exact diagonalization and tensor networks.
Ourwork, rather than focusing on the phase transition, itaims at characterize the sub-diffusive region in its vicin-ity. Specifically, we collect statistics from numerous dis-order realizations and show that not only the mean valueof the smallest commutator, but also the probability dis-tributions are sensitive to the underlying phase. The tails of the distributions are affected by the appearance of rareGriffiths regions, which act as bottlenecks for trans-port and are argued to provide a simple description forthe slow sub-diffusive dynamics in the thermal phase nearthe MBL transition.
We use extreme value theory(EVT) to analyze the tails of the distribution and findthat it is well described by the generalized extreme valuedistribution. Moreover, we discuss how the presence ofrare regions may affect the asymptotic behavior of suchdistributions. Our results demonstrate that finding slowoperators provides access dynamical quantities withoutresorting to the time evolution of a particular state.The paper is organized as follows. In section II wespecify our model. Section III describes the slow opera-tor method, and presents the fundamental ideas of EVTemployed in our analysis. The results for the structureof the slow operators are presented in section IV, the av-erage values of their commutators in the different phasesare shown in section V, and the statistical analysis of theprobability distributions using EVT is presented in sec-tion VI. Finally, in section VII we summarize our findingsand discuss potential extensions of our work.
II. MODEL
We consider the non-integrable spin-1/2 HeisenbergHamiltonian with random magnetic field, H = (cid:88) i J (cid:0) S xi S xi +1 + S yi S yi +1 + S zi S zi +1 (cid:1) + h i S zi , (1) a r X i v : . [ c ond - m a t . d i s - nn ] A p r where S αi are spin-1/2 operators, and the values of thetransverse field h i are randomly chosen from the uniformdistribution in the interval [ − h, h ], where h is the disorderstrength. In the following, we set J = 1. This modelexhibits a phase transition from a thermal to a MBLphase at h c ∼ In the MBL phase the system can be described by LI-OMs, which lead to an emergent integrability. The corre-sponding effective Hamiltonian can be expressed as H = (cid:88) i τ zi + (cid:88) ij J ij τ zi τ zj + (cid:88) n (cid:88) i,j, { k } K ( n ) i { k } j τ zi τ zk ...τ zk n τ zj . (2) The LIOMs, τ zi , commute with each other and withthe Hamiltonian. They are exponentially localized whenwritten in terms of the physical spins { (cid:126)S i } . Moreover,the coefficients J ij and K ( n ) i { k } j decay exponentially withtheir spatial range. III. METHODSA. Slow operator technique
Almost conserved quasi-local operators can identifylong thermalization time scales. A way to compute suchtime scales is to search for slow operators that minimizethe commutator with the Hamilonian (see also Refs.52–54). In particular, we may restrict the search to op-erators O M with finite support on M consecutive sites.Defining L ( O M ) to be the squared norm of the commuta-tor between the Hamiltonian and the (normalized) oper-ator, the problem reduces to solving the variational min-imization λ M := min O M ; tr O M =0 L ( O M ) = min O M ; tr O M =0 (cid:107) [ H, O M ] (cid:107) F (cid:107) O M (cid:107) F , (3)where (cid:107) A (cid:107) F = tr( A † A ) is the Frobenius norm, and O M isrestricted to be traceless. The operator that minimizes(3) is the one evolving the slowest in Heisenberg pictureat time t = 0. Moreover, for a slightly perturbed infinitetemperature state of the form ρ ( O M ) = /Z + (cid:15)O M thethermalization time is lower bounded as t th ≥ / √ λ M . If the system has a conserved extensive quantity, suchas the total energy, it naturally gives rise to slow opera-tors. For a translationally invariant system, the Hamil-tonian commutes with the sum of all translations of thecorresponding local operator, but not with translationsrestricted to a finite window M . Nevertheless, the normof the latter commutator seems to decrease polynomiallywith M . In our particular case of Hamiltonian (1), the totalpolarization in the z direction, (cid:80) i σ zi , is conserved forany value of the disorder strength. The restriction ofthe sum to M sites, Z M = (cid:80) Mi =0 sin (cid:0) πkM (cid:1) σ zi , is thus anatural slow operator for any given support M . It is easyto check that as the support increases, the corresponding commutator decays as L ( Z M ) ∼ /M . This sets areference for comparison with the commutators found bynumerically solving the variational problem, Eq. (3). For small support size M the optimization can besolved via exact diagonalization. This has been used toshow that for model (1) there exists an extensive numberof exponentially localized LIOM in the MBL phase whichare not present in the thermal phase . To reach largersupports the search can be restricted to operators O M of the form of a matrix product operator (MPO). Thenstandard tensor network techniques can be applied tosolve this optimization, and the bond dimension can thenbe systematically increased until convergence is achieved.In this work, we use exact and approximate MPO so-lutions in order to collect information on the minimalcommutators and the corresponding operators. For dif-ferent values of the support M and the disorder strength h , we compute λ M for several configurations of the ran-dom magnetic field (ranging from 2 · for M = 12 to 10 for M = 4). In order to simplify the statistical analysis,we choose independent configurations over M + 2 sites.This is equivalent to choosing non-overlapping windowsover an infinite chain.To solve the problem numerically, the minimizationEq. (3) can be interpreted as an eigenvalue problem. Forfixed support size M and disorder strength h and for aparticular disorder realization, there is an effective op-erator H ( M )eff acting on the d M -dimensional space withsupport on the chosen window, such that the (vectorized)operator that minimizes Eq. (3) corresponds to its eigen-vector with lowest eigenvalue. Since the optimization isrestricted to operators supported on a certain windowof size M , only terms in the Hamiltonian that overlapwith that region will contribute to the commutator. Forour nearest-neighbor Hamiltonian this means that M + 2sites of the Hamiltonian contribute, which we denote by H M +2 . The matrix representation of the correspondingcommutator, acting on vectorized operators, can be writ-ten as C M +2 = H M +2 ⊗ − ⊗ H TM +2 , and the effectiveoperator can be constructed as (cid:104) O M |H ( M )eff | O M (cid:105) := tr (cid:16) [ H M +2 , ˜ O M ] † [ H M +2 , ˜ O M ] (cid:17) = (cid:104) ˜ O M |C M +2 C † M +2 | ˜ O M (cid:105) , (4)where ˜ O M = ⊗ O M ⊗ and the trace is taken over thespace of all M + 2 sites. B. Extreme value theory
The statistical analysis of the smallest commutatorcontains relevant information about the phases of themodel and the critical region, both in the average scalingof λ M with the support size M and in the probability dis-tributions. In particular, the probability density function(PDF) p ( λ M ) and the corresponding cumulative densityfunction (CDF), F ( λ M ) = (cid:82) λ M −∞ p ( x ) dx , will be sensitiveto the presence of rare regions of atypically large disorderwhich emerge in the proximity of the critical point.As we show in this work, these PDFs can be describedand analyzed within the mathematical framework of Ex-treme Value Theory (EVT), a branch of statistics con-cerned with the description of rare events, which is ap-plied to floods, earthquakes or risk in finance and insur-ance, as well as to several branches of physics, includingstatistical mechanics and critical phenomena. In particular, EVT deals with the asymptotic behaviorof the extreme values of a sample, i.e., with the tails ofthe distributions. Let us consider a random variable x governed by a CDF F ( x ). Then EVT ensures that themaxima (equivalently the minima) of samples of F ( x ),properly normalized and centered, will be governed by aCDF G ζ ( y ) = exp (cid:16) − (1 + ζy ) − /ζ (cid:17) , (1 + ζy ) − /ζ > y = ax + b with a, b ∈ R and a >
0. The Fisher-Tippett-Gnedenko theorem (Thm. 1.1.3 in 50) statesthat this single-parameter family, called generalized ex-treme value distribution (GEV), includes all possible lim-iting distributions for the extreme values of a sample ofi.i.d. random variables. The family contains three sub-classes which exhibit quite different behavior, the Gum-bel ( ζ = 0), Fr´echet ( ζ >
0) and Weibull ( ζ <
0) distri-butions. Qualitatively, the deciding criterion for a PDFto belong to the basin of attraction of one family or an-other (i.e. for the extrema of the samples to be describedby the corresponding family) is the form of its tails. In our particular problem, as we try to find the min-imum λ M for a certain configuration, we are effectivelysampling from the left tail of p (Λ M ), which is the distri-bution of eigenvalues Λ M of H ( M )eff . Typically the eigen-values of a matrix are not uncorrelated, and thus GEVis a priori not expected to describe the probability dis-tribution of extreme eigenvalues . Nevertheless, ourresults indicate that the distribution (5) provides indeeda very good description for our data, with the particu-lar form depending only on the asymptotic behavior of p (Λ M ) for small Λ M . IV. STRUCTURE OF THE SLOW OPERATORS
We will first study the structure of the slowest opera-tors. In the strong disorder regime, we expect that theslow operators correspond to some LIOMs, or more pre-cisely to a truncated version of them since our supportsize is fixed. It is worth noticing here that although theslow operator method does not directly target the (trun-cated) LIOMs, as was done in Refs. 39 and 61, in thelocalized regime, truncated LIOMs and their combina-tions are good candidates to attain (exponentially) smallcommutators.To address this question, we can analyze to whichextent the slowest operators are local by examining i − − − − − − − | c ( i ) | i − − − − | c ( i ) | h = 0 . h = 2 . h = 3 . h = 5 . h = 10 . FIG. 1.
Local structure of the slowest operator withsupport M = 40 . The weight of single-site | c ( i ) | contribu-tion is shown as a function of the position i for various disorderstrengths for a single disorder realization. At weak disorder,the local contributions are spread over the whole support andhas the shape of a sine as expected from diffusive dynam-ics. For strong disorder, the operator becomes increasinglylocalized and exhibits exponential tails. The optimization hasbeen performed for MPO bond dimension D = 100. Inset:
For larger support contributions such as | c ( i ) | , we observesimilar tails with considerably smaller values. their spatial structure in different regimes. Thiscan be done by studying their decomposition as asum of tensor products of single-site Pauli matri-ces, O M = 2 − M/ (cid:80) { α j } =0 C α ...α M σ α . . . σ α M M . Inprinciple, for an operator O M found by the mini-mization, either exactly or as a MPO approximation,we can efficiently evaluate any single coefficient C α ...α M = 2 − M/ tr( σ α . . . σ α M M O M ). Nevertheless, dueto the exponential growth of the basis dimension withthe support, already for moderate values of M it is un-feasible to inspect all the individual coefficients. Instead,a more physical quantity for exploring the localized na-ture of the operators is the combined weight of all termswith a fixed range which are supported on a certainsubregion. We define a range- k operator as a product ofPauli matrices that act non trivially on sites i to i + k .Formally, this corresponds to the operator Θ a,b, { α m } i,k = σ . . . σ i − σ ai σ α i +1 . . . σ α k − i + k − σ bi + k σ i + k +1 . . . σ M , where a, b ∈ , , i → i + k ) and { α m } ∈ , , , k terms can be written as | c k ( i ) | := 12 M (cid:88) a,b =1 3 (cid:88) { α m } =0 (cid:12)(cid:12)(cid:12) Tr (cid:16) Θ a,b, { α m } i,k O M (cid:17)(cid:12)(cid:12)(cid:12) . (6)We thus solve the optimization (3) and compute theweights, | c k ( i ) | , along the chain of terms with fixed range k . Even if the support is large, as O M is written as an − h . . . . . / α M − − − E ( λ M ) ( a ) h = 0 . h = 1 . h = 2 . h = 3 . h = 5 . h . . . . / ξ M − − − − E ( λ M ) ( b ) h = 3 . h = 4 . h = 5 . h = 7 . h = 9 . FIG. 2.
Disorder average of the smallest commutator E ( λ M ) . (a) Decay of E ( λ M ) at weak disorder. We observe that theaverage of λ M scales polynomially M − α with an exponent that increases as the transition is approached. Inset:
The inverseexponent 1 /α as a function of the disorder strength h is shown. (b) Decay of E ( λ M ) at strong disorder. For strong disorder,the decay of λ M is compatible with an exponential form, e − M/ξ . Inset:
The inverse length scale 1 /ξ decays as the transitionpoint is approached. MPO, this quantity can be computed efficiently, i.e., witha cost that only scales polynomially in the range k andthe support M . We find clear differences in the structure of the oper-ator depending on the disorder strength. This is illus-trated in Fig. 1 for the MPO, that minimizes Eq. (3) fora particular disorder realization h i = h · r i , where r i arefixed random numbers for each value of h . We have cho-sen a support of M = 40 sites and a relatively large bonddimension of D = 100, in order to ensure that truncationerrors are negligible compared to the effects we observe.The figure shows the different spatial profile of the single-site contributions | c ( i ) | = 2 − M (cid:80) a =1 | Tr ( σ ai O M ) | and two-site contributions | c ( i ) | (Eq.(6) with k = 1)for different values of the disorder strength. For strongdisorder, where the LIOMs are exponentially localized,at least one of the | c ( i ) | is expected to be dominant.Indeed, we observe that for large disorder the operatoris well localized around a single site, with weights thatdecay exponentially around this point. As the disorderstrength decreases, the weights decay slower with dis-tance, and the tails are no longer exponential. Finally forvery small h the profile becomes flat and has the shape ofa sine as expected in the diffusive regime. The single-siteterms shown in the figure always accumulate most of theweight. We found that higher order terms Fig. 1 (in-set), which have smaller contributions, exhibit a similarspatial decay. The observed profile at strong disorder isin good agreement with what one could expect from thedecomposition of the LIOMs in the canonical operatorbasis, what would in principle allow the extraction of alocalization length scale of the LIOM. V. AVERAGE COMMUTATOR
The precise value of the minimum λ M for fixed support M and disorder strength h will depend on the disorderrealization. Hence, the average over realizations E ( λ M )will provide information about the underlying phase andalso indirectly about transport.In the following, we will optimize Eq. (3) over the fam-ily of MPOs with finite bond dimension D = 10. Thisallows us to reach large support sizes M and collect a sig-nificant amount of statistics. The average of the smallestcommutator E ( λ M ) is shown as a function of the supportsize M for different disorder strength h in Fig. 2. Atvery weak disorder h (cid:46) . E ( λ M ) ∼ M − of conserved quantities.When increasing the disorder strength, the decay re-mains a powerlaw E ( λ M ) ∼ M − α , however, with an in-creasing exponent α indicating that the transport be-comes slower due to the presence of disorder, see inset ofFig. 2 (a) which shows the inverse exponent 1 /α . Thisis consistent with the observation of sub-diffusive trans-port on the thermal side of the MBL transition. The longer relaxation times in this regime can be in-terpreted as follows: rare insulating regions with largerthan typical disorder act as bottlenecks for the conven-tional diffusive transport.
Upon approaching theMBL transition at h ∼ we expect the exponent todiverge. Our data indeed shows that the power law de-cay becomes very steep. Due to the finite support sizes M and due to the slow convergence of the algorithm forlarge disorder, we, however, find that the exponent levelsoff at a finite value. For practical purposes, we adopt thestrategy of limiting the maximum computational effortfor each search. Therefore, the result of the algorithmprovides an upper bound for the true minimum λ M ( h ), .
00 0 .
05 0 .
10 0 .
15 0 .
20 0 .
25 0 . λ M P D F M = 12 GEV fit h = 0 . h = 1 . h = 2 . h = 3 . FIG. 3.
Probability density function of λ M . Using themaximum likelihood method we fit the numerical data for λ M with the GEV distribution (5). The shape of the distributionfunction is strongly influenced by the presence of disorder. and the corresponding α is biased towards smaller values.When crossing the MBL transition, an extensive num-ber of exponentially localized conserved quantities, corre-sponding to the LIOMs τ zi in Eq. (2), emerges. There-fore, provided the window M is sufficiently large to cap-ture most of the support of any LIOM, the smallest com-mutator should decay exponentially as λ M ∼ e − M/ξ .Such an exponential decay is supported by our numer-ical simulations for strong disorder, see Fig. 2 (b). Thelength scale ξ is related (in a possibly non-trivial way)to the localization length of the LIOMs. In the inset, weshow that with increasing disorder h , the inverse lengthscale 1 /ξ increases monotonically. The decays of α − and ξ − as the transition is approached are both com-patible with previous estimates of the critical region, seee.g. Refs. 51 and 61. VI. EXTREME VALUE THEORYA. Numerical results
We now study the probability distributions of the slowcommutators λ M , which can be sensitive to rare events.In particular, we show that EVT provides useful math-ematical tools to characterize the tails of the probabil-ity density and cumulative distribution functions. Wefirst illustrate in Fig. 3, that the GEV distribution,obtained from differentiating Eq. (5), provides a verygood fit to the collected data over a range of disorderstrengths. Moreover, we find that the shape of the dis-tribution strongly depends on the disorder value and thatthe peak of the distribution shifts from a finite, relativelylarge value for weak disorder, toward zero as h increases.In between those limits, in a regime where one expectssub-diffusive dynamics, the distribution broadens consid-erably. In the following we will analyze the tails of this distri-bution quantitatively. To this end, we introduce a newvariable, with dimensions of time, as T M = (Λ M ) − / .The minimal eigenvalue λ M thus corresponds to the max-imum of T M , that we call τ M , and has the physical in-terpretation of bounding the thermalization time of theoperator ρ ( O M ) from below. This transformation of vari-ables moves the left tails of p ( λ M ) that arise due torare events, to right tails of p ( τ M ). The correspond-ing survival functions S ( τ M ) = 1 − F ( τ M ) are shownin Fig. 4 on logarithmic scale. We observe that for smalldisorder, the survival function approaches zero super-exponentially fast describing weak tails, due to the smallprobability of large τ M . As the disorder increases, thesurvival function decays with a power-law tail. The valueof the disorder at which the tails change is in agreementwith the crossover from diffusive to sub-diffusive dynam-ics.To make this observation more quantitative, we studythe variation of the shape parameter ζ from the fittedGEV, defined in Eq. (5), with the disorder strength.Our data shows that ζ shifts from values close to zerofor weak disorder, to clearly positive values as the MBLphase is approached. The ζ parameter determines thetype of distribution. In particular, the deviation of ζ fromzero describes the crossover from a peaked distributionto one with polynomial tails. In terms of the GEV fami-lies, this corresponds to a change from a Gumbel ( ζ = 0)to a Fr´echet ( ζ >
0) distribution. The observed qualita-tive change can be explained by an intermediate regimein which atypically slow operators appear for any givensupport M , leading to strong tails in the probability den-sity function of τ M . The value of the disorder strengthat which the shape parameter ζ starts being clearly pos-itive is again consistent with the expected sub-diffusiveregime of the thermal phase. B. Interpretation in terms of Griffiths effects
Although the eigenvalues of the effective operators inour study are not expected to be uncorrelated, the re-sults in the previous section show that GEV does indeeddescribe our findings accurately. Hence, in this sectionwe use EVT arguments to show how the observed dis-tributions can be qualitatively explained in terms of theexistence of rare Griffiths regions.We first consider the effect of rare localized regions inthe thermal phase but close enough to the transition. Inthis situation, given a support size M , the typical valuesof λ M will be polynomial in M . Nevertheless, if rare re-gions are present that (partly) support an exponentiallylocalized operator, they can give rise to exponentiallysmall values of λ M . More concretely, within a fixed win-dow M , such a rare configuration of size (cid:96) < M will oc-cur with an exponentially small probability p ( (cid:96) ) ∼ c (cid:96) , forsome c <
1. This patch can support a localized operator,with localization length ξ , with a correspondingly small τ M − − − S ( τ M ) M = 8( a ) h = 0 . h = 1 . h = 2 . h = 3 . h = 5 . h = 7 . h = 9 . τ M − − − S ( τ M ) M = 12( b ) h = 0 . h = 1 . h = 2 . h = 3 . h = 5 . FIG. 4.
Survival function of τ M = (cid:112) /λ M . We fit the survival function of τ M to the GEV, Eq. (5), which decayssuper-exponentially for weak disorder and polynomially for large disorder. The form of these tails is sensitive to rare Griffithseffects. commutator Λ ∼ e − (cid:96)/ξ (or equivalently T ∼ e (cid:96)/ ξ ).Strictly speaking, this will be detected as the slowest op-erator only if this commutator is smaller than the (poly-nomial) typical commutator for the complementary re-gion of size M − (cid:96) . Since we are interested in the prob-ability for very small commutators, and there is an ex-ponential separation between the scaling of both terms,we may assume that small enough commutators will al-ways come from such rare patches. Evidently, this canonly be true for commutators below a certain ( M and h dependent) threshold. With this assumption, the proba-bility of some very small value of Λ is determined by theprobability of finding a rare patch such that Λ ∼ e − (cid:96)/ξ ,and thus we can identify p ( (cid:96) ) | d(cid:96) | = p ( T ) | dT | . From thatwe find that, for the range of T that correspond to rareregions of length (cid:96) < M , the probability for the largest T values is polynomial, p ( T ) ∼ ξT − ξ | ln( c ) |− . (7)For PDFs with polynomial tails, EVT predicts that theextreme values will be governed by a Fr´echet distributionwith ζ > ζ = 0(see appendix A). This is in agreement with the shapeparameter ζ shown in Fig. 5 where, for weak disorder,we obtain values of ζ very close to zero which correspondto a Gumbel distribution.Beyond the localization transition we expect the typ-ical regions to (partly) support localized operators, and h − . . . . . ζ M = 6 M = 8 M = 10 M = 12 FIG. 5.
Shape parameter ζ of the generalized ex-treme value distribution as a function of the disorderstrength. The values of the shape parameter ζ range from ∼ ζ ∼ / h isincreased. thus give rise to exponentially small values of λ M (cor-respondingly exponentially large values of τ M ), as weexplicitly observed in the average values shown in Sec-tion V. In our data for strong disorder, we have also foundbroad tails of the PDF p ( T ) which are consistent with apower law. Thus, the function describing our data resem-bles the Fr´echet distribution. This might be explainedby the fact that matrix elements in the MBL phase havebeen found to exhibit broad tails. As we approach the transition from the MBL side bylowering the disorder, rare thermal inclusions can appearthat potentially correspond to larger than typical com-mutators. Yet, our method only looks for the smallestcommutator in each given window. Thus, if a support M encloses one such thermal subregion, a competition arisesbetween the values of the commutator for the inclusionand that of the (typically) exponentially localized com-plement. Because we expect that a thermal inclusiongives rise to only polynomially decaying commutators,their value can only be the smallest one when the inclu-sion is sufficiently large in relation to the size M of thesupport. This causes our method to be less sensitive torare thermal regions in the localized side of the transi-tion.In contrast to other numerical studies, where thepresence of Griffiths effects was inferred from averagedobservables, with our method we may directly lo-cate rare regions. In particular, we can obtain the disor-der potential for the eigenvalues λ M that contribute tothe tail of the distribution and analyze the microscopicconfiguration of the random field { h i } in real space. Fol-lowing this procedure, we could, however, not unambigu-ously determine an obvious correlation between strongfluctuations of the field in real space and small commu-tators. It remains an open question whether it is possibleto predict the location of the rare regions from the dis-order landscape using a more direct method than theoptimization. VII. DISCUSSION
We have constructed slow operators with finite supportby minimizing their commutator with the Hamiltonian ofthe system using both exact diagonalization and tensornetwork techniques. In particular, we have consideredthe Heisenberg spin chain with random magnetic field,which displays a dynamical transition from the thermalto the many-body localized phase. The scaling of theminimal commutator with support size provides informa-tion on the localization transition as well as on transportin the system without resorting to a specific initial state.Furthermore, we have demonstrated that the tails ofthe probability distributions are sensitive to rare insu-lating regions in the thermal phase near the many-bodylocalization transition. We have shown that the statisticsof the smallest commutators can be analyzed within themathematical framework of extreme value theory. Inparticular, we have found that the distributions are welldescribed by generalized extreme value functions whoseshape depends on the disorder strength. By extremevalue theory arguments, the observed behavior in thetails can be connected to the appearance of rare, stronglydisordered regions, that give rise to atypically small mini-mal commutators. In particular, the disorder strength atwhich the distribution functions obtain power-law tailsis consistent with the appearance of sub-diffusive trans-port.
We conclude that the slow operator technique com-bined with extreme value theory, constitutes a valuabletool for exploring microscopic mechanisms of the MBLtransition. Further developments may include tailoringthe optimization technique to target explicitly the rare thermal regions on the localized side, which could pro-vide new insights into this less explored aspect of MBLphysics. Another intriguing question would be how cou-pling an MBL system to an external bath wouldchange the structure of slow operators.
ACKNOWLEDGMENTS
We acknowledge financial support from ExQM and IM-PRS (N.P.) and from the Technical University of Mu-nich - Institute for Advanced Study, funded by the Ger-man Excellence Initiative and the European Union FP7under grant agreement 291763 (M.K.). This work wasalso partially funded by the European Union through theERC grant QUENOCOBA, ERC-2016-ADG (Grant no.742102).
Appendix A: Extreme value theory
A given PDF is said to belong to the basin of attrac-tion of one of the extreme values distributions, namelyGumbel, Fr´echet or Weibull, when the extrema are dis-tributed according to the corresponding function. Thevon Mises conditions establish simple criteria to deter-mine whether p ( x ) belongs to one of them. In this ap-pendix we show how the conditions apply to the partic-ular cases of the distributions discussed in section VI B. a. Strong disorder implies a Fr´echet distribution. Rare regions in the thermal phase near the MBL tran-sition may support localized operators. As shown inEq. (7), the corresponding probability distribution is ex-pected to decay as p ( T ) ∼ ξT − ξ | ln c |− . A sufficientcondition for a PDF p ( T ) to belong to the basin ofattraction of the Fr´echet distribution is that the corre-sponding CDF F ( T ) = (cid:82) T −∞ p ( T (cid:48) ) dT (cid:48) satisfies the condi-tion lim T →∞ T F (cid:48) ( T )1 − F ( T ) = 1 ζ , (A1)with ζ > p ( T )1 − F ( T ) = p ( T ) (cid:82) ∞ T p ( t ) dt = 2 ξ | ln c | T , (A2)so that p ( T )1 − F ( T ) = 2 ξ | ln c | T , (A3)and, asymptotically,lim T →∞ T p ( T )1 − F ( T ) = 2 ξ | ln c | > , (A4)which ensures the condition above with ζ − = 2 ξ | ln c | ,and thus implies a limiting Fr´echet distribution. b. Exponentially decaying tails imply a Gumbel dis-tribution. In order to prove that a PDF belongs to thebasin of attraction of the Gumbel distribution it is suffi-cient to check the following condition lim T ↑ T max ddT (cid:18) − F ( T ) F (cid:48) ( T ) (cid:19) = 0 , (A5)We assume the simplest exponential decay for the righttail of the distribution p ( T ) ∼ e − kT . From the corre-sponding cumulative function we get the survival proba- bility, 1 − F ( T ) = (cid:90) ∞ T p ( t ) dt ∝ k e − kT , (A6)so that 1 − F ( T ) p ( T ) ∼ k . (A7)Thus, the derivative vanishes, which ensures Eq.(A5). D.M. Basko, I.L. Aleiner, and B.L. Altshuler, “Metal-insulator transition in a weakly interacting many-electronsystem with localized single-particle states,” Ann. Phys. , 1126 – 1205 (2006). I. Gornyi, A. Mirlin, and D. Polyakov, “Interacting elec-trons in disordered wires: Anderson localization and low- T transport,” Phys. Rev. Lett. , 206603 (2005). Rahul Nandkishore and David A. Huse, “Many-body lo-calization and thermalization in quantum statistical me-chanics,” Annu. Rev. Condens. Matter , 15–38 (2015). Ehud Altman and Ronen Vosk, “Universal dynamics andrenormalization in many body localized systems,” Annu.Rev. Condens. Matter , 383–409 (2015). Vadim Oganesyan and David A. Huse, “Localization ofinteracting fermions at high temperature,” Phys. Rev. B , 155111 (2007). Marko ˇZnidariˇc, Tomas Prosen, and Peter Prelovˇsek,“Many-body localization in the Heisenberg
XXZ magnetin a random field,” Phys. Rev. B , 064426 (2008). John Z. Imbrie, “On many-body localization for quantumspin chains,” J. Stat. Phys. , 998–1048 (2016). Michael Schreiber, Sean S. Hodgman, Pranjal Bordia, Hen-rik P. L¨uschen, Mark H. Fischer, Ronen Vosk, Ehud Alt-man, Ulrich Schneider, and Immanuel Bloch, “Observa-tion of many-body localization of interacting fermions in aquasirandom optical lattice,” Science , 842–845 (2015). S. S. Kondov, W. R. McGehee, W. Xu, and B. De-Marco, “Disorder-induced localization in a strongly corre-lated atomic hubbard gas,” Phys. Rev. Lett. , 083002(2015). J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W.Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Monroe,“Many-body localization in a quantum simulator with pro-grammable random disorder,” Nat. Phys. (AOP) (2016),10.1038/nphys3783. Pranjal Bordia, Henrik P. L¨uschen, Sean S. Hodgman,Michael Schreiber, Immanuel Bloch, and Ulrich Schneider,“Coupling identical one-dimensional many-body localizedsystems,” Phys. Rev. Lett. , 140401 (2016). Jae-yoon Choi, Sebastian Hild, Johannes Zeiher, PeterSchauß, Antonio Rubio-Abadal, Tarik Yefsah, Vedika Khe-mani, David A. Huse, Immanuel Bloch, and ChristianGross, “Exploring the many-body localization transitionin two dimensions,” Science , 1547–1552 (2016). Pranjal Bordia, Henrik L¨uschen, Ulrich Schneider, MichaelKnap, and Immanuel Bloch, “Periodically driving a many-body localized quantum system,” Nat Phys , 460–464. Henrik P. L¨uschen, Pranjal Bordia, Sebastian Scherg, Fa-bien Alet, Ehud Altman, Ulrich Schneider, and ImmanuelBloch, “Evidence for Griffiths-type dynamics near themany-body localization transition in quasi-periodic sys-tems,” (2016), arXiv:1612.07173. Pranjal Bordia, Henrik L¨uschen, Sebastian Scherg, SarangGopalakrishnan, Michael Knap, Ulrich Schneider, and Im-manuel Bloch, “Probing slow relaxation and many-bodylocalization in two-dimensional quasi-periodic systems,”arXiv:1704.03063 (2017). T. Grover, “Certain General Constraints on the Many-Body Localization Transition,” arXiv:1405.1471 (2014). Kartiek Agarwal, Sarang Gopalakrishnan, Michael Knap,Markus M¨uller, and Eugene Demler, “Anomalous diffu-sion and griffiths effects near the many-body localizationtransition,” Phys. Rev. Lett. , 160401 (2015). Sarang Gopalakrishnan, Markus M¨uller, Vedika Khemani,Michael Knap, Eugene Demler, and David A. Huse, “Low-frequency conductivity in many-body localized systems,”Phys. Rev. B , 104202 (2015). A. Chandran, C. R. Laumann, and V. Oganesyan, “Finitesize scaling bounds on many-body localized phase transi-tions,” arXiv:1509.04285 (2015). Sarang Gopalakrishnan, Kartiek Agarwal, Eugene A.Demler, David A. Huse, and Michael Knap, “Griffithseffects and slow dynamics in nearly many-body localizedsystems,” Phys. Rev. B , 134206 (2016). Kartiek Agarwal, Ehud Altman, Eugene Demler, SarangGopalakrishnan, David A. Huse, and Michael Knap,“Rare region effects and dynamics near the many-body localization transition,” Ann. Phys. (2017),10.1002/andp.201600326. Vedika Khemani, S.P. Lim, D.N. Sheng, and David A.Huse, “Critical properties of the many-body localiza-tion transition,” Phys. Rev. X (2017), 10.1103/Phys-RevX.7.021013. Sarang Gopalakrishnan and Rahul Nandkishore, “Mean-field theory of nearly many-body localized metals,” Phys.Rev. B , 224203 (2014). Ronen Vosk, David A. Huse, and Ehud Altman, “The-ory of the many-body localization transition in one-dimensional systems,” Phys. Rev. X , 031032 (2015). Andrew C. Potter, Romain Vasseur, and S. A.Parameswaran, “Universal properties of many-body delo-calization transitions,” Phys. Rev. X , 031033 (2015). Liangsheng Zhang, Bo Zhao, Trithep Devakul, andDavid A. Huse, “Many-body localization phase transition:A simplified strong-randomness approximate renormaliza- tion group,” Phys. Rev. B , 224201 (2016). David A. Huse, Rahul Nandkishore, and VadimOganesyan, “Phenomenology of fully many-body-localizedsystems,” Phys. Rev. B , 174202 (2014). Maksym Serbyn, Z. Papi´c, and Dmitry A. Abanin, “Lo-cal conservation laws and the structure of the many-bodylocalized states,” Phys. Rev. Lett. , 127201 (2013). M. Serbyn, M. Knap, S. Gopalakrishnan, Z. Papi, N.Y.Yao, C.R. Laumann, D.A. Abanin, M.D. Lukin, and E.A.Demler, “Interferometric probes of many-body localiza-tion,” Phys. Rev. Lett. , 147204 (2014). Yasaman Bahri, Ronen Vosk, Ehud Altman, and AshvinVishwanath, “Localization and topology protected quan-tum coherence at the edge of hot matter,” , 7341 (2015). A. Chandran, A. Pal, C. R. Laumann, and A. Scardicchio,“Many-body localization beyond eigenstates in all dimen-sions,” Phys. Rev. B , 144203 (2016). Jens H. Bardarson, Frank Pollmann, and Joel E. Moore,“Unbounded growth of entanglement in models of many-body localization,” Phys. Rev. Lett. , 017202 (2012). V. Ros, M. M¨uller, and A. Scardicchio, “Integrals of mo-tion in the many-body localized phase,” Nuclear PhysicsB , 420 – 465 (2015). Louk Rademaker and Miguel Ortu˜no, “Explicit local inte-grals of motion for the many-body localized state,” Phys.Rev. Lett. , 010404 (2016). Anushya Chandran, Isaac H. Kim, Guifre Vidal, andDmitry A. Abanin, “Constructing local integrals of mo-tion in the many-body localized phase,” arXiv:1407.8480(2014), arXiv: 1407.8480. T. E. O’Brien, Dmitry A. Abanin, Guifre Vidal, andZ. Papi´c, “Explicit construction of local conserved opera-tors in disordered many-body systems,” Phys. Rev. B ,144208 (2016). Rong-Qiang He and Zhong-Yi Lu, “Interaction-inducedcharacteristic length in strongly many-body localized sys-tems,” arXiv:1606.09509 (2016). S. J. Thomson and M. Schir´o, “Time Evolution of Many-Body Localized Systems with the Flow Equation Ap-proach,” ArXiv e-prints (2017), arXiv:1707.06981 [cond-mat.dis-nn]. M. Goihl, M. Gluza, C. Krumnow, and J. Eisert, “Con-struction of exact constants of motion and effective modelsfor many-body localized systems,” ArXiv e-prints (2017),arXiv:1707.05181 [cond-mat.stat-mech]. Stephen Inglis, “Accessing many-body localized statesthrough the generalized gibbs ensemble,” Phys. Rev. Lett. (2016), 10.1103/PhysRevLett.117.120402. A. Chandran, J. Carrasquilla, I. H. Kim, D. A. Abanin,and G. Vidal, “Spectral tensor networks for many-bodylocalization,” Phys. Rev. B , 024201 (2015). David Pekker and Bryan K. Clark, “Encoding the structureof many-body localization with matrix product operators,”Phys. Rev. B , 035116 (2017). Frank Pollmann, Vedika Khemani, J. Ignacio Cirac, andS. L. Sondhi, “Efficient variational diagonalization of fullymany-body localized hamiltonians,” Phys. Rev. B ,041116 (2016). Thorsten B. Wahl, Arijeet Pal, and Steven H. Simon, “Ef-ficient representation of fully many-body localized systemsusing tensor networks,” Phys. Rev. X , 021018 (2017). Hyungwon Kim, Mari Carmen Ba˜nuls, J. Ignacio Cirac,Matthew B. Hastings, and David A. Huse, “Slowest lo-cal operators in quantum spin chains,” Phys. Rev. E , 012128 (2015). F. Verstraete, V. Murg, and J.I. Cirac, “Matrix prod-uct states, projected entangled pair states, and varia-tional renormalization group methods for quantum spinsystems,” Advances in Physics , 143–224 (2008),http://dx.doi.org/10.1080/14789940801912366. Ulrich Schollw¨ock, “The density-matrix renormalizationgroup in the age of matrix product states,” Annals ofPhysics , 96 – 192 (2011), january 2011 Special Issue. Yevgeny Bar Lev, Guy Cohen, and David R. Reichman,“Absence of diffusion in an interacting system of spinlessfermions on a one-dimensional disordered lattice,” Phys.Rev. Lett. , 100601 (2015). Marko ˇZnidariˇc, Antonello Scardicchio, and Vipin KeralaVarma, “Diffusive and Subdiffusive Spin Transport in theErgodic Phase of a Many-Body Localizable System,” Phys.Rev. Lett. , 040601 (2016). Laurens de Haan and Ana Ferreira, “Extreme value the-ory: An introduction,” Springer-Verlag New York (2006),10.1007/0-387-34471-3. David J. Luitz, Nicolas Laflorencie, and Fabien Alet,“Many-body localization edge in the random-field heisen-berg chain,” Phys. Rev. B , 081103 (2015). Marcin Mierzejewski, Peter Prelovˇsek, and Toma ˇz Prosen,“Identifying local and quasilocal conserved quantities inintegrable systems,” Phys. Rev. Lett. , 140601 (2015). M. Mierzejewski, M. Kozarzewski, and P. Prelovsek,“Counting local integrals of motion in disordered spinless-fermion and Hubbard chains,” ArXiv e-prints (2017),arXiv:1708.08931 [cond-mat.str-el]. C.-J. Lin and O. I. Motrunich, “Explicit construction ofquasi-conserved local operator of translationally invariantnon-integrable quantum spin chain in prethermalization,”ArXiv e-prints (2017), arXiv:1709.02135 [cond-mat.stat-mech]. The traceless condition ensures the solution will have nooverlap with the trivially commuting identity. In general we expect that a numerical minimization for anygiven support, will find operators slower than the polariza-tion fluctuations, even in the clean system. . Jean-Philippe Bouchaud and Marc M´ezard, “Universalityclasses for extreme-value statistics,” Journal of Physics A:Mathematical and General , 7997 (1997). Giulio Biroli, Jean-Philippe Bouchaud, and Marc Potters,“Extreme value problems in random matrix theory andother disordered systems,” Journal of Statistical Mechan-ics: Theory and Experiment , P07019 (2007). R´obert Juh´asz, Yu-Cheng Lin, and Ferenc Igl´oi, “Stronggriffiths singularities in random systems and their rela-tion to extreme value statistics,” Phys. Rev. B , 224206(2006). Arvydas Astrauskas, “From extreme values of i.i.d. ran-dom fields to extreme eigenvalues of finite-volume andersonhamiltonian,” Probab. Surveys , 156–244 (2016). A. K. Kulshreshtha, A. Pal, T. B. Wahl, and S. H. Si-mon, “Behaviour of l-bits near the many-body localiza-tion transition,” ArXiv e-prints (2017), arXiv:1707.05362[cond-mat.dis-nn]. The reason for this efficient computation is that the weight | c i ( k ) | can be expressed in terms of the vectorized oper-ator O M as the expectation value of a superoperator withtensor product structure and can then be evaluated withMPS primitives. There will actually be exponentially many conserved quan-tities, which can be constructed by any product or combi-nation of LIOMs τ zi . David Pekker, Bryan K. Clark, Vadim Oganesyan, andGil Refael, “Fixed points of wegner-wilson flows and many-body localization,” Phys. Rev. Lett. , 075701 (2017). Rahul Nandkishore, Sarang Gopalakrishnan, andDavid A. Huse, “Spectral features of a many-body-localized system weakly coupled to a bath,” Phys. Rev.B , 064203 (2014). Mark H Fischer, Mykola Maksymenko, and Ehud Altman,“Dynamics of a many-body-localized system coupled to a bath,” Phys. Rev. Lett. , 160401 (2016). Emanuele Levi, Markus Heyl, Igor Lesanovsky, andJuan P Garrahan, “Robustness of many-body localiza-tion in the presence of dissipation,” Phys. Rev. Lett. ,237203 (2016). Mariya V Medvedyeva, Tomaˇz Prosen, and MarkoˇZnidariˇc, “Influence of dephasing on many-body localiza-tion,” Phys. Rev. B , 094205 (2016). Sarang Gopalakrishnan, K. Ranjibul Islam, and MichaelKnap, “Noise-induced subdiffusion in strongly localizedquantum systems,” Phys. Rev. Lett.119