The eigenstate thermalization hypothesis beyond standard indicators: emergence of random-matrix behavior at small frequencies

Using numerical exact diagonalization, we study matrix elements of a local spin operator in the eigenbasis of two different nonintegrable quantum spin chains. Our emphasis is on the question to what extent local operators can be represented as random matrices and, in particular, to what extent matrix elements can be considered as uncorrelated. As a main result, we show that the eigenvalue distribution of band submatrices at a fixed energy density is a sensitive probe of the correlations between matrix elements. We find that, on the scales where the matrix elements are in a good agreement with all standard indicators of the eigenstate thermalization hypothesis (ETH), the eigenvalue distribution still exhibits clear signatures of the original operator, implying correlations between matrix elements. Moreover, we demonstrate that at much smaller energy scales, the eigenvalue distribution approximately assumes the universal semicircle shape, indicating transition to the random-matrix behavior, and in particular that matrix elements become uncorrelated.


I. INTRODUCTION
Questions of equilibration and thermalization in isolated quantum many-body systems have experienced an upsurge of interest both from the theoretical and the experimental side over the last decades [1][2][3].In this context, the eigenstate thermalization hypothesis (ETH) has been established as a key concept to explain the emergence of thermodynamic behavior, by assuming a certain matrix structure of physical operators O in the eigenbasis of generic Hamiltonians H [4][5][6].Specifically, let O mn = m| O |n denote the matrix element of O within the eigenstates |m and |n of H, then the ETH ansatz reads [3,7] O mn = O( Ē)δ mn + Ω − 1 2 ( Ē)f ( Ē, ω)r mn , where ω = E m −E n is the difference between the eigenenergies E m and E n with mean energy Ē = (E m + E n )/2, O( Ē) and f ( Ē, ω) are smooth functions of their arguments, and Ω( Ē) is the density of states.Furthermore, the r mn = r * nm are conventionally assumed to be (pseudo-)random Gaussian variables with zero mean and unit variance.While the ETH is an assumption and a formal proof is absent, its validity (including the Gaussian distribution of the r nm ) has been numerically confirmed for a variety of models and observables [8][9][10][11][12][13][14][15][16].Generally, the ETH is believed to hold for nonintegrable models and physical (for instance, spatially local) observables.In contrast, the ETH is violated in integrable systems due to their extensive number of integrals of motion [17], as well as in strongly disordered models which undergo a transition to a many-body localized phase in one dimension [18].In these cases, the off-diagonal matrix elements r nm deviate from the Gaussian distribution [16,19].In addition, models exhibiting a weaker violation of the ETH, such as, e.g., models featuring socalled "quantum scars", where rare less entangled states are embedded in an otherwise thermal spectrum, have recently attracted a significant amount of interest (see, e.g., [20,21]).
While the formulation of the ETH in Eq. ( 1) is conventional [3], it is to some degree incomplete with regard to the statistical properties of the O mn .Specifically, for a given H and O, the matrix elements O mn are predetermined, and therefore the notion of (pseudo-)randomness of the r nm needs to be carefully defined.In the spirit of the Bohigas-Giannoni-Schmit conjecture [22], we here advocate the strongest point of view, that beyond a certain energy scale all statistical properties of the off-diagonal matrix elements would match those of a Gaussian random ensemble.The central question of this paper is therefore to what extent the matrix elements O mn can be represented as independently drawn random numbers?
Clearly, all O mn cannot be random in the strict sense as they are constrained by the fact that the observables have to obey various algebraic relations (e.g.O 2 = 1 in case of O being a Pauli matrix acting on an individual spin).Furthermore, correlations between the r mn are necessary to reproduce the growth of certain fourpoint correlation functions in chaotic systems [23][24][25].Likewise, the consistency of relaxation dynamics in local systems also requires the r mn to be correlated [26].We therefore arrive at the important conclusion that the onset of random-matrix behavior has to be limited to matrix elements O mn within a certain energy window specified by the relevant energy scale ∆E RMT .
In this paper, we test the ETH in the case of a local spin operator in the eigenbasis of the paradigmatic spin-1/2 XXZ chain, for which we break integrability by means ´ µ H XXZ 0 5 ´ µ H 1 of (i) an additional next-nearest neighbor interaction or (ii) a single-site magnetic field in the center of the chain.(See Refs.[27][28][29][30] for related studies of the ETH and the emergence of quantum chaos in these models.)Going beyond the "standard" indicators of the ETH, we particularly investigate the existence of the scale ∆E RMT below which random matrix theory (RMT) prevails.To this end, we establish the eigenvalue spectrum of O as a sensitive probe of the correlations between the O mn .While the spectrum of the full spin operator includes only two eigenvalues ±1/2, we particularly focus on the spectrum of band submatrices at a fixed energy density Ē where the O mn are restricted to a narrow band |E n −E m | ≤ ω c .For such band submatrices, we demonstrate that the O mn are in convincing agreement with conventional indicators of the ETH in the following sense: (i) the diagonal matrix elements form a "smooth" function of energy O( Ē), (ii) the off-diagonal matrix elements follow a Gaussian distribution with a variance f 2 ( Ē, ω) that depends smoothly on the mean energy and respective energy difference, and (iii) the ratio between the variances of diagonal and off-diagonal elements for small ω takes on the value predicted by RMT.However, despite (i) -(iii) being satisfied, we find that the eigenvalues of the band submatrices for ω c larger than a certain value ∆E RMT still exhibit clear signatures of the original operator, implying correlations between matrix elements.At the same time, when the bandwidth is sufficiently decreased, the spectrum takes on an approximately semicircular shape, marking the transition where genuine random-matrix behavior occurs.
This paper is structured as follows.In Sec.II, we introduce the models and observables and describe our approach to study the ETH.Our results for "standard" indicators of the ETH are presented in Sec.III A while the existence of the scale ∆E RMT is investigated in Sec.III B. A summary and discussion is given in Sec.IV.

A. Models and observable
In order to test the ETH ansatz (1), we consider different (integrable and nonintegrable) quantum spin chains.A convenient starting point is provided by the one-dimensional XXZ model with open boundary conditions, where L denotes the number of lattice sites, S x,y,z ℓ are spin-1/2 operators at site ℓ, and ∆ is an anisotropy in the z direction.(In the following, we set the anisotropy to ∆ = 1.5.)While H XXZ is integrable in terms of the Bethe ansatz, we break integrability by either an additional next-nearest neighbor interaction of strength ∆ ′ [15,31,32], or by means of a single-site magnetic field h L/2 in the center of the chain [27-30, 33, 34], Note that, although not written explicitly in Eqs. ( 2)-( 4), we furthermore always include a small magnetic field at the first lattice site, h 1 S z 1 with h 1 = 0.1, which breaks the spin-flip and reflection symmetry of the model.
While H XXZ and H 1,2 conserve the total magnetization S z = ℓ S z ℓ , all results presented in this paper are obtained for the largest symmetry subspace which corresponds to S z = 0 and has dimension For L = 18, which is the largest system size we can treat numerically, we have D = 48620.Moreover, our simulations are performed for a representative choice of the integrability-breaking parameters, i.e., ∆ ′ = 1.2 and h L/2 = 1, for which both H 1 and H 2 are robustly nonintegrable (see also Refs.[15,[31][32][33][34] for other parameter choices).
The transition from the integrable XXZ chain to the nonintegrable models H 1 and H 2 can for instance be seen from the level-spacing distribution P (s) which is shown in Fig. 1.While the level spacings follows the Poisson distribution in the integrable case, P (s) matches the Wigner-Dyson distribution for H 1,2 .In this context, let us note that the field h 1 at the edge of the chain does not break integrability of H XXZ [33], while the single impurity h L/2 in the center of the chain induces the onset of chaos [33,34].10) between the variances of diagonal and off-diagonal matrix elements is obtained for regions of size µ shifted along the diagonal.Note that the matrix shown here comprises actual data for the example of H1 and L = 12.(c) We introduce a cutoff frequency ωc, where off-diagonal matrix elements are set to zero, Omn = 0, if |Em − En| > ωc, resulting in a band matrix with relative bandwidth W/D ′ .We study how the distribution of eigenvalues λ ωc 1 , . . ., λ ωc D ′ of the submatrix evolves upon reducing ωc.
For nonintegrable models such as H 1,2 , it is a general expectation that the matrix elements of physical observables O will follow the eigenstate thermalization hypothesis [3].In this paper, we test the ETH for the case of a local spin-1/2 operator acting on the central lattice site of the chain, Specifically, we employ full exact diagonalization to obtain the matrix elements O mn .Note that the indices m and n always refer to the eigenbasis of H.The O mn are real numbers for the chosen operator and the symmetry subspace.

B. Testing the ETH and the onset of RMT
In the following, we introduce the quantities studied in this paper.An accompanying sketch is provided in Fig. 2.

Indicators of diagonal ETH
The ETH ansatz (1) consists of the diagonal part and the off-diagonal part.The diagonal part of the ETH asserts that the function O( Ē) becomes "smooth" in the thermodynamic limit L → ∞.In particular, the eigenstate-to-eigenstate fluctuations O mm − O m+1m+1 should rapidly decrease with the system size L. One way to test this statement is to study the variance σ 2 d ( Ē) of the diagonal matrix elements, where the sum runs over all N Ē eigenstates |m with eigenenergies E m ∈ [ Ē − ∆E/2, Ē + ∆E/2] in a microcanonical energy window around a fixed Ē.For nonintegrable systems including our cases, it has been found that σ 2 d ( Ē) decreases exponentially with increasing L, while the scaling for integrable models is polynomial see, e.g., [3,9,15,35].

Indicators of off-diagonal ETH
Next, in order to test the off-diagonal part of the ETH, we consider matrix elements O mn in a (sufficiently narrow) energy window around a fixed Ē, where Ω( Ē) is approximately constant, which facilitates the analysis of the ω dependence of f ( Ē, ω) and of the distribution of the O mn , see Fig. 2 (a).A useful quantity in this context is the average over matrix elements in a small ω interval, which we denote in this paper by an overline.For instance, the average over |O mn | 2 in an interval of width ∆ω ≪ ω (with fixed Ē) is given by where the sum runs over all N ω matrix elements with versus ω yields the function f 2 ( Ē, ω), cf.Eq. ( 1), except for an overall prefactor [16,36,57].
Assuming the O mn have zero mean, i.e., O mn = 0, (which holds with a very high accuracy), we study the following quantity recently introduced in Ref. [16], which is sensitive to the distribution of r nm , When O mn = 0, the nominator in Eq. ( 9) coincides with the variance of the O mn while the denominator is the squared mean of the folded distribution.In particular, if O mn were to follow the Gaussian distribution, Γ(ω) = π/2.The value we find numerically in Sec.III is very close.To further confirm that the distribution P (O mn ) of the O mn is indeed Gaussian, we plot the histogram of O mn from narrow windows with fixed Ē and ω, see Fig. 5 below.
Next, we consider square-shaped submatrices of O of dimension D ′ < D around a fixed mean energy Ē.In Fig. 2 (b), an example for such a submatrix comprising actual numerical data is shown.As a further check that the O mn are normally distributed, we calculate the ratio Σ 2 (n, µ) between the variances of diagonal and off-diagonal matrix elements for eigenstates in (small) regions Here σ 2 d (n, µ) and σ 2 od (n, µ) are defined analogously to the variance in Eq. (7), see also [37] for details.
For an actual random matrix drawn from the Gaussian orthogonal ensemble (GOE), Σ 2 GOE = 2. Agreement with the GOE was anticipated in [3] and then verified numerically in, e.g., [13,14,38].Our results for Σ 2 (n, µ) in Sec.III are also in agreement with the GOE value.

Indicators of correlations between off-diagonal matrix elements
The indicators of ETH given in Eqs. ( 7) -( 10) have been studied before.In the present work we additionally consider the eigenvalue distribution of the submatrices with dimension D ′ < D around mean energy Ē [see Fig. 2 (b)], and show that it provides a much more sensitive probe of the statistical properties of the O mn .For a full random matrix with all matrix elements being independent, the eigenvalue distribution will follow the celebrated Wigner's semicircle [3,39].In contrast, if there are correlations between the O mn , deviations from the semicircle shape should emerge.Importantly, we find that the eigenvalue spectrum unambiguously shows that the r nm can not be represented as independent Gaussians variables, even though all standard indicators of the ETH are fulfilled, see Sec.III B below.
To identify the scale at which the transition to RMT behavior occurs, we analyze how the eigenvalue distribution depends on the width W of the band, see Fig. 2 (c).Specifically, let ω c denote some cutoff frequency.Then we define the new operator O ωc with matrix elements resulting in a band matrix with the relative bandwidth W/D ′ .Band random matrices have been extensively used in physics to model quantum systems and study their properties [40][41][42][43][44][45].Furthermore, the largest eigenvalues of full (square) and band submatrices have been studied in [26,38,46] in connection with the transition from integrability to chaos as well as relaxation dynamics and thermalization.However, to the best of our knowledge, the full eigenvalue distribution of (band) submatrices of local operators has not been previously considered as a quantity to characterize the presence of correlations between the matrix elements O mn .Provided all matrix elements are independent and identically distributed (except for an overall amplitude which may depend on ω), the eigenvalue distribution of band random matrices is expected to converge towards a semicircle for small W/D ′ [42], although there are corrections at intermediate W/D ′ and the detailed shape is more complicated [47], see Appendix C.
In addition to band submatrices, we also consider the eigenvalue distribution of full submatrices with varying dimension D ′ in Appendix A. One advantage of keeping D ′ fixed and varying W , however, is that the number of eigenvalues remains unchanged and is comparatively large.As shown in Appendix A, the properties of the smaller full submatrices are in fact similar and consistent with our findings for the band submatrices (11).
Given the ordered eigenvalues λ ωc α obtained by diagonalizing O ωc for the cutoff frequency ω c , an important quantity characterizing O ωc is the mean ratio r ωc of adjacent level spacings, where ∆ α = |λ ωc α+1 − λ ωc α | denotes the gap between two adjacent eigenvalues and the averaging is performed over a number (here N r ≈ D ′ /2) of gaps around the center.For a random matrix drawn from the GOE, one expects r GOE ≈ 0.53, while for uncorrelated Poisson distributed eigenvalues, one finds r Poisson ≈ 0.39 [48].In addition to r ωc , the central quantity in this paper is the full eigenvalue distribution P ωc (λ) of the band submatrix O ωc , where δ(•) denotes the delta function, and it is understood that individual peaks are collected in small bins such that P ωc (λ) forms a "continuous" distribution.
Given the corrections to the semicircle distribution for band random matrices with intermediate W/D ′ , a particularly simple and effective scheme to test the randomness of O ωc is to compare the eigenvalue distribution P ωc (λ) with the eigenvalue distribution of the suitably randomized O ωc .Specifically, O ωc is constructed by assigning random signs to the individual matrix elements O ωc mn (while keeping O ωc hermitian), If the matrix elements O ωc nm were random, we expect that the eigenvalue distribution would remain unchanged under this "sign randomization".In contrast, if the matrix elements of O ωc are correlated, these correlations will be erased by the randomization procedure and the eigenvalue distribution of the original and randomized matrices will be different.In order to quantify the difference (and its dependence on ω c ) between P ωc (λ) and the distribution P ωc (λ) of the randomized operator, we introduce −0.5 0 0.5 −0.5 0 0.5 ´ µ H 2 −0.5 0 0.5 ´ µ H 1 12 18 where P ωc (λ) and P ωc (λ) should be understood as the continuous distributions resulting from a binning procedure [see below Eq. ( 13)].If d 2 (ω c ) → 0, both distributions are very similar, which will be interpreted as a further indication that the matrix elements O ωc mn are randomly distributed.

III. RESULTS
Let us now turn to our numerical results for the matrix structure of S z L/2 in the eigenbasis of the two nonintegrable models H 1 and H 2 .The properties of diagonal and off-diagonal matrix elements are discussed in Secs.III A 1 and III A 2 respectively, while Sec.III B presents results for the eigenvalue distribution of band submatrices.Additional results for the integrable XXZ chain can be found in Appendix B. energy density E m /L for H 1,2 and different system sizes L = 14, 16, 18.For energy densities in the center of the spectrum, we find that the "cloud" of matrix elements becomes narrower with increasing L [9,15,27].This finding is in good accord with the ETH prediction that the O mm should form a "smooth" function of energy in the thermodynamic limit L → ∞.At the edges of the spectrum, the scaling with L is significantly slower.Especially for H 1 [Fig.3 (a)] and E m /L 0.2, the O mm are found to fluctuate very strongly.This can be understood as follows.The eigenstates of H 1 with the highest energies are weakly dressed domain-wall states.Consider, for instance, the states Given the distribution of the O mm , we restrict ourselves in the following to eigenstates in an energy window E m /L ∈ [−0.15, −0.05] which is close to the center of the spectrum (shaded area in Fig. 3).As shown in the inset of Fig. 3 (a), the variance σ 2 d ( Ē) of the O mm in this window decays approximately exponentially with L (at least for the system sizes numerically available), indicating that the diagonal part of the ETH is fulfilled.

Off-diagonal matrix elements
Let us now analyze the properties of the off-diagonal matrix elements O mn = m| S z L/2 |n .In view of the previous results for the diagonal matrix elements in Fig. 1, we focus on eigenstates with mean energy density Ē/L ∈ [−0.15, −0.05] as said above.For matrix elements in this window, running averages |O mn | 2 of their absolute squares are shown in Figs. 4 (a) and (b) as a function of ω both for H 1 and H 2 .The data are obtained for frequency intervals of width ∆ω = 10 −2 and system sizes L = 14, 16, 18.Overall, the situation is qualitatively similar for the two models H 1,2 .Namely, |O mn | 2 decays comparatively slowly at low frequencies, while a substantially quicker (presumably superexponential [49,50]) decay can be found at higher ω.In Figs. 4 (a) and (b) the values of |O mn | 2 for different L form smooth curve which collapse on each other when rescaled by the respective Hilbert-space dimension D. Except for a prefactor, these smooth curves correspond to the function f 2 ( Ē, ω) from the ETH ansatz (1).
For a more detailed analysis of |O mn | 2 , Figs. 4 (c) and (d) show a close-up of the low-frequency regime (note that the horizontal and the vertical axis have been rescaled).For both H 1 and H 2 , we observe that |O mn | 2 clearly approaches a nonzero value as ω → 0 with an approximately constant plateau for small ωL 2 (see also the discussion in Ref. [28]).´ µ H 1 0 5 ´ µ H 2 Next, in order to study the distribution of the O mn , Figs. 5 (a) and (b) show the frequency-dependent ratio Γ(ω), defined in Eq. (9).For small ω, we find that Γ(ω) is close to the Gaussian value π/2, while visible deviations appear at higher frequencies.However, these deviations decrease with the increasing system size L, indicating that the O mn follow a Gaussian distribution over a wide range of frequencies if L is sufficiently large.In addition, the full distribution P (O mn ) of the off-diagonal matrix elements is shown in Figs. 5 (c) and (d) for system size L = 18 and frequencies ω = 0.2, 0.4, . . ., 2. For all curves shown, we find that P (O mn ) is indeed well described by the Gaussians with zero mean (see also Refs.[10,16,19]).The width of the Gaussians is found to decrease with increasing ω, which is consistent with the data for |O mn | 2 shown in Figs. 4 (a) and (b).
Let us finally comment on the shaded gray area at low frequencies in Figs. 4 (a), (b) and Figs. 5 (a), (b).This area indicates the frequency range which is covered when considering a square-shaped submatrix in the interval Ē/L ∈ [−0.15, −0.05], cf.Fig. 2. Specifically, since this interval has a width ∆E/L = 0.1, the largest energy difference for L = 18 is ω max = 0.1L = 1.8.Therefore, when studying the eigenvalue distribution of such a submatrix further below, we are probing the region where |O mn | 2 varies comparatively slowly and Γ(ω) ≈ π/2.
To conclude the analysis of the off-diagonal matrix elements, Figs. 6 (a) and (b) show the ratio Σ 2 (n, µ) between the variances of diagonal and off-diagonal matrix ´ µ H 1 ω ´ µ H 2 ´ µ H 2 elements.Specifically, the data are obtained for L = 18 with two different square sizes µ = 100, 1000 and all possible embeddings along the diagonal of the submatrix with dimension D ′ .(Note that for the chosen energy window, we have D ′ ≈ D/4.)For small µ = 100, we find that Σ 2 (n, µ) fluctuates around the GOE value Σ 2 GOE = 2, both for H 1 and H 2 .Increasing the square size to µ = 1000, we observe that the fluctuations of Σ 2 (n, µ) are visibly reduced.For the larger value of µ, Σ 2 (n, µ) is still rather close to Σ 2 GOE for H 1 , while clear deviations between Σ 2 (n, µ) and Σ 2 GOE can be seen in the case of H 2 .
Figures 6 (c) and (d) show the averaged value, We find that Σ 2 (µ) ≈ Σ 2 GOE for small µ, while Σ 2 (µ) monotonously grows with increasing µ (this growth is particularly pronounced in the case of H 2 ).This behavior of Σ 2 (µ) follows from the ω dependence of |O mn | 2 discussed in Fig. 4. Since |O mn | 2 decreases with increasing ω, the variance σ 2 od (n, µ) likewise decreases with increasing µ, simply because matrix elements at higher frequencies are included.Comparing Σ 2 (µ) for different system sizes, we find that Σ 2 (µ) remains closer to Σ 2 GOE for a wider range of µ as L increases (see also [14]).Therefore ´ µ H 1 0 1 ´ µ H 2 ´ µ H 2 in the thermodynamic limit L → ∞ one can expect O nm to approach an actual random matrix drawn from the GOE, at least for a finite region around the diagonal.
To summarize, in this subsection we considered different quantities conventionally considered as standard indicators of ETH.The results presented in Figs. 3 -6 confirm that the matrix structure of the local spin-1/2 operator O = S z L/2 in the eigenbasis of the nonintegrable models H 1,2 is in good agreement with the ETH ansatz (1), at least for the chosen energy window close to the center of the spectrum.Nevertheless, in the next subsection, we will show that the matrix elements O mn within this energy window can not be considered as fully uncorrelated, i.e., the standard indicators in Figs. 3 -6 are not sufficient when it comes to the statistical properties of the O mn .

B. Beyond "standard" indicators: Eigenvalue distribution of band submatrices
We now turn to the eigenvalue distribution for the band submatrices centered around Ē/L = −0.1 with ∆E/L = 0.1.First we discuss the ratio of the adjacent level spacings r ωc defined in (12), which is shown in Fig. 7  ´ µ H 1 0.39 0.53 ´ µ H 2 versus ω c [panels (c) and (d)].We find that the behavior of r ωc is very similar for H 1 and H 2 .Specifically, over a wide range of ω c , r ωc ≈ 0.53 approximately matches the GOE value, while the transition towards the Poissonian value r ωc ≈ 0.39 occurs when the bandwidth becomes too narrow.The crossover from r GOE to r Poisson can be understood from the well-known fact that the eigenstates of a band random matrix with a sufficiently small value of W 2 /D ′ are localized and the eigenvalues become uncorrelated [41].In order to avoid localization effects while studying the eigenvalue distribution P ωc (λ), we restrict our analysis to ω c 0.03 (shaded area in Fig. 7), such that r ωc ≈ r GOE .
In Fig. 8, we show the full spectrum P ωc (λ) for H 1,2 with L = 18 and four exemplary choices of the cutoff frequency ω c .Specifically, we have chosen ω c = 1.8 (i.e. the full nonband submatrix), as well as ω c ≈ 1, ω c ≈ 0.4, and ω c ≈ 0.03, which are all above the transition point of r ωc .In all cases, we compare the spectrum of the bare operator O ωc to the distribution P ωc (λ) of the sign-randomized version O ωc , see Eq. ( 14).As can be clearly seen in Figs. 8 (a) and (b), P ωc (λ) and P ωc (λ) differ strongly for the largest ω c considered.Specifically, while P ωc (λ) closely follows the semicircle law appropriate for random matrices, P ωc (λ) still exhibits pronounced peaks at ±1/2, which is reminiscent to the original spectrum of the spin-1/2 operator.This deviation between 0 0.1 0.2 ´ µ H 1 ¸ωc = 1.8 ´ µ H 2 ¸ωc = 1.8 ´ µ H 1 ¸ωc ≈ 0.03 −0.5 0 0.5 ´ µ H 2 ¸ωc ≈ 0.03 P ωc (λ) and P ωc (λ) illustrates that the matrix elements O mn of a small submatrix cannot automatically be considered as independent random variables, notwithstanding all standard indicators of ETH being in agreement with the Gaussian distribution.This is a main result of the present paper.
Lowering the cutoff frequency to ω c ≈ 1, 0.4 and ω c ≈ 0.03 [see Figs. 8 (c)-(h)], we find that the spectra of the bare and the randomized submatrices become more and more similar.Especially for H 1 , P ωc (λ) and P ωc (λ) are very similar for ω c ≈ 0.4 and virtually indistinguishable from each other for ω c ≈ 0.03.Moreover, the bulk of the spectrum is convincingly described by a semicircular distribution (the width of the semicircle shrinks with ω c ), while small deviations from a perfect semicircle can be observed at the spectral edges.This similarity of P ωc (λ) ´ µ H 1 L = 14, 16, 18 and P ωc (λ) as well as their semicircular shape can be interpreted as an indication that the correlations between the matrix elements are significantly reduced for frequencies around and below ω ∆E RMT ≈ 0.4, i.e., on these smaller scales the O mn can be represented as independent random variables.This is another central result of the present work.Let us emphasize that the full distribution P ωc (λ) is sensitive to the RMT scale ∆E RMT while the mean gap ratio r ωc takes on the GOE value for all ω c considered in Fig. 8.
Comparing properties of H 1 and H 2 , we find that P ωc (λ) and P ωc (λ) still differ visibly for ω c ≈ 0.4 in the case of H 2 , see Fig. 8 (f), but become very similar for the smaller ω c ≈ 0.03, see Fig. 8 (h).This is in accord with Fig. 9 below, which suggests that ∆E RMT ≈ 0.1 is smaller in the case of H 2 .
While it certainly would be desirable to study the eigenvalue distribution P ωc (λ) for even smaller values of ω c in a controlled manner, this is difficult to do numerically as the value of W and the relative bandwidth size W/D ′ become too small.Likewise, if one instead diagonalizes full nonband submatrices with smaller dimension D ′ , see Appendix A, the number of eigenvalues becomes significantly reduced, which complicates the analysis.
Finally, Figs. 9 (c) and (d) show the difference d 2 (ω c ) between the two distributions P ωc (λ) and P ωc (λ).Consistent with our previous observation in Fig. 8, we find that d 2 (ω c ) decreases upon reducing ω c .Moreover, the decrease is slower in the case of H 2 .The minimum of d 2 (ω c ) is reached around the values of ω c which we associate with ∆E RMT (∆E RMT ≈ 0.4 in case of H 1 and ∆E RMT ≈ 0.1 in case of H 2 ), and d 2 (ω c ) remains low for smaller ω c .

IV. DISCUSSION
In this paper we have studied matrix elements of the spin-1/2 operator O = S z L/2 in the eigenbasis of two nonintegrable quantum spin chains.Specifically, we have considered the one-dimensional XXZ model in the pres-ence of two different integrability-breaking perturbations: an additional next-nearest neighbor interaction and a single-site magnetic field in the center.For these models and an energy window close to the center of the spectrum, we have shown that the matrix elements of O are in good agreement with the eigenstate thermalization hypothesis ansatz in the sense that (i) variance of the diagonal matrix elements decreases exponentially with the increasing system size, (ii) the off-diagonal matrix elements follow a Gaussian distribution with a variance that depends smoothly on the energy difference ω, and (iii) the ratio between the variances of diagonal and off-diagonal matrix elements approximately takes on the value predicted by random matrix theory.Overall, our results are in full agreement with previous works [27][28][29][30] and the conventional expectation that for local operators and nonintegrable Hamiltonians the ETH is satisfied.
The central question of this paper was to study to what extent off-diagonal matrix elements O mn can be treated as independently drawn random variables.To this end, we have considered submatrices around a fixed mean energy Ē and restricted the O mn to lie inside a sufficiently narrow band |E n − E m | ≤ ω c .We have established the form of the full eigenvalue distribution to be a sensitive probe to correlations between matrix elements.By comparing the eigenvalue distribution of the band submatrix with its sign-randomized counterpart (14), we have shown that the O mn cannot be considered as independently distributed, even on scales where O mn follow a Gaussian distribution with a variance that varies comparatively slowly with ω, i.e., on the scales where the ETH function f ( Ē, ω) is approximately constant.Specifically, while the spectrum of the sign-randomized matrix closely followed the semicircle law, matching the theoretical expectation for a random matrix, the eigenvalue distribution of the original submatrix was found to exhibit signatures of the spin operator, implying correlations between the matrix elements.When the cutoff frequency ω c is sufficiently reduced, we have found that the eigenvalue distribution of the original and the sign-randomized operator become similar and well described by a semicircle.The energy scale ∆E RMT when this transition occurs marks the onset of validity of the random-matrix behavior.
It should be noted that many important results rooted in the ETH are not sensitive to the statistics of the offdiagonal matrix elements O mn .This includes the central argument that the ETH ensures thermalization [3,6,51], which essentially relies on the exponential smallness of the O mn .At the same time, within the contemporary understanding of ETH, it is often assumed that matrix elements posses additional statistical properties matching the GOE (or some other appropriate Gaussian ensemble), see e.g.[3,10,13,19].In this work, we have advocated that beyond a certain energy scale all statistical properties of the off-diagonal matrix elements would match those of a Gaussian random matrix.We have provided numerical evidence that this onset of random-matrix behavior takes place below a certain energy scale ∆E RMT for specific models and observables.Our results suggest that for frequencies ω < ∆E RMT , the notion of (pseudo-)randomness of the r mn entering the ETH can be interpreted in an even stricter sense.At the same time, we have clearly seen that the scale ∆E RMT , where this transition to genuine random-matrix behavior occurs, is distinctly smaller than the scales on which "standard" indicators of the ETH are fulfilled.
Our work raises a number of straightforward questions.First, we note that our numerical observation ∆E RMT ≪ E τ mirrors the analytical bound ∆E RMT E τ /L established in [26], where E τ is defined as the width of the plateau of f ( Ē, ω) (note that E τ is sometimes referred to as Thouless energy [36]).A natural question would be to establish the scaling of ∆E RMT with the system size L and, in particular, to investigate if (∆E RMT ) −1 can be associated with the time scale of late time chaos at which the dynamics of various observables is captured by RMT [52][53][54][55].
Another direction is to contrast the behavior in chaotic systems with the integrable counterparts.We repeat the analysis of section III for the integrable XXZ model in the Appendix B. One particular observation to point out is that off-diagonal matrix elements in the integrable case also can be regarded as random and independent, although not Gaussian, beyond a certain energy scale.We leave for the future the question of better understanding this transition, and the universal properties of O mn in the integrable case.
Eventually, one avenue of research is to characterize the nature of the correlations between off-diagonal matrix elements for ω > ∆E RMT , and to understand their potential impact on self-averaging properties of the O mn exploited in various works [45,[56][57][58].At the same time, it would be interesting to study the connection between universal properties of O mn at small frequencies and transport, which could be diffusive or ballistic [59].row.In Fig. 10 we show a qualitatively similar result for the Hamiltonian H 1 and full (nonband) submatrices, where the width of the energy window is now chosen as ∆E/L ≈ 0.06 and ∆E/L ≈ 0.012, i.e., narrower than in the main text.Recall that a smaller ∆E a smaller submatrix dimension D ′ .For the examples shown here, we have D ′ ≈ 8000 and D ′ ≈ 1500.While P ωc (λ) still exhibits pronounced features at λ = ±1/2 for the larger ∆E in Fig. 10 (a), we find that the distributions P ωc (λ) and P ωc (λ) are essentially indistinguishable for the smaller ∆E in Fig. 10 (b).Moreover, analogous to the results for the band submatrices in Fig. 8, small deviations from a perfect semicircle law appear at the spectral edges if ∆E is lowered.While we cannot entirely exclude the possibility of finite-size effect, we conclude that our results for band submatrices in the main text, i.e., a semicircular bulk with small deviations at the spectral edges, are not just caused by the finite bandwidth, but are stable features which appear for full nonband matrices as well.
Appendix B: Results for the integrable model HXXZ In Figs. 3 -6 of the main text, we have analyzed the ETH structure of S z L/2 written in the eigenstates of the two nonintegrable models H 1 and H 2 .In Fig. 11, we present analogous data for the integrable model H XXZ .As expected, the results for H XXZ are drastically different compared with the nonintegrable systems H 1,2 .In particular: (i) the width of the distribution of the diagonal matrix elements does not visibly shrink with increasing L, (ii) Γ(ω) = π/2 and is nonconstant and dependent on L, in agreement with [16], (iii) the ratio Σ 2 (n, µ) is orders of magnitude larger compared to Σ 2 GOE = 2, and (iv) the distribution P (O mn ) is clearly non-Gaussian (see also [10,16]).Overall, these results confirm the expectation that the ETH is not satisfied in the case of integrable models.The only quantity which exhibits a similar behavior for H XXZ and H 1,2 is the running aver- age |O mn | 2 shown in Fig. 11 (b).Namely, we find that plotting |O mn | 2 versus ω for system sizes L = 14, 16, 18 yields smooth curves which collapse onto each other when rescaled by the respective Hilbert-space dimension D. This behavior is in agreement with the recent studies in Refs.[16,60].Finally, we consider eigenvalue distribution for the band submatrices in the case of integrable model H XXZ .Figure 12 (a) shows results for the level-spacing ratio r ωc , while Fig. 12 (b) shows the eigenvalue distributions P ωc (λ) and P ωc (λ) of the original and the randomized band submatrices with the cutoff frequency ω c ≈ 0.1.Comparing with the results for the nonintegrable models H 1,2 , the qualitative behavior of both r ωc and P ωc (λ) appears to be very similar.Namely, we find that r ωc exhibits a crossover from r GOE to r Poisson when the cutoff frequency ω c (or bandwidth W ) decreases.Furthermore, P ωc (λ) and P ωc (λ) have the very similar shape for the considered value of ω c , which indicates that O mn within the corresponding band can be considered as independent.Accordingly, as discussed in the Appendix C, eigenvalue distribution is approximately semicircular.Given these results, we conclude that mutual independence of the off-diagonal matrix elements beyond certain energy scale ω < ∆E RMT is also present in the integrable models, raising the question if an appropriate non-Gaussian random matrix theory can capture the universal properties of the O mn in this case.

FIG. 2 .
FIG. 2. (Color online) (a) The ETH ansatz (1) is studied for the spin-1/2 operator S z L/2 written in the eigenbasis of the respective Hamiltonian H.The red shaded area indicates matrix elements around a fixed value of Ē where the density of states is approximately constant, while the gray shaded area indicated a square-shaped submatrix in this energy window.(b) For submatrices with dimension D ′ < D, the ratio Σ 2 (n, µ) defined in Eq. (10) between the variances of diagonal and off-diagonal matrix elements is obtained for regions of size µ shifted along the diagonal.Note that the matrix shown here comprises actual data for the example of H1 and L = 12.(c) We introduce a cutoff frequency ωc, where off-diagonal matrix elements are set to zero, Omn = 0, if |Em − En| > ωc, resulting in a band matrix with relative bandwidth W/D ′ .We study how the distribution of eigenvalues λ ωc 1 , . . ., λ ωc D ′ of the submatrix evolves upon reducing ωc.

∝
FIG. 3. (Color online) Diagonal matrix elements of S z L/2 in the eigenbasis of (a) H1 and (b) H2, for system sizes L = 14, 16, 18.The shaded area indicates the energy window Em/L ∈ [−0.15, −0.05] which is used in the following to further study the properties of off-diagonal matrix elements.For the example of H1, the inset in (a) shows that the variance σ 2 d ( Ē) of the Omm in this window decreases exponentially with increasing L.
A. "Standard" indicators of the ETH1.Diagonal matrix elementsAs a first step, we study the diagonal part of the ETH.To this end, Figs.3 (a) and (b) show the matrix elements O mm = m| S z L/2 |m as a function of the corresponding

2 FIG. 4 .
FIG. 4. (Color online) [(a),(b)] Running averages |Omn| 2 of matrix elements in the energy window Ē/L ∈ [−0.15, −0.05], calculated for system sizes L = 14, 16, 18 and frequency bins of width ∆ω = 0.01.[(c),(d)] Close-up of the low-frequency regime using a bin width of ∆ω = 5 × 10 −4 .Note that both the horizonal and the vertical axis have been rescaled.Panels (a) and (c) show results for H1, while (b) and (d) show data for H2.The shaded area in panels (a) and (b) indicates the ω range which is probed when considering a squareshaped submatrix with eigenstates in an energy interval of width ∆E/L = 0.1.

FIG. 7 .
FIG. 7. (Color online) Mean ratio rω c of adjacent level spacing of the operator O ωc versus [(a),(b)] the scaling parameter W 2 /D ′ , and [(c),(d)] the cutoff frequency ωc.Panels (a) and (c) show data for H1 while (b) and (d) show data for H2.The data is obtained for system sizes L = 16, 18 and the dashed horizontal lines indicate the GOE value rGOE ≈ 0.53 and the Poisson value rPoisson ≈ 0.39.For our analysis of Pω c (λ), we restrict ourselves to ωc 0.03 [shaded area in (c) and (d)], such that rω c ≈ rGOE.

FIG. 9 .
FIG. 9. (Color online) d2(ωc) for (a) H1, and (b) H2.System sizes are chosen as L = 14, 16, 18 and the arrow indicates the direction of increasing L. The dashed vertical lines indicate the approximate location of ∆ERMT below which randommatrix behavior occurs.

FIG. 10 .
FIG. 10. (Color online) Eigenvalue distributions Pω c (λ) and Pω c (λ) for full (nonband) submatrices in the case of the model H1 with L = 18.The matrices are centered around the mean energy Ē/L = −0.1, but the width of the energy window is now smaller compared to the main text, i.e., we choose (a) ∆E/L = 0.06 and (b) ∆E/L = 0.012, which would correspond cutoff frequencies ωc ≈ 1 and ωc ≈ 0.2.
|n 1 and |n 2 are not exact eigenstates).While |n 1 and |n 2 have almost the same energy, one finds that n 1 | S z L/2 |n 1 ≈ 1/2 whereas n 2 | S z L/2 |n 2 ≈ −1/2 such that ETH is not satisfied.We expect the range of energy densities where the ETH applies to increase with L.