Characterizing the many-body localization transition through the entanglement spectrum

TL;DR

This paper investigates the many-body localization transition by analyzing the entanglement spectrum in finite systems, identifying the crossover from criticality to MBL, and establishing bounds on the critical disorder strength.

Contribution

It provides the first numerical observation of the crossover from criticality to MBL using the entanglement spectrum, offering bounds on the critical disorder strength.

Findings

Identified the crossover from criticality to MBL in finite systems.

Placed an upper bound on the critical disorder strength for MBL.

Analyzed the high energy part of the entanglement spectrum to extract universal information.

Abstract

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.

Figures

Equations

$$H=J\sum _{i=1}^L{\vec{S}}_i\cdot {\vec{S}}_{i+1}+h\sum _i({\alpha }_{iz}{S}_i^z+{\alpha }_{ix}{S}_i^x),$$

$$S=-\mathrm{Tr}({\rho }_A\log {\rho }_A);{\rho }_A={\mathrm{Tr}}_B|\Psi ⟩⟨\Psi |,$$

Full text

Introduction

Many body localization (MBL) and the resulting breakdown of statistical mechanics in disordered interacting systems has been the subject of much recent research Gornyi et al. (2005); Basko et al. (2006); Žnidarič et al. (2008); Oganesyan and Huse (2007); Pal and Huse (2010); Nandkishore and Huse (2015). The intensive research has yielded a plethora of insights into the properties of this non-ergodic regime, including its connections of integrability Bardarson et al. (2012); Serbyn et al. (2013); Huse et al. (2014); Ros et al. (2015); Chandran et al. (2016); Geraedts et al. (2017), its response properties Khemani et al. (2015); Gopalakrishnan et al. (2015), and the circumstances under which the phenomenon may arise Nandkishore and Potter (2014); Nandkishore (2014); Gornyi et al. (2016); Nandkishore (2015); Potter and Vasseur (2016); De Roeck et al. (2016); Nandkishore and Gopalakrishnan (2016); De Roeck and Huveneers (2016). However, the quantum phase transition between a ‘thermal’ phase where statistical mechanics is obeyed and a ‘many body localized’ (MBL) phase where it is not continues to be an open problem. This is a dynamical transition which lies outside the usual thermodynamic frameworks, and while it has been attacked with a variety of techniques, from mean field theory Laumann et al. (2014); Gopalakrishnan and Nandkishore (2014) to the strong disorder renormalization group Vosk et al. (2015); Potter et al. (2015); Zhang et al. (2016a); Dumitrescu et al. (2017), and using both analytic arguments Grover (2014); Serbyn and Moore (2016); Chandran et al. (2015) and numerics Serbyn et al. (2015); Luitz et al. (2015); Khemani et al. (2016), (see Parameswaran et al. (2016) for a review), a complete understanding of the transition remains elusive.

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 Grover (2014) 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 Khemani et al. (2016) suggests a different picture, whereby the EEE density is discontinuous at the transition, in the thermodynamic limit. This work Khemani et al. (2016) put forth an appealing picture wherein the MBL and thermal phases are separated in finite size numerics by a ‘quantum critical fan’ [See Fig. 1(a)], just like thermodynamic quantum critical points Sachdev (2011). 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 place a lower bound on the critical disorder strength for the MBL transition 111The transition also occurs as a function of energy density, in this work we will fix the energy density to be in the middle of the density of states.. This lower bound has drifted gradually upwards as numerics on larger system sizes have become available, and while finite size scaling collapses (e.g. Luitz et al. (2015)) 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 Chandran et al. (2015) and are therefore likely incorrect. Existing exact numerics thus cannot exclude the possibility that $h_{t}\rightarrow\infty$ as $L\rightarrow\infty$, i.e. that the MBL phase observed in numerics is nothing more than a finite-size effect, and the correct phase diagram resembles Fig. 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.

In this work we analyze the thermal-MBL transition through a new lens, that of the eigenstate entanglement spectrum (EES) Li and Haldane (2008); Yang et al. (2015); Geraedts et al. (2016); Serbyn et al. (2016); Yang et al. (2017a); Pietracaprina et al. (2016). The EES contains far more information about the pattern of quantum entanglement than the EEE. In earlier work Geraedts et al. (2016), 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 allows us to detect both the thermal-critical and MBL-critical crossovers. Our work thus not only confirms the scenario of a ‘quantum critical fan’ in finite size numerics, but also allows us to upper bound 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:

[TABLE]

where $2S_{\mu}=\sigma^{\mu}$ are the Pauli matrices. Here we set $J=1$ while $h$ controls the disorder strength. The $\alpha_{i\mu}$ are random numbers drawn from a uniform distribution in $[-1,1]$. This model is known to have a thermal-critical transition around $h=2.5$ Geraedts et al. (2016). Eq. (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 spectra. Different sectors are available depending on whether $L$ 222and also $L_{A}$ 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 Eq. (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 Khemani et al. (2016). For $L\leq 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 Pal and Huse (2010); Geraedts et al. (2016). For $L>15$ we use the shift-invert method 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:

[TABLE]

where $|\Psi\rangle\langle\Psi|$ is the density matrix of an eigenstate, and $\rho_{A}$ 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$ 333What we have defined is the Von Neumann entanglement entropy, other entanglement entropies (e.g. Renyi) are defined similarly.. 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 $\rho_{A}$) down to a single number. More information can be potentially found in the entanglement spectrum, $\{\lambda_{i}\}$, which is the set of all eigenvalues of $-\log\rho_{A}$. The entanglement entropy is dominated by ‘entanglement states’ with low ‘entanglement energy’ (i.e. small $\lambda_{i}$). However, in a previous work Geraedts et al. (2016) we showed that entanglement states at high entanglement energy 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.

The quantity we will be interested in is the density of states (EDOS) of the entanglement spectrum. Fig. 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 traced-out region $A$, and $L_{A}$ to be the length of this region in the one-dimentional systems we will study. Then $\lambda_{i}=\log N_{A}\propto L_{A}$ for all $i$. 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 $\lambda_{0}=0$ and all other $\lambda_{i}=\infty$. 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 $\lambda$ values which result from cutting the tails of these l-bits. From this intuition we expect an EDOS which decays exponentially from $\lambda=0$. However in Geraedts et al. (2016) 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 energy. The level spacings in this higher peak obey semi-Poisson statistics, which are indicative of criticality. The location of this ‘high-energy peak’ 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 disorder strength $h$, subsystem size $L_{A}$, and system size $L$. In the main text we study the dependence on system size $L$, which turns out to be the most informative, and allows us to detect the both the thermal-critical and critical-MBL transitions. The scaling with $h$ and $L_{A}$ is discussed in the Appendix. Note that we found that we need $L_{A}\geq 5$ 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\geq 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 $11\leq L\leq 17$. A larger range may be accessible with specifically tailored models (e.g. Zhang et al. (2016b)) or with greater computational resources, and would be an interesting problem for future work.

Results

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 Fig. 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 Khemani et al. (2016) (and anticipated also in Potter et al. (2015)). 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 effects 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 energy peak. In Fig. 3 we plot this peak location as a function of $L$ for a variety of disorder strengths, at $L_{A}=5$. Note that the curves in this plot correspond to moving down vertical lines in Fig. 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. We can also see curves, e.g. at $h=2$ or $h=7$, which have a system size dependence only at small sizes. These curves indicate transitions between the critical regime and either the thermal or MBL phases.

The crossover lengthscale $L_{c}(h)$ for the quantum critical fan is the lengthscale where the dependence of the EDOS peak on $L$ disappears. Note that we observe a crossover lengthscale at weak disorder, which is an increasing function of $h$ (corresponding to the left side of the quantum critical fan), and also a crossover lengthscale at strong disorder. While our resolution is limited by the small range of system sizes we can access numerically, it is clear that $L_{c}(h=6)>17$, $L_{c}(h=7,8,9)=15\pm 1$, and $L_{c}(h=12)\lesssim 12$. The expert reader may wonder if a change in behavior around $L=15$ may be due to a change from full diagonalization to a shift-invert procedure, but we have verified that it is in fact physical, and not a numerical artifact (see Appendix B). This second crossover lengthscale, which is a decreasing function of disorder strength, corresponds to the right side of the quantum critical fan in Fig.1, constituting to our knowledge the first observation of the critical to MBL crossover. Note also that insofar as $h=7$ lies clearly on the right side of the quantum critical fan, we can place a numerical upper bound on the critical disorder strength in this model, $h_{c}<7$, to complement the lower bound $h_{c}\geq 2.5$ from the thermal to critical crossover. While the bounds $2.5<h_{c}<7$ are not particularly tight, the fact that we can place an upper bound at all is, in our view, a major advance - a further tightening of the bounds now requires ‘merely’ the application of greater (or more cleverly applied) computing power. We also 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 Fig.1.

Conclusions: 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 identify 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 allows us to place both a lower and an upper bound on the critical disorder strength for the MBL transition. A further tightening of these bounds, either by the application of greater computational resources or by the refinement of techniques introduced in this paper, or by some other method entirely, would be an interesting challenge for future work. Also interesting would be to apply similar analyses for models with quasi-periodic disorder, which may also be relevant for understanding the MBL transition Khemani et al. (2017), or to examine the dynamical evolution of the entanglement spectrum starting from low entanglement initial conditions, as pioneered by Yang et al. (2017b). These problems too we leave to future work.

Acknowledgements We thank S.A. Parameswaran, R. Vasseur, A.C. Potter, R. N. Bhatt and especially Vedika Khemani for feedback on the manuscript. R.N. is supported by the Foundational Questions Institute (fqxi.org; grant no. FQXi-RFP-1617) through their fund at the Silicon Valley Community Foundation. S.D.G. is supported by Department of Energy BES Grant DE-SC0002140.

Appendix A Dependence of entanglement spectra on subsystem size and disorder strength

In the main text we focussed on the dependence of the entanglement density of states on system size, since this is the quantity which reveals the systems critical properties. In this Appendix we check our understanding of the EDOS by confirming that it behaves correctly as a function of $L_{A}$, the subsystem size, and $h$, the disorder strength.

First we study the dependence on $L_{A}$. As discussed in the main text, in the thermal phase the location of the peak in the EDOS should be $\propto 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 $\sim L_{A}$ 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 $\xi$, where $\xi$ is the localization length. Given exponential localization, a degree of freedom a distance $\sim L_{A}$ away from the cut will have entanglement $\sim\exp(-L_{A}/\xi)$ across the cut, and thus will have entanglement energy $\sim L_{A}/\xi$. In Fig. 4 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 $\propto L_{A}$ as expected [Fig. 4(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[Fig. 4(b)]. Understanding whether these deviations are physical, and to what physics they correspond if so, is an interesting challenge for future work.

We can also 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 $\sim 1/\xi(h)$. Now close to the transition, $\xi\sim(h-h_{c})^{-\nu}$, whereas far from the transition (deep in the locator limit), $\xi\sim 1/\log(h)$. Given that exact diagonalization studies typically see $\nu\approx 0.9$ (Luitz et al. (2015)), 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 Fig. 5 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.

Appendix B Variations between samples vs. variations between eigenstates

In each sample we study, every eigenstate produces its own entanglement spectrum, and therefore its own entanglement density of states. In this Appendix 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 6 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 $\approx 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 shift-invert method to obtain $1000$ eigenstates in the middle of our spectrum, while for smaller system sizes we use full diagonalization. Fig. 6 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 fluctuations between samples. To confirm this expectation in Fig. 7 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.

In the main text we compute 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. Fig. 6 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 Fig. 6 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 Fig. 8 we show the results of extracting peak values from an EDOS averaged over all samples. To determine the values in Fig. 8 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 Fig. 3, and for these reasons we did not use this method in the main text. However we can clearly see that this method gives similar results to those of the main text.