Eigenstate thermalization and quantum chaos in the Holstein polaron model

The eigenstate thermalization hypothesis (ETH) is a successful theory that provides sufficient criteria for ergodicity in quantum many-body systems. Most studies were carried out for Hamiltonians relevant for ultracold quantum gases and single-component systems of spins, fermions, or bosons. The paradigmatic example for thermalization in solid-state physics are phonons serving as a bath for electrons. This situation is often viewed from an open-quantum-system perspective. Here, we ask whether a minimal microscopic model for electron-phonon coupling is quantum chaotic and whether it obeys ETH, if viewed as a closed quantum system. Using exact diagonalization, we address this question in the framework of the Holstein polaron model. Even though the model describes only a single itinerant electron, whose coupling to dispersionless phonons is the only integrability-breaking term, we find that the spectral statistics and the structure of Hamiltonian eigenstates exhibit essential properties of the corresponding random-matrix ensemble. Moreover, we verify the ETH ansatz both for diagonal and off-diagonal matrix elements of typical phonon and electron observables, and show that the ratio of their variances equals the value predicted from random-matrix theory.

Figure 1

Quantum chaos and thermalization (right panel) emerging from coupling of a single electron (upper left) to independent local harmonic oscillators (lower left).

Figure 2

The density of eigenstates $W\left(E\right)$ of the Holstein polaron model using a bin width $\mathrm{\Delta }E=E_{\mathrm{av}}/20$. The system parameters are ${\omega }_0/t_0=1/2, \gamma /t_0=1/\sqrt{2}, L=8$, and $M=3$, which correspond to $\mathcal{D}=L{(M+1)}^L\approx 5\times {10}^5$ eigenstates. The smooth line is a Gaussian function [see Eq. (6)], where $\mathrm{\Gamma }$ is given by Eq. (8).

Figure 3

Level-spacing statistics $\mathcal{P}\left(s\right)$ after the unfolding procedure (see the text for details) for $L=8, M=3, {\omega }_0/t_0=0.5$, and $\gamma /t_0=1/\sqrt{2}$. Numerical results are shown as a histogram for the $k=2\pi /L$ momentum sector. We consider eigenstates within an energy density window given by $\eta =2/3$ in Eq. (11), which for the particular model parameters corresponds to $90\%$ of all eigenstates. Smooth lines are the Wigner-Dyson distribution ${\mathcal{P}}_{\mathrm{WD}}\left(s\right)$, given by Eq. (12).

Figure 4

Statistics of the restricted gap ratio ${\tilde{r}}_{\alpha }$ introduced in Eq. (14). Symbols in (a), (c), and (d) show the average ${\tilde{r}}_{\mathrm{av}}$ as a function of $\gamma /t_0$. We show results at ${\omega }_0/t_0=1/2, M=3$ and three different system sizes $L$ in (a), results at $L=8, M=3$, and three different values of ${\omega }_0/t_0$ in (c), and results at ${\omega }_0/t_0=1/2, M=1$, and for three different system sizes $L$ in (d). Solid lines in (a), (c), and (d) denote ${\tilde{r}}_{\mathrm{GOE}}\approx 0.5307$, which is obtained numerically for asymptotically large systems [92], while the dashed line in (a) denotes the analytical result ${\int }_0^1\tilde{r}{\mathcal{P}}_{\mathrm{GOE}}\left(\tilde{r}\right)\mathrm{d}\tilde{r}=4-2\sqrt{3}\approx 0.5359$ obtained from the GOE distribution ${\mathcal{P}}_{\mathrm{GOE}}\left(\tilde{r}\right)$ for $3\times 3$ matrices, see Eq. (15). (b) Histogram of the distribution $\mathcal{P}\left(\tilde{r}\right)$ at ${\omega }_0/t_0=1/2, \gamma /t_0=1/\sqrt{2}, L=8$, and $M=3$. The solid line is ${\mathcal{P}}_{\mathrm{GOE}}\left(\tilde{r}\right)$ and the dashed line is the distribution $2/{(1+\tilde{r})}^2$ obtained using Poisson statistics [92].

Figure 5

Diagonal matrix elements of observables $\langle \alpha \left|\widehat{O}\right|\alpha \rangle$ for ${\omega }_0/t_0=1/2, \gamma /t_0=1/\sqrt{2}$, and $M=3$ in the quasimomentum sector $k=2\pi /L$. Results are plotted as a function of $E_{\alpha }/E_{\mathrm{av}}$, where the average energy $E_{\mathrm{av}}$ is defined in Eq. (7). Symbols from the back to the front represent results for $L=6$ (red), $L=7$ (dark blue), and $L=8$ (light blue), respectively. Arrows in (a) indicate the trend of the width of the matrix-element distribution as $L$ increases.

Figure 6

Statistics of eigenstate-to-eigenstate fluctuations $z_{\alpha }\left(O\right)$, see Eq. (20), for the same observables and model parameters as in Fig. 5. We show the mean ${\langle z\rangle }_{\eta }\left(O\right)$, see Eq. (21), and the maximum $z_{\max }\left(O\right)$, see Eq. (22). We include eigenstates from the quasimomentum sectors with $0\le k\le \pi$ that belong to the subspace ${\mathcal{Z}}_{\eta }$, defined in Eq. (19). At $\eta =1/3$, the fraction of eigenstates included in ${\mathcal{Z}}_{\eta }$ (out of the total number of eigenstates from the quasimomentum sectors $0\le k\le \pi$) are $37\%, 44\%, 50\%$, and $55\%$ for $L=5,6,7,8$, respectively, while at $\eta =2/3$ the fractions are $70\%, 79\%, 85\%$, and $90\%$, respectively. Solid lines are guides to the eyes, signaling an exponential decay with $L$.

Figure 7

Sketch of a matrix of size ${\mathcal{D}}^{\prime }$ with matrix elements $|\langle \alpha |\widehat{O}|\beta \rangle |$. The shaded region in (a) represents the matrix elements included in the sum in Eq. (23). In (b), each square contains $\mu$ eigenstates, for which we calculate the variance of diagonal matrix elements ${\sigma }_{\mathrm{diag}}^2$, Eq. (25), and the variance of off-diagonal matrix elements ${\sigma }_{\mathrm{offdiag}}^2$, Eq. (26).

Figure 8

Off-diagonal matrix elements of observables $|\langle \alpha |\widehat{O}|\beta \rangle |$ for ${\omega }_0/t_0=1/2, \gamma /t_0=1/\sqrt{2}$, and $M=3$ in the quasimomentum sector $k=2\pi /L$ with the Hilbert-space dimension ${\mathcal{D}}^{\prime }={(M+1)}^L$. The target average energy of pairs of eigenstates is $\overline{\varepsilon }=\overline{E}/E_{\mathrm{av}}=1$. (a) Averages $\overline{|O_{\alpha ,\beta }|}$, defined in Eq. (23), for system sizes $L=5,6,7,8$. Filled symbols are obtained using $\mathrm{\Delta }={10}^{-3}$ in Eq. (23) while open circles behind the filled symbols are obtained using $\mathrm{\Delta }={10}^{-1}$. Lines are fits to $a{\left({\mathcal{D}}^{\prime }\right)}^{-b}$ for $L\ge 6$ and $\mathrm{\Delta }={10}^{-3}$, where the exponent is $b=0.49$ (${\widehat{H}}_{\mathrm{kin}})$, 0.53 (${\widehat{m}}_{q=0}$), 0.58 (${\widehat{N}}_{\mathrm{ph}}$), and 0.58 (${\widehat{T}}_1$). (b) Matrix elements $|\langle \alpha |{\widehat{H}}_{\mathrm{kin}}|\beta \rangle |$ versus $\omega /E_{\mathrm{av}}=(E_{\alpha }-E_{\beta })/E_{\mathrm{av}}$. The dashed line in (b) is a moving average in a window $\delta \omega /E_{\mathrm{av}}=0.1$. Lines in (c) are functions ${\mathcal{F}}_O(\overline{E},\omega )$ defined in Eq. (24), which are moving averages of renormalized off-diagonal matrix elements, averaged in the same window as in (b). Results in (b) and (c) are for $L=8$ and $\mathrm{\Delta }={10}^{-3}$.

Figure 9

Off-diagonal matrix elements of observables $|\langle \alpha |\widehat{O}|\beta \rangle |$ for the electron kinetic energy ${\widehat{H}}_{\mathrm{kin}}$ and the phonon correlator ${\widehat{T}}_1$, using identical system parameters as in Fig. 8. Here, we vary target average energies of pairs of eigenstates $\overline{\varepsilon }=\overline{E}/E_{\mathrm{av}}=1,0.8,0.6,0.4$. (a) Averages $\overline{|O_{\alpha ,\beta }|}$, defined in Eq. (23), using $\mathrm{\Delta }={10}^{-1}$. Results are plotted versus the number of eigenstates $N_{\overline{\varepsilon }}$ within a given microcanonical window. Lines are fits to $a{\left(N_{\overline{\varepsilon }}\right)}^{-b}$ for $L\ge 6$ and $\overline{\varepsilon }=1$, where the exponent is $b=0.45$ (${\widehat{H}}_{\mathrm{kin}})$ and 0.53 (${\widehat{T}}_1$). (b) The function ${\mathcal{F}}_O(\overline{E},\omega )$, defined in Eq. (24), versus $\omega /E_{\mathrm{av}}=(E_{\alpha }-E_{\beta })/E_{\mathrm{av}}$. We use $\mathrm{\Delta }={10}^{-3}, L=8$ and calculate the moving average in a window $\delta \omega /E_{\mathrm{av}}=0.1$.

Figure 10

Variances of matrix elements of observables for ${\omega }_0/t_0=1/2, \gamma /t_0=1/\sqrt{2}$ and $M=3$ in the quasimomentum sector $k=2\pi /L$. Main panels: Ratio of variances ${\mathrm{\Sigma }}_{\alpha ,\mu }^2$, defined in Eq. (27), for four different observables at $L=8$. We fix $\mu =100$ and plot ${\mathrm{\Sigma }}_{\alpha ,\mu =100}^2$ for every $\alpha \in [1,{\mathcal{D}}^{\prime }-\mu ]$. Horizontal dashed lines are averages over those ratios of variances for which $e_{\alpha ,\mu =100}/E_{\mathrm{av}}\in [1-\eta ,1+\eta ]$. We use $\eta =2/3$ such that the set of eigenstates included in the average roughly corresponds to the set of eigenstates analyzed in Fig. 6 (i.e., $90\%$ of all eigenstates). We get $2.05,2.04,2.07,2.03$ in panels (a)–(d), respectively. Insets: Histograms of ${\mathrm{\Sigma }}_{\alpha ,\mu =100}^2$ for lattice sizes $L=7$ and 8, using the same set of $\alpha$ as the one that was used to obtain the horizontal lines in the main panel.

Figure 11

(a) Ratio of variances ${\mathrm{\Sigma }}_{\alpha ,\mu =100}^2$, defined in Eq. (27), for the phonon one-body correlation operator ${\widehat{T}}_1$ at $L=8$ and $M=3$. The data are identical to the ones in Fig. 10, but plotted versus the mean energy density $e_{\alpha ,\mu =100}/E_{\mathrm{av}}$. (b) Average ratio of variances $\overline{{\mathrm{\Sigma }}_{\mu }^2}$, defined in Eq. (28), as a function of $\mu$ for the same model parameters as in Fig. 10. The average is taken over the shaded region in panel (a). The numbers of ratios of variances included in the average are roughly $7.4\times {10}^2, 3.0\times {10}^3$, and $1.2\times {10}^4$ for $L=6,7,8$, respectively.