Comparing many-body localization lengths via non-perturbative construction of local integrals of motion

TL;DR

This paper introduces a non-perturbative method to explicitly construct local integrals of motion for many-body localized systems, enabling detailed analysis of localization lengths and system dynamics.

Contribution

The authors develop a non-perturbative approach to compute complete LIOMs, allowing direct mapping of Hamiltonians and analysis of MBL properties at the microscopic level.

Findings

LIOMs exhibit exponential decay in real space

Localization lengths can be extracted from LIOM weights

Method is effective even at phase transition points

Abstract

Many-body localization (MBL), characterized by the absence of thermalization and the violation of conventional thermodynamics, has elicited much interest both as a fundamental physical phenomenon and for practical applications in quantum information. A phenomenological model, which describes the system using a complete set of local integrals of motion (LIOMs), provides a powerful tool to understand MBL, but can be usually only computed approximately. Here we explicitly compute a complete set of LIOMs with a non-perturbative approach, by maximizing the overlap between LIOMs and physical spin operators in real space. The set of LIOMs satisfies the desired exponential decay of weight of LIOMs in real-space. This LIOM construction enables a direct mapping from the real space Hamiltonian to the phenomenological model and thus enables studying the localized Hamiltonian and the system dynamics. We can thus study and compare the localization lengths extracted from the LIOM weights, their interactions, and dephasing dynamics, revealing interesting aspects of many-body localization. Our scheme is immune to accidental resonances and can be applied even at phase transition point, providing a novel tool to study the microscopic features of the phenomenological model of MBL.

Equations

$$H=\sum _{i=1}^L{h}_i{\sigma }_z^i+\sum _{i=1}^{L-1}{\vec{\sigma }}^i\cdot {\vec{\sigma }}^{i+1},$$

$${\tau }_z^j=\sum _{n=0}^L\sum _k{f}_{n,k}^j{\hat{O}}_{n,k}^j,$$

$$H=\sum _i{\xi }_i{\tau }_z^i+\sum _{ij}{V}_{ij}{\tau }_z^i{\tau }_z^j+\sum _{ijk}{V}_{ijk}{\tau }_z^i{\tau }_z^j{\tau }_z^k+\cdots .$$

$${||{\sigma }_z^j-{\tau }_z^j||}_{W\to 0}\sim \sqrt{||{\sigma }_z^j|{|}^2+||{\tau }_z^j|{|}^2}=2^{(L+1)/2}.$$

$$H_j=\mathrm{Tr}({\tau }_z^{j}H)={\xi }_j+H_j^1+H_j^2+\cdots ,$$

$$H_j^1=V_{j,j+1}{\tau }_z^{j+1}+V_{j-1,j}{\tau }_z^{j-1}+V_{j-1,j,j+1}{\tau }_z^{j-1}{\tau }_z^{j+1}.$$

$$||H_j^l||={\left[\sum _{r=l+1}^{2l+1}{2}^r{e}^{-2r/\kappa }\right]}^{1/2}\approx \frac{{\left(2e^{-2/\kappa }\right)}^{(l+1)/2}}{1-2e^{-2/\kappa }}.$$

$$||\Delta {\tau }_z^j||=\sqrt{\mathrm{Tr}(\Delta {\tau }_z^{j†}\Delta {\tau }_z^j)},$$

Full text

I Introduction

How a many-body quantum system thermalizes –or fails to do so– under its own interaction is a fundamental yet elusive problem. Localization serves as a prototypical example for the absence of thermalization, first studied in the non-interacting single particle regime known as Anderson localization [1; 2], and then revived in the context of interacting systems (many-body localization, MBL) [3]. The existence of MBL as a phase of matter was demonstrated theoretically [4; 5; 6] and numerically [7; 8; 9; 10; 11]. Recently, the MBL phase was observed in cold atoms [12; 13; 14; 15; 16; 17], trapped ions [18; 19] and natural crystals using nuclear magnetic resonances [20]. Most characteristics of MBL, such as area law entanglement [21; 22], Poisson level statistics [8; 9], logarithmic growth of entanglement [23; 7; 24; 25; 26; 16; 27] and power law dephasing [28; 29; 30; 31; 32], can be understood via a phenomenological model that expresses the Hamiltonian in terms of a complete set of local integrals of motion (LIOMs) [25; 21]. However, the explicit computation of LIOMs and their interactions is a challenging task, complicated by the fact that the set of LIOMs is not unique. LIOMs have been calculated by the infinite-time averaging of initially local operators [33; 34], however, the obtained LIOMs does not have binary eigenvalues and thus cannot form complete basis. Binary-eigenvalue LIOMs can be obtained using perturbative treatment of interactions [5; 35; 36; 37; 38], Wegner-Wilson flow renormalization [39], minimizing the commutator with the Hamiltonian [40]. The previous methods either requires strong disorder field strength, or assumes a cutoff of LIOMs in real space, so a complete numerical study of localization lengths is missing.

Here we design and implement a method to compute a complete set of binary LIOMs (i.e., with eigenvalues $\pm 1$) in a non-perturbative way, by maximizing the overlap with physical spin operators. This criterion enables a recursive determination, similar to quicksort, of the LIOMs matrix elements in the energy eigenbasis, without the need to exhaust all the eigenstate permutations, which would be prohibitive for system size $L>5$. We verify that in the MBL phase the LIOMs are exponentially localized in real space, and their interaction strength decays exponentially as a function of interaction range. LIOM localization lengths and interaction localization lengths can be extracted from the two exponential behavior respectively. Deep in the MBL phase, the two localization lengths are well characterized by the inequality derived in Ref. [3]. Near the transition point, that our construction enables exploring, the interaction localization length diverges, while the LIOM localization length remains finite: this should be expected given the constraints imposed by our construction, even if it contradicts the inequality in Ref. [3]. The explicit form of the LIOMs further enables exploring the system evolution, and we show that the LIOMs display a similar dynamics to the physical spin operators, and extract the dynamical localization length from the power law dephasing process [28]. Interestingly, we find that the dynamical localization length is much shorter than would be given by a conjectured relationship to the above two localization lengths [28; 3], suggesting that the dynamics does not only depend on the typical value of LIOMs and their interactions, but also on higher order correlations.

II Algorithm

To understand the construction algorithm, we first review the properties of integrals of motions in the many-body localized phase. LIOMs $\{\tau_{z}^{j}\}$ are diagonal in the Hamiltonian eigenbasis $[H,\tau_{z}^{j}]=0$. A complete set of LIOMs can be related to physical spin operators by a local unitary transformation $\tau_{z}^{j}=U\sigma_{z}^{j}U^{\dagger}$, which implies that (i) half of the eigenvalues of $\tau_{z}^{j}$ are +1 and the other half are -1; (ii) LIOMs are mutually independent (orthonormal) $\mathrm{Tr}(\tau_{z}^{j}\tau_{z}^{k})/2^{L}=\delta_{jk}$; (iii) the weight of $\tau_{z}^{j}$ decays exponentially in real space for localized Hamiltonians. In particular, property (ii) requires that, for any $j$, in either +1 or -1 sector of $\tau_{z}^{j}$, half of the diagonal elements of $\tau_{z}^{k}$ are +1 and the other half are -1 for all $k\neq j$. In another word, the +1 and -1 sectors of $\tau_{z}^{j}$ are effectively two manifolds that represent two instances of a new system with $L-1$ spins, containing all sites except $j$.

With only constraints (i-ii), there are $2^{L}!/L!$ different sets of IOMs among which we want to find the most local one. However, enumerating the $2^{L}!/L!$ different sets, and quantifying the localization of the related $\tau_{z}^{j}$, is numerically prohibitive. Therefore, instead of explicitly demanding the exponential localization, we maximize the overlap of LIOMs and physical spin operators $\mathrm{Tr}(\tau_{z}^{j}\sigma_{z}^{j})$, which enables a systematic and efficient way to find a unique set of LIOMs, and then we verify that these LIOMs are indeed exponentially localized in the MBL phase.

Expanding the IOMs $\tau_{z}^{j}$ in the energy eigenbasis $\{\ket{n}\}$, $n=1,2,\cdots,2^{L}$, as $\tau_{z}^{j}=\sum_{n}a^{j}_{n}|n\rangle\langle n|$, our goal is to find $a_{n}^{j}\in\pm 1$ under the constrains (i-iii). (We thus assume that we have diagonalized the Hamiltonian). The algorithm is reminiscent of quicksort (see Fig. 1):

For all eigenstates $|n\rangle$ and spin $j$, evaluate $s_{n}^{j}=\langle n|\sigma_{z}^{j}|n\rangle$. 2. 2.

For each $j$, sort the eigenstates according to $s_{n}^{j}$, and define candidates $\tilde{\tau}_{z}^{j}=\sum_{n\in S^{j}_{max}}\ket{n}\bra{n}-\sum_{n\in S^{j}_{min}}\ket{n}\bra{n}$, where $S^{j}_{max(min)}$ is the set of eigenstates giving the $2^{L}/2$ largest (smallest) overlaps $s_{n}^{j}$. 3. 3.

For each $j$, compute the overlaps $\langle\tilde{\tau}_{z}^{j}\sigma_{z}^{j}\rangle=\sum_{n\in S^{j}_{max}}s_{n}^{j}-\sum_{n\in S^{j}_{min}}s_{n}^{j}$ and find the site $j_{M}$ that maximizes it. For this site, set $\tau_{z}^{j_{M}}\equiv\tilde{\tau}_{z}^{j_{M}}$. 4. 4.

Consider the two manifolds $\mathbb{S}^{j_{M}}_{\pm}$ corresponding to the $\pm 1$ eigenstates of $\tau_{z}^{j_{M}}$. Each of these manifolds represents two instances of a new system with $L-1$ spins, containing all sites except $j_{M}$. In this new system, perform the same protocol 1-3 to set another LIOM. This results in 4 sectors, each containing $2^{L-2}$ states. 5. 5.

By repeating the previous steps $L-2$ times we finally reduce the dimension of each sector to just 1 and all $a_{n}^{j}$ are assigned.

We note that our scheme does not necessarily find the most local set of $\tau_{z}^{j}$, since once the matrix elements of a LIOM are determined at a given step, the subsequent search for the rest of the LIOMs is restricted to its perpendicular complement to satisfy orthogonality (that is, we are not ensured to find a global optimum). Therefore, we choose to divide sectors using the most local LIOM (largest $\langle\tilde{\tau}_{z}^{j}\sigma_{z}^{j}\rangle$), so that this division sets the least constrains to later divisions. In Fig. 12 of the Appendix/SM we show that this choice indeed gives the most local results among all alternate algorthms we tried. Because we only utilize the overlaps $s_{n}^{j}=\langle n|\sigma_{z}^{j}|n\rangle$ in the computation, the scheme is immune to accidental resonances in the spectrum.

III Results

III.1 Localization of operators and interactions

To test the proposed algorithm and characterize the LIOMs that it finds we consider a prototypical example of an MBL-supporting system, a Heisenberg spin-1/2 chain with random fields,

[TABLE]

where $h_{i}$ is uniformly distributed in $[-W,W]$. It is known [8] that in the thermodynamic limit there is a MBL phase transition at $W_{c}\approx 7\pm 2$. Although this model conserves the total magnetization along $z$, the validity of the algorithm does not depend on this symmetry. To quantitatively check the locality of LIOMs, we decompose them into tensor products of Pauli operators

[TABLE]

where $\hat{O}^{j}_{n,k}$ is a tensor product of Pauli operators whose furthest non-identity Pauli matrix from $j$ is of distance $n$, e.g. $\sigma_{x}^{1}\otimes\sigma_{x}^{2}\otimes\sigma_{y}^{3}\otimes\mathbb{I}^{4}$ is of distance $n=2$ to $j=1$, because $\sigma_{y}^{3}$ is the furthest non-identity Pauli matrix. $k$ labels operators with the same $n$. $f_{n,k}^{j}=\mathrm{Tr}(\tau_{z}^{j}\hat{O}^{j}_{n,k})$ is the weight of $j$-th LIOM on $\hat{O}^{j}_{n,k}$. Figures 2(a) and (b) show the median of $|f_{n.k}^{j}|$ as a function of distance $n$. In the MBL phase, the median weight decays exponentially with distance $n$, while in the ergodic phase it saturates at large $n$.

Because the LIOMs form an orthonormal basis, the Hamiltonian can be decomposed into this basis unambiguously and efficiently:

[TABLE]

For non-interacting models, only the $\xi_{i}$ coefficients are nonzero. We can define the range $r$ of each coupling term $V_{ij\cdots}$ as the largest difference among the indices. For example, the range for 2-body interaction $V_{ij}$ is simply $r=|i-j|$, while for 3-body interactions is $r=\mathrm{max}(|i-j|,|i-k|,|j-k|)$. Figures 2(c) and (d) show the median interaction strength as a function of interaction range. In the MBL phase, the interaction strength decays exponentially. The behavior of two-body interactions $|V_{ij}|$ and three body interactions $|V_{ijk}|,\cdots$ show no significant difference [36; 39] and can be essentially captured by the median of all interaction terms for a given range $V(r)$. We considered the median instead of the mean in order to exclude rare events, i.e., instances where the disorder strength is small in a local region.

To gain more insight into the localization of IOMs and interactions and observe the occurrence of rare events, in Figure 3 we further study the probability distribution of weight $f_{n,k}^{j}$ versus $n$, and the probability distribution of interaction strength $V(r)$ versus $r$ in the localized regime (strong disorder). The distribution of $\log(|f_{n,k}^{j}|)$ can be described by a single Gaussian peak, centered at smaller values of $|f_{n,k}^{j}|$ when the distance $n$ increases, confirming the localization of IOMs. Instead, two peaks can be observed in the distribution of $\log(|V|)$. The left peak shifts to smaller $|V|$ with increasing $r$, while the right peak (larger $|V|$) shows no significant shift. Moreover, the area of the right peak decreases for larger $W$ and smaller $L$. Therefore, we identify the left peak as describing localized cases, the right one as rare events. The exponential localization of the LIOMs and their interactions are usually the two criteria that define the LIOM. In the rare region of low disorder, however, the two requirements cannot be satisfied simultaneously and there is no universal criteria on how to choose LIOMs in this case. Here we require the IOM $\tau_{z}$ to be local by construction, so the presence of a rare region shows up only in the interaction strengths; choosing different criteria for the LIOM construction may lead to different results.

III.2 Localization lengths

From the explicit form of the LIOMs and their interactions, we can extract the LIOM localization length $\xi$, via $|f_{n,k}^{j}|\sim\exp(-n/\xi)$, and interaction localization length $\kappa$, via $|V(r)|\sim\exp(-r/\kappa)$ [3]. In Figure 4 we show $\kappa$ and $\xi$ as a function of disorder strength $W$. The LIOM localization length $\xi$ is extracted using the relation $\mathrm{Tr}(\tau_{z}^{j}\sigma_{z}^{k})\sim\exp(-|k-j|/\xi)$ [33; 36] because calculating $f_{n,k}^{j}$ is numerically demanding (see SM). The interaction localization length $\kappa$ is extracted by fitting the distribution of $\log|V|$ (as in Fig. 3) to two Gaussian peaks and then fitting the localized peak center to a linear function of $r$. Because our method forces $\tau_{z}$ to be local, $\xi$ is always finite, while $\kappa$ diverges around $W=8.1$ [Fig. 4(b)], which agrees with the critical point $W_{c}=7\pm 2$ reported in Ref. [8]. It has been shown in [3] that the two localization lengths satisfy the inequality $\kappa^{-1}\geq(\xi^{-1}-\mathrm{ln}2)/2$. From the numerical results in Fig. 4(c), we find that this inequality is satisfied in the localized phase, except in the vicinity of the phase transition point.

III.3 Non-interacting model: tradeoff of localization

We can better understand why the interaction localization length $\kappa$ diverges at the critical point while the LIOM localization length $\xi$ remains finite by applying our LIOM construction to a non-interacting model $H=\sum_{i=1}^{L}h_{i}\sigma_{z}^{i}+\sum_{i=1}^{L-1}\left(\sigma^{i}_{x}\sigma^{i+1}_{x}+\sigma^{i}_{y}\sigma^{i+1}_{y}\right)$. Due to the lack of interactions, the system is effectively localized for arbitrarily small $W$. This Hamiltonian can be mapped to a free fermionic Hamiltonian via a Jordan-Wigner transformation [41]. The Hamiltonian can be diagonalized by single-particle IOMs $\{\Sigma_{z}^{i}\}$: $H=\sum_{i}\tilde{\xi}_{i}\Sigma_{z}^{i}$, that is, the interaction localization length in the $\{\Sigma_{z}^{i}\}$ basis is zero. However, note that the single-particle IOMs $\{\Sigma_{z}^{i}\}$ can be highly non-local for small $W$. We can instead apply our algorithm to find LIOMs $\{\tau_{z}^{j}\}$ for this model as done for the interacting Hamiltonian and compare $\{\Sigma_{z}^{j}\}$ and $\{\tau_{z}^{j}\}$ (see Fig. 5). For large disorder strength, $W=20$, the Hamiltonian is practically interaction-free even in the $\tau_{z}^{j}$ basis, and indeed the LIOMs $\tau_{z}^{j}$ approach the IOMs, $\tau_{z}^{j}\approx\Sigma_{z}^{j}$. The trade off between the two interaction strength $\kappa$ and $\xi$ becomes evident for small disorder, $W=0.5$, where $\tau_{z}^{j}\neq\Sigma_{z}^{j}$. In this regime, the single-particle IOMs $\Sigma_{z}^{j}$ are delocalized, $\xi\gg 1$, but the Hamiltonian still has no interactions, $\kappa=0$. Instead, the LIOMs obtained by our construction, $\{\tau_{z}^{j}\}$, are localized but they give rise to long-range interactions in the Hamiltonian, $\kappa\gg 1$. For interacting models, it is difficult to obtain IOMs that minimize the interactions in a non-perturbative way. Still, we expect that if one were indeed able to find such a set of IOMs, there would be a similar tradeoff between how local they are (small $\xi$) versus how local the interactions are (small $\kappa$) outside the well-localized phase. Our choice of criterion for constructing LIOMs not only allows a simple and efficient algorithm; by keeping the operators local even when crossing the localization transition, the $\tau_{z}^{j}$ are always well-defined and can be used to explore properties of the system, such as its dynamics, around the localization-delocalization transition point.

III.4 Dephasing Dynamics

Since physical spin and LIOM operators are related by a local unitary transformation, they are expected to exhibit a similar dynamics [Fig. 4(a)]. In particular, the higher order interaction terms in Eq. (3) induce dephasing of the transverse operators by creating an effective magnetic field $H_{\mathrm{eff}}$ at the location of spin $j$ due to all the other spins. The dephasing of the expectation values $\langle\tau_{x}(t)\rangle$ and $\langle\sigma_{x}(t)\rangle$ is closely related to the logarithmic light cone in the MBL phase [28]. It was previously shown that $\langle\langle\sigma_{x}(t)\rangle^{2}\rangle\approx\langle\langle\tau_{x}(t)\rangle^{2}\rangle\propto t^{-\alpha}$, where we took the average of the expectation values over random initial states and disorder realizations. For an initial state given by a product state with each individual spin pointing randomly in the xy plane, $\alpha=2\xi^{\prime}\ln 2$ for bulk spins and $\alpha=\xi^{\prime}\ln 2$ for boundary spins, where $\xi^{\prime}$ is a localization length different from $\xi$ and $\kappa$ [28]. This length $\xi^{\prime}$, that we name dynamical localization length, describes the strength of the contribution to the effective magnetic field felt by spin $j$ due to spins at distance $l$: $H_{\mathrm{eff}}^{l}\sim\exp(-l/\xi^{\prime})$ (see SM). By assuming exponentially decaying interactions, $|V(r)|=\exp(-r/\kappa)$, it was conjectured that $\xi^{\prime-1}\geq\kappa^{-1}+(\ln 2)/2$ [3]. We find instead a much larger dephasing rate [Fig. 4(c)]. To investigate whether this is due solely to our LIOMs construction which does not explicitly enforces an exponentially decaying interaction strength, we artificially generate an Hamiltonian satisfying $|V(r)|\propto\exp(-r/\kappa)$ (see SM). Still, although we indeed find a power law decay, this is even faster than what we observe in Fig. 4(a). We conjecture that the dephasing process cannot be simply described by a mean interaction strength (the model used to justify the relationship to $\kappa$), and higher order correlations may play an important role.

IV Conclusion and Outlook

We provide a novel method to efficiently compute the LIOMs for MBL systems by maximizing the overlap between LIOMs and physical spin operators. The method is non-perturbative and thus immune to resonances in the spectrum, and can be applied at the phase transition point. The only quantity we use in computing the LIOMs and their interactions is the expectation value of physical spin operators on energy eigenstates $\langle n|\sigma_{z}^{j}|n\rangle$. Although we use exact diagonalization here, our scheme is compatible with renormalization group methods and matrix product state representations [42; 38], which can potentially be applied to much larger system and beyond one dimension. We show the power of the constructed LIOMs by extracting the localization length of the LIOMs and the Hamiltonian interactions from their respective exponential decays. We also show that in the MBL phase, the LIOMs and physical spin operators exhibit similar dephasing dynamics, even if it cannot be simply explained by the typical weights of LIOMs and typical interaction strengths.

Appendix A Comparison of LIOMs and physical spin operators

In the main text we defined the overlap $f_{n,k}^{j}$ as a quantifier of the locality of the LIOMs $\tau_{z}^{j}$. Another metric that characterizes the LIOMs as a function of disorder strength is the distance of each $\tau_{z}^{j}$ from the corresponding physical spin-1/2 Pauli operator $\sigma_{z}^{j}$. Indeed, the larger the disorder, the more local are the LIOMs, and therefore the closer to the corresponding Pauli operators. We use the Frobenius norm of the matrix difference between the two operators at the same site [see Fig. 6(a)] to quantify the operator distance. At small disorder strength, the LIOM and physical spin operators are almost perpendicular,

[TABLE]

As the disorder strength increases, the distance decreases, as expected. At strong disorder strength $W>W_{c}\sim 7$, we find that the distance decreases as $1/W$, indicating that the system is in the MBL phase. This result shows that the Frobenius norm distance (or equivalently the trace norm) can be taken as good proxy for the overlap $f_{n,k}^{j}$.

In the main text we state that the LIOM localization length can be extracted from $\mathrm{Tr}(\tau_{z}^{j}\sigma_{z}^{k})$. To confirm this quantitatively, in Fig. 6(b) we compare the weight of first LIOM $f_{n,k}^{1}$ and $\text{Tr}(\tau_{z}^{1}\sigma_{z}^{1+n})$. Both of them show exponential decay with $n$ and the slopes (decay rates) are similar for $n\geq 2$. In numerics, calculating $f_{n,k}^{1}$ is demanding because it is defined in the real space (see SM), while calculating $\text{Tr}(\tau_{z}^{1}\sigma_{z}^{1+n})$ can be done in the energy eigenbasis since the expectation value of $\sigma_{z}^{j}$ on every eigenstate is already obtained during the construction process. Therefore we use $\text{Tr}(\tau_{z}^{j}\sigma_{z}^{1+n})$ with $n\geq 2$ to extract the LIOM localization length $\kappa$ in the main text (Figure 4).

Appendix B Distribution of interaction strengths and rare regions

In the main text we linked the occurrence of rare events in the distribution of interaction strengths to rare regions of the disordered field. We can verify this conjecture by taking a closer look at one particular disorder realization that contains a low-disorder rare region (see Fig. 9 in the SM). To further confirm the connection between a rare region and the rare event peak in the interaction strength distribution, we study a Heisenberg spin chain with disorder field only on part of the chain, i.e. $h_{i}\in[-W,W]$ in the disordered region $i\leq i_{e}$ and $h_{i}=0$ in the disorder-free region $i>i_{e}$. The LIOMs in the disordered region are localized, while the LIOMs in the disorder-free region are delocalized with an exponential tail extending into the disordered region [Fig. 7(a)]. Due to the existence of the disorder-free region, the probability distribution of $\log|V(r)|$ shows a large delocalized peak [blue curve in Fig. 7(b)], which is absent when considering only interactions inside the disordered region. We can further analyze how the occurrence of rare events changes with the system size (see Figure 10 of the SM). We find that for a given interaction range, the area of the delocalized peak gets larger for longer chain, because the frequency of having a local low-disorder region is higher for larger $L$.

Appendix C Dephasing with Artificial Hamiltonian

It has been conjectured that the dephasing rate of $\langle\tau_{x}\rangle$ (and $\langle\sigma_{x}\rangle$) can be related to the interaction localization length via a simple, mean-field model. Using our LIOM construction, we found instead surprising results as shown in Figure 4. Here we want to (i) verify whether assuming an exponentially decaying interaction strength does indeed yield the relationship between localization lengths presented in Ref.[3]; and (ii) determine whether the Hamiltonian approximation with a simpler, exponentially decaying interaction strength is enough to capture the exact dephasing dynamics.

To answer these questions, we consider two artificially generated Hamiltonian: (1) $|V(r)|=\exp(-r/\kappa)$, with each interaction term randomly assigned a plus or minus sign, and (2) $|V(r)|$ randomly sampled from the simulated probability distribution (see Fig. 3(b) in the main text for an example) with a random sign.

The first Hamiltonian exactly satisfies the hypothesis under which the relation between interaction localization length and dynamical localization length was derived in Ref. [3]. Therefore we fit the power law dephasing obtained under this Hamiltonian [see Fig. 8(a)] and extract the dynamical localization length $\xi^{\prime\prime}$ as done in Fig. 4 of the main text. We find that $\kappa^{-1}\approx\xi^{\prime\prime-1}+\ln 2/2$ [Fig. 8(b)], which gives a more stringent relation than the bound $\kappa^{-1}\geq\xi^{\prime\prime-1}-\ln 2/2$ given in Ref. [3]. We can provide a heuristic argument for the relation between $\xi^{\prime\prime}$ and $\kappa$, under the assumption $|V(r)|=\exp(-r/\kappa)$. As described in Ref. [28], the dephasing can be understood as arising from an effective magnetic field at site $j$ due to all other spins $\tau_{z}$. Starting from the phenomenological model in Eq. (3), the effective magnetic field at site $j$ is

[TABLE]

where $H_{j}^{l}$ denotes the magnetic field created by spins within the distance $l$ from spin $j$. For example, the first term is given by

[TABLE]

Similarly, $H_{j}^{l}$ contains interactions of range $l+1,l+2,\cdots,2l+1$. As the interaction strength decays as $|V(r)|=\exp(-r/\kappa)$ and the number of terms grows as $\sim 2^{r}$, the Frobenius norm of $H_{j}^{l}$ is estimated as

[TABLE]

In the last term we assumed that $l\gg 1$ and the system is deep in the MBL phase so that $2\exp(-2/\kappa)<1$. We thus find that $H_{j}^{l}$ also exhibits an exponential decay $H_{j}^{l}\propto\exp(-l/\xi^{\prime\prime})$, with $\xi^{\prime\prime-1}=\kappa^{-1}+\ln 2/2$, yielding the dephasing [28] $\langle\langle\tau_{x}^{L}(t)\rangle^{2}\rangle\sim t^{-\xi^{\prime\prime}\ln 2}$ as shown in Fig. 8.

While we confirm that the dephasing under the approximated Hamiltonian satisfying $|V(r)|=\exp(-r/\kappa)$ follows the predicted relation to $\kappa$, we still find that dephasing under the “real” Hamiltonian is different. The physical spin and LIOMs under the real Hamiltonian in Eq. 1 show similar dephasing as expected. Under either artificial Hamiltonians, however, $\langle\langle\tau_{x}^{L}\rangle^{2}\rangle$ dephases much faster than under the real Hamiltonian, suggesting that the dephasing dynamics cannot be fully captured by the interaction localization length $\kappa$ or even the probability distribution of $|V(r)|$ [Fig. 8(a)]. For instance, in the real system for a given disorder realization the interaction terms may have some correlation, which gives rise to a slower dephasing, but this is not captured by the probability distribution.

Supplementary material

C.1 Rare events in the distribution of interaction strengths and rare regions

To verify that the peak at larger $|V(r)|$ corresponds to the rare regions, we take a closer look at one particular disorder realization that contains a low-disorder rare region (Fig. 9). The rare region leads to a local mix of two LIOMs, and a peak at large interaction strengths, $|V|$, arises in the interaction distribution. In addition to study this link more generally, as done in the Appendix, we can also analyze the behavior as a function of the chain length. Figure 10 shows the probability distribution of $\log|V(r)|$ for $L=14$ and $L=10$. For a given interaction range, the area of the delocalized peak gets larger for longer chain, because the frequency of having a local low-disorder region is higher for larger $L$. There are multiple resonances on the delocalized peak, which are not yet understood.

C.2 Stability

It is known that the MBL phase, contrary to integrable systems, is robust against small perturbations of the Hamiltonian; it is thus desirable that the LIOMs display the same robustness. To show the stability of our construction, i.e. that the LIOMs $\tau_{z}^{j}$ do not change dramatically under small perturbations in the Hamiltonian, we generate an additional disorder field on top of the original Hamiltonian in Eq. (1). We fix the additional disorder field configuration and scale it by $w$. We quantify the deviation of the new LIOM $\tau_{z}^{j,w}$ from the original $\tau_{z}^{j}$ using the Frobenius norm of the matrix difference:

[TABLE]

with $\Delta\tau_{z}^{j}=\tau_{z}^{j,w}-\tau_{z}^{j}$, and plot the median of this distance as a function of $w$ in Fig. 11. At first the distance grows linearly as the perturbation increases, before slowing down to a square root growth for intermediate perturbation, and finally saturating to a constant when the additional disorder dominates. The initial linear growth is expected from a linear expansion of the operators for small perturbation strength. The eventual saturation at $w\gtrsim W$ correspond to the limit where the two operators only overlap at the same (physical) site, thus giving the maximum distance between the two operators. In the intermediate region, the perturbation field does not only modify the eigenstates $|n\rangle$ but might also alter some matrix elements $a_{n}^{j}$, yielding the square root scaling. We note that the overall small deviation $\Delta\tau$ not only demonstrates the stability of our numerical method, but also more broadly the robustness of LIOMs in the MBL phase.

C.3 Comparison to alternative algorithms to find the LIOM set

We compare our proposed algorithm to two similar algorithms that follow however different sector division schemes:

Scheme 1: recursively divide the sector according to the $j_{M}$ that maximizes $\langle\tilde{\tau}_{z}^{j}\sigma_{z}^{j}\rangle$ within the sector. This is the scheme used in the main text.

Scheme 2: divide starting from the leftmost spin (smallest $j$), i.e. divide the $2^{L}$ states into $\mathbb{S}_{+}$ and $\mathbb{S}_{-}$ by sorting $\langle n|\sigma_{z}^{1}|n\rangle$, and then further divide $\mathbb{S}_{\pm}$ each into two sectors by sorting $\langle n|\sigma_{z}^{2}|n\rangle$, etc.

Scheme 3: First, compute and sort $\langle\tilde{\tau}_{z}^{j}\sigma_{z}^{j}\rangle$ to get a sequence $j_{1},j_{2},\cdots j_{L}$ such that $\langle\tilde{\tau}_{z}^{j_{1}}\sigma_{z}^{j_{1}}\rangle>\langle\tilde{\tau}_{z}^{j_{2}}\sigma_{z}^{j_{2}}\rangle>\cdots>\langle\tilde{\tau}_{z}^{j_{L}}\sigma_{z}^{j_{L}}\rangle$. Then divide according to this sequence, i.e. divide the $2^{L}$ states into $\mathbb{S}_{+}$ and $\mathbb{S}_{-}$ by sorting $\langle n|\sigma_{z}^{j_{1}}|n\rangle$, and then further divide $\mathbb{S}_{\pm}$ each into two sectors by sorting $\langle n|\sigma_{z}^{j_{2}}|n\rangle$, etc.

This last scheme differs from scheme 1 starting from the second division: it divides both $\mathbb{S}_{\pm}$ according to $j_{2}$, where $j_{2}$ is chosen such that it has the second largest $\langle\tilde{\tau}_{z}^{j}\sigma_{z}^{j}\rangle$; Scheme 1 instead treats $\mathbb{S}_{\pm}$ individually as two instances of a new system with $2^{L-1}$ eigenstates, so $j_{M}$ in $\mathbb{S}_{+}$ may differ from $j_{M}$ in $\mathbb{S}_{-}$ and $j_{M}$ is chosen such that it has largest $\langle\tilde{\tau}_{z}^{j}\sigma_{z}^{j}\rangle$ within the new system.

Scheme 1, our chosen algorithm, gives the most local results, especially for smaller disorder $W$ (see Fig. 12). For the same disorder $W$, scheme 2 gives less local results for $\tau_{z}^{j}$ with larger $j$, because the sector division for large $j$ is constrained by the sectors of small $j$, which is what motivated us to start with the most local spin. The LIOMs generated using scheme 3 are almost as local as using scheme 1, however with respect to the interactions, scheme 3 gives larger interaction strength than scheme 1, even if it still shows exponential decay with a comparable decay rate. These two observations suggest that scheme 3 has higher frequency of delocalized events than scheme 1.

Among the three schemes we have shown and many others we have tried, scheme 1 gives the best result, and we believe it successfully captures the localized cases due to the similarity with scheme 3. However, we cannot exclude the possibility that there might be another scheme which gives even lower frequency of delocalized cases.

One could also define LIOMs by maximizing their overlap with the corresponding physical spin operators, $\tau_{z}=\tilde{\tau}_{z}$, i.e. $a_{n}^{j}=-1$ for the $2^{L-1}$ eigenstate $|n\rangle$ with larger $\langle n|\sigma_{z}^{j}|n\rangle$ and $a_{n}^{j}=1$ for others. By requiring maximum overlap, $\{\tilde{\tau}_{z}^{j}\}$ are in principle more local than $\{\tau_{z}^{j}\}$ which are described in the main text, but they are not mutually independent and thus cannot form a complete basis. Numerical results for the two constructions are presented in Fig. 13, showing that the LIOMs generated by the two methods display no significant difference, even for moderate disorder strength $W=10$.

C.4 Computational complexity

The computational complexity of our LIOMs construction is set by the diagonalization, which is $O(2^{3L})$ (Fig. 14). We note that there are several methods to reduce this complexity and obtain an approximate diagonalization in the localized phase [42; 38].

Here we then analyze only the computation complexity of the other steps of the algorithm, which are particular to our scheme:

• •

The complexity of evaluating $\langle n|\sigma_{z}^{j}|n\rangle$ for all $j$ and $n$ is $O(L2^{2L})$ because $\sigma_{z}^{j}$ is sparse.

• •

The complexity of the recursion step (sorting eigenstates and dividing into sectors) is $O(L^{3}2^{L})$: for a sector containing $2^{N}$ states, sorting the $\langle n|\sigma_{z}^{j}|n\rangle$ for each $j$ is $N2^{N}$. For all $j$ is thus $N^{2}2^{N}$. There are $2^{(L-N)}$ such sectors. Total complexity is $\sum_{N=1}^{L}N^{2}2^{N}2^{L-N}\sim L^{3}2^{L}$.

• •

Assigning $a_{n}^{k}$ has a cost $O(L2^{L})$

• •

The decomposition of $H$ to LIOM basis is $O(L2^{L})$, because only diagonal elements are nonzero.

Figure 14 confirms that for large $L$, most of the time is spent on diagonalization.

The only operation that could lead to a complexity higher than $O(2^{3L})$ is computing LIOMs in the physical spin basis, which is required for calculating $f_{n}^{k}$. Transferring each LIOM from the eigenbasis to the physical spin basis is $O(2^{3L})$, so total is $O(L2^{3L})$. This high cost is the reason why we studied an alternate metric for the operator localization (the overlap with the physical spin operator), see A.