Thouless energy and multifractality across the many-body localization transition

Thermal and many-body localized phases are separated by a dynamical phase transition of a new kind. We analyze the distribution of off-diagonal matrix elements of local operators across this transition in two different models of disordered spin chains. We show that the behavior of matrix elements can be used to characterize the breakdown of thermalization and to extract the many-body Thouless energy. We find that upon increasing the disorder strength the system enters a critical region around the many-body localization transition. The properties of the system in this region are: (i) the Thouless energy becomes smaller than the level spacing, (ii) the matrix elements show critical dependence on the energy difference, and (iii) the matrix elements, viewed as amplitudes of a fictitious wave function, exhibit strong multifractality. This critical region decreases with the system size, which we interpret as evidence for a diverging correlation length at the many-body localization transition. Our findings show that the correlation length becomes larger than the accessible system sizes in a broad range of disorder strength values and shed light on the critical behavior near the many-body localization transition.

Figure 1

Mean ratio of energy level spacings $\langle r\rangle$ in the Ising model (9) as a function of disorder strength $\mathrm{\Gamma }$. The ratio interpolates between its values for the Wigner-Dyson (top dashed line) and Poisson distributions (bottom dashed line). The crossing between the curves drifts towards stronger disorder (smaller values of $\mathrm{\Gamma }$) with increasing the system size, similar to the $XXZ$ spin chain [10].

Figure 2

Spectral function $f^2\left(\omega \right)$ in the $XXZ$ model. Left: At weak disorder ($W=0.5$), $f^2\left(\omega \right)$ for different system sizes collapse onto each other for both average (solid lines) and typical (dashed) curves. Right: For disorder $W=2.75$, the collapse breaks down for dashed lines, which correspond to the typical spectral function. The Thouless energy, that was visible in the left panel at small $\omega$, now cannot be resolved, and solid lines collapse onto each other in the full range.

Figure 3

Plotting $f^2\left(\omega \right)$ as a function of $\omega /\mathrm{\Delta }$ reveals its behavior at very low energies. Top row [(a)–(c)]: $XXZ$ model. (a) For weak disorder $W=1$, the plateau in $f^2\left(\omega \right)$ for $\omega <E_{\text{Th}}$ is followed by the power-law decay. (b) For intermediate disorder $W=2$, the plateau fails to fully develop even for $L=20$ spins, which signals that $E_{\text{Th}}\sim \mathrm{\Delta }$. (c) In the MBL phase $W=5, f^2\left(\omega \right)$ decays as a power law even for $\omega <\mathrm{\Delta }$. Bottom row [(d)–(f)]: Ising model. Average $f^2\left(\omega \right)$ as a function of $\omega /\mathrm{\Delta }$ shows a developed plateau only for $\mathrm{\Gamma }=0.95$ (d). For stronger disorder [$\mathrm{\Gamma }=0.7$, (e)] the concave part is visible only for $L=14$ spins. Finally, in the MBL phase at $\mathrm{\Gamma }=0.4$ (f) the spectral function decays as a power law for all system sizes.

Figure 4

Thouless energy for the $XXZ$ model at various disorders. The growth of $E_{\text{Th}}/\mathrm{\Delta }$ with $L$ slows down as the value of the disorder is increased towards the MBLT. Already for $W=2.5$ we have $E_{\text{Th}}<\mathrm{\Delta }$ even for $L=20$. For larger disorders $E_{\text{Th}}$ is too small to be reliably determined.

Figure 5

Powers ${\varphi }_{\mathrm{av}}$ and ${\varphi }_{\mathrm{typ}}$ controlling the decay of $f^2\left(\omega \right)$ and ${\left[f^2\left(\omega \right)\right]}_{\text{typ}}$. Top: $XXZ$ model where the powers saturate to one near $W_{*}\approx 2$. At larger disorder, ${\varphi }_{\mathrm{av}}$ decreases again, while ${\varphi }_{\mathrm{typ}}$ grows above one. Bottom: Similar results for the Ising model, where the powers start to systematically deviate from each other for $\mathrm{\Gamma }\le 0.7$. This coincides with the disappearance of the last traces of the Thouless plateau in $f^2\left(\omega \right)$ in the panel Fig. 3.

Figure 6

The multifractal spectrum ${\tau }_q$ in Eq. (15) as a function of the disorder strength, $W$ for the $XXZ$ model (top panel) and $\mathrm{\Gamma }$ for the Ising model (bottom panel). The spectrum evolves from being very close to that of a metal (red lines) to the “frozen” fractal spectrum deep in the MBL phase (black line). Inset shows ${\tau }_q$ near $q\le 1$.

Figure 7

Illustration of the fits quality of $f^2\left(\omega \right)$ in the $XXZ$ model for $W=1$ (top panel) and $W=2$ (bottom panel). Note that Eq. (5) fails to capture the onset of decay of $f^2\left(\omega \right)$ at weak disorder in the top panel, thus we fitted the data with the power law (straight dashed lines). Dashed lines denote fits used to extract the power $\varphi$, whereas solid lines were used to determine $E_{\text{Th}}$ and $f^2\left(0\right)$.

Figure 8

Dependence of $f^2\left(0\right)$ in the $XXZ$ model on $L$ for different values of the disorder. For weak disorder, the data is consistent with the power-law dependence of $E_{\text{Th}}$ on the system size (the bottom dashed line corresponds to diffusive behavior). Upon increasing the disorder, $f^2\left(0\right)$ gets enhanced, but still grows slower that $1/\mathrm{\Delta }$ which is illustrated by the top dot-dashed line.

Figure 9

Average and typical exponents $\varphi$ governing the decay of $f^2\left(\omega \right)$ in the $XXZ$ model for the operator ${\sigma }_1^x{\sigma }_2^x+{\sigma }_1^y{\sigma }_2^y$. The results are similar to Fig. 5 where the operator was ${\sigma }_1^z$.

Figure 10

The finite size flow of ${\tau }_q$ in the $XXZ$ model (for operator $\widehat{O}={\sigma }_1^z$). The solid lines show ${\tau }_q$ using data for $L=10...16$, while dashed lines are extracted using $L=8...14$. Note that while ${\tau }_q$ always flows to metallic asymptote with increasing system size, this flow becomes hardly noticeable for $W\ge 2$.

Figure 11

Scaling dimension for the operator $\widehat{O}={\sigma }_1^x{\sigma }_2^x+{\sigma }_1^y{\sigma }_2^y$, in the $XXZ$ model.

Figure 12

Distribution of participation ratios $P_2$ on the ergodic side of the MBLT in the $XXZ$ model. For weak disorder, $W=0.5$ (left), the distribution becomes increasingly narrower for larger systems. However, already at disorder $W=1$ (middle) the distribution looks scale invariant, with the power law decay of probability to have $P_2>P_2^{\mathrm{typ}}$. Finally, at even stronger disorder $W=2$ (right), the distribution consistently broadens with the system size.