{'fraction_non_alphanumeric': 0.031556039173014146, 'fraction_numerical': 0.01074537540805223, 'mean_word_length': 4.642555640828856, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Almost-conserved operators in nearly many-body localized systems

We construct almost conserved local operators, that possess a minimal commutator with the Hamiltonian of the system, near the many-body localization transition of a one-dimensional disordered spin chain. We collect statistics of these slow operators for different support sizes and disorder strengths, both using exact diagonalization and tensor networks. Our results show that the scaling of the average of the smallest commutators with the support size is sensitive to Griffiths effects in the thermal phase and the onset of many-body localization. Furthermore, we demonstrate that the probability distributions of the commutators can be analyzed using extreme value theory and that their tails reveal the difference between diffusive and sub-diffusive dynamics in the thermal phase.


I. INTRODUCTION
Isolated interacting quantum systems can undergo an unconventional transition from an ergodic thermal phase in which the system approaches thermal equilibrium irrespective of the initial conditions to a many-body localized (MBL) phase in which it does not. 1-4 There exists firm evidence, from numerical, 5,6 rigorous mathematical, 7 and experimental work, [8][9][10][11][12][13][14][15] that the MBL phase is realized in strongly disordered and interacting quantum systems. Recently, the MBL transition between the thermal and the MBL phase has been started to be investigated using general arguments, [16][17][18][19][20][21][22] mean-field theory, 23 and renormalization-group schemes [24][25][26] . Since the MBL transition is inherently dynamical and can occur at arbitrary energy density, it is beyond the scope of a conventional thermodynamic description.
A possible strategy to explore MBL is to directly study properties of the Hamiltonian. In particular, a phenomenological description of the Hamiltonian in the MBL phase using an extensive set of local integrals of motion (LIOMs) has been proposed. 27,28 This approach explains, for instance, the dephasing 29-31 and entanglement dynamics 6,32 in the localized phase. Moreover, LIOMs have been constructed explicitly for some models using analytical techniques, 33 exact diagonalization [34][35][36][37][38][39] , stochastic methods, 40 and tensor networks. [41][42][43][44] Upon reducing the disorder, LIOMs become more and more extended in space and eventually cease to exist at the MBL transition.
Here, we take another route and study the dynamics of a disordered and interacting Heisenberg spin chain by constructing almost conserved local operators with finite support, which minimize the commutator with the system Hamiltonian. 45 Deep in the MBL phase the slow operators resemble LIOMs. However, our procedure can be directly extended to the thermal phase as well. The slow operators can potentially be used to study dynamical properties, because the value of the commutator gives a lower bound on the thermalization times of the corre-sponding operator. 45 We compute the slowest operators for different support sizes and disorder strengths, using both exact diagonalization and tensor networks. 46,47 Our work, rather than focusing on the phase transition, it aims at characterize the sub-diffusive region in its vicinity. Specifically, we collect statistics from numerous disorder realizations and show that not only the mean value of the smallest commutator, but also the probability distributions are sensitive to the underlying phase. The tails of the distributions are affected by the appearance of rare Griffiths regions, 17,21 which act as bottlenecks for transport and are argued to provide a simple description for the slow sub-diffusive dynamics in the thermal phase near the MBL transition. 17,48,49 We use extreme value theory (EVT) to analyze the tails of the distribution and find that it is well described by the generalized extreme value distribution. 50 Moreover, we discuss how the presence of rare regions may affect the asymptotic behavior of such distributions. Our results demonstrate that finding slow operators provides access dynamical quantities without resorting to the time evolution of a particular state.
The paper is organized as follows. In section II we specify our model. Section III describes the slow operator method, and presents the fundamental ideas of EVT employed in our analysis. The results for the structure of the slow operators are presented in section IV, the average values of their commutators in the different phases are shown in section V, and the statistical analysis of the probability distributions using EVT is presented in section VI. Finally, in section VII we summarize our findings and discuss potential extensions of our work.

II. MODEL
We consider the non-integrable spin-1/2 Heisenberg Hamiltonian with random magnetic field, where S α i are spin-1/2 operators, and the values of the transverse field h i are randomly chosen from the uniform distribution in the interval [−h, h], where h is the disorder strength. In the following, we set J = 1. This model exhibits a phase transition from a thermal to a MBL phase at h c ∼ 4. 5,51 In the MBL phase the system can be described by LI-OMs, which lead to an emergent integrability. The corresponding effective Hamiltonian can be expressed as 27,28 The LIOMs, τ z i , commute with each other and with the Hamiltonian. They are exponentially localized when written in terms of the physical spins { S i }. Moreover, the coefficients J ij and K (n) i{k}j decay exponentially with their spatial range.

A. Slow operator technique
Almost conserved quasi-local operators can identify long thermalization time scales. A way to compute such time scales is to search for slow operators that minimize the commutator with the Hamilonian 45 (see also . In particular, we may restrict the search to operators O M with finite support on M consecutive sites. Defining L(O M ) to be the squared norm of the commutator between the Hamiltonian and the (normalized) operator, the problem reduces to solving the variational minimization is the Frobenius norm, and O M is restricted to be traceless. 55 The operator that minimizes (3) is the one evolving the slowest in Heisenberg picture at time t = 0. Moreover, for a slightly perturbed infinite temperature state of the form ρ(O M ) = 1/Z + O M the thermalization time is lower bounded as t th ≥ 1/ √ λ M . 45 If the system has a conserved extensive quantity, such as the total energy, it naturally gives rise to slow operators. For a translationally invariant system, the Hamiltonian commutes with the sum of all translations of the corresponding local operator, but not with translations restricted to a finite window M . Nevertheless, the norm of the latter commutator seems to decrease polynomially with M . 45 In our particular case of Hamiltonian (1), the total polarization in the z direction, i σ z i , is conserved for any value of the disorder strength. The restriction of the sum to M sites, Z M = M i=0 sin πk M σ z i , is thus a natural slow operator for any given support M . It is easy to check that as the support increases, the corresponding commutator decays as L(Z M ) ∼ 1/M 2 . 45 This sets a reference for comparison with the commutators found by numerically solving the variational problem, Eq. (3). 56 For small support size M the optimization can be solved via exact diagonalization. This has been used to show that for model (1) there exists an extensive number of exponentially localized LIOM in the MBL phase which are not present in the thermal phase 36 . To reach larger supports the search can be restricted to operators O M of the form of a matrix product operator (MPO). Then standard tensor network techniques can be applied to solve this optimization, and the bond dimension can then be systematically increased until convergence is achieved.
In this work, we use exact and approximate MPO solutions in order to collect information on the minimal commutators and the corresponding operators. For different values of the support M and the disorder strength h, we compute λ M for several configurations of the random magnetic field (ranging from 2·10 3 for M = 12 to 10 5 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 windows over an infinite chain.
To solve the problem numerically, the minimization Eq. (3) can be interpreted as an eigenvalue problem. For fixed support size M and disorder strength h and for a particular disorder realization, there is an effective operator H (M ) eff acting on the d 2M -dimensional space with support on the chosen window, such that the (vectorized) operator that minimizes Eq. (3) corresponds to its eigenvector with lowest eigenvalue. Since the optimization is restricted to operators supported on a certain window of size M , only terms in the Hamiltonian that overlap with that region will contribute to the commutator. For our nearest-neighbor Hamiltonian this means that M + 2 sites of the Hamiltonian contribute, which we denote by H M +2 . The matrix representation of the corresponding commutator, acting on vectorized operators, can be written as , and the effective operator can be constructed as whereÕ M = 1 ⊗ O M ⊗ 1 and the trace is taken over the space of all M + 2 sites.

B. Extreme value theory
The statistical analysis of the smallest commutator contains relevant information about the phases of the model and the critical region, both in the average scaling of λ M with the support size M and in the probability distributions. In particular, the probability density function (PDF) p(λ M ) and the corresponding cumulative density function (CDF), F (λ M ) = λ M −∞ p(x)dx, will be sensitive to the presence of rare regions of atypically large disorder which emerge in the proximity of the critical point.
As we show in this work, these PDFs can be described and analyzed within the mathematical framework of Extreme Value Theory (EVT), a branch of statistics concerned with the description of rare events, 50 which is applied to floods, earthquakes or risk in finance and insurance, as well as to several branches of physics, including statistical mechanics and critical phenomena. [57][58][59] In particular, EVT deals with the asymptotic behavior of the extreme values of a sample, i.e., with the tails of the distributions. Let us consider a random variable x governed by a CDF F (x). Then EVT ensures that the maxima (equivalently the minima) of samples of F (x), properly normalized and centered, will be governed by a CDF 50 for y = ax + b with a, b ∈ R and a > 0. The Fisher-Tippett-Gnedenko theorem (Thm. 1.1.3 in 50) states that this single-parameter family, called generalized extreme value distribution (GEV), includes all possible limiting distributions for the extreme values of a sample of i.i.d. random variables. The family contains three subclasses which exhibit quite different behavior, the Gumbel (ζ = 0), Fréchet (ζ > 0) and Weibull (ζ < 0) distributions. Qualitatively, the deciding criterion for a PDF to belong to the basin of attraction of one family or another (i.e. for the extrema of the samples to be described by the corresponding family) is the form of its tails. 50 In our particular problem, as we try to find the minimum λ M for a certain configuration, we are effectively sampling from the left tail of p(Λ M ), which is the distribution of eigenvalues Λ M of H (M ) eff . Typically the eigenvalues of a matrix are not uncorrelated, and thus GEV is a priori not expected to describe the probability distribution of extreme eigenvalues 58,60 . Nevertheless, our results indicate that the distribution (5) provides indeed a very good description for our data, with the particular 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 operators. In the strong disorder regime, we expect that the slow operators correspond to some LIOMs, or more precisely to a truncated version of them since our support size is fixed. It is worth noticing here that although the slow operator method does not directly target the (truncated) LIOMs, as was done in Refs. 39 and 61, in the localized regime, truncated LIOMs and their combinations are good candidates to attain (exponentially) small commutators.
To address this question, we can analyze to which extent the slowest operators are local by examining their spatial structure in different regimes. This can be done by studying their decomposition as a sum of tensor products of single-site Pauli matri- In principle, for an operator O M found by the minimization, either exactly or as a MPO approximation, we can efficiently evaluate any single coefficient Nevertheless, due to the exponential growth of the basis dimension with the support, already for moderate values of M it is unfeasible to inspect all the individual coefficients. Instead, a more physical quantity for exploring the localized nature of the operators is the combined weight of all terms with a fixed range which are supported on a certain subregion. We define a range-k operator as a product of Pauli matrices that act non trivially on sites i to i + k. Formally, this corresponds to the operator Θ where a, b ∈ 1, 2, 3 to ensure Pauli matrices at the edges (thus imposing that the support is strictly i → i + k) and {α m } ∈ 0, 1, 2, 3.
The weight of these range-k terms can be written as We thus solve the optimization (3) and compute the weights, |c k (i)| 2 , along the chain of terms with fixed range k. Even if the support is large, as O M is written as an MPO, this quantity can be computed efficiently, i.e., with a cost that only scales polynomially in the range k and the support M . 62 We find clear differences in the structure of the operator depending on the disorder strength. This is illustrated in Fig. 1 for the MPO, that minimizes Eq. (3) for a particular disorder realization h i = h · r i , where r i are fixed random numbers for each value of h. We have chosen a support of M = 40 sites and a relatively large bond dimension of D = 100, in order to ensure that truncation errors are negligible compared to the effects we observe. The figure shows the different spatial profile of the singlesite contributions |c 0 (i)| 2 = 2 −M 3 a=1 |Tr (σ a i O M ) | 2 and two-site contributions |c 1 (i)| 2 (Eq.(6) with k = 1) for different values of the disorder strength. For strong disorder, where the LIOMs are exponentially localized, at least one of the |c 0 (i)| 2 is expected to be dominant. Indeed, we observe that for large disorder the operator is well localized around a single site, with weights that decay exponentially around this point. As the disorder strength decreases, the weights decay slower with distance, and the tails are no longer exponential. Finally for very small h the profile becomes flat and has the shape of a sine as expected in the diffusive regime. The single-site terms shown in the figure always accumulate most of the weight. We found that higher order terms Fig. 1 (inset), which have smaller contributions, exhibit a similar spatial decay. The observed profile at strong disorder is in good agreement with what one could expect from the decomposition of the LIOMs in the canonical operator basis, what would in principle allow the extraction of a localization 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 disorder realization. Hence, the average over realizations E(λ M ) will provide information about the underlying phase and also indirectly about transport.
In the following, we will optimize Eq. When increasing the disorder strength, the decay remains a powerlaw E(λ M ) ∼ M −α , however, with an increasing exponent α indicating that the transport becomes slower due to the presence of disorder, see inset of Fig. 2 (a) which shows the inverse exponent 1/α. This is consistent with the observation of sub-diffusive transport on the thermal side of the MBL transition. 17,48,49 The longer relaxation times in this regime can be interpreted as follows: rare insulating regions with larger than typical disorder act as bottlenecks for the conventional diffusive transport. 17,24,25 Upon approaching the MBL transition at h ∼ 4, 51 we expect the exponent to diverge. Our data indeed shows that the power law decay becomes very steep. Due to the finite support sizes M and due to the slow convergence of the algorithm for large disorder, we, however, find that the exponent levels off at a finite value. For practical purposes, we adopt the strategy of limiting the maximum computational effort for each search. Therefore, the result of the algorithm provides an upper bound for the true minimum λ M (h), and the corresponding α is biased towards smaller values. When crossing the MBL transition, an extensive number of exponentially localized conserved quantities, corresponding to the LIOMs τ z i in Eq.
(2), emerges. 63 Therefore, provided the window M is sufficiently large to capture most of the support of any LIOM, the smallest commutator should decay exponentially as λ M ∼ e −M/ξ . Such an exponential decay is supported by our numerical simulations for strong disorder, see Fig. 2 (b). The length scale ξ is related (in a possibly non-trivial way) to the localization length of the LIOMs. In the inset, we show that with increasing disorder h, the inverse length scale 1/ξ increases monotonically. The decays of α −1 and ξ −1 as the transition is approached are both compatible with previous estimates of the critical region, see e.g. Refs. 51 and 61.

A. Numerical results
We now study the probability distributions of the slow commutators λ M , which can be sensitive to rare events. In particular, we show that EVT provides useful mathematical tools to characterize the tails of the probability density and cumulative distribution functions. We first illustrate in Fig. 3, that the GEV distribution, obtained from differentiating Eq. (5), provides a very good fit to the collected data over a range of disorder strengths. Moreover, we find that the shape of the distribution strongly depends on the disorder value and that the peak of the distribution shifts from a finite, relatively large value for weak disorder, toward zero as h increases. In between those limits, in a regime where one expects sub-diffusive dynamics, the distribution broadens considerably.
In the following we will analyze the tails of this distribution quantitatively. To this end, we introduce a new variable, with dimensions of time, as T M = (Λ M ) −1/2 . The minimal eigenvalue λ M thus corresponds to the maximum of T M , that we call τ M , and has the physical interpretation of bounding the thermalization time of the operator ρ(O M ) from below. This transformation of variables moves the left tails of p(λ M ) that arise due to rare events, to right tails of p(τ M ). The corresponding survival functions S(τ M ) = 1 − F (τ M ) are shown in Fig. 4 on logarithmic scale. We observe that for small disorder, the survival function approaches zero superexponentially fast describing weak tails, due to the small probability of large τ M . As the disorder increases, the survival function decays with a power-law tail. The value of the disorder at which the tails change is in agreement with the crossover from diffusive to sub-diffusive dynamics.
To make this observation more quantitative, we study the variation of the shape parameter ζ from the fitted GEV, defined in Eq. (5), with the disorder strength. Our data shows that ζ shifts from values close to zero for weak disorder, to clearly positive values as the MBL phase is approached. The ζ parameter determines the type of distribution. In particular, the deviation of ζ from zero describes the crossover from a peaked distribution to one with polynomial tails. In terms of the GEV families, this corresponds to a change from a Gumbel (ζ = 0) to a Fréchet (ζ > 0) distribution. The observed qualitative change can be explained by an intermediate regime in which atypically slow operators appear for any given support M , leading to strong tails in the probability density function of τ M . The value of the disorder strength at which the shape parameter ζ starts being clearly positive is again consistent with the expected sub-diffusive regime of the thermal phase. 49

B. Interpretation in terms of Griffiths effects
Although the eigenvalues of the effective operators in our study are not expected to be uncorrelated, the results in the previous section show that GEV does indeed describe our findings accurately. Hence, in this section we use EVT arguments to show how the observed distributions can be qualitatively explained in terms of the existence of rare Griffiths regions.
We first consider the effect of rare localized regions in the thermal phase but close enough to the transition. In this situation, given a support size M , the typical values of λ M will be polynomial in M . Nevertheless, if rare regions are present that (partly) support an exponentially localized operator, they can give rise to exponentially small values of λ M . More concretely, within a fixed window M , such a rare configuration of size < M will occur with an exponentially small probability p( ) ∼ c , for some c < 1. This patch can support a localized operator, with localization length ξ, with a correspondingly small commutator Λ ∼ e − /ξ (or equivalently 59 T ∼ e /2ξ ). Strictly speaking, this will be detected as the slowest operator only if this commutator is smaller than the (polynomial) typical commutator for the complementary region of size M − . Since we are interested in the probability for very small commutators, and there is an exponential separation between the scaling of both terms, we may assume that small enough commutators will always come from such rare patches. Evidently, this can only be true for commutators below a certain (M and h dependent) threshold. With this assumption, the probability of some very small value of Λ is determined by the probability of finding a rare patch such that Λ ∼ e − /ξ , and thus we can identify p( )|d | = p (T ) |dT |. From that we find that, for the range of T that correspond to rare regions of length < M , the probability for the largest T values is polynomial, For PDFs with polynomial tails, EVT predicts that the extreme values will be governed by a Fréchet distribution with ζ > 0 (see appendix A). For very weak disorder, instead, the same rare configurations of the field will not support an exponentially localized operator, so that, even if a large effective disorder on some region may give rise to commutators that are below the typical one, they still decay at most polynomially. Consequently, the argument above does not apply and we expect the probability of the smallest values to decay faster (even exponentially). For a PDF with such properties, EVT predicts a Gumbel distribution with ζ = 0 (see appendix A). This is in agreement with the shape parameter ζ shown in Fig. 5 where, for weak disorder, we obtain values of ζ very close to zero which correspond to a Gumbel distribution.
Beyond the localization transition we expect the typical regions to (partly) support localized operators, and The values of the shape parameter ζ range from ∼ 0 (corresponding to a Gumbel distribution) to ζ ∼ 1/2 (corresponding to Fréchet distribution) as the disorder h is increased.
thus give rise to exponentially small values of λ M (correspondingly exponentially large values of τ M ), as we explicitly observed in the average values shown in Section V. In our data for strong disorder, we have also found broad tails of the PDF p(T ) which are consistent with a power law. Thus, the function describing our data resembles the Fréchet distribution. This might be explained by the fact that matrix elements in the MBL phase have been found to exhibit broad tails. 18,64 As we approach the transition from the MBL side by lowering the disorder, rare thermal inclusions can appear that potentially correspond to larger than typical commutators. Yet, our method only looks for the smallest commutator in each given window. Thus, if a support M encloses one such thermal subregion, a competition arises between the values of the commutator for the inclusion and that of the (typically) exponentially localized complement. Because we expect that a thermal inclusion gives rise to only polynomially decaying commutators, their value can only be the smallest one when the inclusion is sufficiently large in relation to the size M of the support. This causes our method to be less sensitive to rare thermal regions in the localized side of the transition.
In contrast to other numerical studies, where the presence of Griffiths effects was inferred from averaged observables, 17,48,49 with our method we may directly locate rare regions. In particular, we can obtain the disorder potential for the eigenvalues λ M that contribute to the tail of the distribution and analyze the microscopic configuration of the random field {h i } in real space. Following this procedure, we could, however, not unambiguously determine an obvious correlation between strong fluctuations of the field in real space and small commutators. It remains an open question whether it is possible to predict the location of the rare regions from the disorder landscape using a more direct method than the optimization.

VII. DISCUSSION
We have constructed slow operators with finite support by minimizing their commutator with the Hamiltonian of the system using both exact diagonalization and tensor network techniques. In particular, we have considered the Heisenberg spin chain with random magnetic field, which displays a dynamical transition from the thermal to the many-body localized phase. The scaling of the minimal commutator with support size provides information on the localization transition as well as on transport in the system without resorting to a specific initial state.
Furthermore, we have demonstrated that the tails of the probability distributions are sensitive to rare insulating regions in the thermal phase near the many-body localization transition. We have shown that the statistics of the smallest commutators can be analyzed within the mathematical framework of extreme value theory. 50 In particular, we have found that the distributions are well described by generalized extreme value functions whose shape depends on the disorder strength. By extreme value theory arguments, the observed behavior in the tails can be connected to the appearance of rare, strongly disordered regions, that give rise to atypically small minimal commutators. In particular, the disorder strength at which the distribution functions obtain power-law tails is consistent with the appearance of sub-diffusive transport. 17,48,49 We conclude that the slow operator technique combined with extreme value theory, constitutes a valuable tool for exploring microscopic mechanisms of the MBL transition. Further developments may include tailoring the optimization technique to target explicitly the rare thermal regions on the localized side, which could provide new insights into this less explored aspect of MBL physics. Another intriguing question would be how coupling an MBL system to an external bath 65-69 would change the structure of slow operators.

ACKNOWLEDGMENTS
We acknowledge financial support from ExQM and IM-PRS (N.P.) and from the Technical University of Munich -Institute for Advanced Study, funded by the German Excellence Initiative and the European Union FP7 under grant agreement 291763 (M.K.). This work was also partially funded by the European Union through the ERC grant QUENOCOBA, ERC-2016-ADG (Grant no. 742102).

Appendix A: Extreme value theory
A given PDF is said to belong to the basin of attraction of one of the extreme values distributions, namely Gumbel, Fréchet or Weibull, when the extrema are distributed according to the corresponding function. The von Mises conditions 50 establish simple criteria to determine whether p(x) belongs to one of them. In this appendix we show how the conditions apply to the particular cases of the distributions discussed in section VI B.
a. Strong disorder implies a Fréchet distribution. Rare regions in the thermal phase near the MBL transition may support localized operators. As shown in Eq. (7), the corresponding probability distribution is expected to decay as p(T ) ∼ 2ξT −2ξ| ln c |−1 . A sufficient condition 50 for a PDF p(T ) to belong to the basin of attraction of the Fréchet distribution is that the corresponding CDF F (T ) = with ζ > 0. From Eq. (7) we obtain so that and, asymptotically, which ensures the condition above with ζ −1 = 2ξ| ln c |, and thus implies a limiting Fréchet distribution.
b. Exponentially decaying tails imply a Gumbel distribution. In order to prove that a PDF belongs to the basin of attraction of the Gumbel distribution it is sufficient to check the following condition 50 We assume the simplest exponential decay for the right tail of the distribution p(T ) ∼ e −kT . From the corresponding cumulative function we get the survival proba-bility, so that Thus, the derivative vanishes, which ensures Eq.(A5).