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 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.


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 [1,2]. This topic became viral after quantum-gas experiments [3][4][5] 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 [6][7][8][9] or the exceptions such as integrable models [10][11][12] or many-body localization [13][14][15].
An additional impetus to the theoretical community was given by the work of Rigol and coworkers [16] who demonstrated that the eigenstate thermalization hypothesis (ETH), pioneered by Deutsch [17] and Srednicki [18], provides a relevant framework to describe statistical properties of eigenstates of lattice Hamiltonians, applicable to ongoing experiments with ultracold atoms on optical lattices [7,13,14]. 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 [19][20][21][22].
So far, the ETH has been verified for a wide number of lattice models such as nonintegrable spin-1/2 chains [23][24][25][26][27][28][29][30][31][32][33], ladders [26,[34][35][36] and square lattices [37][38][39], interacting spinless fermions [40,41], Bose-Hubbard [26,42] and Fermi-Hubbard chains [43], dipolar hard-core bosons [44], quantum dimer models [45] and Fibonacci anyons [46]. 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 thermal-ization 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.
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 [69] 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 energylevel 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. 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 [27,38], i.e., we show that eigenstateto-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 Gaussian orthogonal ensemble.
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 operators are not extensive but of order 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 [17]. Recently, a large number of studies addressed the influence of integrability-breaking static impurities on quantum ergodicity and thermalization (see, e.g., Refs. [70][71][72][73]), and showed that an O(1) integrability-breaking term is enough to induce quantum-chaotic statistics of energy levels [74][75][76][77][78]. Here, we provide numerical evidence that an 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 It consists of the electron kinetic energy operator whereĉ 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Ĥ whereb j is the phonon annihilation operator at site j and ≡ 1, and the local coupling of phonons to the electron densityn j =ĉ † jĉj iŝ We use periodic boundary conditions,ĉ L+j ≡ĉ 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 Holsteinpolaron model, denoted by {E α } and {|α }, 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., N = jĉ † jĉ j = 1. The Hilbert-space dimension of each k sector is D = (M + 1) L , and hence the total Hilbert-space dimension is 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 We define the distribution of the Hamiltonian eigenvalues {E α }, i.e., the density of eigenstates, as In Fig. 2, we plot W (E) for a system with L = 8, M = 3 and the model parameters γ/t 0 = 1/ √ 2, ω 0 /t 0 = 1/2. The density of eigenstates is close to a Gaussian function which is shown in the same figure (solid line). The average energy E av and the energy variance Γ 2 are given by and The latter result is consistent with the expectation that the relative standard deviation Γ/L vanishes as 1/ √ L when L → ∞. On the other hand, if L is kept fixed and M → ∞, the relative standard deviation Γ/M does not vanish but goes to a constant Γ/M → ω 0 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 randommatrix theory, in particular, to the Gaussian orthogonal ensemble (GOE), which represents the appropriate symmetry class for the Holstein-polaron model, i.e., models with time-reversal symmetry. Historically, randommatrix theory was shown to provide a relevant framework to describe spectral properties of quantum systems whose classical counterparts are chaotic [79]. As a consequence, all quantum systems with spectral statistics identical to the ones in random matrices are called quantum chaotic [80]. 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 and the eigenstates exhibit ETH properties [20].
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 γ → 0 and the small-polaron limit t 0 → 0), where the quantum-chaotic properties are expected to be absent. Exploring the parameter regimes in which the model exhibits quantumchaotic 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 [20]. In the Holsteinpolaron model, the ground-state energy density in the limit L → ∞ is E min /E av = 0 and the highest-eigenstate energy density is E max /E av = 2. We introduce a control parameter η ∈ [0, 1) to consider eigenstates in a finite energy-density window, such that and In our model the corresponding set of eigenstates within the target energy-density window is where {k} denotes the momentum sectors from which eigenstates are selected. It is well known [20,81] 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π/L in the set of eigenstates Z {k} η in Eq. (11) and fix η = 2/3. We study quantum-chaos indicators by analyzing the statistics of neighboring energy-level spacings δ α = E α+1 − E α . We perform an unfolding of the energy levels in Z {k} η such that the mean level spacing is one. The unfolding procedure is carried out by introducing a cumulative spectral function G(E) = α Θ(E − E α ), where Θ is the unit step function. In the practical analysis, we find it useful to smoothen it by fitting a polynomial Level-spacing statistics P(s) after the unfolding procedure (see the text for details). (a) P(s) for the Holsteinpolaron model, Eq. (1), for L = 8, M = 3, ω0/t0 = 0.5 and (13), for L = 24 and 8 particles. In both cases, numerical results are shown as histograms for the k = 2π/L momentum sector. In (a) we consider eigenstates within an energy density window given by η = 2/3 in Eq. of degree six g 6 (E) to G(E). The spectral analysis is then performed on the unfolded neighboring level spacings s α = g 6 (E α+1 ) − g 6 (E α ). We verified that using polynomials g n (E) of different degrees has a negligible influence on the results.
The distribution P(s) of the unfolded neighboring level spacings s α for the Holstein-polaron model at ω 0 /t 0 = 1/2 and γ/t 0 = 1/ √ 2 is shown as histograms in Fig. 3(a). We compare the distribution P(s) to the Wigner-Dyson distribution (smooth lines in Fig. 3) which was derived for an ensemble of 2 × 2 random matrices and, despite not being exact, is now considered as a prototypical distribution of the neighboring level spacings in the GOE [20]. We observe a reasonable agreement of P(s) with P WD (s). In Fig. 3(b), we also show the distribution P(s) of the so-called t-t -V -V model, which is a model of hard-core bosons on a one-dimensional lattice with periodic boundaries, We impose the hard-core boson constraint (â l ) 2 = (â † l ) 2 = 0. Note that the nonzero values of the nextnearest neighbor hopping t and the interaction strength V are crucial to break integrability of the model. In our calculations, we set t = t = 1 and V = V = 1.1. For these parameters, this model has been shown [23,80,82] to be quantum chaotic and to exhibit ETH behavior. The comparison of the results shown in Figs. 3(a) and 3(b) unveils that the level-spacing distribution P(s) of the Holstein-polaron model is, already in finite systems with a Hilbert-space dimensions of the order ∼ 10 4 −10 5 , almost indistinguishable from the P(s) of the prototypical quantum-chaotic model introduced in Eq. (13). While the agreement may appear surprising since the only integrability-breaking term in the Holstein-polaron model isĤ eph (4), which is of the order O(1), our observations are qualitatively consistent with previous studies of level statistics in spin-1/2 chains with integrabilitybreaking terms of the same order [74][75][76][77][78].
The numerical data shown in Fig. 3(a) 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 P(s) we study a single number (defined below) derived from the distribution.
The authors of Ref. [84] considered the ratio of two consecutive level spacings (shortly, the gap ratio) and introduced the restricted gap ratiõ A convenient property ofr α is that no unfolding is necessary to eliminate the influence of a finite-size dependence through the local density of states. We studyr av = r η , which represents the average value of the restricted gap ratior α within the set of eigenstates Z k η introduced in Eq. (11). A numerically accurate prediction for the average ofr α in the GOE is r GOE ≈ 0.5307 [83]. 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 isr uncorr ≈ 0.3863. The fact thatr GOE >r uncorr can be understood from the level repulsion of the GOE, which suppresses the presence of small values ofr α in the probability distribution P(r) compared to the spectra with uncorrelated level spacings (dashed line), as illustrated in Fig. 4(b) [20]. Figure 4(a) showsr av at ω 0 /t 0 = 1/2 as a function of γ/t 0 for different system sizes L. It exhibits two important features: (i) when increasing L, the fluctuations above the averages decrease, and at L = 8, the averages approach the GOE predictionr GOE with high accuracy; (ii) deviations fromr GOE are enhanced in the limits γ/t 0 → 0 and γ/t 0 → ∞. However, when L increases the onset of deviations is shifted to smaller γ/t 0 in the regime γ/t 0 → 0, and to larger γ/t 0 in the regime γ/t 0 → ∞. This suggests that eventually,r av →r GOE for all nonzero and finite γ/t 0 in the thermodynamic limit, which is in agreement with the finite-size dependencies observed in other quantum-chaotic Hamiltonians [39].
Before exploring the behavior ofr av in other parameter regimes, we verify thatr av does not only agree with predictions of the GOE, but that the entire distribution P(r) ofr α does. We plot P(r) in Fig. 4(b) at ω 0 /t 0 = 1/2 and γ/t 0 = 1/ √ 2. The data show an excellent agreement with the analytical GOE prediction P GOE (r) = 27 4r which is exact for 3 × 3 matrices and very accurate for asymptotically large systems [83]. As a side remark, note thatr GOE = 0.5307 introduced above is a numerical result for asymptotically large systems, which on the third digit differs from the result for 3 × 3 matrices 1 0r P GOE (r)dr = 4 − 2 √ 3 ≈ 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 ω 0 /t 0 on quantum-chaotic properties of the model. Figure 4(c) compares results forr av versus γ/t 0 at ω 0 /t 0 = 1/2, 1 and 4. As a general trend, fluctuations ofr av increase with ω 0 /t 0 . This is, in particular, evident for ω 0 /t 0 = 4, where the agreement withr GOE is weaker. In the antiadiabatic regime of the Holstein-polaron model, ω 0 > 4t 0 , the phonon energy becomes larger than the electron bandwidth and gaps in the spectrum increase. As a consequence, some values ofr α in Eq. (15) become very small and the overall value ofr 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 → ∞, 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 → ∞, we observe the onset of quantum chaos (Sec. III) and the validity of the ETH (to be studied in Sec. IV) even when L → ∞ and M fixed. We hence conjecture that the statements made for the latter case remain valid when M → ∞. 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 ω 0 /t 0 = 1/2, γ/t 0 = 1/ √ 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 [85], can be written as whereĒ = (E α + E β )/2 is the average energy of the pair of eigenstates |α and |β , ω = E α − E β is the corresponding energy difference, and S(Ē) is the thermodynamic entropy at energyĒ. The requirement for O(Ē) and f O (Ē, ω) is to be smooth functions of their arguments, while R α,β are random numbers with zero mean and unit variance. When the ETH ansatz is applicable, O(Ē) from the diagonal part is the microcanonical average of an observableÔ, while the offdiagonal part can be used, e.g., to prove the fluctuation-dissipation relation for a single eigenstate [20].  . We include eigenstates from the quasimomentum sectors with 0 ≤ k ≤ π that belong to the subspace Zη, defined in Eq. (20). At η = 1/3, the fraction of eigenstates included in Zη (out of the total number of eigenstates from the quasimomentum sectors 0 ≤ k ≤ π) are 37%, 44%, 50% and 55% for L = 5, 6, 7, 8, respectively, while at η = 2/3 the fractions are 70%, 79%, 85% and 90%, respectively. Solid lines are guides to the eyes, signaling an exponential decay with L.
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Ĥ kin , defined in Eq. (2), and the quasimomentum occupation operator 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. In the phonon sector, we study the average phonon number per siteN ph =Ĥ ph /(ω 0 L), whereĤ ph was introduced in Eq. (3), and the nearest-neighbor offdiagonal matrix element of the phonon one-body correlation matrix, with the underlying operator The observables in both sectors are chosen such that one observable is part of the Hamiltonian (1) [Ĥ kin andN ph ] while the other one is not [m q andT 1 ].

A. Diagonal matrix elements of observables
We first focus on the eigenstate-expectation values α|Ô|α , i.e., on the diagonal part of the ETH ansatz in Eq. (17). Figures 5(a)-5(d) show results for observables in the k = 2π/L quasimomentum sector for three different system sizes L = 6, 7, 8, plotted versus the eigenstate energy density E α /E 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, eigenstateexpectation values become smooth and sharp functions of the energy density in the limit L → ∞.
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 this end, we use the most restrictive criterion [27] 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 Z {k} η , where η 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 ≤ k ≤ π and denote this set of eigenstates as Eigenstates from quasimomentum sectors k < 0 are excluded from Z η since the eigenstate-expectation values at k and −k are identical due to time-reversal symmetry. We introduce a measure of eigenstate-to-eigenstatefluctuations of diagonal expectation values and study two statistical properties of z α (O), the mean and the maximum. The mean is defined as where ||Z η || is the number of eigenstates in Z η , and the maximum is The statistics of eigenstate-to-eigenstate fluctuations, Eq. (21), was first studied in Ref. [27]. 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. [38]), defined by the parameter η in Eq. (11). By fixing η, 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 → ∞. For example, at η = 2/3 and L = 8 the fraction f is around 0.9 (see the caption of Fig. 6 for details). We show results for the set of eigenstates Z η using η = 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.

B. Offdiagonal matrix elements of observables
We now turn our focus to offdiagonal matrix elements α|Ô|β of observables and test the second part of the ETH ansatz in Eq. (17). We consider the symmetry sector k = 2π/L, which includes D = (M + 1) L eigenstates. We focus on matrix elements of pairs of eigenstates with a similar average energyĒ. An example of the set of eigenstates for whichĒ ≈ E av is sketched as a shaded region in Fig. 7(a). Figure 8(b) shows results for the matrix elements | α|Ĥ kin |β | of the kinetic energy as a function of ω/E av (we only plot results for ω = E α − E β > 0). We restrict the eigenstates to a narrow energy window |(E α + E β )/(2E av ) − 1| < ∆/2 using ∆ = 10 −3 . The results reveal strong fluctuations of | α|Ĥ kin |β |, which is in agreement with the presence of the random term R α,β in the ETH ansatz (17). The moving average shown in the same panel indicates that the averaged function is almost flat at small ω/E av and then rapidly decays at larger ω/E av , consistent with results observed in the two-dimensional transverse field Ising model [39] and in the hard-core boson model with dipolar interactions [44].
The overall prefactor of the matrix elements | α|Ô|β | depends on system size. To extract that prefactor we define the average offdiagonal matrix element where the normalization N equals the number of elements included in the sum andε =Ē/E av is the target average energy density of the pair of eigenstates. The ETH ansatz predicts the prefactor to be ∼ e −S(Ē)/2 , for which the leading term atε = 1 scales as ∼ 1/ √ D . Figure 8(a) shows our numerical results for |O α,β | atε = 1 as a function of D . We use ∆ = 10 −3 and 10 −1 in Eq. (24) and show that results for different ∆ are almost identical. The data points can, for all observables, accurately be fitted using the function a(D ) −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 correc- tions to the scaling controlled by Eq. (17), as discussed in Ref. [29], more data points would be needed.
The universal properties of the offdiagonal matrix elements | α|Ô|β | encoded in f O (Ē, ω) can be studied after randomness and system-size dependent prefactors are removed. We define which is a moving average of the renormalized offdiagonal matrix elements of observables. The function F O (Ē, ω) is, up to a constant prefactor, identical to the universal function f O (Ē, ω) used in the ETH ansatz (17). We therefore plot F O (Ē, ω) in Fig. 8(c) for different observables as a function of ω/E av . The results reveal that the offdiagonal matrix elements of both electron and phonon observables exhibit a similar scaling with ω. We complement our previous results by calculating the offdiagonal matrix elements of the pairs of eigenstates at average energy densitiesε < 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 (17), the scaling of | α|Ô|β | with the system size is governed by e −S(Ē)/2 . In the context of Eq. (24), S(Ē) is the logarithm of the number of microstates Nε in the energy density window [ε − ∆/2,ε + ∆/2]. In Fig. 9(a), we plot | α|Ô|β | versus Nε and observe a reasonably good collapse of the data for differentε. We fit the results to a(Nε) −b and obtain b close to 0.5, as predicted by the ETH ansatz.
On the other hand, the dependence of f O (Ē, ω) on the average energyĒ is not described by the ETH ansatz and rarely studied in the literature (see Ref. [20] for a discussion on how f O (Ē, ω) is related to the fluctuationdissipation relation). Our results in Fig. 9(b) suggest that F O (Ē, ω) in the Holstein polaron model is roughly independent of the pair energy averageε =Ē/E av . We note that a recent study ad-hoc assumed a slowly varying f O (Ē, ω) withĒ in a nonintegrable spin chain [31]. Our results lend numerical support to this assumption. The study of f O (Ē, ω) and itsĒ-dependence in other models remains an open problem.

C. 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π/L. We define the vari- ance of a set of consecutive diagonal matrix elements where α determines the first matrix element and µ determines the number of matrix elements in the set. The average O α,α µ = µ −1 α+µ−1 ρ=α ρ|Ô|ρ is defined in a microcanonical window, i.e., it is an equal-weight average over µ consecutive eigenstates, starting at the index α. Similarly, the variance of the submatrix of offdiagonal matrix elements is defined as where O α,β µ = (µ 2 − µ) −1 α+µ−1 ρ,ρ =α ρ |Ô|ρ and ρ = ρ . Note that the offdiagonal matrix elements can be complex as a consequence of diagonalizing the Hamiltonian in the translationally invariant basis. Both variances are  Fig. 10(d), but plotted versus the mean energy density eα,µ=100/Eav. (b) Average ratio of variances Σ 2 µ , defined in Eq. (29), as a function of µ 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 × 10 2 , 3.0 × 10 3 and 1.2 × 10 4 for L = 6, 7, 8, respectively. sketched in Fig. 7(b) for a given α and µ. We associate the mean energy e α,µ to every microcanonical window defined by (α, µ), where e α,µ = H α,α µ .
In our calculations, we fix µ and calculate the variances for every α in the interval [1, D − µ]. Our main focus is on the ratio of variances defined as Within random-matrix theory one can show [20] that the ratio of variances takes a universal value Σ 2 GOE (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 and if the answer is affirmative, at which energy densities this behavior is found.
The numerical extraction of Σ 2 α,µ obeying predictions from the random-matrix theory is a nontrivial task, and the results strongly depend on the parameters α and µ as well as on the Hilbert-space dimension D . In the main panels of Figs. 10(a)-10(d), we fix µ = 100 and set L = 8, M = 3, which yields D = 65536 and hence µ D . Our choice of µ is consistent with the one used in a recent study of a two-dimensional spin-1/2 system [39]. We compute Σ 2 α,µ=100 (O) for four different observables as a function of the eigenstate index α and observe Σ 2 α,µ=100 (O) ≈ 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 α,µ=100 /E av ∈ [1 − η, 1 + η] and η = 2/3. The resulting averages are shown as horizontal dashed lines in Fig. 10 and yield, for all four observables, almost perfect agreement with Σ 2 GOE (O). The insets of Figs. 10(a)-10(d) show the distributions of Σ 2 α,µ=100 (O) for all α and two system sizes L = 7 and 8. They reveal that, by increasing the system size, the distributions become narrower and their mean approaches 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 superposition of base kets in some simple basis [82], eigenstates away from the middle of the spectrum are no longer entirely random. Figure 11(a) shows the same data for Σ 2 α,µ=100 as in Fig. 10(d), but plotted versus the mean energy density e α,µ=100 /E av . It reveals a broad region away from e α,µ=100 /E av = 1 where Σ 2 α,µ=100 ≈ 2. Finally, we study the dependence of Σ 2 α,µ on µ and the system size, and we only consider eigenstates α away from the middle of the spectrum. To this end, we consider a subset of ratios of variances S µ = {Σ 2 α,µ ; e α,µ /E av ∈ [1/3, 2/3]}, which corresponds to a shaded region in Fig. 11(a). The average Σ 2 α,µ in this region is defined as and plotted in Fig. 11(b). The results reveal that when L increases, a larger number of eigenstates µ can be included in the variances to observe Σ 2 α,µ → Σ 2 GOE . Remarkably, this suggest that Σ 2 α,µ = Σ 2 GOE in the limit L → ∞, which is, to our knowledge, the first observation of this property for Hamiltonian eigenstates well below the middle of the spectrum.

V. CONCLUSIONS
In this work, 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 actually is a canonical example of the ETH, to a quantitative degree that is observed in single-component models as well [39]. In particular, we showed that: (i) the statistics of neighboring energy-level spacings agree with predictions of the Gaussian orthogonal ensemble; (ii) the eigenstate-toeigenstate 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 one is to demonstrate that quantum ergodicity and thermalization in a manybody 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 [69]. The second one is to demonstrate that an integrability breaking term of the order O(1) is sufficient to induce perfect ETH behavior. [In our case, the eigenstate expectation value of the electron-phonon coupling operator is O(1) in every eigenstate.] This extends previous studies [74][75][76][77][78] that showed that an O(1) integrability breaking term renders the spectral level statistics to be quantum chaotic.
There are several interesting extensions to our work. First of all, one may study generic conditions for an integrability breaking term of the order 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. These questions are left for future work.