Characterizing the many-body localization transition through the entanglement spectrum

We numerically explore the many body localization (MBL) transition through the lens of the {\it entanglement spectrum}. While a direct transition from localization to thermalization is believed to obtain in the thermodynamic limit (the exact details of which remain an open problem), in finite system sizes there exists an intermediate `quantum critical' regime. Previous numerical investigations have explored the crossover from thermalization to criticality, and have used this to place a numerical {\it lower} bound on the critical disorder strength for MBL. A careful analysis of the {\it high energy} part of the entanglement spectrum (which contains universal information about the critical point) allows us to make the first ever observation in exact numerics of the crossover from criticality to MBL and hence to place a numerical {\it upper bound} on the critical disorder strength for MBL.

Many investigations have focused on the (Von Neuman) eigenstate entanglement entropy (EEE) as an order parameter for the phase transition. An early analytic paper argued [31] that if the EEE density evolved in a smooth fashion across the MBL transition, then the critical point must exhibit thermal entanglement. However, a careful parsing of results from numerical exact diagonalization [37] suggests a different picture, whereby the EEE density is discontinuous at the transition, in the thermodynamic limit. This work [37] put forth an appealing picture wherein the MBL and thermal phases are separated in finite size numerics by a 'quantum critical fan' (see figure 1(a)), just like thermodynamic quantum critical points [39]. However, exact numerical investigations to date have only seen the left side of the fan i.e. the crossover from the thermal phase to the quantum critical regime. As a result, not only is the global structure of the phase diagram still unresolved, but since the disorder strength h t for the thermal to critical crossover increases with system size, existing exact numerics are only able to Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. place a lower bound on the critical disorder strength for the MBL transition 4 . This lower bound has drifted gradually upwards as numerics on larger system sizes have become available, and while finite size scaling collapses (e.g. [36]) suggest a finite critical disorder strength, it is possible that the sizes accessed are not large enough to be in the scaling regime, a perspective reinforced by the observation that the critical exponents extracted from such analysis violate analytic bounds [33] and are therefore likely incorrect. Existing exact numerics thus cannot exclude the possibility that h t  ¥ as L  ¥, i.e. that the MBL phase observed in numerics is nothing more than a finite-size effect, and the correct phase diagram resembles figure 1(b). A direct observation of the MBL to critical crossover using exact numerical methods is thus highly desirable, both to confirm the global structure of the phase diagram, and to enable the extraction of a numerical upper bound on the critical disorder strength for the MBL transition. (Of course, it remains possible that rare region effects could turn the critical point obtained from such an extrapolation into an avoided critical point, as happens in e.g. [40].) In this work we analyze the thermal-MBL transition through a new lens, that of the eigenstate entanglement spectrum (EES) [41][42][43][44][45][46]. The EES contains far more information about the pattern of quantum entanglement than the EEE (and it may be experimentally observable [47]). In earlier work [43], we revealed a rich and universal structure in the EES of the MBL phase, and pointed out that the high energy part of the EES (the part typically discarded in numerical analyses) contains robust information about the critical regime. We now parse the system size dependence of the high energy EES-and find this provides information about both the thermalcritical and MBL-critical crossovers. Our work thus not only confirms the scenario of a 'quantum critical fan' in finite size numerics, but also suggests an upper bound on the critical disorder strength for the MBL transition.

Model
In this work we perform exact diagonalization on a spin one half Heisenberg model of length L with random fields in the x and z directions and periodic boundary conditions: ]. This model is known to have a thermal-critical transition around h 2.5 = [43]. Equation (1) is similar to the 'standard model' used in studies of many-body localization which has a random field in only one direction. We use two random fields to break the total S z conservation. Without breaking this symmetry, we have to choose a total S z sector for every entanglement  [37] is given in (a), in the thermodynamic limit it has a direct transition from the thermal to MBL phases, but at the finite sizes accessible numerically there is an intermediate 'critical' regime. Figure (b) shows an alternative phase diagram in which there is no MBL phase in the thermodynamic limit, though finite sized systems may be glassy. Previous numerical investigations have only observed transition between the thermal and critical phases, and therefore cannot determine which of these phase diagrams are correct. Our examination of the eigenstate entanglement spectrum tentatively allows us to see both sides of the quantum critical fan, thus ruling out scenario (b). The symbols correspond to where we have detected a transition out of the critical fan using the entanglement spectrum. spectra. Different sectors are available depending on whether L 5 is even or odd, and this leads to significant even/ odd effects in our data. Breaking the symmetry by introducing a field along the x axis eliminates these effects, allowing us to compare even and odd L and therefore doubling our resolution in system size.
We use exact diagonalization to obtain the eigenstates of equation (1) for system sizes in the range L = 11-17, for 200−500 realizations of disorder. Extensive disorder averaging is necessary, since sample to sample fluctuations are known to be large close to the transition [37]. For L 15  we compute all the eigenstates of the system. However since the entanglement properties depend on the energy density we use only the middle 1/3 of these eigenstates in our calculations, similar to [5,43]. For L 15 > we use the shift-invert method (see e.g. [48] to obtain 1000 eigenstates at energy densities approximately halfway between the top and bottom of the spectrum).

Entanglement density of states
Having obtained the eigenstates the next task is to extract from them some useful information. Entanglement has long been understood to be a useful quantity to probe MBL. Historically work focuses on the entanglement entropy, defined as: where YñáY | |is the density matrix of an eigenstate, and A r is the 'reduced density matrix' for subregion A obtained by tracing out all the degrees of freedom in some region B which is the complement to A 6 . The subregion A is composed of L A consecutive sites. Entanglement entropy exhibits area law behavior in the MBL phase and volume law behavior in the thermal phase, making it a useful probe of the physics of the transition. A potential drawback of it is that it distills a lot of information (e.g. the entire matrix A r ) down to a single number. More information can be potentially found in the entanglement spectrum, i l { }, which is the set of all eigenvalues of log A r -. The entanglement entropy is dominated by 'entanglement states' with low 'entanglement eigenvalue' (i.e.small i l ). However, in a previous work [43] we showed that entanglement states at high entanglement eigenvalue have universal structure in the MBL phase, and appear to carry a memory of the critical point. Motivated by this observation, we focus in the present work on the high energy part of the entanglement spectrum to tease out insights into the MBL transition. Such information can only be obtained by studying the spectrum directly since these states' contribution to the entanglement entropy is small. Note that this procedure produces a lot of data: for each realization of disorder we get 1000  eigenstates, and for each eigenstate we get an entanglement spectra with a number of eigenvalues equal to the square root of the size of the Hilbert space.
We begin by reviewing the relevant results of [43] on the properties of entanglement spectra in the ETH and MBL phases. The quantity we will be most interested in is the density of states (EDOS) of the entanglement spectrum, which was found to be qualitatively different in the two phases. Figure 2 shows examples of the EDOS in the ETH (a) and MBL (c) phases. To understand the ETH density of states, the prototypical example is a system at infinite temperature. In this case every eigenvalue of the reduced density matrix is equal. Let us define N A to be the size of the Hilbert space of the region A (the complement to which is traced out), and L A to be the length of this region in the one-dimensional systems we will study. Then Since the systems we study are not at infinite temperature we instead get a density of states narrowly distributed around L A . In the MBL phase a model wavefunction could be a pure product state which has 0 0 l = and all other i l = ¥. Away from this ideal limit we believe that the eigenstates should be products of 'l-bits', which decay exponentially in space. Therefore there will be some non-zero λ values which result from cutting the tails of these l-bits. From this intuition we expect an EDOS which decays exponentially from 0 l = . However in [43] we showed that the majority of entanglement eigenvalues are not in this exponentially decaying part but are instead in another peak which is at much higher entanglement eigenvalue. The level spacings in this higher peak obey semi-Poisson statistics, which are indicative of criticality.
We study the location of this 'high-energy peak', which turns out to be a very sensitive probe of the critical properties of the system. We have studied the scaling of the 'high energy peak' with system size L, subsystem size L A , and disorder strength h. We discuss these in turn. Note that we found that we need L 5 A  in order to have good statistics. Since for any location of entanglement cut, L A is the smaller of the two subregions, it then follows that we need L 11  . Meanwhile, the largest system sizes we can access given available numerical resources are L=17, and thus our results are limited to the range L 11 1 7   . A larger range may be accessible with specifically tailored models (e.g. [49]) or with greater computational resources, and would be an interesting problem for future work.

System size dependence
Deep in either the thermal or MBL phases, we expect the EDOS to be independent of L. In the thermal phase, our intuition from the infinite-temperature case tells us that the EDOS should depend only on the number of degrees of freedom in the traced-out region A, i.e. it depends on L A but not on L. In the MBL phase the correlation length is very small. Increasing L for fixed L A means adding degrees of freedom far from the region A, and these should not affect the entanglement between A and the rest of the system. In figure 2 we show how the EDOS changes with system sizes for a few different disorder strengths. Deep in the thermal (a) and MBL phases (c) the EDOS peak position changes little with L, as expected. But at intermediate values of h (b), which we might expect to lie in the critical region, we see significant changes in the EDOS peak position as a function of L.
The dependence of the EDOS in the critical region has a natural explanation in terms of the scenario outlined in [37] (and anticipated also in [27]). In that work it was argued that entanglement in the critical regime is dominated by sparse 'fractal' networks of long range resonances. Adding degrees of freedom far away from the entanglement cut alters the structures of these resonant networks, and hence alters the EDOS. As we increase system size at subcritical disorder, tuning through the critical to thermal crossover, these 'resonant networks' densely fill in, so that EDOS becomes insensitive to system size as discussed above. Meanwhile, as we increase system size at supercritical disorder, tuning through the critical to MBL crossover, then eventually system size exceeds the size of the resonances, such that adding degrees of freedom far from the cut no longer affects the resonant networks across the entanglement cut, or the EDOS. The key diagnostic of criticality is thus a dependence of the EDOS on total system size L, at constant h and L A .
We can quantify the motion of the density of states with system size by tracking the position of the top of the high entanglement eigenvalue peak. In figure 3 we plot this peak location as a function of L for a variety of disorder strengths, at L 5 A = . Note that the curves in this plot correspond to moving down vertical lines in figure 1. At small or large h the curves are flat, since at the system sizes we can access these disorder strengths are entirely in either the thermal or MBL phases. At intermediate disorder, e.g. h=3, we see curves which decrease steadily, implying that at the system sizes we can access these curves are all in the critical regime. In figure 4 we show the results of fitting the data in figure 3 to straight lines, which shows a non-zero slope in the critical region which seems to decay to zero in both the small-and large-disorder limits.
The data in figure 4 was the result of fitting all the data to a straight line, but if the picture in figure 1(a) is correct then the slope of such lines should decrease as L is increased (though it is possible that the range of L that we study is too narrow to observe this). If such a change occured we could imagine defining a 'crossover lengthscale' L h c ( ) for the quantum critical fan as the lengthscale where the dependence of the EDOS peak on L disappears. The observation of such a crossover at weak disorder sets a lower bound on the location of the , and in the MBL phase with h=8 (c), for L 5 A = . In the thermal phase the EDOS is peaked around a value proportional to L A and independent of L. In the MBL phase the EDOS has two peaks, one which exponentially decays from zero and the other at higher entanglement energies. In the critical phase the system moves between the two cases as a function of L. The EDOs shown we obtained by averaging over the middle third of the eigenstates in 128−512 realizations of disorder.  We can clearly detect L h c ( ) at weak disorder, and we plot its observed values as green circles in figure 1, finding as expected that it moves to large L as h is increased. Interestingly, we also observe a flattening of the curves at h 7 » , corresponding to L 15 c » , and we plot these results as green squares in figure 1. ) . This second crossover lengthscale, which is a decreasing function of disorder strength, corresponds to the right side of the quantum critical fan in figure 1, constituting to our knowledge the first observation of the critical to MBL crossover. We believe that the apparent 'uptick' in the peak position at large h and largest L is a numerical artifact (in this regime the curves are 'flat within error bars'), although it is interesting to speculate about a possible second crossover, which could modify the picture of the phase diagram in figure 1.

Dependence of entanglement spectra on subsystem size
We now explore how the peak position varies as a function of subsystem size L A , at fixed system size L. In the thermal phase the location of the peak in the EDOS should be L A µ . In the MBL phase, we also expect the peakʼs location to scale as L A , by the following argument: most degrees of freedom live a distance L Ã away from the entanglement cut. (This is just dimensional analysis-L A is the only lengthscale characterizing the cut.) However, in the eigenstates degrees of freedom are typically only entangled within a lengthscale ξ, where ξ is the localization length. Given exponential localization, a degree of freedom a distance L Ã away from the cut will have entanglement ) across the cut, and thus will have entanglement eigenvalue L A x . In figure 5 we show the peak location, as a function of L A for fixed L=17 and several disorder strengths. We see that in the thermal phase the behavior is L A µ as expected ( figure 5(a)). In the MBL phase the peak location scales approximately as L A , but there are small deviations from perfect linear behavior in the deep MBL regime ( figure 5(b)). Understanding whether these deviations are physical, and to what physics they correspond if so, is an interesting challenge for future work.

Dependence on disorder strength
We now consider the dependence on the disorder strength h. Note that our toy model above suggests that in the thermal phase, the peak location should be independent of h (depending only on L A ), whereas in the MBL phase, peak location should be L h A x ( ). Now close to the transition, h h c x~-n -( ) , whereas far from the transition (deep in the locator limit), h 1 log x~( ). Given that exact diagonalization studies typically see 0.9 n » [36], we should expect to see the peak location be (i) independent of h in the thermal regime (ii) increasing roughly linearly with h close to the MBL side of the transition, and (iii) increasing much more slowly with h in the deep MBL regime. In figure 6 we show the numerically obtained dependence of peak position on h for several system sizes. Note that the behavior is broadly consistent with the predictions of the toy model detailed above. Figure 5. The location of the maximum of the EDOS (excluding the peak around 0 l = which occurs in the MBL phase), as a function of L A for several h values and L=16. As discussed in the text, the peak location grows L A µ in all phases. The data has been broken into multiple panels to make the small h behavior easy to see.

Statistical fluctuations: variations between samples versus variations between eigenstates
In this subsection, we investigate the dominant source of the statistical fluctuations. In each sample we study, every eigenstate produces its own entanglement spectrum, and therefore its own entanglement density of states. Here we compare the density of states of eigenstates coming from the same sample to those from different samples. In addition to fundamental interest, it is important to understand these difference because they inform how we perform our data analysis. Figure 7 shows such a comparison. To obtain this figure, we created a microcanonical average EDOS out of 1000 eigenstates of a single sample. At L=14, where this data was taken, we focus 5000 » eigenstates (note that we restrict ourselves to the middle of the spectrum so that we do not have to worry about mobility edges), and therefore we can make five such different EDOS curves. We then repeat this procedure for multiple disorder realizations, with data coming from the same realization of disorder plotted in the same color. The plot shows that fluctuations between the EDOS obtained by microcanonical averaging over different energy windows within the same sample (different groups of 1000 eigenstates) are small compared to fluctuations between microcanonically averaged EDOS from different samples.
These results inform our numerical procedures in a number of ways. For large system sizes we use the shiftinvert method to obtain 1000 eigenstates in the middle of our spectrum, while for smaller system sizes we use full diagonalization. Figure 7 suggests that these methods should give very similar results, since the difference between EDOS obtianed from microcanonical averaging over different windows is much smaller than the  We can see that distinct microcanonical averages of the EDOS from the same sample are much more similar than microcanonical averages from different samples, implying that sample-to-sample variation is more important than variation between microcanonical energy windows (of 1000 states) within the same sample.
fluctuations between samples. To confirm this expectation in figure 8 we show data for 256 samples, where each point is the result of averaging over 16 samples. The first four points (64 samples) were obtained using full diagonalization while the remainder we taken using shift-invert. If shift-invert gave different results we would expect larger error bars on the points taken with shift-invert, but we do not observe this.
Previously we computed the location of the EDOS peak by first computing an EDOS which is averaged over all eigenstates in a given sample, and then averaging those peak locations. We could imagine other ways of performing this analysis. We could compute the peak location for each eigenstate, and average those peak locations over eigenstates and samples. Figure 7 tells us that both of these averaging methods should give the same results, since the EDOS for all eigenstates in a given sample and the final error bar will be dominated by sample-to-sample fluctuations in any case. We could also imagine finding the peak of an EDOS which is averaged over all eigenstates and samples. The values of the peaks extracted this way will be different: in figure 7 for example one can see that the results of averaging the peaks of the different samples will not necessarily be the same as the peak of the average distribution. In figure 9 we show the results of extracting peak values from an EDOS averaged over all samples. To determine the values in figure 9 we fit the data near the top of the curve to a Gaussian distribution, and the peak location is the maximum of this fit. Using such a method it is less clear how to assign an error bar to our data, and furthermore the data seems more noisy than that in figure 3. However we can clearly see that this method gives similar results to those of figure 3.

Discussion
We have examined the MBL transition through the lens of the entanglement spectrum. In particular, we have shown that system size dependence of the high energy peak in the entanglement spectrum is a sensitive measure of criticality, which allows us to study both the critical to thermal crossover (which many studies have seen before), and the critical to MBL crossover (which has never before been observed to our knowledge). As such, not only does our work justify a picture in which at finite size the MBL and thermal regimes are separated by a 'critical' fan, but it also suggests both a lower and an upper bound on the critical disorder strength for the MBL transition. (Of course, while observation of a critical fan suggests a critical point between the bounds we place herein, it remains possible that e.g. rare region effects could render the critical point avoided in thermodynamically large samples.) We provided two pieces of evidence suggesting that the entanglement spectrum can detect the critical-MBL crossover. The first was that at large disorder the peak location is independent of L, which is a qualitatively different behavior than the decrease of the peak location with L observed in the critical region. The second evidence was that for a given h we can observe a crossover lengthscale L h c ( ) where the peak location versus L curves become flat. If these results are correct this implies a critical-MBL crossover. However there is still work to be done to improve our results. Though figure 4 seems to show a curve which goes to zero at large h, we cannot rule out that it only asymptotically approaches zero. It may be desirable to apply greater computation resources (or new techniques) to this problem in an effort to more definitively establish whether the slope indeed vanishes. It is also possible that other systems, while exhibiting the same universal properties, might yield cleaner results at the system sizes we are able to study. We note that the data we present required thousands of hours of CPU time, so improving our data by obtaining more samples, though possible, would require significant computational resources. A similar problem exists with our putative observation of L h c ( ). It is possible that the flattening we observe is statistical noise, and higher-quality data could help us to establish this 7 . Higher quality data might also provide a tighter bound on the location of the thermodynamic transition point (recall that our current data suggests h 2 7 c < < ). Also interesting would be to apply similar analyses for models with quasi-periodic disorder, which may also be relevant for understanding the MBL transition [50], or to examine the dynamical evolution of the entanglement spectrum starting from low entanglement initial conditions, as pioneered by [51]. Figure 9. Similar to figure 3, but using a different method of determining the average peak location (error bars not shown). Instead of finding a peak location for each sample and averaging those, we determine one EDOS by averaging over all samples and plot its peak value. Though the results of the two methods are qualitatively the same we find that the data using the full EDOS is noisier, and moreover it is harder to estimate error bars.