Chain breaking and Kosterlitz-Thouless scaling at the many-body localization transition

Despite tremendous theoretical efforts to understand subtleties of the many-body localization (MBL) transition, many questions remain open, in particular concerning its critical properties. Here we make the key observation that MBL in one dimension is accompanied by a spin freezing mechanism which causes chain breakings in the thermodynamic limit. Using analytical and numerical approaches, we show that such chain breakings directly probe the typical localization length, and that their scaling properties at the MBL transition agree with the Kosterlitz-Thouless scenario predicted by phenomenological renormalization group approaches.

In such a puzzling context, several progresses have been made to build an analytical theory, able to describe the ergodic-MBL transition [54][55][56][57][58][59][60][61][62]. The most successful description, based on the so-called "avalanche" scenario [63,64], proposes a phenomenological renormalization group (RG) treatment, working in the MBL regime where large insulating blocks compete with small ergodic inclusions. In this framework, recent works [65][66][67] found a RG flow of the Kosterlitz-Thouless (KT) form. The MBL phase is described as a line of fixed points with a vanishing density of ergodic inclusions and a finite typical localization length ζ, which controls the spatial extension of the l-bits [68][69][70]. When the delocalization transition is approached from the MBL side, the typical localization length ζ grows and reach a finite critical value ζ c = (ln 2) −1 at the transition point h c , with a singular behavior [71] As argued in Ref. [66], a diverging length scale ∼ exp(c/ √ h − h c ) should control finite-size effects, thus explaining some limitations in the numerics.
In this Letter, we show that one can overcome such finite size constraints by measuring spin polarization, a simple local observable of the S = 1/2 Hamiltonian Eq. (1). Our main finding is that the extreme statistics of the polarizations directly probes the typical length scale ζ of MBL physics. At the thermodynamic limit, the MBL regime is driven by a chain breaking mechanism, signalled by a complete freezing of some sites with | S z i | − 1/2 ∼ L −γ . The freezing exponent γ > 0, which controls the corrections, can be analytically related to the typical localization length by γ ∝ ζ −1 . Building on Refs. [36,45,49], a careful scaling analysis of our ED data reveals that they are best described by a volumic scaling variable N /Λ in the delocalized regime, where N is the Hilbert-space size and Λ a non-ergodicty volume, while a logarithmic scaling variable (ln L)/λ dictates the behavior in the MBL regime, with λ −1 ∝ ζ −1 − ζ −1 c following the KT singularity Eq. (2). arXiv:2004.02861v2 [cond-mat.dis-nn] 7 Apr 2020 Distribution of local polarizations and extreme value statistics-Let us start the discussion by looking at the local polarizations m z = S z , computed by ED simulations performed at infinite temperature for the Hamiltonian Eq. (1). As already noticed in Refs. [72][73][74][75], the histograms P (m z ) display distinct features across the two regimes, as illustrated in Fig. 1 (a) for 3 representative values of the disorder strength (h = 1, 3, 7). At weak disorder, we expect from the eigenstate thermalization hypothesis (ETH) [76][77][78] to observe Gaussian distributions, peaked around m z = 0 and shrinking with increasing system size, as clearly visible for h = 1. At h = 3, a more complex form emerges with deviations from gaussianity [79,80]: strongly polarized (m z ±1/2) sites appear, but their density shrinks down with system size (as evidenced by the decrease of magnetization variance, see inset). At strong disorder h = 7, the density of strongly polarized sites no longer shrinks down: ETH is violated and the distribution is U-shaped, almost free of finite-size effects.
To quantify this effect, we introduce the deviation from perfect polarization δ = 1 2 − |m z |. While δ > 0 for finite systems, δ becomes arbitrarily small for large disorder, as evidenced in Fig. 1 (a) at h = 7 where one observes [81] with an exponent γ ≥ 0 related to extreme value statistics [82], as we discuss now. In each finite sample, we define for the most polarized site whose distributions are shown at h = 3 and h = 7 in Fig. 1 (b). There, we explicitly report two distinct trends with increasing system sizes L. In the ETH regime (h = 3), some weight is transferred with increasing L towards small − ln δ min , with a peak developing at − ln 2. The opposite is observed in the MBL regime (h = 7), where weight is moved towards large values of − ln δ min . In both cases we qualitatively spot that rare events of the competing phase (rare insulating bottenecks for the delocalized phase vs. thermal bubbles in the MBL regime) become less and less relevant upon increasing L.
Let us now turn to a more quantitative analysis of the extreme polarization finite size scaling. At strong disorder, assuming that the magnetizations along the chain are independently drawn from the distribution P (m z ) of Fig. 1 (a), the deviation from perfect polarization verifies δmin 0 P (x)dx ∼ 1/L [82], which for a power-law tail Eq. (3), yields Such a scaling has dramatic consequences since for any γ = 0, we expect δ min (L) → 0 when L → ∞, meaning spin freezing, and therefore chain breaking at the thermodynamic limit.
Extreme polarizations: numerical results-In order to numerically check the power-law scaling Eq. (S2), and probe the chain breaking mechanism at the microscopic level, the most polarized site is recorded for each finite-length sample. Typical values are computed from shift-invert ED simulations at infinite temperature, and disorder averaging is performed over a large number of realizations, at least 10 3 samples. Results as a function of chain lengths L are shown in Fig. 2 (a) for a wide range of disorder strengths. As expected from the extreme value argument Eq. (S2), at strong disorder we observe a power-law decay with L, indicating a chain breaking at the thermodynamic limit. In contrast, for weak disorder δ min does not vanish with L, but instead tends to 1/2, showing "healing" in a way similar to Kane-Fisher physics [83,84]. Such radically different behaviors correspond to the two phases of the model, as we argue below. In Fig. 2 (b), we show the h-dependence of the freezing exponent γ obtained from power-law fits to the form   Eq. (S2) over 4 different fitting windows. Finite-size scaling towards perfect polarization is governed by a disorder-dependent freezing exponent γ, which shows a logarithmic divergence at large h, a behavior explained below by analytical arguments. Conversely, at weak disorder we observe a finite-size (possibly non-monotonous) crossover, and ultimately γ → 0 as expected from the ETH. Before making a quantitative scaling analysis of ED data (Fig. 3), we first provide an analytical understanding of the chain breaking mechanism.
Analytical derivation at large disorder-Using the Jordan-Wigner transformation, we can rewrite Eq. (1) as interacting spinless fermions in a random potential The first sum, describing free fermions, can be diagonal- The single-particle orbitals φ k i are exponentially localized for any h = 0, and the interacting Hamiltonian rewrites in this "Anderson basis" as 4 terms: In the strong disorder limit h 1, the second line can be neglected [81]. The number operators n k then form a complete basis of local operators commuting with H and among themselves: in other words, the Hamiltonian is written in l-bit form [55,68,70,[85][86][87], and the n k are good approximation of the l-bits in the strong disorder limit [88]. The fermion density and thus the spin polarization at real-space position i is then simply the sum of contributions coming from occupied orbitals: Building on previous ideas [89], we are now ready to understand the spin freezing mechanism. If we approximate the Anderson orbitals by simple exponential func- is the localization center of orbital k, and ξ ∼ (ln h) −1 the localization length at strong disorder [89], we expect the maximal density n i max to occur in the middle of the longest region of max consecutive sites that are occupied by an orbital. At half-filling, a configuration with consecutive occupied sites appears with probability 2 − , which for a finite chain of length L 1 yields max ≈ ln L/ ln 2. The most polarized site i 0 corresponds to the site with maximal n i0 ≈ 1 (resp. minimal n i0 ≈ 0) fermionic density, yielding S z i0 ≈ +1/2 (resp. S z i0 ≈ −1/2). As a result, the sum of exponentials Eq. (8) gives δ min ∼ exp(− max /(2ξ)), which naturally leads to the power-law decay Eq. (S2) with a freezing exponent At large disorder, numerical data perfectly agree with a logarithmic growth γ(h) ∝ ln h (panel (b) of Fig. 2), which validates our analytical description of the freezing mechanism. Within this description, 1/γ is identified with the localization length of l-bits deep in the bulk of the largest non-thermal region. Being far away from rare thermal inclusions, the inverse freezing exponent provides an estimate of the typical l-bit extension, i.e. the typical localization length ζ: 1/γ ∼ ζ. As we further elaborate below, this has decisive consequences for our understanding of the critical behavior. Already in Fig. 2 (b), when the transition is approached from the MBL side, a singular behavior develops for γ(h) close to h c , in agreement with Eq. (2), followed by a jump to zero in the ergodic phase. We now provide an explicit description of this critical behavior.
Scaling analysis: KT behavior-We have performed a very careful analysis of our ED data for δ typ min (L, h) in order to address the transition between ergodic and MBL phases. The best data scalings [81], shown in Fig. 3 (a), are obtained using two distinct scaling functions: In the delocalized phase (h < h c ), we obtain a volumic scaling, in agreement with our recent multifractality analysis [36] (see also [45] for the Anderson transition on random graphs), and with ETH predictions [76], given that the density of state scales with Hilbert space           size. The scaling variable is the ratio between the Hilbert space size N ≈ 2 L / √ L and a disorder-dependent nonergodicity volume Λ, which diverges exponentially at criticality: ln Λ ∼ (h c − h) −ν d with ν d ≈ 0.5, see Fig. 3

(d).
Here the fit of the scaling function to the data is equally good for h c = 3.8 and h c = 4.2, as visible in Fig. 3 (c) where the chi-squared statistic divided by the number of degrees of freedom of the fit χ 2 /N df is minimum [81].
In contrast, the MBL regime (h > h c ) is best described by the function g ln L λ displayed in Fig. 3 (a). This is a direct consequence (see [36,49,81]) of the power-law decay Eq. (S2), observed at criticality and in the MBL phase, which yields an MBL scaling function g ∝ (ln L)/λ for large ln L λ. As an outcome, the scale λ is directly related to the freezing exponent, and therefore with the typical localization length, such that The larger h c ≥ 4.5, the better the goodness of fit for the MBL scaling, see Fig. 3 (c). Indeed, by definition the logarithmic scaling perfectly describes data with algebraic behavior such as δ typ min in the MBL phase. To correctly estimate h c , we therefore need to sum the goodness of fit from both ergodic and MBL regimes, see Fig. 3 (c). Then a bootstrap analysis of the data [81] gives a critical disorder strength h c = 4.2(5). This quite large relative error of ≈ 10% reflects the difficulties inherent to this ergodic-MBL transition.
However, this very number is not decisive for our theoretical understanding of the critical point [81]. As shown for h c = 4.2 in Fig. 3 (e), we get the following singularity for the inverse typical localization length: with ν loc = 0.52(3) after a bootstrap analysis, and γ c = 1.7(1) (see panel (b) of Fig. 3). This behavior quantitatively agrees with the recently proposed [66,67] KT mechanism, see Eq.
(2). Moreover, the scaling variable (ln L)/λ implies that finite size effects in real-space are formally controlled by an exponentially diverging length scale exp(λ) ∼ exp(h − h c ) −ν loc , thus confirming a KT scenario, compatible with the Harris bound.
Consequences and discussion-Our key-result is that the extreme value statistics in the MBL regime gives a direct access to the typical localization length ζ, and therefore the typical l-bits extension of the MBL states. Maximally polarized sites correspond to entanglement bottlenecks: nearly frozen spins are almost disentangled from the rest of the system, with an entanglement entropy cutting such sites given by S = −δ ln δ to leading order [92]. Consequently, δ controls the leading eigenvalues of the entanglement spectrum, whose observed power-law distribution [93,94] is also governed by the freezing exponent γ.
Our analysis further indicates that the entire MBL regime, including the critical point, witnesses a chain breaking mechanism at the thermodynamic limit, with finite size effects controlled by a power-law behavior Eq. (S2). The power-law exponent, related to the typical localization length γ ∝ ζ −1 , diverges logarithmicaly at strong disorder, and displays a jump at the transition, with a singular critical behavior Eq. (12) in perfect agreement with a KT mechanism Eq. (2).
Besides probing the transition universality class, spin freezing has deep consequences for MBL physics. Firstly, it provides a simple picture accounting for the absence of thermalization at the thermodynamic limit. Another important aspect concerns the recently discussed Hilbertspace fragmentation [95][96][97][98] of the MBL regime, which here is expected to naturally emerge from spin freezing upon increasing system sizes. In addition, our results are very encouraging for the development of a perturbative decimation method [99] discarding the strongly polarized sites. Coupled to exact methods this could provide quantitative results at strong disorder for system sizes much larger than currently accessible, and thus potentially useful beyond one dimension. Supplemental Material for "Chain breaking and Kosterlitz-Thouless scaling at the many-body localization transition" In this supplemental material we provide additional informations regarding: 1. The power-law divergence of the distributions P ( 1 2 − |m z |), in particular finite-size effects. 2. The analytical derivation in connection with free-fermions. 3. The finite-size scaling analysis of the MBL transition.

S1. POWER-LAW DIVERGENCE OF THE DISTRIBUTIONS: FINITE-SIZE EFFECTS
In the MBL regime, the deviation form perfect polarization δ = 1 2 − |m z | displays power-law tails in the distributions, according to P (δ) ∝ δ −1+ 1 γ (δ → 0), which corresponds to an exponential tail for ln δ: (S1) Histograms of ED data for L = 8, 10, 12, . . . , 22 are shown in Fig. S1 (a) in the MBL regime (h = 7). In the left panel, P (δ) is displayed in log-log scale, and a power-law behavior is visible, but does not correspond to the extreme value tail which we aim at exploring in this work. In order to better see such a tail, we show P (ln δ) in the right panel where an exponential tail is clearly visible, thus targeting much smaller values of δ, and thus the extreme values. However, the apparent prefactor 1/γ extracted from the exponential tail shows important finitesize effects, as reported in the inset of Fig. S1 (a) where we have compared this estimate with the one extracted from the extreme value scaling which has almost no finite-size effects. Considering the limited system sizes available to shift-invert ED techniques for the interacting problem, we have also explored this effect for the non-interacting case of many-body freefermions, described by the random-field XX model in the S z tot = 0 sector. Using standard free-fermion techniques, one can easily access to large systems, typically L ∼ 10 3 without particular numerical effort. Results are shown in Fig. S1 (b) for the very same system sizes as for the MBL case, with the following additional lengths L = 24,32,48,64,96,128,192,256,384,512,1024. Disorder strength h = 4 ensures that non-interacting fermions have short localization lengths. Interestingly, we clearly observe a slow convergence of the tail exponent 1/γ, as compared to the extreme value Eq. (S2) which, like for the MBL case, shows almost no finite-size effects. Therefore, we clearly see that focusing on the extreme value δ min (L) offers the most direct way to get the asymptotic value γ, as compared to the distribution tail. The spin-1 2 XXZ Hamiltonian can be decomposed in two terms: H 0 is the random-field XX model, Eq. (S3), which corresponds to free-fermions (using the Jordan-Wigner transformation) in a random potential The interaction terms read where C is an irrelevant constant. The non-interacting part Eq. (S5) is diagonalized by b k = L i=1 φ k i c i (φ k being the single particle orbitals, Anderson localized for any finite disorder in 1D), yielding In this Anderson basis, the interaction term Eq. (S6) is where we have assumed real orbitals. As proposed in the main text, we decompose the interacting part in 4 terms: The first two terms V (1,2) are diagonal, and can be interpreted as a first approximation for the so-called l-bit Hamiltonian, while V (3,4) are off-diagonal.

s2. Strong disorder limit
Below, we show that off-diagonal terms can be ignored at strong disorder. In such a limit, the Hamiltonian remains diagonal in the Anderson basis, i.e.
k,l n k n l . (S9) In this case, the maximum of the fermionic density can be readily obtained. We first use a simplified expression for the non-interacting orbitals, assuming the following exponential form The particle density at a given site i given by is maximal in the middle of a region of max consecutive occupied orbitals. Indeed, for i ∼ max /2 we have where we have neglected the contribution of orbitals outside of the occupied region of size max . This yields for the deviation δ min = 1 − n max the following finite-size scaling (S13) In the above expression we have used the fact that at strong disorder ξ ≈ (2 ln h) −1 . Indeed, a simple perturbative expansion of any wavefunction away from its localization center yields amplitudes vanishing ∼ h −2r , where r is the distance to the localization center.
We now provide simplified expressions for various terms in Eq. (S8) in the limit h 1. In particular we justify why it is safe to ignore off-diagonal terms which vanish at large h.
s3. Diagonal terms In the large system size limit we have for the one-body term The two-body contributions, assuming constant ξ k (this is justified at strong disorder where the distribution P (ξ) is strongly peaked) reads where r = |i k 0 − i l 0 | ≥ 1 is the distance between two orbitals k and l. Therefore the two-body interaction, which reads k =l V (2) k,l n k n l ≈ ∆ k n k n k+1 + 2 h 2 n k n k+2 (S16) + 3 h 4 n k n k+3 + · · · , is clearly dominated by the nearest neighbor repulsion ∆.

s4. Off-diagonal terms
We first discuss the three-body terms of the form From Eq. (S7) we see that the amplitude can be randomly positive or negative. Nevertheless, one can easily anticipate that V k,l,m will be dominated by maximally overlaping orbitals, more precisely by three nearest-neighbor orbitals k, l, m, such that |i k where each orbital k is labelled such that φ k i ∼ exp − |i−k| ξ , meaning that its localization center i k 0 = k. In the more generic case, we have Thus according to the observables considered, there is a so-called " linear " scaling in y = y c F (L/ξ) where L is the linear size of the random graph (ie its diameter), or a "volumic" scaling in y = y c G(N /Λ) where N denotes the volume of the graph (i.e. the number of sites) and Λ is a characteristic volume diverging exponentially at the transition. Indeed, on the random graphs considered, the volume increases exponentially with the system size L and the linear or volumic scalings are not equivalent (contrary to the case of finite dimension). Different scalings have been found on both sides of the transition and two critical exponents exist in the localized phase: the average localization length diverges as ξ ∼ (W −W c ) −1 ( where W is the disorder strength) while the typical localization length follows a behavior: e. a critical exponent 1/2) identical to the critical behavior recently predicted by phenomenological RGs for the MBL transition [65][66][67].

Recently, the MBL transition has been analyzed in
the Hilbert space where critical and MBL regions have been shown to host multifractal many-body states occupying a sub-volumic Hilbert space region N D2 , with D 2 < 1 an h-dependent fractal dimension [36]. This type of behavior is associated with a linear scaling. In the ETH phase on the contrary, the scaling is found to be volumic, indicating that the many-body states are ergodic in the Hilbert space at the thermodynamic limit. Supporting the random network analogy, the linear scaling of the MBL region is associated to a diverging length scale ξ ∼ (W − W c ) −0.7 , while on the ETH side, the volumic scaling is associated to a diverging volume Λ ∼ e a|h−hc| −0.5 .
These observations invite us to consider non-standard scaling hypotheses: scaling in N /Λ, in L/ξ or even in ln L/λ, which can be different on either side of the transition. We will explain below the preferred scaling hypotheses and then describe the method we used to validate or not these hypotheses. Finally we will describe the results of these analyzes.
The phenomenological RG studies suggest very different critical and MBL behaviors for "average" observables, dominated by rare events (thermal bubbles), and for "typical" observables, controlled by the typical localization length ζ [63][64][65][66][67]. An idea to access the typical quantities is to consider observables that are the least affected by rare thermal bubbles. Maximally polarized sites, being found in the bulk of localized regions, well separated from thermal inclusions, make up for a good candidate to probe typical properties. This is why we focused on δ typ min .

s1. Scalings in the critical and MBL regions
The power-law behavior δ typ min ∼ L −γ ∼ exp(−γ ln L) observed at the largest sizes and predicted analytically at strong disorder suggests a scaling in ln L/λ, since with g(X) ≈ A0 − A 1 X for large ln L λ. Moreover, the observed power-law N D2 for the eigenstates multifractality on the Hilbert space [36] (see also [45,49]), which has been shown to be associated with the scaling ln N /ξ indicates, with L playing the role of N , that the same hypothesis is also valid here.
Note that we have observed a smooth dependence of the prefactor A(h) on the algebraic law for δ typ min ≈ A(h)L −γ(h) . Such dependence could be described by an irrelevant correction to the logarithmic scaling in the form g(ln L/λ) + h(ln L/λ)/ ln L with h(0) = 0 and h(X) ≈ CX for large ln L λ, C a constant. Due to the limited range of variation of L, we did not take into account that correction.
We can therefore expect that this observable obeys a volumic scaling with N /Λ, similarly to the participation entropies on the Hilbert space [36].

s2. Controlled finite-size scaling approach
To quantitatively test the compatibility of the numerical data with these different scaling hypotheses, we have adapted the controlled finite-size scaling approach [90,91] as detailed below: • We first assume a value of h c and consider the data for the observable ln[δ typ min /δ typ min (h c )] along with their uncertainties given by σ 2 + σ 2 c , where σ (σ c , respectively) is the standard deviation of ln[δ typ min ] (ln[δ typ min (h c )]). We aim at asserting whether these data can be fitted by a function of the form: where S can be either ln L, L or N , ρ ∼ (h − h c ) close to h c , and ν is a critical exponent to be determined by the fitting procedure. Importantly, we will fit the data for h > h c independently from the data for h < h c , as they can follow different scalings. The value χ of the chi-squared statistic for the best fit divided by the number of degrees of freedom N df is plotted as a function hc. χ/N df should be of order one for an acceptable fit. Volumic scaling N /Λ prevails in the delocalized regime, while the data in the localized regime are compatible with a scaling as ln L/λ. The best value of hc, corresponding to the minimum of χ/N df , corresponds to hc = 3.8 and hc = 4.2 for the delocalized data, while the scaling of the localized data does not allow to determine hc (see text). If we combine the localized and delocalized data, the optimal hc is hc = 4.2 (see Fig. 3(c)).
• A controlled finite-size scaling analysis [90,91] consists in a Taylor development of the scaling function around h = h c : with The orders of expansion have been set to n = 5 and m = 3. The total number of parameters to be determined in the fit is N p = n+m+1 (including ν).
We assessed the goodness of the fit by calculating the value χ of the chi-squared statistic for the best fit divided by the number of degrees of freedom N df = N D − N p where N D is the number of data, which should be of order one for an acceptable fit.
• We then tested systematically different values of h c and different scaling hypotheses, i.e. different choices S = ln L, L or N for the localized h > h c and delocalized h < h c phases. The best value of h c corresponds to the minimum of χ/N df , taking into account only the localized or delocalized data, or all the data by considering χ tot /N tot df ≡ (χ loc + χ deloc )/(N loc df + N deloc df ).
• To determine the uncertainties on h c and ν, we have used a bootstrap procedure. From the data, we generated 100 synthetic data sets by sampling randomly from normal distributions centered on the true data and with standard deviations given by the σs. We then fitted these data sets and calculated the χ/N df . The best value of h c for each synthetic data set was determined and the uncertainty corresponds to the standard deviation of these h c s. The fluctuations of ν within the data sets was also recorded.

s3. Results
We present here additional results to the ones reported in the letter. The analysis by finite-size scaling has been done without excluding any size or value of the disorder. The sizes are between L = 8 and L = 22 and the disorder between h = 0.4 and h = 90.
The figure S3 represents χ/N df as a function of h c in the localized (h > h c ) and delocalized (h < h c ) regimes. This quantitative test of the different scaling hypotheses confirm our expectation that volumic scaling N /Λ prevails in the delocalized regime, while the data in the localized regime are compatible with a scaling as ln L/λ.
It is not surprising to see this ln L/λ scaling in the localized phase favors large values of h c , since it is by construction good for algebraic data, which is the case for large values of h.
On the delocalized regime, the volumic scaling has a optimum around h c = 3.8 and h c = 4.2, with a value of χ deloc /N deloc df = 1.9 which is systematically reduced if we exclude the lowest values of h.
If we combine the localized and delocalized data, the optimal h c is h c = 4.2 with a χ tot /N tot df = 1.9 which is remarkably low given the considerable ranges of variation of h.
This scaling procedure assumes a power law divergence of the characteristic volume Λ in the delocalized phase. This is done for practical reasons. We have not been able to implement an exponential divergence in this procedure, but we have added non-linear corrections for ρ to allow for a more complex behavior than just a strict power-law divergence. The best fit gives a very large value of the critical exponent ν d ≈ 8. which indicates that the divergence of Λ close to h c is most probably an exponential divergence. Indeed, the figure 3(d) shows that the data for Λ are compatible with Λ ∼ e a(hc−h) −0.5 . Such a behavior has been observed for the scaling of the participation entropy in the Hilbert space [36], and is known to control the Anderson transition on random graphs in the delocalized phase [45,49]. Importantly, the volumic scaling implies that the whole delocalized phase obeys ETH at the thermodynamic limit.
The bootstrap procedure gives for the delocalized data (h < h c ) an average h c = 3.96 and a standard deviation of σ hc = 0.22. For the localized data, h c = 6.22, σ hc = 0.64 which confirms that the scaling as ln L/λ is not able determine the critical value of h c . However, combining the localized and delocalized data gives h c = 4.46, σ hc = 0.32. We therefore combine the delocalized and total estimations to give the estimation h c = 4.2(5).
On the localized phase, the divergence of λ ∼ (h − h c ) −ν loc corresponds to values of the critical exponent ν loc ≈ 0.5 which drifts slowly as h c increases (see figure 3(b)). The uncertainty on h c translates into a larger uncertainty on ν loc than that of the fluctuations of ν loc found from the bootstrap procedure. We find ν loc = 0.52 (3). Importantly, this behavior is compatible with the predictions from phenomenological RGs [65][66][67], with a characteristic length scale controlling finitesize effects in the MBL phase diverging exponentially as e b(h−hc) −0.5 close to the transition. Together with the algebraic critical behavior δ typ min (h c ) ∼ L −γc and the relation found between γ and the typical localization length ζ, Eq.(9) of the letter, this also confirms that the typical localization length reaches a finite value at the transition with a square-root singularity, as predicted by phe-nomenological RGs [63][64][65][66][67].

s4. Non-parametric finite-size scaling
In parallel with the controlled finite-size scaling approach whose results are presented above and in the main text, we have performed a non-parametric finite-size scaling. This approach does not assume a parametrization for the scaling functions f and g (see Eq. (10) of the main text) and scaling parameters Λ(h) and λ(h), but instead finds the best collapse of the data by optimizing an objective function which quantifies how close the data points are from being collapsed onto a unique curve. The objective function we chose to optimize is Spearman's rank correlation coefficient, R, which is a measure of how close a cloud of points are from forming the graph of a monotonic function. This approach did not allow us to quantitatively determine the best value of h c nor the uncertainties on the critical exponent ν loc . However, assuming h c = 4.2, we obtain scaling results consistent with the parametric approach, as Fig. S4 shows. In particular, we find a power-law divergence of the MBL scaling parameter with an exponent ν loc ≈ 0.44, close to the previously reported ν loc = 0.52. In the ETH region, we observe a behavior compatible with an exponentially divergent non-ergodicity volume with exponent ν d = 0.5.