Multifractality and self-averaging at the many-body localization transition

Finite-size effects have been a major and justifiable source of concern for studies of many-body localization, and several works have been dedicated to the subject. In this paper, however, we discuss yet another crucial problem that has received much less attention, that of the lack of self-averaging and the consequent danger of reducing the number of random realizations as the system size increases. By taking this into account and considering ensembles with a large number of samples for all system sizes analyzed, we find that the generalized dimensions of the eigenstates of the disordered Heisenberg spin-1/2 chain close to the transition point to localization are described remarkably well by an exact analytical expression derived for the non-interacting Fibonacci lattice, thus providing an additional tool for studies of many-body localization.

An eigenstate is multifractal when it is extended, but covers only a finite fraction of the available physical space. Multifractality is characterized by the so-called generalized dimensions D q , for fully delocalized states D q = 1, for multifractal states 1 < D q < 0, and for localized states D q = 0. In the thermodynamic limit, all eigenstates of one-dimensional (1D) noninteracting systems with uncorrelated random onsite disorder are exponentially localized in configuration space for any disorder strength. It is at higher dimensions that the delocalization-localization transition takes place and this happens at a single critical point, where the eigenstates are multifractal. In contrast, if interactions are added to these systems, the delocalization-localization transition happens already in 1D and for finite disorder strengths, fractality exists even in the MBL phase. For these interacting systems, it is still under debate whether before the MBL phase there is a single critical point or an extended phase where the eigenstates are multifractal [24,[27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42]. In fact, even the very existence of the MBL phase has now gone under debate [43][44][45]. One of the reasons why it is so hard to settle these disputes is the presence of serious finite-size effects. Recent large-scale numerical studies [39,41] of the disordered spin-1/2 Heisenberg chain, where Hilbert space dimensions of sizes ∼ 3 × 10 6 have been reached, did not question the transition to a localized phase, but were not entirely conclusive with respect to the existence of an extended nonergodic phase or a single critical point, although the latter is strongly advocated in Ref. [39].
In this work, we consider the same Heisenberg model and emphasize another problem that has not received as much attention as finite-size effects, but is also crucial for studies of disordered systems, that of lack of self-averaging. This issue becomes particularly alarming as the system approaches the transition to the MBL phase [38,46,47]. If a quantity is non-self-averaging, the number of samples used in statistical analysis cannot be reduced as the system size increases [46,[48][49][50][51][52][53][54][55][56][57][58]. This reduction is a very common procedure due to the limited computational resources when dealing with exponentially large Hilbert spaces, but it may lead to wrong results. We show that when the disorder strength of the spin model gets larger than the interaction strength and it moves away from the strong chaotic (thermal) regime, the fluctuations of the moments of the energy eigenstates increase as the system size grows, exhibiting strong lack of self-averaging. Decreasing the number of random realizations in this case may affect the analysis of the structures of the eigenstates, including the results for the generalized dimensions.
The various challenges faced by the numerical studies of the MBL is a great motivator for theoretical works, which, however, have difficulties of their own. The current trend is to focus on phenomenological renormalization group approaches [59][60][61][62][63][64][65][66] that aim at improving our understanding of the MBL transition in 1D systems with quenched randomness, without providing microscopic details. Some of these studies suggest that the transition is characterized by a finite jump of the inverse localization length. Similarly, numerical studies indicate that the generalized dimensions jump at the critical point [39], and a connection between these two jumps was proposed in [47].
Our contribution to those theoretical efforts is to show that an exact analytical expression for the generalized dimensions derived for the 1D non-interacting Fibonacci lattice [67,68] matches surprisingly well our numerical results for the disordered spin-1/2 Heisenberg chain in the vicinity of the MBL critical point. This expression provides an additional tool in the construction of effective models for the MBL transition. Its derivation is based on a renormalization group map of the transfer matrices used to investigate the wave functions of the Fibonacci model [67,68]. Our 1D lattice system has L interacting spin-1/2 particles subjected to on-site magnetic fields. It is described by the Hamiltonian where S x,y,z k are spin-1/2 operators, the coupling strength was set equal to 1, h k are random numbers from a flat distribution in [−h, h], h being the disorder strength, and periodic boundary conditions, S x,y,z L+1 = S x,y,z 1 , are imposed. Since H (1) conserves the total spin in the z-direction, S z = S z k , we work in the largest subspace corresponding to S z = 0, which has dimension N = L!/(L/2)! 2 . The model is integrable when h = 0 and chaotic, that is, it shows level statistics similar to those from full random matrices [69], when h chaos ≤ h < h c . The value of h chaos for the transition from integrability to chaos and of the critical point h c for the transition from delocalization to the MBL phase are not yet known exactly. Our focus here is on the second transition, and for that, some works estimate that 3 < h c < 4 [10,29,[70][71][72] and others that h c > 4 [73,74].
Multifractality and ensemble size.-To obtain the generalized dimensions D q , we perform scaling analysis of the generalized inverse participation ratios, which are defined as IPR α q := N k=1 | φ k |ψ α | 2q , where q can take, in principle, any real value, |ψ α is an eigenstate of the Hamiltonian (1), and |φ k represents a physically relevant basis. Since we study localization in the configuration space, |φ k is a state where the spins point up or down in the z-direction, such as | ↑↓↑↓ . . . . We average the generalized inverse participation ratios, IPR q , over ensembles with n samples that include 0.02N eigenstates with energy close to the middle of the spectrum and n/(0.02N ) random realizations, and then extract the generalized dimensions using Multifractality holds when D q is a nonlinear function of q.
In practice, D q is obtained from the slope of the linear fit of ln IPR q versus ln N . In Fig. 1, we show some representative examples of the scaling analysis for different values of h and q, and also for ensembles of different sizes n, varying from n = 10 2 to n = 3 × 10 4 . The symbols are numerical data and the solid lines are the corresponding fitting curves.
In the chaotic region, for example when h = 1, the scaling of ln IPR q with ln N is independent of the size of the ensemble, with all points and lines for a given L coinciding and leading to D q ∼ 1. This is shown in Fig. 1 (a) for q = 1.2 and it holds for all other values of q that we studied, 0.1 ≤ q ≤ 3. In contrast, when h > 1, the numerical points strongly depend on the number of samples used, as seen from Fig. 1 Fig. 1 (f). Notice that this dependence becomes more evident for the larger system sizes. In the particular cases of Figs. 1 (b)-(e), where 1 < h < ∼ h c , the fittings lead to larger slopes when the ensemble sizes are smaller. For these smaller n's, the values of D q would get even larger if we would neglect the smallest system sizes when doing the fittings. These results illustrate the danger of reducing the number of samples as the system size increases.
We verify in Figs. 1 (b)-(f) that the convergence of our numerical points happens for ensembles with n > ∼ 2×10 4 . Indeed the points for n = 2×10 4 and n = 3×10 4 are nearly indistinguishable, so in all of our subsequent studies, we use n = 3×10 4 for all L's. It may be, however, that for system sizes larger than the ones considered here, convergence would require even larger ensembles.
Self-averaging.-The fluctuations of the values of IPR q bring us to the discussion of self-averaging. A given quantity O is self-averaging when its relative variance decreases as the system size increases [48][49][50][51][52][53][54][55][56]. This implies that in the thermodynamic limit, the result for a single sample agrees with the average over the whole ensemble of samples.
In quantum many-body systems, the eigenstates can spread over the many-body Hilbert space, which is exponentially large in L, so we study the scaling of R IPRq with N [46,57], If ν < 0, then IPR q is self-averaging and one can reduce the number of samples for the average as the system size increases. This cannot be done when ν ∼ 0, and it is even worse in the extreme scenario where ν > 0 and the relative fluctuations increase as the system size grows. In the main panels of Fig. 2, we show the dependence of R IPRq on the disorder strength h for different values of q and each line represents one system size. It is clear that deep in the chaotic region, the relative variance decreases as the system size grows, implying self-averaging of IPR q . This is also illustrated in the insets, where ν < 0 for h chaos ≤ h < ∼ 1, which is consistent with Fig. 1 (a), where the scaling of IPR q does not depend on the number of samples.
There is, however, a turning point at h > ∼ 1, where ν suddenly jumps above zero and R IPRq grows significantly with system size. As seen in Figs. 2 (a)-(d), this is particularly bad in the region preceding the MBL phase, For h > ∼ 4, where the system should already be in the MBL phase, the relative variance R IPRq continues to grow with system size, but ν is close to zero and the curves for L = 16 and L = 18 are not far from each other.
The results in Fig. 1 and Fig. 2 show that extra care needs to be taken when performing scaling analysis away from the chaotic region, not only due to finite-size effects, but also due to the lack of self-averaging. No matter how large the system size is, large numbers of samples are required and may even need to be increased as L grows.
One can reduce the fluctuations of the generalized inverse participation ratios by using their logarithm, known as participation Rényi entropies. In fact, using a toy model, it was shown in [46] that in the MBL phase, R IPR2 grows with system size, while R − ln IPR2 decreases with L. However, for 1 < ∼ h < ∼ 4, even though we observe a reduction of the fluctuations, ln IPR q remains non-selfaveraging and we still have ν ∼ 0 [75]. We indeed verified that the plots shown in Fig. 1 remain similar if instead of ln IPR q , we use ln IPR q .
Multifractality and analytical expression for D q .-After taken the necessary precautions for performing the scaling analysis of the generalized inverse participation ratios, as discussed in Fig. 1, we now proceed with the study of how D q depends on q and h.
In Refs. [67,68], an exact analytical expression was derived for the structure of the eigenstate at the center of the spectrum of the off-diagonal version of the Fibonacci model in the thermodynamic limit, leading to the generalized dimensions where σ = ( √ 5 + 1)/2 is the golden mean and λ(h) = (h + 1) 2 + (h + 1) 4 + 4h 2 1/2 /(2h) is the maximum eigenvalue of the transfer matrix [67,68]. For the Fibonacci model, h denotes the ratio between its two hopping constants, which are arranged in a Fibonacci sequence.
In the case of our interacting spin model, we use Eq. (4) as an ansatz. Since in this case, the eigenstates are extended for h chaos ≤ h < ∼ 1, while Eq. (4) predicts a monotonic decrease of D Fibonacci q for h < 1, we compare our results with the expression for Θ(h−1)D Fibonacci q +Θ(1−h), where Θ is the Heaviside step function. We find that this expression matches the numerical values of D q for the spin chain extremely well for disorder strengths in the vicinity of the critical value, 3 < h c < 5.
The inverse participation ratio, IPR 2 , is the most commonly used quantity in studies of localization, so we start by analyzing D 2.0 and D 2.4 as a function of the disorder strength. They are shown in Fig. 3 (b) and Fig. 3 (c), respectively. The numerical results agree very well with the analytical expression for D Fibonacci q for all h's. However, for smaller q's, the agreement is not as good, as illustrated in Fig. 3 (a), the same happening for larger q's, as seen in Fig. 3 (d).
We cannot say whether the agreement of the curves for D q vs h with D Fibonacci q for q away from 2 would improve or get even worse if larger system sizes were considered. If it would improve, that would point to the existence of an extended phase of multifractal eigenstates before the MBL phase and described by the analytical expression of the Fibonacci lattice. Large-scale numerical studies [39,41] indicate that if such a phase exists, it should appear for h > 2 [41], and it may as well be a single point [39].
We stress, however, that the most relevant and less controversial information provided by Fig. 3 (a) and Fig. 3 (d) is that the numerical points for different values of q's cross the curve of the analytical expression of D Fibonacci Fig. 3 (a) and Fig. 3 (d)]. This indicates that at least in the vicinity of (or right at) the critical point, the generalized dimensions of the disordered spin model is indeed extremely well described by Eq. (4).
The bottom panels of Fig. 3 give further support for this observation. There, we plot D q as a function of q for different values of the disorder strength. For 1 < h < h c , as illustrated by Fig. 3 (e), there is no good agreement between the numerical points and D Fibonacci q . The same happens for h > h c , as seen in Fig. 3 (h), although the mismatch in this case is not as large. However, for h ∼ h c , as shown in Fig. 3 (f) and Fig. 3 (g), the agreement is extremely good.
Conclusions.-Our analysis of the disordered spin-1/2 Heisenberg chain calls attention to the strong lack of selfaveraging of the generalized inverse participation ratios for a range of disorder strengths that precedes the critical point of the MBL transition. This implies that in theoretical and experimental studies of this region, one should not decrease the number of samples as the system size increases. We notice also that the logarithm of the generalized inverse participation ratios can be used to reduce fluctuations, but it still does not lead to selfaveraging in that region.
Our studies indicate a strong relationship between multifractality, 0 < D q < 1, and the lack of self-averaging of the generalized inverse participation ratios, ν ≥ 0. Multifractality reflects the fragmentation of the Hilbert space [76], and this fragmentation, in turn, leads to the sample-to-sample fluctuations associated with the absence of self-averaging of IPR q . The latter should then hint at the existence of multifractal states.
The comparison between our numerical results for the generalized dimensions of the disordered spin chain and the analytical expression for D q derived for the offdiagonal version of the Fibonacci model shows remarkable agreement in the vicinity of the MBL transition. This connection is useful for theoretical efforts seeking to adequately describe the critical point and may serve as a reference for studies of transport behavior. It should also motivate additional numerical studies to verify whether the agreement holds in a finite region or only at a single critical point.
LFS was supported by the NSF grant No. DMR-1936006. E.J.T.-H. is grateful to LNS-BUAP for their supercomputing facility.
In this supplemental material (SM), we show that the application of a logarithmic transformation to the generalized inverse participation ratios of the energy eigenstates reduces the size of the fluctuations, but is not enough to achieve self-averaging, specially around the critical point. This implies that reducing the numbers of statistical data to compute averages as the system size increases remains a problem also in this case.
We also present results for the errors δD q involved in the computation of the generalized dimensions. We find that the errors obtained for the linear fit of ln IPR q versus ln N and those for ln IPR q versus ln N are very similar. They are larger in the vicinity of the critical point and get significantly larger for all values of h > 1 as one decreases the size of the ensembles, while in the chaotic region, h ∼ 1, we have that δD q ∼ 0 for all numbers of samples considered. This reinforces our claims that outside the chaotic region, we should not reduce the size of the ensembles as the system size increases, especially in the vicinity of the critical point.

LOGARITHM OF THE GENERALIZED PARTICIPATION RATIOS
An alternative approach to compute the generalized dimensions D q consists in using the logarithm of the generalized participation ratios, which corresponds to the so-called participation Rényi entropies. In this case, instead of doing the linear fit of ln IPR q versus ln N , as in the main text, we study That is, instead of computing the logarithm of the averaged generalized inverse participation ratios, which is an arithmetic mean, we now compute the average of the logarithm of the generalized participation ratios, which is a geometric mean. The logarithmic transformation is commonly applied to reduce the fluctuations of a set of data and it has a significant effect on the tails of the distributions.
In Fig. S1, we show ν for R − ln IPRq ∝ N ν . Comparing it with the insets of the Fig. 2 in the main text, we see two main differences. One is that in the localized phase, − ln IPR q becomes self-averaging, in agreement to what was discussed in Ref. [42]. The other difference is that the values of ν around the critical point become much smaller, but it is still not negative, so the lack of selfaveraging persists.

ERRORS
As show in Fig. S2, there is no significant difference between the error δD q obtained by extracting the generalized dimensions D q from the scaling of ln IPR q with ln N [δD The results in all four panels are very similar. In the chaotic region, h chaos ≤ h < ∼ 1, δD q is close to zero for all ensemble sizes considered. The errors increase monotonically for 1 < ∼ h < ∼ 2, but remain almost independent of the number of samples. It is in the vicinity of the critical point, 2 < ∼ h < ∼ 4.5, that the errors become clearly larger as the number of samples gets decreased. For h > ∼ 4.5, the errors still depend on the number of samples, but they are smaller that in the preceding region, specially for the ensembles with 3 × 10 4 samples for which a sudden drop is seen at h ∼ 2.8.