Statistics of correlation functions in the random Heisenberg chain

TL;DR

This paper investigates the probability distributions of eigenstate correlation functions in the random-field Heisenberg chain, revealing distinct statistical features across ergodic, thermal, and many-body localized phases, and providing theoretical explanations for these behaviors.

Contribution

It provides a comprehensive analysis of correlation function distributions in the MBL transition, including a perturbation theory that explains the observed distribution features.

Findings

Gaussian distributions at weak disorder consistent with ETH

Heavy tails indicating anomalous thermalization in the thermal phase

Asymmetry and multimodality in MBL phase correlation distributions

Abstract

Ergodic quantum many-body systems satisfy the eigenstate thermalization hypothesis (ETH). However, strong disorder can destroy ergodicity through many-body localization (MBL) -- at least in one dimensional systems -- leading to a clear signal of the MBL transition in the probability distributions of energy eigenstate expectation values of local operators. For a paradigmatic model of MBL, namely the random-field Heisenberg spin chain, we consider the full probability distribution of eigenstate correlation functions across the entire phase diagram. We find gaussian distributions at weak disorder, as predicted by pure ETH. At intermediate disorder -- in the thermal phase -- we find further evidence for anomalous thermalization in the form of heavy tails of the distributions. In the MBL phase, we observe peculiar features of the correlator distributions: a strong asymmetry in $S_i^z S_{i+r}^z$ correlators skewed towards negative values; and a multimodal distribution for spin-flip correlators. A quantitative quasi-degenerate perturbation theory calculation of these correlators yields a surprising agreement of the full distribution with the exact results, revealing, in particular, the origin of the multiple peaks in the spin-flip correlator distribution as arising from the resonant and off-resonant admixture of spin configurations. The distribution of the $S_i^zS_{i+r}^z$ correlator exhibits striking differences between the MBL and Anderson insulator cases.

Figures

Equations

$$H=\sum _{i=0}^{L-1}\frac{J}{2}\left(S_i^{+}{S}_{i+1}^{-}+S_i^{-}{S}_{i+1}^{+}\right)+\Delta S_i^z{S}_{i+1}^z-h_i{S}_i^z.$$

$$⟨n|S_i^{+}{S}_{i+r}^{-}/2+h.c.|n⟩=⟨n|F_{i,i+r}|n⟩,$$

$$⟨n|S_i^z{S}_{i+r}^z{|n⟩}_c=⟨n|S_i^z{S}_{i+r}^z|n⟩-⟨n|S_i^z|n⟩⟨n|S_{i+r}^z|n⟩.$$

$$D_{\mathrm{KL}}\equiv -\int \mathrm{d}xP(x)\log \left(\frac{Q(x)}{P(x)}\right)$$

$$H/W=\frac{1}{W}\sum _i{S}_i^z{S}_{i+1}^z+{\tilde{h}}_i{S}_i^z+\frac{1}{W}\sum _i{F}_{i,i+1}=H_0+\frac{1}{W}V.$$

$$\begin{array}{rl}⟨\tilde{n}|F_{i,i+r}|\tilde{n}⟩& =\frac{1}{W}\sum _{k_0}⟨k_0|F_{i,i+r}|n_0⟩\frac{⟨n_0|V|k_0⟩}{E_{k_0}-E_{n_0}}\\ & +(k_0\leftrightarrow n_0).\end{array}$$

$$⟨\tilde{n}|F_{i,i+1}|\tilde{n}⟩=\frac{1}{2W\left(E_{n_0}-E_{k_0}\right)}=\frac{1}{2W\left({\tilde{h}}_i-{\tilde{h}}_{i+1}\right)}.$$

$$P(x)=\frac{4W|x|-1}{16W^2|x{|}^3},|x|\ge 1/4W.$$

$$\mathrm{D}=\text{span}\left(|n_0⟩,\{|i_0⟩:|i_0⟩=F_{i,i+1}|n_0⟩\text{and}{E}_{i_0}\approx E_{n_0}\}\right).$$

$$|\psi ⟩=a|\dots {\sigma }_{i-1}{\sigma }_i{\sigma }_{i+1}{\sigma }_{i+2}\dots ⟩+b|\dots {\tau }_{i-1}{\tau }_i{\tau }_{i+1}{\tau }_{i+2}\dots ⟩.$$

$$\begin{array}{rl}& 4⟨S_i^z{S}_{i+1}^z{⟩}_c=|a{|}^2{\sigma }_i{\sigma }_{i+1}-\left(|a{|}^2-1\right){\tau }_i{\tau }_{i+1}\\ & -\left[|a{|}^2{\sigma }_i-\left(|a{|}^2-1\right){\tau }_i\right]\left[|a{|}^2{\sigma }_{i+1}-\left(|a{|}^2-1\right){\tau }_{i+1}\right].\end{array}$$

$$\begin{gathered} {\sigma }_i{\sigma }_{i+1}{\tau }_i{\tau }_{i+1}\in \{0011,1100\}: \\ 4⟨S_i^z{S}_{i+1}^z{⟩}_c=1-(2|a{|}^2-1{)}^2>0. \end{gathered}$$

$$\begin{gathered} {\sigma }_i{\sigma }_{i+1}{\tau }_i{\tau }_{i+1}\in \{0110,1001\}: \\ 4⟨S_i^z{S}_{i+1}^z{⟩}_c=(2|a{|}^2-1{)}^2-1<0. \end{gathered}$$

$$H_0$$

$$V$$

$$H_{\mathrm{eff}}\equiv \mathbb{P}H\Omega =H_{\text{eff}}^{(0)}+H_{\text{eff}}^{(1)}+H_{\text{eff}}^{(2)}+\dots$$

$${\Omega }^{(0)}$$

$${\Omega }^{(1)}$$

$${\Omega }^{(2)}$$

$$⟨k|SA|m⟩\equiv \frac{⟨k|A|m⟩}{E_0^m-E_0^k}.$$

$$⟨m|H_{\text{eff}}^{(0)}|m^{\prime }⟩=⟨m|H|m^{\prime }⟩.$$

$$|{\Psi }_{\lambda }⟩=\sum _{|m⟩\in \text{D}}⟨m⟩v_{\lambda }^{\text{D}}|m⟩=\sum _{m\in \text{D}}{v}_m^{\lambda }|m⟩.$$

$$⟨{\Psi }_{\lambda }|O|{\Psi }_{\lambda }⟩=\sum _{m,m^{\prime }\in \text{D}}{v}_{m^{\prime }}^{\lambda ,\ast }{v}_m^{\lambda }⟨m^{\prime }|O|m⟩.$$

$$⟨m|H_{\mathrm{eff}}^{(1)}|m^{\prime }⟩=\sum _{k\in \overline{\text{D}}}\frac{⟨m|H|k⟩⟨k|V|m^{\prime }⟩}{E_0^{m^{\prime }}-E_0^k}+\dots$$

$$H_{\text{eff}}=H_{\text{eff}}^{(0)}+H_{\text{eff}}^{(1)}\in {\mathbb{C}}^{\text{dim}(\text{D})\times \text{dim}(\text{D})}.$$

$$\ket{\Psi_{\lambda}}=\ket{v_{\lambda}}+\sum_{\tiny\begin{array}[]{c}k\in\overline{\mathrm{D}}\\ m\in\mathrm{D}\end{array}}v^{\lambda}_{m}\frac{\langle k|V|m\rangle}{E_{0}^{m}-E_{0}^{k}}\ket{k}.$$

$$⟨{\Psi }_{\lambda }|O|{\Psi }_{\lambda }⟩=\sum _{m,m^{\prime }\in \mathrm{D}}{v}_m^{\lambda }{v}_{m^{\prime }}^{\lambda \ast }⟨m^{\prime }|O|m⟩$$

$$\displaystyle+\sum_{\tiny\begin{array}[]{c}k\in\overline{\mathrm{D}}\\ m,m^{\prime}\in\mathrm{D}\end{array}}v_{m}^{\lambda}v_{m^{\prime}}^{\lambda\hskip 1.42262pt*}\langle m^{\prime}|{O}|k\rangle\frac{\langle k|V|m\rangle}{E_{0}^{m}-E_{0}^{k}}$$

$$\displaystyle+\sum_{\tiny\begin{array}[]{c}k\in\overline{\mathrm{D}}\\ m,m^{\prime}\in\mathrm{D}\end{array}}v_{m}^{\lambda}v_{m^{\prime}}^{\lambda\hskip 1.42262pt*}\langle k|{O}|m\rangle\frac{\langle m^{\prime}|V|k\rangle}{E_{0}^{m^{\prime}}-E_{0}^{k}}$$

$$⟨m|H_{\mathrm{eff}}^{(2)}|m^{\prime }⟩=\sum _{k,k^{\prime }\in \overline{\text{D}}}⟨m|H|k⟩\frac{⟨k|V|k^{\prime }⟩⟨k^{\prime }|V|m^{\prime }⟩}{(E_0^{m^{\prime }}-E_0^k)(E_0^{m^{\prime }}-E_0^{k^{\prime }})}$$

$$\displaystyle-\sum_{\tiny\begin{array}[]{c}k\in\overline{\mathrm{D}}\\ l\in\mathrm{D}\end{array}}\langle m|H|k\rangle\frac{\langle k|V|l\rangle\langle l|V|m^{\prime}\rangle}{(E_{0}^{l}-E_{0}^{k})(E_{0}^{m^{\prime}}-E_{0}^{k})}.$$

$$H_{\text{eff}}=H_{\text{eff}}^{(0)}+H_{\text{eff}}^{(1)}+H_{\text{eff}}^{(2)}.$$

$$\begin{split}\Omega^{(2)}\ket{v_{\lambda}}&=\sum_{\tiny\begin{array}[]{c}k,k^{\prime}\in\overline{\mathrm{D}}\\ m\in\mathrm{D}\end{array}}\ket{k}\frac{\langle k|V|k^{\prime}\rangle\langle k^{\prime}|V|m\rangle}{(E_{0}^{m}-E_{0}^{k})(E_{0}^{m}-E_{0}^{k^{\prime}})}v^{\lambda}_{m}\\ &-\sum_{\tiny\begin{array}[]{c}k\in\overline{\mathrm{D}}\\ l,m\in\mathrm{D}\end{array}}\ket{k}\frac{\langle k|V|l\rangle\langle l|V|m\rangle}{(E_{0}^{l}-E_{0}^{k})(E_{0}^{m}-E_{0}^{k})}v^{\lambda}_{m}.\end{split}$$

Full text

I Introduction

Over the last decade, the Eigenstate Thermalization Hypothesis (ETH) Peres (1984); Deutsch (1991); Srednicki (1994, 1999); Rigol et al. (2008); D’Alessio et al. (2016); Borgonovi et al. (2016) has become the essential framework for reconciling quantum dynamics with statistical mechanics. In its simplest form, ETH posits that expectation values of local observables in energy eigenstates are smooth functions of the energy eigenvalue in the thermodynamic limit. This provides a mechanism for thermalization in isolated quantum systems. ETH can be understood in terms of random matrix theory: ergodic quantum systems are essentially well described by random matrix ensembles at least where local observables are concerned. This leads to the above smoothness condition, and also predicts the correct scaling of statistical deviations from it with system size Beugeling et al. (2014); D’Alessio et al. (2016); Khaymovich et al. (2019); Ikeda and Ueda (2015). In particular, the expectation values of local observables have gaussian distributions — the distribution shape is an important characteristic of ETH behavior Khaymovich et al. (2019).

It turns out that some disordered interacting systems can avoid thermalization if disorder is strong enough. Such a nonequilibrium phase of matter is called the Many Body Localized (MBL) phase Anderson (1958); Fleishman and Anderson (1980); Georgeot and Shepelyansky (1998); Basko et al. (2006); Gornyi et al. (2005); Oganesyan and Huse (2007); Žnidarič et al. (2008); Berkelbach and Reichman (2010); Pal and Huse (2010); Schreiber et al. (2015); Luitz et al. (2015); Nandkishore and Huse (2015); Altman and Vosk (2015); Imbrie (2016a); Abanin and Papić (2017); Alet and Laflorencie (2018); Pietracaprina and Laflorencie (2019). In this phase, transport is completely halted and the system becomes a perfect insulator. In particular ETH is not valid Pal and Huse (2010); Luitz et al. (2015). The current theoretical understanding of the MBL phase relies on the emergence of integrability via a complete set of local integrals of motion (LIOM) Serbyn et al. (2013); Huse et al. (2014); Imbrie (2016a); Ros et al. (2015); Imbrie (2016b). For instance, this theory accounts for the failure of thermalization, the area law of entanglement entropy in infinite temperature eigenstates Bauer and Nayak (2013) and the logarithmic growth of entanglement entropy after a quench Žnidarič et al. (2008); Bardarson et al. (2012).

Even though the existence of the MBL phase is by now well established in one dimensional systems, in both theory Basko et al. (2006); Imbrie (2016b, a) and experiments Schreiber et al. (2015); Smith et al. (2016), the nature of the localization-delocalization transition remains an active area of research. One outstanding question is the universality of anomalous thermalization Luitz and Bar Lev (2016); Roy et al. (2018), characterized by sub-diffusive transport, close to the localization transition coming from the ergodic side Bar Lev and Reichman (2014a); Agarwal et al. (2015); Potter et al. (2015); Vosk et al. (2015); Luitz et al. (2016); Lev et al. (2017); Bordia et al. (2017); Luitz and Lev (2017); Žnidarič and Ljubotina (2018); Kozarzewski et al. (2018); Schulz et al. (2018); Lezama et al. (2019). Moreover, there is evidence that distributions of diagonal matrix elements of the local (globally conserved) density develop heavy tails in this anomalous thermal phase Luitz (2016). It has been suggested that the latter is connected to the sub-diffusive transport and, in addition, could be described by a modified version of ETH Luitz and Bar Lev (2016); Roy et al. (2018). However it is not clear whether power law tails in the distribution of local operators are a general feature of the sub-diffusive regime.

In this work, we consider the probability distributions of local correlation functions in mid-spectrum energy eigenstates to determine their features in the ergodic as well as in the MBL phase. While the gaussian shape of these distributions is a central property of pure ETH, their behavior is equally important to characterize non-ergodic phases, in particular the MBL phase. Considering the Heisenberg model with random on-site fields, we present and analyze the distributions for two-point operators: spin-flip and $S_{i}^{z}S_{i+r}^{z}$ operators. Due to the U(1) symmetry of the XXZ model, there are no other non-vanishing two-point correlators. Furthermore, we carry out quantitative quasi-degenerate perturbation-theory calculations (around the limit of infinite disorder) to explain various features of the distributions in the MBL phase.

The energy eigenstate distributions of spin-flip and $S_{i}^{z}S_{i+r}^{z}$ correlators considered in this paper are gaussian for small disorder strengths but acquire significant weight in the tails already for intermediate disorder $W\approx 2<W_{c}\approx 3.7$, $W_{c}$ being the critical disorder strength to enter the MBL phase. Despite the heavy tails in the thermal regime, the variance of the distribution falls off with increasing system size for fixed $W$ up to the critical value $W_{c}$. Within the MBL regime, $W>W_{c}$, the variation of the distribution with increasing system size is negligible and the distribution has features not present in the thermal regime. In particular, the spin-flip operator distribution exhibits a sharp peak at zero with smaller satellite peaks on each side and further small peaks at the edge of the distribution, $\pm 1/2$. Perturbation theory captures the form of the large disorder distribution quantitatively.

Perturbative methods to describe localization-delocalization phenomena in condensed matter physics have a long history dating back to Anderson’s seminal work and continuing today to address questions relating to MBL Anderson (1958); Ros et al. (2015); Scardicchio and Thiery (2017); Basko et al. (2006); Gornyi et al. (2005); Imbrie (2016a). In the context of MBL, two of the main questions were to systematically construct the local integrals of motion that are thought to characterize the MBL phase and to estimate the transition point between MBL and the thermal phase. Both can be achieved by computing perturbative corrections to the mutually commuting occupation numbers at infinite disorder under the constraint that the corrections themselves continue to commute Ros et al. (2015). This has to be done to infinite order within some suitable approximation to capture possible delocalization. To make sense of the perturbation theory, as in the case of Anderson localization, there are resonances that lead to naive divergences coming from states close in energy that are mixed by hopping in the latter case and interactions in the former. In both cases, the divergences may be resolvable giving the perturbation theory a finite radius of convergence. Resolving these divergences amounts to diagonalizing the resonating configurations exactly.

The perturbation theory discussed in this paper is an expansion in the hopping part of the Hamiltonian around the infinite disorder limit. We carry out the expansion to low orders to be quantitative at large disorder for our finite system and to capture the main qualitative features for smaller disorder within the MBL regime. In the spirit of earlier works, we deal with resonances in the non-degenerate perturbation theory by diagonalizing exactly on the resonant subspaces.

In Section II we present the model and the local operators whose correlation functions we study. Section III focuses on the spin-flip operators across the whole phase diagram: first to nearest neighbor, then the further neighbor spin-flip operator distributions. The form of the spin-flip operator distributions in the MBL regime are rationalized within perturbation theory in Section III.3. We then turn to the $\langle S_{i}^{z}S_{i+r}^{z}\rangle$ correlators (Section IV) and the corresponding connected correlators (Section IV.2), in both cases showing the development of heavy tails at $W\approx 2.0$ and the evolution of these distributions into the MBL regime. We highlight the distinctive form of the distributions in the Anderson localized phase and the difference with the corresponding MBL distributions (Section IV.3). Finally, we compute the distribution using quantitative perturbation theory showing, once again, that it captures well the form of the distributions at strong disorder (Section IV.4).

II Background Material

II.1 Model

We study the canonical XXZ model with random fields $h_{i}$ along the $z$ direction,

[TABLE]

This model – which is widely studied in the context of MBL Žnidarič et al. (2008); Pal and Huse (2010); Berkelbach and Reichman (2010); Bardarson et al. (2012); Luca and Scardicchio (2013); Serbyn et al. (2013); Bauer and Nayak (2013); Nanduri et al. (2014); Bar Lev and Reichman (2014b); Luitz et al. (2015); Nandkishore and Huse (2015); Agarwal et al. (2015); Bar Lev et al. (2015); Bera et al. (2015); Torres-Herrera and Santos (2015); Luitz et al. (2016); Serbyn and Moore (2016); Singh et al. (2016); Enss et al. (2017); De Tomasi et al. (2017); Bera et al. (2017); Lezama et al. (2017); Alet and Laflorencie (2018); Herviou et al. (2019); Sierant and Zakrzewski (2019); Šuntajs et al. (2019); Serbyn et al. (2017); Maksymov et al. (2019) – can be mapped to a spinless fermion model with nearest neighbor hopping $J/2$, interaction term $\Delta$ and on-site potential $h_{i}$. In this paper, periodic boundary conditions are set, $J=1$ is fixed throughout the paper and the fields $h_{i}$ are distributed uniformly in $[-W,W]$ with disorder strength $W$. We focus mainly on the isotropic point $\Delta=1$ (interacting spinless fermions) of the parameter space. However in Section IV.3 we compare also to results for various $\Delta$, including the point $\Delta=0$ (free spinless fermions). The operators $S^{\alpha}_{i}=\sigma^{\alpha}_{i}/2$ are spin 1/2 operators, with $\alpha=0,x,y,z$ and $i$ is the site index.

The total magnetization $M=\sum_{i=0}^{L-1}S^{z}_{i}$ along the $z$ direction is conserved. We therefore focus on the largest magnetization sector $M=0$ for even system sizes $L$, corresponding to the Hilbert space dimension $\text{dim}(\mathcal{H})=\text{binom}(L,\lfloor{L/2\rfloor})$. For each disorder realization $\{h_{0},\dots,h_{L-1}\}$, we obtain $\gtrsim 50$ eigenstates closest to the energy target $(E_{\text{max}}+E_{\text{min}})/2$ ($E_{\text{min}}$ being the ground state energy and $E_{\text{max}}$ the highest energy of the sample) using a state-of-the-art shift-invert code Luitz et al. (2015); Pietracaprina et al. (2018). We consider the probability distributions of various eigenstate expectation values of local operators, i.e. the diagonal matrix elements of these operators in the eigenbasis of the Hamiltonian. Our results are histograms over at least $10^{3}$ disorder realizations for each system size $L$ and disorder strength $W$, we also calculate the correlators for all sites $i\in[0,L-1]$ to improve the statistics, since the average over disorder is translation invariant. The mid-spectrum states of this model are known to exhibit two dynamical phases Pal and Huse (2010); Luitz et al. (2015): at low disorder ($W\lesssim 3.7$) they obey the ETH, while at strong disorder ($W\gtrsim 3.7$) all eigenstates are many-body localized (MBL).

II.2 Operators

In previous works in the context of many-body localization and the MBL transition, the distributions of local operators were considered, mostly focussing on distributions of diagonal or off-diagonal matrix elements of simple local observables such as the local magnetization (or number density in the language of spinless fermions) $\bra{n}S_{i}^{z}\ket{n}$, where $\ket{n}$ is a central eigenstate of the Hamiltonian Luitz (2016); Luitz and Bar Lev (2016). In this work, we consider more complicated operators given by two point correlation functions. First, we consider the correlators

[TABLE]

i.e., the matrix elements of spin-flip operators $F_{i,i+r}$. We also consider diagonal two-point correlators, namely $\bra{n}S^{z}_{i}S^{z}_{i+r}\ket{n}$ and its ‘connected’ version

[TABLE]

For $r=1$, the first expression above corresponds to the kinetic energy density, while the second expression is the interaction energy density in the language of spinless fermions. The connected correlator $\bra{n}S^{z}_{i}S^{z}_{i+r}\ket{n}_{c}$ was previously considered in Ref. Pal and Huse, 2010.

III Eigenstate expectation values of $S^{+}_{i}S^{-}_{i+r}+h.c.$

In Fig. 1 we show the probability distribution of eigenstate expectation values of $F_{i,i+r}=S^{+}_{i}S^{-}_{i+r}/2+h.c.$ for a system of size $L=20$ and different disorder strengths $W$.

III.1 Nearest neighbor flip

We start by considering the special case $r=1$, where the operator $F_{i,i+1}$ corresponds to the kinetic energy per bond. In the thermal phase at weak disorder, we expect this operator to follow ETH and be distributed according to a gaussian distribution, which is true to very good precision.

We observe that at weak disorder ($W\lesssim 2$) and $r=1$, the mean of the distribution is slightly negative, a result of the fact that the eigenstates of the Hamiltonian we consider are in the center of the spectrum, and correspond to high but finite temperature due to the asymmetry of the density of states (cf. Appendix B for an analysis of the energy dependence). States corresponding to strictly infinite temperature correspond to energies given by $\mathrm{Tr}(H_{M=0})/\text{dim}(\mathcal{H}_{M=0})=-\frac{L}{4L-4}$, where ${H}_{M}$ is the Hamiltonian matrix in the zero magnetization sector. Such states have a zero mean for traceless operators like $F_{i,i+r}$. Zero mean distributions are recovered at intermediate disorder where the asymmetry of the spectrum is less pronounced and the energy of the eigenstates we consider is indeed close to $-\frac{1}{4}$ for large $L$.

At intermediate disorder $W\approx 2$, we observe the development of heavy tails in the distribution, very similar to the situation for the distribution of $\bra{n}S_{i}^{z}\ket{n}$ studied in Ref. Luitz, 2016, confirming that the presence of such tails appears to be a generic feature at intermediate disorder in the thermal phase. We note that heavy tails are also observed in the $S_{i}^{z}S_{i+r}^{z}$ correlation function studied in Sec. IV.

At strong disorder $W>W_{c}$ in the MBL phase, we find a strikingly different distribution of the spin-flip operator expectation values $F_{i,i+1}$; it features a pronounced central peak at zero, accompanied by two minima adjacent to it, which are framed by two satellite peaks, before the probability density $p(\bra{n}F_{i,i+1}\ket{n})$ decays towards the edges of its domain $\left[-\frac{1}{2},\frac{1}{2}\right]$. We have found that this intriguing shape persists at strong disorder and can be explained using perturbation theory, a discussion of which we postpone to the end of this section.

In Fig. 2, we analyze the system size dependence of the probability density of the nearest neighbor flip operator $F_{i,i+1}$ over the whole range of disorder strengths, comparing distributions for sizes $L=12,14,16,18,20$. At the weakest disorder $W=0.4$, we find gaussian probability distributions, with the variance decreasing exponentially in system size $L$, as expected from ETH (cf. Fig. 3). At intermediate disorder $W=2.0$, the distribution is no longer gaussian, but the variance still decreases exponentially with size. It appears that the heavy tails, deviating from the gaussian shape, persist even at large system size, following the same phenomenology observed for the $S_{i}^{z}$ operator in Ref. Luitz, 2016; Luitz and Bar Lev, 2016; Roy et al., 2018.

To quantify departures from gaussianity, we compute the excess kurtosis $\kappa=(\mu_{4}/\sigma^{4})-3$, ($\mu_{4}$ being the $4$th central moment of the distribution) and the Kullback-Leibler divergence defined by

[TABLE]

where $Q(x)$ is the reference gaussian distribution and $P(x)$ is the computed distribution of the correlator, where the gaussian is defined by the mean and variance of $P(x)$. Results are shown in Fig. 4. Both quantities indicate that the distribution is quantitatively gaussian for $W\lesssim 1.5$ and that they become strikingly less gaussian with a peak at about $W=2$ that increases with system size. Beyond the peak for larger disorder both measures increase smoothly with little system size dependence, indicating strongly nongaussian distributions in the MBL phase.

In the MBL phase at strong disorder, there is no discernible system size dependence of the distribution (Figs. 2, lower panels, and 3), showing a pronounced maximum at zero, framed by two symmetric satellite peaks, which seem to get closer to each other at stronger disorder.

III.2 Long distance flip

The flip operator of distant spins $F_{i,i+r}$, with $r>1$ is not a term of the Hamiltonian and could therefore behave differently. We have verified that this is so by examining the energy dependence of the mean of the distribution which is constant over a large range of energies for $r>1$, linear for $r=1$, cf. Appendix B. For this reason, the mean of the $r>1$ distribution is close to zero at weak disorder. At intermediate disorder, the distribution again shows heavy tails, and, in general, the variance decreases with longer distance between the operators, which we attribute to decreasing long distance correlations. Most interestingly, at strong disorder in the MBL phase and at long distance, the peculiar satellite peaks of the distribution at $r=1$ disappear, leading to simple, yet heavy tails. Additionally, we see that the standard deviation of both correlation functions at larger distances decreases as function of disorder (Fig. 3) and stay constant at $r=1$. In the limit $W\rightarrow\infty$ the spins are uncorrelated so both standard deviations will go to zero. In this range of disorder, the localization length is big enough for allowing correlations at $r=1$, hence we expect the standard deviation to start decreasing only at large enough disorder.

The absence of satellite peaks for $r>1$, as well as most of the other features in this and the preceding subsection, can be understood through perturbation theory in $1/W$, as we describe in the next subsection.

III.3 Perturbation theory analysis

We have postponed the discussion of the peculiar features of the distribution of the nearest neighbor flip operator $F_{i,i+r}$ in the MBL phase – a topic to which we now turn.

In Fig. 5 we have a closer look at its distribution for different (strong) disorder strengths. While the qualitative features (central and satellite peaks) are independent of disorder and apparently characteristic of MBL, there is a quantitative evolution: the satellite peaks become sharper and move towards zero as the disorder $W$ is increased (inset in Fig. 5). Furthermore, at very strong disorder $W>10$, additional peaks at $-\frac{1}{2}$ and $\frac{1}{2}$ develop, which are not present at weaker disorder $W\lesssim 6$ (cf. Fig. 2).

As a first step towards a more quantitative analysis, we consider the drift of the position of the satellite peaks as a function of disorder. The lower left panel of Fig. 5 shows the estimated peak positions, which are consistent with a $1/W$ dependence, suggesting a perturbative analysis.

At very strong disorder, it is natural to treat the kinetic term of the Hamiltonian as a perturbation of order $1/W$. Noting that the eigenstates of $H/W$ are equal to those of $H$, we cast the Hamiltonian in the form

[TABLE]

The scaled fields, $\tilde{h}_{i}$, are now distributed uniformly in a fixed range $[-1,1]$. The eigenstates of the unperturbed Hamiltonian $H_{0}$ are product states and eigenstates of all $S_{i}^{z}$ operators and can therefore be enumerated by their eigenvalues. The eigenenergies of $H_{0}$ for each eigenstate can be easily calculated using these quantum numbers.

Naive perturbation theory produces divergences when the spacing between unperturbed energy levels goes to zero. Such divergences – dubbed resonances – are unphysical and are resolved by admixing clusters of nearly degenerate states. Resonances are of great importance in disordered systems and become increasingly so as the system size increases. In order to incorporate the effect of resonances from the large disorder limit, we carry out a mixed degenerate and non-degenerate perturbation theory for the operator $F_{i,i+r}$. Details of the perturbation theory are given in Appendix A. In addition to the rather general discussion given in the appendix, we note here various peculiarities of the perturbative calculation of $\bra{\tilde{n}}F_{i,i+r}\ket{\tilde{n}}$ which simplify our task. In order to obtain a matrix element $\bra{\tilde{n}}F_{i,i+r}\ket{\tilde{n}}$ of an eigenstate $\ket{\tilde{n}}$ of the perturbed Hamiltonian in perturbation theory, we start with an eigenstate $\ket{n_{0}}$ of $H_{0}$. The matrix element $\bra{n_{0}}F_{i,i+r}\ket{n_{0}}$ is the zeroth order contribution and is identical to zero because $F_{i,i+r}$ is off-diagonal in the $z$ basis. It therefore contributes to the prominent peak of the distribution of this matrix element at zero. More precisely, for $r=1$, states with $|\dots 00\dots\rangle$ or $|\dots 11\dots\rangle$ on the sites $i$ and $i+1$ yield a zero contribution at zeroth and first order in perturbation theory. This accounts for half of the states so we expect the fraction of such states to tend to $1/2$ as $W$ increases and this is indeed what is found (Fig. 5 lower right panel). For $r>1$, one must go to higher order in perturbation to obtain any non-vanishing contribution so the central peak is significantly higher. To understand the satellite peaks, we have to go to first order in perturbation theory (cf. e.g. Fig. 5).

The eigenstate $\ket{n_{0}}$ is connected to a set of other eigenstates $\{\ket{k_{0}}\}$ of $H_{0}$ by the perturbation $V$. By this, we mean that the $\bra{n_{0}}V\ket{k_{0}}\neq 0$ for all $\ket{k_{0}}$ in this set, while matrix elements of $V$ with all other states vanish. Let us first deal with the case in which all energies $E_{k_{0}}$ are sufficiently different from $E_{n_{0}}$, such that in nondegenerate perturbation theory the denominators $1/(E_{k_{0}}-E_{n_{0}})$ do not diverge. In this case, we obtain for the matrix element $\bra{\tilde{n}}F_{i,i+r}\ket{\tilde{n}}$ up to first order in $1/W$:

[TABLE]

From this, we can now understand several features of the distribution: if $\ket{n_{0}}$ does not have eigenvalues of $S_{j}^{z}$ with opposite sign on sites $j=i$ and $j=i+r$, the matrix element $\bra{\tilde{n}}F_{i,i+r}\ket{\tilde{n}}$ vanishes. This implies that due to the incompatibility of $V$ and $F_{i,i+r}$ for $r>1$, to first order in $1/W$ all matrix elements vanish, which explains the different behavior of $F_{i,i+1}$ and $F_{i,i+r},r\geq 2$. If the spins on sites $i$, $i+1$ have opposite $S_{j}^{z}$ eigenvalues (i.e. they are “flippable”), there is only one nonvanishing term in the sum, in Eq. (4). The matrix element $\bra{n_{0}}F_{i,i+1}\ket{k_{0}}=1/2$ in this case giving

[TABLE]

Since the on-site fields have a uniform distribution bounded by $\tilde{h}_{i}\in[-1,1]$, the expression in Eq. 5 can be computed. We first note that the lower bound on the matrix element is $1/4W$ as the maximum difference in the fields is $\pm 2$. Now, rewriting $x=\bra{\tilde{n}}F_{i,i+1}\ket{\tilde{n}}$ as a random variable that takes values in the range $[-\infty,-1/4W]\cap[1/4W,\infty]$, its probability distribution is:

[TABLE]

The maxima of $P(x)$ are located at $x=\pm 3/8W$ To summarize, first order nondegenerate perturbation theory explains the presence and weight of the central peak at zero, the presence and location of the satellite peaks at $O(1/W)$, as well as the local minima separating these peaks. The satellite peak stems therefore from the admixture of states which change their energy maximally upon flipping of two neighboring spins. Fig. 6 shows plots of the exact result (black) together with the distribution Eq. (6) (red dashed) showing that the formula captures the exact distribution very well (the delta peak at zero with weight $\frac{1}{2}$ is not shown in Fig. 6). To examine the agreement in more detail, Fig. 5 shows the close correspondence between the analytical calculation of the satellite peak location and the exact result at least for larger values of disorder. We notice that the perturbation theory produces a higher central peak. This is caused by the missing weight around the central maximum (see inset Fig. 5) due to the lower bound in magnitude of the matrix elements $\bra{\tilde{n}}F_{i,i+1}\ket{\tilde{n}}$ up to first order (Eq. 5).

The distribution, Eq. (6), does not reproduce the small peaks at the edge of the domain of the distribution, close to $\pm\frac{1}{2}$. To understand the origin of these peaks, we come back to the consideration of the case that the eigenenergy of the state with flipped spins $E_{k_{0}}$ is close to the energy of $E_{n_{0}}$, in which case we have a “resonance” and nondegenerate perturbation theory breaks down. In this case, we have to use quasi degenerate perturbation theory and include $\ket{n_{0}}$ and its flipped partner $\ket{k_{0}}$ in the model space of quasi-degenerate states. Due to the constraint by the matrix elements of $F_{i,i+1}$ in the model space, this is the only state which contributes to the model space. The mixing of these two states leads to the emergence of the peaks at $\pm\frac{1}{2}$ of the distribution. To see this, we consider the zeroth order mixing of quasi-degenerate states through a single spin flip term in the Hamiltonian. This generates pairs of admixed states of the form $\alpha|\ldots 10\ldots\rangle+\beta|\ldots 01\ldots\rangle$. The flip-flop operator expectation value is ${\rm Re}(\alpha\beta^{\star})$.Since the perturbation maximally mixes these quasidegenerate states $\alpha=\pm\beta$ and this accounts for the $\pm 1/2$ peaks.

Our numerical treatment of the exact mixing by quasi-degenerate perturbation theory up to second order in $1/W$ captures also corrections to these features quantitatively and we show the full distributions obtained from it as colored solid histograms in Fig. 6. The perturbation theory can be carried out for much larger system sizes than treatable in shift-invert diagonalization and show no visible system size dependence at strong disorder as shown in Fig. 6. We conclude that the parts of the distribution of $\bra{n}F_{i,i+1}\ket{n}$ close to zero are due to off-resonant mixing of flippable and not flippable states, while the edges of the distribution close to $\pm\frac{1}{2}$ reveal the effect of resonances. It should be noted that we do not compute $\bra{\tilde{n}}F_{i,i+r}\ket{\tilde{n}}$ at distances $r>1$ because low order contributions are trivial and higher orders in perturbation theory make the numerical implementation hard to deal with. Considering this, we restrict our perturbation theory computations to operators with $r=1$.

IV Eigenstate Expectation Values of $S^{z}_{i}S^{z}_{i+r}$

IV.1 $\left\langle S^{z}_{i}S^{z}_{i+r}\right\rangle$ Correlators

We now turn to the $S^{z}_{i}$ correlation function. Fig. 7 shows the probability distribution of energy eigenstate expectation values of $S^{z}_{i}S^{z}_{i+r}$ for a system of size $L=20$, $r=1,2,3,4$ and for different disorder strengths $W$. For weak disorder $W\lesssim 1.2$, the distributions are gaussian in accordance with ETH and the variance of the distribution increases with disorder strength. As remarked in Section III.1 and, similarly to the spin-flip correlators studied there, heavy tails are apparent for disorder strength $W=2$ in the thermal regime. Again, similarly to the spin-flip correlators, the gaussian mean is displaced from zero and the reason for this displacement is the same as in that case (cf. Appendix B). As one expects in the ETH regime, the variance of the distribution falls off inversely in the Hilbert space dimension (exponential in $L$), which is visible for the case of the connected correlator in Fig. 3 by the equidistant spacing of the standard deviations for different system sizes on the semilogarithmic scale.

For strong disorder, deep in the MBL regime, the distribution is qualitatively different. The central peak is still present but is obscured by a very broad distribution that extends out to the tails where there are more pronounced peaks. There are again negligible differences between the distributions for different $L$ within the MBL regime. The presence of the outer peaks is simply explained from the strong disorder limit where eigenstates of the Hamiltonian are also eigenstates of the local $S^{z}_{i}S^{z}_{i+r}$ operators with eigenvalue $\pm 1/4$. The fact that the main new feature of the distribution appears in the large $W$ limit suggests that perturbation theory might be as successful as it was for the spin-flip correlators. We address this question in Section IV.4.

In order to remove the trivial contribution to the correlation function coming from $\bra{n}S_{i}^{z}\ket{n}$ expectation values, we discuss, in the following section, the connected correlation function.

IV.2 Connected Correlators

Fig. 8 shows probability distributions of the connected correlation function for $L=20$ and for different values of $W$ and Fig. 9 shows distributions for different system sizes. For small $W$, the expectation is that ETH is obeyed and the figures demonstrate that, at least for $W\lesssim 1.2$, the distributions are gaussian (Fig. 8) while the finite size scaling is consistent with random matrix theory (cf. Fig 3). In the ETH regime, there is little variation in the distributions for different $r$ for a given system size – one merely observes that the mean of the distribution shifts towards zero from $r=1$ to $r>1$ as discussed above due to the different energy dependence of different $r$ operators (cf. Appendix B).

For larger values of $W$, a sharp peak forms at zero and persists deep into the MBL phase while the distributions further depart from gaussianity by acquiring a distinctive asymmetry with higher weight for negative values of the correlator. For $r=1$, the left-hand-side of the distribution acquires a shoulder down to $-1/4$ while the positive side tapers off towards $+1/4$. For larger $r$ the shoulder is rounded on the left side, so that the asymmetry is less pronounced. Our analysis in Appendix C confirms heavy tails on either side at strong disorder.

In common with other distributions of matrix elements of local operators there is little apparent variation between different system sizes in the MBL regime. In contrast, within the ETH regime for significant values of disorder as exemplified by the $W=2$ data, the central width of the distribution narrows for larger system sizes while weight at the tails remains.

IV.3 Anderson Insulator vs MBL

To understand the asymmetry of the distribution of the connected correlator $\langle S_{i}^{z}S_{i+r}^{z}\rangle_{c}$ for small distances $r$, it is useful to compare to the noninteracting limit. In Eq. 1, $\Delta=0$ corresponds to the case of an Anderson insulator of noninteracting spinless fermions. Fig. 10 shows the connected $S_{i}^{z}S_{i+r}^{z}$ correlator for $\Delta=0.0,0.1,0.2,0.5,1.0$ and for $W=6$ and $r=1,2,3,4$. We see that for distance $r=1$ the distribution of negative correlations has little sensitivity to the value of $\Delta$ while the weight of positive correlations is exactly zero for the Anderson insulator, giving a clear signature where MBL differs from the non-interacting case albeit one where the asymmetry between positive and negative weights persists to $\Delta=1$.

We can understand the vanishing positive weight for the Anderson insulating case for arbitrary $r$ through a straightforward application of Wick’s theorem since the Anderson case $\Delta=0$ is a free spinless fermion model with a gaussian action. In fermionic language, the $\hat{S}^{z}_{i}\hat{S}^{z}_{i+r}$ operator takes the form $\left(c_{i}^{\dagger}c_{i}-\frac{1}{2}\right)\left(c_{i+r}^{\dagger}c_{i+r}-\frac{1}{2}\right)$. Using Wick’s theorem for $\Delta=0$, we obtain in any eigenstate of the Hamiltonian: $\langle c_{i}^{\dagger}c_{i}c_{j}^{\dagger}c_{j}\rangle=\langle c_{i}^{\dagger}c_{i}\rangle\langle c_{i}^{\dagger}c_{i}\rangle-\langle c_{i}^{\dagger}c_{j}\rangle\langle c_{j}^{\dagger}c_{i}\rangle$. It follows that the connected correlator is $\langle n|\hat{S}^{z}_{i}\hat{S}^{z}_{i+r}|n\rangle_{c}=-\left|\langle n|c_{i}^{\dagger}c_{i+r}|n\rangle\right|^{2}$, which is necessarily $\leq 0$. This leads to the extreme asymmetry of the connected correlator distribution at $\Delta=0$.

We see in Fig. 10 that the asymmetry decreases as $r$ increases. This can be understood perturbatively, as discussed in the next section.

IV.4 Perturbation Theory

As in the case of the spin-flip correlator, the distributions of the $S^{z}$ correlator are well reproduced at large $W$ by the perturbation theory discussed in Sec. III.3 and Appendix A. This is demonstrated in Fig. 11 for the distributions of connected correlators with $r=1$. Perturbation theory also provides qualitative physical explanations for the features reported in the last three subsections, as we elaborate below.

In the infinite disorder limit $W\to\infty$, the correlators $\langle S^{z}_{i}S^{z}_{i+r}\rangle$ and $\langle S^{z}_{i}\rangle$ are respectively $\pm 1/4$ and $\pm 1/2$ and the connected correlator simply vanishes, contributing to the sharp peak of the distribution at zero. Indeed, the numerical data in Fig. 7, for large values of $W$, shows peaks at the extreme values of the matrix elements. Similarly the connected correlator distribution in Fig. 8 exhibits a smooth decay of the central peak at zero towards finite values of the connected correlator.

Inspecting the general expressions for matrix elements in mixed degenerate and nondegenerate perturbation theory presented in appendix A, we observe that the zeroth order term coming from mixing within the model space should account for some degree of broadening of the peaks. The first order terms in $1/W$ trivially vanish, because the $S_{i}^{z}S_{i+r}^{z}$ operator does not connect states in the model space to those outside it, and this is why we proceed to compute the second order contribution. Fig. 11 shows that perturbation theory to second order (yellow) compares very well to the exact finite size results for $L=18$ and large disorder $W=12,\dots 16$ whereas the zeroth order results (orange) fail to capture the smooth falloff of the distribution to the left of the central peak. We also show results for larger system sizes in Fig. 11, which are not reachable otherwise and show that at large disorder the distributions are essentially converged.

We show below that, while zeroth order (i.e., quasi-degenerate) perturbation theory does not fully account for the distribution shape of the connected correlators, it does capture the asymmetry.

Zeroth order perturbation theory mixes quasi-degenerate eigenstates of ${H}_{0}$ connected by ${V}=\sum_{i}F_{i,i+1}$ – in other words states connected by flippable spins $|\ldots 01\ldots\rangle$ or $|\ldots 10\ldots\rangle$ where the ellipses denote some spin configurations. This means, starting from an eigenstate $\ket{n_{0}}$ of $H_{0}$, we can expect mixing with the states $\{F_{i,i+1}\ket{n_{0}}\}$, which are quasi-degenerate 111In practice, we set a cutoff value for the energy difference which we still treat as quasi degenerate of the order of $1/W$. We checked that choosing $2/W$ or larger values does not significantly change the results with $\ket{n_{0}}$. This means that the model space is then spanned by

[TABLE]

Let us now try to understand why the distribution of the connected correlator is asymmetric, starting from the case $r=1$. For simplicity, we consider the case of a two dimensional model space, yielding a state (with $|b|^{2}=1-|a|^{2}$) of the form:

[TABLE]

The connected correlator is then given by

[TABLE]

Inspecting this expression shows that most combinations of spin configurations $\sigma_{i}$,$\sigma_{i+1}$,$\tau_{i}$,$\tau_{i+1}$ yield vanishing connected correlators and these contribute to the central peak. The spin configurations on $i$, $i+1$ that yield non-vanishing contributions are:

[TABLE]

[TABLE]

This means we obtain two cases for positive correlators and two for negative correlators. Evidently the case with a flippable pair $\sigma_{i}\neq\sigma_{i+1}$ (yielding a negative correlator) appears at first order in $V$, since the two states are directly connected through $V$ and are included in the model space if they are quasi-degenerate, independent of the state of the neighboring spins $\sigma_{i-1}$ and $\sigma_{i+2}$. The case of an aligned pair $\sigma_{i}=\sigma_{i+1}$ (yielding a positive correlator) is connected to its flipped partner state $\tau_{i}=-\sigma_{i}$ and $\tau_{i+1}=-\sigma_{i+1}$ only in second order of $V$, including a constraint on the neighboring spins $\sigma_{i-1}=-\sigma_{i}$ and $\sigma_{i+2}=-\sigma_{i+1}$. We note that in addition, an intermediate state with one spin flip has to be quasi-degenerate, which is an additional constraint. For simplicity we have left this state out of the discussion.

From these arguments, we conclude that the case of admixed states which yield negative correlations is much more probable than the case yielding positive correlations due to their appearance at different orders in $V$ and, additionally, owing to constraints which reduce the number of possibilities giving obtaining positive correlations.

We now understand that for the case $r=1$, negative correlations are more probable than positive ones to zeroth order in degenerate perturbation theory for two reasons: negative correlations need only one order in $V$, while the appearance of positive weight requires two applications of $V$ and only a specific set of spin configurations can lead to positive correlations, thus reducing their likelihood.

The same set of arguments can now be generalized to the case $r=2$. We see that in this case, we always need to apply $V$ twice to get nonzero (both positive and negative) correlations, however there are more possibilities of having a flippable pair $i,i+2$ (necessary for negative correlations), compared to the possibilities of getting a mixture of $\ket{\dots 0\text{{x}}0\dots}$ and $\ket{\dots 1\text{{x}}1\dots}$ (necessary for positive correlations), since in this case the flippability depends on the state x of the middle spin $i+1$. Therefore, also in the case $r=2$, the distribution of the connected correlator is skewed towards negative correlations. In the case of longer distances $r>2$, these constraints become increasingly weak (while requiring an order of $V^{r}$ to get nonzero correlators), leading to more and more symmetric distributions.

V Conclusions

We have presented the exact energy eigenstate distributions of spin-flip and $S_{i}^{z}S_{i+r}^{z}$ correlators in the disordered XXZ chain across the many-body localization transition. While – at very weak disorder – we find gaussian distributions to very high precision, the distributions depart from gaussianity at intermediate disorder – still well inside the thermal regime – through the appearance of heavy tails that persist into the MBL regime. The presence of these tails correlates to the appearance of sub-diffusive behavior in transport properties observed in previous studies Luitz (2016); Luitz and Bar Lev (2016); Roy et al. (2018). In the entire thermal regime, the variance of the correlator distributions falls off with increasing Hilbert space dimension as one should expect for operators obeying ETH but significant weight remains in the tails of the distribution and measures of departures from gaussianity including the Kullback-Leibler divergence and the kurtosis show a peak for $W<W_{c}$ that sharpens with system size. The system size dependence of local operator distributions is negligible inside the MBL regime where ETH fails.

For large disorder, we have carefully investigated the distinctive forms of the correlator distributions, unraveling various features of the distributions. We find that strong disorder perturbation theory can reproduce the full distributions in the MBL phase. We note that our semianalytical perturbation theory scheme should be applicable to other models and could provide information about the effect of resonances in different systems.

For the spin-flip correlator, the distributions are highly structured with a central peak at zero, a pair of neighboring satellite peaks with disorder strength dependent positions at $\pm 3/8W$ and further maxima at the edge of the distribution at $\pm 1/2$. All these features are perfectly captured by a quantitative strong disorder perturbation theory that also gives insight into their origins. In particular, (i) the central peak comes from eigenstates where the eigenstate carries no pairs of spins that are flippable by the spin-flip operator, this accounts for $1/2$ of all states at strong disorder (ii) the satellite peaks at $\pm 3/8W$ arise from flipped pairs of spins that are maximally pinned by the random field and therefore maximally off resonant (iii) the $\pm 1/2$ peaks arise from resonances - strongly admixed quasi-degenerate states. These extremal peaks can only be captured by quasi-degenerate perturbation theory. Overall, mixed quasi-degenerate and degenerate perturbation theory unifies all contributions and yields an unbiased result, matching the full exact distribution almost perfectly.

The $S_{i}^{z}S_{j}^{z}$ correlator distribution is more complicated to analyze since we have to go to second order in $1/W$ in our perturbative treatment. Our analysis reveals that for short distances $r=|i-j|$, the correlator is predominantly negative in the MBL phase, since eigenstates are biased to contain mixtures of flippable neighboring pairs. This leads to distributions skewed towards negative weights, most strongly so for noninteracting Anderson Insulators, where no weight for positive correlators is present due to Wick’s theorem. Therefore, the $S_{i}^{z}S_{j}^{z}$ correlator distribution reveals a strikingly different behavior generated by interactions in the MBL case compared to the noninteracting model.

Acknowledgements.

We would like to thank Jeff Rau for useful discussions. We furthermore acknowledge PRACE for awarding access to HLRS’s Hazel Hen computer based in Stuttgart, Germany under grant number 2016153659. Our code is based on the PETSCBalay et al. (2014a, 1997, b) and SLEPcHernandez et al. (2005) libraries and uses MUMPSAmestoy et al. (2001, 2006).

Appendix A Perturbation Theory

In this section, we give details of the Rayleigh-Schrödinger perturbation theory that is used to obtain the leading contributions in $1/W$ from the transverse exchange term $V$ of the Hamiltonian to the energy eigenstate expectation values of $S_{i}^{z}S_{i+r}^{z}$ and spin-flip correlators (see main text). We intend for this appendix and the accompanying part of the main text to be self-contained on the method but we refer the interested reader to Lindgren Lindgren (1974) for further details.

The Hamiltonian from Eq. (1) is rewritten $H/W=H_{0}+(1/W)V$ where

[TABLE]

and $\tilde{h}_{i}\in[-1,1]$ are normalized normal fields.

We wish to compute matrix elements of operator ${O}$ in the eigenstate basis where the eigenstates are computed to some order in perturbation theory in $1/W$. The eigenstates of interest here are taken to be highly excited states from the middle of the spectrum. We start from the eigenstates $\ket{n}$ of $H_{0}$, which are product states in the $S^{z}$ basis, their eigenenergies $E_{0}^{n}$ is trivially obtained from $H_{0}$. To compute the corrections to the eigenvectors $\ket{n}$ perturbatively in $1/W$ we should bear in mind that the energies of the product states are essentially randomly distributed and will mix strongly under the perturbation $V$ if the states are quasidegenerate. Loosely speaking, we would like to organize the perturbation theory so that states coupled by the pairwise spinflips in the perturbation $V$ that are separated by an energy $\Delta E>1/W$ are treated via non-degenerate perturbation theory while clusters of quasi-degenerate states (with mutual energy differences $\lesssim 1/W$) are treated via degenerate perturbation theory. We therefore use the formalism of a mixed quasi-degenerate and nondegenerate perturbation theory as described in Ref. Lindgren, 1974.

We discuss how to organize this calculation up to order $1/W^{2}$, keeping in mind our main goal: the computation of energy eigenstate expectation values of local operators.

Let us focus on a single $H_{0}$ eigenstate $\ket{n}$ in one disorder realization that has eigenvalue $E_{0}^{n}$. We compute all other states $\{\ket{m}\}$ connected to it by $V$ and $V^{2}$, i.e. $\bra{n}V\ket{m}\neq 0$ or $\bra{n}V^{2}\ket{m}\neq 0$. Now, we define the model space, D, to consist of all states in this set ${\ket{m}}$ (including the “parent state” $\ket{n}$), and with $|E_{k}-E_{n}|<\alpha/W$. We have experimented with the threshold $\alpha$ to include states in the model space and the results presented in the paper are essentially identical for $\alpha\in[1,3]$. Let $\mathbb{P}$ be the projector onto the model space D and $\mathbb{Q}=1-\mathbb{P}$ the projector on its complement $\overline{\text{D}}$.

Suppose $\ket{\Psi_{n}}$ is an exact eigenstate of the full Hamiltonian that has a nonvanishing projector into the model space $\ket{\Psi^{D}_{\lambda}}=\mathbb{P}\ket{\Psi_{\lambda}}$. We introduce the inverse wave operator $\Omega=\Omega^{(0)}+\Omega^{(1)}+\Omega^{(2)}+\ldots$ that is expanded in powers of $V$. The action of $\Omega$ is to rotate a state in the model space into the full space eigenfunction $\ket{\Psi_{\lambda}}=\Omega\ket{\Psi^{0}_{\lambda}}$. So now if we compute an eigenstate $\ket{v^{D}_{\lambda}}$ of the effective Hamiltonian

[TABLE]

this is nothing but the projection of the exact wavefunction onto the model space. The eigenfunctions of $H_{\rm eff}$ are generally non-orthogonal and the effective Hamiltonian non-Hermitian (albeit with real eigenvalues). To obtain the eigenstates of $H$ we lift them out of the model space by acting with $\Omega$: $\ket{\Psi_{\lambda}}=\Omega\ket{v^{D}_{\lambda}}$. Note that $\ket{v_{\lambda}^{D}}\in\text{D}$ is a vector of dimension $\text{dim}(\text{D})\ll\text{dim}(\mathcal{H})$, whereas $\ket{\Psi_{\lambda}}\in\mathcal{H}$ is a vector in the entire Hilbert space $\mathcal{H}$.

The problem is now to compute $\Omega^{(n)}$. One may show that Lindgren (1974)

[TABLE]

where $S$ is defined by

[TABLE]

We now spell out the procedure for computing eigenstate matrix elements of local operators to first and second order in perturbation theory for mid-spectrum states deep in the MBL phase.

The required steps are the following:

Select a random “parent state” $\ket{n}$, which is an eigenstate of the unperturbed Hamiltonian $H_{0}$. Its energy $E_{0}^{n}$ will typically lie in the middle of the spectrum of $H_{0}$. 2. 2.

Generate the “family” of states $\ket{m}$ connected to $\ket{n}$ by the perturbation $V$ (i.e. by neighboring pairwise spin flips). 3. 3.

Compare the energies of all states $\ket{m}$ with the parent state $\ket{n}$, include $\ket{n}$ and all states $\ket{m}$ which are quasidegenerate (energy difference $<\alpha/W$) with $\ket{n}$ in the model space D. 4. 4.

For each state $\ket{m}$ which is added to D, its family has to be created and energy differences have to be checked again, possibly including more states in the model space D. States which are well separated from the model space states are included in the complement $\overline{\text{D}}$. 5. 5.

Once this iterative process stops (for the spin-flip operator expectation values further constraints can be used, which simplify this), the effective Hamiltonian $H_{\text{eff}}\in\mathbb{C}^{\text{dim}(\text{D})\times\text{dim}(\text{D})}$ is calculated (see below). 6. 6.

Next the eigenstates of $H_{\text{eff}}$ are computed and we pick one eigenstate $\ket{v_{\lambda}^{\text{D}}}$ at random. 7. 7.

In the next step, the eigenstate of $H_{\text{eff}}$ is promoted to the full Hilbert space, using the perturbation expansion of the wave operator $\Omega$ up to the required order. This step can be skipped using the expressions derived below to directly obtain the eigenstate expectation values of the operators we are after.

We note that since we are dealing with central eigenstates of a quantum spin chain with a dense spectrum, there is no separation of clusters of eigenvalues of $H_{0}$ from other parts of the spectrum. Therefore in principle the procedure above is not guaranteed to yield a model space which is of smaller dimension than the full Hilbert space. It turns out, however, that for the spin-flip correlator, additional selection rules (i.e. the removal of terms yielding zero contributions) guarantee a small model space. For the $S_{i}^{z}S_{i+r}^{z}$ correlators this is not true, and we typically find much larger model spaces, which are however still significantly smaller than the full Hilbert space.

Our procedure is designed to yield the minimal set of states necessary to get the non-vanishing contributions to matrix elements to a given order in perturbation theory.

A.1 Zeroth Order Perturbation Theory

The zeroth order effective Hamiltonian is nothing but the projection $H_{\text{eff}}^{(0)}=\mathbb{P}H\mathbb{P}$ into the model space. Its matrix elements are therefore $\ket{m},\ket{m^{\prime}}\in\text{D}$:

[TABLE]

An eigenstate $\ket{v_{\lambda}^{D}}$ of $H_{\text{eff}}^{(0)}$ with eigenvalue $\lambda$ then yields the corresponding eigenstate in the full Hilbert space

[TABLE]

From this, we may now calculate the zeroth order contribution to the expectation value of an operator $O$:

[TABLE]

A.2 First Order Perturbation Theory

The first order effective Hamiltonian is constructed as follows

[TABLE]

where ${\ket{k}}$ lives in the complement space $\overline{\mathrm{D}}$ and $\ket{m}$, $\ket{m^{\prime}}$ live in the model space $\mathrm{D}$.

To first order, the total effective Hamiltonian is then given by

[TABLE]

The eigenstates $\ket{v_{\lambda}}=\sum_{\mathrm{D}}v^{\lambda}_{m}\ket{m}$ of $H_{\text{eff}}$ with eigenvalue $\lambda$ can now be lifted into the full space by acting with $\Omega$ to first order: $(\Omega^{(0)}+\Omega^{(1)})\ket{v_{\lambda}}$. Thus:

[TABLE]

Then the diagonal matrix elements of operator ${O}$ are

[TABLE]

and we have omitted the one term that appears to order $1/W^{2}$.

How do we decide which states should go into the model space? In order to keep the perturbation consistently of order $1/W$ we should, in principle, keep all states with nonzero amplitude onto $V\ket{n}$ in the model space where $\ket{n}$ is the state chosen initially. However, the expression for the matrix element simplifies matters. We consider operators ${O}$ that are either diagonal in the configuration basis or which flip a pair of spins amounting to a single term in $V$. Then we are justified in restricting the model space to $\ket{n}$ and ${O}\ket{n}$ if the latter lies within an energy window of $1/W$ from $\ket{n}$.

A.3 Second Order Perturbation Theory

We now discuss the perturbative corrections to eigenstates to second order in $V$. The first step in constructing the eigenstates is to include terms to second order in the effective Hamiltonian.

[TABLE]

As above, the total effective Hamiltonian up to second order is given by the sum of all lower order contributions as

[TABLE]

Once we have computed the eigenstates of the effective Hamiltonian to this order, we again lift them into the full Hilbert space by operating on them with $\Omega^{(0)}+\Omega^{(1)}+\Omega^{(2)}$. The first two terms are reported above, the second order contribution is:

[TABLE]

Diagonal matrix elements may now be computed. In full,

[TABLE]

Appendix B Energy dependence of local operators

In the ergodic regime, the validity of the ETH implies that energy eigenstate expectation values of local operators are equivalent to the thermal average at the temperature corresponding to the energy eigenvalue:

[TABLE]

where $E_{n}$ is the energy density of the state $\ket{n}$ and $\beta(E_{n})$ is the temperature chosen such that

[TABLE]

In Fig. 12 the expectation value of the spin-flip correlator and the connected $z$ correlation at distance $r=1,2,3,4$ are plotted as function of energy. The corresponding thermal average matches the mean value of $\bra{n}{O}\ket{n}$ .

We note that the thermal expectation values correspond to the mean of the distributions we consider in the main text. Since the correlators for $r=1$ are terms of the Hamiltonian, they exhibit a significant slope in the middle of the spectrum, leading to the observed sensitivity to the energy target. In the $M=0$ sector $\mathrm{Tr}(H_{M=0})/\text{dim}(\mathcal{H}_{M=0})=-L/(4L-4)$ which sets the infinite temperature limit.

Appendix C Heavy tails in the MBL phase

Some common features of the matrix elements distributions of correlation functions in the MBL phase are the sharp peaks at zero and the presence of tails. In Fig. 13 and 14 such tails are highlighted for different distances using a doubly logarithmic scale. They both show heavy tails, which seem consistent with a power law behavior. Since the connected $S_{i}^{z}S_{i+r}^{z}$ correlator is asymmetric, we show the positive part of the distribution separately from the negative part in Fig. 14 (the sign of the negative part is flipped to show them on the same plot).

Remarkably, the right and left hand tails of the connected $z$ correlation seem to follow the same power law dependence despite their asymmetry, which is clearly visible in this representation.