Probing Off-diagonal Eigenstate Thermalization with Tensor Networks

Energy filter methods in combination with quantum simulation can efficiently access the properties of quantum many-body systems at finite energy densities [Lu et al. PRX Quantum 2, 020321 (2021)]. Classically simulating this algorithm with tensor networks can be used to investigate the microcanonical properties of large spin chains, as recently shown in [Yang et al. Phys. Rev. B 106, 024307 (2022)]. Here we extend this strategy to explore the properties of off-diagonal matrix elements of observables in the energy eigenbasis, fundamentally connected to the thermalization behavior and the eigenstate thermalization hypothesis. We test the method on integrable and non-integrable spin chains of up to 60 sites, much larger than accessible with exact diagonalization. Our results allow us to explore the scaling of the off-diagonal functions with the size and energy difference, and to establish quantitative differences between integrable and non-integrable cases.


I. INTRODUCTION
Since the early days of quantum mechanics, the emergence of thermalization behavior in isolated quantum systems has been a fundamental and intriguing question [1,2].But it has been only in recent years that, thanks to the high levels of control and isolation of ultracold atomic experiments, it has become possible to explore the quantum thermalization phenomenon experimentally [3][4][5], which has rekindled the attention to the topic.
Tensor network (TN) methods offer possibilities to numerically explore quantum systems of sizes much larger than the ones allowed by ED.The most successful TN methods target equilibrium states (ground or low energy eigenstates, or thermal equilibrium [38][39][40]).But probing ETH requires investigating states at finite energy, which typically are not efficiently described by a TN ansatz [41].Nevertheless, a new algorithm has been recently proposed that precisely allows studying an ensemble of eigenstates at finite energy density [42,43].The method simulates the effect of a narrow energy filter operator through its expansion as a sum of evolution operators.By (quantum or classically) simulating each of these evolutions and post-processing the data, it is possible to approximate the properties of a microcanonical ensemble over an extensive region of the spectrum.Classically simulating the evo-lution with tensor networks imposes a limit on the width of the accessible filters, but, as demonstrated in [43] it suffices to efficiently access the microcanonical values for spin chains up to 80 sites, thus providing a way to probe the diagonal part of the ETH ansatz.
In this paper, we generalize the applications of the energy filter method to probe the more challenging off-diagonal part of ETH.More concretely, our method computes a (broadened) filter spectral function for a given (local) operator.This function can be interpreted as an average of matrix elements over eigenstate pairs selected by two spectral filters, one selecting the average energy and the other the energy difference.We demonstrate how, by simulating finite time evolutions with standard tensor network routines, in the spirit of [43], we can obtain the effect of two combined filters.Using our method we compute and compare the spectral functions for several Ising spin chains, including integrable and non-integrable clean systems, as well as a disordered one, for much larger system sizes than allowed by exact diagonalization.For the region we can reliably probe, we obtain convergence of the spectral functions with system size, and are able to discriminate qualitatively distinct features, such as a different scaling with energy difference.
The rest of the paper is structured as follows.Section II provides a brief review of ETH and the filter ensemble.In section III we present the method used to compute the spectral function numerically with TNS algorithms.Section IV describes the three different spin models studied, and collects our numerical results.In particular our results capture the asymptotic behavior of off-diagonal matrix elements, and exhibit qualitative differences in the energy and energy difference dependence between the integrable and generic cases.We also study how in the latter case, the fluctuation dissipation relation is fulfilled by the filter ensemble for large enough systems.

A. Eigenstate Thermalization Hypothesis (ETH)
Consider a many-body Hamiltonian with spectral decomposition Ĥ = α E α |α⟩ ⟨α|.Given a physical observable Ô, ETH predicts that its matrix elements O αβ ≡ ⟨α| Ô |β⟩ in the energy eigenbasis obey the following form [7,9] ( where we define the energy variables Ē ≡ (E α + E β )/2 and ω ≡ E β − E α .R αβ is a random variable with zero mean and unit variance.The thermodynamic entropy S(E) can be defined as the logarithm of the number of available states in the microcanonical window, DoS(E)∆E.In the literature it is nevertheless common to drop the dependence on the window width, which contributes only a small constant [44], and use as definition S(E) = ln DoS(E).O( Ē) and f O ( Ē, ω) are smooth functions of their arguments.While Eq. 1 is probably the most common one for the ETH ansatz, other expressions exist that use a different entropy factor [10], which results in a slightly different definition of f O ( Ē, ω).In our case, we write the off-diagonal term as e −[S(Eα)+S ETH is a sufficient condition for quantum thermalization [9]: starting from a (sufficiently narrow in energy) out-ofequilibrium state, if (1) is satisfied, the time-averaged expectation value of physical observables will relax to its microcanonical average O(E) in the limit of infinitely long time, with the fluctuations around the thermalization value controlled by the off-diagonal matrix elements.Furthermore, it is believed that ETH holds for generic non-integrable systems and fewbody operators, and multiple numerical studies have been conducted to verify its validity in such scenarios [23][24][25][26][27][28][29].Violations of ETH can be observed in integrable systems [30][31][32][33][34][35][36], as well as strongly disordered models [20][21][22]37], due to their extensive number of (quasilocal) integrals of motions.
The off-diagonal structure function |f ( Ē, ω)| 2 appearing in the ETH ansatz is related to dynamic properties, and also determines the fluctuation-dissipation theorem (FDT) of nonequilibrium states.Recent works have focused on studying some of its properties in generic and non-generic systems.The statistics of off-diagonal matrix elements in integrable spin chains were analyzed in [19,25,[34][35][36] and its dependence on energy difference ω and system size N in [9,[24][25][26][27][33][34][35][36][37].On the other hand, the off-diagonal matrix elements in disordered models can display spectral properties deviating from the ETH prediction, as for instance shown in [20,45,46].Finally, the standard ETH does not make explicit predictions for the correlation between matrix elements, which is nevertheless related to quantum chaos, and has been a focus of recent works [47][48][49][50][51][52].

B. The filter ensemble
The most direct way to verify ETH numerically is to diagonalize the many-body Hamiltonian and analyze the exact energy eigenstates in a targeted energy window.This is, however, limited to small systems of the order of 20 spins, as the dimension of the Hilbert space increases exponentially with the size.
A potential workaround is to study, instead of the properties of individual eigenstates, those of an ensemble that is narrow in energy.This strategy has been recently followed in [42,43,53] to investigate finite energy properties using the filter ensemble where g σ (x) is the Gaussian function, The filter ensemble ρσ (E) is diagonal in the energy eigenbasis, and it is centered at E, with σ being the energy width.The expectation values for ρσ (E) will converge to the microcanonical ones in the limit of small σ, such that the filter ensemble can be used to probe the diagonal part of the ETH.In [43] it was argued that, in a system fulfilling ETH, a width σ = o( √ N ) should be sufficient for the l.h.s. of (4) to converge to the microcanonical values for intensive quantities (see also [54]).

III. THE SPECTRAL FUNCTION OF FILTER ENSEMBLES
In this work, we are interested in the application of filter methods to probe the off-diagonal structure of the ETH.In particular, we aim at extracting information about the offdiagonal structure function |f O (E, ω)| 2 .The main quantity in our study will be the spectral function of the filter ensemble, defined as follows.The autocorrelation for some local observable Ô in a given state ρ can be computed as Its Fourier transform yields the spectral function If the state is diagonal in the energy basis, we denote its matrix elements ⟨α| ρ |α⟩.This is the case for a single eigenstate or for the Gibbs ensemble, but also for the filter ensemble defined in Eq. ( 2).In such case, the spectral function can be written as i.e. S ρ O (ω) is an average over squared matrix elements O αβ = ⟨α| Ô |β⟩ between energy eigenstates with fixed energy difference ω, weighted by the probability of the eigenstates |α⟩ in the distribution defined by ρ.In the following we focus on the filter ensemble ρσ (E), for which the corresponding probabilities are ⟨α| ρ |α⟩ ∝ g σ (E − E α ), for the function g σ defined in Eq. ( 3), with the normalization factor specified in Eq. (2).
In order to compute this quantity, we introduce a generalized version of the spectral function, where we replace the δ function in Eq. ( 6) by a Gaussian of width σ ω , which we will approximate by a second filter acting on the energy difference E β − E α .This results in a broadened spectral function S ′ρσ(E) O (t) that performs an average of squared matrix elements |O αβ | 2 for states α in the support of ρ σ (E) and states β with energies around E α + ω, as graphically illustrated in Fig. 1.The function can be written as For a local and bounded Hamiltonian, the density of states converges weakly to a Gaussian distribution in the thermodynamic limit, with a width of √ N σ 0 , where σ 0 is a constant independent of the system size [55,56].If the filters are narrow enough compared to

√
N σ 0 , we can consider the density of states almost constant within the peak (while still much wider than the level spacing), and we can express S ′ρσ(E) O in terms of an average of matrix elements as where the average is taken over pairs of eigenstates with energies around E and E + ω, respectively.

Relation with ETH
Since the filter does not select precise energy differences, the sum will in general also pick up a contribution from the (much larger) diagonal matrix elements.In order to study the off-diagonal part of ETH, we thus choose observables for which the microcanonical value, and thus the diagonal contribution, vanishes, i.e.O(E) = 0.For these observables, the ETH ansatz predicts and thus Therefore, S ′ρσ(E) O (ω) is directly related to the function f O and should also be a smooth function if ETH holds.The finite filter widths imply multiplicative corrections to Eq. ( 10) of order O(σ 2 /N 2 + σ 2 ω ), as explicitly shown in appendix B. Since the filter strategy can be used to determine the density of states (see also [53,57]), we can extract the value of |f O | from the computed spectral function.But for Hermitian Ô, a better strategy is making use of the fact that f O (E, ω) = f O (E, −ω), which allows us to eliminate the exponential factor by combining the spectral functions at different arguments in a single function where the last part holds for ETH, as it follows from (10).The function V O (E, ω) can be understood as the variance of the matrix elements within the filter, and can obviously be computed for any system, fulfilling ETH or not, but in the latter case, it is not ensured to be a smooth function.
Notice that our method is also applicable to operators with nonzero microcanonical values.In such case, the microcanonical value can be approximated using the method in [43,53], and subsequently subtracted from our calculations, to extract the off-diagonal component.
Finally, it is worth noticing that the regularized correlators introduced in [58], which can be related to filtered functions (see appendix A 3), provide a similar strategy to extract the function |f O |.

The filter method
Our numerical strategy is based on the TN simulation of the filter operators presented in [43].While the details are discussed in [42,59,60], we sketch here the main steps for the sake of clarity.It is convenient to approximate the Gaussian by a cosine function as where M = ⌊(α/σ) 2 ⌋ 2 ( ⌊. ..⌋ 2 indicating the nearest even integer) and α is a rescaling factor introduced to ensure that the range of the argument ξ/α is smaller than the period of the cosine function π.
Using the binomial expansion, the cosine power can be written as a sum of M + 1 complex exponentials.The number of terms in this sum can be reduced to O(x √ M ), by introducing a small error controlled by x = O(1), yielding where t m = 2m/α and c Using the expansion (13) for both Gaussian functions in Eq. ( 7) we obtain the following expression for the generalized spectral function S where in the second line we have simply used the trace and product of operators to rewrite the sums over the spectrum.Notice that both Gaussian filters have independent widths, corresponding to two different expansion parameters M and M ω .Also the normalization factor α could in principle be different for each filter.Nevertheless, it is convenient for the numerics to use a common value of α, so we choose the largest of both.Notice that a different pair of filters could be used resulting in an average of matrix elements that probes the structure of the off-diagonal matrix elements using these filters.We briefly introduce other possibilities in App. A.

Cosine filter parameters
The cosine filter for the ensemble is determined by the three parameters (σ, α, x).Similarly, the triplet (σ ω , α ω , x ω ) determines the cosine filter approximating the second Gaussian function in Eq. (7).The flexibility in parameter selection provides greater control over the performance of the approximation.
As mentioned above, for simplicity, we choose the same rescaling factor for both of the filters.In practice we find α = α ω > E max −E min , with E max (E min ) being the highest (lowest) energy in the spectrum of H, to be enough to ensure the proper bound of both cosine filter arguments in the regime we study.In order to determine the suitable value, we estimate E max (E min ) using a variational MPS optimization, and choose α > 1.1(E max − E min )/π.This fixes the time step in the filter expansion ∆t = 2/α to be the same for both filters, such that we can use common values for t m and t n in Eq. ( 14).
The longest time in the filter expansion scales as The smallest accessible widths are fixed by the longest times that we can reliably simulate using the TNS algorithms given the available computational resources and the finite precision of the numerical estimates.In App.A we discuss in detail how we choose t max and x for the specific models under study, to achieve the smallest filter widths while keeping the numerical error under control.Based on our error analysis results, we employ filter widths of σ = O( √ N ) and σ ω = O(1).

Tensor network simulations
Each of the terms in Eq. ( 14) can be evaluated numerically with TNS techniques similar to the ones employed in [43].In particular, here we need to compute the trace expressions Tr[e iHtm ] and Tr[e iHtm Ô(t n ) Ô † ].To calculate the latter, we approximate e iH(tm+tn)/2 Ôe −iHtn/2 using matrix product operators (MPO) [38][39][40] In our simulation this is achieved by the time-evolving block decimation (TEBD) algorithm [38][39][40].Evaluating the trace Tr[e iHtm Ô(t n ) Ô † ] corresponds to computing the inner product of the two MPOs, which can be done efficiently.This strategy of evolving the MPO on both physical indices, splitting the evolution between both operators and evaluating an inner product has been shown to extend the time one can reach with a limited bond dimension [61,62].
Similarly, in order to obtain Tr[e iHtm ], we compute an MPO approximation to the evolution operator at times t m .Evaluating the trace is then equivalent to a contraction with the vectorized identity operator.
We have studied several spin models for different system sizes, up to N = 60.In all cases, we appoximated the evolution operators by a second order Trotter expansion with small Trotter step 0.01.The results shown in the following were obtained with maximal bond dimension D = 600, which found to be enough to guarantee convergence for the studied timescales (see app.A for more detailed description of the numerical errors).

IV. PROBING ETH WITH SPECTRAL FUNCTIONS:
NUMERICAL RESULTS

A. Setup
We benchmark the method on a quantum Ising chain with open boundary conditions, selectively including a disordered field and an integrability-breaking next-to-nearest-neighbor term, J sets the energy scale and in the following, we fix it to J=1.
The model is integrable when J 2 = 0, when it can be mapped to free fermions.The transverse field includes a homogeneous component h, and potentially a disordered one r i .For simulations of the disordered model, the values of r i are sampled from the uniform distribution in the interval [−r, r].
We focus on three different sets of parameters.The first one, (J 2 , g, r) = (0.0, 1.05, 0) is the transverse field Ising model with uniform potential and thus integrable.The second one is (J 2 , g, r) = (0.2, 1.05, 0), which corresponds to a nonintegrable case, for which ETH is expected to hold.We finally consider also a disordered case, (J 2 , g, r) = (0.2, 0, 3.0), where the on-site field takes random values.
The observable we focus on is the longitudinal magnetization of the central site Ô = σ z N/2 .Notice that Ĥ is invariant under the transformation F = N j=1 σ x j , which flips all the spins in the chain.As σ z N/2 anti-commutes with F, the expectation value of σ z N/2 in eigenstates of both H and F is 0, that is the diagonal elements are automatically zero and only offdiagonal elements contribute.In the integrable system, where the eigenstates are characterized by free excitations, σ z N/2 is a many-particle operator in terms of these excitations, and thus the majority of off-diagonal elements are non-zero [19].

B. Effect of the filter widths
The spectral function obtained with our method depends on the filter parameters described above.Thus, first of all, we need to analyze the effect of the filter widths σ and σ ω on the results.To isolate the influence of the filter ensemble width σ, we study the dependence of the autocorrelator C ρσ O (t) (independent of σ ω ) on this parameter.Figure 2   with mean energy density E/N = 0.5, using filters of varying width σ = q √ N , with proportionality factor q = 0.5, 0.2 and 0.1.We observe that the results are converged for σ ≤ 0.2 √ N .We observe a similar convergence for all the size and energy ranges studied in this work, thus in the following we choose values of σ within that range.
Additionally, from Fig. 2 it becomes evident that the various studied models exhibit widely different time dependence of the autocorrelator.In the clean systems (Figs.2(a) and 2(b)), the autocorrelator decays exponentially with t (notice the logarithmic scale of the vertical axis), whereas for the disordered system (Fig. 2(c)) it remains significant at long times.
Next, to study the effect of the width of the energydifference filter, we fix σ = 0.2

√
N and vary σ ω .In the clean systems, the spectral function is a smooth function.In contrast, for the disordered case, the function exhibits multiple peaks, at large values of the energy difference.This is a feature observed in localized systems [20,45,46].The results shown in Fig. 3(c) correspond to a particular realization of the disorder, but the figure is qualitatively similar for other realizations, with the positions of the peaks varying.
It is also worth noticing that in the disordered case the error originated by the sum truncation in (13) is more significant, so we need to choose a larger x ω factor.We select x ω ≥ 3 which, together with the upper bound on the simulated time t max ≤ 20 imposed by the truncation error, means we can reach σ ω ≥ 0.3 (see App. A).
The previous arguments allow us to fix the filter widths σ and σ ω to suitable values in the following studies.Specifically, we set the ensemble width to σ = 0.2

√
N for all systems.Thus, for the clean systems, we choose σ ω = 0.1, and for the disordered model, σ ω = 0.3.

C. The off-diagonal matrix elements
To study the off-diagonal matrix elements, we analyze the function V O (E, ω) in Eq. (11).As discussed in section III, Just as the spectral functions, V O (E, ω) appears smooth in the clean systems [Fig.5(a,b)], whereas it exhibits multiple sharp peaks for the disordered system [Fig.5(c)].This feature is also visible for the non-interacting (J 2 = 0) version of the same disordered chain and in the interacting case becomes more pronounced for stronger disorder, which seems to relate it to localization, even though the case shown in Fig. 5(c) corresponds to moderate disorder strength.
Even though the function V O (E, ω) is smooth for both clean cases, significant differences are visible.In the integrable system, V O (E, ω) exhibits a symmetry V O (E, ω) = V O (−E, −ω).This property follows from the particle-hole symmetry in the integrable model, patent in its fermionic formulation, as explicitly shown in appendix C. In the nonintegrable system there is no such symmetry, and V O (E, ω) displays very different functional dependence on ω at high and low energies E, as seen in fig. 5

(b).
A remarkable advantage of our method is that it can be applied to much larger systems than exact diagonalization, thus making it possible to address the system size scaling of these features.To study the system size dependence of V O (E, ω), we plot the function in Fig. 4 at a fixed mean energy density E/N = 0.6 for system sizes up to N = 60, for the three models introduced above.In all cases, we obtain convergence of the function with the system size, showing that the method allows us to observe the asymptotic behavior.The figures show that for ω > 5, V O (E, ω) decreases very fast with ω.However, the form of the decay is qualitatively different.As can be appreciated in the logarithmic plots shown in the insets, for the non-integrable case [Fig.4(b)], the decay is exponential, as expected from ETH [9].In contrast, for the integrable case, the decay of V O (E, ω) is faster, compatible with exp(−k 1 ω 2 ), [see inset of Fig. 4(a)].The exp(−k 1 ω 2 ) decay of V O (E, ω) at high frequency regime was also found in [34][35][36] for other integrable models, such as XXZ chain and hard-core bosons.

D. Numerical test of the Fluctuation-Dissipation Theorem
The fluctuation-dissipation theorem is a general property of systems in thermodynamic equilibrium.For quantum systems in the Gibbs state ρβ = e −β Ĥ /Tr[e −β Ĥ ], the spectral function automatically satisfies the Kubo-Martin-Schwinger (KMS) condition [63,64], which can be expressed This is a sufficient and necessary condition for the fluctuationdissipation theorem (FDT) to hold.But the KMS condition can hold in more general situations.Beyond thermal equilibrium, it has been shown to hold in particular for some non-equilibrium initial states [65,66].More relevant for our case, it holds also for individual eigenstates [26], and for ensembles narrow in energy [9,67] in the thermodynamic limit when ETH is valid.We thus expect the filter ensemble to also fulfill FDT in the general case, at least in the limit of large systems.Using the filter strategy we can actually probe the validity of the relation in finite systems as a function of size.
We can define the following indicator function that will test the KMS condition (hence the FDT) for a generic ensemble [67], If the FDT holds we expect the function to be independent of ω, and equal to the inverse microcanonical temperature corresponding to the mean energy of the ensemble ρ.For a generic system (fulfilling ETH), we can compute the value of ( 18) in the filter ensemble by expanding the spectral function around E as where β ≡ ∂ E S(E) is the inverse temperature at energy E in the microcanonical ensemble.We obtain for the indicator function We thus expect the function to be approximately equal to the inverse temperature, as KMS requires, with some correction.The major correction term ∂ E ln |f O (E, ω)| 2 scales as O(1/N ) (see App. B for more details).Thus we expect that, in the generic case, FDT indeed holds for the filter ensemble in the thermodynamic limit.
We have computed this indicator function using our generalized spectral function for the non-integrable Ising chain, up to N = 60 sites.To be able to compare different system sizes, we fix a value of the (reference) inverse temperature β 0 .This corresponds to a value of the mean energy E 0 fulfilling β 0 = ∂ E S(E 0 ), which we can determine from the results of the last section.To be specific, we obtain DoS(E) with filters and then extract ∂ E S(E) = ∂ E ln DoS(E) using finite derivatives.Then, we compute the indicator function Eq. ( 18) from the ratio between the values of the spectral function at ω and −ω at the corresponding energy E 0 .To check the validity of our approach, we plot in Fig. 6 the value obtained for this indicator function as a function of 2 , which we obtained in an analogous way to ∂ E S(E 0 ).The figure shows this dependence for two values of the reference inverse temperature, β 0 = 0 and 0.2, in the range ω ∈ [ −5, 5].The plots show a near-perfect linear dependence, passing near the point (0, β 0 ), as predicted by Eq. (20).
Our findings are consistent with those of the ED studies in [26,67], but extend the results to significantly larger system sizes.

V. CONCLUSIONS
Using energy filter operators offers an alternative strategy to study the properties of the off-diagonal matrix elements of observables in the energy eigenbasis, and thus probe the ETH, fundamental ingredient in the theoretical understanding of quantum thermalization.In this work we have explored this possibility, with focus on the spectral function of the filter ensemble S ρσ(E) O (ω).We have shown that this quantity corresponds to an average of such matrix elements, and gives access to the off-diagonal function in the ETH ansatz, with corrections that depend on the filter width and the system size.
We have shown that this quantity can be simulated classically using TNS techniques, in a generalization of methods presented in [43] for the microcanonical averages.In particular, this strategy allows addressing much larger systems than exact diagonalization, and allows us to observe convergence in the system size and thus to identify the asymptotic features of the off-diagonal matrix elements.It thus provides a powerful tool to explore the ETH ansatz.
In order to test the strategy, we have applied it to several operators in Ising chains including different terms, such that they span from integrable, non-integrable generic (ETH) to ergodicity breaking behaviors.In particular, we have shown that the spectral functions S ρσ(E) O (ω) for a non-integrable, generic instance and for a disordered chain exhibit clear qualitative and quantitative differences.
We have also shown that the filter spectral function provides a way to probe the validity of the FDT for the filter ensemble.In the limit of vanishing filter width this would converge to a probe of the relation for energy eigenstates, so far realized with ED for small systems [67].Our numerical results show good agreement with FDT in the generic non-integrable chain up to 60 sites.The filter width σ For fixed x, the filter width σ inversely depends on the max(t m ).As observed from Figure 7, the cutoff time max(t m ) becomes shorter for larger systems.Here we provide a theoretical explanation of the dependence of max(t m ) on system sizes.For a traceless, local and bounded Hamiltonian, the density of states DoS(E) converges weakly to Gaussian distribution [55,56] in the thermodynamic limit with width proportional to where d is the local Hilbert space and σ 0 is a constant independent of the system size., and therefore the argument above also applies here.As for the time dependence of |Tr[ Ô(t n ) Ô † ]|, we find it varies for different models.We thus simulate values of t n as large as possible, until the quantity becomes too small, or the truncation error in our MPO approximation becomes significant.
In Fig. 8, we compare the |Tr[ Ô(t n ) Ô † ]| evaluated with bond dimension D = 400 and D = 600.For the clean systems, the results of different bond dimensions start to deviate when t n ≈ 10 for all system sizes, which indicates the bond dimension is saturated.Therefore we only reserve the simulation result before t n = 10.While for the disordered system, the difference between different bond dimensions is not significant, indicating slow entanglement growth.For efficiency consideration, we cut off the simulation at t n = 20.where σ 0 is the constant determining the width of DoS(E) in Eq.A12.With that one could obtain Let's estimate the scale of the correction.Notice that Combining all the derivatives, one obtains One can in addition replace Ē with its expression in Eq. (B4), Altogether we achieve The indicator function of FDT defined in the mean text is To obtain the indicator function, we expand ln S Notice that the terms with an even power of ω cancel out in the indicator function, The first correction term ∂ E ln f O ( Ē, ω) 2 scales as O(1/N ).As the ensemble width scales as σ = O( √ N ) in our simulation, the second correction term is also O(1/N ).All the corrections vanish in the thermodynamic limit, meaning that β ρσ(E) F DT converges to the thermal β in the thermodynamic limit.Notice that a similar derivation was in [67], which is in line with our result.

The generalized spectral function
In this section, we are going to analyze the numerical error of replacing the δ function in the spectral function with a Gaussian filter with width σ ω .The generalized spectral function for an ensemble ρ can be written as where K is a positive constant.Then we have The error is of the order O(σ ω ).
2. For the observable O fulfilling ETH, S ρ O (ω) should be a smooth function of ω.Using the property of the Gaussian function, we have The relative error is of the order O(σ 2 ω ).In general, we expect the smoother S ρ O (ω) is, the smaller this error would be.This Hamiltonian features an inherent particle-hole symmetry: for any given eigenstate |{n µ }⟩, if we define S as the operation that interchanges particles and holes, then S |{n µ }⟩ remains an eigenstate with an energy of the opposite sign.This results in the symmetry in the matrix elements of the observable Ô = σ z N/2 .To be concrete, one could check that the matrix elements of The symmetry in matrix elements immediately implies V O (E, ω) = V O (−E, −ω).

Figure 1 .
Figure 1.A graphical illustration of Eq. 7. The heatmap shows the matrix elements |O αβ | 2 in the eigenstate basis.The x-axis and y-axis are the eigenenergies Eα and E β .The shadowed stripes indicate the filters on energy Eα and energy difference E β − Eα.The summation occurs over the matrix elements within the intersection of two shadowed stripes.

Figure 3 .
Figure 3. Generalized spectral functions S ′ρσ (E) O (ω) for Ô = σ z N/2 , as a function of the energy difference ω, for a filter ensemble of width σ = 0.2 √ N , at energy E/N = 0.5, with N = 40 and varying σω = 0.6, 0.3, 0.1.The panels show the results for (a) integrable, (b) non-integrable, (c) disordered model.Notice that in the disordered system we impose σω ≥ 0.3 to ensure the smallness of the numerical error (see App. A).
shows, for each set of Hamiltonian parameters, the autocorrelator as a function of time for a chain of N = 40 sites and an ensemble
Tr[e iHt ] = dEe iEt DoS(E) is the Fourier transform of DoS(E) into the time domain.Therefore Tr[e iHt ] is also Gaussian, with width proportional to 1/ √ N .As a result, the time for Tr[e iHt ] to fall below 10 −5 of the initial value also scales as O(1/ √ N ).According to Eq. (15), this corresponds to a filter width σ = O( √ N ).