Many-body localization and thermalization: Insights from the entanglement spectrum

We study the entanglement spectrum in the many-body localizing and thermalizing phases of one- and two-dimensional Hamiltonian systems and periodically driven “Floquet” systems. We focus on the level statistics of the entanglement spectrum as obtained through numerical diagonalization, finding structure beyond that revealed by more limited measures such as entanglement entropy. In the thermalizing phase the entanglement spectrum obeys level statistics governed by an appropriate random matrix ensemble. For Hamiltonian systems this can be viewed as evidence in favor of a strong version of the eigenstate thermalization hypothesis (ETH). Similar results are also obtained for Floquet systems, where they constitute a result “beyond ETH” and show that the corrections to ETH governing the Floquet entanglement spectrum have statistical properties governed by a random matrix ensemble. The particular random matrix ensemble governing the Floquet entanglement spectrum depends on the symmetries of the Floquet drive and therefore can depend on the choice of origin of time. In the many-body localized phase the entanglement spectrum is also found to show level repulsion, following a semi-Poisson distribution (in contrast to the energy spectrum, which follows a Poisson distribution). This semi-Poisson distribution is found to come mainly from states at high entanglement energies. The observed level repulsion occurs only for interacting localized phases. We also demonstrate that equivalent results can be obtained by calculating with a single typical eigenstate or by averaging over a microcanonical energy window, a surprising result in the localized phase. This discovery of new structure in the pattern of entanglement of localized and thermalizing phases may open up new lines of attack on many-body localization, thermalization, and the localization transition.

Figure 1

Results from numerical exact diagonalization on the one-dimensional spin chain [Eq. (1)] showing density of states (a) and level statistics (b) for a system of $N=14$ spins for $h_x=h_z=0.25$ and $h_y=0$ corresponding to GOE (blue dots) or $h_x=h_z=h_y=0.25$ corresponding to GUE (red dots). Even without any average over disorder, we clearly see the good agreement between energy level statistics in the thermal phase and the predictions of random matrix theory for the appropriate ensemble (solid lines). The level statistics clearly differ from the Poisson distribution (brown line). We also show in the left panel the region of the spectrum (shaded area) that we will consider for the entanglement spectrum analysis (we average over the states in the middle third of the spectrum).

Figure 2

Density of states (a) and level statistics (b) for a system of $N=14$ spins for $h_x=h_z=12$ and $h_y=0$ corresponding to the orthogonal class (red dots) or $h_x=h_z=h_y=12$ corresponding to the unitary class (blue dots). Even without any average over disorder, we clearly see the good agreement with the Poisson distribution (brown line). We also show in the left panel the region of the spectrum (shaded area) that we consider for the entanglement spectrum analysis (we average over the middle third of states).

Figure 3

(a) Density of states for the entanglement spectrum in the thermalizing phase of [Eq. (1)] (with $h_x=h_z=0.25,h_y=0)$ and for a variety of different system sizes. We work with the standard entanglement cut which cuts the system into two equal halves. Note that the density of states has qualitatively the same behavior whether it is extracted from a single entanglement spectrum $\left(N=16,18\right)$ or an average of many entanglement spectra $\left(N=12,14\right)$. We have obtained qualitatively similar data for the unitary Hamiltonian. (b) Level spacing distributions for the entanglement spectrum in the thermal phase of Eq. (1) $\left(h_x=h_z=0.25\right)$ for chains obeying GOE $\left(h_y=0\right)$ and GUE $\left(h_y=0.25\right)$. The solid lines show the corresponding predictions for the energy spacings distributions, and we see that there is good agreement; i.e., the level spacings of the entanglement spectrum follow the same distribution as the level spacings of the energy spectrum. This is true even for $N=18$, where we show data from the entanglement spectrum of a single typical state.

Figure 4

Entanglement spectrum properties with the standard entanglement cut which cuts the system into two equal halves. Data are for $h_x=h_z=12$ and $h_y=0$ (orthogonal class). (a) Density of states for the entanglement spectrum in the MBL phase $\left(h_x=h_z=12,h_y=0\right)$. Note that the density of states has qualitatively the same behavior whether it is extracted from a single entanglement spectrum $\left(N=16,18\right)$ or an average of many entanglement spectra $\left(N=12,14\right)$. (b) Entanglement level splitting distribution averaged over the middle third of the spectrum $\left(N=14\right)$ of states (red dots), and for a single eigenstate $(N=18$, blue dots). We note that averaging over a narrow energy window gives the same level statistics as computing the spectrum for a single eigenstate.

Figure 5

A comparison of the ES in the localized phase $(h_x=h_z=12,h_y=0$ for the orthogonal Hamiltonian; $h_x=h_y=h_z=12$ and 24 for unitary) to the semi-Poisson form $P\left(s\right)\sim sexp(-s)$ (i.e., $\alpha =\gamma =1)$. The data in the orthogonal case are the same as those in Fig. 4. The $N=18$ data are noisier because we have not averaged over eigenstates.

Figure 6

(a) Evolution of the fitted values of $\gamma$ and $\alpha$ as we tune disorder strength for an orthogonal Hamiltonian with $h_x=h_z,h_y=0$, with $N=14$. The fit is to an ES averaged over states in the middle third of the spectrum. For each exponent we show two curves, each corresponding to a different realization of disorder. There are two sources of error in such measurement: variations between states in the same realization of disorder (ROD) and variations between ROD. The error bars represent the error within an ROD, while the error between RODs can be estimated from the differences between the two curves and is discussed in more detail in the Appendix. (b) Same but in the unitary class $h_x=h_y=h_z$.

Figure 7

(a) The density-of-states data plotted in Fig. 4, but only for $N=14$. There are two peaks, once around $\xi =0$ (which contains only $2\%$ of the states), and a much larger one around $\xi =30$. The red points correspond to the “low entanglement energy” part of the spectrum, where we expect that nonuniversal behavior is encoded, while the blue “high entanglement energy” points have the universal distribution of level spacings. (b) Level spacing distributions for the low and high entanglement energy states. The high entanglement energy states follow a semi-Poisson distribution with $\alpha =\gamma =1$, while the low entanglement energy states do not.

Figure 8

(a) Entanglement spectra for a single eigenstate in the middle of the energy spectrum (the 8000th eigenstate in a system with $N=14)$, for various values of $h_x=h_z$ and with $h_y=0$. (b) Density of states for the same system, averaged over all states in the middle third of the spectrum. We see that most of the states are concentrated in a peak above the entanglement gap, and this peak moves to higher entanglement energies as disorder is increased.

Figure 9

Entanglement spacings for a Hamiltonian with $Jx=J_y=1,h_x=h_y=0,h_z=6$, and varying $J_z$. This system is always in a localized phase, but at $J_z=0$ it is a noninteracting system. We see that in this case (red symbols) the entanglement spectrum does not display the level repulsion observed elsewhere in this work; instead, it seems to follow a Poisson distribution. Data were taken for $N=18$ restricted to the $S_z=0$ sector, averaging over the middle third of the energy spectrum. When calculating the entanglement spectra we only considered entanglement states with $S_z=1$ on the left half of the system. The blue points show that even small interactions return the spectrum to semi-Poisson statistics $\alpha =\gamma =1$. One complication when taking this data is that the lack of interactions leads to a large number of very small eigenvalues of the reduced density matrix. In order to obtain good statistics, one must throw out these levels. For the data shown we threw out all eigenvalues less than ${10}^{-24}$.

Figure 10

The two-dimensional square lattice used for the transverse field Ising model. Here we depict the case that we have considered, namely $N=18$ spin-$\frac{1}{2}$ on a square torus of $N_x=6$ sites in one direction and $N_y=3$ in the perpendicular direction. When we perform the ES, the system is cut into a part $A$ (red dotted region) of $N_A=8=4\times 2$ sites and a part $B$ (green dotted region) made of the 10 remaining sites. Note that the cut breaks both translation symmetries and the inversion symmetry.

Figure 11

(a) Energy spectrum density of states and (b) energy level statistics for the system of $N=18$ spins depicted in Fig. 10. Here we have used $h_x=h_z=-J=1$.We look at a typical momentum sector not invariant under the inversion symmetry $\left[\right(\pi /3,2\pi /3)$ black dots] and a typical momentum sector invariant under the inversion symmetry $(\pi ,0)$. In the right panel, we clearly see that level statistics is given by the GOE for the $(\pi /3,2\pi /3)$ sector while it is somehow closer to the Poisson distribution for $(\pi ,0)$ (the same is true in other inversion-symmetric sectors). We also show in the left panel the region of the spectrum (shaded area) that we consider for the entanglement spectrum analysis, setting $i_{\min }=6000$ and $i_{\max }=9500$ for $(\pi /3,2\pi /3)$ $(\simeq 24\%$ of all states in this sector) and $i_{\min }=3000$ and $i_{\max }=4500$ for $(\pi ,0)$ and ${\lambda }_{\mathrm{I}}=\pm 1$ $(\simeq 20\%$ of all states in this sector).

Figure 12

Density of states (a) and level statistics (b) for the entanglement spectra of the system with $N=18$ spins. The system geometry and the bipartition are shown in Fig. 10. Here we consider to the Hamiltonians of Eq. (7) with $h_x=h_z=-J=1$. We look at the same momentum sectors as those of Fig. 12. The average over the entanglement spectra is done over the 3500 states $[1500$ for each sector of $(\pi ,0)]$ located in the bulk of the spectrum as shown in Fig. 11.

Figure 13

(a) The $r$ level spacing parameter for the energy eigenstates for the 2D model. The gray lines are a guide for the eye, showing the predicted values for the Poisson $(r\approx 0.386$ [19]) and GOE $(r\approx 0.530$ [34]) cases. We see a transition from GOE to Poisson. For the $N=15$ data, only the 200 states near the middle of the energy spectrum were used, while the $N=12$ data include all states. (b) Energy level spacing distributions for a system with $N=12$ at representative points in the ETH and MBL phases, which follow the correct distributions.

Figure 14

(a) Entanglement density of states for the 2D case. As in the 1D case, we see a small peak at $\xi =0$, as well as a larger peak at higher entanglement entropies. The data are normalized such that each curve has an equal area. This hides the behavior at large $N$ and small $\xi$, so in the inset we plot data normalized such that each curved area is proportional to the number of states, and the small $\xi$ peak is visible. Data are taken for $h_{\mathrm{random}}=8$. (b) Entanglement level spacings for the 2D case, which follow a strict semi-Poisson distribution $\alpha =\gamma =1$. Data are taken for $h_{\mathrm{random}}=8$. We see that the data match well the expected distribution shown by the solid green line.

Figure 15

(a) Ratio of adjacent gaps $r_{\mathrm{energy}}$ for the Floquet energies of the two-bang models where $H_1$ and $H_2$ are given by Eqs. (10) and (11), respectively. The calculations were done for $N=14$ spins for a single realization of disorder. The Floquet energy data shows the same behavior of that in Ref. [38]. Gray lines show the predicted values for the Poisson $(r\approx 0.386$ [19]) and COE $(r\approx 0.530$ [34]) cases. (b) Level statistics (right) for a system of $N=14$ spins for $h=6$. Even without any average over disorder, we clearly see the good agreement with either the Poisson distribution for $\tau =0.25$ or the COE for $\tau =2$ and $\tau =6$.

Figure 16

Density of states (a) and level statistics (b) for the entanglement spectra of a system with $N=14$ spins keeping $N_A=7$ spins and $h=6$. Here we consider evolution times $\tau =0.25$ (black dots), 2 (blue dots), and 6 (red dots). The average over the entanglement spectra is done over the all the 4096 states. The origin of time is chosen according to Eq. (9). Considering only a fraction of the states does not change the picture. The level statistics in the thermal phase are well described by the GUE. In the localized phase $\left(\tau =0.25\right)$, there is a good fit to the strict semi-Poisson form $\alpha =\gamma =1$ (green line). We have also fitted the level spacing distribution to the interpolating distribution considered in Ref. [31]. The fitting parameters are $\gamma \simeq 1.68\left(6\right)$ and $\alpha =1.31\left(2\right)$ (red line). Errors bars are from the fits alone and do not include averaging over disorder or finite size effects.

Figure 17

Demonstration of how the entanglement level statistics depends on the origin of time. The blue dots are the entanglement spacings for a time-reversal-symmetric origin of time; they satisfy a GOE distribution, which is the same as the Floquet energy spacings. The green dots are for the TR-breaking drive of Eq. (9); they follow a GUE distribution. All data were taken for $N=14,h=6$. (Inset) The level spacing ratio for various system sizes, which takes the GOE value only when $\varepsilon =0$ (the TR-symmetric point).

Figure 18

Level spacing ratios $r_{\mathrm{ES}}$ for the entanglement spectrum for the same model as that shown in Fig. 15, but with a time-reversal-symmetric origin of time. The gray line shows the predicted value for the GOE.

Figure 19

This plot was obtained by extracting the values of $\alpha$ and $\gamma$ for 256 energy eigenstates at a time. The horizontal axis shows the index of the first state on the set of 256. We see that the at values at small and large energy are different from that in the middle of the spectrum. The region which we average over in the main text is shown in gray; we can see the exponents are relatively constant in this region. The different colored points correspond to different realizations of disorder; we see that changing the disorder realization can affect the extracted exponents, at least for the finite-size systems available to our numerics. Data were taken with $N=14,h_x=h_z=12,h_y=0$.

Figure 20

Fitted values of (a) $\gamma$ and (b) $\alpha$ [as defined in Eq. (5)] for the ES, for multiple sizes and averaged over multiple realizations of disorder (16 RODs for $N=12$, 4 for $N=14,16)$. Data were taken for the GOE case with $h_x=h_z=h,h_y=0$. For $N=12,14$ the entanglement spectra in each disorder realization was averaged over the one-third of energy levels in the middle of the spectrum; for $N=16$ the spectra were averaged over 1000 energy levels near the middle of the energy spectrum. $\alpha$ and $\gamma$ were extracted for each realization of disorder, and error bars are the variance of the different values. (c),(d) Same as (a),(b) but in the unitary symmetry class $h_x=h_y=h_z$. In this case the entanglement eigenvalues are smaller so that at $h>6$ we run into problems with machine precision at $N=16$.