Time evolution of correlation functions in quantum many-body systems
Álvaro M. Alhambra Jonathon Riddell Luis Pedro García-Pintos







aalhambra@perimeterinstitute.ca 
Perimeter Institute for Theoretical Physics, Waterloo, ON N2L 2Y5, Canada


riddeljp@mcmaster.ca
Department of Physics & Astronomy, McMaster University
	1280 Main St.  W., Hamilton ON L8S 4M1, Canada.


Corresponding author
lpgp@umd.edu
Department of Physics, University of Massachusetts, Boston, MA 02125, USA
Joint Center for Quantum Information and Computer Science, NIST/University of Maryland, College Park, Maryland 20742, USA
Joint Quantum Institute, NIST/University of Maryland, College Park, Maryland 20742, USA



We give rigorous analytical results on the temporal behavior of two-point correlation functions –also known as dynamical response functions or Green’s functions– in closed many-body quantum systems. We show that in a large class of  translation-invariant  models the correlation functions factorize at late times ⟨ A(t) B⟩_β→⟨ A ⟩_β⟨ B ⟩_β, thus proving that dissipation emerges out of the unitary dynamics of the system. We also show that  for systems with a generic spectrum  the fluctuations around this late-time value are bounded by the purity of the thermal ensemble, which generally decays exponentially with system size.
For auto-correlation functions we provide an upper bound on the timescale at which they reach the factorized late time value. Remarkably, this bound is only a function of local expectation values, and does not increase with system size. We give numerical examples that show that this bound is a good estimate in non-integrable models, and argue that the timescale that appears can be understood in terms of an emergent fluctuation-dissipation theorem. 
Our study extends to further classes of two point functions such as the symmetrized ones and the Kubo function that appears in linear response theory, for which we give analogous results. 



Time evolution of correlation functions in quantum many-body systems
    Luis Pedro García-Pintos
    December 30, 2023
====================================================================




Two-point correlation functions –or dynamical response/Green's functions– are the central object of the theory of linear response <cit.>, and appear in the characterization of a wide range of non-equilibrium and statistical phenomena in the study of quantum many-body systems and condensed matter physics <cit.>.
This includes different types of scattering and spectroscopy experiments <cit.>, quantum transport <cit.>, and fluctuation-dissipation relations <cit.>.
They have also appeared in the characterization of topological <cit.> and crystalline ordering <cit.>, of quantum-chaotic systems <cit.> and of different notions of ergodicity in quantum and classical systems <cit.>. 
 

Here we study the time evolution of such correlation functions in isolated systems evolving under unitary dynamics. More precisely, we focus on functions of the form 

    C^AB(t)   ≡⟨ A(t) B ⟩_β =ρ  A(t) B,

where the evolution is generated by a time-independent Hamiltonian H, ρ≡  e^-β H / Z_β is a thermal state at inverse temperature β with partition function Z_β, and A(t) = e^i HtAe^-i Ht is the evolved observable in the Heisenberg picture. Both A and B are usually taken to be either local (such as a single-site spin) or extensive operators (such as a global current or magnetization).

Two-point correlation functions have been widely studied before, mostly through numerical methods such as exact diagonalization <cit.>, QMC <cit.> and tensor networks <cit.>, and analytically 
 for specific models, e.g.  <cit.>. Also, a number of experimental schemes to measure it directly have been proposed <cit.>, which manage to circumvent the obstacle of having to measure two non-commuting observables on a single system.
Here, we give rigorous analytical results on their dynamical behavior with as few assumptions on the Hamiltonian as possible. 



Our results apply to most  translation-invariant  non-integrable Hamiltonians, in which the degeneracy of the energy spectrum is small.

First, for arbitrary local observables A and B we 

prove that, for late times, the following signature of dissipation occurs in a large class of translation-invariant  models 

    ⟨ A(t) B ⟩_βt→∞⟶⟨ A ⟩_β⟨ B ⟩_β.

Moreover, we show that the fluctuations around the late-time value are in fact
 bounded by the effective dimension d^-1_eff≡ρ^2 of the ensemble, which decays quickly with system size.
 
For the 
 case of auto-correlation functions, when A=B, we also derive an upper bound on the timescale at which the factorization of Eq. (<ref>) happens, which, remarkably, is independent of the size of the system.  
 We provide numerical evidence showing that the bound is in fact a good estimate even for moderate system sizes, and becomes tighter as the size increases. 

Our study can be extended to a large class of 2-point correlation functions. 
For instance, for 
 symmetric correlation functions C_s^A(t) ≡⟨ A(t) A ⟩_β + ⟨ A A(t) ⟩_β, 
 we find that  
  evolution is dominated by a timescale which is at most of 
  order  t^2 ∼⟨ A^2⟩_β/⟨ [A,H][H,A] ⟩_β. We argue that this can be interpreted in terms of a fluctuation-dissipation theorem that arises from the unitary dynamics of the system. Finally, we consider the timescales of evolution of the Kubo correlation function that appears in linear response theory <cit.>, which dictates the response of a system at equilibrium to a perturbation in its Hamiltonian. 



Late-time behaviour — 
We now show the rigorous formulation of the late-time factorization of  2-point functions.
First, we need the following definition.


	A state ρ on an Euclidean lattice ℤ^D has finite correlation length ξ >0 if it holds that
	
    max_X∈ M,Y ∈ N|ρ X⊗ Y -ρ Xρ Y|/|| X |||| Y ||≤ e^-dist(M,N)/ξ,

	where M,N are regions on the lattice separated by a distance of at least  dist(M,N), and X,Y are arbitrary operators with support on each region.
	

This condition is generic of thermal states at finite temperature away from a phase transition. It has been 
proven at least for 1D systems <cit.> and arbitrary models above a threshold temperature <cit.>.
 In order to prove factorization at late times,
we focus on systems on states that show clustering of correlations, and 
whose Hamiltonians are k-local, i.e. which can be written as H=∑_j h_j, where h_j couples at most k  closest neighbors.

 

Given that evolution is unitary and the system is finite-dimensional, limits such as lim_t→∞ C^AB(t) are not well-defined.  Hence, we consider the 

  late-time behaviour 
  
  under infinite-time averages of the correlation functions
  lim_T→∞∫_0^T dt/T C^AB(t). 
With these considerations, our first main result is the following. 



	Let H be a k-local, translation-invariant, non-degenerate Hamiltonian on a D-dimensional Euclidean lattice of N sites, and let [ρ,H]=0 be an equilibrium ensemble (such as a thermal state) of finite correlation length ξ>0. Let A,B be local observables with support on  at most 𝒪(N^1/D+1-ν) sites, with ν>0. Then
	
    lim_T→∞∫_0^T C^AB(t) dt/T=   ρ Aρ B
       +𝒪( ξ^2 D/D+1log^2 (N) N^-2/D+1)


 This guarantees that all operators supported on a region with size scaling like any function smaller than 𝒪(N^1/D+1), satisfy the assumptions of the theorem.
The proof, found in Appendix <ref>, relies on a weak form of the Eigenstate Thermalization Hypothesis (ETH)
shown in <cit.>, which is itself based on previous works on large deviation theory for lattice models <cit.>. 
This shows that, in fact, any model obeying the weak ETH and without too many degeneracies will display identical factorization of correlation functions at long times <cit.>.

 Note that we assume that the energy spectrum is non-degenerate, which is accurate for systems without non-trivial symmetries or extensive number of conserved quantities. In particular, non-integrable systems usually display Wigner-Dyson statistics in their fine-grained spectrum, which imply level repulsion <cit.>.

This factorization of the correlation function can be thought of as a signature of the emergence of dissipation due to unitary dynamics, since the lack of correlations at different times indicates the loss of information about an initial perturbation of B at time t=0, as reflected in the observable A at time t <cit.>.


Fluctuations around late-time value —

For most times, the 2-point correlation function is in fact close to its late-time average, with small fluctuations around the equilibrium value.
 In order to prove this, one needs the extra assumption that the energy gaps are non-degenerate, which is again reasonable in non-integrable systems with connected Hamiltonians <cit.>,
where it is generally expected to hold as random perturbations are sufficient to lift degeneracies in energy gaps <cit.>. 
  

Let us define C^AB_∞=lim_T→∞∫_0^T dt/T C^AB(t), and the average fluctuations around the late-time value as

    σ_C^2=lim_T →∞∫_0^T dt/T( C^AB(t)-C^AB_∞)^2.

The following result puts an upper bound on average fluctuations.


	Let H=∑_j E_j |j⟩⟨j| be a Hamiltonian with non-degenerate energy gaps, such that
	
    E_j-E_k=E_m-E_l   ⇔ j=m  ,k=l,

and let [ρ,H]=0. It holds that  





    σ_C^2 ≤    A Bmax_j ≠ k{| A_kj B_jk| }ρ^2  ,

where A_kj, B_jk are matrix elements in the energy eigenbasis,  A = ∑_jk A_jk|j⟩⟨k|.


The proof can be found in Appendix <ref>. It follows the same steps as the main result in <cit.>. 
Here, we also find that the purity ρ^2 of the equilibrium ensemble plays a key role. For a microcanonical ensemble (𝕀 / d)^2=1/d, so the RHS of Eq. (<ref>) is expected to decay exponentially with system size in most situations of interest. Also, notice that for a thermal state ρ_β^2≤ 1/Z_β. Moreover, the ETH predicts that | A_kj B_jk|∼ 1/d <cit.>. 
 


Timescales of equilibration — 
Theorems <ref> and <ref> combined imply that correlation functions of the form ⟨ A(t) B ⟩_β are, for most times t ∈ [0,∞],

 close to the uncorrelated average ⟨ A ⟩_β⟨ B ⟩_β, for a wide class of  translation-invariant systems.
It is expected that the timescale at which this happens may depend on a number of factors, such as the distance between A and B. If the operators are far apart on the lattice the correlations are limited by the Lieb-Robinson bound <cit.>, and timescales associated with ballistic (∝ N^1/D) or diffusive (∝ N^2/D) processes may play a role. However, for the autocorrelation function C^A(t) ≡⟨ A(t)A ⟩_β, we can show that  equilibration to the late-time value occurs in a short timescale, independent of system size. 

There may also be further effects at larger timescales, such as the Thouless time <cit.>, and for those effects our result limits their relative size.









Let us define ρ = ∑_j ρ_jj|j⟩⟨j|, so that ρ_jj and A_jk are the matrix elements of ρ and A in the energy basis. We can then write

    C^A(t)/C^A(0)   = ∑_j kρ_jj |A_jk|^2/C^A(0) e^-i(E_j - E_k)t≡∑_α v_α e^-iG_α t,

where we denote each pair of levels {i,j} by a Greek index, and the corresponding energy gaps by G_α≡ E_j-E_k (notice that both E_j-E_k and E_k-E_j appear in the sum).   The normalized distribution v_α≡ρ_jj |A_jk|^2/C^A(0) is central to our proofs, since it contains all the relevant information about the state, observable and Hamiltonian, and determines which frequencies contribute to the dynamics of the autocorrelation function.  Based on it, we define the following functions.


	Given a normalized distribution p_α over  energy gaps G_α, we define ξ_p(x) as the maximum weight that fits an interval of energy gaps with width x:
	
    ξ_p(x) ≡max_λ∑_α : G_α∈[ G_λ, G_λ + x ]  p_α.

We also define

    a(ϵ) ≡ξ_p(ϵ) /ϵσ_G,     δ(ϵ) ≡ξ_p(ϵ),

where σ_G = √(∑_α p_α G_α^2 - ( ∑_α p_α G_α)^2 ) is  the standard deviation of the distribution p_α over the energy gaps G_α.  


 The important point behind these definitions is that, for a sufficiently smooth and unimodal probability distribution, one can find an ϵ small enough such that a(ϵ)∼𝒪(1) and δ(ϵ)≪ 1.
In the following theorem, the relevant probability distribution is given by v_α. 
Our main result regarding the timescales of correlation functions, proven in Appendix <ref>,  is:


For any Hamiltonian H and state ρ such that [H,ρ] = 0, and any observable A, the autocorrelation function C^A(t) = ρ A(t) A satisfies

    1/T∫_0^T  |C^A(t) - C^A_∞|^2  /(C^A(0))^2 dt     ≤ 3 π( a(ϵ)/σ_G1/T + δ(ϵ) )

 for any ϵ>0. Here, a(ϵ) and δ(ϵ) are as in Definition <ref> for the normalized distribution v_α≡ρ_jj |A_jk|^2/C^A(0), and 
 σ_G is given by

    σ_G^2    =  1 / C^A(0) ρ [A ,H] [H,A ]  - ρ [H,A ] A ^2/(C^A(0))^2.



Theorem [<ref>] provides an upper bound of T_eq≡3π  a(ϵ)/σ_G on the timescales under which autocorrelation functions approach their steady state value.  
To see this note that, if for a given T the RHS of Eq. (<ref>) is small, C^A(t) must have spent a significant amount of time  during the interval [0,T] near the late-time value C^A_∞ .


 The crucial point is that   for distributions v_α that are uniformly spread over many values of the gaps G_α, one can always find an ϵ such that δ≪ 1. In that case, the right hand side of Eq. (<ref>) becomes small on timescales 𝒪(T_eq). 
As discussed in <cit.> and in Appendix <ref>, if one further assumes smooth unimodal distributions, typically  one also finds that a ∼𝒪(1). In that case,  the timescale is governed by 1/σ_G. 
 Given that σ_G is a combination of expectation values of local observables, it does not scale with the system of the system.
In fact, a result of <cit.> shows that a timescale of order 1/σ_G provides a lower bound to the timescale of equilibration, which strongly suggests that our upper bound is tight when the conditions of a∼𝒪(1) and δ≪ 1 hold.

As a prime example, for local operators in non-integrable lattice models, in which (as per the ETH) | A_jk| are uniformly distributed around a peak at zero energy gap <cit.>, one should be able to choose ϵ such that a ∼𝒪(1) and δ≪ 1. In Fig. <ref> we numerically show that this is indeed the case in a non-integrable Ising model.


Theorem <ref> does not make assumptions on the specifics of the Hamiltonian, the observable or the state,  making it completely general. 
However, we do not expect the correlation functions to equilibrate well in all cases, as in some scenarios a(ϵ) and δ(ϵ) will be large  no matter what value of ϵ is chosen, in which case the RHS of Eq.  (<ref>) may not become small within reasonable timescales.  This can happen, for instance, in models with highly degenerate energy spectrum. 
 To illustrate this, in Appendix <ref> we compute  these parameters in an integrable model, where we see that the gap degeneracies of the model negatively affect the quantities a(ϵ) and δ(ϵ), making the estimated equilibration timescales longer.
 
 


Symmetric correlation functions — 
The previous results can be extended to other correlation functions, such as

    C^A_s(t) ≡1/2ρ{ A , A(t) }  = C^A(t) + C^A(t)^*/2.

Along the same lines of Theorem [<ref>], in Appendix <ref> we prove the following.


For any Hamiltonian H and state ρ such that [H,ρ] = 0, and any observable A, the time correlation function C_s^A(t) = ρ{A, A(t)} satisfies

    1/T∫_0^T  |C^A_s(t) - C^A_s,∞|^2  /(C^A_s(0))^2 dt     ≤ 3 π( a(ϵ)/σ_G1/T + δ(ϵ) ),

 for any ϵ>0. Here,  a(ϵ) and δ(ϵ) are as in Definition  [<ref>] for the normalized distribution v_α^s ≡ρ_jj + ρ_kk/2|A_jk|^2/C^A_s(0), and  

    σ_G^2    =  1 / C^A_s(0) ρ [A_0,H] [H,A_0].


Thus an upper bound for the equilibration timescale is

    T_eq   =   3 π  a(ϵ)  √(C^A_s(0))/√(ρ [A ,H] [H,A ]),

where again  
 we expect that for small enough ϵ, a(ϵ)∼𝒪(1) and δ≪ 1 for the same reasons as before. 
The denominator in T_eq can be seen as an “acceleration” of the symmetric autocorrelation function. 
Eq. (<ref>) can in fact be written as

    T_eq =     3 π  a(ϵ)  √(C^A_s(0))/√(| d^2 C^A_s(t) /dt^2|_0    | ).



Such timescale turns out to be similar to that of a short-time analysis. A Taylor expansion gives  

    C^A_s(t)    =C^A_s(0) (  1 -    1/2 C^A_s(0) d^2 C^A_s(t) /dt^2|_0   t^2  )  +𝒪(t^3).

For early times, the above expression decays on a timescale τ =  √(2)/3 π a(ϵ) T_eq, 
identical to our upper bound Eq. (<ref>)
up to a prefactor. 

The timescale of Eq. (<ref>) suggests an interpretation in terms of an emergent fluctuation-dissipation theorem. Consider i) T_eq to be the timescale of dissipation of unitary dynamics, meaning that ⟨ A(t) A ⟩_β⟶⟨ A ⟩_β⟨ A ⟩_β occurs, and ii) C^A_s(0)=ρ A^2 as a measure of the  fluctuations of A. Then, Eq. (<ref>) gives a proportionality relation between the strength of the fluctuations and the timescale of equilibration, in a similar spirit to what was found in <cit.> using random matrix theory arguments.




Linear response and the Kubo correlation function —
As a further application of our methods, we study the evolution of a quantum system under a perturbation of its Hamiltonian. Let the system start in a thermal state, such that ρ∝ e^-β (H+λ A). Subsequently, the Hamiltonian is slightly perturbed by λ A,  so that the evolved state is ρ_t = e^-i t Hρ e^i t H. 

It was shown by Kubo <cit.> that, to leading order in λ, the expectation value of A satisfies ρ A(t) = C_Kubo(t)ρ A, where




for thermal initial  states ρ  the Kubo correlation function can be written as

    C_Kubo(t) ∝∑_j≠ ke^-β E_k-e^-β E_j/E_j-E_k| A_jk|^2 e^i t(E_j-E_k).



Equilibration of ρ A(t) is then equivalent to equilibration of the function (t), for which we prove in Appendix <ref>

that the following holds.



For any Hamiltonian H, thermal state ρ∝ e^-β (H+λ A), and any observable A, the Kubo correlation function C_Kubo satisfies

    1/T∫_0^T  |C_Kubo(t) - C_Kubo,∞ |^2  /C_Kubo(0)^2 dt     ≤ 3 π( a(ϵ)/σ_G1/T + δ(ϵ) ),

 for any ϵ>0. Here,  a(ϵ) and δ(ϵ) are as in Definition  [<ref>] for the normalized distribution w_α≡e^-β E_k-e^-β E_j/E_j-E_k| A_jk|^2 /C_Kubo(0), and 

    σ_G^2    =  1 /(0)  [A,ρ] [A,H].


This again implies an upper bound T_eq =     3 π  a(ϵ)/σ_G on the equilibration timescale of , 
and therefore on the time to return to thermal equilibrium after a perturbation of the system Hamiltonian by A. 
 The distribution w_α plays the same role as v_α and v^s_α before. If w_α is smoothly distributed and unimodal (which we expect for local observables in non-integrable models)  then a ∼𝒪(1) and δ(ϵ)≪ 1 holds (see Appendix <ref>).








































Simulations —
We test Theorem [<ref>] in a spin model governed by the Hamiltonian

    H = ∑_j=1^L(γσ_j^X + λσ_j^Z   ) + J∑_j=1^L-1σ_j^Zσ_j+1^Z+ α∑_j=1^L-2σ_j^Zσ_j+2^Z,

where σ_j^Z and σ_j^X are the Pauli spin operators along Z and X directions for spin j, and
we take open boundary conditions.
The field and interaction coefficients (γ,λ, J, α ) characterize the model. 
 We focus on a case corresponding to a system satisfying ETH by choosing (γ,λ, J, α ) = (0.8,0.5,1,1) <cit.>,
 
  and study the autocorrelation functions of the observable A = σ_L/2^x. For simplicity we set β = 1 in our numerics, though no significant changes were observed for β∈ [0.1,5].
Figure <ref> depicts the functions a(ϵ) and δ(ϵ) that appear in Theorem <ref>,  confirming that there exist regions of ϵ such that δ≪ 1, ensuring equilibration occurs, and a ∼ 0.4. Importantly, this is increasingly the case as the size of the system grows.  
 
 
Figure <ref> compares the two sides in bound (<ref>), 
showing that dynamics obtained from the upper bound differs from the actual dynamics by roughly an order of magnitude.
Thus, the general, model-independent bounds obtained from Theorem <ref> provide remarkably good estimates of the actual (simulated) dynamics.
Note that the estimate becomes increasingly better as the size of the system increases. 
This discrepancy could, however, be a finite-size effect, which is also suggested by the lower bound obtained in <cit.>.
Details of the simulations can be found in <ref>.
 




 


Discussion — 
We derived analytic results on the dynamical behavior of 2-point correlation functions in quantum systems. These include conditions that imply that time-correlation functions factorize for long times, as well as easy-to-estimate upper bounds on the timescales under which such process occurs
 which hold regardless of details of the model under consideration.
Remarkably, 
our numerical findings show that the
derived upper bounds can correctly estimate the actual dynamics of the system to within an order of magnitude, and
become increasingly better
  estimates as the size of the system increases. 



 




We used techniques previously applied in the context 
of equilibration of quenched quantum systems <cit.>, for which finding rigorous estimates on the timescales is a largely open problem <cit.>. This connection is not surprising, specially considering that previous works <cit.> have argued that in some situations one can approximate the out of equilibrium dynamics with the autocorrelation functions covered here. 

Given the importance of time-correlation functions in the analysis of a wide range of problems in many-body physics –for instance, in transport phenomena– we anticipate that our results will be useful in the description of closed system dynamics, whose study has surged in recent times due to enormous experimental advances in settings such as cold atoms or ion traps <cit.>.



The authors acknowledge useful discussions with Anurag Anshu, Beni Yoshida, Charlie Nation and Jens Eisert. We are also thankful to Yichen Huang for pointing out an improvement to the error term in Theorem <ref> from a previous version (see <cit.>). This research was supported in part by the Perimeter Institute for Theoretical Physics. Research at
Perimeter Institute is supported by the Government of
Canada through the Department of Innovation, Science
and Economic Development and by the Province of Ontario
through the Ministry of Research, Innovation and
Science. This research was also supported by NSERC and 
enabled in part by support provided by (SHARCNET) (www.sharcnet.ca) and Compute/Calcul Canada (www.computecanada.ca).
LPGP acknowledges support from the John Templeton Foundation, UMass Boston Project No. P20150000029279,
DOE Grant No. DE-SC0019515, AFOSR MURI project "Scalable Certification of Quantum Computing Devices and Networks", DoE ASCR Quantum Testbed Pathfinder program (award No. DE-SC0019040), DoE BES QIS program (award No. DE-SC0019449), DoE ASCR FAR-QC (award No. DE-SC0020312), NSF PFCQC program, AFOSR, ARO MURI, ARL CDQI, and NSF PFC at JQI.






§ LATE TIME BEHAVIOUR OF TWO-POINT FUNCTIONS



 §.§ Proof of late-time equilibration


First, we state the key result from <cit.> that we use.  It says that expectation values for single eigenstates, of the form ⟨E_k| A |E_k⟩, are close to the ensemble average ⟨ A ⟩_β with very high probability. In contrast, the strong form of the ETH states that the above happens for all eigenstates within an energy window. We reproduce the proof of <cit.> (which itself builds on <cit.>), with the difference that our version holds for lattices of dimension larger than 1.


	Let H be a translation-invariant, non-degenerate Hamiltonian with N sites on a D-dimensional lattice, ρ an equilibrium ensemble [ρ,H]=0 with finite correlation length ξ, and A some observable with support on a connected region of at most 𝒪(N^1/(D+1)-ν) sites, with ν>0. Then, for any δ > 0,
	
    Pr_|E_k⟩∈ρ (|⟨E_k| A |E_k⟩-ρ A|≥δ )≤exp(-c δ N^1/D+1ξ^-D/D+1),

	where c>0 is a constant, and |E_k⟩∈ρ  indicates that the eigenstates are sampled from the equilibrium distribution ρ.


	We show a bound Pr_|E_k⟩∈ρ ( ⟨E_k| A |E_k⟩-ρ A≥δ ), and Pr_|E_k⟩∈ρ ( ρ A- ⟨E_k| A |E_k⟩≥δ ) will follow analogous steps. Notice that since the Hamiltonian is translation-invariant and non-degenerate we can write ⟨E_k| A |E_k⟩=⟨E_k|A̅/N  |E_k⟩,
	where A̅=∑_i A_i is the extensive observable built out of translations of A. Define A̅'=A̅-⟨A̅⟩_β. Then, using Markov's inequality and e^⟨ψ| A |ψ⟩≤⟨ψ| e^A |ψ⟩  we can write
	
    Pr_|E_k⟩∈ρ ( ⟨E_k| A |E_k⟩-ρ A≥δ )    = Pr_|E_k⟩∈ρ(e^λ⟨E_k|A̅' |E_k⟩≤ e^λδ N) 
       ≤ e^- λδ N ∑_k ρ_kk e^λ⟨E_k|A̅' |E_k⟩
       ≤ e^- λδ N ρ e^λA̅'.

	Now let us write the spectral decomposition e^λA̅'=∑_l e^λ a_lΠ_l, with Π_l the projectors into the subspace with eigenvalue e^λ a_l, and write the average as
	
    ρ e^λA̅'=∑_a_l ≤δ N /2   e^λ a_lρΠ_l + ∑_a_l > δ N /2  e^λ a_lρΠ_l.

	The first term is upper bounded by e^λδ N/2. For the second, we write
	
    ∑_a_l > δ N /2  e^λ a_lρΠ_l =∑_j ∑_δ N/2+(j+1) ≥ a_l ≥δ N/2+j e^λ a_lρΠ_l≤∑_j e^λ (δ N  /2+j+1)ρΠ_≥δ N/2+j,

	where Π_≥ x denotes the projector on the subspace with a_l ≥ x. The main result of <cit.> states that
	
    ρΠ_≥δ N/2+j≤ e^D k N^α /ξexp(- c'/D ξ(N (δ/2+j/N)^2 ξ)^1/D+1).

	Here, k N^α is the number of sites on which an individual local term of observable A̅ has support.
	Since j is at most 𝒪(N), we can choose some λ =𝒪((N ξ)^-D/D+1) such that, for some constant c_1>0,
	
    ∑_a_l > δ N /2  e^λ a_lρΠ_l≤poly(N) e^D k N^α /ξexp(-c_1 N^1/D+1ξ^-D/D+1).

	This way, the dominant contribution of Eq. (<ref>) is the first term. Plugging the bounds back in Eq. (<ref>) results in the following, for some constant c>0 and large enough N,
	
    Pr_|E_k⟩∈ρ ( ⟨E_k| A |E_k⟩-ρ A≥δ )    ≤ e^-λδ N/2(1 + e^-λδ N/2poly(N) e^D k N^α /ξexp(-c_1 N^1/D+1ξ^-D/D+1)) 
       ≤exp(-c δ N^1/D+1ξ^-D/D+1),

	where the last line follows from the fact that the second term in the RHS of Eq.(<ref>) is subleading (much smaller than 1) as long as α<1/D+1. This sets the support of of the observable A ot be within a region of size at most 𝒪(N^1/(D+1)-ν) for any ν>0. 

With it, we are now ready to prove the result on late-time factorization of correlation functions.
1
	
	Let H be a k-local, translation-invariant and non-degenerate Hamiltonian on a D-dimensional Euclidean lattice of N sites, and let [ρ,H]=0 be an equilibrium ensemble (such as a thermal state) of finite correlation length ξ>0. Let A,B be local observables with support on at most 𝒪(N^1/D+1-ν) sites, with ν>0. Then
	
    lim_T→∞∫_0^T ρ A(t) Bdt/T=ρ Aρ B+𝒪(ξ^2 D/D+1log^2 (N) N^-2/D+1)



	First let us write
	
    ρ=∑_k ρ_k k|E_k⟩⟨E_k|  ,   A=∑_k j A_k j|E_k⟩⟨E_j|   ,    B=∑_k j B_k j|E_k⟩⟨E_j|,

	from which we have
	
    ρ A(t) B=∑_k j A_k j B_j kρ_kk e^i t (E_j-E_k).

	Since the Hamiltonian is non-degenerate by assumption, it holds that 
	
    lim_T→∞∫_0^T e^i t (E_j-E_k)d t/T=δ_k,j,

	so that the limit becomes
	
    lim_T→∞∫_0^T ρ A(t) Bdt/T=∑_k A_k k B_k kρ_kk.

	Now let us define A_k k-ρ A=Δ_k,A and B_k k-ρ B=Δ_k,B , so that <cit.> 
	
	
    ∑_k A_k k B_k kρ_kk   =ρ Aρ B+∑_k ρ_k k(ρ AΔ_k,B+ρ BΔ_k,A+Δ_k,AΔ_k,B)
        =ρ Aρ B+∑_k ρ_k kΔ_k,AΔ_k,B

	Let us define Δ≡ K log N/ N^1/D+1, and split the sum over energies of the error term as
	
    ∑_k ρ_k kΔ_k,AΔ_k,B=∑_k ∈𝒮ρ_k kΔ_k,AΔ_k,B+∑_k ∉𝒮ρ_k kΔ_k,AΔ_k,B.

	where 𝒮≡{k: |Δ_k,A|,|Δ_k,B|≤Δ} (that is, the set of k for which both errors are small).
	Notice that the first term is smaller than Δ^2 by definition. On the other hand, the second term can be bounded as
	
    ∑_k ∉𝒮ρ_k kΔ_k,AΔ_k,B   ≤max_k'Δ_k',AΔ_k',B∑_k ∉𝒮ρ_k k
       ≤max_k'Δ_k',AΔ_k',B(∑_|Δ_k,A|≥Δρ_k k+∑_|Δ_k,B|≥Δρ_k k )
       ≤ 2max_k'Δ_k',AΔ_k',Bexp(-c Δ N^1/D+1ξ^-D/D+1)
       ≤ 2 || A |||| B ||( 1/N)^c Kξ^-D/D+1.

	The third line follows from Lemma <ref>, and the fourth from |Δ_k,A|≤|| A ||,|Δ_k,B|≤|| B ||. The constant K is arbitrary, so we can choose it such that c Kξ^-D/D+1 = 2/D+1. In that case the dominant contribution to Eq. (<ref>) is that of the first term, and hence ∑_k ρ_k kΔ_k,AΔ_k,B =𝒪(Δ^2), so that
	
    ∑_k A_k k B_k kρ_kk=ρ Aρ B+𝒪(Δ ^2),

	completing the proof.





 §.§ Proof of fluctuations around late-time value



2
	
	Let H=∑_i E_i |i⟩⟨i| be a Hamiltonian with non-degenerate energy gaps, such that
	
    E_j-E_k=E_m-E_l   ⇔ j=m  ,k=l,

	and let [ρ,H]=0. It holds that  
	
	
	
	
	
    σ_C^2 ≤    A Bmax_j ≠ k{| A_kj B_jk| }ρ^2  ,

	where A_kj, B_jk are off-diagonal matrix elements in the energy eigenbasis,  A = ∑_jk A_jk|j⟩⟨k|.



	Let us expand in the energy eigenbasis.
	
	
	
	
	
	
	
	
	
    σ_C^2    =lim_T →∞∫_0^T dt/T (C(t)^AB-C_∞^AB)^2 
       = ∑_j≠ k∑_l≠ mρ_jjρ_ll A_jkB_kj A_lm B_ml∫_0^T dt/T e^-i t (E_j-E_k+E_l-E_m)
       = ∑_j≠ kρ_jjρ_kk A_jkA_kj B_jkB_kj
       ≤∑_j≠ kρ_jjρ_kk| A_jk A_kj B_jkB_kj|
    	
       ≤max_j ≠ k{| A_kj B_jk| }∑_j≠ kρ_jjρ_kk | A_jkB_kj|

	In the second to the third line we used the assumption of non-degenerate energy gaps. 
	We now use the Cauchy-Schwarz inequality twice.
	
    σ_C^2    ≤max_j ≠ k{| A_kj B_jk| }√(∑_j≠ kρ_jjρ_kk | A_jk|^2 )√(∑_j≠ kρ_jjρ_kk | B_kj|^2 )
       ≤max_j ≠ k{| A_kj B_jk| }√(∑_j,kρ_jjρ_kk | A_jk|^2 )√(∑_j,kρ_jjρ_kk | B_kj|^2 )
       =   max_j ≠ k{| A_kj B_jk| }√(ρ A ρ A)√(ρ B ρ B)
       ≤max_j ≠ k{| A_kj B_jk| }√(A^2 ρ^2 )√(B^2 ρ^2 )
       ≤A Bmax_j ≠ k{| A_kj B_jk| }ρ^2 .

	
	The last inequality follows from the fact that for positive operators PQ≤|| P ||Q. 
	




§ DYNAMICS OF TWO-TIME CORRELATION FUNCTIONS





 §.§ Preliminaries


Here we prove an intermediate lemma required for Appendix <ref>. This is similar to results found in  <cit.>, improving on them by a numerical factor. Let us first recall the following definition from the main text. 

2
	
	Given a normalized distribution p_α over G_α, we define ξ_p(x) as the maximum weight that fits an interval of energy gaps with width x:
	
    ξ_p(x) ≡max_λ∑_α : G_α∈[ G_λ, G_λ + x ]  p_α.

   
With it, we have a general upper bound that will be central to the later proofs.


	
	Let f(t) be a positive function of the form 
	
    f(t) = ∑_α,βp_α p_β e^i(G_α- G_β)t,

	such that the p_α  form a discrete probability distribution ∑_α p_α = 1. The uniform average of such a function is upper bounded by, 
	
    ⟨ f(t) ⟩_T ≡∫_0^T dt/T f(t)  ≤ 3 πξ_p(1/T),

	where ξ_p(x) is as in Definition <ref>. 


	Let us first define the uniform and Gaussian probability density functions, as 
	
    p_T(t) =  {[        1/T  t∈ [0,T],;          0 otherwise. ].
 
	and
	
    p_G(t)  = 	1/√(2πα^2T^2) e^-(t-T/2)^2/2α^2T^2,  t∈ R,

	respectively, 
	where the mean of the Gaussian is written as μ = T/2 and the standard deviation as σ = α T. The parameter α>0 is free for us to choose.
	We will first show that for any positive function f(t)  ≥ 0, the uniform average  of a function on the interval t∈ [0,T] can be bounded tightly by the Gaussian average by, 
	
    ⟨  f(t) ⟩_T ≤γ⟨ f(t) ⟩_G_T,

	where γ = √(2π)α e^1/8α^2, and we denote the uniform and Gaussian averages over distributions p_T(t) and p_G(t) as ⟨  f(t) ⟩_T and ⟨ f(t) ⟩_G_T respectively. 
	
	To do this, note that due to the choice of μ = T/2 the two distributions have identical means. To find the smallest possible γ such that (<ref>) holds we set that p_T(t=T) = γ p_G(t=T). 
	This then implies,
	
    1/T = γ1/√(2 π)α T e^-1/8 α^2γ = √(2π)α e^1/8α^2,

	which guarantees that p_T(t) ≤γ p_G(t)  ∀ t. Since over the interval [0,T] p_T(t) is constant and p_G(t)≥ p_G(T),it follows that
	
    ⟨  f(t) ⟩_T ≤γ⟨ f(t) ⟩_G_T.

	With this we can now bound the uniform average of f(t) as
	
    ⟨ f(t) ⟩_T = 1/T∫_0^T ∑_α, βp_α p_β e^i(G_α- G_β)tdt ≤γ/√(2πα^2T^2)∫_-∞^∞∑_α, βp_α p_β e^i(G_α- G_β)t e^-(t-T/2)^2/2α^2T^2 dt.

	Interchanging the integral with the sum and then integrating gives the characteristic function of a Gaussian distribution. This allows us to write, 
	
    ⟨ f(t) ⟩_T ≤γ∑_αβ p_α p_β e^ - σ^2 Δ G^2/2,

	where Δ G = G_α-G_β. To further bound Eq. (<ref>) we may introduce the function, 
	
    g(x) = {[            1 if  x∈[0,1),;            0   otherwise. ].

	This function satisfies
	
    r^-x^2≤∑_n=0^∞ r^-n^2 g(|x|-n)

	where r > 1. Applying this bound to Eq. (<ref>) with r =  e^α/2 and x = Δ G^2 T^2, it follows that 
	
    ⟨ f(t) ⟩_T ≤γ∑_n=0^∞ r^-n^2∑_α,β p_α p_β g(Δ G^2 T^2-n).

	The presence of g(x) allows to restrict the summation in Eq. (<ref>). With the condition ( Δ G^2 T^2-n ) ∈ [0,1), the two following intervals can be identified for the variable G_β given a particular value of G_α,
	
    G_β∈ I_+ =  [ G_α+ √(n)/T, G_α+√(n+1)/T), 
     G_β∈ I_- = ( G_α -  √(n+1)/T, G_α+√(n)/T].
 
	With this we rewrite the summation in the bound as,
	
    ⟨ f(t) ⟩_T ≤γ∑_n=0^∞ r^-n^2∑_αp_α( ∑_G_β∈ I_+ p_β + ∑_G_β∈ I_- p_β).

	Given that √(n)≤ n for all integer n, one can expand the length of the interval, so that
	
    ⟨ f(t) ⟩_T ≤γ∑_n=0^∞ r^-n^2∑_αp_α( ∑_G_β∈ I'_+ p_β + ∑_G_β∈ I'_- p_β) ,

	where we define
	
    I'_+ =  [ G_α+ n/T, G_α+n+1/T), 
    
    	I'_- = ( G_α -  n+1/T, G_α+n/T].

	It is straightforward to see that ∑_G_β∈ I'_+ p_β + ∑_G_β∈ I'_- p_β≤ 2 ξ_p(1/T), which allows us to further bound the summation by using Def. <ref>,
	
    ⟨ f(t) ⟩_T ≤ 2 γξ_p(1/T) ∑_n=0^∞ r^-n^2.

	To complete the proof we find the α that minimizes
	
    κ(α) = γ∑_n=0^∞ r^-n^2 = γ( Θ_3(0,r^-1)+1  ),

	which occurs for α≈ 0.6347. This gives κ≈ 2.8637<3. Finally,
	
    ⟨ f(t) ⟩_T ≤ 3 πξ_p(1/T),
 
	which completes the proof. 
	
	
	
	











 §.§ Proof of equilibration timescales of correlation functions









3
	For any Hamiltonian H and state ρ such that [H,ρ] = 0, and any observable A, the time correlation function C^A(t) = ρ A(t) A satisfies
	
    1/T∫_0^T  |C^A(t) - C^A_∞|^2  /(C^A(0))^2 dt     ≤ 3 π( a(ϵ)/σ_G1/T + δ(ϵ) )

	for any ϵ>0. Here, a(ϵ) and δ(ϵ) are as in Definition <ref> for the normalized distribution v_α≡ρ_jj |A_jk|^2/C^A(0), and 
	σ_G is given by
	
    σ_G^2    =  1 / C^A(0) ρ [A ,H] [H,A ]  - ρ [H,A ] A ^2/(C^A(0))^2.




	Following Lemma <ref>, we have that
	
    1/T∫_0^T  |C_t - C_∞|^2  /(C^A(0))^2 dt    ≤ 3 πξ_v( 1T),

	where ξ_v( 1T) is taken with respect to the distribution v_α, which is normalized:
	
    ∑ v_α = ∑_jkρ_jj |A_jk|^2/C^A(0) = ρ  A^2/C^A(0) = 1.

	
	It was shown in <cit.> (Proposition 5) that the function ξ_p(x) satisfies 
	
    ξ_p(x) ≤a(ϵ)/σ_G x + δ(ϵ),

	with a(ϵ) = ξ_p(ϵ) /ϵσ_G and δ(ϵ) = ξ_p(ϵ), and where ε∈ (0,∞).
	
	Finally, the standard deviation of the distribution p_v is 
	
    σ_G^2    =  ∑_α p_α G_α^2 - ( ∑_α p_α G_α)^2  
        = ∑_jkρ_jj |A_jk|^2 / C^A(0)  (E_j - E_k)^2   - ( ∑_jkρ_jj |A_jk|^2 / C^A(0)  (E_j - E_k)   )^2 
        =  1 / C^A(0) ρ [A ,H] [H,A ] - ρ [H,A ] A ^2/(C^A(0))^2,

	which completes the proof.






 §.§ Proof of equilibration timescales of symmetric correlation functions


4
	For any Hamiltonian H and state ρ such that [H,ρ] = 0, and any observable A, the time correlation function C_s^A(t) = ρ{A, A(t)} satisfies
	
    1/T∫_0^T  |C^A_s(t) - C^A_s,∞|^2  /(C^A_s(0))^2 dt     ≤ 3 π( a(ϵ)/σ_G1/T + δ(ϵ) ),

	for any ϵ>0. Here, a(ϵ) and δ(ϵ) are as in Definition <ref>
	
	for the normalized distribution v_α^s ≡ρ_jj + ρ_kk/2|A_jk|^2/C^A_s(0), and  
	
    σ_G^2    =  1 / C^A_s(0) ρ [A_0,H] [H,A_0].

	



	The symmetric correlation function is defined as
	
    C^A_s(t) ≡1/2ρ{ A , A(t) }  = C^A(t) + C^A(t)^*/2.

	The equivalent of Eq.<ref> 
	becomes
	
    C^A_s(t)/C^A_s(0)    = ∑_j kρ_jj + ρ_kk/2|A_jk|^2/C^A_s(0)  e^-i(E_j - E_k)t,

	The proof of Theorem <ref>
	
	is identical as the previous proof, with the symmetrized distribution v_α^S ≡ρ_jj + ρ_kk/2 |A_jk|^2.
	In this case the variance of the normalized distribution v_α^S becomes
	
    σ_G^2    =  ∑_α p_α^S G_α^2 - ( ∑_α p_α^S G_α)^2   
        = ∑_jkρ_jj + ρ_kk/2 |A_jk|^2 / C^A_s(0)  (E_j - E_k)^2   - ( ∑_jkρ_jj + ρ_kk/2 |A_jk|^2 / C^A_s(0)  (E_j - E_k)   )^2   
        =  1 / C^A_s(0) ρ [A ,H] [H,A ].



Equilibration then occurs within a timescale

    T_eq   =   κπ  a(ϵ)  √(C^A_s(0))/√(ρ [A ,H] [H,A ]),

The denominator in T_eq can be identified as an “acceleration” of the symmetric autocorrelation function. 
Indeed, 

    d^2 C^A_s(t) /dt^2 = - ∑_j k (E_j - E_k)^2 ρ_jj+ρ_kk/2 |A_jk|^2 e^-i(E_j - E_k)t.

Then, the equilibration timescale is

    T_eq =     κπ  a(ϵ)  √(C^A_s(0))/√(| d^2 C^A_s(t) /dt^2|_0    | ).





  §.§.§ Short-time evolution of symmetric correlation functions


The symmetric autocorrelation functions is given by

    C^A_s(t)    = ∑_j kρ_jj + ρ_kk/2  |A_jk|^2 e^-i(E_j - E_k)t.

Taking the Taylor expansion of C^A_s(t),

    C^A_s(t)    = C^A_s(0)  - ( i ∑_j k (E_j - E_k  ) 
    ρ_jj + ρ_kk/2      |A_jk|^2 ) t 
       - (  ∑_j k (E_j - E_k)^2  ρ_jj + ρ_kk/2  |A_jk|^2  ) t^2/2 +𝒪(t^3) 
       =C^A_s(0)- d^2 C^A_s /dt^2|_0  t^2 /2  +𝒪(t^3) 
       =C^A_s(0) (  1 -    1/2 C_0 d^2 C^A_s /dt^2|_0   t^2  )  +𝒪(t^3).

For early times, the above expression decays on a timescale 

    τ≡√(2)√(C^A_s(0))/√(| d^2 C^A_s(t) /dt^2|_0    | ).

























































































 §.§ Proof of equilibration timescales of Kubo correlation functions



5
	For any Hamiltonian H, thermal state ρ, and any observable A, the Kubo correlation function C_Kubo satisfies
	
    1/T∫_0^T  |C_Kubo(t) - C_Kubo,∞ |^2  /C_Kubo(0)^2 dt     ≤ 3 π( a(ϵ)/σ_G1/T + δ(ϵ) ),

	for any ϵ>0. Here, a(ϵ) and δ(ϵ) are as in Definition <ref> 
	for the normalized distribution w_α≡e^-β E_k-e^-β E_j/E_j-E_k| A_jk|^2 /C_Kubo(0), and 
	
    σ_G^2    =  1 /(0)  [A,ρ] [A,H].




	The Kubo correlation function can be written as
	
    C_Kubo(t) ∝∑_j≠ ke^-β E_k-e^-β E_j/E_j-E_k| A_jk|^2 e^i t(E_j-E_k),

	with the proportionality constant defined by C_Kubo(0) = 1.
	We can then write
	
    C_Kubo(t)-C_Kubo(∞)/C_Kubo(0)   =∑_j ≠ ke^-β E_k-e^-β E_j/E_j-E_k| A_jk|^2 e^i t(E_j-E_k)/∑_j,ke^-β E_k-e^-β E_j/E_j-E_k| A_jk|^2
       ≡∑_α≠ 0 w_α e^-i t G_α/∑_α w_α,

	where we define w_α≡e^-β E_k-e^-β E_j/E_j-E_k| A_jk|^2≥ 0.
	Given that w_α≥ 0, we can perform similar calculations as for Theorem <ref>, 
	albeit with a different probability distribution. 
	Thus, we also have
	
    ⟨(C_Kubo(t)-C_Kubo(∞)/C_Kubo(0))^2 ⟩_T ≤ 3 πξ_w (1/T).

	Defining the normalized distribution q_α≡ w_α/(0). The variance is 
	
    σ_G^2    =  ∑_α q_α G_α^2 - ( ∑_α q_α G_α)^2  
        = ∑_jk1/(0) e^-β E_k-e^-β E_j/E_j-E_k| A_jk|^2  (E_j - E_k)^2   - ( ∑_jk1/(0) e^-β E_k-e^-β E_j/E_j-E_k| A_jk|^2 (E_j - E_k)   )^2 
        = 1/(0)∑_jk(e^-β E_k-e^-β E_j)   | A_jk|^2 (E_j - E_k)    
        =  1 /(0)  [A,ρ] [A,H] ,

	completing the proof.





 §.§ Scaling of a and δ



The proofs of Theorems (3-5)
rely on the fact that the function 
ξ_p(x),
defined for any normalized distribution p_α as  the maximum distribution that fits an interval x 

    ξ_p(x) ≡max_β∑_α : G_β∈[ G_β, G_β + x ]  p_α,

satisfies


    ξ_p(x) ≤a(ϵ)/σ_G x + δ(ϵ),

which was shown in <cit.> (Proposition 5). Here a(ϵ) ≡ξ_p(ϵ) /ϵσ_G and δ(ϵ) ≡ξ_p(ϵ), where

σ_G is the standard deviation of the distribution p_α.
The function a(ϵ) ends up in the bound of the equilibration timescales, as T_eq =     π  a(ϵ)/σ_G, while δ(ϵ) governs the long time behavior in the bounds. 

Given that ξ_p(x) characterizes how much of the distribution p_α fits an interval x, the value of a in Eq. (<ref>) depends on how well 1/σ_G serves to characterize the region where the distribution p_α is supported.
Roughly speaking, whenever 1/σ_G is a good estimate of the width of such small region, then one expects a ∼𝒪(1). This is well illustrated when considering a unimodal distribution (e.g. a Gaussian). In such case, the fraction of the distribution that fits an interval x is roughly x times the width 1/σ_G of the window where the distribution is supported, and ξ_p(x) ∼ x/σ_G, so that a ∼𝒪(1).
Multimodal distributions violate such condition, as for them the standard deviation does not characterize the regions in which the distribution has considerable support.
At the same time, δ in Eq. (<ref>) carries information of the fine structure of p_α, indicating the scale at which the distribution can no longer be coarse-grained to a continuous distribution. The only way that δ≪ 1 fails is for distributions that are not smooth, in which a small region of width ϵ is significantly populated. Thus, for distributions that are smooth in a coarse-grained sense, and approximately unimodal, one expects to be able to find a small enough ϵ  such that δ(ϵ )  ≪ 1 and a(ϵ) ∼𝒪(1). 


In summary, the problem of proving fast equilibration timescales in our approach can thus be linked to knowing whether the relevant distribution p_α is `approximately unimodal'.
We argue that for an strongly interacting many-body system this will typically be the case. 
Consider for instance the case of Theorem <ref>,
where the relevant distribution is given by v_α≡ρ_jj |A_jk|^2/C^A(0).

The large number of energy gaps present in a many-body system implies a dominance of small gaps over larger ones, which favors that, on a coarse-grained sense, the distribution over gaps shows a decay as the size |G_α| of the gap increases. This is reinforced by the tendency of off-diagonal matrix elements |A_jk| of local observables to decay as the levels considered are further apart.
Existing numerical results on off-diagonal matrix elements of local observables in non-integrable models are consistent with all the requirements listed here <cit.>.
The present arguments suggest distributions v_α that decay for larger values of |G_α|, and are therefore unimodal, and also smoothly distributed. This is confirmed in the simulations in Appendix <ref> in a non-integrable model on Fig. <ref> (left), and to a somewhat lesser extent in an integrable model too on Fig. <ref> (right).




§ SIMULATIONS


To calculate the function given in Definition <ref> in the main text

exactly one needs to find the maximum sum of p_α such that α:G_α∈ [G_λ, G_λ+x]. This calculation scales quite unfavourably with system size. If we have N energies, the number of intervals one must probe is quadratic in N. For each x ∼𝒪(10^-1) the intervals near the center are quite dense, making the entire algorithm for one choice of x approximately  scale like  𝒪(N^3). For this reason, we exactly diagonalize the Hamiltonian given in Eq. <ref>, 
and numerically approximate ξ(x). This is done by means of a Monte Carlo scheme where we randomly select intervals defined by G_λ using a normal distribution defined by μ_G = ∑_α p_α G_α and σ_G given in Definition <ref>. 


Figure <ref> (left) depicts the accuracy of this scheme. Unsurprisingly L=6 is exactly calculated and is not visible on the except at one location. The other cases show the approximation scheme performs better at larger system sizes. Despite this improvement, the accuracy of the scheme roughly puts us accurate to the fourth digit in all cases, making this scheme more than accurate enough. 
Quantities such as μ_G and σ_G can be calculated exactly given the exact diagonalization. However the left hand side of Eq. <ref>
in the main text has a time order complexity of 𝒪(N^4), making it again extremely difficult to calculate exactly. To get around this, we simply define a grid t_k = Δ k where k=0,1,2 … and average over the values calculated of |C^A(t_k)-C_∞^A|^2 as, 

    c(k)_T = 1/k+1∑_i=0^k |C^A(t_i)-C_∞^A|^2.









Figure <ref> (right)  shows the forward error of this scheme, showing an expected first order accuracy in time. Since we have an accuracy which is satisfactory for the scales we are comparing with the bound which tends have a roughly 10^-1 disagreement between the two sides of the bound. Finally to the optimal choice of a(ϵ) and δ(ϵ). 
For the plots present in Figure <ref> in the main text
we simply took the smallest δ(ϵ) available to minimize the resolution of our bound and picked the corresponding a(ϵ). This choice has an obvious issue in the L=8 case but begins to be more favourable in the larger system sizes. Looking at Figure <ref>
we see the value a(ϵ) can grow quite quickly due to finite size effects, making the prefactor outside the 1/T term quite large.




 §.§ Integrable models

Next, this section provide an example of showing how our bounds on timescales are affected in integrable models, highlighting the negative effect of degeneracies. Suppose we choose to define our Hamiltonian of Eq. (<ref>)
in the main text with parameters H = (-0.5,0, -0.5, 0 ). This corresponds to an Ising model with a transverse field. The issue in general with this model comes from investigating the behaviour of the corresponding δ(ϵ) and the fact that the frequencies G_α are very degenerate, meaning this function will not necessarily decay to zero as we take ϵ→ 0. 



In Figure <ref> we see the issue emerging with the bound found in Theorem <ref>. 
The degeneracy of the G_α terms cause the decrease in δ(ϵ) to happen in discrete steps triggered by calculating δ(ϵ) in a small enough region to differentiate two degenerate values of G_α which are close. Thus at small ϵ we still expect our resolution of equilibrium to be quite large. This slow decay of δ(ϵ) also causes a(ϵ) to become quite large very quickly, as one needs δ(ϵ) to be roughly linear for a(ϵ) to be reasonably small. This suggests that perhaps alternative approaches are required to bound the equilibration of two point time correlation functions in integrable models. 




 §.§ Distribution of v_α≡ρ_jj |A_jk|^2/C_0


Finally, we show the distributions of v_α≡ρ_jj |A_jk|^2/C_0 and comment on the differences between the integrable case  and the case that obeys the ETH. To proceed we define a coarse grained version p̅_α of v_α≡ρ_jj |A_jk|^2/C_0, where we define n bins and bins, 
b_1 = [G_min, G_min+ Δ G], b_2 = [ G_min + Δ G, G_min + 2Δ G ], …, where Δ G = G_max-G_min/n. The coarse grained probability is then obtained by summing the associated probabilities, p̅_α = ∑_G_β∈ b_β v_β. 




The result is given in Figure <ref>. The ETH case approaches a unimodal distribution quicker than the integrable case, however both distributions appear favorable in the coarse grained probabilities. Note, increasing the number of bins significantly did not significantly affect the shape of the curve. 


