Eigenstate thermalization and quantum chaos in the Holstein polaron model

TL;DR

This paper investigates whether the Holstein polaron model, a minimal electron-phonon system, exhibits quantum chaos and thermalization properties consistent with the eigenstate thermalization hypothesis (ETH) using exact diagonalization.

Contribution

It demonstrates that despite its simplicity, the Holstein polaron model shows spectral and eigenstate properties indicative of quantum chaos and obeys ETH, extending understanding to solid-state systems.

Findings

Spectral statistics match random-matrix ensemble predictions

Eigenstates exhibit properties consistent with quantum chaos

ETH holds for both diagonal and off-diagonal matrix elements

Abstract

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 offdiagonal matrix elements of typical phonon and electron observables, and show that the ratio of their variances equals the value predicted from random-matrix theory.

Figures

Equations

$$\hat{H}={\hat{H}}_{\mathrm{kin}}+{\hat{H}}_{\mathrm{ph}}+{\hat{H}}_{\mathrm{eph}}.$$

$${\hat{H}}_{\mathrm{kin}}=-t_0\sum _{j=1}^L\left({\hat{c}}_j^{†}{\hat{c}}_{j+1}^{}+{\hat{c}}_{j+1}^{†}{\hat{c}}_j^{}\right),$$

$${\hat{H}}_{\mathrm{ph}}=\hslash {\omega }_0\sum _{j=1}^L{\hat{b}}_j^{†}{\hat{b}}_j^{},$$

$${\hat{H}}_{\mathrm{eph}}=\gamma \sum _{j=1}^L\left({\hat{b}}_j^{†}+{\hat{b}}_j^{}\right){\hat{n}}_j.$$

$$W(E)={\mathcal{D}}^{-1}\sum _{\alpha =1}^{\mathcal{D}}\delta (E-E_{\alpha }).$$

$$g(E)=E_{\mathrm{av}}/(\sqrt{2\pi }\Gamma )e^{-(E-E_{\mathrm{av}}{)}^2/(2{\Gamma }^2)},$$

$$E_{\mathrm{av}}=⟨\hat{H}⟩={\mathcal{D}}^{-1}\mathrm{Tr}\{\hat{H}\}=L\frac{{\omega }_{0}M}{2}$$

$${\Gamma }^2=⟨{\hat{H}}^2⟩-⟨\hat{H}{⟩}^2=L\frac{{\omega }_0^{2}M(M+2)}{12}+2t_0^2+{\gamma }^{2}M.$$

$$\frac{E_{\mathrm{av}}-E_{\alpha }}{E_{\mathrm{av}}-E_{\min }}<\eta \text{mbox}if{E}_{\alpha }<E_{\mathrm{av}}$$

$$\frac{E_{\alpha }-E_{\mathrm{av}}}{E_{\max }-E_{\mathrm{av}}}<\eta \text{mbox}if{E}_{\alpha }>E_{\mathrm{av}}.$$

$${\mathcal{Z}}_{\eta }^{\{k\}}=\{|\alpha ⟩;E_{\alpha }/E_{\mathrm{av}}\in [1-\eta ,1+\eta ]\},$$

$${\mathcal{P}}_{\mathrm{WD}}(s)=\frac{\pi s}{2}{e}^{-\pi s^2/4},$$

$$r_{\alpha }=\frac{{\delta }_{\alpha }}{{\delta }_{\alpha -1}},$$

$${\tilde{r}}_{\alpha }=\frac{\min \{{\delta }_{\alpha },{\delta }_{\alpha -1}\}}{\max \{{\delta }_{\alpha },{\delta }_{\alpha -1}\}}=\min \{r_{\alpha },r_{\alpha }^{-1}\}.$$

$${\mathcal{P}}_{\mathrm{GOE}}(\tilde{r})=\frac{27}{4}\frac{\tilde{r}+{\tilde{r}}^2}{(1+\tilde{r}+{\tilde{r}}^2{)}^{5/2}},$$

$$⟨\alpha |\hat{O}|\beta ⟩=O(\bar{E}){\delta }_{\alpha ,\beta }+e^{-S(\bar{E})/2}{f}_O(\bar{E},\omega )R_{\alpha ,\beta },$$

$${\hat{m}}_q=\frac{1}{L}\sum _{j,l=1}^L{e}^{-i(l-j)q}{\hat{c}}_j^{†}{\hat{c}}_l$$

$${\hat{T}}_1=\frac{1}{L}\sum _{j=1}^L\left({\hat{b}}_j^{†}{\hat{b}}_{j+1}+H.c.\right).$$

$${\mathcal{Z}}_{\eta }=\bigcup _{0\le k\le \pi }{\mathcal{Z}}_{\eta }^{\{k\}}.$$

$$z_{\alpha }(O)=⟨\alpha +1|\hat{O}|\alpha +1⟩-⟨\alpha |\hat{O}|\alpha ⟩$$

$$⟨z{⟩}_{\eta }(O)=||{\mathcal{Z}}_{\eta }|{|}^{-1}\sum _{|\alpha ⟩\in {\mathcal{Z}}_{\eta }}|z_{\alpha }(O)|,$$

$$z_{\max }(O)=\underset{|\alpha ⟩\in {\mathcal{Z}}_{\eta }}{\max }|z_{\alpha }(O)|.$$

$$\overline{|O_{\alpha ,\beta }|}=\frac{1}{\mathcal{N}}\sum _{\begin{array}{c}\alpha ,\beta ;\alpha \ne \beta \\ |(E_{\alpha }+E_{\beta })/(2E_{\mathrm{av}})-\bar{\epsilon }|<\Delta /2\end{array}}|⟨\alpha |\hat{O}|\beta ⟩|,$$

$${\mathcal{F}}_O(\bar{E},\omega )=\mathrm{MA}\left(\frac{|⟨\alpha |\hat{O}|\beta ⟩|}{\overline{|O_{\alpha ,\beta }|}}\right),$$

$$[{\sigma }_{\mathrm{diag}}^{(\alpha ,\mu )}(O){]}^2=⟨(O_{\alpha ,\alpha }{)}^2{⟩}_{\mu }-⟨O_{\alpha ,\alpha }{⟩}_{\mu }^2,$$

$$[{\sigma }_{\mathrm{offdiag}}^{(\alpha ,\mu )}(O){]}^2=⟨|O_{\alpha ,\beta }{|}^2{⟩}_{\mu }-|⟨O_{\alpha ,\beta }{⟩}_{\mu }{|}^2,$$

$${\Sigma }_{\alpha ,\mu }^2(O)=\frac{[{\sigma }_{\mathrm{diag}}^{(\alpha ,\mu )}(O){]}^2}{[{\sigma }_{\mathrm{offdiag}}^{(\alpha ,\mu )}(O){]}^2}.$$

$$\overline{{\Sigma }_{\mu }^2}(O)=||{\mathcal{S}}_{\mu }|{|}^{-1}\sum _{{\Sigma }_{\alpha ,\mu }^2\in {\mathcal{S}}_{\mu }}{\Sigma }_{\alpha ,\mu }^2(O)$$

Full text

I Introduction

Understanding whether and how an isolated quantum many-body system approaches thermal equilibrium after being driven far from equilibrium has been a tremendous theoretical challenge since the birth of quantum mechanics vonneumann_29 ; Goldstein2010 . This topic became viral after quantum-gas experiments Bloch2008 ; Langen2015 ; Eisert2015 and other quantum simulators such as ion traps became available that are, to a good approximation, closed quantum systems. Using these platforms, a series of experiments on nonequilibrium dynamics was conducted, focusing either on the generic case of thermalizing systems Trotzky2012 ; Kaufman2016 ; Neill2016 ; tang_kao_18 or the exceptions such as integrable models Kinoshita2006 ; gring_kuhnert_12 ; langen15a or many-body localization Schreiber2015 ; Choi2016 ; Smith2016 .

An additional impetus to the theoretical community was given by the work of Rigol et al. rigol_dunjko_08 who demonstrated that the eigenstate thermalization hypothesis (ETH), pioneered by Deutsch deutsch_91 and Srednicki srednicki_94 , provides a relevant framework to describe statistical properties of eigenstates of lattice Hamiltonians, applicable to ongoing experiments with ultracold atoms on optical lattices Kaufman2016 ; Choi2016 ; Schreiber2015 . As a crucial consequence, if the ETH is satisfied in a given quantum many-body system, it implies thermalization of local observables after a unitary time evolution rigol_srednicki_12 ; dalessio_kafri_16 ; mori_ikeda_18 ; deutsch_18 .

So far, the ETH has been verified for a wide number of lattice models such as nonintegrable spin-1/2 chains rigol_09a ; rigol_santos_10 ; steinigeweg_herbrych_13 ; beugeling_moessner_14 ; kim_ikeda_14 ; luitz_16 ; luitz_barlev_16 ; garrison_grover_18 ; richter_gemmer_18 ; khaymovich_haque_18 ; hamazaki_ueda_19 , ladders beugeling_moessner_14 ; steinigeweg_khodja_14 ; khodja_steinigeweg_15 ; beugeling_moessner_15 and square lattices fratus_srednicki_15 ; mondaini_fratus_16 ; mondaini_rigol_17 , interacting spinless fermions rigol09 ; neuenhahn_marquardt_12 , Bose-Hubbard sorg14 ; beugeling_moessner_14 and Fermi-Hubbard chains mondaini_rigol_15 , dipolar hard-core bosons khatami_pupillo_13 , quantum dimer models lan_powell_17 and Fibonacci anyons chandaran_schulz_16 . In these examples, mostly, direct two-body interactions in systems of either spins, fermions or bosons are responsible for rendering the system ergodic. In condensed matter systems, however, the presence of phonons is ubiquitous. In fact, the textbook example of thermalization in solid-state physics is the phonons being a bath for the electrons. There, the distinction between bath and system is not made by a real-space bipartition that is often considered in ETH studies, but by tracing out one type of physical degree of freedom, namely the phonons. Moreover, in metals, direct electron-electron interactions are often irrelevant while the coupling to phonons leads to effective interactions responsible for thermalization and transport.

A second timely motivation to study nonequilibrium problems in electron-phonon coupled systems stems from the experimental advances with time-resolved spectroscopy giannetti_capone_16 . Such experiments have stimulated an increased interest in exploring nonequilibrium dynamics of isolated electron-phonon coupled systems theoretically ku2007 ; vidmar11 ; fehske11 ; vidmar11c ; golez12b ; li13 ; hohenadler13 ; kemper13 ; sentef13 ; murakami_werner_15 ; dorfner_vidmar_15 ; werner_eckstein_15 ; brockt_dorfner_15 ; kogoj_vidmar_16 ; kogoj_mierzejewski_16 ; huang_chen_17 ; hashimoto_ishihara_17 ; chen_borrelli_17 ; huang_hoshina_19 . Among them, systems with a single excited electron coupled to phonons (the so-called polaron case) provide a platform for accurate numerical simulations bonca99 ; fehske2007 ; vidmar11 ; dorfner_vidmar_15 ; brockt_dorfner_15 , allowing us to demonstrate, for certain nonequilibrium protocols, equilibration golez12b ; dorfner_vidmar_15 and indications for thermalization kogoj_vidmar_16 .

In this paper, we explore quantum ergodicity of a polaron system by studying indicators of quantum chaos and the validity of the ETH, which implies thermalization for generic far-from-equilibrium initial states. We focus on the one-dimensional Holstein polaron model, which describes properties of a single electron locally coupled to dispersionless phonons, as sketched in Fig. 1. Using an exact-diagonalization analysis, we demonstrate that despite the apparent simplicity of the model, it exhibits clear manifestations of quantum chaos and the Hamiltonian eigenstates in the bulk of the spectrum obey the ETH. Therefore, the model appears to encompass minimal ingredients for thermalization in condensed matter systems including both electrons and phonons. Note that dispersionless phonons do not match the standard properties expected for a thermal reservoir devega_alonso_17 since they exhibit, e.g., infinite autocorrelation times. The results of the present study, however, suggest that the presence of phonons with a dispersion is not a requirement for thermalization and quantum ergodicity in many-body quantum systems.

We first verify that the statistics of neighboring energy-level spacings is well described by the Wigner-Dyson distribution, and that the averages of restricted gap ratios agree with predictions of the Gaussian orthogonal ensemble (GOE). We then study matrix elements of observables in both the electron and phonon sector. For the diagonal matrix elements we use the most restrictive criterion to test the ETH kim_ikeda_14 ; mondaini_fratus_16 , i.e., we show that eigenstate-to-eigenstate fluctuations vanish for all eigenstates in a finite energy-density window. We then present an extensive analysis of the offdiagonal matrix elements and extract the universal function at different energy densities. Finally, we study variances of both diagonal and offdiagonal matrix elements in narrow microcanonical windows and find that their ratio approaches a universal number, which matches predictions of the GOE.

An interesting aspect of the Holstein polaron model is that the electron-phonon coupling is the only mechanism that breaks integrability of the uncoupled electronic and bosonic subsystems and gives rise to quantum ergodicity. Since there is only a single electron in the system, expectation values of the electron-phonon coupling energy are not extensive but of order ${\cal O}(1)$ in every Hamiltonian eigenstate. The question of the minimal perturbation strength to render an integrable system quantum chaotic is highly nontrivial and can be traced back to the pioneering work on ETH by Deutsch deutsch_91 . Recently, a large number of studies addressed the influence of integrability-breaking static impurities on quantum ergodicity and thermalization (see, e.g., Refs. bertini_fagotti_16 ; fagotti_17 ; bastianello_deluca_18b ; bastianello_deluca_18 ; bastianello_19 ), and showed that an ${\cal O}(1)$ integrability-breaking term is enough to induce quantum-chaotic statistics of energy levels santos_04 ; barisic_prelovsek_09 ; santos_mitra_11 ; torresherrera_santos_14 ; brenes_mascarenhas_18 . Here, we provide numerical evidence that an ${\cal O}(1)$-integrability-breaking term, introduced by an itinerant electron coupled to noninteracting phonons, is sufficient to observe perfect ETH properties of all Hamiltonian eigenstates in the bulk of the spectrum.

The paper is organized as follows. We introduce the Holstein polaron model and its basic spectral properties in Sec. II. In Sec. III we study quantum-chaos indicators in the statistics of the neighboring level spacings for different parameter regimes of the model. We then focus on the analysis of the ETH in Sec. IV and explore both diagonal and offdiagonal matrix elements of observables. We conclude in Sec. V.

II The Holstein polaron model

The one-dimensional Holstein polaron model on a lattice with $L$ sites is described by the Hamiltonian

[TABLE]

It consists of the electron kinetic energy operator

[TABLE]

where $\hat{c}_{j}$ is the electron annihilation operator acting on site $j$ (in the Holstein model, a spin-polarized gas of electrons is considered). Next, the phonon-energy operator reads

[TABLE]

where $\hat{b}_{j}$ is the phonon annihilation operator at site $j$ and $\hbar\equiv 1$, and the local coupling of phonons to the electron density $\hat{n}_{j}=\hat{c}_{j}^{{\dagger}}\hat{c}_{j}^{\phantom{{\dagger}}}$ is

[TABLE]

We use periodic boundary conditions, $\hat{c}_{L+j}\equiv\hat{c}_{j}$ and we set the lattice spacing to $a=1$. Throughout the paper we use $t_{0}$ as the energy unit and set $t_{0}=1$ in numerical calculations.

We employ full exact diagonalization to numerically obtain all eigenvalues and eigenstates of the Holstein polaron model, denoted by $\{E_{\alpha}\}$ and $\{|\alpha\rangle\}$, respectively. We truncate the number of local bosonic degrees of freedom, with $M$ representing the maximal number of phonons per site. We exploit translational invariance to split the Hamiltonian into $L$ distinct sectors with quasimomenta $k$. Each sector contains one electron, i.e., $\langle\hat{N}\rangle=\langle\sum_{j}\hat{c}_{j}^{\dagger}\hat{c}_{j}\rangle=1$. The Hilbert-space dimension of each $k$ sector is ${\cal D}^{\prime}=(M+1)^{L}$, and hence the total Hilbert-space dimension is ${\cal D}=L(M+1)^{L}$.

Since we use exact diagonalization, it is necessary to deal with the finite cutoff $M$ in the local phonon number. This introduces a dependence on $M$ in some quantities. The approach taken in this paper is to establish, for a fixed $M$, quantum-chaotic properties of the model and the validity of the ETH for eigenstates that fall within a finite window of energy densities. We will argue that the qualitative conclusions do not depend on the choice of $M$.

We define the distribution of the Hamiltonian eigenvalues $\{E_{\alpha}\}$, i.e., the density of eigenstates, as

[TABLE]

In Fig. 2, we plot $W(E)$ for a system with $L=8$, $M=3$ and the model parameters $\gamma/t_{0}=1/\sqrt{2}$, $\omega_{0}/t_{0}=1/2$. The density of eigenstates is close to a Gaussian function

[TABLE]

which is shown in the same figure (solid line). The average energy $E_{\rm av}$ and the energy variance $\Gamma^{2}$ are given by

[TABLE]

and

[TABLE]

The latter result is consistent with the expectation that the relative standard deviation $\Gamma/L$ vanishes as $1/\sqrt{L}$ when $L\to\infty$. On the other hand, if $L$ is kept fixed and $M\to\infty$, the relative standard deviation $\Gamma/M$ does not vanish but goes to a constant $\Gamma/M\to\omega_{0}\sqrt{L/12}$. We study finite-size dependencies by fixing $M$ and making $L$ larger.

III Quantum-chaos indicators

We first study statistical properties of the spectrum of the Holstein polaron model. A standard approach is to compare these properties with predictions of the random-matrix theory, in particular, to the GOE, which represents the appropriate symmetry class for the Holstein polaron model, i.e., models with time-reversal symmetry. Historically, random-matrix theory was shown to provide a relevant framework to describe spectral properties of quantum systems whose classical counterparts are chaotic bohigas_giannoni_84 . As a consequence, all quantum systems with spectral statistics identical to the ones in random matrices are called quantum chaotic santos_rigol_10a . Since there is no classical counterpart of the Holstein polaron model, we refer to the model as being quantum chaotic if its spectral statistics matches those of the GOE. Note that the level-spacing statistics (to be studied below) can not distinguish between a completely ergodic system and a system with a nonzero number of nonergodic eigenstates whose fraction vanishes in the thermodynamic limit kormos_collura_17 ; turner_michailidis_18a ; turner_michailidis_18b ; ho_choi_19 ; khemani_laumann_19 ; robinson_james_19 ; lin_motrunich_19 ; iadecola_znidaric_19 . In Sec. IV we show that the system under investigation here is indeed completely ergodic, i.e., all eigenstates away from the edges of the spectrum exhibit ETH properties dalessio_kafri_16 .

In finite systems, quantum-chaotic properties are usually observed most unambiguously at intermediate parameter regimes, i.e., when all parameters are of the same order. Moreover, the Holstein polaron model has two integrable limits (the noninteracting limit $\gamma\to 0$ and the small-polaron limit $t_{0}\to 0$), where the quantum-chaotic properties are expected to be absent. Exploring the parameter regimes in which the model exhibits quantum-chaotic properties is one of the main goals of this section.

In the parameter regime where the model is quantum chaotic, this manifests itself for eigenstates at a nonzero energy density above the ground state and below the highest energy eigenstate dalessio_kafri_16 . In the Holstein polaron model, the ground-state energy density in the limit $L\to\infty$ is $E_{\rm min}/E_{\rm av}=0$ and the highest-eigenstate energy density is $E_{\rm max}/E_{\rm av}=2$. We introduce a control parameter $\eta\in[0,1)$ to consider eigenstates in a finite energy-density window, such that

[TABLE]

and

[TABLE]

In our model the corresponding set of eigenstates within the target energy-density window is

[TABLE]

where $\{k\}$ denotes the momentum sectors from which eigenstates are selected. It is well known santos_rigol_10b ; dalessio_kafri_16 that quantum-chaos indicators can only be observed in statistics of Hamiltonian eigenvalues of a single symmetry sector. We therefore limit our analysis in this section to the quasimomentum sector $k=2\pi/L$ in the set of eigenstates ${\cal Z}_{\eta}^{\{k\}}$ in Eq. (11) and fix $\eta=2/3$.

We study quantum-chaos indicators by analyzing the statistics of neighboring energy-level spacings $\delta_{\alpha}=E_{\alpha+1}-E_{\alpha}$. We perform an unfolding of the energy levels in ${\cal Z}_{\eta}^{\{k\}}$ such that the mean level spacing is one brody_flores_81 ; santos_rigol_10a . The unfolding procedure is carried out by introducing a cumulative spectral function ${\cal G}(E)=\sum_{\alpha}\Theta(E-E_{\alpha})$, where $\Theta$ is the unit step function. In the practical analysis, we find it useful to smoothen it by fitting a polynomial of degree six $g_{6}(E)$ to ${\cal G}(E)$. The spectral analysis is then performed on the unfolded neighboring level spacings $s_{\alpha}=g_{6}(E_{\alpha+1})-g_{6}(E_{\alpha})$. We verified that using polynomials $g_{n}(E)$ of different degrees has a negligible influence on the results.

The distribution ${\cal P}(s)$ of the unfolded neighboring level spacings $s_{\alpha}$ for the Holstein polaron model at $\omega_{0}/t_{0}=1/2$ and $\gamma/t_{0}=1/\sqrt{2}$ is shown as a histogram in Fig. 3. We compare the distribution ${\cal P}(s)$ to the Wigner-Dyson distribution (smooth line in Fig. 3)

[TABLE]

which was derived for an ensemble of $2\times 2$ random matrices and, despite not being exact, is now considered as a prototypical distribution of the neighboring level spacings in the GOE dalessio_kafri_16 . We observe a reasonable agreement of ${\cal P}(s)$ with ${\cal P}_{\rm WD}(s)$. The agreement may appear surprising since the only integrability-breaking term in the Holstein polaron model is $\hat{H}_{\rm eph}$ (4), which is intensive. Moreover, our observations are qualitatively consistent with previous studies of level statistics in spin-1/2 chains with integrability-breaking terms of the same order santos_04 ; barisic_prelovsek_09 ; santos_mitra_11 ; torresherrera_santos_14 ; brenes_mascarenhas_18 .

The numerical data shown in Fig. 3 are for the largest numerically available system at $M=3$, and in a regime where the model parameters are of the same order. We now extend our analysis to other parameter regimes and system sizes. Instead of the full distribution ${\cal P}(s)$ we study a single number (defined below) derived from the distribution.

The authors of Ref. oganesyan_huse_07 considered the ratio of two consecutive level spacings (shortly, the gap ratio)

[TABLE]

and introduced the restricted gap ratio

[TABLE]

A convenient property of $\tilde{r}_{\alpha}$ is that no unfolding is necessary to eliminate the influence of a finite-size dependence through the local density of states.

We study $\tilde{r}_{\rm av}=\langle\tilde{r}\rangle_{\eta}$, which represents the average value of the restricted gap ratio $\tilde{r}_{\alpha}$ within the set of eigenstates ${\cal Z}_{\eta}^{k}$ introduced in Eq. (11). A numerically accurate prediction for the average of $\tilde{r}_{\alpha}$ in the GOE is $\tilde{r}_{\rm GOE}\approx 0.5307$ atas_bogomolny_13 . This value is different from the average value of uncorrelated energy levels (relevant for spectra of integrable Hamiltonians with a Poisson distribution of nearest level spacings), which is $\tilde{r}_{\rm uncorr}\approx 0.3863$. The fact that $\tilde{r}_{\rm GOE}>\tilde{r}_{\rm uncorr}$ can be understood from the level repulsion of the GOE, which suppresses the presence of small values of $\tilde{r}_{\alpha}$ in the probability distribution ${\cal P}(\tilde{r})$ compared to the spectra with uncorrelated level spacings (dashed line), as illustrated in Fig. 4(b) dalessio_kafri_16 .

Figure 4(a) shows $\tilde{r}_{\rm av}$ at $\omega_{0}/t_{0}=1/2$ as a function of $\gamma/t_{0}$ for different system sizes $L$. When increasing $L$, the fluctuations of the averages decrease, and at $L=8$, the averages approach the GOE prediction $\tilde{r}_{\rm GOE}$ with high accuracy. Deviations from $\tilde{r}_{\rm GOE}$ are the largest in the limits $\gamma/t_{0}\to 0$ and $\gamma/t_{0}\to\infty$, which is in agreement with observations in other quantum-chaotic Hamiltonians mondaini_rigol_17 .

Before exploring the behavior of $\tilde{r}_{\rm av}$ in other parameter regimes, we verify that $\tilde{r}_{\rm av}$ does not only agree with predictions of the GOE, but that the entire distribution ${\cal P}(\tilde{r})$ of $\tilde{r}_{\alpha}$ does. We plot ${\cal P}(\tilde{r})$ in Fig. 4(b) at $\omega_{0}/t_{0}=1/2$ and $\gamma/t_{0}=1/\sqrt{2}$. The data show an excellent agreement with the analytical GOE prediction

[TABLE]

which is exact for $3\times 3$ matrices and very accurate for asymptotically large systems atas_bogomolny_13 . As a side remark, note that $\tilde{r}_{\rm GOE}=0.5307$ introduced above is a numerical result for asymptotically large systems, which on the third digit differs from the result for $3\times 3$ matrices $\int_{0}^{1}\tilde{r}{\cal P}_{\rm GOE}(\tilde{r}){\rm d}\tilde{r}=4-2\sqrt{3}\approx 0.5359$. Our numerical results in Fig. 4(a) are accurate enough to resolve this difference.

We next study the influence of the phonon energy $\omega_{0}/t_{0}$ on quantum-chaotic properties of the model. Figure 4(c) compares results for $\tilde{r}_{\rm av}$ versus $\gamma/t_{0}$ at $\omega_{0}/t_{0}=1/2$, 1 and 4. As a general trend, fluctuations of $\tilde{r}_{\rm av}$ increase with $\omega_{0}/t_{0}$. This is, in particular, evident for $\omega_{0}/t_{0}=4$, where the agreement with $\tilde{r}_{\rm GOE}$ is weaker. In the antiadiabatic regime of the Holstein polaron model, $\omega_{0}>4t_{0}$, the phonon energy becomes larger than the electron bandwidth and gaps in the spectrum increase. As a consequence, some values of $\tilde{r}_{\alpha}$ in Eq. (14) become very small and the overall value of $\tilde{r}_{\rm av}$ decreases. Since finite-size effects are expected to increase with increasing gaps, we do not study quantum-chaotic properties in the antiadiabatic regime using the existing numerical method.

Finally, we study the influence of the local phonon cutoff $M$ on the restricted gap ratio. Figure 4(d) shows results for $M=1$ (hard-core bosons) at $L=12$, 14 and 16, which leads to identical Hilbert-space dimensions as for the results for $M=3$ and $L=6$, 7, 8, respectively, shown in Fig. 4(a). We observe increased fluctuations when $M$ decreases. Still, fluctuations decrease in both cases when fixing $M$ and increasing $L$, suggesting that in the limit $L\to\infty$, the model exhibits quantum-chaotic properties even in the case of hard-core bosons.

While strictly speaking, the thermodynamic limit corresponds to sending both $L,M\to\infty$, we observe the onset of quantum chaos (Sec. III) and the validity of the ETH (to be studied in Sec. IV) even when $L\to\infty$ and $M$ fixed. A plausible conjecture, which remains to be verified in future work, is that our key observations remain valid when $M\to\infty$. In the following section, we study the ETH properties of the Holstein polaron model using $M=3$, which turned out to be the best compromise between allowing for relatively large local phonon fluctuations and still large enough lattice sizes.

IV Eigenstate thermalization

The previous section studied quantum-chaos indicators and showed that for a wide range of model parameters they agree with predictions of the GOE. We now focus on a specific set of model parameters given by $\omega_{0}/t_{0}=1/2$, $\gamma/t_{0}=1/\sqrt{2}$, $M=3$, and test the validity of the ETH.

The ansatz for expectation values of observables in eigenstates of the Hamiltonian, known also as the ETH ansatz srednicki_99 , can be written as

[TABLE]

where $\bar{E}=(E_{\alpha}+E_{\beta})/2$ is the average energy of the pair of eigenstates $|\alpha\rangle$ and $|\beta\rangle$, $\omega=E_{\alpha}-E_{\beta}$ is the corresponding energy difference, and $S(\bar{E})$ is the thermodynamic entropy at energy $\bar{E}$. The requirement for $O(\bar{E})$ and $f_{O}(\bar{E},\omega)$ is to be smooth functions of their arguments, while $R_{\alpha,\beta}$ are random numbers with zero mean and unit variance. When the ETH ansatz is applicable, $O(\bar{E})$ from the diagonal part is the microcanonical average of an observable $\hat{O}$, while the offdiagonal part can be used, e.g., to prove the fluctuation-dissipation relation for a single eigenstate dalessio_kafri_16 .

We study four different observables, two in the electron sector and two in the phonon sector. In the electron sector, we study the electron kinetic energy $\hat{H}_{\rm kin}$, defined in Eq. (2), and the quasimomentum occupation operator

[TABLE]

at $q=0$. Note that the expectation values of both electron observables are intensive quantities, which is a consequence of having only a single electron in the Holstein polaron model. The kinetic energy $\hat{H}_{\rm kin}$ is the sum over all $\hat{m}_{q}$. Therefore, if ETH behavior is observed for $\hat{m}_{q=0}$ it will not immediately imply ETH behavior for $\hat{H}_{\rm kin}$ (and vice versa). Nevertheless, both quantities are related and it may not be surprising to see ETH behavior in both cases.

In the phonon sector, we study the average phonon number per site $\hat{N}_{\rm ph}=\hat{H}_{\rm ph}/(\omega_{0}L)$, where $\hat{H}_{\rm ph}$ was introduced in Eq. (3), and the nearest-neighbor offdiagonal matrix element of the phonon one-body correlation matrix, with the underlying operator

[TABLE]

The observables in both sectors are chosen such that one observable is part of the Hamiltonian (1) [$\hat{H}_{\rm kin}$ and $\hat{N}_{\rm ph}$] while the other one is not [$\hat{m}_{q}$ and $\hat{T}_{1}$].

IV.1 Diagonal matrix elements of observables

We first focus on the eigenstate-expectation values $\langle\alpha|\hat{O}|\alpha\rangle$, i.e., on the diagonal part of the ETH ansatz in Eq. (16). Figures 5(a)-5(d) show results for observables in the $k=2\pi/L$ quasimomentum sector for three different system sizes $L=6,7,8$, plotted versus the eigenstate energy density $E_{\alpha}/E_{\rm av}$. In all cases we observe clear manifestations of ETH, namely, eigenstate-expectation values are only functions of the energy density. In particular, results for different $L$ suggest that the fluctuations decrease with $L$ and as a consequence, eigenstate-expectation values become smooth and sharp functions of the energy density in the limit $L\to\infty$.

The key step in validating the ETH for the diagonal matrix elements of observables is to quantify finite-size fluctuations of the results in Fig. 5. We study whether (and how) eigenstate-to-eigenstate fluctuations in the bulk of the spectrum vanish when $L\to\infty$. To this end, we use the most restrictive criterion kim_ikeda_14 defined below.

The ETH is expected to be satisfied for the same set of eigenstates in the bulk of the spectrum for which quantum-chaos indicators were found to exhibit GOE behavior in Sec. III. In Eq. (11), we defined these sets of eigenstates ${\cal Z}_{\eta}^{\{k\}}$, where $\eta$ defines the interval of energy densities and $\{k\}$ denotes the set of quasimomentum sectors, from which eigenstates are selected. For the analysis of eigenstate-expectation values, we include eigenstates from quasimomentum sectors $0\leq k\leq\pi$ and denote this set of eigenstates as

[TABLE]

Eigenstates from quasimomentum sectors $k<0$ are excluded from ${\cal Z}_{\eta}$ since the eigenstate-expectation values at $k$ and $-k$ are identical due to time-reversal symmetry. We introduce a measure of eigenstate-to-eigenstate-fluctuations of diagonal expectation values

[TABLE]

and study two statistical properties of $z_{\alpha}(O)$, the mean and the maximum. The mean is defined as

[TABLE]

where $||{\cal Z}_{\eta}||$ is the number of eigenstates in ${\cal Z}_{\eta}$, and the maximum is

[TABLE]

The statistics of eigenstate-to-eigenstate fluctuations, Eq. (20), was first studied in Ref. kim_ikeda_14 . The number of eigenstates analyzed in that work was a finite fraction of the total number of eigenstates, implying that those eigenstates correspond to infinite temperature in the thermodynamic limit. Here, in contrast, we study the statistics of fluctuations in a window of finite energy densities (similar to Ref. mondaini_fratus_16 ), defined by the parameter $\eta$ in Eq. (11). By fixing $\eta$, the fraction $f$ of eigenstates involved in the statistics relative to the total number of eigenstates increases with $L$ and eventually approaches 1 when $L\to\infty$. For example, at $\eta=2/3$ and $L=8$ the fraction $f$ is around 0.9 (see the caption of Fig. 6 for details).

Figures 6(a)-6(d) show the finite-size scaling of $z_{\rm max}$ and $\langle z\rangle_{\eta}$ for the same observables as in Figs. 5(a)-5(d). We show results for the set of eigenstates ${\cal Z}_{\eta}$ using $\eta=1/3$ and $2/3$. Remarkably, in all cases, we observe an exponential decay of fluctuations with the system size $L$, while in fact the electronic density decreases from $1/5$ to $1/8$. These results support the validity of the ETH in the Holstein polaron model, i.e., they suggest that all eigenstates in the bulk of the spectrum obey the ETH.

IV.2 Offdiagonal matrix elements of observables

We now turn our focus to offdiagonal matrix elements $\langle\alpha|\hat{O}|\beta\rangle$ of observables and test the second part of the ETH ansatz in Eq. (16). We consider the symmetry sector $k=2\pi/L$, which includes ${\cal D}^{\prime}=(M+1)^{L}$ eigenstates. We focus on matrix elements of pairs of eigenstates with a similar average energy $\bar{E}$. An example of the set of eigenstates for which $\bar{E}\approx E_{\rm av}$ is sketched as a shaded region in Fig. 7(a).

Figure 8(b) shows results for the matrix elements $|\langle\alpha|\hat{H}_{\rm kin}|\beta\rangle|$ of the kinetic energy as a function of $\omega/E_{\rm av}$ (we only plot results for $\omega=E_{\alpha}-E_{\beta}>0$). We restrict the eigenstates to a narrow energy window $|(E_{\alpha}+E_{\beta})/(2E_{\rm av})-1|<\Delta/2$ using $\Delta=10^{-3}$. The results reveal strong fluctuations of $|\langle\alpha|\hat{H}_{\rm kin}|\beta\rangle|$, which is in agreement with the presence of the random term $R_{\alpha,\beta}$ in the ETH ansatz (16). The moving average shown in the same panel indicates that the averaged function is almost flat at small $\omega/E_{\rm av}$ and then rapidly decays at larger $\omega/E_{\rm av}$, consistent with results observed in the two-dimensional transverse field Ising model mondaini_rigol_17 and in the hard-core boson model with dipolar interactions khatami_pupillo_13 .

The overall prefactor of the matrix elements $|\langle\alpha|\hat{O}|\beta\rangle|$ depends on system size. To extract that prefactor we define the average offdiagonal matrix element

[TABLE]

where the normalization ${\cal N}$ equals the number of elements included in the sum and $\bar{\varepsilon}=\bar{E}/E_{\rm av}$ is the target average energy density of the pair of eigenstates. The ETH ansatz predicts the prefactor to be $\sim e^{-S(\bar{E})/2}$, for which the leading term at $\bar{\varepsilon}=1$ scales as $\sim 1/\sqrt{{\cal D}^{\prime}}$. Figure 8(a) shows our numerical results for $\overline{|O_{\alpha,\beta}|}$ at $\bar{\varepsilon}=1$ as a function of ${\cal D}^{\prime}$. We use $\Delta=10^{-3}$ and $10^{-1}$ in Eq. (23) and show that results for different $\Delta$ are almost identical. The data points can, for all observables, accurately be fitted using the function $a({\cal D}^{\prime})^{-b}$. We find that $b$ is very close to the predicted value $b=1/2$ (see the caption of Fig. 8 for the actual numerical values of $b$). However, if one also aims at resolving eventual subleading corrections to the scaling controlled by Eq. (16), as discussed in Ref. luitz_barlev_16 , more data points would be needed.

The universal properties of the offdiagonal matrix elements $|\langle\alpha|\hat{O}|\beta\rangle|$ encoded in $f_{O}(\bar{E},\omega)$ can be studied after randomness and system-size dependent prefactors are removed. We define

[TABLE]

which is a moving average of the renormalized offdiagonal matrix elements of observables. The function ${\cal F}_{O}(\bar{E},\omega)$ is, up to a constant prefactor, identical to the universal function $f_{O}(\bar{E},\omega)$ used in the ETH ansatz (16). We therefore plot ${\cal F}_{O}(\bar{E},\omega)$ in Fig. 8(c) for different observables as a function of $\omega/E_{\rm av}$. The results reveal that the offdiagonal matrix elements of both electron and phonon observables exhibit a similar scaling with $\omega$.

We complement our previous results by calculating the offdiagonal matrix elements of the pairs of eigenstates at average energy densities $\bar{\varepsilon}<1$. In Fig. 9, we show results for two observables (the electron kinetic energy and the phonon correlator), while results for the other two observables studied in Fig. 8(c) exhibit a similar behavior (not shown here). According to the ETH ansatz (16), the scaling of $|\langle\alpha|\hat{O}|\beta\rangle|$ with the system size is governed by $e^{-S(\bar{E})/2}$. In the context of Eq. (23), $S(\bar{E})$ is the logarithm of the number of microstates $N_{\bar{\varepsilon}}$ in the energy density window $[\bar{\varepsilon}-\Delta/2,\bar{\varepsilon}+\Delta/2]$. In Fig. 9(a), we plot $|\langle\alpha|\hat{O}|\beta\rangle|$ versus $N_{\bar{\varepsilon}}$ and observe a reasonably good collapse of the data for different $\bar{\varepsilon}$. We fit the results to $a(N_{\bar{\varepsilon}})^{-b}$ and obtain $b$ close to 0.5, as predicted by the ETH ansatz.

On the other hand, the dependence of $f_{O}(\bar{E},\omega)$ on the average energy $\bar{E}$ is not described by the ETH ansatz and rarely studied in the literature (see Ref. dalessio_kafri_16 for a discussion on how $f_{O}(\bar{E},\omega)$ is related to the fluctuation-dissipation relation). Our results in Fig. 9(b) suggest that ${\cal F}_{O}(\bar{E},\omega)$ in the Holstein polaron model is roughly independent of the pair energy average $\bar{\varepsilon}=\bar{E}/E_{\rm av}$. We note that a recent study ad hoc assumed a slowly varying $f_{O}(\bar{E},\omega)$ with $\bar{E}$ in a nonintegrable spin chain richter_gemmer_18 . Our results lend numerical support to this assumption. The study of $f_{O}(\bar{E},\omega)$ and its $\bar{E}$-dependence in other models remains an open problem.

IV.3 Variances of matrix elements of observables

Finally, we study fluctuations of both diagonal and offdiagonal matrix elements of observables, focusing on the quasimomentum sector $k=2\pi/L$. We define the variance of a set of consecutive diagonal matrix elements

[TABLE]

where $\alpha$ determines the first matrix element and $\mu$ determines the number of matrix elements in the set. The average $\langle O_{\alpha,\alpha}\rangle_{\mu}=\mu^{-1}\sum_{\rho=\alpha}^{\alpha+\mu-1}\langle\rho|\hat{O}|\rho\rangle$ is defined in a microcanonical window, i.e., it is an equal-weight average over $\mu$ consecutive eigenstates, starting at the index $\alpha$. Similarly, the variance of the submatrix of offdiagonal matrix elements is defined as

[TABLE]

where $\langle O_{\alpha,\beta}\rangle_{\mu}=(\mu^{2}-\mu)^{-1}\sum_{\begin{subarray}{c}\rho,\rho^{\prime}=\alpha\end{subarray}}^{\alpha+\mu-1}\langle\rho^{\prime}|\hat{O}|\rho\rangle$ and $\rho\neq\rho^{\prime}$. Note that the offdiagonal matrix elements can be complex as a consequence of diagonalizing the Hamiltonian in the translationally invariant basis. Both variances are sketched in Fig. 7(b) for a given $\alpha$ and $\mu$. We associate the mean energy $e_{\alpha,\mu}$ to every microcanonical window defined by ($\alpha,\mu$), where $e_{\alpha,\mu}=\langle H_{\alpha,\alpha}\rangle_{\mu}$.

In our calculations, we fix $\mu$ and calculate the variances for every $\alpha$ in the interval $[1,{\cal D}^{\prime}-\mu]$. Our main focus is on the ratio of variances defined as

[TABLE]

Within random-matrix theory one can show dalessio_kafri_16 that the ratio of variances takes a universal value $\Sigma_{\rm GOE}^{2}(O)=2$ for generic local observables in the GOE. Since the ETH generalizes random-matrix theory to Hamiltonians of real physical systems, one may ask whether the same result can be obtained for matrix elements of Hamiltonian eigenstates prosen_99 and if the answer is affirmative, at which energy densities this behavior is found.

The numerical extraction of $\Sigma_{\alpha,\mu}^{2}$ obeying predictions from the random-matrix theory is a nontrivial task, and the results depend on the parameters $\alpha$ and $\mu$ as well as on the Hilbert-space dimension ${\cal D}^{\prime}$. In the main panels of Figs. 10(a)-10(d), we fix $\mu=100$ and set $L=8$, $M=3$, which yields ${\cal D^{\prime}}=65536$ and hence $\mu\ll{\cal D}^{\prime}$. Our choice of $\mu$ is consistent with the one used in a recent study of a two-dimensional spin-1/2 system mondaini_rigol_17 . We compute $\Sigma_{\alpha,\mu=100}^{2}(O)$ for four different observables as a function of the eigenstate index $\alpha$ and observe $\Sigma_{\alpha,\mu=100}^{2}(O)\approx 2$ for the majority of eigenstates for all observables.

We quantify the excellent agreement with predictions of random-matrix theory, observed in Figs. 10(a)-10(d), by computing the averages over those ratios of variances for which $e_{\alpha,\mu=100}/E_{\rm av}\in[1-\eta,1+\eta]$ and $\eta=2/3$. The resulting averages are shown as horizontal dashed lines in Fig. 10 and yield, for all four observables, almost perfect agreement with $\Sigma_{\rm GOE}^{2}(O)$. The insets of Figs. 10(a)-10(d) show the distributions of $\Sigma_{\alpha,\mu=100}^{2}(O)$ for the same set of $\alpha$ as the one that was used to obtain the horizontal lines in the main panel, and two system sizes $L=7$ and 8. They reveal that, by increasing the system size, the distributions become narrower and their means approach the random-matrix theory result.

It is remarkable that the agreement of the ratio of variances for Hamiltonian eigenstates with the GOE predictions carries over to eigenstates away from the middle of the spectrum. Namely, in contrast to eigenstates in the middle of the spectrum, which can be well approximated by pure states that are random superpositions of base kets in some simple basis vidmar_rigol_17 , eigenstates away from the middle of the spectrum are no longer entirely random. Figure 11(a) shows the same data for $\Sigma_{\alpha,\mu=100}^{2}$ as in Fig. 10(d), but plotted versus the mean energy density $e_{\alpha,\mu=100}/E_{\rm av}$. It reveals a broad region away from $e_{\alpha,\mu=100}/E_{\rm av}=1$ where $\Sigma_{\alpha,\mu=100}^{2}\approx 2$.

Finally, we study the dependence of $\Sigma_{\alpha,\mu}^{2}$ on $\mu$ and the system size, and we only consider eigenstates $\alpha$ away from the middle of the spectrum. To this end, we consider a subset of ratios of variances ${\cal S}_{\mu}=\{\Sigma_{\alpha,\mu}^{2}\;;\;e_{\alpha,\mu}/E_{\rm av}\in[1/3,2/3]\}$, which corresponds to a shaded region in Fig. 11(a). The average $\Sigma_{\alpha,\mu}^{2}$ in this region is defined as

[TABLE]

and plotted in Fig. 11(b). The results reveal that when $L$ increases, a larger number of eigenstates $\mu$ can be included in the variances to observe $\Sigma_{\alpha,\mu}^{2}\to\Sigma_{\rm GOE}^{2}$. Remarkably, this suggests that $\Sigma_{\alpha,\mu}^{2}=\Sigma_{\rm GOE}^{2}$ in the limit $L\to\infty$, for eigenstates both in the middle and away from the center of the spectrum.

V Conclusions

In this paper, we addressed the question whether the hallmark features of the ETH can be observed in a paradigmatic condensed matter model that includes both electron and phonon degrees of freedom. We studied the Holstein polaron model, which is a minimal model describing a single electron that locally couples to dispersionless phonons. In spite of the apparent simplicity of the model, we showed that it represents an example where ETH is very well fulfilled, to a quantitative degree that is observed in single-component models as well mondaini_rigol_17 . In particular, we showed that: (i) the statistics of neighboring energy-level spacings agree with predictions of the GOE; (ii) the eigenstate-to-eigenstate fluctuations of the diagonal matrix elements of observables decay exponentially for all eigenstates in the bulk of the spectrum; (iii) the ’universal’ part of the offdiagonal matrix elements of all four observables under investigation exhibits a similar decay as a function of the energy difference of pairs of eigenstates, and exhibits an almost negligible dependence on the average energy of the pair of eigenstates; and (iv) the ratio of variances of fluctuations of diagonal to offdiagonal matrix elements of observables fulfills predictions of the random-matrix theory. The latter is a robust feature of Hamiltonian eigenstates also away from the middle of the spectrum.

Our results establish two important aspects of quantum many-body systems: The first is to demonstrate that quantum ergodicity and thermalization in a many-body system of electrons and phonons does not require phonons to exhibit properties of a bath in isolation, as usually assumed in studies of open quantum systems devega_alonso_17 . The second is to demonstrate that an integrability breaking term of the order ${\cal O}(1)$ is sufficient to induce perfect ETH behavior. [In our case, the eigenstate expectation value of the electron-phonon coupling operator is ${\cal O}(1)$ in every eigenstate.] This extends previous studies santos_04 ; barisic_prelovsek_09 ; santos_mitra_11 ; torresherrera_santos_14 ; brenes_mascarenhas_18 that showed that an ${\cal O}(1)$ integrability breaking term renders the spectral level statistics to be quantum chaotic.

There are several interesting extensions to our paper. First, one may study generic conditions for an integrability breaking term of the order ${\cal O}(1)$ to induce ETH behavior. Of particular interest may be systems where the uncoupled electron-phonon system is not as highly degenerate as in the Holstein model. This can be achieved, e.g., by considering dispersive or disordered phonons. Second, extending the analysis to finite electronic densities would be interesting, but also numerically much more challenging. Third, exact diagonalization studies are obviously limited to small systems. Therefore, complementary analytical insights into the problem of ETH in electron-phonon models would be desirable. Finally, concrete studies of quench problems in the polaron case with an extensive quench energy are lacking. This could be carried out by using, e.g., matrix-product-state (MPS) methods schollwock2005density ; schollwock2011density ; paeckel_koehler_19 (see Refs. znidaric_scardicchio_16 ; kloss_barlev_18 ; doggen_schindler_18 ; hallam_morley_19 for examples of recent developments). References brockt_dorfner_15 ; schroeder_chin_16 ; wall_safavinaini_16 ; kloss_reichman_19 discuss MPS methods that are specifically designed for electron-phonon coupled systems. These questions are left for future work.

Acknowledgments

We acknowledge useful discussions with J. Bonča, I. de Vega, M. Goldstein, D. J. Luitz, S. Kehrein, M. Mierzejewski, P. Prelovšek, M. Rigol and R. Steinigeweg. This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Project No. 207383564 via Research Unit FOR 1807, and via CRC 1073 under Project No. 217133147. D.J. acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy No. EXC-2111-390814868. L.V. acknowledges support from the Slovenian Research Agency (ARRS), Research Core Funding No. P1-0044. Part of this research was conducted at KITP at the University of California at Santa Barbara. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958.