# Statistics of systemwide correlations in the random-field XXZ chain: Importance of rare events in the many-body localized phase

Jeanne Colbois, Fabien Alet, Nicolas Laflorencie

► To cite this version:

Jeanne Colbois, Fabien Alet, Nicolas Laflorencie. Statistics of systemwide correlations in the random-field XXZ chain: Importance of rare events in the many-body localized phase. Physical Review B, 2024, 110 (21), pp.214210. 10.1103/PhysRevB.110.214210. hal-04757906

HAL Id: hal-04757906

https://hal.science/hal-04757906v1

Submitted on 20 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Statistics of systemwide correlations in the random-field XXZ chain: Importance of rare events in the many-body localized phase

Jeanne Colbois,¹,²,³,* Fabien Alet,¹,† and Nicolas Laflorencie¹,‡

¹Laboratoire de Physique Théorique, Université de Toulouse, CNRS, UPS, France

²MajuLab, CNRS-UCA-SU-NUS-NTU International Joint Research Unit, Singapore

³Department of Materials Science and Engineering and Centre for Quantum Technologies, National University of Singapore, Singapore

(Dated: January 10, 2025)

Motivated by recent debates around the many-body localization (MBL) problem, and in particular its stability against systemwide resonances, we investigate long-distance spin-spin correlations across the phase diagram of the random-field XXZ model, with a particular focus on the strong disorder regime. Building on state-of-the-art shift-invert diagonalization techniques, we study the high-energy behavior of transverse and longitudinal correlation functions, computed at the largest possible distance, for a broad range of disorder and interaction strengths. Our results show that while transverse correlations display a fairly stable exponential decay over the entire XXZ phase diagram, longitudinal correlations exhibit markedly different behavior, revealing distinct physical regimes. More precisely, we identify an intermediate disorder region where standard observables show well-converged MBL behavior [J. Colbois et al., Phys. Rev. Lett. 133, 116502 (2024)] while the distributions of longitudinal correlations reveal unexpected fat-tails towards large values. These rare events strongly influence the average decay of longitudinal correlations, which we find to be algebraic in a broad region inside the supposed MBL phase, whereas the typical decay remains mostly exponential. At stronger disorder and weaker interactions, this intermediate regime is replaced by a more conventional exponential decay with short correlation lengths for both typical and average correlators, as expected for standard localization. Our findings shed light on the systemwide instabilities and raise important questions about the impact of such rare but large long-range correlations on the stability of the MBL phase. Finally, we discuss the possible fate of the intermediate region in the context of recent perspectives in the field.

## CONTENTS

<table><tr><td>I. Introduction</td><td>2</td></tr><tr><td>A. Instabilities</td><td>2</td></tr><tr><td>B. Spin correlations</td><td>2</td></tr><tr><td>C. Main results</td><td>3</td></tr><tr><td>D. Paper organization</td><td>3</td></tr></table>

<table><tr><td>II. Two limiting cases: ETH and Anderson</td><td>4</td></tr><tr><td>A. The many-body Anderson insulator</td><td>4</td></tr><tr><td>1. Exponential decay of the correlation functions</td><td>4</td></tr><tr><td>2. Histograms and correlation map between C<sup>xx</sup><sub>L/2</sub> and C<sup>zz</sup><sub>L/2</sub></td><td>5</td></tr><tr><td>B. Deep in the ergodic phase</td><td>6</td></tr><tr><td>1. ETH and random states</td><td>6</td></tr><tr><td>2. ED results deep in the ergodic regime</td><td>6</td></tr><tr><td>C. Summary of the two limits</td><td>8</td></tr></table>

<table><tr><td>III. Correlations across the phase diagram</td><td>8</td></tr><tr><td>A. Overview</td><td>8</td></tr><tr><td>1. Ergodic regime (h = 3.5)</td><td>8</td></tr><tr><td>2. Strong disorder (h = 6.5, 8.5)</td><td>9</td></tr><tr><td>B. Average vs. typical correlations</td><td>10</td></tr><tr><td>1. Midchain correlation functions</td><td>10</td></tr><tr><td>2. Correlation lengths</td><td>10</td></tr></table>

<table><tr><td>IV. Distributions of correlations, extreme statistics, rare events</td><td>11</td></tr><tr><td>A. Focus on large correlations: general overview</td><td>11</td></tr><tr><td>B. Quantitative analysis of the fat-tail regime</td><td>12</td></tr><tr><td>1. Δ = 1: horizontal scan</td><td>13</td></tr><tr><td>2. h = 6: vertical scan</td><td>14</td></tr><tr><td>C. Discussion</td><td>15</td></tr></table>

<table><tr><td>V. Summary, discussions and outlooks</td><td>15</td></tr><tr><td>Acknowledgments</td><td>17</td></tr><tr><td>A. Many-body Anderson insulator</td><td>17</td></tr><tr><td>1. Transformation to free fermions</td><td>17</td></tr><tr><td>2. Spin-spin correlations</td><td>17</td></tr><tr><td>3. Link between longitudinal and transverse correlations</td><td>18</td></tr></table>

<table><tr><td>B. Random state</td><td>18</td></tr><tr><td>C. Variance</td><td>19</td></tr><tr><td>D. Constant interactions</td><td>19</td></tr><tr><td>E. Constant fields</td><td>21</td></tr><tr><td>1. Fat-tail regime</td><td>21</td></tr><tr><td>2. Weak interaction limit</td><td>22</td></tr></table>

## References

<footnote>* physics.jeannecolbois@gmail.com</footnote>

<footnote>† fabien.aret@cnrs.fr</footnote>

<footnote>‡ nicolas.laflorencie@cnrs.fr</footnote>
# I. INTRODUCTION

Over the past decades, the many-body localization (MBL) phenomenon has attracted a great deal of interest, leading to a significant effort deployed to understand the interacting version of an Anderson insulator [1–10]. Among the crucial points that have been discussed recently, the interpretation of finite-size and/or finite-time numerical simulations plays a central role [11–19]. Indeed, our current understanding of the MBL problem has been strongly shaped by the results of simulations of relatively small interacting quantum disordered systems.

Nevertheless, it is quite widely accepted [20–33] that in one dimension, strong enough disorder favors a stable MBL phase or regime, on the length scales (a few tens of spins) probed by these simulations. However, after almost two decades of very hard work, there is still no clear consensus on some aspects, not least being the very existence of the MBL phase in the thermodynamic limit [34–37], even for the few models that have been intensively studied, such as the random-field Heisenberg chain [19, 21, 24]. One should also mention the few attempts to mathematically prove the existence of the MBL as a new state of matter preventing thermal equilibrium [38–41]. Interestingly, they have all focused on a certain class of models of quantum spin 1/2 chains with no symmetry, in contrast to the by-now standard U(1) symmetric random-field Heisenberg chain. In addition to these theoretical results, experimental evidences of localization have been reported by several groups [42–50], reinforcing the idea that the MBL phase can exist, but also raising the same questions about finite size and finite time limitations.

## A. Instabilities

One important aspect in understanding MBL stability is to identify the relevant mechanisms that could destabilize the interacting insulator. This knowledge would then allow a more complete theoretical modeling of the elusive MBL phenomenon, a problem that has recently benefited from important advances such as the theory of avalanches [51–53] and of many-body resonances [54–59]. Attempts at observing or modeling thermalization through avalanches, as triggered by an ergodic inclusion in an otherwise strongly localized background, have taken several forms. Idealized setups modeling the ergodic inclusion as a random matrix [60–63] are consistent with theoretical predictions [51, 52], even for small systems. Other setups [19, 64], assuming the avalanche has started and modeling its action through a weakly-coupled bath, find a critical disorder much larger than those previously obtained using standard observables [24, 29–32, 65], while with a similar setup Ref. [66] reports a larger range of stability of MBL. Investigations of more realistic models of local thermal inclusions [67–69] reach different conclusions, some finding no or little evidence of effective thermalization by the inclusion [68, 69], some concluding in favor of thermalization albeit in a different model [67]. Remarkably, these ideas have been tested experimentally [50], but again on small chains, making the interpretation in terms of avalanches quite delicate.

Another type of MBL instabilities actively discussed are many-body resonances which, loosely speaking, correspond to the hybridization by interaction terms of “localized” many-body eigenstates differing by spin flips located in a certain spatial range in real-space [39, 54, 55, 57, 70]. The “localized” many-body eigenstates that participate the resonance can be adiabatically connected to the large disorder limit (where a local integral of motions picture apply [71–73]), and have an energy difference exponentially small with the resonance range [57]. Long-range resonances imply hybridization between states that differ by an extensive amount (a fraction of the system size), and are extremely close in energy. Such many-body resonances are at the heart of the formal proofs of stability of MBL [38, 39, 41] and are thought to destabilize the MBL phase. Several prior works found indirect signs of (long-range) many-body resonances in microscopic MBL models [70, 74–77] while more direct approaches [19, 54–56, 58, 78] have also been designed to identify resonances – most of these coming with a certain level of modeling [57, 59].

We finally note that these two destabilization mechanisms need not be unrelated, and Ref. [58] discusses that avalanches primarily propagate through strong, rare, (near-)resonances. The statistics of long-range resonances, including rare events, certainly play a crucial role when approaching the ergodic transition, as been argued e.g. in Refs. [79, 80] borrowing ideas from glassy physics. Nevertheless, a good understanding of systemwide resonances is not yet fully achieved, especially at the microscopic level. Our work aims to present a contribution to this topic through a comprehensive study of simple objects, namely long-range spin correlation functions in eigenstates, which have the potential to capture long-range resonances and allow us to provide a finer insight into the MBL phase diagram, as we show below.

## B. Spin correlations

In the vast literature on MBL and its lattice simulations in U(1) symmetric models, only a handful of papers focused on pairwise correlation functions in eigenstates [21, 76, 77, 81–84], despite the fact that they can be routinely measured in several experimental platforms [43, 44, 47]. In practice, the main efforts have focused on the two-site quantum mutual information [19, 77], an entanglement estimate that directly depends on single-site and two-body correlations [85].

In the present work, we aim to further explore the phase diagram of a prototypical model hosting MBL with the help of long-distance correlations, extending in several ways our recent work [65]. The random-field spin-1/2 XXZ model is defined on a periodic chain of length L, by the Hamiltonian

$$ \mathcal{H}_{\Delta} = \sum_{i=1}^{L} \left( S_{i}^{x} S_{i+1}^{x} + S_{i}^{y} S_{i+1}^{y} + \Delta S_{i}^{z} S_{i+1}^{z} + h_{i} S_{i}^{z} \right), \quad (1.1) $$

where $\Delta \in [0, 1]$ controls the interaction (the Ising anisotropy), and the random fields are drawn from a uniform distribution, $h_i \in [-h, h]$. Crucially, $\mathcal{H}_\Delta$ preserves the total magnetization $S_{\text{tot}}^z = \sum_{i=1}^L S_i^z$. We highly rely on this U(1) symmetry.
Throughout this paper, we consider eigenstates of Eq. (1.1) in the middle of the spectrum for the largest sector of fixed total magnetization. As we focus on even sizes L, we always consider $S_{tot}^z = 0$. Eigenstates are obtained through shift-invert exact diagonalization (ED) [29] on lattices up to length L = 20. For each value of $\Delta$, $h$, we consider between 1000 and 5000 realizations $\{h_i\}$ of disorder.

Before presenting the main results of the current paper, let us recall previous findings obtained in Ref. [65] on the phase diagram for mid-spectrum eigenstates, see Fig. 1. First, we have compiled evidence showing that the Anderson localized (AL) line is unstable for any finite interaction $\Delta$, leading to an ergodic regime when $h < h^* \sim 2-3$, whereas it connects to the MBL regime for sufficiently large disorder. The ergodic to MBL transition line in the disorder $h$ — interaction $\Delta$ plane, has been obtained through an infinite-size extrapolation ($L \to \infty$, symbols in Fig. 1) of finite-size data (up to $L = 21$) of four standard observables: participation entropy (PE), entanglement entropy (EE), extreme magnetizations (EM) and gap ratio (GR). Second, beyond this line, at stronger disorder, we identified a regime of instabilities where, in particular, the spin-spin correlation functions at long-distance ($L/2$) show a non-trivial response upon increasing the system length L.

This last point is the main motivation of the present work, where we propose a comprehensive study of the spin-spin correlation functions as a probe across the phase diagram of the random-field XXZ chain model Eq. (1.1) at high energy, shown in Fig. 1. We denote as follows the connected spin-spin correlation functions

$$C_{ij}^{\alpha\alpha} = \langle S_i^{\alpha} S_j^{\alpha} \rangle - \langle S_i^{\alpha} \rangle \langle S_j^{\alpha} \rangle, \quad (1.2)$$

with $\alpha = x, y, z$ and where the expectation values are calculated for many-body mid-spectrum eigenstates. The XXZ Hamiltonian Eq. (1.1) commutes with $S_{tot}^z$, which implies that the $\alpha = x$ and $\alpha = y$ components, transverse with respect to the random field, are equivalent by symmetry, and that $\langle S_i^x \rangle = 0$, which is not true for the longitudinal component ($\langle S_i^z \rangle \neq 0$).

Throughout this work, we mostly use the absolute value of the correlations [21, 81], so that we can define the typical correlations at distance $r = |i - j|$

$$C_{r,\text{typ}}^{\alpha\alpha} := \exp(\ln |C_r^{\alpha\alpha}|), \quad (1.3)$$

where the bar denotes averaging over positions, eigenstates as well as disorder realizations. Considering the absolute value also in the *average* correlations

$$C_{r,\text{avg}}^{\alpha\alpha} := |C_r^{\alpha\alpha}| \quad (1.4)$$

facilitates the comparison with the typical behavior, and is helpful for the study of the longitudinal correlations, which are in majority negative in the ergodic regime, and fully negative in the non-interacting limit [76]. Because of the exponential decay of various correlations throughout the phase diagram, we often display directly their logarithm.

As a probe for long-range resonances, we focus particularly on the midchain correlations, between spins at the maximal distance $r = L/2$ in periodic chains. We consider their full

probability distributions $P(C_{L/2}^{\alpha\alpha})$, and, to emphasize the typical behavior, $P(\ln|C_{L/2}^{\alpha\alpha}|)$. The main focus is therefore on the decay of the typical and average midchain correlations for different system sizes L, rather than on their decay in the bulk as a function of r for a given system size.

### C. Main results

Figure 1 shows the phase diagram of the random-field XXZ chain Eq. (1.1), across which we have studied in great detail the different behaviors of the spin-spin correlation functions at long distances $C_{L/2}^{\alpha\alpha}$. We thus provide a deep extension of our initial results [65], which can be summarized as follows.

(i) The transverse correlation functions $C_{L/2}^{xx}$ show no significant change throughout the phase diagram: they always show an exponential decay, both for their average and for their typical behavior. While the origin of such exponential decays is very different for weak and strong disorder, $C_{L/2}^{xx}$ is not able to distinguish between the different regimes.

(ii) This is in sharp contrast to the longitudinal correlators $C_{L/2}^{zz}$, which follow distinctly different patterns across the phase diagram: an exponential decay for both typical and average in the strong disorder limit, identified as MBL, precedes a marked difference between typical and average correlators at intermediate disorder, before finally giving rise to the thermal phase characterized by midchain correlations decaying as a power law (due to magnetization conservation) at lower disorder.

(iii) These different behaviors result from a qualitative change in the probability distribution of the longitudinal correlators. Between the strong and weak disorder limits, we show the existence of a rather broad intermediate regime characterized by a fat tail (presumably power-law) for the large correlations, leading to an average dominated by rare events. In this intermediate region, however, the usual observables (spectral statistics, entanglement, participation, extreme magnetization) show well-converged MBL features [31, 65].

Our results are strengthened by a detailed study of the two limiting cases of the thermal phase (at low disorder, finite interaction) and the many-body Anderson localized phase (at $\Delta = 0$). In the main text, results along several scans are discussed: the AL line at $\Delta = 0$; the random-field Heisenberg model (RFHM) line at $\Delta = 1$; and a vertical cut at $h = 6$. Results along other horizontal ($\Delta = 0.25, 0.5, 0.75$) and vertical ($h = 4, 6$) lines are discussed in Appendix.

### D. Paper organization

The rest of the paper is organized as follows. We first discuss in Section II the two limit cases: on the one hand, the many-body Anderson insulator ($\Delta = 0$), and on the other hand, deep
<figure><img src="image_2.png" /></figure>

in the ergodic regime, where the eigenstate thermalization hypothesis (ETH) [86–88] applies. In particular, we study the joint behavior of the longitudinal and transverse correlation functions in these two regimes, which will serve as a guide for our initial exploration of the phase diagram. Section III provides an overview of the phase diagram through the lens of midchain correlators following two representative cuts: one horitzontal along the highly studied random-field Heisenberg model line at $\Delta = 1$, and another one along a vertical cut in the phase diagram at $h = 6$. Typical and average correlations and the related correlation lengths are discussed in detail. Their striking differences motivate a focus on the large correlations, which we discuss in Section IV. Our results are summarized and put in perspective in Section V in relation with the recent literature on instabilities and resonances of the MBL phase, along with a brief conclusion.

Our paper is supported by several appendices. Appendix A provides further details regarding the free fermions results. Appendices B and C discuss respectively the distributions for a random vector in the $S_{tot}^z = 0$ sector and the variances in the interacting case. Appendix D contains a similar analysis as in Sections III and IV for horizontal cuts at $\Delta = 0.25$ and $\Delta = 0.5$. Finally, Appendix E provides further insight into vertical cuts at $h = 4, 6$ and $10$.

## II. TWO LIMITING CASES: ETH AND ANDERSON

To set the stage for the discussion, we present the main features of the spin-spin correlations and their probability distributions, highlighting the differences between two well-controlled limits. We first consider the “many-body” Anderson-localized phase ($\Delta = 0$, Sec. II A), before going deep into the ergodic phase ($\Delta = 1, h = 1$, Sec. II B), where the ETH holds.

FIG. 1. High energy phase diagram of the random-field XXZ chain Eq. (1.1) as obtained for mid-spectrum eigenstates from shift-invert diagonalization. The extrapolated transition points out of the ergodic phase (blue region), as previously obtained [65] using four standard observables, are shown together with new markers (green dots, not extrapolated), dividing the phase diagram into three qualitatively different regimes with respect to the midchain longitudinal correlations $C_{L/2}^{zz}$. Between the red and blue regions, our work highlights the importance of a new intermediate regime (green region) where anomalously large correlation events appear in fat-tailed distributions and dominate the average decay of $C_{L/2}^{zz}$. The estimates of the positions of these green markers are discussed in Appendices D and E. The top insets show the three qualitatively different behaviors for the distributions $P(-\ln|4C_{L/2}^{zz}|)$ for three representative points, indicated by colored stars: (a) self-averaging in the ergodic regime; (b) fat-tail behavior for large correlation events ($\ell^z \to 0$) in the intermediate regime; (c) log-normal-like for strong disorder, where both the average and the typical function decay equally exponentially.

### A. The many-body Anderson insulator

#### 1. Exponential decay of the correlation functions

In the absence of interaction $\Delta = 0$, the XXZ Hamiltonian becomes equivalent (after Jordan-Wigner fermionization) to free fermions in a random potential (half-filling corresponding to $S_{tot}^z = 0$). Using Wick's theorem, it is straightforward [89, 90] to express the various correlation functions in terms of the Anderson localized fermionic orbitals $\phi_m(i)$ (see Appendix A for free-fermion calculations). For instance one easily sees that the longitudinal spin-spin correlations, equivalent to the density-density correlations for fermions, are always negative (for any $i \neq j$):

$$C_{ij}^{zz} = -\left|\langle c_i^\dagger c_j \rangle\right|^2 = -\left|\sum_m \phi_m^*(i)\phi_m(j)\right|^2, \quad (2.1)$$

where $m$ labels the occupied AL orbitals. In contrast, the transverse correlations can be positive or negative

$$C_{ij}^{xx} = \frac{1}{2} \langle S_i^+ S_j^- \rangle \quad (2.2)$$

$$= \frac{1}{2} \langle c_i^\dagger (-1)^{\hat{\phi}_{ij}} c_j \rangle \quad (\text{with } \hat{\phi}_{ij} = \sum_{\ell=i}^{j-1} c_\ell^\dagger c_\ell). \quad (2.3)$$

In this non-interacting, Anderson-localized case, the correlations are expected to decay exponentially with the separation $r = |i - j|$, as shown by ED results in Fig. 2 (a) for the typical correlations Eq. (1.3). We note that the rounding of $\ln|C_r^{\alpha\alpha}|$ observed at large $r$ and fixed size $L$ is an effect of the periodic boundary conditions. Importantly, in this case the bulk decay of the typical correlations for $1 < r < L/2 - 2$:

$$C_{r, \text{typ}}^{\alpha\alpha} \propto \exp\left(-\frac{r}{\xi_{\text{bulk,typ}}^{\alpha}}\right) \quad (2.4)$$
<figure><img src="image_3.png" /><figcaption>FIG. 2. Typical correlations and associated lengths in the XX chain (Anderson localized insulator) at $\epsilon = 0.5$, averaging over $5 \cdot 10^4$ disorder realizations. Small sizes are taken in order to facilitate the comparison with the interacting model. (a) Exponential decay of the typical correlations at $h = 5$. Errors are smaller than the data points. The bulk correlation lengths $\xi_{\text{bulk,typ}}^z = 0.448 \pm 0.005$ and $\xi_{\text{bulk,typ}}^x = 0.868 \pm 0.006$ agree with the midchain correlation lengths $\xi_{\text{typ}}^z = 0.449 \pm 0.003$ and $\xi_{\text{typ}}^x = 0.894 \pm 0.009$. (b) Evolution of the midchain typical correlation length compared to the many-body localization length $\xi_{\text{loc}}$ obtained from the Lyapunov exponent. Dotted lines indicate the fits, Eqs. (2.6), (2.7). At weak disorder, finite-size effects limit the typical midchain correlation length.</figcaption></figure>

agrees with the decay of the midchain typical correlations taken at $r = L/2$ as a function of system size $L$

$$C_{L/2, \text{typ}}^{\alpha\alpha} \propto \exp\left(-\frac{L}{2\xi_{\text{typ}}^{\alpha}}\right). \qquad (2.5)$$

It is instructive to compare these correlation lengths with a many-body version of the Anderson localization length, which can be computed from the Lyapunov exponent associated with single-particle orbitals [53, 91, 92]. The single-particle localization length depends on both disorder strength and single-particle energy. To obtain an overall localization length $\xi_{\text{loc}}$ characteristic of many-body eigenstates at a given many-body energy target (e.g. $\epsilon = 0.5$ [93]), we have to average over the density of states [53, 92]. Fig. 2 (b) shows the disorder dependence of $\xi_{\text{typ}}^{\alpha}$ and $\xi_{\text{loc}}$. Inspired by the fact that the latter is well described by the form [92, 94]

$$\xi_{\text{loc}}(h) = \frac{1}{\ln[1 + (h/h_0)^2]}, \quad h_0 \approx 1.13 \qquad (2.6)$$

we numerically fit the typical correlation lengths and obtain similar behaviors for $h \gtrsim 1.5$

$$\xi_{\text{typ}}^{x}(h) \approx 2\xi_{\text{typ}}^{z}(h) \approx \frac{1}{\ln[1 + (h/h'_{0})^{2}]}, \quad h'_{0} \approx 1.72. \qquad (2.7)$$

The estimated midchain correlation lengths plotted in Fig. 2 (b) deviate slightly from their expected behavior for $h \lesssim 1.5$. For $\xi_{\text{typ}}^z$, this is a very simple finite-size effect due to the small sizes considered here, of the same order as those reachable in the interacting problem. We note that this occurs for $\xi_{\text{typ}}^z \gtrsim 3$. For $\xi_{\text{typ}}^x$, the finite-size effects are more subtle: the correlation length shows a maximum as a function of the disorder strength, drifting to zero disorder as the size increases. For stronger

random fields, the fits are good and the ratio of the transverse to longitudinal correlation lengths is $\approx 2$.

There is a simple argument to explain the factor of 2: for strong enough disorder, the local densities are very close to 1 or 0, and therefore in most cases we expect the Jordan-Wigner phase factors $\hat{\phi}_{ij} = \sum_{\ell=1}^{j-1} c_{\ell}^{\dagger} c_{\ell}$ to contribute almost exclusively to the overall sign in Eq. (2.3), leading to

$$C_{ij}^{xx} \approx \pm \frac{1}{2} \langle c_i^{\dagger} c_j \rangle \qquad (2.8)$$

$$\approx \pm \frac{1}{2} \sqrt{|C_{ij}^{zz}|}, \qquad (2.9)$$

and thus to the observed proportionality between the two components of the correlation lengths

$$\xi_{\text{typ}}^{x} \approx 2\xi_{\text{typ}}^{z}. \qquad (2.10)$$

The average midchain correlations also decay exponentially in the non-interacting case (not shown here, but see Sec. III B), yielding a midchain "average" correlation length:

$$\overline{|C_{L/2}^{\alpha\alpha}|} \propto \exp\left(-\frac{L}{2\xi_{\text{avg}}^{\alpha}}\right) \qquad (2.11)$$

We note two important differences with the typical correlations stemming from the fact that the average correlations are dominated by the large values: the average is more noisy, and the longitudinal and transverse correlation lengths are not related by a simple constant ratio. Instead we observe that $\xi_{\text{avg}}^x < 2\xi_{\text{avg}}^z$, as further discussed in Appendix A.

## 2. Histograms and correlation map between $C_{L/2}^{xx}$ and $C_{L/2}^{zz}$

Keeping in mind our objective of exploring long-range correlations, we focus on the mid-chain correlation functions ($|i - j| = L/2$). In the panels (a) and (b) of Fig. 3, we consider the "correlations of the correlations" at disorder $h = 3.5$, exposing how transverse and longitudinal responses are correlated at the level of individual eigenstates. Most of the weight essentially accumulates along the red dashed line

$$\ln |C_{L/2}^{zz}| = 2 \ln |C_{L/2}^{xx}| + 2 \ln 2, \qquad (2.12)$$

in agreement with Eq. (2.9). This behavior does not change with increasing system sizes, as shown in Fig. 3(b). Although Eq. (2.12) is the typical behavior, which is not rigorously valid for each correlation at the individual eigenstate level, the results in Fig. 3 (a-b) clearly reinforce the above arguments of the irrelevance of the Jordan-Wigner string, leading to Eqs. (2.8) and (2.9), and thus the factor of two in the correlations lengths Eq. (2.10). A simple physical interpretation is that the random magnetic field strongly pins and constrains the spins along the field direction (whether aligned with or against the field), and the dominant quantum fluctuations occur in the transverse channel, perpendicular to the random field direction.

When studied independently, the individual distributions of the xx and zz midchain correlations, plotted in Fig. 3(c,d) for
<figure><img src="image_4.png" /><figcaption>FIG. 3. Distribution of the midchain correlations in the XX chain (Anderson localized) at $\epsilon = 0.5$, for $h = 3.5$, taken over $5 \cdot 10^4$ disorder realizations. The system size color code is the same in all panels. (a) Joint distribution for $L = 18$. The heatmap encodes $P(\ln|C^{xx}|, \ln|C^{zz}|)$, in logarithmic scale. (b) Evolution with increasing system sizes of the joint distribution. For readability we only plot the contour at half the maximum probability (shown in black on panel (a)) (c) and (d) : distribution of the natural logarithm of the longitudinal ($z$) and transverse ($x$) correlations. The bright red dotted lines indicate the maximal value of the correlations, $\ln|C^{\alpha\alpha}| \le -\ln(4)$. The size dependence of the variance of the log distribution is shown in inset of panel (b), with the square root growth marked by a light gray line.</figcaption></figure>

$h = 3.5$, $\Delta = 0$, show very interesting features as system size is increased. Here we focus on small lengths ($L = 8, \dots, 24$) to favor the comparison with the interacting problem later on. Except for the already discussed factor of two, both families of distributions behave similarly:

(i) They are upper bounded by $\ln|C_{L/2}^{\alpha\alpha}| \le -\ln(4)$.

(ii) They drift toward negative values linearly with $L$ as $\overline{\ln|C_{L/2}^{\alpha\alpha}|} \sim -L/2\xi^\alpha$.

(iii) Their standard deviation increases, with a slower power-law compatible with

$$\sigma_{\ln|C_{L/2}^{\alpha\alpha}|} \sim \sqrt{L}, \qquad (2.13)$$

as shown in the inset of Fig. 3(b) (see also Appendix A). We note that the behavior Eq. (2.13) and the presence of logarithmic corrections has been conjectured by Fisher and Young [95] for the related many-body AL regime of the random transverse-field Ising chain. These distributions are not perfectly well captured by (even truncated) simple log-normal distributions [96].

## B. Deep in the ergodic phase

### 1. ETH and random states

In the ergodic phase, a typical “infinite temperature” eigen-state (in the middle of the many-body spectrum, $\epsilon = 0.5$) should be well described by a featureless random state, with no spatial dependence of two-point correlations, namely $C_{ij}^{\alpha\alpha}$ is independent of $|i-j|$. The finite-size effects in expectation values should be entirely ascribed to the finite-dimension $N_H$ of the available Hilbert space, and for some quantities the existence of an underlying sum rule due to a conservation law. In the remainder of this paper, the size of the available Hilbert space $N_H$ will be the sizes of the $S_{tot}^z = 0$ magnetization sector. This leads to markedly different behavior as a function of system size $L$ for longitudinal and transverse correlations, clearly visible in the two leftmost panels (a, b) of Fig. 4, which we now discuss.

a. **Longitudinal correlations**— The conservation of the total magnetization $\sum_{i=1}^L S_i^z = S_{tot}^z$ can be rewritten as follows

$$\sum_{i,j} \langle S_i^z S_j^z \rangle = \frac{L}{4} + \sum_{i \neq j} \langle S_i^z S_j^z \rangle = (S_{tot}^z)^2. \qquad (2.14)$$

For the $S_{tot}^z = 0$ sector (for even $L$), when associated to the featureless property of a random vector, one expects on average

$$\overline{\langle S_i^z S_j^z \rangle} \approx -\frac{1}{4(L-1)} \quad \forall i \neq j. \qquad (2.15)$$

While Eq. (2.15) applies to the average disconnected correlations, it is also valid for the typical, connected correlations, up to corrections decaying exponentially with system size. The absence of dependence on $r$ and the decay of the typical mid-chain correlations with $L$ is readily checked in Fig. 4(a) against exact numerics for fully random states in the $S_{tot}^z = 0$ sector, averaged over hundreds of samples [97].

b. **Transverse correlations**— In contrast to the above slow algebraic decay, the transverse component is suppressed much faster with $L$. Since for a random state the correlations do not depend on the distance, $\langle S_i^\dagger S_j^\dagger \rangle$ is given by a sum of $N_H/4$ variables of order $1/N_H$ with random signs, where $N_H$ is the size of the Hilbert space. We therefore expect

$$\overline{\langle S_i^x S_j^x \rangle} \propto \frac{1}{\sqrt{N_H}} \quad \forall i \neq j, \qquad (2.16)$$

a behavior that is checked for $S_{tot}^z = 0$ random vectors in Fig. 4(b), where we observe

$$C_{L/2, \text{typ}}^{xx} \approx \frac{0.137}{\sqrt{N_H}}, \qquad (2.17)$$

with $N_H = \frac{L!}{(L/2)!^2} \sim 2^L/\sqrt{L}$ for the largest ($S_{tot}^z = 0$) sector.

### 2. ED results deep in the ergodic regime

We now discuss the more realistic situation of the XXZ chain with parameters $h = 1$, $\Delta = 1$, i.e. deep in the ergodic
<figure><img src="image_5.png" /><figcaption>FIG. 4. Typical decay of the spin-spin correlations of a random vector in the $S_{tot}^z = 0$ sector compared to results deep in the ergodic phase at $\Delta = 1, h = 1$. (a) Longitudinal correlations for the random vector. The gray line corresponds exactly to Eq. (2.15). (b) Transverse correlations in the random vector. The gray line is given by Eq. (2.17) (c) Longitudinal correlations in the ergodic phase. The gray line is given by Eq. (2.18). (d) Transverse correlations in the ergodic phase. The gray line is a fit to a square-root decay with the $S_{tot}^z = 0$ sector's size. (e) Joint distribution of the midchain correlations for $L = 18$, with the white contour indicating the height at half-maximum and the plus sign indicating the locations of the maxima. (f) Evolution of the contour at half-maximum of the joint midchain correlations distribution with increasing system sizes. The colored plus signs indicate the locations of the maxima of the distribution. Typical and average values of $|C_{L/2}^{\alpha\alpha}|$ are shown by triangles and crosses respectively. In panels (e) and (f), the AL line is given by Eq. (2.12) while the ETH line is given by Eqs. (2.17) and (2.15).</figcaption></figure>

regime, where we expect the above ETH behavior to apply.

<figure><img src="image_6.png" /><figcaption>FIG. 5. Distributions of the midchain spin-spin correlations deep in the ergodic phase, at $\Delta = 1, h = 1$. (a)-(b) Longitudinal correlations. The majority of the correlations are negative. The variance of the distribution of the logarithm clearly decreases with increasing system sizes. (c)-(d) Transverse correlations. The distributions of $|C_{L/2}^{xx}|$ are gaussian and result in a traveling-wave behavior for the distribution of $\ln|C_{L/2}^{xx}|$. In the bottom panels, the bright red dotted lines indicate the maximal value of the correlations, $\ln|C^{\alpha\alpha}| \le -\ln(4)$. The arrows indicate the typical value of $|C_{L/2}^{\alpha\alpha}|$, while the crosses indicate the average.</figcaption></figure>

Figs. 4 (c, d) show the relatively good agreement with the analytical forms Eq. (2.15) and Eq. (2.16), despite some small discrepancies at short distance.

In order to minimize this non-universal residual effect, we focus on midchain correlations, with a leading term that indeed scales as a power-law for $zz$, shown by a gray line in Fig. 4(c):

$$C_{L/2, \text{typ}}^{zz} \approx \frac{0.27}{L + \ell_0}, \quad (2.18)$$

where the shift $\ell_0 \approx 6.2$ is due to the small size dependence at short distance. The transverse (xx) typical correlator is better described by the random vector results, the gray line corresponding to $C_{L/2, \text{typ}}^{xx} = 0.17/\sqrt{N_H}$ in Fig. 4(d).

In order to emphasize the contrasting behavior of the correlators in the ergodic and Anderson-localized phase, we show in a heatmap the joint probability distribution $P(\ln|C^{xx}|, \ln|C^{zz}|)$ for $L = 18$ in Fig. 4 (e), and the size dependence of its half-maximum in Fig. 4 (f), which can be compared with Fig. 3 (a, b) for the AL case. Clearly, the distribution of correlations between the two components $C_{L/2}^{xx}$ and $C_{L/2}^{zz}$ take a different (triangular-like) shape, with a much smaller variance rapidly decaying with system size. The typical, average and maximal value rapidly converge towards each other, and follow roughly the ETH behavior (from random states, displayed in blue, long-dashed lines in the Figure), up to the correction due to $\ell_0$. This behavior contrasts with the convergence towards the AL line ($\ln|C_{L/2}^{zz}| = 2 \ln|C_{L/2}^{xx}| + 2 \ln 2$), see Fig. 3 (a,b), and the red dashed line in Fig. 4 (f).

In addition, we also show in Fig. 5 the individual distributions $P(C_{L/2}^{\alpha\alpha})$ (panels a,c) and $P(\ln|C_{L/2}^{\alpha\alpha}|)$ (panels b, d). They can be compared with the AL case, displayed in Fig. 3 (c-d), and with random states (see Fig. S2 in Appendix B). This set of plots showcases not only the strong differences between AL and ergodic behaviors, but also between the two different components (transverse and longitudinal) inside the ergodic phase.
In the top panels of Fig. 5, one immediately remarks the strong sharpening with system size of the peaks of $P(C_{L/2}^{\alpha\alpha})$ around their mean values. This is typical of the self-averaging nature of the ergodic phase. More precisely we find (see Appendix C) that the standard deviation of $C_{L/2}^{\alpha\alpha}$ decreases rapidly with the Hilbert space size as

$$ \sigma_{C_{L/2}^{\alpha\alpha}} \sim \frac{1}{\sqrt{N_H}}. \qquad (2.19) $$

Furthermore, in agreement with the random state predictions of Eq. (2.15), we find that the average longitudinal correlation $\overline{C_{L/2}^{zz}}$ takes mostly negative values, in contrast with the transverse correlation $C_{L/2}^{xx}$ which has a symmetric distribution around zero.

The lower panels of Fig. 5 show the probability distributions of the logarithm of the correlators $P(\ln|C_{L/2}^{\alpha\alpha}|)$ where the difference between the longitudinal and transverse components is much more visible. In particular, the standard deviation for the $zz$ component decreases very fast

$$ \sigma_{\ln|C_{L/2}^{zz}|} \sim \frac{L}{\sqrt{N_H}}, \qquad (2.20) $$

as a simple consequence of ETH Eq. (2.19) and the conservation of the total magnetization yielding Eq. (2.15). In contrast, $P(\ln|C_{L/2}^{xx}|)$ shows a clear traveling wave behavior with a constant standard deviation

$$ \sigma_{\ln|C_{L/2}^{xx}|} \sim O(1), \qquad (2.21) $$

which is also a consequence of ETH, since both $\sigma_{C_{L/2}^{xx}}$ and $C_{L/2}^{xx}$ decay $\sim 1/\sqrt{N_H}$.

### C. Summary of the two limits

To summarize the main results in the two limiting cases (many-body AL and deep ergodic case, as captured by random states), we report in Table I the different scaling behaviors as a function of chain length $L$ or Hilbert space size $N_H$ for midchain correlators, related correlation lengths, and variance of probability distribution. We find that the midchain correlators (at $r = L/2$) are particularly useful for several reasons: (i) in the ETH phase, they naturally capture the strongly different decays (exponential versus algebraic) of the $C^{xx}$ and $C^{zz}$ correlators; (ii) in this same phase, they are less affected by finite-size effects that influence short-distance correlators (at small $r$); (iii) in the AL phase they perfectly capture the correlation lengths $\xi^x$ and $\xi^z$, equally well if not better than bulk correlators.

In the next sections, we explore in detail how the correlations and associated correlation lengths evolve as a function of interaction and disorder strengths. In Ref. [65], we described as a smoking gun the fact that the AL hierarchy $\xi_{typ}^x \approx 2\xi_{typ}^z$ evolves to eventually switch, such that $\xi_{typ}^x \approx 1/\ln 2 \ll \xi_{typ}^{z,eff} \to \infty$ in the ergodic phase. Here, our discussion goes much deeper into the distributions, the statistics of rare events, and provides further insight by exploring the subtle differences between typical and average correlations.

<table><caption>TABLE I. Summary of the different behaviors for the transverse (xx) and longitudinal (zz) components of the midchain correlation functions $C_{L/2}^{\alpha\alpha}$ for the ergodic (ETH) and Anderson localized (AL) limits. $N_H$ is the dimension of the $S_{tot}^z = 0$ magnetization sector.</caption><thead><tr><th></th><th>ETH</th><th>Many-body AL</th></tr></thead><tbody><tr><td>Mid-Chain Correlations Functions</td><td>$\overline{\ln|C_{L/2}^{xx}|} \sim -\ln\sqrt{N_H}$<br>$\overline{C_{L/2}^{zz}} \sim -\frac{1}{4(L-1)}$</td><td>$\overline{\ln|C_{L/2}^{xx}|} \sim -\frac{L}{2\xi^x}$<br>$\overline{\ln|C_{L/2}^{zz}|} \sim -\frac{L}{2\xi^z}$</td></tr><tr><td>Correlation lengths</td><td>$\xi_{typ}^x \sim O(1)$, $\xi_{typ}^{z,eff} \to \infty$</td><td>$\xi_{typ}^x = 2\xi_{typ}^z \sim O(1)$</td></tr><tr><td>Standard Deviations</td><td>$\sigma_{\ln|C_{L/2}^{xx}|} \sim O(1)$<br>$\sigma_{\ln|C_{L/2}^{zz}|} \sim \frac{L}{\sqrt{N_H}}$</td><td>$\sigma_{\ln|C_{L/2}^{xx}|} \sim \sqrt{L}$<br>$\sigma_{\ln|C_{L/2}^{zz}|} \sim \sqrt{L}$</td></tr></tbody></table>

## III. CORRELATIONS ACROSS THE PHASE DIAGRAM

### A. Overview

We now consider the fate of both transverse and longitudinal correlation functions over the entire XXZ phase diagram shown in Fig. 1. This is first illustrated in Fig. 6 with a horizontal cut performed at $\Delta = 1$ (the RFHM) for three representative disorder strengths: one point in the ergodic regime ($h = 3.5$) and two points of stronger disorder, presumably belonging to the MBL phase ($h = 6.5$ and $h = 8.5$), at least according to the extrapolation of standard estimates [31, 65].

#### 1. Ergodic regime ($h = 3.5$)

We begin with a discussion of Fig. 6 (a), which shows spin-spin correlation data for $h = 3.5$. Although less deep in the ergodic regime than the previous case ($h = 1$ in Fig. 4 and Fig. 5), these data are consistent with ergodic behavior, and nicely show a crossover to the ETH limit. This is best seen in the lower panel, where the contour at the half-maximum of the joint distribution of the xx and zz midchain correlators $P(\ln|C_{L/2}^{xx}|, \ln|C_{L/2}^{zz}|)$ clearly evolves toward ETH behavior with increasing system size. This trend can also be observed in the typical correlation functions (top panels), where the clear exponential decay of the transverse component $C^{xx}(r)$ differs from the transient behavior of $C^{zz}(r)$, which signals a slow change towards ETH.

The histograms of the midchain correlations $P(\ln|C_{L/2}^{\alpha\alpha}|)$, shown in the middle panels, reveal additional features compared to $h = 1$ in Fig. 5 (c-d). Indeed, the contrasted behaviors of the two components are different. Both distributions broaden with increasing system sizes, while at the same time an increasingly pronounced peak slowly develops for the $z$
<figure><img src="image_7.png" /><figcaption>FIG. 6. Overview of the spin-spin correlations along the RFHM line ($\Delta = 1$) for three representative disorder strengths (a) $h = 3.5$, (b) $h = 6.5$ and (c) $h = 8.5$, for even system sizes from $L = 10$ to $L = 20$. First row: Decay of the typical longitudinal (left) and transverse (right) correlations with the distance $r = |i - j|$. Second row: Probability distribution $P(\ln|C_{L/2}^{\alpha\alpha}|)$. The short vertical lines denote horizontally the typical (triangles) and average (crosses) values of $|C_{L/2}^{\alpha\alpha}|$. The bright red dotted lines indicate the maximal value of the correlations, $\ln|C_{L/2}^{\alpha\alpha}| \le -\ln(4)$. Third row: Evolution of the contour at half-maximum of the joint distribution $P(\ln|C_{L/2}^{xx}|, \ln|C_{L/2}^{xx}|)$ with increasing system sizes. The colored plus signs indicate the locations of the maxima of the distribution. Typical and average values of $|C_{L/2}^{\alpha\alpha}|$ are shown by triangles and crosses respectively. For all disorder strengths, the typical and average values are bounded by the two asymptotic AL and ETH behaviors. The former is represented by a dashed red line corresponding to Eq. (2.12), and the latter by a dashed blue line given by Eqs. (2.17) and (2.15).</figcaption></figure>

component near the average value, signaling the onset of the ETH behavior.

Furthermore, typical $\overline{\ln|C_{L/2}^{\alpha\alpha}|}$ and average $\overline{\ln|C_{L/2}^{\alpha\alpha}|}$ correlations are indicated by down arrows and crosses on the histogram (middle panels), as well as in the $x-z$ correlation map (lower panel). For $h = 3.5$, they show similar trends towards ETH as $L$ increases. The distinction between typical and average behaviors will be discussed below in Sec. III B.

## 2. Strong disorder ($h = 6.5, 8.5$)

For stronger disorder, the situation differs both qualitatively and quantitatively from the ergodic regime, as shown in Fig. 6 (b, c). In the upper panels, typical midchain correlations decay exponentially for both components (discussed further in Sec. III B), although $C_r^{zz}$ shows more rounding in the bulk. Furthermore, the histograms of midchain correlations show a very different behavior compared to both ergodic and AL regimes. While the transverse correlations change only slightly

with $h$, the longitudinal component $C_{L/2}^{zz}$ exhibits several remarkable properties.

(i) Together with a shift towards large negative values, a clear broadening is observed with increasing $L$ (middle panels) [98]. However, unlike the AL results in Fig. 3 (c), a strong asymmetry develops with an accumulation of the weight in the right tail (large correlations). This effect is most pronounced at $h = 6.5$.

(ii) The joint correlations distributions (bottom panels) are very different from the AL case (Fig. 3 (b)). The contours at half maximum encompass a much broader area instead of accumulating on the line given by the typical decay.

(iii) The bottom panels show that average and typical behaviors are different at both quantitative and qualitative levels. This is particularly noticeable at $h = 6.5$, where there is a curvature for the average correlations, but not for the typical.

We explore this last and very significant effect in more detail in the next section, before further discussing important features of the distributions in Sec. IV A.
## B. Average vs. typical correlations

We now address the distinction between average and typical correlations along two representative scans across the phase diagram, namely varying disorder $h$ with fixed interaction $\Delta = 1$ (RFHM), and then varying interaction $\Delta$ for a fixed disorder strength $h = 6$. Other cuts are shown in Appendix D ($\Delta = 0.25, 0.5, 0.75$) and E ($h = 4$, and more details on $h = 6$).

### 1. Midchain correlation functions

a. $\Delta = 1$: horizontal scan— Fig. 7 displays along the RFHM line the evolution with increasing disorder strength ($h = 4.5, \dots, 8.5$) of the average $\ln|C_{L/2}^{\alpha\alpha}|$ and the typical $\overline{\ln|C_{L/2}^{\alpha\alpha}|}$ midchain correlations for both longitudinal ($\alpha = z$) and transverse ($\alpha = x$) components. As already noticed, the transverse correlations remain very well described by a purely exponential decay, with the additional observation that average correlations decay with a larger correlation length: $C_{L/2}^{xx}$ is not sensitive to the transition, and behaves essentially the same in all regimes.

In sharp contrast, the longitudinal correlations $C_{L/2}^{zz}$ show very different properties as $h$ varies. Starting from strong disorder, we first find that average and typical correlations appear to decrease exponentially over the whole range of $L$ for $h \ge 8$. For $8 \ge h \ge 6$, some deviations start to become very visible for the average which no longer follows a pure exponential decay over the full range. Indeed, exponential fits (dashed lines) do not display the same slope depending on the fitting range in the left panels of Fig. 7. The discrepancy becomes increasingly important as $h$ decreases, while at the same time the typical correlations instead seem to remain mostly better described by an exponential decay, at least for $h \ge 6$, where the various fits agree reasonably well. A more systematic analysis in terms of the effective longitudinal correlation length (dependent on the fitting window) is given below in Sec. III B 2.

b. $h = 6$: vertical scan— The results along a vertical cut also display interesting features, as shown in Fig. 8 for increasing interaction ($\Delta = 0, \dots, 0.75$) at $h = 6$. As already observed in Ref. [65], for such a disorder the transverse correlations are practically insensitive to the interaction strength, while the longitudinal component is remarkably enhanced by any $\Delta \ne 0$: this is easily seen by comparing the left and right panels of Fig. 8. Moreover, the distinction between average and typical values is particularly relevant for the $zz$ component. Indeed, as can be seen in the left panels, there is a striking contrast between $\Delta = 0$ (top panel) and an interaction as weak as $\Delta = 0.1$, where the average does not decay as a pure exponential over the whole range, while the typical value does. If we increase $\Delta$ further, the situation remains qualitatively the same, as we now analyze quantitatively by focusing on the correlation lengths.

<figure><img src="image_8.png" /><figcaption>FIG. 7. Evolution of longitudinal (zz, left) and transverse (xx, right) midchain correlations. The typical $\overline{\ln|C_{L/2}^{\alpha\alpha}|}$ (triangles) and average $\overline{\ln|C_{L/2}^{\alpha\alpha}|}$ (crosses) are plotted against $L$ for $\Delta = 1$ and various disorder strengths (shown on plots). A linear scale is used to highlight possible exponential decays. Dashed lines show exponential fits. For $zz$ (left), two fits are shown: light green for the four smallest sizes, and dark blue for the four largest, qualitatively revealing fitting instabilities. For $xx$ (right), we only show one fit.</figcaption></figure>

### 2. Correlation lengths

We first observe that for midchain transverse correlation functions $C_{L/2}^{xx}$, the associated transverse correlation length $\xi^x$ is practically featureless, as it remains $O(1)$ over the entire parameter space. The lengths $\xi^x$ are very stable with system size, and barely depend on the averaging process. This reinforces the result that the midchain transverse correlations always decay exponentially with $L$ over the entire phase diagram. However, it is important to realize that the origin of these exponential decays is very different in the ergodic regime,
<figure><img src="image_9.png" /><figcaption>FIG. 8. Same as for Fig. 7, but for a fixed value of disorder $h = 6$ and various values of interaction (from top to bottom). The first panel corresponds to the non-interacting many-body AL insulator.</figcaption></figure>

where it comes from the simple ETH argument, see Eq. (2.16), while at strong disorder it is a signature of localization.

The longitudinal length turns out to be much more sensitive to the parameters $\Delta$ and $h$, the system size $L$, and the averaging procedure. This can be clearly observed for example at $\Delta = 1$ (Fig. 9 left), where we see that as $h$ decreases, $\xi_{typ}^z$ starts to increase with $L$ when $h \le 5-6$, while $\xi_{avg}^z$ also shows instability with increasing $L$, but already at $h \le 7-8$. A similar effect occurs at $h = 6$ (Fig. 9 right) as a function of the interaction $\Delta$, where one sees that the enhancement with $L$ is much more pronounced for the average length than for the typical. The two lower panels of Fig. 9 provide another clear view of the strong difference between the average and typical correlation lengths, by showing how their ratio evolves. Despite the rather large error bars (due to the fact that we are taking a ratio of two estimates already averaged over the disorder), we can easily see that the ratio increases with $L$

<figure><img src="image_10.png" /><figcaption>FIG. 9. Correlation lengths associated to the typical (top) and average (middle) decays, and their ratio (bottom) in the longitudinal (zz, green) and transverse (xx, orange) directions. Left: $h$ dependence at fixed interaction $\Delta = 1$. Right: $\Delta$ dependence at $h = 6$. Labels denote the various system sizes $L$ used for exponential fits to extract correlation lengths. The top panels were already presented in Ref. [65].</figcaption></figure>

below $h \approx 7$ for $\Delta = 1$ and above $\Delta \approx 0.1$ for $h = 6$, which signals a growing discrepancy between average and typical behavior.

All these observations reinforce the idea that there might be a regime where the decrease of the average correlations is not well captured by a stable exponential decay. This calls for a closer inspection of the distributions of the longitudinal correlators $C_{L/2}^{zz}$, focusing in particular on certain rare events with large values.

## IV. DISTRIBUTIONS OF CORRELATIONS, EXTREME STATISTICS, RARE EVENTS

### A. Focus on large correlations: general overview

To emphasize the large values of correlators, we consider for both components ($\alpha = x, z$) the strictly positive quantities

$$ \ell^{\alpha} = -\ln|4C_{L/2}^{\alpha\alpha}|, \qquad (4.1) $$

and examine their distributions. Compared to the histograms $P(\ln|C_{L/2}^{\alpha\alpha}|)$ plotted in Fig. 6, there is a clear advantage in showing $P(\ell^{\alpha})$ instead as we directly focus on the large correlation events $|C_{L/2}^{\alpha\alpha}| \approx 1/4$ corresponding to the $\ell^{\alpha} \to 0$
<figure><img src="image_11.png" /><figcaption>FIG. 10. Evolution of the distributions $P(\ell^z)$ along the RFHM line ($\Delta = 1$) for a wide range of disorder strengths $h \in [2, 20]$. The histograms, collected across all eigenstates (from full diagonalization) and several thousands of disordered samples, are shown for various chain lengths $L = 8, 10, 12, 14, 16$. As the disorder strength increases, two clear tendencies are observed: there is a first broadening, together with the development of a left fat tail visible for $h$ approximately in the range [4, 10], followed by a gradual vanishing of the fat tail, which then becomes more and more Gaussian-like as $L$ increases. For the largest values of $h$, a peak can be seen at the right-hand limit of the support ($\ell^z \approx 35$), corresponding to vanishingly small correlations that have reached machine precision.</figcaption></figure>

regime.

**Longitudinal correlators** — Fig. 10 gives a broad overview of the evolution of $P(\ell^z)$ for different sizes $L = 8 \cdots 16$ along the RFHM line ($\Delta = 1$) in the phase diagram with $h = 2, 4, \cdots, 20$, thus highlighting three different regimes. Note that to improve the statistics, the histograms are constructed over all the eigenstates, rather than only for the mid-spectrum states. For weak disorder $h = 2$ [panel (a)], the distributions narrow with increasing $L$, in agreement with previous observations [21]. Then, as the disorder strength increases, a new pattern emerges in a regime where otherwise standard observables indicate a rather well-converged MBL-type behavior for these system sizes [31, 65]: the distributions tend to broaden slightly, but more strikingly, a fat tail (seemingly power-law) appears for small $\ell^z$ and then persists for disorder strengths roughly in the range $h \sim 4-10$, see panels (b-e) (this behavior is best seen for the largest $L$). Finally, as the disorder is further increased, this peculiar regime is followed by a more conventional behavior where the distributions become closer to a normal distribution (in the $\ell^z$ variable, with an exponential-like tail not emphasized here) as $L$ increases, see panels (f-h). Note that the observation of this "fat tail" behavior for $P(\ell^z)$ was not reported in Ref. [21], which instead mentioned log-normal distributions once $h > 4$ in the putative MBL phase.

**Comparison with transverse correlator** — Fig. 11 compares $P(\ell^x)$ and $P(\ell^z)$ for three representative disorder strengths ($h = 5, 6.5, 8.5$) in the regime where the fat tail emerges for $\ell^z$ (here, only eigenstates close to the middle of the spectrum are kept). $P(\ell^x)$ is clearly suppressed much faster with system size

and does not exhibit a tail for low values of $\ell^x$. In this regime, $P(\ell^x)$ is qualitatively similar to the Gaussian-like behavior observed for $P(\ell^z)$ at stronger disorders (compare to Fig. 10 (f-h)). The distribution of $\ell^x$ translates into a behavior much more similar to the many-body Anderson insulator (Sec. II A), that is, a standard localized regime in which the statistical analysis is not dominated by rare events.

The different shapes of $P(\ell^z)$, and their evolution with increasing system size, provide a qualitative explanation for the differences between average and typical observed in Sec. III B, a point which we discuss more quantitatively below.

## B. Quantitative analysis of the fat-tail regime

We quantify the importance of the tail in the distribution of $|C_{L/2}^{zz}|$ (gray shaded area in Fig. 11) by the weight

$$W_* = \int_0^{\ell_*^z} P(\ell^z) d\ell^z, \qquad (4.2)$$

where $\ell_*^z$ is a size-independent $O(1)$ cutoff. For a log-normal distribution, assuming that the standard deviation (in $\ell^z$) grows slowly enough [99], $W_*$ would vanish exponentially with $L$, such that average and typical correlations would be similar. The situation is quite different in the presence of a fat tail, which leads to a slower decay of both $W_*(L)$ and the average correlation function, as we discuss below for two representative cuts.
<figure><img src="image_12.png" /><figcaption>FIG. 11. Distributions $P(\ell^z)$ (left) and $P(\ell^x)$ (right) for three disorder strengths $h = 5, 6.5, 8.5$ along the RFHM line ($\Delta = 1$). The histograms, collected from mid-spectrum eigenstates are shown for $L = 10, 12, 14, 16, 18, 20$. The fat-tail regime for the longitudinal component $\ell^z$ contrasts with the quickly vanishing tail for the transverse one. The region indicated by the shaded gray area corresponds to $\ell_*^z < 3.22$ ($|C_{L/2}^{zz}| > 10^{-2}$), and is used to compute the weight $W_*$ in Eq. (4.2), shown in Fig. 12.</figcaption></figure>

1. $\Delta = 1$: horizontal scan

Focusing on the regime $4.5 \le h \le 8.5$, we plot in Fig. 12 both the weight $W_*$ Eq. (4.2), computed for $|C_{L/2}^{zz}| > 0.01$ (corresponding to $\ell^z < 3.22$ [100], shaded area in Fig. 11), and the average correlation, formally defined as

$$ \overline{|C_{L/2}^{zz}|} = \frac{1}{4} \int_{0}^{\infty} P(\ell^z) \exp(-\ell^z) d\ell^z. \quad (4.3) $$

These two quantities show qualitatively similar trends, clearly demonstrating that the average is controlled by the left tail of $P(\ell^z)$, i.e. the weight of large correlations. For most values of $h$ in this regime, we observe a sub-exponential decay with $L$, presumably better described by a power law, as shown in the log-log plots of Fig. 12 (right panels). On the other hand, the data for $h = 8.5$ seem to more reasonably obey a pure exponential decay with $L$, as shown by the semi-log plots in Fig. 12 (left panels).

To make this analysis more systematic across the entire

<figure><img src="image_13.png" /><figcaption>FIG. 12. Average midchain longitudinal correlations $|C_{L/2}^{zz}|$ (top) and large correlation weight $W_*$, Eq. (4.2) (bottom), both plotted against $L$ for $\Delta = 1$ and various disorder strengths $h$, as indicated on the plot. The semi-log scale (left) indicates that only the data at $h = 8.5$ can be reasonably described by an exponential decay, while the log-log scale (right) shows a plausible power-law decay for most of the other data, see Fig. 13 for a quantitative analysis.</figcaption></figure>

horizontal cut ($\Delta = 1, h \le 10$), we fit the decay of both $W_*(L)$ and $|C_{L/2}^{zz}|$ by exponential forms

$$ W_*(L) \sim \exp\left(-\frac{L}{2\lambda_w}\right) \quad \text{and} \quad \overline{|C_{L/2}^{zz}|} \sim \exp\left(-\frac{L}{2\xi_{\text{avg}}^z}\right), \quad (4.4) $$

as well as using algebraic laws

$$ W_*(L) \sim L^{-\eta_w} \quad \text{and} \quad \overline{|C_{L/2}^{zz}|} \sim L^{-\eta_{\text{avg}}^z}. \quad (4.5) $$

As previously done in Sec. IIIB, the length scales $\lambda_w$, $\xi_{\text{avg}}^z$, and the power-law exponents $\eta_{w, \text{avg}}$ are estimated from fits performed over 3 different sliding windows ($L \in [10, 16]; [12, 18]; [14, 20]$) so that one can estimate the (in)stability of the presumed forms Eqs. (4.4, 4.5).

The results are shown in Fig. 13 where we observe a strong similarity between $W_*(L)$ and $|C_{L/2}^{zz}|$, such that $\lambda_w \approx \xi_{\text{avg}}^z$ and $\eta_w \approx \eta_{\text{avg}}^z$. In addition, the study of fit stability with the sliding windows allows one to distinguish three qualitatively different regimes: (i) For strong disorder (above $h \sim 8$) the data are best described by an exponential decay, stable with $L$ and with rather short characteristic length scales; (ii) in the other limit for $h \lesssim 5$, the exponential fits are clearly not stable with $L$, while the power-law fits show a slow $L$-dependent crossover for the exponent $\eta$ possibly at weak disorder towards the value $\eta = 1$ expected for ETH (see Sec. IIB and Fig. 4, this is more clearly seen for $\eta_{\text{avg}}^z$); (iii) finally, the intermediate region $5 \le h \le 8$ is not correctly described by an exponential
<figure><img src="image_14.png" /><figcaption>FIG. 13. $\Delta = 1$. Characteristic lengths $\lambda_w$ and $\xi_{avg}^z$ (bottom) and power-law exponents $\eta_{w,avg}$ (top), obtained by fitting the longitudinal tail weight $W_*$ and correlator $|C_{L/2}^{zz}|$ to respectively exponential and power-law forms Eq. (4.4), and Eq. (4.5)). These estimates, obtained over different sliding windows ($L \in [10, 16]$; $[12, 18]$; $[14, 20]$), are plotted against the disorder strength $h$, highlighting three qualitatively distinct regimes (see text). The insets quantify the qualities of the fitting forms, defined in Eq. (4.6).</figcaption></figure>

decay, and is best characterized by a rather stable algebraic decay, with exponents $2 < \eta < 4$. To confirm this visual estimate of stability of the fits, we compute their coefficient of determination obtained by comparing the residuals of the data $y$ with respect to the fit $f$ from those with respect to the average:

$$R^2 := 1 - \frac{\sum_i (y_i - f_i)^2}{\sum_i (y_i - \bar{y})^2}. \qquad (4.6)$$

$R^2 = 1$ ($R^2 \ll 1$) indicates a perfect (bad) fit. Despite the noise inherent to the original data, the plateau close to 1 for $R^2$ for the power-law fits displayed in the insets in Fig. 13 highlights the existence of an intermediate regime where the power-law fits indeed perform better.

## 2. $h = 6$: vertical scan

An identical analysis is performed for a vertical scan, examining how $P(\ell^z)$ evolves for a fixed disorder $h = 6$ while $\Delta$ is varied. Fig. 14 shows how the weight $W_*(L)$ and the average correlations $|C_{L/2}^{zz}|$ decay with $L$ in the range $\Delta \in [0, 1]$. The results of a quantitative analysis of the decays, fitted with exponential Eq. (4.4) and algebraic Eq. (4.5) forms, are shown in Fig. 15, similar to the analysis in Fig. 13 for the RFHM.

Although the error bars are quite large, especially for the power-law exponents $\eta_{w,avg}$ (top panels), we can still conclude that at sufficiently low interaction, roughly below $\Delta \sim 0.1$, we do not observe signs of fat tails, with both $W^*(L)$ and $|C_{L/2}^{zz}|$

<figure><img src="image_15.png" /><figcaption>FIG. 14. Same as Fig. 12 but for a vertical scan at $h = 6$ plotted for $\Delta \in [0, 1]$. The semi-log scale (left) indicates that the data for $\Delta < 0.1$ can be described by an exponential decay, while with the log-log scale (right) we observe a reasonable power-law decay for $\Delta \gtrsim 0.1$. Fig. 15 shows a quantitative analysis of the fitting parameters.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIG. 15. Same as Fig. 13 but for a vertical scan at $h = 6$ with $\Delta \in [0, 1]$. The power-law regime appears to be stable above $\Delta \sim 0.1$, while below that an exponential decay seems to describe the available data slightly better.</figcaption></figure>

compatible with exponential decay, at least in the range of system sizes accessible to our simulations. This regime would then be adiabatically connected to the many-body Anderson insulator. However, as the interaction strength increases, we observe a qualitative change towards power-law decay at and above $\Delta \sim 0.1$, signaling the onset of the intermediate fat-tail regime. The results in Appendix E support this picture even further.
## C. Discussion

With the above analysis, we have clearly identified a new regime where the midchain longitudinal correlations follow a fat-tailed distribution for $\ell^z \to 0$, presumably power-law. Along the RFHM line ($\Delta = 1$), this roughly occurs for $5 \le h \le 8$, i.e. in a regime where otherwise standard eigen-state observables indicate MBL physics [31, 65]. An important consequence is that in this region, typical and average correlations have qualitatively different decays:

$$|C_{L/2}^{zz}| \propto L^{-\eta_{\text{avg}}} \qquad (4.7)$$

$$\exp\left(\frac{\ln|C_{L/2}^{zz}|}{L}\right) \propto \exp\left(-\frac{L}{2\xi_{z}^{\text{typ}}}\right), \qquad (4.8)$$

suggesting that certain types of rare events with large correlation dominate the mean, leading to a behavior clearly distinct from that of the typical value. The other scan for which we have also presented a detailed analysis at $h = 6$ similarly displays such an intermediate window, readily observed above $\Delta \sim 0.1$ (for lower values of $\Delta < 0.1$, this regime is not observed within our system sizes). Importantly, we did not find trace of a similar regime for the transverse correlator $C_{L/2}^{xx}$, which instead displays more conventional behaviors (Gaussian-like in the ergodic phase, log-normal like at larger disorder).

The fact that typical and average longitudinal correlations show such different trends suggests the absence of self-averaging for this particular observable. Such unconventional behavior reminds us of the case of Infinite Randomness Fixed Point (IRFP) criticality [101–103], for which we briefly recall the peculiar properties of the correlation functions at distance $r$. They are broadly distributed, with a fat tail describing rare but large correlation events associated with long-range singlets (aka the random singlet phase). As a striking consequence, the average decays algebraically $C_{\text{avg}}(r) \sim r^{-\eta}$, since it is dominated by $O(1)$ correlations occurring with a probability $\sim r^{-\eta}$, see Refs. [104, 105] for numerical evidence. At the same time, the typical correlation, controlled by a scale-invariant part of the distribution [95], decays much faster, like a stretched exponential $C_{\text{typ}}(r) \sim \exp(-A\sqrt{L})$. Although IRFP physics is usually associated with a strong disorder type renormalization group flow [106], it is not clear that such a framework provides a good description of MBL physics [107–109]. For example, Refs. [110–112] concluded for $\mathbb{Z}_2$-symmetric models that at high energies, the non-interacting IRFP is destabilized by any weak interactions towards a delocalized ergodic phase, contrasting with the zero-temperature properties where instead IRFP remains stable against finite interactions in the ground state [113]. For the present U(1)-symmetric XXZ model, the non-interacting limit is the many-body Anderson insulator, which is not described by IRFP physics, but rather by a standard finite-disorder fixed point, which has been argued [114] to be analogous to the off-critical (Griffiths) regime of the random transverse-field Ising chain [103, 104]. Therefore, we do not expect finite interaction to emulate any kind of IRFP physics, but our numerical results strongly suggest that rare events may still play a key role.

Remarkably, qualitative differences between average and typical decay have also been reported in apparently quite different situations. In the context of classical spin glasses, Ref. [115] discussed the relevance of large deviations of correlation functions in random magnets using the example of the random-field Ising model on the Bethe lattice, where they found rare events with anomalously strong correlation functions driven by a few pairs of strongly correlated spins. In a series of advances on the Anderson transition problem in random graphs, Refs. [116–118] proposed that rare events dominate the physics in the localized regime. In particular, they numerically found qualitatively different behavior for the decays of single-particle two-point correlation functions over the random graph, a behavior that may be related to our observations, although the correspondence between the many-body problem and random graphs is known to be highly non-trivial [119–126] and deserves more careful investigation.

## V. SUMMARY, DISCUSSIONS AND OUTLOOKS

We summarize here our main results and connect them with the relevant literature. The basis for our comprehensive study of correlators was laid in section II with the limiting cases, where their behavior is well understood. The ergodic regime at low disorder shows different natures for the midchain correlators: algebraic decay with system size $L$ for the longitudinal $C_{L/2}^{zz}$, due to the magnetization preserving nature of $H_\Delta$ in Eq. (1.1), versus the exponentially decaying transverse component $C_{L/2}^{xx}$. In the non-interacting many-body Anderson case, both correlators decay exponentially with associated typical length scales that follow the localization length and respect a hierarchy $\xi_{\text{typ}}^x > \xi_{\text{typ}}^z$.

Section III contrasts the exponential decay of transverse correlators observed throughout the phase diagram with the richer longitudinal behavior. While for sufficiently strong disorder $C_{L/2}^{zz}$ also decays exponentially, upon decreasing disorder we find a quite unconventional intermediate regime with a broad and asymmetric distribution of systemwide correlation $P(\ln|C_{L/2}^{zz}|)$ that appears well before the ergodic phase. In this interesting regime, we report very pronounced differences between average and typical midchain correlators. In particular, the average is no longer well described by a stable exponentially decay, while the typical is (note that the latter will eventually also show signs of instability as the disorder is further reduced, as previously reported in Ref. [65]).

By focusing in Section IV on large correlation events using the logarithmic variable $\ell^z = -\ln(4|C_{L/2}^{zz}|)$, we can explore the large deviation statistics in this intermediate disorder regime, where we reveal the existence of a fat-tailed probability distribution (presumably a power law) for strong correlators $|C_{L/2}^{zz}|$ (small values of $\ell^z$), an effect that is clearly visible in Fig. 10. These rare but strong correlations are at the origin of the non-trivial qualitative difference observed between typical and average decays, as also reported in Fig. 9, where the effective correlation lengths, $\xi_{\text{typ}}^z(L)$ and $\xi_{\text{avg}}^z(L)$, show very contrasting behaviors. This intermediate "fat-tail" regime, in which
systemwide correlations are dominated by rare events, exhibits a rather large extension in the disorder-interaction plane (green region in the phase diagram, Fig. 1), typically nearly as large as the ergodic region. Beyond this regime, for larger disorder, we report more conventional exponentially suppressed long-range correlations, consistent with a more standard real-space localization.

The existence of a region of the phase diagram, or a regime in the parameter range, presumably inside the MBL phase, with unusual (not-strictly localized) properties has been advocated in some previous works, albeit not always at the same location or not delimited by the same markers [127]. An earlier exact diagonalization analysis [81] argued on the existence of a critical regime for $h \in [3, 4]$ based on the power-law decay of the average correlator $C_r^{zz}$, which we instead attribute to the crossover to the ergodic regime in this disorder range, as already visible for the typical correlator represented in Fig. 6 (a). Ref. [128] considered the distribution of sizes of entanglement clusters across the phase diagram and found that the probability of small 2-sites entanglement clusters remains finite for disorder values where we observe the fat-tails, but finds that the two sites are at relatively short distances. Similarly, Ref. [84] obtains relatively local correlation clusters. This contrasts with the long-range nature of the enhanced correlations that we find. Interestingly, Ref. [77] analyzes the probability distribution of (the log of) the quantum mutual information (QMI) between two distant spins (closely related to $C^{zz}$ in case of magnetization conservation) for the RFHM. One of the main results is a stretched exponential behavior (as a function of distance $r$) for the typical QMI of open chains in an intermediate disorder regime, whereas our typical data at maximal distance $r = L/2$ are best fit with an exponential decay (albeit in some cases with a growing correlation length). Ref. [77] also studies events associated to large QMI, but to the best of our understanding does not report power-law tails in the range where we find them, but rather a maximum around $h \approx 3$ (identified as the transition). We note nevertheless the positive skewness in the distribution of the log of the QMI, in the range where we see fat tails. Anomalous dynamical behaviors (in dynamics after a quench, or in transport) have also been reported [12, 14, 18, 70, 129–132] in a similar disorder regime where we find fat-tails for strong correlators, and it would be interesting to understand if these two observations could be related.

Our finding of the enhancement of strong long-distance correlations could be interpreted as a direct witness of long-range many-body resonances, extensively argued [19, 55–58, 70] to be the final precursor to the ergodic phase, and as such complements other (direct or indirect) supports for the identification of long-range many-body resonances in microscopic MBL models [19, 54–56, 58, 59, 70, 74–78, 133].

Considering the logarithmic variable $\ell_z$ effectively zooms on anomalous strong long-distance correlators. With a similar objective but a different approach, Ref. [80] proposes a way to put a biased emphasis on long-distance many-body resonances and also find an intermediate regime (between the ergodic and deeply localized) where such rare long-distance resonances dominate the physical behavior (albeit in a different microscopic model).

In a related spirit, Ref. [19] considers the extreme statistics of the QMI (the maximal QMI taken over all eigenstates at the largest distance for open chains) as a marker of the appearance of at least one systemwide resonance, identified at a disorder strength $h_{sw} \approx 8$ for the RFHM. This value is apparently close to the one where, for $\Delta = 1$, we begin to observe power laws for $P(\ell^z)$, but we note that taking the maximum QMI over the entire many-body spectrum should provide an upper bound on the extent of this regime. Our results suggest, at least for periodic systems free of edge effects, that this upper bound should be larger.

Our study does not shed light on the exact anatomy [120, 123] of the resonant states that give rise to these strong distant correlators, even though we note that that simultaneous strong $C_{L/2}^{zz}$ and vanishing $C_{L/2}^{xx}$ suggests the involvement of (at least two) spin flips in the bulk due to the total magnetization conservation in our samples. Further work is needed to identify the form and range distribution of these strong correlation events, as well as to make connections with phenomenological approaches to many-body resonances [55–57, 59, 134], avalanche models [19, 51–53, 57, 58, 134], and other renormalization group / effective approaches to the transition [52, 94, 117, 135–139].

Although our results are effectively limited to finite systems, one can nevertheless try to elaborate possible scenarios for the intermediate regime in larger chains, and in particular how it relates to its surrounding ergodic and MBL phases.

1) The fat tail of large systemwide correlations may be an early sign of a regime that eventually becomes ergodic at larger scales, a scenario that falls into the phenomenology of the prethermal MBL [59]. However, it remains to be understood why in this regime many-body spectroscopy adheres so strongly to Poisson statistics, as also raised by Biroli et al. [80].

2) As an alternative scenario, one could imagine that the observed power-law decay of the average systemwide correlations is only a transient finite-size effect before an asymptotic exponential decay signaling a stable conventional MBL phase at larger $L$. Although we cannot rule it out, this scenario seems unlikely from our data, which instead show that the boundary of this regime drifts toward larger disorder with increasing $L$.

3) Finally, as a third possibility, one may expect the power-law decay of the average correlations to persist at larger sizes, giving rise to an intermediate stable phase where ergodicity is broken but real-space localization is not total.

This clearly paves the way for interesting future research aimed at exploring these different scenarios. In that respect, tensor network approaches targeting excited states [81, 140–142] may help clarify the fate of the average systemwide correlations for larger systems. Their bias towards less entangled states could, however, be a significant limitation [84, 141].

As a final remark, we believe that further studies of long-distance correlators in microscopic models (with or without magnetization conservation, including in dynamical setups) would provide very useful insights, since they are natural, experimentally accessible [50], and interpretable probes of many-body resonances. This real-space perspective on localization could form a basis for deeper understanding, and reveal the key role of extreme statistics in the stability of MBL.
# ACKNOWLEDGMENTS

We thank G. Biroli, W. Chen, G. Lemarié and M. Tarzia for very fruitful discussions. This work has been partly supported by the EUR grant NanoX No. ANR-17-EURE0009 in the framework of the "Programme des Investissements d'Avenir", is part of HQI initiative (www.hqi.fr) and is supported by France 2030 under the French National Research Agency award number ANR- 22-PNCQ-0002, and also benefited from the support of the Fondation Simone et Cino Del Duca. JC acknowledges support from the Singapore Ministry of Education Academic Research Fund Tier 2 (grant MOE-T2EP50222-0005). We acknowledge the use of HPC resources from CALMIP (grants 2023-P0677) and GENCI (project A0150500225), as well as of the PETSc [143, 144], SLEPc [145, 146], MUMPS [147, 148] and Strumpack [149] sparse linear algebra libraries.

# Appendix

In these appendices, we provide details for several points of the main text. Appendix A deals with the free-fermion calculations. Appendix B gives more details on correlations in random vectors modeling ETH. Appendix C provides a brief discussion of the variance of the correlation distributions in the interacting case. Finally, Appendices D and E show the results of a similar analysis as in the main text, but for other cuts in the phase diagram. Appendix D gives the results for $\Delta = 0.25, 0.5$ and $0.75$ while Appendix E provides further results at constant fields $h = 4$ and $6$. We discuss the evaluation of the boundary of the intermediate regime (green dots in Fig. 1), and of the associated error bars.

## Appendix A: Many-body Anderson insulator

### 1. Transformation to free fermions

The XXZ Hamiltonian Eq. (1.1) can be re-written in terms of spinless fermions in a random potential through the Jordan-Wigner transformation [150]:

$$n_i = S_i^z + 1/2, \quad c_i = \exp\left(i\pi \sum_{j=1}^{i-1} n_j\right) S_i^{-}, \quad (S1)$$

with the term in the exponent corresponding to the Jordan-Wigner string. In one dimension, this string cancels out for nearest-neighbor hopping, yielding

$$H_f = \sum_i \left[ \frac{1}{2} (c_i^\dagger c_{i+1} + c_{i+1}^\dagger c_i + 2\Delta n_i n_{i+1}) + h_i n_i \right] + H_B + C \quad (S2)$$

where the boundary term $H_B$ comes from the Jordan-Wigner string, and $C$ is a constant energy shift. In the Anderson case,

$\Delta = 0$ and the fermionic Hamiltonian is quadratic and can be exactly diagonalized by new fermionic operators of the form

$$b_m = \sum_{i=1}^{L} \phi_m(i) c_i. \quad (S3)$$

Thus we have $H_f = \sum_{m=1}^{L} \epsilon_m b_m^\dagger b_m$ with $\epsilon_m$ the single-particle energy associated with the $m^{th}$ exponentially localized orbital $\phi_m$. The conservation of the magnetization corresponds to a conserved number of spinless fermions, with the $S_{tot}^z = 0$ sector of the spin chain corresponding to $L/2$ occupied orbitals $\{m_\mu\}_{\mu=1}^{L/2}$. To form the many-body wavefunction, we randomly select the occupied orbitals in such a way as to target many-body states at energy density above the ground state $\epsilon = (E - E_{min})/(E - E_{max}) = 0.5$ [151]. For a given disorder realization, the many-body energy is given by $E = \sum_{j=1}^{L/2} \epsilon_{mj}$, with $E_{min}$ ($E_{max}$) corresponding to the ground (top) state.

### 2. Spin-spin correlations

Both the longitudinal and transverse correlations can be expressed in terms of fermions and then evaluated using Wick's theorem, as detailed in particular in Refs. [89, 152]. For $i = j$, we have $\langle S_i^\alpha S_j^\alpha \rangle = \frac{1}{4}$. For $i \neq j$ the longitudinal correlations,

$$\langle S_i^z S_j^z \rangle = \langle (c_i^\dagger c_i - 1/2)(c_j^\dagger c_j - 1/2) \rangle \quad (S4)$$

$$= \langle c_i^\dagger c_i c_j^\dagger c_j \rangle - \frac{1}{2} \langle n_i + n_j \rangle + 1/4 \quad (S5)$$

so that the connected $zz$ correlations read

$$C_{ij}^{zz} = \langle c_i^\dagger c_i c_j^\dagger c_j \rangle - \langle n_i \rangle \langle n_j \rangle = -\langle c_i^\dagger c_j \rangle \langle c_j^\dagger c_i \rangle, \quad (S6)$$

where in the last equality we used Wick's decoupling. This is readily expressed in terms of the occupied orbitals as in Eq. (2.1) and results in purely negative $C^{zz}$ [76].

Ref. [152] shows how to evaluate the transverse correlations. We need to define $A_i = c_i^\dagger + c_i$, $B_i = c_i^\dagger - c_i$. Both these operators intervene when using Wick's theorem, with the only non-zero contraction given by:

$$\langle A_i B_j \rangle = -\langle B_j A_i \rangle = \delta_{ij} - 2 \left( \sum_{\mu=1}^{L/2} \phi_{m_\mu}^*(j) \phi_{m_\mu}(i) \right). \quad (S7)$$

Defining $G_{ij} := \langle B_i A_j \rangle$, we can obtain the transverse correlations from the determinant:

$$\langle S_i^\dagger S_j \rangle = \frac{1}{2} \begin{vmatrix} G_{i,i+1} & G_{i,i+2} & \cdots & G_{i,j} \\ G_{i+1,i+1} & G_{i+1,i+2} & \cdots & G_{i+1,j} \\ \vdots & \vdots & \ddots & \vdots \\ G_{j-1,i+1} & G_{j-1,i+2} & \cdots & G_{j-1,j} \end{vmatrix} \quad (S8)$$

which is easy to evaluate numerically.
### 3. Link between longitudinal and transverse correlations

In Sec. II A, we argued that $\xi_{typ}^x \approx 2\xi_{typ}^z$ but also noticed that $\xi_{avg}^x < 2\xi_{avg}^z$. Here, we use the simplified picture of simple log-normal distributions to explain this observation.

Fig. 3 shows the joint distributions of the logarithm $P(\ln|C_{L/2}^{xx}|, \ln|C_{L/2}^{zz}|)$. The maximal probability lies on the line corresponding to $\ln|C_{L/2}^{zz}| = 2 \ln|C_{L/2}^{xx}| + 2 \ln 2$, implying a factor two between the longitudinal and transverse typical correlation lengths. Yet, the distribution is quite broad, compatible with a standard deviation of the form (see also Appendix C)

$$ \sigma(\ln|C_{L/2}^{\alpha\alpha}|) \sim \beta_{\alpha}\sqrt{(L + o(\ln L))}. \quad (S9) $$

To see how this can impact the average behavior, we first rescale the correlation functions using the Pauli matrices, such that $C^{\alpha\alpha} = 4C^{\alpha\alpha}$. Then, using $\ell^\alpha = -\ln|C_{L/2}^{\alpha\alpha}|$, we can write $C_{typ}^{\alpha\alpha} = \exp(-\ell^\alpha)$ and $C_{avg}^{\alpha\alpha} = \exp(-\ell^\alpha)$ (the subscripts L/2 have been dropped to simplify the notations). For sufficiently strong $h$, we assume that each sample/eigenstate expectation yields

$$ \ell^z = 2\ell^x, \quad (S10) $$

obviously giving the expected result for the typical correlations

$$ C_{typ}^{zz} = \exp(-2\ell^x) = [\exp(-\ell^x)]^2 = (C_{typ}^{xx})^2. \quad (S11) $$

However, for the average we have

$$ C_{avg}^{zz} = \overline{\exp(-2\ell^x)} = \overline{[\exp(-\ell^x)]^2} \quad (S12) $$

$$ = (C_{avg}^{xx})^2 + \sigma_{|C_{avg}^{xx}|}^2 \neq (C_{avg}^{xx})^2 \quad (S13) $$

where we have used the fact that $(C^{xx})^2 = (C_{avg}^{xx})^2 + \sigma_{|C_{avg}^{xx}|}^2$.

If we assume a log-normal $P(C^{\alpha\alpha})$ (implying that $\ell^\alpha$ follows a normal distribution), we obtain

$$ \sigma_{\ell^\alpha}^2 = \ln \left( 1 + \frac{\sigma_{C^{\alpha\alpha}}^2}{(C_{avg}^{\alpha\alpha})^2} \right). \quad (S14) $$

Ignoring the subleading log corrections, we use

$$ \sigma_{\ell^\alpha}^2 = \beta_{\alpha}^2 L, \quad (S15) $$

(in the large-$L$ limit) and invert Eq. (S14), to obtain

$$ C_{avg}^{zz} = (C_{avg}^{xx})^2 \exp(\beta_x^2 L). \quad (S16) $$

Finally, this gives the following relation between the inverse average correlation lengths from midchain correlations

$$ \frac{1}{\xi_{avg}^x} \approx \frac{1}{2\xi_{avg}^z} + \beta_x^2. \quad (S17) $$

Thus the factor 2 present for the typical correlation lengths is replaced by a slowly varying term for the average correlation lengths, as shown in Fig. S1. Recall that Eq. (S17) is valid for log-normal distributions, which do not perfectly describe the transverse and longitudinal correlations in AL. Despite this limitation, the agreement in Fig. S1 is somewhat reasonable.

<figure><img src="image_17.png" /><figcaption>FIG. S1. Relation between the longitudinal and transverse correlation lengths at $\Delta = 0$. While the typical lengths are related by a factor 2 when the disorder is sufficiently strong, there is a further disorder-dependent correction for the average, qualitatively captured by Eq. (S17) despite the approximations.</figcaption></figure>

## Appendix B: Random state

At finite interaction, the limit of weak disorder, i.e. deep in the ergodic regime, was discussed in Sec. II B. In Fig. 5, we presented the distributions of the correlations at $\Delta = 1, h = 1$. For reference, in Fig. S2 we show the distributions of the longitudinal and transverse correlations and of their logarithms for random vectors sampled from the $S_{tot}^z = 0$ sector. Panels (a) and (b) clearly show Gaussian distributions with a reducing variance. Additionally, it is easy to see that longitudinal correlations are shifted to negative values due to the sum rule and that they are almost all negative. Panels (c) and (d) of Fig. S2 illustrate the strong contrast created by this simple sum rule

<figure><img src="image_18.png" /><figcaption>FIG. S2. Distributions of the midchain spin-spin correlations for a random vector in the zero magnetization sector, to compare to Fig. 5.</figcaption></figure>
for the distributions of the logarithms. The transverse correlations have self-similar distributions which only get shifted with increasing system sizes with a constant variance. In contrast, the longitudinal correlations get more and more peaked around the typical value, which decays as a power-law of the system size (see Fig. 4).

Thus, the random vector captures very well the most important effects present at $\Delta = 1, h = 1$ (Fig. 5). The only slight difference is that in Fig. 5, $P(C_{L/2}^{zz})$ has more weight on positive correlations than in the random vector case. This results in the tail at small correlations for $P(\ln |C_{L/2}^{zz}|)$, which however reduces with increasing system sizes.

## Appendix C: Variance

We now mention results for the standard deviations of the logarithm of the longitudinal and transverse correlations. The limiting cases are shown in Fig. S3. In the non-interacting AL insulator, fitting the standard deviation growth with simply a power-law yields $\sigma \sim L^{-0.43}$ at $h = 5$. However, this behavior can be ascribed to additive logarithmic corrections to the variance [95], as stated in Eq. (S9). Fig. S3(a) shows that the numerical results are compatible with such corrections. At $h = 5$, we find:

$$ \sigma(\ln |C_{L/2}^{\alpha\alpha}|) = \beta_{\alpha}\sqrt{(L + 2 \ln L)} \quad \beta_{z} \approx 2\beta_{x} \approx 0.943 \quad (S1) $$

At weak disorder and strong interaction in the ergodic regime, $\sigma(\ln |C_{L/2}^{xx}|)$ is constant as expected from ETH (see the distributions in Fig. 5), whereas $\sigma(\ln |C_{L/2}^{zz}|)$ reaches a maximum and decays with increasing system sizes. More precisely, Fig. S3(b) shows that the decay is roughly related to the square-root of the Hilbert space dimension, see also Section II B).

<figure><img src="image_19.png" /><figcaption>FIG. S3. (a) Standard deviation of the distributions of the log longitudinal and transverse correlations at $\Delta = 0$, and fits to square-root growth with logarithmic corrections. (b) Decay of the variance of the longitudinal correlations, deep in the ETH regime ($h = 0.5, \Delta = 1$).</figcaption></figure>

Fig. S4(a) presents the standard deviation of the logarithm of the longitudinal correlations at $\Delta = 1$. For sufficiently strong disorder, both $\sigma(\ln |C_{L/2}^{zz}|)$ and $\sigma(\ln |C_{L/2}^{xx}|)$ (not shown) grow as a power-law of the system size. The effective exponents in Fig. S4(b) show a very clear square-root behavior for the longitudinal correlations, and a slightly lower power for the transverse correlations, compatible with the logarithmic corrections present also in the non-interacting results.

<figure><img src="image_20.png" /><figcaption>FIG. S4. (a) Standard deviation of the distributions of the log longitudinal correlations at $\Delta = 1$. For disorders $h \ge 3.9$, the standard deviation grows as a power-law with system size $L$. (b) Exponent of the power-law fit for even sizes $L = 8$ to $L = 20$ for both components. At strong disorder for the longitudinal correlations, it is compatible with the expected square-root increase. The transverse standard deviation can be accounted for a logarithmic corrections to a square-root growth.</figcaption></figure>

Taken together with the evolution of the typical correlations, these results are compatible with those of Ref. [21] for the random-field Heisenberg chain, where their $\sigma$ stands for $\sigma(\ln |C^{zz}|)/\sqrt{\ln |C^{zz}|}$. However, our findings in the main text highlight the essential role played by the fat-tail at strong longitudinal correlations, an effect that is not captured by the scaling of $\sigma(\ln |C^{zz}|)$.

## Appendix D: Constant interactions

In the main text, we have focused on the correlation functions in the random-field Heisenberg model. Here, we provide an analysis of other constant interaction lines, at $\Delta = 0.25, 0.5$ and $0.75$. We focus on the typical and average decays as well as on the weight of large correlations (fat-tails).

Fig. S5 provides an overview of the evolution of the typical longitudinal and transverse correlation lengths. It can be compared to Figs. 1 and 3 of Ref. [65] as well as Fig. 9 of the present work. Altogether, these figures support the following important results (see also main text):

(i) The typical longitudinal length grows quickly with interactions, contrasting with the almost constant transverse correlation length.

(ii) Throughout the presumably MBL phase (as given by standard observables), this results in a change of ordering from $\xi_{typ}^z < \xi_{typ}^x$ at strong disorder and weak interactions (similar to the AL limit) to $\xi_{typ}^z > \xi_{typ}^x$ at weaker disorder and stronger interactions (Fig. S5).

(iii) An instability of the longitudinal typical length shows up in roughly the region where $\xi_{typ}^z \sim \xi_{typ}^x$.

We suggested in the main text that this instability could be a weak signature of rare, large correlations, that are much more manifest in the average correlations. Here, we provide additional evidence supporting this scenario throughout the phase diagram. Despite the large errors, Fig. S6 signals a region of power-law decay of the average correlations, supporting further our results in Sec. IV.
<figure><img src="image_21.png" /><figcaption>FIG. S5. Overview of the midchain correlation lengths throughout the phase diagram (compare to Fig. 9). The lengths associated with the midchain average correlations are shown in Fig. S6. The middle panel was already shown in [65].</figcaption></figure>

Using these results to trace out the intermediate regime in Fig. 1 is not easy for the following reason. At very strong disorder ($h \gtrsim 9$), the average correlations fluctuate a lot. Importantly, this purely numerical effect shows up already in the Anderson localized case for strong disorder $h \gtrsim 8$ – in some cases, the decay can seem algebraic even though a pure exponential is expected (not shown here). The limit on small correlations imposed by machine precision together with the low probability of large correlations imply that sampling effects become much more important. Upon introducing interactions, the challenging region starts from stronger disorders (as the correlations grow). Hence, this difficulty does not undermine our results at disorders below $h \sim 8-9$ in the phase diagram.

From Fig. 1, we extract the boundary and corresponding errors as follows:

1. For $\Delta = 0.25$ in Fig. S6(a), the fits to power-law decays of both the average correlations and of the weight

<figure><img src="image_22.png" /><figcaption>FIG. S6. Overview of the midchain correlations throughout the phase diagram (compare to Fig. 13). Left column: decay of the weight of large correlations fitted to a power-law (top panels) and exponential (bottom panels). Right column: power-law (top) and exponential (bottom) fits to the decay of the average $|C_{L/2}^{zz}|$.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>FIG. S7. Evolution of the distributions $P(\ell^z)$ along a constant field line $h = 6$ for interactions from $\Delta = 0.01$ to $\Delta = 1$. As in Fig. 11, and unlike Fig. 10, the histograms are collected across about 60 eigenstates near $\epsilon = 0.5$ (and less than 1% of the spectrum, for small sizes) and thousands of disordered samples. Small values of $\ell^z$ correspond to very large correlations (close to the maximal possible value), while large values of $\ell^z$ correspond to extremely small correlations.</figcaption></figure>

of large correlations are consistently good until $h \le 6$ (presumably even $h \le 6.5$), while $h \ge 7$ is on the edge with an exponential fit seeming slightly better. The error on the right is estimated in light of the fact that the decay at $h = 8$ is clearly exponential, suggesting that the change of behavior occurs before.

2. For $\Delta = 0.5$ (panel (b)), the situation is much more delicate, with a lot of fluctuations in the longitudinal correlations. However, the $R^2$ (Eq. (4.6)) suggests that the correlations decay is slightly better described by a power-law at least until $h \le 6.5$, which provides the lower bound on the error in Fig. 1, while most points beyond this are compatible with an exponential decay. The very large error in the phase diagram corresponds to the unexpected, probably accidental, agreement with a power-law decay at $h = 9, \Delta = 0.5$.

3. For $\Delta = 0.75$ in panel (c), the situation is again much clearer, with a power-law decay until $h \le 7$, an intermediate situation for $h \sim 8$ followed by an exponential decay at $h \sim 9$: we place the upper limit for the end of the power-law regime before this value.

In our view, despite the challenges, these results present a consistent picture for the phase diagram. The next section

supports this even further by providing a similar analysis, this time for fixed disorder strength and variable interaction.

## Appendix E: Constant fields

### 1. Fat-tail regime

In Sec. IV B 2, we have presented a quantitative analysis of the decays of the average correlations and the weights of the fat-tails. Here, we provide some further evidence, qualitative for $h = 6$, and quantitative for $h = 4$.

First, Fig. S7 gives a qualitative overview of the distributions at constant fields. The distributions for weak interactions [panels (a)-(c)] exhibit a strikingly different behavior from strong interactions [panels (f)-(h)]. In the first panels, the large correlations are suppressed very quickly with increasing system sizes, in agreement with the exponential decay of the associated weights Figs. 14 and 15. Instead, at strong interactions, a very clear power-law tail develops at small $\ell^z$, i.e. large correlations, again agreeing with the power-law decay of the weights of large correlations shown in Fig. 14. In the intermediate values of interactions, $0.1 \le \Delta \le 0.25$, there are some signatures of the presence of a fat-tail regime. While it
is difficult to identify the precise limiting value of the inter-
actions from only a qualitative inspection of the distributions,
the integrated weight in Fig. 15 indicates clearly that all this
region also corresponds to a power-law decay. Taken together,
the results at *h* = 6 do provide ample evidence for the exist-
ence of a wide power-law regime for the decay of the average
correlations.

<figure><img src="image_24.png" /></figure>

FIG. S8. Same as Fig. 15 but for *h* = 4. For this field, we note that for the finite sizes presented here, all standard observables show non-ergodic behavior for any *Δ* ≤ 0.5. Note that in Fig. 1, the extrapolated critical disorder out of the ergodic phase for *Δ* = 0.25 is of order *h*<sub>c</sub> = 3.7±0.6. This suggests that above *Δ* ~ 0.25, the power law decay of the average correlations may be due to the presence of the ergodic phase.

To further establish that a truly *exponential* decay of average
longitudinal correlations can be found at weaker (but non-zero)
interactions, we complement our results at *h* = 6 with evidence
at *h* = 4, shown in Figure S8. For *Δ* ≤ 0.05, this plot reveals a
consistent exponential decay of the average correlations, with
the top panels clearly ruling out a power-law decay. In contrast,
for *Δ* ≥ 0.1, the results are consistent with a power-law.

## 2. Weak interaction limit

We finally point out an intriguing effect of the weak in-
teractions on the joint distributions of mid-chain correlation,
$P(\ln|C_{L/2}^{xx}|, \ln|C_{L/2}^{zz}|)$. In Sec. II A and in particular Fig. 3, we
have seen that the majority of the weight accumulates along
the “AL line” given by Eq. (2.12). In Sec. III and Fig. 6,
we have seen very clearly that even at strong interactions and
strong disorder (in the presumably MBL regime), the corre-
lations instead have much broader distributions, noticeable in
the contours at half-maximum. The typical correlations also
follow a very different trend as the one shown by the AL line.
In Fig. S9, we put the emphasis on how this behavior arises at
weak but non-zero interaction.

For very weak interactions (orange, $\Delta = 0.01$), the contour
at half maximum for $L = 12$ accumulates on the AL line.

However, the relation between transverse and longitudinal cor-
relations (orange triangles) exhibits a slightly different slope.
Upon increasing the system size, eventually, the contour also
starts detaching from the AL accumulation line. This shows
the extreme sensitivity to vanishingly small but finite $\Delta$ on the
Anderson insulator.

At $\Delta = 0.25$ (blue), the slope of the typical correlations is
markedly different, and the contours at half-maximum cover
a much broader region. They clearly detach from the AL
accumulation line for increasing system sizes. For $L = 12$, the
most probable point ("+" symbol) is still on the AL line, but
similarly to what happens in Fig. 6 at $\Delta = 1$, $h = 8.5$, it jumps
away for a large enough size.

<figure><img src="image_25.png" /></figure>

FIG. S9. Contours at half-maximum weight of $\mathcal{P}(\ln|C_{L/2}^{xx}|, \ln|C_{L/2}^{zz}|)$ at $h = 6$, for sizes and interaction strengths as indicated in the legend, with the location of maximal probability indicated by a "+" sign. The down triangles illustrate the typical correlations for even sizes $L = 8, \dots, 20$.
[1] L. Fleishman and P. W. Anderson, Interactions and the Anderson transition, Phys. Rev. B 21, 2366 (1980).

[2] B. L. Altshuler, Y. Gefen, A. Kamenev, and L. S. Levitov, Quasi-particle Lifetime in a Finite System: A Nonperturbative Approach, Phys. Rev. Lett. 78, 2803 (1997).

[3] P. Jacquod and D. L. Shepelyansky, Emergence of Quantum Chaos in Finite Interacting Fermi Systems, Phys. Rev. Lett. 79, 1837 (1997).

[4] I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Interacting Electrons in Disordered Wires: Anderson Localization and Low-T Transport, Phys. Rev. Lett. 95, 206603 (2005).

[5] D. M. Basko, I. L. Aleiner, and B. L. Altshuler, Metal-insulator transition in a weakly interacting many-electron system with localized single-particle states, Annals of Physics 321, 1126 (2006).

[6] V. Oganesyan and D. A. Huse, Localization of interacting fermions at high temperature, Phys. Rev. B 75, 155111 (2007).

[7] R. Nandkishore and D. A. Huse, Many-Body Localization and Thermalization in Quantum Statistical Mechanics, Annu. Rev. Condens. Matter Phys. 6, 15 (2015).

[8] F. Alet and N. Laflorencie, Many-body localization: An introduction and selected topics, Comptes Rendus Physique 19, 498 (2018).

[9] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Colloquium: Many-body localization, thermalization, and entanglement, Rev. Mod. Phys. 91, 021001 (2019).

[10] P. Sierant, M. Lewenstein, A. Scardicchio, L. Vidmar, and J. Zakrzewski, Many-body localization in the age of classical computing, arXiv:2403.07111 (2024).

[11] E. V. H. Doggen, F. Schindler, K. S. Tikhonov, A. D. Mirlin, T. Neupert, D. G. Polyakov, and I. V. Gornyi, Many-body localization and delocalization in large quantum chains, Phys. Rev. B 98, 174202 (2018).

[12] F. Weiner, F. Evers, and S. Bera, Slow dynamics and strong finite-size effects in many-body localization with random and quasiperiodic potentials, Phys. Rev. B 100, 104204 (2019).

[13] R. K. Panda, A. Scardicchio, M. Schulz, S. R. Taylor, and M. Žnidarič, Can we study the many-body localisation transition?, EPL 128, 67003 (2020).

[14] T. Chanda, P. Sierant, and J. Zakrzewski, Time dynamics with matrix product states: Many-body localization transition of large systems revisited, Phys. Rev. B 101, 035148 (2020).

[15] P. Sierant, D. Delande, and J. Zakrzewski, Thouless Time Analysis of Anderson and Many-Body Localization Transitions, Phys. Rev. Lett. 124, 186601 (2020).

[16] D. A. Abanin, J. H. Bardarson, G. De Tomasi, S. Gopalakrishnan, V. Khemani, S. A. Parameswaran, F. Pollmann, A. C. Potter, M. Serbyn, and R. Vasseur, Distinguishing localization from chaos: Challenges in finite-size systems, Annals of Physics 427, 168415 (2021).

[17] E. V. Doggen, I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Many-body localization in large systems: Matrix-product-state approach, Annals of Physics 435, 168437 (2021).

[18] P. Sierant and J. Zakrzewski, Challenges to observation of many-body localization, Phys. Rev. B 105, 224203 (2022).

[19] A. Morningstar, L. Colmenarez, V. Khemani, D. J. Luitz, and D. A. Huse, Avalanches and many-body resonances in many-body localized systems, Phys. Rev. B 105, 174205 (2022).

[20] M. Žnidarič, T. Prosen, and P. Prelovsek, Many-body localization in the Heisenberg XXZ magnet in a random field, Phys. Rev. B 77, 064426 (2008).

[21] A. Pal and D. A. Huse, Many-body localization phase transition, Phys. Rev. B 82, 174411 (2010).

[22] J. H. Bardarson, F. Pollmann, and J. E. Moore, Unbounded Growth of Entanglement in Models of Many-Body Localization, Phys. Rev. Lett. 109, 017202 (2012).

[23] A. D. Luca and A. Scardicchio, Ergodicity breaking in a model showing many-body localization, EPL 101, 37003 (2013).

[24] D. J. Luitz, N. Laflorencie, and F. Alet, Many-body localization edge in the random-field Heisenberg chain, Phys. Rev. B 91, 081103(R) (2015).

[25] S. Bera, H. Schomerus, F. Heidrich-Meisner, and J. H. Bardarson, Many-Body Localization Characterized from a One-Particle Perspective, Phys. Rev. Lett. 115, 046603 (2015).

[26] Y. Bar Lev, G. Cohen, and D. R. Reichman, Absence of Diffusion in an Interacting System of Spinless Fermions on a One-Dimensional Disordered Lattice, Phys. Rev. Lett. 114, 100601 (2015).

[27] M. Serbyn, Z. Papić, and D. A. Abanin, Criterion for many-body localization-delocalization phase transition, Phys. Rev. X 5, 041047 (2015).

[28] J. Gray, S. Bose, and A. Bayat, Many-body localization transition: Schmidt gap, entanglement length, and scaling, Phys. Rev. B 97, 201105(R) (2018).

[29] F. Pietracaprina, N. Macé, D. J. Luitz, and F. Alet, Shift-invert diagonalization of large many-body localizing spin chains, Sci-Post Physics 5, 045 (2018).

[30] N. Macé, F. Alet, and N. Laflorencie, Multifractal Scalings Across the Many-Body Localization Transition, Phys. Rev. Lett. 123, 180601 (2019).

[31] P. Sierant, M. Lewenstein, and J. Zakrzewski, Polynomially filtered exact diagonalization approach to many-body localization, Phys. Rev. Lett. 125, 156601 (2020).

[32] N. Laflorencie, G. Lemarié, and N. Macé, Chain breaking and Kosterlitz-Thouless scaling at the many-body localization transition in the random-field Heisenberg spin chain, Phys. Rev. Research 2, 042033(R) (2020).

[33] S. Roy and D. E. Logan, Fock-space anatomy of eigenstates across the many-body localization transition, Phys. Rev. B 104, 174201 (2021).

[34] J. Suntajs, J. Bonca, T. Prosen, and L. Vidmar, Quantum chaos challenges many-body localization, Phys. Rev. E 102, 062144 (2020).

[35] D. Sels and A. Polkovnikov, Dynamical obstruction to localization in a disordered spin chain, Phys. Rev. E 104, 054105 (2021).

[36] D. Sels and A. Polkovnikov, Thermalization of dilute impurities in one-dimensional spin chains, Phys. Rev. X 13, 011041 (2023).

[37] A. Weisse, R. Gerstner, and J. Sirker, Operator growth in disordered spin chains: Indications for the absence of many-body localization, arXiv:2401.08031 (2024).

[38] J. Z. Imbrie, On Many-Body Localization for Quantum Spin Chains, J Stat Phys 163, 998 (2016).

[39] J. Z. Imbrie, Diagonalization and many-body localization for a disordered quantum spin chain, Phys. Rev. Lett. 117, 027201 (2016).

[40] C. Yin, R. Nandkishore, and A. Lucas, Eigenstate localization in a many-body quantum system, Phys. Rev. Lett. 133, 137101 (2024).

[41] W. D. Roeck, L. Giacomin, F. Huveneers, and O. Prosniak, Absence of normal heat conduction in strongly disordered interacting quantum chains, arXiv:2408.04338 (2024).
[42] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch, Observation of many-body localization of interacting fermions in a quasirandom optical lattice, *Science* **349**, 842 (2015).

[43] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, Many-body localization in a quantum simulator with programmable random disorder, *Nature Physics* **12**, 907 (2016).

[44] J.-y. Choi, S. Hild, J. Zeiher, P. Schauß, A. Rubio-Abadal, T. Yef-sah, V. Khemani, D. A. Huse, I. Bloch, and C. Gross, Exploring the many-body localization transition in two dimensions, *Science* **352**, 1547 (2016).

[45] H. P. Lüschen, P. Bordia, S. Scherg, F. Alet, E. Altman, U. Schneider, and I. Bloch, Observation of Slow Dynamics near the Many-Body Localization Transition in One-Dimensional Quasiperiodic Systems, *Phys. Rev. Lett.* **119**, 260401 (2017).

[46] P. Roushan, C. Neill, J. Tangpanitanon, V. M. Bastidas, A. Megrant, R. Barends, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. Fowler, B. Foxen, M. Giustina, E. Jeffrey, J. Kelly, E. Lucero, J. Mutus, M. Neeley, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. White, H. Neven, D. G. Angelakis, and J. Martinis, Spectroscopic signatures of localization with interacting photons in superconducting qubits, *Science* **358**, 1175 (2017).

[47] A. Lukin, M. Rispoli, R. Schittko, M. E. Tai, A. M. Kaufman, S. Choi, V. Khemani, J. Léonard, and M. Greiner, Probing entanglement in a many-body-localized system, *Science* **364**, 256 (2019).

[48] M. Rispoli, A. Lukin, R. Schittko, S. Kim, M. E. Tai, J. Léonard, and M. Greiner, Quantum critical behaviour at the many-body localization transition, *Nature* **573**, 385 (2019).

[49] Q. Guo, C. Cheng, Z.-H. Sun, Z. Song, H. Li, Z. Wang, W. Ren, H. Dong, D. Zheng, Y.-R. Zhang, and others, Observation of energy-resolved many-body localization, *Nature Physics* **17**, 234 (2020).

[50] J. Léonard, S. Kim, M. Rispoli, A. Lukin, R. Schittko, J. Kwan, E. Demler, D. Sels, and M. Greiner, Probing the onset of quantum avalanches in a many-body localized system, *Nature Physics* **19**, 481 (2023).

[51] W. De Roeck and F. Huveneers, Stability and instability towards delocalization in many-body localization systems, *Phys. Rev. B* **95**, 155129 (2017).

[52] T. Thiery, F. Huveneers, M. Müller, and W. De Roeck, Many-Body Delocalization as a Quantum Avalanche, *Phys. Rev. Lett.* **121**, 140601 (2018).

[53] P. J. D. Crowley and A. Chandran, Avalanche induced coexisting localized and thermal regions in disordered chains, *Phys. Rev. Research* **2**, 033262 (2020).

[54] B. Villalonga and B. K. Clark, Eigenstates hybridize on all length scales at the many-body localization transition, *arXiv:2005.13558* (2020).

[55] S. J. Garratt, S. Roy, and J. T. Chalker, Local resonances and parametric level dynamics in the many-body localized phase, *Phys. Rev. B* **104**, 184203 (2021).

[56] S. J. Garratt and S. Roy, Resonant energy scales and local observables in the many-body localized phase, *Phys. Rev. B* **106**, 054309 (2022).

[57] P. J. D. Crowley and A. Chandran, A constructive theory of the numerically accessible many-body localized to thermal crossover, *SciPost Phys.* **12**, 201 (2022).

[58] H. Ha, A. Morningstar, and D. A. Huse, Many-body resonances in the avalanche instability of many-body localization, *Phys. Rev. Lett.* **130**, 250405 (2023).

[59] D. M. Long, P. J. D. Crowley, V. Khemani, and A. Chandran,

Phenomenology of the prethermal many-body localized regime, *Phys. Rev. Lett.* **131**, 106301 (2023).

[60] D. J. Luitz and others, How a Small Quantum Bath Can Thermalize Long Localized Chains, *Phys. Rev. Lett.* **119**, 150602 (2017).

[61] J. Šuntajs and L. Vidmar, Ergodicity breaking transition in zero dimensions, *Phys. Rev. Lett.* **129**, 060602 (2022).

[62] K. Pawlik, P. Sierant, L. Vidmar, and J. Zakrzewski, Many-body mobility edge in quantum sun models, *Phys. Rev. B* **109**, L180201 (2024).

[63] J. Šuntajs, M. Hopjan, W. De Roeck, and L. Vidmar, Similarity between a many-body quantum avalanche model and the ultra-metric random matrix model, *Phys. Rev. Res.* **6**, 023030 (2024).

[64] D. Sels, Bath-induced delocalization in interacting disordered spin chains, *Phys. Rev. B* **106**, L020202 (2022).

[65] J. Colbois, F. Alet, and N. Laflorencie, Interaction-driven instabilities in the random-field xxz chain, *Phys. Rev. Lett.* **133**, 116502 (2024).

[66] A. Scocco, G. Passarelli, M. Collura, P. Lucignano, and A. Russomanno, Thermalization propagation front and robustness against avalanches in localized systems, *arXiv* (2024), 2407.20985.

[67] J. C. Peacock and D. Sels, Many-body delocalization from embedded thermal inclusion, *Phys. Rev. B* **108**, L020201 (2023).

[68] L. Colmenarez, D. J. Luitz, and W. De Roeck, Ergodic inclusions in many-body localized systems, *Phys. Rev. B* **109**, L081117 (2024).

[69] T. Szołdra, P. Sierant, M. Lewenstein, and J. Zakrzewski, Catching thermal avalanches in the disordered xxz model, *Phys. Rev. B* **109**, 134202 (2024).

[70] S. Gopalakrishnan, M. Müller, V. Khemani, M. Knap, E. Demler, and D. A. Huse, Low-frequency conductivity in many-body localized systems, *Phys. Rev. B* **92**, 104202 (2015).

[71] D. A. Huse, R. Nandkishore, and V. Oganesyan, Phenomenology of fully many-body-localized systems, *Phys. Rev. B* **90**, 174202 (2014).

[72] M. Serbyn, Z. Papić, and D. A. Abanin, Local Conservation Laws and the Structure of the Many-Body Localized States, *Phys. Rev. Lett.* **111**, 127201 (2013).

[73] V. Ros, M. Müller, and A. Scardicchio, Integrals of motion in the many-body localized phase, *Nucl. Phys. B* **891**, 420 (2015).

[74] S. D. Geraedts, R. Nandkishore, and N. Regnault, Many-body localization and thermalization: Insights from the entanglement spectrum, *Phys. Rev. B* **93**, 174202 (2016).

[75] V. Khemani, S. P. Lim, D. N. Sheng, and D. A. Huse, Critical Properties of the Many-Body Localization Transition, *Phys. Rev. X* **7**, 021013 (2017).

[76] L. Colmenarez, P. A. McClarty, M. Haque, and D. J. Luitz, Statistics of correlation functions in the random Heisenberg chain, *SciPost Phys.* **7**, 064 (2019).

[77] B. Villalonga and B. K. Clark, Characterizing the many-body localization transition through correlations, *arXiv:2007.06586* (2020).

[78] J. A. Kjäll, Many-body localization and level repulsion, *Phys. Rev. B* **97**, 035163 (2018).

[79] C. Monthus, Many-Body-Localization Transition in the Strong Disorder Limit: Entanglement Entropy from the Statistics of Rare Extensive Resonances, *Entropy* **18**, 122 (2016).

[80] G. Biroli, A. K. Hartmann, and M. Tarzia, Large-deviation analysis of rare resonances for the many-body localization transition, *Phys. Rev. B* **110**, 014205 (2024).

[81] S. P. Lim and D. N. Sheng, Many-body localization and transition by density matrix renormalization group and exact diagonalization studies, *Phys. Rev. B* **94**, 045111 (2016).
[82] A. K. Kulshreshtha, A. Pal, T. B. Wahl, and S. H. Simon, Approximating observables on eigenstates of large many-body localized systems, Phys. Rev. B 99, 104201 (2019).

[83] V. K. Varma, A. Raj, S. Gopalakrishnan, V. Oganesyan, and D. Pekker, Length scales in the many-body localized phase and their spectral signatures, Phys. Rev. B 100, 115136 (2019).

[84] K. Hémery, F. Pollmann, and A. Smith, Identifying correlation clusters in many-body localized systems, Phys. Rev. B 105, 064202 (2022).

[85] G. De Tomasi, S. Bera, J. H. Bardarson, and F. Pollmann, Quantum mutual information as a probe for many-body localization, Phys. Rev. Lett. 118, 016804 (2017).

[86] J. M. Deutsch, Quantum statistical mechanics in a closed system, Phys. Rev. A 43, 2046 (1991).

[87] M. Srednicki, Chaos and quantum thermalization, Phys. Rev. E 50, 888 (1994).

[88] M. Rigol, V. Dunjko, and M. Olshanii, Thermalization and its mechanism for generic isolated quantum systems, Nature 452, 854 (2008).

[89] E. Lieb, T. Schultz, and D. Mattis, Two soluble models of an antiferromagnetic chain, Annals of Physics 16, 407 (1961).

[90] P. Pfeuty, The one-dimensional Ising model with a transverse field, Annals of Physics 57, 79 (1970).

[91] C. A. Müller and D. Delande, Disorder and interference: localization phenomena, arXiv:1005.0915 (2024).

[92] J. Colbois and N. Laflorencie, Breaking the chains: Extreme value statistics and localization in random spin chains, Phys. Rev. B 108, 144206 (2023).

[93] Here $\epsilon = (E - E_{\min}) / (E_{\max} - E_{\min}) \in [0, 1]$ defines the (many-body) energy density above the ground-state energy $E_{\min}$, with $E$ the total energy and $E_{\max}$ the maximal energy.

[94] A. C. Potter, R. Vasseur, and S. A. Parameswaran, Universal Properties of Many-Body Delocalization Transitions, Phys. Rev. X 5, 031033 (2015).

[95] D. S. Fisher and A. P. Young, Distributions of gaps and end-to-end correlations in random transverse-field Ising spin chains, Phys. Rev. B 58, 9131 (1998).

[96] The large correlations are captured correctly but not the weak-correlations tail.

[97] More than 1500 samples for $L = 18$ and below, and more than 600 for $20 \le L \le 24$.

[98] One can notice in Fig. 6 for $h = 8.5$ and $L = 20$ the weakest correlations start to reach the numerical precision when $\ln|C| < -36$.

[99] It is sufficient that $\sigma$ grows slowly enough compared to $\overline{l^z}$. Here we have $\sigma \sim \sqrt{L}$, see Appendix C.

[100] The precise value of $\ell_*^z$ is not crucial, as long as it remains $O(1)$ and not too small so that we can numerically resolve the weight $W_*$.

[101] D. S. Fisher, Random transverse field Ising spin chains, Phys. Rev. Lett. 69, 534 (1992).

[102] D. S. Fisher, Random antiferromagnetic quantum spin chains, Phys. Rev. B 50, 3799 (1994).

[103] D. S. Fisher, Critical behavior of random transverse-field Ising spin chains, Phys. Rev. B 51, 6411 (1995).

[104] A. P. Young and H. Rieger, Numerical study of the random transverse-field Ising spin chain, Phys. Rev. B 53, 8486 (1996).

[105] N. Laflorencie, H. Rieger, A. W. Sandvik, and P. Henelius, Crossover effects in the random-exchange spin-1/2 antiferromagnetic chain, Phys. Rev. B 70, 054430 (2004).

[106] F. Igloi and C. Monthus, Strong disorder RG approach of random systems, Phys. Rep. 412, 277 (2005).

[107] D. Pekker, G. Refael, E. Altman, E. Demler, and V. Oganesyan, Hilbert-Glass Transition: New Universality of Temperature-

Tuned Many-Body Dynamical Quantum Criticality, Phys. Rev. X 4, 011052 (2014).

[108] K. Slagle, Y.-Z. You, and C. Xu, Disordered XYZ spin chain simulations using the spectrum bifurcation renormalization group, Phys. Rev. B 94, 014205 (2016).

[109] Y.-Z. You, X.-L. Qi, and C. Xu, Entanglement holographic mapping of many-body localized system by spectrum bifurcation renormalization group, Phys. Rev. B 93, 104205 (2016).

[110] R. Sahay, F. Machado, B. Ye, C. R. Laumann, and N. Y. Yao, Emergent Ergodicity at the Transition between Many-Body Localized Phases, Phys. Rev. Lett. 126, 100604 (2021).

[111] S. Moudgalya, D. A. Huse, and V. Khemani, Perturbative instability towards delocalization at phase transitions between MBL phases, arXiv:2008.09113 (2020).

[112] N. Laflorencie, G. Lemarié, and N. Macé, Topological order in random interacting Ising-Majorana chains stabilized by many-body localization, Phys. Rev. Research 4, L032016 (2022).

[113] N. Chepiga and N. Laflorencie, Resilient Infinite Randomness Criticality for a Disordered Chain of Interacting Majorana Fermions, Phys. Rev. Lett. 132, 056502 (2024).

[114] H. Javan Mard, J. A. Hoyos, E. Miranda, and V. Dobrosavljević, Strong-disorder renormalization-group study of the one-dimensional tight-binding model, Phys. Rev. B 90, 125141 (2014).

[115] F. Morone, G. Parisi, and F. Ricci-Tersenghi, Large deviations of correlation functions in random magnets, Phys. Rev. B 89, 214202 (2014).

[116] I. García-Mata, O. Giraud, B. Georgeot, J. Martin, R. Dubertrand, and G. Lemarié, Scaling Theory of the Anderson Transition in Random Graphs: Ergodicity and Universality, Phys. Rev. Lett. 118, 166801 (2017).

[117] I. García-Mata, J. Martin, R. Dubertrand, O. Giraud, B. Georgeot, and G. Lemarié, Two critical localization lengths in the Anderson transition on random graphs, Phys. Rev. Res. 2, 012020(R) (2020).

[118] I. García-Mata, J. Martin, O. Giraud, B. Georgeot, R. Dubertrand, and G. Lemarié, Critical properties of the anderson transition on random graphs: Two-parameter scaling theory, kosterlitz-thouless type flow, and many-body localization, Phys. Rev. B 106, 214202 (2022).

[119] M. Tarzia, Many-body localization transition in hilbert space, Phys. Rev. B 102, 014208 (2020).

[120] S. Roy and D. E. Logan, Fock-space correlations and the origins of many-body localization, Phys. Rev. B 101, 134202 (2020).

[121] S. Roy and D. E. Logan, Localization on certain graphs with strongly correlated disorder, Phys. Rev. Lett. 125, 250402 (2020).

[122] F. Pietracaprina and N. Laflorencie, Hilbert-space fragmentation, multifractality, and many-body localization, Annals of Physics 435, 168502 (2021).

[123] S. Roy and D. E. Logan, The fock-space landscape of many-body localisation, arXiv:2408.12685 (2024).

[124] J.-N. Herre, J. F. Karcher, K. S. Tikhonov, and A. D. Mirlin, Ergodicity-to-localization transition on random regular graphs with large connectivity and in many-body quantum dots, Phys. Rev. B 108, 014203 (2023).

[125] K. Tikhonov and A. Mirlin, From anderson localization on random regular graphs to many-body localization, Annals of Physics 435, 168525 (2021), special Issue on Localisation 2020.

[126] T. Scoquart, I. V. Gornyi, and A. D. Mirlin, Role of fock-space correlations in many-body localization, Phys. Rev. B 109, 214203 (2024).

[127] The fat-, power-law, tails that we report for the probability $P(\ell^z)$ should not be confused with the “heavy-tails” of Ref. [76].
which finds power-law behavior for the RFHM for $P(|C_r^{\alpha\alpha}|)$ for intermediate values of $|C|$, at strong disorder (deep in the MBL regime), and for small distances $r \le 4$ and both directions $\alpha = x, z$.

[128] L. Herviou, S. Bera, and J. H. Bardarson, Multiscale entanglement clusters at the many-body localization phase transition, Phys. Rev. B 99, 134205 (2019).

[129] O. S. Barišić, J. Kokalj, I. Balog, and P. Prelovšek, Dynamical conductivity and its fluctuations along the crossover to many-body localization, Phys. Rev. B 94, 045126 (2016).

[130] P. Prelovšek, M. Mierzejewski, O. Barišić, and J. Herbrych, Density correlations and transport in models of many-body localization, Annalen der Physik 529, 1600362 (2017).

[131] S. Nandy, F. Evers, and S. Bera, Dephasing in strongly disordered interacting quantum wires, Phys. Rev. B 103, 085105 (2021).

[132] F. Evers, I. Modak, and S. Bera, Internal clock of many-body delocalization, Phys. Rev. B 108, 134204 (2023).

[133] R. Ghosh and M. Žnidarič, Resonance-induced growth of number entropy in strongly disordered systems, Phys. Rev. B 105, 144203 (2022).

[134] P. J. D. Crowley and A. Chandran, Mean-field theory of failed thermalizing avalanches, Phys. Rev. B 106, 184208 (2022).

[135] T. Thiery, M. Müller, and W. D. Roeck, A microscopically motivated renormalization scheme for the MBL/ETH transition, arXiv:1711.09880 (2017).

[136] P. T. Dumitrescu, A. Goremykina, S. A. Parameswaran, M. Serbyn, and R. Vasseur, Kosterlitz-Thouless scaling at many-body localization phase transitions, Phys. Rev. B 99, 094205 (2019).

[137] A. Morningstar and D. A. Huse, Renormalization-group study of the many-body localization transition in one dimension, Phys. Rev. B 99, 224205 (2019).

[138] A. Morningstar, D. A. Huse, and J. Z. Imbrie, Many-body localization near the critical point, Phys. Rev. B 102, 125134 (2020).

[139] S. Roy, Spectral multifractality and emergent energyscales across the many-body localisation transition, arXiv:2404.07975 (2024).

[140] V. Khemani, F. Pollmann, and S. Sondhi, Obtaining Highly Excited Eigenstates of Many-Body Localized Hamiltonians by the Density Matrix Renormalization Group Approach, Phys. Rev. Lett. 116, 247204 (2016).

[141] X. Yu, D. Pekker, and B. K. Clark, Finding Matrix Product State Representations of Highly Excited Eigenstates of Many-Body Localized Hamiltonians, Phys. Rev. Lett. 118, 017201 (2017).

[142] B. Villalonga, X. Yu, D. J. Luitz, and B. K. Clark, Exploring one-particle orbitals in large many-body localized systems, Phys. Rev. B 97, 104406 (2018).

[143] S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, J. Faibussowitsch, W. D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan, J. Sarich, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, and J. Zhang, PETSc/TAO Users Manual, Tech. Rep. ANL-21/39 - Revision 3.21 (Argonne National Laboratory, 2024).

[144] S. Balay, W. D. Gropp, L. C. McInnes, and B. F. Smith, Efficient management of parallelism in object oriented numerical software libraries, in Modern Software Tools in Scientific Computing, edited by E. Arge, A. M. Bruaset, and H. P. Langtangen (Birkhäuser Press, 1997) pp. 163–202.

[145] V. Hernandez, J. E. Roman, and V. Vidal, SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems, ACM Trans. Math. Software 31, 351 (2005).

[146] J. E. Roman, C. Campos, L. Dalcin, E. Romero, and A. Tomas, SLEPc Users Manual, Tech. Rep. DSIC-II/24/02 - Revision 3.21 (D. Sistemes Informàtics i Computació, Universitat Politècnica de València, 2024).

[147] P. Amestoy, I. S. Duff, J. Koster, and J.-Y. L'Excellent, A fully asynchronous multifrontal solver using distributed dynamic scheduling, SIAM Journal on Matrix Analysis and Applications 23, 15 (2001).

[148] P. Amestoy, A. Buttari, J.-Y. L'Excellent, and T. Mary, Performance and Scalability of the Block Low-Rank Multifrontal Factorization on Multicore Architectures, ACM Transactions on Mathematical Software 45, 2:1 (2019).

[149] P. Ghysels, X. S. Li, F.-H. Rouet, S. Williams, and A. Napov, An efficient multicore implementation of a novel hss-structured multifrontal solver using randomized sampling, SIAM Journal on Scientific Computing 38, S358 (2016).

[150] P. Jordan and E. P. Wigner, Uber das paulische aquivalenzverbot, in The Collected Works of Eugene Paul Wigner, The Collected Works of Eugene Paul Wigner No. A / 1, edited by A. S. Wightman (Springer Berlin Heidelberg, 1993) pp. 109–129.

[151] M. Hopjan, G. Orso, and F. Heidrich-Meisner, Detecting delocalization-localization transitions from full density distributions, Phys. Rev. B 104, 235112 (2021).

[152] P. Henelius and S. M. Girvin, Numerical study of the random dimerized XX spin-1/2 chain, Phys. Rev. B 57, 11457 (1998)