Multifractality and its role in anomalous transport in the disordered XXZ spin-chain

The disordered XXZ model is a prototype model of the many-body localization transition (MBL). Despite numerous studies of this model, the available numerical evidence of multifractality of its eigenstates is not very conclusive due severe finite size effects. Moreover it is not clear if similarly to the case of single-particle physics, multifractal properties of the many-body eigenstates are related to anomalous transport, which is observed in this model. In this work, using a state-of-the-art, massively parallel, numerically exact method, we study systems of up to 24 spins and show that a large fraction of the delocalized phase flows towards ergodicity in the thermodynamic limit, while a region immediately preceding the MBL transition appears to be multifractal in this limit. We discuss the implication of our finding on the mechanism of subdiffusive transport.

connectivity between them is mediated by the Hamiltonian (cf. discussion in [33]). Since the structure of this graph is rather involved it is normally approximated by either the Bethe lattice [34] or random-regular graphs (RRG, see also review by Imbrie et al. [35]). The first proposal of an intermediate nonergodic extended phase sandwiched between the deeply ergodic and insulating (MBL) phases appeared almost 20 years ago [34]. This phase, colloquially dubbed by Altshuler a "bad metal" [36], was defined as a phase where the eigenfunctions are extended over the Hilbert space, but cover only N γ states, where γ < 1 and N is the Hilbert space dimension. Whether such an intermediate phase, with multifractal eigenfunctions, exists for the Anderson localization problem on the Bethe lattice or RRGs, is still an ongoing debate. Large scale studies on random regular graphs (RRGs) suggest that this phase disappears in the thermodynamic limit [37][38][39][40][41], although there is also no consensus here [37][38][39][42][43][44][45][46][47]. In addition, for weak disorder where all researchers agree that eigenfunctions are ergodic on RRGs, subdiffusion has been recently observed [48,49]. Notwithstanding, while Anderson localization on graphs and MBL are related, it is not clear whether results from RRGs apply for MBL.
Multifractal properties of eigenstates of systems which exhibit MBL where examined in a number of studies [50][51][52][53][54]. The outcome is however rather inconclusive, mostly due to presence of severe finite size effects. While Ref. [50] suggests that there is no intermediate multifractal phase, Refs. [51,55] argue in favor of a stable intermediate phase.
In Refs. [52,56] multifractal properties of matrix elements of local operators are studied and found to be multifractal, though Ref. [52] argued that the intermediate phase shrinks to the MBL critical point in the thermodynamic limit, as it occurs in the standard Anderson transition. There is therefore a need for a large-scale, numerical study, which attempts to resolve these discrepancies, and shed light whether multifractality is related to the anomalous dynamical properties of the delocalized phase. Two multifractal moments of the disordered XXZ model were studied in Refs. [50,53], and suggest that the extended phase is ergodic. While we see similar behavior of the relevant moments, in our work we find them insufficient to unveil possible nonergodic behavior, which becomes only apparent at higher moments. Our analysis thus allows us to locate a region in the extended phase which appears to be nonergodic. Our study also provides, for the first time, the presentation of the multifractal spectrum. We are able to identify a large portion of the delocalized phase, where anomalous transport was previously observed, but which is consistent with a transient multifractality. Our results support the existence of multifractality in a region which precedes the MBL transition, although we cannot say whether this region shrinks to the critical point when the system size is increased (cf. [52]).

II. MODEL
In this work we analyze the properties of the disordered XXZ chain, which is given by the Hamiltonian, wereŜ z i , is the z−projection of the spin-1/2 operator,Ŝ ± i are the corresponding lowering and raising operators, J xy and J z are inter-spin couplings and h i are random magnetic fields taken to be uniformly distributed in the interval h i ∈ [−W, W ] . This model conserves the z−projection of the total spin, and serves as the prototypical model of the MBL transition, which for infinite temperature occurs for W ∼ 3.7 [50,57,58]. For W 3.7 the system is in a MBL phase, with a completely suppressed transport of all globally conserved quantities [4], while for W 3.7 it exhibits an anomalous transport with a dynamical exponents which depends on the disorder strength [11][12][13][14][15]. Since one of the major purposes of this work is to study the connection between anomalous transport and multifractality, we limit this study to disorder strengths sufficiently far from the MBL transition, W ≤ 3.

III. RESULTS
Multifractal analysis requires the calculation of the eigenstates of (1) in a certain energy density window and for a large number of disorder realizations. Since full diagonalization becomes overwhelmingly expensive for system sizes L 18, and access to large system sizes is essential, we utilize the shift-invert technique [59], which transforms the spectrum of the Hamiltonian such that the states of interest are moved to the lowest (highest) energies in the transformed spectrum and become tractable by Krylov space methods. The most commonly used spectral transformation for this purpose is (H − σI) −1 , where the explicit inversion of the shifted Hamiltonian can be avoided and replaced by the solution of a set of linear equations using the Gauss algorithm. We use the massively parallel strumpack library [60,61] to extract about 50 eigenstates in the middle of the many-body spectrum, where the density of states is at its maximum. The largest system size we consider is L = 24, which corresponds to a Hilbert space dimension of N = 2 704 156. We repeat this procedure for 100 − 15 000 realizations of the disordered magnetic field h i in (1). Overall, the total number of calculated eigenstate coefficients in the computational basis for each disorder strength and system size is 10 8 − 10 10 , which in most cases allows us to reach statistical errors smaller than the symbol size.

A. Distributions of eigenstate coefficients
The high energy states of ergodic systems are well approximated by eigenstates of random-matrices drawn from a Wigner-Dyson ensemble of matrices [62] which shares the same temporal symmetry as the Hamiltonian. Specifically, eigenstates of real ergodic Hamiltonians, which are time-inversion invariant, are well described by eigenstates of matrices drawn from the Gaussian Orthogonal Ensenble (GOE), suggesting that the elements of eigenstates, |β , written in a certain basis |n , are almost independent random variables, normally distributed according to, where we defined the random variable, x ≡ | n|β | and N is the Hilbert space dimension. This assertion, known as Berry's conjecture [63], was verified numerically in several single-and many-body ergodic systems (see Ref. [64] for a review). Multifractal eigenstates, on the contrary, do not satisfy Berry's conjecture, but are distributed according to, where f (α) is a function called the spectrum of fractal dimensions [3], depending on the only variable α taken to be α ≡ − ln x 2 / ln N (see Section III C). The distribution of eigenstate coefficients for our model (1) has been studied by two of us in Ref. [22], and was found to exhibit significant deviations from Berry's conjecture for 0.4 ≤ W ≤ 1.8, hinting that the underlying eigenstates are multifractal. In this work we study these distributions in detail, focusing on their flow towards the thermodynamic limit. In both above cases we focus on eigenstate coefficients in the computational basis, where the basis states |n are labeled by the eigenvalues of the localŜ z i operators. To give equal weight to small and large values of the eigenstate coefficients, the bins of the histogram are equally spaced on a logarithmic scale, which we achieve by calculating the histogram of α ≡ − ln | n|β | 2 / ln N using bins of equal size. The histogram of the wavefunction coefficients | n|β | and the corresponding probability density can then be straightforwardly inferred. In Fig. 1 we show the result of this calculation for a number of disorder strengths in the extended phase, and a range of system sizes. We compare these distributions with the normal distribution of GOE, (2), and see a visible departure for all disorder strengths, similarly to Fig. 2 in Ref. [22]. The departure is especially apparent in the head of the distribution, indicating an excess in small values of the eigenstate elements compared to GOE, and the tails of the distribution, indicating an excess in large values of the eigenstate elements. The departure becomes more prominent with the strength of the disorder.
While at first glance the rescaled distributions look collapsed, a more detailed examination by zooming into various parts of the distribution, shows a noticeable, yet slow, flow towards the (Gaussian) GOE distribution in most parts of the distribution. In what follows we examine this flow in detail, by considering the moments of the distribution and its multifractal spectrum (3).  In the multifractal analysis one defines the standard inverse participation ratio (IPR) I β 2 and its generalizations,

B. Moments of the distributions: inverse participation ratios
which measure how many "sites" in the Hilbert space (a site here is a certain basis state |n ) the wavefunction occupies [3], the generalized IPR is directly related to the corresponding q Rényi entropies S q = ln I β q /(1 − q) [65]. For eigenstates extended over the entire basis, such as for eigenstates extracted from GOE, x ≡ | β|n | ∼ N −1/2 giving, I β q = n N −q = N −(q−1) and τ q = q − 1 for q > −0.5 (the average of I β q diverges otherwise). Eigenstates which occupy a finite number of configurations |n which doesn't scale with the Hilbert space dimension N will have τ q = 0 for q > 0 (and −∞ otherwise). The parameter q is used to tune the weight in the average from large to small values. Under the assumption that the eigenstate coefficients are statistically independent, the IPRs are related to the moments of P (x), since one can write, Using the definition of τ q in (4), the normalization of the wavefunction, which gives, I β q=1 = 1, and the fact that L ∈ [14,24] L ∈ [18,24]  n 1 = N , which gives I β q=0 = N one can show that in the limit of N → ∞, τ q is a monotonically increasing and concave function of q, namely τ q > 0 and τ q < 0 [3].
To evaluate the τ q we calculate the IPRs for each eigenfunction and a range 0 ≤ q ≤ 4. We then average I β q over the nearby in energy eigenstates, as also different disorder realizations, and obtain I N q . The finite-size average τ avg q is the given by, Since the relation (4) is only expected to hold asymptotically, in Fig. 2 we plot τ avg q (N ) vs 1/ ln N and extrapolate to N → ∞ using a linear function in 1/ ln N . [66] The extrapolated values are then plotted as a function of q and compared to the prediction of GOE, τ q = q − 1 (dashed black line). While the scaling of τ avg q (N ) with respect to 1/ ln N is mostly linear indicating a high quality of the extrapolation, for larger values of q a departure from the linear dependence is apparent, suggesting that the data is still far from being asymptotic. To quantify the curving of the data, we extrapolate to N → ∞ using a sliding fit window of system sizes, which are shown in the legend. The error bars in the extrapolated data are estimated using a bootstrap fitting procedure, quantifying the statistical errors in τ avg q (N ) (which are in all cases smaller than the symbol size). From the extrapolated data we see that the average, τ avg q ∼ q − 1 for all values of q up to some q * (W, L), which depends on the strength of the disorder. While for weak disorder q * spans the entire range of q considered here, for W → W c we see that q * → 1 . On the other hand we also see that, q * (W, L) increases when higher weight in the extrapolation is given to the larger system sizes, suggesting that the departure from the q − 1 could be a finite size effect, though we cannot rule out a saturation of the form lim L→∞ q * (W, L) = q * (W ), which will indicate residual multifractality at large moments q > q * (W ).
We also study the typical τ typ q (N ), which is defined as The advantage of this measure is that it suppresses the weight of the outliers. The results of the same analysis for τ typ q as described above for τ avg q are presented in Fig. 3, and are qualitatively identical to the analysis of τ avg q . Quantitatively q * is pushed to much larger values [67], almost entirely eliminating the departure from the q − 1 line [16,24] L ∈ [20,24] [16,24] L ∈ [20,24] 0 0.05 0.10 16,24] L ∈ [20,24] 0 0.05 0.10 16,24] L ∈ [20,24]  for all W < 2.6. This can be viewed as another indication that q * is dominated by outliers and is likely to flow to infinity for larger system sizes.
Another advantage of τ typ q is that unlike τ avg q it doesn't diverge for q < −0.5, and thus allows to study the behavior of the small values of the eigenstate coefficients. We recall that these values are of particular interest given their abundance compared to the Gaussian distribution (see Fig. 1) . In Fig. 4 we repeat the analysis done in Fig. 3 for q < 0 (for technical reasons we use a different set of data here, which includes less samples).
We estimate the value of τ typ q for GOE eigenstates based on the behavior of f GOE avg (α) for α > 1, which can be calculated analytically based on (2) and (3), Since the typical f typ (α) is determined by histogram counts growing with the system size [3], it coincides with f avg (α) for f avg (α) ≥ 0 and tends to −∞ (zero counts) otherwise. Thus, we can evaluate τ typ q from f GOE av (α) using a truncated Legendre transform [3], which is designated in Fig. 4 by the dashed black lines. We note that similarly to q > 0, τ typ q for q < 0 appears to flow to the predictions of GOE.
This conclusion is in apparent contradiction to the behavior of the distributions in Fig. 1, where an excess of zeros of the eigenstates compared to GOE prediction is clearly visible for W ∼ 1, and does not appear to vanish in the N → ∞. The discrepancy must follow from the prefactor in the definition of τ typ q , I N eigenstates from a Gaussian probability distribution in Eq. (2), while fixing the norm of the eigenstate. This procedure is very efficient and allows us to study the same Hilbert space dimensions as we do for the XXZ model, at a negligible computational cost. The results of the evaluation of I N q typ / I N ,GOE q typ can be seen in Fig. 5. For W < 1 we indeed see that for the system sizes we have the ratio flows away from 1. We strongly suspect that this apparent divergence from GOE, is a finite size effect, which has to do with the proximity of an integrable point (for W = 0 the XXZ model is integrable). We leave the examination of this effect to future studies. For W > 1 we see an apparent flow towards 1, with clearest evidence for W = 2.2. For W = 2.6 and 3.0, the finite size behavior is non-monotonic (highlighting the importance of the use of large system sizes), and appears to flow towards 1, though for these disorder strengths it is less apparent.
To summarize this section, we have seen that while a naïve examination of the distributions of the eigenstates elements in Fig. 1, shows apparent convergence to a non-Gaussian, and thus multifractal distribution, a more detailed analysis looking on the moments of the distribution, shows a slow but clear flow towards the predictions of GOE, in τ typ q , τ avg q and even directly in the finite size generalized IPRs compared to thier GOE values I N q typ / I N ,GOE q typ . In the next section we will complement this analysis, by examining an additional multifractal measure -the multifractal spectrum.

C. Multifractal spectrum
In this section we analyze the multifractal spectrum, f (α), of the eigenstates of (1), which appeared in (3), but we repeat it here for convenience, with α ≡ − ln x 2 / ln N and x ≡ | n|β |, where |n is a basis state, and |β are eigenstates of (1) [3]. The multifractal spectrum, f (α) is the fractal dimension of the set x ∈ N −α/2 , N −(α+dα)/2 , namely the probability for x to be in   Figure 6. Odd columns. f (α, N ) vs 1/ ln N for different α values (see legend). Dashed colored lines correspond to an extrapolation to N → ∞, and the stars on the y-axis correspond to the infinite size GOE prediction. Even columns. The finite size multifractal spectrum, f (α, N ) calculated using Eq. (15) for L = 12, 16, 20, 24 (darker colors, correspond to larger system sizes) and disorder strengths W = 1.0, 1.8, 2.2, 2.6, the width of the lines correspond to statistical errors. Red dashed line corresponds to the infinite size GOE prediction, fGOE (α) according to (14), and black dashed lines show the upper bounds on f (α) according to (13). this interval is given by, Using (10) and the relation (5) to I q one can see that, namely in the limit N → ∞, τ q and f (α) are related via a Legendre transform [3]. Here inf is an infimum. We note however that while τ q is a concave function the definition (10) above allows for f (α) to be non-concave a feature which we will utilize in our analysis below. Similarly to the restrictions on τ q , described above Eq. (5), from normalization of the probability distribution, P N (x) dx = 1 and the wavefunction x 2 P N (x) dx = N −1 in the limit N → ∞ one can derive, that For the Gaussian distribution of GOE eigenstates (2), using the definition (8), one can obtain the finite size correction, which in the limit N → ∞ gives the already mentioned result (8). Here A is a normalization constant being a slow (at most logarithmic) function of N . Note that this multifractal spectrum differs from f GOE (1) = 1, f GOE (α = 1) → −∞ in Ref. [3] because the latter is written for wavefunction envelopes, while the raw numerical eigenstates contain de Broglie-like oscillations corresponding to the increased statistics of zeros (large values of α > 1). While there are methods to remove these superfluous zeros (see, e.g., [44,45]), since the statistics of the zeros does not affect q > 0 moments of the eigenfunctions we don't consider such methods in this work.
To obtain the finite size multifractal spectrum we numerically compute a histogram of α with 0 ≤ α ≤ 4 (we used 50, 100, and 200 bins, and verified that our results don't change with respect to the bin number, not shown), then using (11) yields, which is presented in the even columns of Fig. 6, for various disorder strengths 1 ≤ W ≤ 2.6. Here we assumed A to be a constant (not a logarithmic function of N ) like in the GOE case as we focus on the small disorder amplitudes.
In the odd columns of Fig. 6 we extrapolate the data to N → ∞, with the same procedure used in several random matrix models (see, e.g., [44,45]). For sufficiently strong disorder or α sufficiently far from 1, our finite size data shows a nonlinear behavior in 1/ ln N (like in GOE case (8)) indicating the importance of using large system sizes in the the determination of the asymptotic behavior. Even with the state-of-the-art system sizes we use here, the extrapolation procedure is not justified due to nonlinearity of the data for some values of α. Nevertheless, for W < 2, the extrapolation works fairly well, and similarly to the moments analysis in the previous section, supports a flow towards the predictions of GOE. The extrapolation is not entirely satisfactory for W = 2.2 and 2.6, thus we cannot rule out multifractal behavior in this case. We also note that in the calculation of the histograms of α, the sampling for our system sizes starts yielding zero counts (for all used bin sizes) for α 2.5 and the majority of eigenvectors. In this interesting regime (corresponding to the excess of zeros in the wavefunction histogram, Fig. 1), the (low probability) contribution to the histogram seems to stem from the distribution over disorder realizations, rather than from representative eigenstates. This leads to large fluctuations, also visible in the errorbars (shaded area).  Figure 7. Upper row. The finite size multifractal spectrum f (α, N ) for small values of α and L = 12, 16, 20, 24 (darker colors, correspond to larger system sizes). Dashed black lines indicate the upper bounds on f (α) according to (13). The insets show f (α, N ) (solid black lines) for the maximal available size (L = 24) together with its multifractally symmetric counterpart according to Eq. (16) (dashed red lines). Lower row. Same as the upper row, but for the tilted multifractal spectrum, f (α, N ) − qα, and q = 1.5.
The finite size multifractal spectrum f (α, N ) is also useful to understand the deviation from the (q − 1) GOE line we observe for both τ avg q and τ typ q for some value q * (W, L) (see Figs. 2 and 3). By looking on the direction of the flow of f (α, N ) with the system size in Fig. 6 for all W ≤ 2.2 and α close to the maximum, corresponding to q < q * , f (α, N ) increases with N , while for small α, corresponding to q > q * , the spectrum f (α, N ) decreases with N . Moreover the crossing points of f (α, N ) at two adjacent N values increases towards α = 1 with increasing N . For W = 2.6 the situation is drastically different, since there is no downward finite-size flow at small α, but instead f (α, N ) appears to saturate there. This is also visible in the extrapolation curves for small α = 0.1 in Fig. 6. To emphasize this point in Fig. 7 we plot f (α) − qα, the supremum of which gives −τ q (see (9) for example). For W = 2.6 the left local maximum at α 0.1 does not appear to flow with system size, while the right local maximum at α α 0 drifts upward. Nevertheless this upward flow is bounded from above by the normalization conditions (13) (shown by black dashed lines in the figure) and thus the right local maximum at α α 0 cannot overcome the one at α 0.1 even in thermodynamic limit. While it is possible that there is a very slow downward flow of the left local maximum, which will eventually restore GOE, we do not see it within the available system sizes. Further support for possible multifractality at W = 2.6 can be obtained by examining the well-known symmetry of multifractal spectrum, which can be analytically derived for wavefunction envelopes in multifractal states of various models [3], In the insets of Fig. 7 we test this symmetry for the maximal available system size L = 24. To suppress the effects of zeros of the eigenstates we only examine the symmetry in the regime where the tail of f (α) is significantly above its ergodic value (3 − α) /2 (see Eq. (8)), which for our data occurs for W ≥ 2.6 (see Fig. 6). In this range of disorder strengths the multifractal symmetry (16) is satisfied (insets of Fig. 7), while in the complementary range, W ≤ 2.2 the symmetry doesn't apply. This indicates a possible multifractal phase for W > 2.2.

IV. SUMMARY AND DISCUSSION
In this work we have conducted a detailed large-scale numerical study of multifractal properties of eigenstates on the delocalized side of the many-body localization transition (MBL). This phase is known to have a number of anomalous dynamical features, such as subdiffusive transport, sublinear entanglement entropy growth and suppressed spreading of information [23]. For the single-particle case, suppressed relaxation and dynamics are often associated with spatial sparseness of the underlying eigenstates [31,32]. A natural question to ask is whether a similar relation exists also in the many-body case, namely if sparseness of the eigenstates in Hilbert space implies slow relaxation and suppressed transport of local observables. In this work we answer this question in the negative, by identifying a large fraction of the delocalized phase which is consistent with ergodicity, while still showing a clear signature of subdiffusion and slow relaxation in both numerical and experimental data. We reach this conclusion by a careful analysis of the finite size flows of eigenstate coefficient distributions, moments of these distributions and their spectrum of multifractal dimensions. Our analysis focuses on the computational basis, where the basis states are labeled by the eigenvalues of the localŜ z i operators. This is the natural basis for the XXZ chain in the context of MBL since it is compatible with the disorder and is naturally linked to the hopping problem in Hilbert space. To the best of our knowledge, the multifractal spectrum of the disordered XXZ chain has not been studied before due to severe finite-size behavior and non-monotonic behavior (see Figs. 5 and 6 for example), which hindered reliable extrapolation to the thermodynamic limit. All measures we study give a coherent picture of a steady flow towards the predictions of GOE for disorder strengths W 2.6. A slight discrepancy compared to the GOE prediction, is observed for W < 1, as an excess of small values of the eigenstate coefficients compared to GOE. Although it is consistent with a so-called weak ergodic phase, where the wavefunction occupies a finite, but tiny fraction of the Hilbert space, observed in several single-particle and many-body systems [69][70][71] , we argue that this discrepancy is a result of proximity to an integrable point at W = 0, and should -if this is indeed the case -disappear either in the thermodynamic limit, or if integrability is broken. We leave the verification of this prediction to a future study. For larger disorder strengths, our analysis becomes unreliable, due to slowing down of finite-size flows. While we cannot rule out a slow residual flow to GOE, we don't observe it within our range of accessible system sizes. Moreover, for W = 2.6 the multifractal spectrum perfectly satisfies one of its basic symmetries, which would be consistent with multifractality of the eigenstates in this region. Given the immense numerical cost of our calculations, we could only compute the spectrum in a limited range of disorder strengths across the delocalized phase. Combined with the slow finite size flows at stronger disorder, we cannot determine whether the region consistent with multifractality shrinks to the critical point when the system size is increased, as was claimed in Ref. [52]. It would be interesting to study this important question in more detail in the future.
One of the central outcomes of our study suggests that the previously observed anomalous dynamics is not related to multifractality of many-body eigenstates. Since multifractal features are generically basis dependent, one can wonder whether the outcome of our study changes with the change of the basis. For GOE eigenstates the multifractal features do not depend on the basis, since the GOE distribution is invariant under orthogonal transformation [62]. However, this is only true for most bases (for example the eigenstate basis is clearly not a good basis to study multifractality).
The family of bases generated by locally exciting the eigenstates of the system, and which can be directly tied to relaxation of local observables, is a good candidate for such purpose [52].

ACKNOWLEDGMENTS
The authors acknowledge fruitful discussions with Nicolas Macé, Fabien Alet and Nicolas Laflorencie. This research was supported by the Israel Science Foundation (grant No. 527/19). We acknowledge PRACE for awarding access to HLRS's Hazel Hen computer based in Stuttgart, Germany under grant number 2016153659. Our code is based on the Petsc [72,73], Slepc [74,75] and Strumpack [60,61]  Note that for the diagonal matrix elements α|Ŝ z i |α the distribution for each disorder realization has a (slightly) nonzero mean, which we subtracted here (cf. discussion in Ref. [76] In this appendix, we turn our attention to the analysis of the distributions of matrix elements of the local magne-tizationŜ z i in the eigenbasis of the Hamiltonian. Similarly to our analysis of the eigenstates coefficients in the main text, for each disorder realization we consider 50 eigenstates |α of the HamiltonianĤ with an eigenvalue E α closest to middle of the many-body spectrum (E max − E min ) /2. These eigenstates correspond roughly to infinite temperature. We note however, that since we study the microcanonical ensemble with S tot z = 0 here, where Tr Sz=0 H = 0, in each disorder realization there is a slightly different "effective temperature", which we correct by subtracting the mean of the diagonal matrix elements (computed over the extracted eigenstates) for each disorder realization (cf. discussion in Ref. [76] and in particular Appendix B therein).
We complement our previous work in Refs. [22,77], by calculating the distribution of the matrix elements α|Ŝ z i |β of the localŜ z i in the eigenbasis {|α } of the Hamiltonian, with a massively improved statistics and one additional system size (L = 24). We also add logarithmic binning of the histograms, a direct distribution of diagonal matrix elements rather than their differences as well as a direct comparison of diagonal α = β and offdiagonal α = β matrix elements distributions. Fig. 8 shows the results for the logarithmically binned probability density of diagonal α|Ŝ z i |α and offdiagonal α|Ŝ z i |α | matrix elements ofŜ z i for disorder strengths W = 1.0, 1.8, 2.2, 2.6, which are well on the delocalized side of the phase diagram, for system sizes L = 10, 12 . . . , 22, 24. The logarithmic binning highlights the maximum of the distributions, where the matrix elements are closeset to zero. At weak disorder W 1.0 we note that both diagonal and offdiagonal matrix elements assume a distribution very close to Gaussian as predicted by ETH [78]. Furthermore, the prediction from random matrix theory that distributions of diagonal and offdiagonal matrix elements should be directly related [79,80] is verified to very high precision (dashed red lines in the lower panels of Fig. 8 are diagonal distributions for L = 22, 24 (renormalized by √ 2 in order to take account of convolution of two gaussian distributions) in comparison to offdiagonal distributions shown in color). It is clear however (as was shown in Ref. [80]), that for stronger disorder W 1.8, this correspondence is violated.
For the diagonal matrix elements the shape of the maximum appears to be Gaussian (flat on a logarithmic scale), however the tails of the distribution deviate from Gaussian distributions at disorder strengths W 1.0. The double logarithmic scale reveals a long straight tail, particularly well developed for W = 1.8 and W = 2.2, which seems to be consistent with a power law tail over more then one decade.For the offdiagonal matrix elements the tail seems to decay faster than a power law. In contrast to the diagonal matrix elements, there is a significant excess weight at small values of the offdiagonal matrix elements α|Ŝ z i |β , which seems to scale to zero for W = 2.2 but survives up to at least L = 24 for W = 2.6. The scaling of the matrix element distribution variance inversely with Hilbert space dimension is well visible at weak disorder W = 1.0 in the (almost) equidistant distributions for both diagonal and offdiagonal matrix elements and was analyzed in detail in Refs. [22,77]. Increasingly strong deviations from this scaling are observed at stronger disorder, which were connected to subdiffusive transport [22].