Time evolution of correlation functions in quantum many-body systems

TL;DR

This paper rigorously analyzes the time evolution of two-point correlation functions in quantum many-body systems, showing how they factorize at late times and providing bounds on the timescales involved, with implications for understanding thermalization.

Contribution

It provides the first rigorous bounds on the timescales for correlation functions to reach equilibrium in quantum many-body systems, extending to various types of two-point functions.

Findings

Correlation functions factorize at late times in translation-invariant models.

Fluctuations around late-time values are bounded by the purity of the thermal state.

The timescale for auto-correlation functions to reach equilibrium is size-independent and bounded by local expectation values.

Abstract

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 $\langle A(t) B\rangle_\beta \rightarrow \langle A \rangle_\beta \langle B \rangle_\beta$, 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.

Equations

$$C^{AB}(t)$$

$$⟨A(t)B{⟩}_{\beta }\stackrel{t\to \infty }{⟶}⟨A{⟩}_{\beta }⟨B{⟩}_{\beta }.$$

$$\underset{X\in M,Y\in N}{\max }\frac{|\mathrm{Tr}\left(\rho X\otimes Y\right)-\mathrm{Tr}\left(\rho X\right)\mathrm{Tr}\left(\rho Y\right)|}{||X||||Y||}\le e^{-\frac{\text{dist}(M,N)}{\xi }},$$

$$\underset{T\to \infty }{\lim }{\int }_0^T{C}^{AB}(t)\frac{\text{d}t}{T}=$$

$$+\mathcal{O}\left({\xi }^{\frac{2D}{D+1}}{\log }^2(N)N^{-\frac{2}{D+1}}\right)$$

$${\sigma }_C^2=\underset{T\to \infty }{\lim }{\int }_0^T\frac{\text{d}t}{T}{\left(C^{AB}(t)-C_{\infty }^{AB}\right)}^2.$$

$$E_j-E_k=E_m-E_l\iff j=m,k=l,$$

$${\sigma }_C^2\le$$

$$\frac{C^A(t)}{C^A(0)}$$

$${\xi }_p(x)\equiv \underset{\lambda }{\max }\sum _{\alpha :G_{\alpha }\in \left[G_{\lambda },G_{\lambda }+x\right]}{p}_{\alpha }.$$

$$a(ϵ)\equiv \frac{{\xi }_p(ϵ)}{ϵ}{\sigma }_G,\delta (ϵ)\equiv {\xi }_p(ϵ),$$

$$\frac{1}{T}{\int }_0^T\frac{|C^A(t)-C_{\infty }^A{|}^2}{{\left(C^A(0)\right)}^2}dt$$

$${\sigma }_G^2$$

$$C_s^A(t)\equiv \frac{1}{2}\mathrm{Tr}\left(\rho \left\{A,A(t)\right\}\right)=\frac{C^A(t)+C^A(t{)}^{\ast }}{2}.$$

$$\frac{1}{T}{\int }_0^T\frac{|C_s^A(t)-C_{s,\infty }^A{|}^2}{(C_s^A(0){)}^2}dt$$

$${\sigma }_G^2$$

$$T_{eq}$$

$$\displaystyle T_{eq}=\frac{3\pi\,a(\epsilon)\,\sqrt{C^{A}_{s}(0)}}{\sqrt{\left|\frac{d^{2}C^{A}_{s}(t)}{dt^{2}}\big{|}_{0}\right|}}.$$

$$C_s^A(t)$$

$$C_{\text{Kubo}}(t)\propto \sum _{j\ne k}\frac{e^{-\beta E_k}-e^{-\beta E_j}}{E_j-E_k}|A_{jk}{|}^2{e}^{it(E_j-E_k)}.$$

$$\frac{1}{T}{\int }_0^T\frac{|C_{\text{Kubo}}(t)-C_{\text{Kubo},\infty }{|}^2}{C_{\text{Kubo}}(0{)}^2}dt$$

$${\sigma }_G^2$$

$$H=\sum _{j=1}^L\left(\gamma {\sigma }_j^X+\lambda {\sigma }_j^Z\right)+J\sum _{j=1}^{L-1}{\sigma }_j^Z{\sigma }_{j+1}^Z+\alpha \sum _{j=1}^{L-2}{\sigma }_j^Z{\sigma }_{j+2}^Z,$$

$${\text{Pr}}_{|E_k⟩\in \rho }(|⟨E_k|A|E_k⟩-\mathrm{Tr}\left(\rho A\right)|\ge \delta )\le \exp (-c\delta N^{\frac{1}{D+1}}{\xi }^{-\frac{D}{D+1}}),$$

$${\text{Pr}}_{|E_k⟩\in \rho }(⟨E_k|A|E_k⟩-\mathrm{Tr}\left(\rho A\right)\ge \delta )$$

$$\le e^{-\lambda \delta N}\sum _k{\rho }_{kk}{e}^{\lambda ⟨E_k|{\bar{A}}^{\prime }|E_k⟩}$$

$$\le e^{-\lambda \delta N}\mathrm{Tr}\left(\rho e^{\lambda {\bar{A}}^{\prime }}\right).$$

$$\mathrm{Tr}\left(\rho e^{\lambda {\bar{A}}^{\prime }}\right)=\sum _{a_l\le \delta N/2}{e}^{\lambda a_l}\mathrm{Tr}\left(\rho {\Pi }_l\right)+\sum _{a_l>\delta N/2}{e}^{\lambda a_l}\mathrm{Tr}\left(\rho {\Pi }_l\right).$$

$$\sum _{a_l>\delta N/2}{e}^{\lambda a_l}\mathrm{Tr}\left(\rho {\Pi }_l\right)=\sum _j\sum _{\delta N/2+(j+1)\ge a_l\ge \delta N/2+j}{e}^{\lambda a_l}\mathrm{Tr}\left(\rho {\Pi }_l\right)\le \sum _j{e}^{\lambda (\delta N/2+j+1)}\mathrm{Tr}\left(\rho {\Pi }_{\ge \delta N/2+j}\right),$$

$$\mathrm{Tr}\left(\rho {\Pi }_{\ge \delta N/2+j}\right)\le e^{Dk{N}^{\alpha }/\xi }\text{exp}\left(-\frac{c^{\prime }}{D\xi }{\left(N(\delta /2+j/N{)}^2\xi \right)}^{1/D+1}\right).$$

$$\sum _{a_l>\delta N/2}{e}^{\lambda a_l}\mathrm{Tr}\left(\rho {\Pi }_l\right)\le \text{poly}(N)e^{Dk{N}^{\alpha }/\xi }\text{exp}\left(-c_1{N}^{\frac{1}{D+1}}{\xi }^{-\frac{D}{D+1}}\right).$$

$${\text{Pr}}_{|E_k⟩\in \rho }(⟨E_k|A|E_k⟩-\mathrm{Tr}\left(\rho A\right)\ge \delta )$$

Full text

Two-point correlation functions –or dynamical response/Green’s functions– are the central object of the theory of linear response Kubo (1957), 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 Rickayzen (1980). This includes different types of scattering and spectroscopy experiments Zhu (2005), quantum transport Zwanzig (1965); Zotos et al. (1997), and fluctuation-dissipation relations Srednicki (1999); Khatami et al. (2013); D’Alessio et al. (2016). They have also appeared in the characterization of topological Nussinov and Ortiz (2008) and crystalline ordering Watanabe and Oshikawa (2015), of quantum-chaotic systems Luitz and Lev (2016); Gharibyan et al. (2019) and of different notions of ergodicity in quantum and classical systems Cornfeld et al. (2012); Venuti and Liu (2019).

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

[TABLE]

where the evolution is generated by a time-independent Hamiltonian $H$, $\rho\equiv e^{-\beta H}/Z_{\beta}$ is a thermal state at inverse temperature $\beta$ with partition function $Z_{\beta}$, and $A(t)=e^{iHt}Ae^{-iHt}$ 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 Steinigeweg and Gemmer (2009), QMC Starykh et al. (1997) and tensor networks Sirker and Klümper (2005); Barthel et al. (2009); Feiguin and Fiete (2010); Karrasch et al. (2012); Dargel et al. (2012); Barthel (2013); Bohrdt et al. (2017), and analytically for specific models, e.g. Korepin et al. (1997); Its et al. (1993); Stolze et al. (1995); Reyes and Tsvelik (2006); Perk and Au-Yang (2009); Bertini et al. (2019). Also, a number of experimental schemes to measure it directly have been proposed Romero-Isart et al. (2012); Knap et al. (2013); Buscemi et al. (2013); Gessner et al. (2014); Uhrich et al. (2017), 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

[TABLE]

Moreover, we show that the fluctuations around the late-time value are in fact bounded by the effective dimension $d^{-1}_{\text{eff}}\equiv\operatorname{\textnormal{Tr}}\left({\rho^{2}}\right)$ 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. (2) 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)\equiv\langle A(t)A\rangle_{\beta}+\langle AA(t)\rangle_{\beta}$, we find that evolution is dominated by a timescale which is at most of order $t^{2}\sim\frac{\langle A^{2}\rangle_{\beta}}{\langle[A,H][H,A]\rangle_{\beta}}$. 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 Kubo (1957); Khatami et al. (2013), 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.

[TABLE]

where $M,N$ are regions on the lattice separated by a distance of at least $\textnormal{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 Araki (1969) and arbitrary models above a threshold temperature Kliesch et al. (2014). 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=\sum_{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\rightarrow\infty}C^{AB}(t)$ are not well-defined. Hence, we consider the late-time behaviour under infinite-time averages of the correlation functions $\lim_{T\rightarrow\infty}\int_{0}^{T}\frac{\text{d}t}{T}C^{AB}(t)$. With these considerations, our first main result is the following.

[TABLE]

This guarantees that all operators supported on a region with size scaling like any function smaller than $\mathcal{O}(N^{\frac{1}{D+1}})$, satisfy the assumptions of the theorem. The proof, found in Appendix A.1, relies on a weak form of the Eigenstate Thermalization Hypothesis (ETH) shown in Brandão et al. (2019), which is itself based on previous works on large deviation theory for lattice models Anshu (2016); Mori (2016). 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 Biroli et al. (2010).

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 D’Alessio et al. (2016).

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$ Kubo (1957).

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 D’Alessio et al. (2016); Gogolin and Eisert (2016), where it is generally expected to hold as random perturbations are sufficient to lift degeneracies in energy gaps Cohen-Tannoudji et al. (2006).

Let us define $C^{AB}_{\infty}=\lim_{T\rightarrow\infty}\int_{0}^{T}\frac{\text{d}t}{T}C^{AB}(t)$, and the average fluctuations around the late-time value as

[TABLE]

The following result puts an upper bound on average fluctuations.

[TABLE]

*and let $[\rho,H]=0$. It holds that *

[TABLE]

where $A_{kj},B_{jk}$ are matrix elements in the energy eigenbasis, $A=\sum_{jk}A_{jk}\ket{j}\!\bra{k}$.

The proof can be found in Appendix A.2. It follows the same steps as the main result in Short (2011). Here, we also find that the purity $\operatorname{\textnormal{Tr}}\left({\rho^{2}}\right)$ of the equilibrium ensemble plays a key role. For a microcanonical ensemble $\operatorname{\textnormal{Tr}}\left({(\mathbbm{1}/d)^{2}}\right)=1/d$, so the RHS of Eq. (44) is expected to decay exponentially with system size in most situations of interest. Also, notice that for a thermal state $\operatorname{\textnormal{Tr}}\left({\rho_{\beta}^{2}}\right)\leq 1/Z_{\beta}$. Moreover, the ETH predicts that $|A_{kj}B_{jk}|\sim 1/d$ Deutsch (2018).

Timescales of equilibration — Theorems 1 and 2 combined imply that correlation functions of the form $\langle A(t)B\rangle_{\beta}$ are, for most times $t\in[0,\infty]$, close to the uncorrelated average $\langle A\rangle_{\beta}\langle B\rangle_{\beta}$, 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 Lieb and Robinson (1972); Huang et al. (2018), and timescales associated with ballistic ($\propto N^{1/D}$) or diffusive ($\propto N^{2/D}$) processes may play a role. However, for the autocorrelation function $C^{A}(t)\equiv\langle A(t)A\rangle_{\beta}$, 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 Dymarsky (2018); Schiulaz et al. (2019), and for those effects our result limits their relative size.

Let us define $\rho=\sum_{j}\rho_{jj}\ket{j}\!\bra{j}$, so that $\rho_{jj}$ and $A_{jk}$ are the matrix elements of $\rho$ and $A$ in the energy basis. We can then write

[TABLE]

where we denote each pair of levels $\{i,j\}$ by a Greek index, and the corresponding energy gaps by $G_{\alpha}\equiv E_{j}-E_{k}$ (notice that both $E_{j}-E_{k}$ and $E_{k}-E_{j}$ appear in the sum). The normalized distribution $v_{\alpha}\equiv\frac{\rho_{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.

[TABLE]

We also define

[TABLE]

where $\sigma_{G}=\sqrt{\sum_{\alpha}p_{\alpha}G_{\alpha}^{2}-\left(\sum_{\alpha}p_{\alpha}G_{\alpha}\right)^{2}}$ is the standard deviation of the distribution $p_{\alpha}$ over the energy gaps $G_{\alpha}$.

The important point behind these definitions is that, for a sufficiently smooth and unimodal probability distribution, one can find an $\epsilon$ small enough such that $a(\epsilon)\sim\mathcal{O}(1)$ and $\delta(\epsilon)\ll 1$. In the following theorem, the relevant probability distribution is given by $v_{\alpha}$. Our main result regarding the timescales of correlation functions, proven in Appendix B.2, is:

[TABLE]

for any $\epsilon>0$. Here, $a(\epsilon)$ and $\delta(\epsilon)$ are as in Definition 2 for the normalized distribution $v_{\alpha}\equiv\frac{\rho_{jj}|A_{jk}|^{2}}{C^{A}(0)}$, and $\sigma_{G}$ is given by

[TABLE]

Theorem [3] provides an upper bound of $T_{eq}\equiv\frac{3\pi\,a(\epsilon)}{\sigma_{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. (11) 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}_{\infty}$ .

The crucial point is that for distributions $v_{\alpha}$ that are uniformly spread over many values of the gaps $G_{\alpha}$, one can always find an $\epsilon$ such that $\delta\ll 1$. In that case, the right hand side of Eq. (11) becomes small on timescales $\mathcal{O}(T_{eq})$. As discussed in (García-Pintos et al., 2017) and in Appendix B.5, if one further assumes smooth unimodal distributions, typically one also finds that $a\sim\mathcal{O}(1)$. In that case, the timescale is governed by $1/\sigma_{G}$. Given that $\sigma_{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 Kim et al. (2015) shows that a timescale of order $1/\sigma_{G}$ provides a lower bound to the timescale of equilibration, which strongly suggests that our upper bound is tight when the conditions of $a\sim\mathcal{O}(1)$ and $\delta\ll 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 Beugeling et al. (2015); Mondaini and Rigol (2017), one should be able to choose $\epsilon$ such that $a\sim\mathcal{O}(1)$ and $\delta\ll 1$. In Fig. 1 we numerically show that this is indeed the case in a non-integrable Ising model.

Theorem 3 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(\epsilon)$ and $\delta(\epsilon)$ will be large no matter what value of $\epsilon$ is chosen, in which case the RHS of Eq.  (11) may not become small within reasonable timescales. This can happen, for instance, in models with highly degenerate energy spectrum. To illustrate this, in Appendix C we compute these parameters in an integrable model, where we see that the gap degeneracies of the model negatively affect the quantities $a(\epsilon)$ and $\delta(\epsilon)$, making the estimated equilibration timescales longer.

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

[TABLE]

Along the same lines of Theorem [3], in Appendix B.3 we prove the following.

[TABLE]

for any $\epsilon>0$. Here, $a(\epsilon)$ and $\delta(\epsilon)$ are as in Definition  [2] for the normalized distribution $v_{\alpha}^{s}\equiv\frac{\rho_{jj}+\rho_{kk}}{2}\frac{|A_{jk}|^{2}}{C^{A}_{s}(0)}$, and

[TABLE]

Thus an upper bound for the equilibration timescale is

[TABLE]

where again we expect that for small enough $\epsilon$, $a(\epsilon)\sim\mathcal{O}(1)$ and $\delta\ll 1$ for the same reasons as before. The denominator in $T_{eq}$ can be seen as an “acceleration” of the symmetric autocorrelation function. Eq. (16) can in fact be written as

[TABLE]

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

[TABLE]

For early times, the above expression decays on a timescale $\tau=\frac{\sqrt{2}}{3\pi a(\epsilon)}T_{eq}$, identical to our upper bound Eq. (17) up to a prefactor.

The timescale of Eq. (16) 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 $\langle A(t)A\rangle_{\beta}\longrightarrow\langle A\rangle_{\beta}\langle A\rangle_{\beta}$ occurs, and ii) $C^{A}_{s}(0)=\operatorname{\textnormal{Tr}}\left({\rho A^{2}}\right)$ as a measure of the fluctuations of $A$. Then, Eq. (16) gives a proportionality relation between the strength of the fluctuations and the timescale of equilibration, in a similar spirit to what was found in Nation and Porras (2019) 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 $\rho\propto e^{-\beta(H+\lambda A)}$. Subsequently, the Hamiltonian is slightly perturbed by $\lambda A$, so that the evolved state is $\rho_{t}=e^{-itH}\rho e^{itH}$.

It was shown by Kubo Kubo (1957) that, to leading order in $\lambda$, the expectation value of $A$ satisfies $\operatorname{\textnormal{Tr}}\left({\rho A(t)}\right)~{}=~{}C_{\text{Kubo}}(t)\operatorname{\textnormal{Tr}}\left({\rho A}\right)$, where for thermal initial states $\rho$ the Kubo correlation function can be written as

[TABLE]

Equilibration of $\operatorname{\textnormal{Tr}}\left({\rho A(t)}\right)$ is then equivalent to equilibration of the function $\operatorname{C_{\textnormal{Kubo}}}(t)$, for which we prove in Appendix B.4 that the following holds.

[TABLE]

for any $\epsilon>0$. Here, $a(\epsilon)$ and $\delta(\epsilon)$ are as in Definition  [2] for the normalized distribution $w_{\alpha}\equiv\frac{e^{-\beta E_{k}}-e^{-\beta E_{j}}}{E_{j}-E_{k}}\frac{|A_{jk}|^{2}}{C_{\textnormal{Kubo}}(0)}$, and

[TABLE]

This again implies an upper bound $T_{eq}=\frac{3\pi\,a(\epsilon)}{\sigma_{G}}$ on the equilibration timescale of $\operatorname{C_{\textnormal{Kubo}}}$, and therefore on the time to return to thermal equilibrium after a perturbation of the system Hamiltonian by $A$. The distribution $w_{\alpha}$ plays the same role as $v_{\alpha}$ and $v^{s}_{\alpha}$ before. If $w_{\alpha}$ is smoothly distributed and unimodal (which we expect for local observables in non-integrable models) then $a\sim\mathcal{O}(1)$ and $\delta(\epsilon)\ll 1$ holds (see Appendix B.5).

Simulations — We test Theorem [3] in a spin model governed by the Hamiltonian

[TABLE]

where $\sigma_{j}^{Z}$ and $\sigma_{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 $(\gamma,\lambda,J,\alpha)$ characterize the model. We focus on a case corresponding to a system satisfying ETH by choosing $(\gamma,\lambda,J,\alpha)=(0.8,0.5,1,1)$ Kaneko et al. (2019), and study the autocorrelation functions of the observable $A=\sigma_{\frac{L}{2}}^{x}.$ For simplicity we set $\beta=1$ in our numerics, though no significant changes were observed for $\beta\in[0.1,5]$. Figure 1 depicts the functions $a(\epsilon)$ and $\delta(\epsilon)$ that appear in Theorem 3, confirming that there exist regions of $\epsilon$ such that $\delta\ll 1$, ensuring equilibration occurs, and $a\sim 0.4$. Importantly, this is increasingly the case as the size of the system grows.

Figure 2 compares the two sides in bound (11), 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 3 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 Kim et al. (2015). Details of the simulations can be found in C.

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 Short and Farrelly (2012); Malabarba et al. (2014); García-Pintos et al. (2017), for which finding rigorous estimates on the timescales is a largely open problem Wilming et al. (2017); de Oliveira et al. (2018); Dymarsky (2019); Reimann (2019). This connection is not surprising, specially considering that previous works Srednicki (1999); Richter et al. (2019); Richter and Steinigeweg (2019) 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 Bloch et al. (2008); Schneider et al. (2012).

Appendix A Late time behaviour of two-point functions

A.1 Proof of late-time equilibration

First, we state the key result from Brandão et al. (2019) that we use. It says that expectation values for single eigenstates, of the form $\bra{E_{k}}A\ket{E_{k}}$, are close to the ensemble average $\langle A\rangle_{\beta}$ 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 Brandão et al. (2019) (which itself builds on Mori (2016)), with the difference that our version holds for lattices of dimension larger than $1$.

Lemma 1** (Proposition 7 Brandão et al. (2019)).**

Let $H$ be a translation-invariant, non-degenerate Hamiltonian with $N$ sites on a $D$-dimensional lattice, $\rho$ an equilibrium ensemble $[\rho,H]=0$ with finite correlation length $\xi$, and $A$ some observable with support on a connected region of at most $\mathcal{O}(N^{1/(D+1)-\nu})$ sites, with $\nu>0$. Then, for any $\delta>0$,

[TABLE]

where $c>0$ is a constant, and $\ket{E_{k}}\in\rho$ indicates that the eigenstates are sampled from the equilibrium distribution $\rho$.

Proof.

We show a bound $\text{Pr}_{\ket{E_{k}}\in\rho}(\bra{E_{k}}A\ket{E_{k}}-\operatorname{\textnormal{Tr}}\left({\rho A}\right)\geq\delta)$, and $\text{Pr}_{\ket{E_{k}}\in\rho}(\operatorname{\textnormal{Tr}}\left({\rho A}\right)-\bra{E_{k}}A\ket{E_{k}}\geq\delta)$ will follow analogous steps. Notice that since the Hamiltonian is translation-invariant and non-degenerate we can write $\bra{E_{k}}A\ket{E_{k}}=\bra{E_{k}}\bar{A}/N\ket{E_{k}}$, where $\bar{A}=\sum_{i}A_{i}$ is the extensive observable built out of translations of $A$. Define $\bar{A}^{\prime}=\bar{A}-\langle\bar{A}\rangle_{\beta}$. Then, using Markov’s inequality and $e^{\bra{\psi}A\ket{\psi}}\leq\bra{\psi}e^{A}\ket{\psi}$ we can write

[TABLE]

Now let us write the spectral decomposition $e^{\lambda\bar{A}^{\prime}}=\sum_{l}e^{\lambda a_{l}}\Pi_{l}$, with $\Pi_{l}$ the projectors into the subspace with eigenvalue $e^{\lambda a_{l}}$, and write the average as

[TABLE]

The first term is upper bounded by $e^{\lambda\delta N/2}$. For the second, we write

[TABLE]

where $\Pi_{\geq x}$ denotes the projector on the subspace with $a_{l}\geq x$. The main result of Anshu (2016) states that

[TABLE]

Here, $kN^{\alpha}$ is the number of sites on which an individual local term of observable $\bar{A}$ has support. Since $j$ is at most $\mathcal{O}(N)$, we can choose some $\lambda=\mathcal{O}\left((N\xi)^{-\frac{D}{D+1}}\right)$ such that, for some constant $c_{1}>0$,

[TABLE]

This way, the dominant contribution of Eq. (27) is the first term. Plugging the bounds back in Eq. (26) results in the following, for some constant $c>0$ and large enough $N$,

[TABLE]

where the last line follows from the fact that the second term in the RHS of Eq.(31) is subleading (much smaller than $1$) as long as $\alpha<\frac{1}{D+1}$. This sets the support of of the observable $A$ ot be within a region of size at most $\mathcal{O}(N^{1/(D+1)-\nu})$ for any $\nu>0$. ∎

With it, we are now ready to prove the result on late-time factorization of correlation functions.

Theorem 1**.**

Let $H$ be a $k$-local, translation-invariant and non-degenerate Hamiltonian on a $D$-dimensional Euclidean lattice of $N$ sites, and let $[\rho,H]=0$ be an equilibrium ensemble (such as a thermal state) of finite correlation length $\xi>0$. Let $A,B$ be local observables with support on at most $\mathcal{O}(N^{\frac{1}{D+1}-\nu})$ sites, with $\nu>0$. Then

[TABLE]

Proof.

First let us write

[TABLE]

from which we have

[TABLE]

Since the Hamiltonian is non-degenerate by assumption, it holds that

[TABLE]

so that the limit becomes

[TABLE]

Now let us define $A_{kk}-\operatorname{\textnormal{Tr}}\left({\rho A}\right)=\Delta_{k,A}$ and $B_{kk}-\operatorname{\textnormal{Tr}}\left({\rho B}\right)=\Delta_{k,B}$ , so that Huang (2019)

[TABLE]

Let us define $\Delta\equiv K\log N/N^{\frac{1}{D+1}}$, and split the sum over energies of the error term as

[TABLE]

where $\mathcal{S}\equiv\{k:|\Delta_{k,A}|,|\Delta_{k,B}|\leq\Delta\}$ (that is, the set of $k$ for which both errors are small). Notice that the first term is smaller than $\Delta^{2}$ by definition. On the other hand, the second term can be bounded as

[TABLE]

The third line follows from Lemma 1, and the fourth from $|\Delta_{k,A}|\leq||A||,|\Delta_{k,B}|\leq||B||$. The constant $K$ is arbitrary, so we can choose it such that $cK\xi^{-\frac{D}{D+1}}=\frac{2}{D+1}$. In that case the dominant contribution to Eq. (40) is that of the first term, and hence $\sum_{k}\rho_{kk}\Delta_{k,A}\Delta_{k,B}=\mathcal{O}(\Delta^{2})$, so that

[TABLE]

completing the proof. ∎

A.2 Proof of fluctuations around late-time value

Theorem 2**.**

Let $H=\sum_{i}E_{i}\ket{i}\bra{i}$ be a Hamiltonian with non-degenerate energy gaps, such that

[TABLE]

*and let $[\rho,H]=0$. It holds that *

[TABLE]

where $A_{kj},B_{jk}$ are off-diagonal matrix elements in the energy eigenbasis, $A=\sum_{jk}A_{jk}\ket{j}\!\bra{k}$.

Proof.

Let us expand in the energy eigenbasis.

[TABLE]

In the second to the third line we used the assumption of non-degenerate energy gaps. We now use the Cauchy-Schwarz inequality twice.

[TABLE]

The last inequality follows from the fact that for positive operators $\operatorname{\textnormal{Tr}}\left({PQ}\right)\leq||P||\operatorname{\textnormal{Tr}}\left({Q}\right)$.

∎

Appendix B Dynamics of two-time correlation functions

B.1 Preliminaries

Here we prove an intermediate lemma required for Appendix B.2. This is similar to results found in  Malabarba et al. (2014); García-Pintos et al. (2017), improving on them by a numerical factor. Let us first recall the following definition from the main text.

Definition 2**.**

Given a normalized distribution $p_{\alpha}$ over $G_{\alpha}$, we define $\xi_{p}(x)$ as the maximum weight that fits an interval of energy gaps with width $x$:

[TABLE]

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

Lemma 2**.**

Let $f(t)$ be a positive function of the form

[TABLE]

such that the $p_{\alpha}$ form a discrete probability distribution $\sum_{\alpha}p_{\alpha}=1$. The uniform average of such a function is upper bounded by,

[TABLE]

where $\xi_{p}(x)$ is as in Definition 2.

Proof.

Let us first define the uniform and Gaussian probability density functions, as

[TABLE]

and

[TABLE]

respectively, where the mean of the Gaussian is written as $\mu=\frac{T}{2}$ and the standard deviation as $\sigma=\alpha T$. The parameter $\alpha>0$ is free for us to choose. We will first show that for any positive function $f(t)\geq 0$, the uniform average of a function on the interval $t\in[0,T]$ can be bounded tightly by the Gaussian average by,

[TABLE]

where $\gamma=\sqrt{2\pi}\alpha e^{\frac{1}{8\alpha^{2}}}$, and we denote the uniform and Gaussian averages over distributions $p_{T}(t)$ and $p_{G}(t)$ as $\langle f(t)\rangle_{T}$ and $\langle f(t)\rangle_{G_{T}}$ respectively.

To do this, note that due to the choice of $\mu=\frac{T}{2}$ the two distributions have identical means. To find the smallest possible $\gamma$ such that (60) holds we set that $p_{T}(t=T)=\gamma p_{G}(t=T)$. This then implies,

[TABLE]

which guarantees that $p_{T}(t)\leq\gamma p_{G}(t)\,\,\forall t$. Since over the interval $[0,T]$ $p_{T}(t)$ is constant and $p_{G}(t)\geq p_{G}(T)$,it follows that

[TABLE]

With this we can now bound the uniform average of $f(t)$ as

[TABLE]

Interchanging the integral with the sum and then integrating gives the characteristic function of a Gaussian distribution. This allows us to write,

[TABLE]

where $\Delta G=G_{\alpha}-G_{\beta}$. To further bound Eq. (64) we may introduce the function,

[TABLE]

This function satisfies

[TABLE]

where $r>1$. Applying this bound to Eq. (64) with $r=e^{\frac{\alpha}{2}}$ and $x=\Delta G^{2}T^{2}$, it follows that

[TABLE]

The presence of $g(x)$ allows to restrict the summation in Eq. (67). With the condition $\left(\Delta G^{2}T^{2}-n\right)\in[0,1)$, the two following intervals can be identified for the variable $G_{\beta}$ given a particular value of $G_{\alpha}$,

[TABLE]

With this we rewrite the summation in the bound as,

[TABLE]

Given that $\sqrt{n}\leq n$ for all integer $n$, one can expand the length of the interval, so that

[TABLE]

where we define

[TABLE]

It is straightforward to see that $\sum_{G_{\beta}\in I^{\prime}_{+}}p_{\beta}+\sum_{G_{\beta}\in I^{\prime}_{-}}p_{\beta}\leq 2\xi_{p}\left(\frac{1}{T}\right)$, which allows us to further bound the summation by using Def. 2,

[TABLE]

To complete the proof we find the $\alpha$ that minimizes

[TABLE]

which occurs for $\alpha\approx 0.6347$. This gives $\kappa\approx 2.8637<3$. Finally,

[TABLE]

which completes the proof.

∎

B.2 Proof of equilibration timescales of correlation functions

Theorem 3**.**

For any Hamiltonian $H$ and state $\rho$ such that $[H,\rho]=0$, and any observable $A$, the time correlation function $C^{A}(t)=\operatorname{\textnormal{Tr}}\left({\rho A(t)A}\right)$ satisfies

[TABLE]

for any $\epsilon>0$. Here, $a(\epsilon)$ and $\delta(\epsilon)$ are as in Definition 2 for the normalized distribution $v_{\alpha}\equiv\frac{\rho_{jj}|A_{jk}|^{2}}{C^{A}(0)}$, and $\sigma_{G}$ is given by

[TABLE]

Proof.

Following Lemma 2, we have that

[TABLE]

where $\xi_{v}\left(\tfrac{1}{T}\right)$ is taken with respect to the distribution $v_{\alpha}$, which is normalized:

[TABLE]

It was shown in García-Pintos et al. (2017) (Proposition 5) that the function $\xi_{p}(x)$ satisfies

[TABLE]

with $a(\epsilon)=\frac{\xi_{p}(\epsilon)}{\epsilon}\sigma_{G}$ and $\delta(\epsilon)=\xi_{p}(\epsilon)$, and where $\varepsilon\in(0,\infty)$. Finally, the standard deviation of the distribution $p_{v}$ is

[TABLE]

which completes the proof. ∎

B.3 Proof of equilibration timescales of symmetric correlation functions

Theorem 4**.**

For any Hamiltonian $H$ and state $\rho$ such that $[H,\rho]=0$, and any observable $A$, the time correlation function $C_{s}^{A}(t)=\operatorname{\textnormal{Tr}}\left({\rho\{A,A(t)\}}\right)$ satisfies

[TABLE]

for any $\epsilon>0$. Here, $a(\epsilon)$ and $\delta(\epsilon)$ are as in Definition 2 for the normalized distribution $v_{\alpha}^{s}\equiv\frac{\rho_{jj}+\rho_{kk}}{2}\frac{|A_{jk}|^{2}}{C^{A}_{s}(0)}$, and

[TABLE]

Proof.

The symmetric correlation function is defined as

[TABLE]

The equivalent of Eq.8 becomes

[TABLE]

The proof of Theorem 4 is identical as the previous proof, with the symmetrized distribution $v_{\alpha}^{S}\equiv\frac{\rho_{jj}+\rho_{kk}}{2}|A_{jk}|^{2}$. In this case the variance of the normalized distribution $v_{\alpha}^{S}$ becomes

[TABLE]

∎

Equilibration then occurs within a timescale

[TABLE]

The denominator in $T_{eq}$ can be identified as an “acceleration” of the symmetric autocorrelation function. Indeed,

[TABLE]

Then, the equilibration timescale is

[TABLE]

B.3.1 Short-time evolution of symmetric correlation functions

The symmetric autocorrelation functions is given by

[TABLE]

Taking the Taylor expansion of $C^{A}_{s}(t)$,

[TABLE]

For early times, the above expression decays on a timescale

[TABLE]

B.4 Proof of equilibration timescales of Kubo correlation functions

Theorem 5**.**

For any Hamiltonian $H$, thermal state $\rho$, and any observable $A$, the Kubo correlation function $C_{\textnormal{Kubo}}$ satisfies

[TABLE]

for any $\epsilon>0$. Here, $a(\epsilon)$ and $\delta(\epsilon)$ are as in Definition 2 for the normalized distribution $w_{\alpha}\equiv\frac{e^{-\beta E_{k}}-e^{-\beta E_{j}}}{E_{j}-E_{k}}\frac{|A_{jk}|^{2}}{C_{\textnormal{Kubo}}(0)}$, and

[TABLE]

Proof.

The Kubo correlation function can be written as

[TABLE]

with the proportionality constant defined by $C_{\text{Kubo}}(0)=1$. We can then write

[TABLE]

where we define $w_{\alpha}\equiv\frac{e^{-\beta E_{k}}-e^{-\beta E_{j}}}{E_{j}-E_{k}}|A_{jk}|^{2}\geq 0$. Given that $w_{\alpha}\geq 0$, we can perform similar calculations as for Theorem 3, albeit with a different probability distribution. Thus, we also have

[TABLE]

Defining the normalized distribution $q_{\alpha}\equiv w_{\alpha}/\operatorname{C_{\textnormal{Kubo}}}(0)$. The variance is

[TABLE]

completing the proof. ∎

B.5 Scaling of $a$ and $\delta$

The proofs of Theorems (3-5) rely on the fact that the function $\xi_{p}(x)$, defined for any normalized distribution $p_{\alpha}$ as the maximum distribution that fits an interval $x$

[TABLE]

satisfies

[TABLE]

which was shown in García-Pintos et al. (2017) (Proposition 5). Here $a(\epsilon)\equiv\frac{\xi_{p}(\epsilon)}{\epsilon}\sigma_{G}$ and $\delta(\epsilon)\equiv\xi_{p}(\epsilon)$, where $\sigma_{G}$ is the standard deviation of the distribution $p_{\alpha}$. The function $a(\epsilon)$ ends up in the bound of the equilibration timescales, as $T_{eq}=\frac{\pi\,a(\epsilon)}{\sigma_{G}}$, while $\delta(\epsilon)$ governs the long time behavior in the bounds.

Given that $\xi_{p}(x)$ characterizes how much of the distribution $p_{\alpha}$ fits an interval $x$, the value of $a$ in Eq. (101) depends on how well $1/\sigma_{G}$ serves to characterize the region where the distribution $p_{\alpha}$ is supported. Roughly speaking, whenever $1/\sigma_{G}$ is a good estimate of the width of such small region, then one expects $a\sim\mathcal{O}(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/\sigma_{G}$ of the window where the distribution is supported, and $\xi_{p}(x)\sim x/\sigma_{G}$, so that $a\sim\mathcal{O}(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, $\delta$ in Eq. (101) carries information of the fine structure of $p_{\alpha}$, indicating the scale at which the distribution can no longer be coarse-grained to a continuous distribution. The only way that $\delta\ll 1$ fails is for distributions that are not smooth, in which a small region of width $\epsilon$ 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 $\epsilon$ such that $\delta(\epsilon)\ll 1$ and $a(\epsilon)\sim\mathcal{O}(1)$.

In summary, the problem of proving fast equilibration timescales in our approach can thus be linked to knowing whether the relevant distribution $p_{\alpha}$ 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 3, where the relevant distribution is given by $v_{\alpha}\equiv\frac{\rho_{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_{\alpha}|$ 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 Beugeling et al. (2015); Mondaini and Rigol (2017). The present arguments suggest distributions $v_{\alpha}$ that decay for larger values of $|G_{\alpha}|$, and are therefore unimodal, and also smoothly distributed. This is confirmed in the simulations in Appendix C in a non-integrable model on Fig. 5 (left), and to a somewhat lesser extent in an integrable model too on Fig. 5 (right).

Appendix C Simulations

To calculate the function given in Definition 2 in the main text exactly one needs to find the maximum sum of $p_{\alpha}$ such that $\alpha:G_{\alpha}\in[G_{\lambda},G_{\lambda}+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\sim\mathcal{O}(10^{-1})$ the intervals near the center are quite dense, making the entire algorithm for one choice of $x$ approximately scale like $\mathcal{O}(N^{3})$. For this reason, we exactly diagonalize the Hamiltonian given in Eq. 22, and numerically approximate $\xi(x)$. This is done by means of a Monte Carlo scheme where we randomly select intervals defined by $G_{\lambda}$ using a normal distribution defined by $\mu_{G}=\sum_{\alpha}p_{\alpha}G_{\alpha}$ and $\sigma_{G}$ given in Definition 2.

Figure 3 (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 $\mu_{G}$ and $\sigma_{G}$ can be calculated exactly given the exact diagonalization. However the left hand side of Eq. 11 in the main text has a time order complexity of $\mathcal{O}(N^{4})$, making it again extremely difficult to calculate exactly. To get around this, we simply define a grid $t_{k}=\Delta k$ where $k=0,1,2\dots$ and average over the values calculated of $|C^{A}(t_{k})-C_{\infty}^{A}|^{2}$ as,

[TABLE]

Figure 3 (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(\epsilon)$ and $\delta(\epsilon)$. For the plots present in Figure 2 in the main text we simply took the smallest $\delta(\epsilon)$ available to minimize the resolution of our bound and picked the corresponding $a(\epsilon)$. 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 1 we see the value $a(\epsilon)$ can grow quite quickly due to finite size effects, making the prefactor outside the $1/T$ term quite large.

C.1 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. (22) 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 $\delta(\epsilon)$ and the fact that the frequencies $G_{\alpha}$ are very degenerate, meaning this function will not necessarily decay to zero as we take $\epsilon\to 0$.

In Figure 4 we see the issue emerging with the bound found in Theorem 3. The degeneracy of the $G_{\alpha}$ terms cause the decrease in $\delta(\epsilon)$ to happen in discrete steps triggered by calculating $\delta(\epsilon)$ in a small enough region to differentiate two degenerate values of $G_{\alpha}$ which are close. Thus at small $\epsilon$ we still expect our resolution of equilibrium to be quite large. This slow decay of $\delta(\epsilon)$ also causes $a(\epsilon)$ to become quite large very quickly, as one needs $\delta(\epsilon)$ to be roughly linear for $a(\epsilon)$ 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.

C.2 Distribution of $v_{\alpha}\equiv\frac{\rho_{jj}|A_{jk}|^{2}}{C_{0}}$

Finally, we show the distributions of $v_{\alpha}\equiv\frac{\rho_{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 $\bar{p}_{\alpha}$ of $v_{\alpha}\equiv\frac{\rho_{jj}|A_{jk}|^{2}}{C_{0}}$, where we define $n$ bins and bins, $b_{1}=[G_{min},G_{min}+\Delta G],b_{2}=[G_{min}+\Delta G,G_{min}+2\Delta G],\dots$, where $\Delta G=\frac{G_{max}-G_{min}}{n}$. The coarse grained probability is then obtained by summing the associated probabilities, $\bar{p}_{\alpha}=\sum_{G_{\beta}\in b_{\beta}}v_{\beta}$.

The result is given in Figure 5. 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.