Avalanches and many-body resonances in many-body localized systems

We numerically study both the avalanche instability and many-body resonances in strongly disordered spin chains exhibiting many-body localization (MBL). Finite-size systems behave like MBL within the MBL regimes, which we divide into the asymptotic MBL phase and the finite-size MBL regime; the latter regime is, however, thermal in the limit of large systems and long times. In both Floquet and Hamiltonian models, we identify some landmarks within the MBL regimes. Our first landmark is an estimate of where the MBL phase becomes unstable to avalanches, obtained by measuring the slowest relaxation rate of a finite chain coupled to an infinite bath at one end. Our estimates indicate that the actual MBL-to-thermal phase transition occurs much deeper in the MBL regimes than has been suggested by most previous studies. Our other landmarks involve systemwide many-body resonances: We find that the effective matrix elements producing eigenstates with systemwide many-body resonances are enormously broadly distributed. This broad distribution means that the onset of such resonances in typical samples occurs quite deep in the MBL regimes, and the first such resonances typically involve rare pairs of eigenstates that are farther apart in energy than the minimum gap. Thus we find that the resonance properties define two landmarks that divide the MBL regimes of finite-size systems into three subregimes: (i) at strongest randomness, typical samples do not have any eigenstates that are involved in systemwide many-body resonances; (ii) there is a substantial intermediate subregime where typical samples do have such resonances but the pair of eigenstates with the minimum spectral gap does not, so the size of the minimum gap agrees with expectations from Poisson statistics; and (iii) in the weaker randomness subregime, the minimum gap is larger than predicted by Poisson level statistics because it is involved in a many-body resonance and thus subject to level repulsion. Nevertheless, even in this third subregime, all but a vanishing fraction of eigenstates remain nonresonant and the system thus still appears MBL in most respects. Based on our estimates of the location of the avalanche instability, it might be that the MBL phase is only part of subregime (i) and the other subregimes are entirely in the thermal phase, even though they look localized in most respects, so are in the finite-size MBL regime.

Figure 1

Sketch of the MBL phase diagram. Top: In the infinite size limit, the transition out of the MBL phase is believed to occur due to an instability toward the formation of thermalizing avalanches. The disorder strength at which this occurs marks the $L=\infty$ limit of one of our landmarks, denoted ${\mathcal{L}}_{\text{avch}}$. Bottom: At accessible finite system sizes $L$ and/or finite times $t$, we observe the significantly larger MBL regimes. The finite-$L$ crossover in the mean spectral gap ratio $\langle r\rangle$ from random matrix to Poisson statistics is one convenient landmark, denoted ${\mathcal{L}}_r$, marking the crossover from thermal behavior into the MBL regimes; $\langle r\rangle$ probes the level repulsion between all neighboring energy levels. The MBL regimes contain a few additional landmarks that probe the behavior of rare many-body resonances in eigenstates. At strong disorder exceeding a threshold marked by ${\mathcal{L}}_{\text{swr}}$, there are no systemwide resonances (swr) in any of the eigenstates of typical samples. For weaker disorder, a small number of eigenstates display systemwide resonances, but the minimum gap (mg) in the spectrum is nonresonant and shows negligible level repulsion. The level repulsion in the minimum gap appears at even weaker disorder within the MBL regime at a landmark denoted by ${\mathcal{L}}_{\text{mg}}$. Any of these finite-system landmarks that occur in the finite-size MBL regime will drift toward the true transition as $L$ and $t$ are taken to infinity, and some might drift past the transition and end up within the asymptotic MBL phase.

Figure 2

Mean-level spacing ratio $\langle r\rangle$ as a function of disorder. The finite-size crossings drift to larger disorder as $L$ is increased. For the largest $L$, the crossings are at $\alpha \cong 5.9$ for the Floquet model left and at $W\cong 3.1$ for the Hamiltonian model right. Error bars are 68% bootstrap confidence intervals here and in all other figures, but are too small to see here. The Floquet model does not have time-reversal invariance and has CUE-level statistics in its thermal phase, while the Hamiltonian model does have time-reversal invariance and hence GOE statistics, thus the difference in the thermal values of ${\langle r\rangle }_{\text{CUE}}\approx 0.5996$ versus ${\langle r\rangle }_{\text{GOE}}\approx 0.5307$ [77].

Figure 3

Example spectra of the two superoperators for $L=5$ spins. Left: Floquet superoperator with $\alpha =3.33$. The spectrum (blue) is strictly contained inside the unit circle (gray). The slowest mode (red) is the eigenvalue with the largest modulus less than 1. Right: Liouvillian of the random field Heisenberg model coupled to a bath. The spectrum lies in the complex left half plane with nonpositive real parts. The slowest mode (red) is the eigenvalue whose real part is negative and is closest to zero.

Figure 4

Typical slowest decay rate as a function of disorder strength. The curves are the (scaled) median of the distribution of $1/\tau$ over realizations for the Floquet and Lindblad superoperators. The crossings occur when the scaling of the decay rates with system size is $1/\tau \propto 4^{-L}$ and thus indicate the location of ${\mathcal{L}}_{\text{avch}}$. Smaller sizes $L<7$ are calculated using full diagonalization, bigger sizes are computed using shift-invert diagonalization with target eigenvalue $\lambda =0$ (Lindblad) and $\sigma =1$ (Floquet). For the open Floquet system, we used at least $20000,20000,10000,4000,1000,500$ disorder realizations for $L=3,4,5,6,7,8$ respectively. For the Lindblad operator, we collected $5000,10000,10000,8000,9000,8000,1000$ disorder realizations for $L=3,4,5,6,7,8,9$. Error bars are 68% bootstrap confidence intervals.

Figure 5

Typical decay rate as a function of system size computed using perturbation theory in the weak bath-coupling limit. The curves are the (scaled) 80th percentile of the distribution of $1/\tau$ over realizations for the Floquet and Lindblad superoperators. System sizes are $L=5-12$ and $L=5-14$ for Floquet and Lindblad set ups, respectively. At least 5000 disorder realizations were used for each set of parameters except for those computed using quadruple precision and $L=14$ (Lindblad) or $L=12$ (Floquet), for which 1500–2000 realizations were used. Error bars are 68% bootstrap confidence intervals.

Figure 6

Distributions over realizations of the slowest decay rate $1/\tau$ of the Floquet superoperator for disorder strengths $\alpha =3.33,6.25,12.5,25.0$. This is part of the data shown in Fig. 4.

Figure 7

The Bloch sphere representing the 2D subspace spanned by two eigenstates and the rotation to the corresponding two demixed states. A Bloch sphere can be defined using any two orthogonal states. Here we associate the poles to the two eigenstates $|\alpha \rangle$ and $|\beta \rangle$, and the points that represent the demixed states $|a\rangle$ and $|b\rangle$ are rotated away from the poles at polar and azimuthal angles $\theta$ and $\varphi$. Note that orthogonal states in the 2D subspace represented by this Bloch sphere are antipodal on the surface of the sphere.

Figure 8

Demixing of two neighboring Floquet eigenstates. $Z$ magnetizations are shown for two neighboring eigenstates (light blue and red) of one realization of the Floquet model and the corresponding two demixed states (dark blue and red). The pair of neighboring eigenstates is chosen such that it is the one with the most significant rotation angle $\theta$ on the Bloch sphere [see Eq. (16)], while also requiring that the resonance spans the system. Top: Floquet model at $\alpha =12$. In a strongly localized system we find a rare two-state resonance and are able to undo the resonance by demixing the eigenstates into highly magnetized states. Bottom: Floquet model at $\alpha =1$. In the thermal phase, resonances involve many states and so a pair of eigenstates in isolation are not able to be successfully demixed into highly magnetized states.

Figure 9

Distributions of matrix elements for end-to-end processes in the MBL regime. These distributions are collected over such pairs of adjacent eigenstates (see main text), and extracted via the method detailed in Sec. 5. Left: The Floquet circuit model. The data was collected from ${10}^4,{10}^4,1.5\times {10}^3$ realizations at $\alpha =10$ and $L=8,10,12$, respectively. Right: The Hamiltonian model. The data was collected from ${10}^4$ realizations at $W=6$ and $L=10,12,14$. See Appendix pp2 for a discussion of numerical errors in finite-precision arithmetic.

Figure 10

End-to-end resonances: Mutual information and the maximum entanglement bottleneck. The finite-size crossings are estimates of the location of ${\mathcal{L}}_{\text{swr}}$, which we estimate to be at $\alpha >13$ and $W>8.5$ for $L$ larger than we can access. The dashed lines indicate the maximal value of $ln2$. The Floquet data is on the left and the Hamiltonian on the right. The Floquet data makes use of the full spectrum. The number of realizations used to compute statistics are $32\times {10}^3$ for $L\le 10, 5\times {10}^3$ for $L=12$, and $3.5\times {10}^2$ for $L=14$. The Hamiltonian data is calculated from the center fifth of the spectrum in the ${\sum }_i{Z}_i=0$ sector. The number of realizations used to compute statistics are $64\times {10}^3$ for $L\in \{10,12\}, 12\times {10}^3$ for $L=14$, and $5\times {10}^2$ for $L=16$. Top: The median over realizations of the maximum over eigenstates of the quantum mutual information between end sites $I^{1,L}$. The median is used to target typical realizations. We have dropped data points at high $\alpha$ and $W$ that are affected by finite numerical precision. Bottom: The median over realizations of the maximum over eigenstates of $S_{\text{min}}$, the minimum entanglement entropy over cuts.

Figure 11

Typical maximum end-to-end $G$ ratio. The crossings are estimates of ${\mathcal{L}}_{\text{swr}}$. med($·$) denotes the median. The median is used to target typical realizations, and the maximum is taken over pairs of neighboring eigenstates that are systemwide (differ in the signs of their magnetizations on both end sites). Again the Floquet data is on the left and the Hamiltonian data is on the right. The system sizes $L$, disorder strengths $\alpha$ and $W$, and number of realizations are all the same as in Fig. 10.

Figure 12

Typical mean end-to-end matrix element $|U_{ab}|$, scaled by $4^L$ (left) and $2^L$ (right), for different system sizes $L$ as a function of $\alpha$ for the Floquet circuit. The crossings are estimates of ${\mathcal{L}}_{\text{swr}}$ (left) and ${\mathcal{L}}_r$ (right). The median is taken over realizations and the mean is taken over systemwide pairs of neighboring eigenstates. As explained in the main text, the reason we scale by $4^L$ on the left is to make a quantity that is proportional to the number of resonant systemwide pairs of neighboring eigenstates in a realization. On the right, we scale the same data by $2^L$ so the finite-size crossing occurs when the number of resonances is $\sim 2^L$, i.e., an $O\left(1\right)$ fraction of states.

Figure 13

Level repulsion of the smallest gaps and scaling of the typical matrix element. Floquet data is on the left and Hamiltonian data is on the right, as in the above figures. The crossings are estimates of ${\mathcal{L}}_{\text{mg}}$, which we estimate to be at $\alpha >7.9$ and $W>5.7$ as $L$ increases. Top: The average minimum gap (Floquet) and minimum gap ratio (Hamiltonian). Data is scaled and shifted so the value for a Poisson (uncorrelated) spectrum is 0 for both quantities. There is not an important difference between average and typical values for these quantities. Bottom: The log of the typical scaled matrix element (restricted to neighboring pairs of states but not to systemwide pairs). $D$ is the number of states in a realization: $2^L$ for Floquet and $\frac{1}{5}\left(\genfrac{}{}{0pt}{}{L}{L/2}\right)$ for the Hamiltonian model because we restrict to the ${\sum }_i{Z}_i=0$ sector and take the middle fifth of eigenstates. As explained in the main text, we scale the matrix elements by $D^2$ because the minimum gap is $\sim D^{-2}$. We have dropped data points at high $\alpha$ and $W$ that are affected by finite numerical precision.

Figure 14

Cumulative probability density function (CDF) of the minimal gap for different system sizes $L$ and disorder strengths $\alpha$ for the Floquet unitary circuit. The dashed lines show the expected CDF $P[{\text{min}}_n\left({\delta }_n\right)<\delta )]=1-exp(4^L\delta /\left(2\pi \right))$ for the case of completely uncorrelated eigenvalues. The difference between the dashed and solid curves signals the development of level repulsion in the small gaps.

Figure 15

Summary of estimates of landmarks in the MBL regime for the Floquet circuit. The $L$-dependent location of a landmark, ${\alpha }_c\left(L\right)$, is obtained via the condition $Q({\alpha }_c,L)=Q({\alpha }_c,L-\delta L)$, where $Q$ is a relevant quantity whose finite-size crossing indicates the landmark and $\delta L$ is the increment by which the system size is increased for a series of data. The landmarks shown are (1) ${\mathcal{L}}_{\text{mg}}$, where the minimal gap switches from Poissonian to non-Poissonian. ${\alpha }_c\left(L\right)$ is estimated from the data shown in the top left of Fig. 13 with $\delta L=2$. 2) ${\mathcal{L}}_r$, the crossover to the thermal regime, where all states exhibit significant level repulsion. ${\alpha }_c\left(L\right)$ was estimated from data shown on the left of Fig. 2 with ${\delta }_L=2$. 3) ${\mathcal{L}}_{\text{swr}}$, the landmark that marks the edge of the deep-MBL subregime where no systemwide resonances are present. ${\alpha }_c\left(L\right)$ comes from the crossings in the top left of Fig. 10 with $\delta L=2$. 4) ${\mathcal{L}}_{\text{avch}}$, an estimate for the position of the avalanche instability from the perturbative decay rate shown in Fig. 5, using $\delta L=1$.

Figure 16

Distribution, over realizations, of the slowest rate of thermalization in Lindblad dynamics for the open XXX model. This is part of the data shown in Fig. 4.

Figure 17

Left: Distribution of decay rate $\tau$ in the open Hamiltonian setup with $L=14$ and $W=18$ at weak coupling as described in Sec. 3g. Dashed (continues) lines in both panels show the decay rate computed using quadruple (double) precision. Right: 80th percentile of $-{log}_{10}\tau$ as function of disorder $W$ at fixed system sizes. In both plots, 1000 disorder realizations are used for every set of parameters. At each disorder realization, the decay rate is computed twice, once using double precision and the other using quadruple precision. Error bars are 68% bootstrap confidence interval.