Exploration of the stability of many-body localization in $d>1$

Recent work by De Roeck et al. [Phys. Rev. B 95, 155129 (2017)] has argued that many-body localization is unstable in two and higher dimensions due to a thermalization avalanche triggered by rare regions of weak disorder. To examine these arguments, we construct several models of a finite ergodic bubble coupled to an Anderson insulator of noninteracting fermions. We first describe the ergodic region using a Gaussian orthogonal ensemble random matrix and perform an exact diagonalization study of small systems. The results are in excellent agreement with a refined theory of the thermalization avalanche that includes transient finite-size effects, lending strong support to the avalanche scenario. We then explore the limit of large system sizes by modeling the ergodic region via a Hubbard model with all-to-all random hopping: the combined system, consisting of the bubble and the insulator, can be reduced to an effective Anderson impurity problem. We find that the spectral function of a local operator in the ergodic region changes dramatically when coupling to a large number of localized fermionic states; this occurs even when the localized sites are weakly coupled to the bubble. In principle, for a given size of the ergodic region, this may arrest the avalanche. However, this back-action effect is suppressed and the avalanche can be recovered if the ergodic bubble is large enough. Thus, the main effect of the back-action is to renormalize the critical bubble size.

Figure 1

Schematic illustration of an ergodic bubble inside an Anderson insulator. The red dots represent sites in the bubble, while the blue dots are the positions of the localized states of the insulator. The blue shaded region within a distance of a localization length $\xi$ from the bubble is strongly coupled to it.

Figure 2

The two geometries under consideration. In both cases, the coupling between the first LIOM and the bubble is set to be $V_1=1.0$, independent of the localization length $\xi$.

(a) One-dimensional geometry: the other LIOMs, indexed by $\alpha \ge 2$, are located at distances $r_{\alpha }=(\alpha -1)a$ in units of $a=1$ from the ergodic bubble. The coupling strengths are $V_{\alpha }=V_1{e}^{-r_{\alpha }/\xi }$ and $V_{\alpha }=V_1{e}^{-\sqrt{r_{\alpha }/\xi }}$ for the 1D insulator and 1D insulator with subexponentially decaying wave functions, respectively. (b) Two-dimensional geometry: LIOMs are arranged in concentric layers around the ergodic bubble such that the $n\mathrm{th}$ layer (for $n\ge 2$) is at a distance $r_n=(n-1)a$ from the bubble; this layer contains $n$ LIOMs indexed by $\alpha =\frac{n(n-1)}{2}+1,\cdots ,\frac{n(n+1)}{2}$ and their coupling strength is taken to be $V_{\alpha }=V_1{e}^{-r_n/\xi }$.

Figure 3

The spectral function $\rho \left(\omega \right)$ of a local operator acting on the $M\mathrm{th}$ LIOM in a fixed disorder realization $\left\{{\epsilon }_{\alpha }\right\}$ and averaged over $N_A$ eigenstates in the middle of the band. $\rho \left(\omega \right)$ is plotted as histogram with a bin width equal to twice the many-body level spacing $2{\delta }_E$. The spectral function exhibits two peaks located at $\pm 2{\epsilon }_M$. For a large localization length $\xi$, the width $\gamma$ of these peaks is much larger than the level spacing, indicating that the LIOM is hybridized with the ergodic bubble. For a small localization length, the LIOM is not hybridized and the width of the peaks is limited by the level spacing.

Figure 4

Exact diagonalization (ED) results for the spin model defined in Sec. 3 using an ergodic bubble of $N=8$ spins coupled to $M=1,\cdots ,6$ LIOMs. We plot the ratio $\gamma /{\delta }_E$ between the width $\gamma$ of the spectral function peaks of a local operator acting on the $M\mathrm{th}$ LIOM (see Fig. 3) and the many-body level spacing ${\delta }_E$ as a function of the number of LIOMs added. Different curves and colors (online) correspond to different localization lengths in the bulk: the solid curves represent numerical ED results, whereas the dashed lines represent the theoretical expectations from Eqs. (6) and (7) based on ETH assumptions. (a) One-dimensional geometry with exponentially decaying couplings between the bubble and LIOMs. (b) One dimensional geometry with a stretched exponential decay of couplings. (c) Two-dimensional geometry with an exponential decay of couplings. In (b) and (c) the failure of the avalanche is due to the subcritical bubble size. In these cases, the avalanche can be restored by increasing the bubble size $N$ or by increasing the bare bath LIOM coupling $V_1$, as verified in the inset of (b) for the one-dimensional model with stretched exponentials.

Figure 5

Increase of the effective bath size due to adding the $M\text{th}$ LIOM as measured by the change of the Thouless parameter [see Eq. (9)]. (a) One-dimensional geometry with exponentially decaying couplings. The different curves and colors (online) correspond to different localization lengths: as described in the main text, a transition occurs at ${\xi }_c=2/log2\approx 2.9$. The inset shows how the distribution $p\left(\mathcal{G}\right)$ for $\xi =0.6$ evolves with the number of added spins $M$. The distribution eventually collapses to the left and broadens, signaling localization. (b) One-dimensional geometry with stretched exponential decay of interactions. (c) Two-dimensional geometry with exponentially decaying couplings. Inset: evolution of the distribution $p\left(\mathcal{G}\right)$ for $\xi =4.0$ in the one-dimensional geometry with exponentially decaying interactions. The distribution moves to the right at a constant rate, signaling thermalization.

Figure 6

The evolution of the bubble spectral function $\rho \left(\omega \right)$ for (a), (b) $V=0.5$ and (c), (d) $V=1.5$ as a function of $M$, the number of LIOMs coupled to a Hubbard bubble of $N=30$ sites. Each curve (color) corresponds to a different $M$. The inset in (a) shows the 2D geometry in which concentric layers of LIOMs are coupled to the bubble. Here, $t=1$, $U=2$, $W=2$, $\xi =10$, $n_L=1$, and $T=2$. In (b), for the weak coupling $V=0.5$, $\rho \left(\omega \right)$ does not change much even at large $M$. For comparison, we also show $\rho \left(\omega \right)$ for the noninteracting case ($U=0$). For the stronger coupling $V=1.5$, $\rho \left(\omega \right)$ changes drastically in (c) and (d) when large numbers of LIOMs are coupled to the bubble. In (d), $\rho \left(\omega \right)$ for the interacting case becomes identical to the noninteracting one over an interval $-W<\omega <W$ for large $M$.

Figure 7

(a)–(c) The imaginary part of the three different contributions, ${\mathrm{\Sigma }}_t$ (due to the hopping), ${\mathrm{\Sigma }}_U$ (due to the interactions in the bubble), and ${\mathrm{\Sigma }}_V$ (due to the bubble-LIOM coupling), in the self-energy of the bubble for $V=1.5$ at three different values of $M$. We find that ${\mathrm{\Sigma }}_V$ dominates the other two self-energies at large $M\gtrsim 300$. (d) and (e) compare the typical values (see the main text for the definition) of the three contributions for $V=1.5$ and 0.5, respectively. For $V=1.5$ in (d), ${\mathrm{\Sigma }}_V$ crosses ${\mathrm{\Sigma }}_t$ and ${\mathrm{\Sigma }}_U$ as a function of $M$, signaling the onset of the strong back-action in the spectral function, as shown in Fig. 6.

Figure 8

(a)–(d) Correspond to the spectral functions of a bath spin, whereas (e)–(h) correspond to the spectral functions of a LIOM in a 2D geometry with exponentially decaying couplings. The curves in (a)–(f) have been averaged over $N_A$ eigenstates in the middle of the spectrum and over many disorder realizations. The curves in (g) and (h) correspond to a given disorder realization.

Figure 9

The integral over all frequencies for the typical (geometrically averaged over eigenstates and disorder realizations) spectral function of a local operator acting on a bath spin: $\mathcal{I}={\int }_{-\infty }^{\infty }{\rho }_{\mathrm{typ}}\left(\omega \right)d\omega$. We plot $\mathcal{I}$ as a function of the LIOMs added and each curve (color) corresponds to a different localization length. (a) Corresponds to a $d=1$ geometry with exponentially decaying couplings. (b) Corresponds to a $d=1$ geometry with stretched exponentials. (c) Corresponds to a $d=2$ geometry with exponentially decaying couplings.

Figure 10

The entanglement entropy $S$ of the farthest LIOM coupled to the ergodic bubble as a function of $M$. Each curve (color) corresponds to a different localization length and we have averaged over $N_A$ eigenstates in the middle of the spectrum and over many disorder realizations. (a) Corresponds to a $d=1$ geometry with exponentially decaying couplings. (b) Corresponds to a $d=1$ geometry with stretched exponentials. (c) Corresponds to a $d=2$ geometry with exponentially decaying couplings.

Figure 11

The entanglement entropy $S$ of the farthest LIOM coupled to the ergodic region as a function of the localization length $\xi$. Each curve (color) corresponds to a different total number $M$ of LIOMs (i.e., a different system size). (a) Corresponds to a $d=1$ geometry with exponentially decaying couplings. (Inset) A zoom-in around the crossing point ${\xi }_c\sim 2.9$ between the different curves and this value of ${\xi }_c$ is in very good agreement with the one obtained in Refs. [9, 21]. (b) Corresponds to a $d=1$ geometry with stretched exponentials. (Inset) A zoom-in around $\xi \sim 0.7$ between the curves, showing that there is no crossing (i.e., it is not a transition). (c) Corresponds to a $d=2$ geometry with exponentially decaying couplings.