Equilibration time scales in closed many-body quantum systems

There is currently a renewed interest in the derivation of statistical mechanics from the kinematics and dynamics of a closed quantum system [1]. In this approach, instead of assuming a priori that the system is in some mixed state, such as e.g. a micro-canonical ensemble, one describes it at all times using a pure state $| \psi (t)\rangle$. One then seeks to show that, under reasonable conditions, the system behaves as if it were described by a statistical ensemble. In this way the use of statistical mechanics can be justified without introducing additional external degrees of freedom, such as e.g. thermal 'baths'.

A central part of this programme has been to understand the process of equilibration, i.e., how a constantly-evolving closed quantum system can behave as if relaxing to a stable equilibrium. The main insight relies on the fact [2–4] that, if measurements are limited to small subsystems or restricted sets of observables, then 'typical' pure states of large quantum systems are essentially indistinguishable from thermal states. It can then be shown [5, 6] that under very general conditions on the Hamiltonian and nearly all initial states, the system will eventually equilibrate, in the sense that an (again, restricted) set of relevant physical quantities will remain most of the time very close to fixed, 'equilibrium' values. For example, given some observable A and a system of finite but arbitrarily large size, if its expectation value $\langle A(t)\rangle$ equilibrates, then it must do so around the infinite time average (see section 5.1 of [1])

$$\begin{eqnarray}&&\overline{A}=\mathop{\mathrm{lim}}\limits_{T\to \infty }\displaystyle \frac{1}{T}{\int }_{0}^{T}\langle A(t)\rangle \,{\rm{d}}t.\end{eqnarray} \tag{ 1 }$$

If the infinite-time average fluctuation of $\langle A(t)\rangle$ around $\bar{A}$ is small, then we say that the observable A equilibrates.

One major open question is to understand the time scale at which equilibration occurs in a given system, and in particular its scaling with respect to system parameters such as its size (number of degrees of freedom). Various authors have tackled this question, e.g. [7–19], producing upper bounds that imply finite-time equilibration in various contexts (see also the 'supplementary information' section in [17] for a brief survey of the literature).

Several of these results [8–13, 17, 18] are again obtained in a typicality framework: they estimate, in various different senses, the average equilibration time of evolutions. While in many cases these calculated averages can have an impressive correspondence to experimentally measured equilibration times [17, 18], this approach also has some inherent weaknesses. First of all, it is generally believed that many specific physical conditions that are realisable in Nature or in the lab can be very far from 'typical' (for example, the actual Hamiltonians in Nature tend to have a locality structure that may be absent from most members of a mathematically generated ensemble [18]). In addition, by averaging, one loses information about the physical properties that are relevant to the equilibration time scale of any one specific evolution.

Bounds on equilibration times without taking averages have also been obtained, but only for certain restricted classes of evolutions or observables, e.g. in Malabarba et al [11], Farrelly [16], Goldstein et al [20], Monnai [7], and Santos et al [19].

Finally, there are a few works by the Bristol group [14, 15] that derive general and rigorous bounds on the equilibration times of arbitrary observables, systems and initial states, without any ensemble averaging. However, the bound in [15] scales with the inverse of the minimum energy difference (gap) in the system's spectrum. For physically realistic systems it therefore increase exponentially with the system size, and cannot be a good estimate of the actual equilibration timescale—in particular since equilibration would then not occur in the thermodynamic limit (see section 7 for details). In contrast, a bound derived in [14] can be independent of the system's size, however only in a regime that requires it to be initially in a nearly completely mixed state, failing to give a physically reasonable estimate in the case of a closed system in a pure initial state.

In this work, we seek to identify the properties of a closed quantum system which are relevant for the equilibration timescale of a given (arbitrary) observable. Our approach is more heuristic than rigorous, but it allows us to estimate a timescale which, under reasonable circumstances, depends only weakly on the system size, and thus seems to capture the relevant physics. The main insight we rely on is that equilibration is due primarily to a process of dephasing between different Fourier components of the dynamical evolution—a point that was briefly made in a classic [21], but that has apparently not been fully appreciated by the current community.

Although our argument does not result in a rigorous bound such as those in [14, 15], nor in a definite average evolution such as in [17, 18], we are able to discuss how the equilibration timescale of a given observable A depends on the physical properties of the system. Specifically, we find that the coherences of the observable of interest in the energy basis, $\langle {E}_{i}| A| {E}_{j}\rangle$, play a fundamental role. More specifically, the equilibration time depends critically on the range of energy gaps E_i − E_j for which these coherences have non-negligible values. In particular, if this range remains roughly constant as the system size increases, the same will be true for the equilibration time. As we discuss below, this indeed happens for many observables of interest in many-body systems. We illustrate these results with numerical simulations of a spin chain, finding a reasonable qualitative agreement. It should be noted that similar heuristic methods and conclusions have also recently been proposed in simultaneous, independent research by Wilming et al [22].

This point of view also gives a new understanding of some existing results. For instance, it allows us to identify the reason for the limitations of bounds and estimates such as those in [15, 20], which, while rigorous, can vastly overestimate the time scale at which equilibration occurs in realistic systems. In section 7, we argue that the reason for this behaviour is that these estimates ultimately rely on the wrong physical mechanism of equilibration, disregarding the crucial role played by dephasing.

Our main findings can be summarised as follows:

• In section 3 we discuss qualitatively why dephasing is the underlying mechanism of equilibration in closed quantum systems.
• In section 4 we develop a formalism based on the coarse-graining of functions in frequency space, which allows us to apply basic tools from Fourier transform theory, such as uncertainty relations, to equilibration related questions.
• In section 5 we determine the relevant energy scales that govern the equilibration time for a given observable, namely the energy fluctuations of the initial state and the bandwidth of the matrix of the observable in the energy basis. In particular, we give an independent proof of the fact [23] that local observables have banded matrices when written in the energy basis of a short-ranged spin Hamiltonian.
• In section 6 we illustrate our results with a numerical simulation of the XXZ model.
• In section 7 we discuss some implications of our results. In section 7.1, we reinterpret existing results from the point of view of our dephasing framework. In section 7.2 we discuss implications for the fields of quantum chaos and integrability. In particular we present two models with identical eigenbases but different level statistics that have indistinguishable dynamics over realistic time-scales.

Let us consider a closed system whose state is described by a vector in a Hilbert space of dimension d_T and whose Hamiltonian has a spectral representation

$$\begin{eqnarray}&&H=\displaystyle \sum _{k=1}^{{d}_{E}}{E}_{k}{P}_{k},\end{eqnarray} \tag{ 2 }$$

where E_k are its energies and P_k the projectors onto its eigenspaces. Note that the sum runs over d_E ≤ d_T terms, since some eigenspaces can be degenerate.

We denote the initial state by $| \psi (0)\rangle$. If the Hamiltonian has degenerate energies, we choose an eigenbasis of H such that $| \psi (0)\rangle$ has non-zero overlap with only one eigenstate $| {E}_{k}\rangle$ for each distinct energy. Choosing units such that ℏ = 1, the state at time t is then given by

$$\begin{eqnarray}&&| \psi (t)\rangle =\displaystyle \sum _{k}{c}_{k}{{\rm{e}}}^{-{{\rm{i}}{E}}_{k}t}| {E}_{k}\rangle ,\end{eqnarray} \tag{ 3 }$$

with ${c}_{k}\equiv \langle {E}_{k}| \psi (0)\rangle$. It is clear that $| \psi (t)\rangle$ evolves in the subspace spanned by $\{| {E}_{k}\rangle \}$ as if it were acted on by the non-degenerate Hamiltonian $H^{\prime} ={\sum }_{k}{E}_{k}| {E}_{k}\rangle \langle {E}_{k}|$. In this case, if the system equilibrates, the equilibrium state must be

$$\begin{eqnarray}&&\omega =\displaystyle \sum _{k}| {c}_{k}{| }^{2}| {E}_{k}\rangle \langle {E}_{k}| .\end{eqnarray} \tag{ 4 }$$

In this article, following a number of authors, [5, 6, 14, 15] we will study equilibration by focusing on observables. The idea is that a system can be considered in equilibrium if all experimentally relevant (typically, coarse-grained) observables A have equilibrated. In other words, we will focus on understanding how the expectation value $\langle A(t)\rangle$ approaches its equilibrium value $\mathrm{Tr}(A\omega )$. Note that, in order to even talk about equilibration time scales, we must assume that such a condition holds, i.e., that this observable sooner or later equilibrates.

Let us introduce the time signal of A, given the initial state $| \psi (0)\rangle$, as the distance of $\langle A(t)\rangle$ from equilibrium at time t

$$\begin{eqnarray}g(t) & := & \,\displaystyle \frac{1}{{{\rm{\Delta }}}_{A}}(\langle \psi (t)| A| \psi (t)\rangle -\mathrm{Tr}(A\omega ))\\ & = & \displaystyle \frac{1}{{{\rm{\Delta }}}_{A}}\displaystyle \sum _{i\ne j}({c}_{j}^{* }{A}_{{ji}}{c}_{i}){{\rm{e}}}^{-{\rm{i}}({E}_{i}-{E}_{j})t},\end{eqnarray} \tag{ 5 }$$

where ${A}_{{ij}}:= \langle {E}_{i}| A| {E}_{j}\rangle$ are the matrix elements of A in the energy eigenbasis, and Δ_A = a_max − a_min is the range of possible outcomes, being ${a}_{\max (\min )}$ the largest (smallest) eigenvalue of A. The denominator Δ_A is introduced to make the time signal dimensionless and satisfying $| g(t)| \leqslant 1$. Note that the time signals of two observables A and A' = b(A − a₀) are identical for ${a}_{0},b\in {\mathbb{R}}$.

We can conveniently rewrite the time signal as

$$\begin{eqnarray}&&g(t)=\displaystyle \sum _{\alpha \in { \mathcal G }}{v}_{\alpha }{{\rm{e}}}^{{\rm{i}}{G}_{\alpha }t},\end{eqnarray} \tag{ 6 }$$

where $\alpha \in { \mathcal G }=\{(i,j):i,j\in \{1,\,\ldots \,{d}_{E}\},i\ne j\}$ labels each energy gap G_α = (E_j − E_i) appearing in the system's spectrum, and where

$$\begin{eqnarray}&&{v}_{\alpha }={v}_{(i,j)}=\displaystyle \frac{{c}_{j}^{* }{A}_{{ji}}{c}_{i}}{{{\rm{\Delta }}}_{A}}.\end{eqnarray} \tag{ 7 }$$

We will refer to the complex number v_α as the amplitude of the corresponding gap G_α, and to its normalised square modulus ${q}_{\alpha }:= | {v}_{\alpha }{| }^{2}/{\sum }_{\beta }| {v}_{\beta }{| }^{2}$ as the relevance of G_α. Note that the set of relevances form a probability distribution over ${ \mathcal G }$.

A physical interpretation of this normalisation factor can be given as follows [5, 6, 24]: note that, if the system has non-degenerate gaps, then the time-averaged fluctuations of the time signal

$$\begin{eqnarray}&&{\langle | g{| }^{2}\rangle }_{T}:= \displaystyle \frac{1}{T}{\int }_{0}^{T}{\rm{d}}t| g(t){| }^{2}\end{eqnarray} \tag{ 8 }$$

satisfy the limit

$$\begin{eqnarray}&&{\langle | g{| }^{2}\rangle }_{\infty }:= \mathop{\mathrm{lim}}\limits_{T\to \infty }{\langle | g{| }^{2}\rangle }_{T}=\displaystyle \sum _{\alpha }| {v}_{\alpha }{| }^{2}.\end{eqnarray} \tag{ 9 }$$

In other words, this quantity gives the infinite time average of (the square of) the deviation of $\langle A(t)\rangle$ from its equilibrium value. We consider that the observable A equilibrates if this quantity is small, in the sense that

$$\begin{eqnarray}&&{\langle | g{| }^{2}\rangle }_{\infty }\lesssim {\left(\displaystyle \frac{\delta A}{{{\rm{\Delta }}}_{A}}\right)}^{2},\end{eqnarray} \tag{ 10 }$$

where δA is the experimentally available resolution. Furthermore, note that Δ_A/δA quantifies the the amount of different possible outcomes of a measurement and hence δA/Δ_A is expected to be much much smaller than 1. From now on, without loss of generality and for the sake of simplicity we will consider that Δ_A = 1 .

In [5, 6, 24], sufficient conditions for equilibration in the sense defined above are given by bounding the average distance from equilibrium equation (9). In particular, in [24] it is shown that for any Hamiltonian with non-degenerate gaps,

$$\begin{eqnarray}&&\displaystyle \sum _{\alpha }| {v}_{\alpha }{| }^{2}\leqslant \displaystyle \frac{1}{{d}_{{\rm{eff}}}},\end{eqnarray} \tag{ 11 }$$

where the effective dimension is defined as ${d}_{{\rm{eff}}}:= 1/({\sum }_{k}| {c}_{k}{| }^{4})=1/\mathrm{Tr}({\omega }^{2})$, which roughly speaking tells us how many eigenstates of the Hamiltonian participate in the superposition of the initial state. Thus, a large effective dimension, which is usually the case in many body systems (see appendix A and [25]), is sufficient to guarantee that condition (10) will be satisfied. Concerning the assumption of the Hamiltonian having non-degenerate gaps, it is shown in [15] that as long as there are not exponentially many degeneracies the argument stays the same.

Let us consider a situation in which the initial state is out of equilibrium, i.e. the time signal of a given operator is initially significantly larger than the equilibrium value: $| g(0)| \gg \sqrt{{\langle | g{| }^{2}\rangle }_{\infty }}$. For this to happen, the phases of the complex numbers v_α in the time signal equation (6) need to be highly synchronised. This case is presented pictorially in figure 1 (left) where the v_α's are depicted as points in the complex plane.

Zoom In Zoom Out Reset image size

Result 1 Given a time signal $g(t)={\sum }_{\alpha }{v}_{\alpha }{{\rm{e}}}^{{\rm{i}}{G}_{\alpha }t}$ with ${v}_{\alpha }=| {v}_{\alpha }| {{\rm{e}}}^{{\rm{i}}{\theta }_{\alpha }}\in {\mathbb{C}}$ being the initial amplitude of the gap ${G}_{\alpha }\in {\mathbb{R}}$, a necessary condition for the system to be initially out of equilibrium, i.e., $| g(0)|$ significantly larger than the typical equilibrium fluctuation $\sqrt{{\langle | g{| }^{2}\rangle }_{\infty }}$, is that the initial phases θ_α are not isotropically distributed but significantly synchronised. More precisely, we quantify the distance from equilibrium as

$$\begin{eqnarray}&&| g(0){| }^{2}-{\langle | g{| }^{2}\rangle }_{\infty }=2\displaystyle \sum _{\alpha \lt \beta }| {v}_{\alpha }| | {v}_{\beta }| \cos ({\theta }_{\alpha }-{\theta }_{\beta }),\end{eqnarray} \tag{ 12 }$$

which becomes negligible when the phases θ_α are isotropically distributed.

Equation (12) follows from a straightforward calculation

$$\begin{eqnarray}&&| g(0){| }^{2}=\left(\displaystyle \sum _{\alpha }{v}_{\alpha }\right){\left(\displaystyle \sum _{\beta }{v}_{\beta }\right)}^{* }={\langle | g{| }^{2}\rangle }_{\infty }+2\displaystyle \sum _{\alpha \lt \beta }| {v}_{\alpha }| | {v}_{\beta }| \cos ({\theta }_{\alpha }-{\theta }_{\beta }).\end{eqnarray} \tag{ 13 }$$

To see that isotropic randomly distributed phases give $| g(0){| }^{2}\simeq {\langle | g{| }^{2}\rangle }_{\infty }$, let ${v}_{\alpha }=| {v}_{\alpha }| {{\rm{e}}}^{{\rm{i}}{\theta }_{\alpha }}$ be a set of independent random complex variables with an isotropic probability distribution ${p}_{\alpha }(r,\theta )={p}_{\alpha }(r)=\tfrac{1}{2\pi }\delta (r-{r}_{\alpha })$, i.e. the random variable v_α has fixed modulus r_α and a random phase θ_α uniformly distributed around the circle. Then, the variance of the random variable ${\sum }_{\alpha }{v}_{\alpha }$ is

$$\begin{eqnarray}&&\mathrm{var}\left(\displaystyle \sum _{\alpha }{v}_{\alpha }\right)=\displaystyle \sum _{\alpha }\mathrm{var}({v}_{\alpha })=\displaystyle \sum _{\alpha }\langle | {v}_{\alpha }{| }^{2}\rangle =\displaystyle \sum _{\alpha }| {v}_{\alpha }{| }^{2},\end{eqnarray} \tag{ 14 }$$

where we have used the fact that the variance of a sum of independent random variables is the sum of variances and the first moments $\langle {v}_{\alpha }\rangle =0$. That is, if the phases of v_α are uniformly distributed, then the typical initial value of the time signal is $| g(0)| =| {\sum }_{\alpha }{v}_{\alpha }| \simeq {({\sum }_{\alpha }| {v}_{\alpha }{| }^{2})}^{1/2}$.

In contrast, an out-of-equilibrium initial state, for a distribution of gaps with a non-zero dispersion and a large enough system, will evolve to an equilibrium state, represented by an isotropic cloud of points centred at the origin of the complex plane as shown in figure 1 (right). This mechanism for equilibration, or more generally the vanishing amplitude of a signal over time, is usually called dephasing and is a well-known feature in many different fields of physics [19, 26–32],

The time-scale necessary for the sum in equation (6) to dephase can be estimated with the following simple argument. Suppose for simplicity that all terms have the same phase at t = 0. One way to lower bound the time needed for these phases to spread around the whole range (0, 2π) is to use the time it takes for the difference between the fastest and slowest phase to differ by 2π: tG_max − tG_min = 2π which leads to T_eq ∼ 2π/(G_max − G_min). However, the slowest and fastest gaps are not necessarily very relevant to the sum, that is, they can have a relatively small amplitude v_α. For this reason, a better estimate is obtained by replacing the denominator G_max − G_min with 2σ_G, where we define the gap dispersion σ_G to be simply the standard deviation of the G_α when weighted by their respective relevances q_α, i.e.

$$\begin{eqnarray}&&{\sigma }_{G}^{2}:= \displaystyle \sum _{\alpha }{q}_{\alpha }{({G}_{\alpha }-{\mu }_{G})}^{2}.\end{eqnarray} \tag{ 15 }$$

Since for each gap G_{(i, j)} = E_j − E_i there is also G_(j,i) = −G_(i,j), then by symmetry the average gap is μ_G = 0, and also $| {v}_{(i,j)}| =| {v}_{(j,i)}|$. By the argument made above, we can thus expect that, at least in cases where the distribution G_α is of a unimodal type, the equilibration time can be estimated by

$$\begin{eqnarray}&&{T}_{{\rm{eq}}}\sim \pi /{\sigma }_{G}.\end{eqnarray} \tag{ 16 }$$

Note that similar estimates are made in the various fields where dephasing is relevant, e.g. in [29, 31, 32]. In the next section we will give a more detailed justification for this estimate using standard tools from Fourier transform theory.

The physical properties that control the dispersion σ_G, and thus the equilibration time ${T}_{{\rm{eq}}}$, can be identified once we note that the probabilities q_α are proportional to the $| {v}_{\alpha }{| }^{2}$. Since ${v}_{\alpha }\propto {c}_{i}^{* }{c}_{j}\langle {E}_{i}| A| {E}_{j}\rangle$ (equation (7)), then σ_G is determined by: (i) the probability distribution $| {c}_{i}{| }^{2}$ for the energies; (ii) the matrix-elements of the observable A in the energy basis; and (iii) the distribution of the values of the gaps G_α = E_j − E_i themselves.

In section 5 we identify the energy scales of the observable and the initial state which are relevant to determine the dispersion of gaps σ_G. Let us now anticipate the requirements on the initial state and the observable in order that σ_G does not diverge in the thermodynamic limit, i.e. when the system size $n\to \infty$ and the equilibration time ${T}_{{\rm{eq}}}$ does not vanish. We can identify two regimes where this will happen:

• First, there are observables A that equilibrate in finite time regardless of the details of the initial state. Note that, for large many-body systems, the variance of the density of states usually scales as $\sqrt{n}$, and for generic states the width of the distribution of $| {c}_{i}{| }^{2}$ is as wide as the energy spectrum, which increases with n. This implies that, unless restrictions are placed on A, the gap variance σ_G will have the same scaling and ${T}_{{\rm{eq}}}$ will vanish in the thermodynamic limit. In order to avoid this problem, the variance of the distribution of v_α should not increase with n. This requires that the distribution $\langle {E}_{i}| A| {E}_{j}\rangle$ should decrease (or become null) for large values of E_j − E_i—in other words, A needs to be banded in the energy basis.
• Conversely, suppose the initial state has support over an energy range that does not scale with n. This can happen, for example, in the case of a so-called 'local quench' [33], when a local subsystem of fixed dimension is excited regardless of the full size of the system. In this case, σ_G will also at worst be independent of n, and so even observables with long-range coherence between very different energies, which would otherwise equilibrate quickly, will now take a finite time ${T}_{{\rm{eq}}}$.

Finally, it is also important to mention that the gap dispersion σ_G may actually decay to zero with n, leading to equilibration times that become very long. As an example of this situation consider two subsystems of increasing size n interacting through a spatially localised border of fixed size. The coherence in the energy basis of the interaction Hamiltonian is bounded by the operator norm of such an interaction, matching with the intuition that the stronger the interaction between systems, the faster the relaxation, and vice versa. By rescaling the global Hamiltonian, we can see that the interaction terms become relatively weaker as n grows and thereby the equilibration slower.

In this section we give further substance to the above heuristic argument to estimate the equilibration time-scale by means of Fourier transform techniques.

Let us first give a general idea of our approach. Suppose the time-signal g(t) decayed more or less steadily to zero, and stayed there. If so, then a good estimate for the equilibration time scale would be given by a few multiples of the standard deviation Δt defined by

$$\begin{eqnarray}&&{\rm{\Delta }}{t}^{2}:= \displaystyle \frac{1}{{\parallel g\parallel }_{2}}\int {\rm{d}}t| g(t){| }^{2}{(t-{\mu }_{t})}^{2}\end{eqnarray} \tag{ 17 }$$

with ${\mu }_{t}:= \int {\rm{d}}t| g(t){| }^{2}t/{\parallel g\parallel }_{2}$ (see figure 2).

Zoom In Zoom Out Reset image size

In such a case, following the spirit of our previous heuristic example, it would also be tempting to estimate the order of magnitude of Δt by taking the inverse of the spectral variance Δω of the signal. Indeed, this can be justified if we recall the standard uncertainty principle of Fourier analysis [34]

$$\begin{eqnarray}&&{\rm{\Delta }}t\cdot {\rm{\Delta }}\omega \geqslant 1/2.\end{eqnarray} \tag{ 18 }$$

Of course, this is a lower bound, that is saturated exactly only in the case of a Gaussian spectrum. However, it will be nearly saturated (Δt · Δω = c₁ , where c₁ is a constant of order 1), when the spectrum is unimodal and without long tails. In this case, we can also expect the time signal to decrease to a very small value after a time c₂Δt, for some small multiple c₂ again of order 1. Taking both multiples together, we expect a good estimate for the equilibration time to be of order $\sim \tfrac{{c}_{1}{c}_{2}}{2{\rm{\Delta }}\omega }$.This is satisfied by equation (16), which we will therefore continue to take as our estimate.

Unfortunately, for finite systems, our initial assumption of steady decay does not apply. The time signal has recurrences, that is, a long time after the dephasing has occurred and the system has equilibrated, the phases get again aligned (synchronised) in the complex plane, and the signal regains strength (figure 2). In order to avoid this problem, in the next subsections we introduce a coarse-grained version of the signal spectrum, which dampens out the recurrences. This allows us to exploit the uncertainty principle to estimate the equilibration time-scales, as described above. Under some mild conditions, we show that the equilibration time-scale estimated provided by this procedure coincides with the one previously given by the heuristic argument of points dephasing in the complex plane.

4.1. The frequency signal

We define the frequency signal, $\tilde{g}(\omega )$, as the Fourier transform of the time signal g(t)

$$\begin{eqnarray}&&\tilde{g}(\omega ):= { \mathcal F }[g](\omega )=\displaystyle \frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{\infty }g(t){{\rm{e}}}^{-{\rm{i}}\omega t}{\rm{d}}{t}\end{eqnarray} \tag{ 19 }$$

which roughly speaking tells us the relevance with which every frequency contributes to the time signal. When both the time and frequency signals are square-integrable ($g,\tilde{g}\in {L}^{2}$ ), the standard uncertainty principle of equation (18) applies, where

$$\begin{eqnarray}&&{\rm{\Delta }}{\omega }^{2}:= \displaystyle \frac{1}{{\parallel \tilde{g}\parallel }_{2}}\int {\rm{d}}\omega | \tilde{g}(\omega ){| }^{2}{(\omega -{\mu }_{\omega })}^{2}\end{eqnarray} \tag{ 20 }$$

with ${\mu }_{\omega }:= \int {\rm{d}}\omega | \tilde{g}(\omega ){| }^{2}\omega /{\parallel \tilde{g}\parallel }_{2}$.

However, in the case of time signals such as in equation (6), $g(t)\,\not\hspace{-2pt}{\in }\,{L}^{2}$, as can be seen from equation (9) (the integral diverges proportional to T). The same is true for the frequency signal, since

$$\begin{eqnarray}&&\tilde{g}(\omega )={ \mathcal F }[g](\omega )=\displaystyle \sum _{\alpha }{v}_{\alpha }\delta (\omega -{G}_{\alpha }),\end{eqnarray} \tag{ 21 }$$

where δ(x) is the Dirac delta distribution. Hence, the uncertainty principle in equation (18) cannot be directly applied.

It is worth noting here that, due to the finite range of energies present in our system, there is an asymmetry between the uncertainties in time and frequency of the signal g(t). On the one hand, the uncertainty in frequency can still be well-defined. To see how, recall first that, for g(t) ∈ L², it is possible to write the moments of a frequency signal in terms of the corresponding time signal and its derivatives, e.g.

$$\begin{eqnarray}&&\langle {\omega }^{2}\rangle =-\displaystyle \frac{{\int }_{-\infty }^{\infty }{g}^{* }(t)g^{\prime\prime} (t){\rm{d}}{t}}{{\int }_{-\infty }^{\infty }| g(t){| }^{2}{\rm{d}}{t}}.\end{eqnarray} \tag{ 22 }$$

In our case, although each of these integrals diverges, their ratio does have a well-defined limit, in the sense that

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{T\to \infty }-\displaystyle \frac{{\int }_{-T}^{T}{g}^{* }(t)g^{\prime\prime} (t){\rm{d}}{t}}{{\int }_{-T}^{T}| g(t){| }^{2}{\rm{d}}{t}}=\displaystyle \frac{\displaystyle \sum _{\alpha }| {v}_{\alpha }{| }^{2}{G}_{\alpha }^{2}}{\displaystyle \sum _{\alpha }| {v}_{\alpha }{| }^{2}}.\end{eqnarray} \tag{ 23 }$$

Taking then this limit as the appropriate definition of $\langle {\omega }^{2}\rangle$ in this case, and noting that, in the same sense, $\langle \omega \rangle =0$, we obtain that ${\rm{\Delta }}\omega =\sqrt{\langle {\omega }^{2}\rangle }$ is indeed precisely equal to the gap dispersion σ_G defined in equation (15).

On the other hand, though, the value of Δt diverges, even when taking limits in the same sense above. This can be understood physically due to the previously mentioned recurrences in the time signal. Indeed, g(t) is a quasi-periodic function that experiences, over an infinitely large time interval, an infinite number of recurrences to a value arbitrary close to its initial one [35, 36].

4.2. The coarse-grained signal

We now define the notion of coarse-graining, in which we introduce a microscopic energy scale below which the fine-grained details of the spectrum are washed out. As we show below, this is done by replacing the discrete spectrum present in equation (21) by a suitable smooth version.

As previously mentioned, such coarse graining of the frequency signal dampens the time signal g(t), removing the recurrences seen in figure 2 and making g(t) and $\tilde{g}(\omega )$ belong to L². We shall see later that it will also allow us to exploit certain existing statements concerning the shapes of energy densities and density of states of realistic initial states and short-ranged local Hamiltonians [37], which will justify our assumption of quasi-saturating the uncertainty bound.

An important issue in the coarse-graining is obviously the choice of the energy scale . This will be discussed later in detail but we can already understand that in order to remove the recurrences, the discreteness of the frequency signal has to be removed, which implies an much larger than the separation between consecutive gaps.

Mathematically, the coarse-graining is accomplished by convolving the frequency signal with an appropriate window function h(x), which is only nonzero over an interval of size O(). In our case, we find it convenient to choose h(x) = C N(x), where

$$\begin{eqnarray}&&{N}_{\sigma }(x):= \displaystyle \frac{1}{\sqrt{2\pi }\sigma }{{\rm{e}}}^{-\tfrac{{x}^{2}}{2{\sigma }^{2}}}\end{eqnarray} \tag{ 24 }$$

is the normalised Gaussian distribution centred at the origin and with standard deviation σ, and C a constant that is determined below.

With this spirit, the -coarse-grained version of the frequency signal is defined as

$$\begin{eqnarray}&&{\tilde{g}}_{\epsilon }(\omega ):= ({h}_{\epsilon }* \tilde{g})(\omega )={\int }_{-\infty }^{\infty }{\rm{d}}\omega ^{\prime} {h}_{\epsilon }(\omega ^{\prime} -\omega )\tilde{g}(\omega ^{\prime} ).\end{eqnarray} \tag{ 25 }$$

If $\tilde{g}(\omega )$ is given by equation (21), then the -coarse-grained frequency signal is

$$\begin{eqnarray}&&{\tilde{g}}_{\epsilon }(\omega )=\displaystyle \sum _{\alpha }{v}_{\alpha }\ {h}_{\epsilon }(\omega -{G}_{\alpha }).\end{eqnarray} \tag{ 26 }$$

In other words, coarse-graining corresponds to widening each Dirac-δ in the original spectral function into a Gaussian of O() width (figure 3). Note that, in doing this, we remove fine details of the spectrum such as the level statistics. Furthermore, unlike $\tilde{g}(\omega )$, ${\tilde{g}}_{\epsilon }(\omega )$ is square integrable and lies in L² :

$$\begin{eqnarray}&&\int {\rm{d}}\omega | {\tilde{g}}_{\epsilon }(\omega ){| }^{2}=\displaystyle \sum _{\alpha ,\beta }{v}_{\alpha }{v}_{\beta }^{* }\ ({h}_{\epsilon }\ast {h}_{\epsilon })({G}_{\alpha }-{G}_{\beta })={C}^{2}\displaystyle \sum _{\alpha ,\beta }{v}_{\alpha }{v}_{\beta }^{* }\ {N}_{\sqrt{2}\epsilon }({G}_{\alpha }-{G}_{\beta })\lt \infty ,\end{eqnarray} \tag{ 27 }$$

where we have used the fact that $({h}_{\epsilon }\ast {h}_{\epsilon })(x)={C}^{2}{N}_{\sqrt{2}\epsilon }(x)$.

Zoom In Zoom Out Reset image size

A coarse-grained frequency signal defines a coarse-grained time signal given by

$$\begin{eqnarray}&&{g}_{\epsilon }(t)={{ \mathcal F }}^{-1}[{\tilde{g}}_{\epsilon }](t)={{ \mathcal F }}^{-1}[{h}_{\epsilon }](t)\cdot g(t)=\displaystyle \frac{C}{\sqrt{2\pi }}{{\rm{e}}}^{-\tfrac{1}{2}{\epsilon }^{2}{t}^{2}}g(t),\end{eqnarray} \tag{ 28 }$$

where we have used the convolution theorem for Fourier transforms. The constant C is fixed by imposing that the time signal is not affected by coarse-graining the frequency signal in time scales t ≪ ⁻¹ i.e. g(0) = g(0). This leads to $C=\sqrt{2\pi }$ and

$$\begin{eqnarray}&&{h}_{\epsilon }(\omega )=\sqrt{2\pi }{N}_{\epsilon }(\omega )=\displaystyle \frac{1}{\epsilon }{{\rm{e}}}^{-\tfrac{{\omega }^{2}}{2{\epsilon }^{2}}},\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{g}_{\epsilon }(t)={{\rm{e}}}^{-\tfrac{1}{2}{\epsilon }^{2}{t}^{2}}g(t).\end{eqnarray} \tag{ 30 }$$

Even if the time signal g(t) equilibrates after some time, we know that it must eventually have recurrences. To determine the equilibration time-scale from g(t), we need that ⁻¹ is much greater than the equilibration timescale, but much smaller than the recurrence timescale. In this way, we ensure that the coarse-grained time signal will be indistinguishable from the original one during the equilibration process, but unlike the latter will then decay to zero.

4.3. The coarse-grained density of relevant gaps

We now focus on the properties of the variance of gap values with respect to the coarse-grained frequency signal ${\tilde{g}}_{\epsilon }(\omega )$, i.e. 

$$\begin{eqnarray}&&{\rm{\Delta }}{\omega }_{\epsilon }^{2}:= \displaystyle \frac{1}{{\parallel {\tilde{g}}_{\epsilon }\parallel }_{2}}\int {\rm{d}}\omega | {\tilde{g}}_{\epsilon }(\omega ){| }^{2}{(\omega -{\tilde{\mu }}_{\omega })}^{2}\end{eqnarray} \tag{ 31 }$$

with ${\tilde{\mu }}_{\omega }:= \int {\rm{d}}\omega | {\tilde{g}}_{\epsilon }(\omega ){| }^{2}\omega /{\parallel {\tilde{g}}_{\epsilon }\parallel }_{2}$. We will refer to Δω as the 'dispersion of relevant gaps', and to the weight

$$\begin{eqnarray}&&| {\tilde{g}}_{\epsilon }(\omega ){| }^{2}=\displaystyle \sum _{\alpha ,\beta }{v}_{\alpha }{v}_{\beta }^{* }{h}_{\epsilon }(\omega -{G}_{\alpha }){h}_{\epsilon }(\omega -{G}_{\beta })\end{eqnarray} \tag{ 32 }$$

as the density of relevant gaps.

Our goal is to show that, under a wide range of choices of and of physically relevant circumstances: (i) Δω² is very close to the gap dispersion σ_G² of the original signal, as defined in equation (15), and at the same time (ii) the inverse Δω⁻¹ is a good estimate for the equilibration time.

Note first that, by construction, for every gap G_α = E_j − E_i in equation (6) (and, by extension, equation (26)), its negative ${G}_{\bar{\alpha }}:= {E}_{i}-{E}_{j}$ also appears, with $| {v}_{\bar{\alpha }}| =| {v}_{\alpha }|$. Hence the 'average gap' ${\tilde{\mu }}_{\omega }$ vanishes for any and the variance Δω² is equivalent to $\langle {\omega }^{2}\rangle$, i.e.

$$\begin{eqnarray}&&{\rm{\Delta }}{\omega }_{\epsilon }^{2}=\displaystyle \frac{\int {\omega }^{2}| {\tilde{g}}_{\epsilon }(\omega ){| }^{2}{\rm{d}}\omega }{\int | {\tilde{g}}_{\epsilon }(\omega ){| }^{2}{\rm{d}}\omega }.\end{eqnarray} \tag{ 33 }$$

After some straightforward manipulation, using equations (29) and (32), we obtain

$$\begin{eqnarray}&&{\rm{\Delta }}{\omega }_{\epsilon }^{2}=\displaystyle \frac{{\sum }_{\alpha }| {v}_{\alpha }{| }^{2}({G}_{\alpha }^{2}+{\epsilon }^{2}/2)+\tfrac{1}{4}{\sum }_{\alpha \ne \beta }{v}_{\alpha }{v}_{\beta }^{* }[{({G}_{\alpha }+{G}_{\beta })}^{2}+2{\epsilon }^{2}]{{\rm{e}}}^{-\tfrac{{({G}_{\alpha }-{G}_{\beta })}^{2}}{4{\epsilon }^{2}}}}{{\sum }_{\alpha }| {v}_{\alpha }{| }^{2}+{\sum }_{\alpha \ne \beta }{v}_{\alpha }{v}_{\beta }^{* }\,{{\rm{e}}}^{-\tfrac{{({G}_{\alpha }-{G}_{\beta })}^{2}}{4{\epsilon }^{2}}}}.\end{eqnarray} \tag{ 34 }$$

It can be easily checked that, if $\epsilon \to 0$, this expression indeed reduces to equation (15). More specifically, Δω² will be very close to σ_G² for all $\epsilon \ll \min ({G}_{\alpha }-{G}_{\beta })$. Indeed, in this limit the Gaussian window function h(ω) becomes negligibly thin with respect to the smallest separation between gaps, and the coarse-grained spectrum ${\tilde{g}}_{\epsilon }(\omega )$ resembles the original discrete spectrum g(ω). Precisely for this reason, however, this limit is of little use to our goals. Another way of putting this is that, for such small values of the coarse-grained time signal (equation (30)) does not have time to decay before the recurrence timescale of the original signal, which is of order $\min {({G}_{\alpha }-{G}_{\beta })}^{-1}$.

To make further progress at this point, it is necessary to assume some features about the amplitudes v_α. Otherwise, a fine tuning between phases and modulus of v_α can make $| {\tilde{g}}_{\epsilon }(\omega )|$ have an arbitrary behaviour, preventing any general statement concerning the equilibration of g(t).

Inspired by [25], we will adopt a weak-typicality point of view: let us assume that the evolution we are considering is drawn from an ensemble for which the v_α's are describable by some smooth functions plus stochastic fluctuations. Note that we do not assume a uniform ensemble over all states (or any one specific ensemble), as is the case in most typicality studies [8–13, 17, 18], merely one for which the resulting distribution over the v_α's has some very general features which are described below. In the spirit of statistical physics, we basically replace complexity by apparent randomness. In most situations, the description of the gap relevances v_α in terms of a smooth function plus stochastic fluctuations is a consequence of the energy level populations c_i and the matrix-elements of the observable A_ij having this same behaviour. In the next section we discuss under which conditions this is indeed the case.

In the following we show that the process of coarse-graining removes the fluctuations and makes the density of relevant gaps $| {\tilde{g}}_{\epsilon }(\omega ){| }^{2}$ have a smooth behaviour.

(Coarse-grained frequency signal).

Result 2 Let us consider the amplitudes v_α of the gaps G_α to be described by

$$\begin{eqnarray}&&{v}_{\alpha }=v({G}_{\alpha })+\delta {v}_{\alpha },\end{eqnarray} \tag{ 35 }$$

where $v(\omega )=r(\omega ){{\rm{e}}}^{{\rm{i}}\theta (\omega )}$, with r(ω) and θ(ω) two functions with a Lipshitz constant upper bounded by some K ≪ ⁻¹, and where δv_α are independent random variables that average to zero $\langle \delta {v}_{\alpha }\rangle =0$ and $\langle \delta {v}_{\alpha }\delta {v}_{\beta }\rangle ={\gamma }^{2}({G}_{\alpha }){\delta }_{\alpha \beta }$. The variance γ²(ω) is a function that represents the strength of the fluctuations and also has a Lipshitz constant upper bounded by K. Then, with a high probability $\mathrm{erf}(m)\geqslant 1-\exp (-{m}^{2})$, the density of relevant gaps fulfils the following bound

$$\begin{eqnarray}&&| {\tilde{g}}_{\epsilon }(\omega )-\sqrt{2\pi }v(\omega ){\rho }_{\epsilon }(\omega )| \leqslant {\rho }_{\epsilon }(\omega )\left({c}_{1}K\epsilon +{\pi }^{1/4}m\displaystyle \frac{\gamma (\omega )}{\sqrt{\epsilon {\rho }_{\epsilon }(\omega )}}\right),\end{eqnarray} \tag{ 36 }$$

where c₁ > 0 is a constant, ${\rho }_{\epsilon }(\omega )=(\rho * {N}_{\epsilon })(\omega )$ is the coarse-grained density of gaps, with $\rho (\omega )={\sum }_{\alpha }\delta (\omega -{G}_{\alpha })$ being the density of gaps. The coarse-grained density of gaps ρ(ω) describes how many gaps are -close to the frequency ω. That is, the process of coarse-graining washes out the fluctuations δv_α turning the coarse-grained frequency signal into a smooth function

$$\begin{eqnarray}&&{\tilde{g}}_{\epsilon }(\omega )\simeq \sqrt{2\pi }v(\omega ){\rho }_{\epsilon }(\omega ),\end{eqnarray} \tag{ 37 }$$

where the meaning of the approximation is made precise in equation (36).

The proof is tedious and non-illuminating so we give it in appendix B.

Note that the error in (36) has two components. The term K is due to the variation of the 'smooth' functions with Lipshitz constant smaller than K within an interval of width O(). The second component $\gamma (\omega )/\sqrt{\epsilon {\rho }_{\epsilon }(\omega )}$ is the fluctuation that shrinks according to the central limit theorem with $\sqrt{\epsilon {\rho }_{\epsilon }(\omega )}$, where ρ(ω) is the number of gaps G_α within an interval O().

In order for equation (36) to be meaningful and the error bound small, an optimisation over is needed. It is easy to see that for the error to be small, the parameter should be much larger than the spacing between consecutive gaps and much smaller than the inverse of the Lipshitz constants of the continuous functions v(ω) and γ(ω), i.e.

$$\begin{eqnarray}&&{G}_{\alpha +1}-{G}_{\alpha }\ll \epsilon \ll {K}^{-1}.\end{eqnarray} \tag{ 38 }$$

Here 'consecutive' refers to gaps ordered by size, and we use the informal notation 'G_α+1' as a shorthand for 'the gap immediately larger than G_α'. Note that, in practice, the lower bound above should be applied only for consecutive relevant gaps, i.e., gaps that whose amplitudes v_α make a non-negligible contribution to the sum in equation (6).

For many-body systems, the number of gaps increases exponentially in the system size n, and energy differences between consecutive gaps shrink exponentially to zero in n. If for example these gaps all have roughly equal (exponentially small) relevances q_α, then can also be taken exponentially to zero, as long as this is done at a slower rate than the difference in gaps. This makes the time signal and its coarse-grained version indistinguishable in any realistic time-scale.

In sum, if the gap relevances v_α can be described by a continuous part plus a fluctuating part, as in equation (35), then the coarse-grained density of relevant gaps is a smooth function given by

$$\begin{eqnarray}&&| {\tilde{g}}_{\epsilon }(\omega ){| }^{2}\simeq 2\pi | v(\omega ){| }^{2}{\rho }_{\epsilon }{(\omega )}^{2}.\end{eqnarray} \tag{ 39 }$$

This will be particularly useful in the following section, where we apply these ideas in the case of many body systems described by short-ranged Hamiltonians.

It is worth noting that the factorisation in equation (37) is also automatically obtained if one assumes, as is often done, that in the thermodynamic limit $n\to \infty$ the discrete gap spectrum may be replaced by a smooth continuous gap density. i.e. taking

$$\begin{eqnarray}&&g(t)=\displaystyle \sum _{\alpha }{v}_{\alpha }{{\rm{e}}}^{{\rm{i}}{G}_{\alpha }t}\to \int v(\omega ){{\rm{e}}}^{{\rm{i}}\omega t}{\rho }_{G}(\omega ){\rm{d}}\omega .\end{eqnarray} \tag{ 40 }$$

Comparing with the definition in equation (19), we see that in this case the frequency signal is again $\tilde{g}(\omega )=\sqrt{2\pi }\,v(\omega ){\rho }_{G}(\omega )$. This indicates that the assumptions we have made concerning the smoothness of v_α are not severe. However, the point of attaining this relation via coarse-graining, while maintaining a finite dimension n, is that it allows us to control how the equilibration timescale scales with the system size, and in particular to understand how this scaling depends on the energy scaling of the relevant observable. We turn to this question in the next section.

In this section, we focus on the particular but relevant case of short-range Hamiltonians and initial states that have a finite correlation length. For such systems, we express the density of relevant gaps in terms of the energy density of the initial state and the function that describes the matrix-elements of the observable. By doing so, we identify the energy scales that determine the dispersion of gaps. We find that there are mainly two relevant energy scales: the energy fluctuations of the initial state and the bandwidth of the matrix of the observable A in the Hamiltonian eigenbasis. In the case of systems globally out of equilibrium, only those observables that are banded in the Hamiltonian basis can be observed out of equilibrium for a non-negligible time.

Local Hamiltonian. Let us define a short-ranged or local Hamiltonian of a spin lattice system, i.e., acting on a Hilbert space ${ \mathcal H }={\bigotimes }_{x\in V}{{ \mathcal H }}_{x}$ with $\dim ({{ \mathcal H }}_{x})=d$, as

$$\begin{eqnarray}&&H=\displaystyle \sum _{u\in { \mathcal E }}{h}_{u},\end{eqnarray} \tag{ 41 }$$

where the locality structure is given by a graph $(V,{ \mathcal E })$ with a vertex set V and edge set ${ \mathcal E }$. The number of terms of the Hamiltonian is denoted by $n=| { \mathcal E }|$. We consider systems for which it is possible to define a sequence of Hamiltonians H_n of different sizes. This becomes trivial in the case of translational invariant systems and regular lattices, but also includes systems with disorder, and defects. The reason for introducing such a sequence of Hamiltonians H_n is that it allows us to define the thermodynamic limit. For simplicity, the subindex n is not explicitly written from now on.

Energy density of the initial state. The energy density f(E) of an initial state $| \psi (0)\rangle ={\sum }_{i}{c}_{i}| {E}_{i}\rangle$ is defined by

$$\begin{eqnarray}&&f(E):= \displaystyle \sum _{i=1}^{{d}_{E}}| {c}_{i}{| }^{2}\delta (E-{E}_{i}),\end{eqnarray} \tag{ 42 }$$

and, what is more relevant for us, its coarse-grained version reads ${f}_{\epsilon }(E):= ({N}_{\epsilon }* f)(E)$.

For the case of local Hamiltonians, the coarse-grained energy density of states with a finite correlation length has been proven to approach a Gaussian as the size n of the system increases (lemma 8 in [37]). More specifically, for all E,

$$\begin{eqnarray}&&{\int }_{-\infty }^{E}{\rm{d}}{E}^{\prime} {f}_{\epsilon }(E^{\prime} )\to {\int }_{-\infty }^{E}{\rm{d}}{E}^{\prime} \displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_{E}}{{\rm{e}}}^{-\tfrac{1}{2}\tfrac{{(E^{\prime} -{\mu }_{E})}^{2}}{2{\sigma }_{E}^{2}}},\end{eqnarray} \tag{ 43 }$$

where the Gaussian has mean ${\mu }_{E}=\langle \psi (0)| H| \psi (0)\rangle$ and standard deviation ${\sigma }_{E}^{2}\,=\,\langle \psi (0)| {(H-{\mu }_{E})}^{2}| \psi (0)\rangle$, and where the difference in these expressions falls, at worst, essentially as $1/\sqrt{n}$. In what follows, we only make use of the energy density inside integrals, so in practice equation (43) allows us to replace the energy density ${f}_{\epsilon }(E^{\prime} )$ by the corresponding Gaussian, with vanishing error.

Initial states that are globally out of equilibrium, e.g. globally quenched, have energy densities with mean and standard deviation that scale in the system size as ${\sigma }_{E}\propto \sqrt{n}$ and μ_E ∝ n.

If the system is at criticality and the correlations decay in a power law, this Gaussian shape cannot be guaranteed anymore. In any case, the energy fluctuations can still scale as ${\sigma }_{E}\propto \sqrt{n}$ as long as the power m of the decay is sufficiently fast, i.e. m > D + 1 where D is the spatial dimension of the lattice (see appendix A for details).

Matrix-elements of an observable in the energy basis. Taking again a weakly stochastic approach, in the same spirit of the argument used in result 2, we constrain ourselves to observables with off-diagonal matrix-elements in the Hamiltonian basis that can be described by a 'continuous' function S(E, ω) plus some fluctuations δA_ij

$$\begin{eqnarray}&&| {A}_{{ij}}{| }^{2}=S\left(\displaystyle \frac{{E}_{i}+{E}_{j}}{2},\displaystyle \frac{{E}_{i}-{E}_{j}}{2}\right)+\delta {A}_{{ij}},\end{eqnarray} \tag{ 44 }$$

where this choice of writing the arguments of the function S(E, ω) will be shown to be convenient in the following section. As in the previous section, the Lipshitz constant of S(E, ω) is assumed to be bounded by K.

Note that the so called eigenstate thermalisation hypothesis (ETH) [21] can be seen as a particular case of this assumption (44). One popular version of the ETH [38] is an ansatz on the matrix-elements of an observable A in the Hamiltonian eigenbasis,

$$\begin{eqnarray}&&{A}_{{ij}}=\langle {E}_{i}| A| {E}_{j}\rangle ={ \mathcal A }(E){\delta }_{{ij}}+\rho {(E)}^{-1}{f}_{{\rm{ETH}}}(E,\omega ){R}_{{ij}},\end{eqnarray} \tag{ 45 }$$

where E = (E_i + E_j)/2 and ω = (E_i − E_j)/2. The functions ${ \mathcal A }$ and f_ETH (E, ω) are smooth functions of their arguments, and R_ij are complex numbers randomly distributed, each with zero mean and unit variance. The essential idea is that both in (44) and (45) the off-diagonal matrix-elements of the observable can be described by a smooth function plus fluctuations that vanish when coarse-graining.

Note that in our case, in contrast to ETH, we do not assume anything about the diagonal elements of the observable in the energy basis. This is due to the fact that we are not concerned about what is the equilibrium state of the system (whether is thermal or not), but only about how long the relaxation process takes.

Now that we have introduced the energy density f(E) and the function S(E, ω) that describes the observable, we are ready to express the density of relevant gaps in terms of these functions:

Result 3 Given a Hamiltonian and an initial state with populations $| {c}_{i}{| }^{2}$, let A be an observable whose matrix-elements in the Hamiltonian eigenbasis, $| {A}_{{ij}}{| }^{2}$ , can be described by means of the smooth function S(E, ω) plus some fluctuations as in equation (44). Then, up to errors O( K), the density of relevant gaps can be written as

$$\begin{eqnarray}&&| {\tilde{g}}_{\epsilon }(\omega ){| }^{2}=2\pi {\rho }_{\epsilon }(\omega )\ {\int }_{-\infty }^{\infty }{\rm{d}}{{Ef}}_{\epsilon }(E-\omega /2){f}_{\epsilon }(E+\omega /2)S(E,\omega ).\end{eqnarray} \tag{ 46 }$$

Furthermore, if the Hamiltonian is local and the initial state has a finite correlation length, such that the coarse-grained energy density f(E) is the Gaussian (43), then

$$\begin{eqnarray}&&| {\tilde{g}}_{\epsilon }(\omega ){| }^{2}=2\pi \,{N}_{\sqrt{2}{\sigma }_{E}}(\omega )S(\omega ){\rho }_{\epsilon }(\omega ),\end{eqnarray} \tag{ 47 }$$

where the function S(ω) is defined as

$$\begin{eqnarray}&&S(\omega ):= \int {\rm{d}}E\,{N}_{\tfrac{{\sigma }_{E}}{\sqrt{2}}}(E-{\mu }_{E})S(E,\omega ),\end{eqnarray} \tag{ 48 }$$

with μ_E and σ_E the mean and standard deviation of the energy density.

Proof. Let us introduce the density

$$\begin{eqnarray}&&{R}_{\epsilon }(\omega ):= {\int }_{-\infty }^{\infty }{\rm{d}}E\,{f}_{\epsilon }(E-\omega /2){f}_{\epsilon }(E+\omega /2)S(E,\omega )\end{eqnarray} \tag{ 49 }$$

and plug in it the energy density ${f}_{\epsilon }(E):= {\sum }_{i\,=\,1}^{{d}_{E}}| {c}_{i}{| }^{2}{N}_{\epsilon }(E-{E}_{i})$. A straightforward calculation leads to

$$\begin{eqnarray}&&{R}_{\epsilon }(\omega )=\displaystyle \sum _{i,j}{N}_{\sqrt{2}\epsilon }(\omega -({E}_{i}-{E}_{j}))| {c}_{i}{| }^{2}| {c}_{j}{| }^{2}\int {\rm{d}}E\,S(E,\omega ){N}_{\tfrac{\epsilon }{\sqrt{2}}}\left(E-\displaystyle \frac{{E}_{i}+{E}_{j}}{2}\right),\end{eqnarray} \tag{ 50 }$$

where we have used that ${N}_{\epsilon }(E+\omega -{E}_{i}){N}_{\epsilon }(E-\omega -{E}_{j})={N}_{\tfrac{\epsilon }{\sqrt{2}}}\left(E-\tfrac{{E}_{i}+{E}_{j}}{2}\right){N}_{\sqrt{2}\epsilon }(\omega -({E}_{i}-{E}_{j}))$. By considering that S (E, ω) is the smooth description of the matrix elements $| {A}_{{ij}}{| }^{2}$, we get, up to errors of order K,

$$\begin{eqnarray}&&{R}_{\epsilon }(\omega )=\displaystyle \sum _{\alpha }| {v}_{\alpha }{| }^{2}{N}_{\sqrt{2}\epsilon }(\omega -{G}_{\alpha })=| v(\omega ){| }^{2}{\rho }_{\sqrt{2}\epsilon }(\omega ).\end{eqnarray} \tag{ 51 }$$

Putting equations (49) and (51) together with result 3 implies (46). Now for initial states with a Gaussian energy density as in equation (43), simple algebra leads to the identity

$$\begin{eqnarray}&&{f}_{\epsilon }(E-\omega /2){f}_{\epsilon }(E+\omega /2)={N}_{\tfrac{{\sigma }_{E}}{\sqrt{2}}}(E-{\mu }_{E}){N}_{\sqrt{2}{\sigma }_{E}}(\omega ).\end{eqnarray} \tag{ 52 }$$

Plugging (52) in (46) completes the proof.

The function S(ω) describes the average magnitude of the off diagonal matrix-elements $| {A}_{{ij}}{| }^{2}$ at a distance ω = E_i − E_j from the diagonal.

Result 3 and in particular equation (47) show that the density of relevant gaps $| {\tilde{g}}_{\epsilon }(\omega ){| }^{2}$ for local Hamiltonians and initial states with decaying correlations decomposes in the product of two densities: S(ω)ρ(ω) and ${N}_{\sqrt{2}{\sigma }_{E}}(\omega )$. The density ${N}_{\sqrt{2}{\sigma }_{E}}(\omega )$ is a Gaussian with standard deviation controlled by the energy fluctuations of the initial state σ_E , and S(ω)ρ(ω) is the density of the off diagonal elements of the observable A.

The dispersion of relevant gaps σ_G, which we expect to estimate the equilibration time, is then controlled by the smallest of the standard deviations of these two densities. In the case that the system is globally out of equilibrium, e.g., a global quench, the variance of the energy density of the initial state is extensive with the system size, and ${\sigma }_{E}\propto \sqrt{n}$. This implies the following statement:

Result 4 Given a local Hamiltonian, and an initial state with clustering of correlations let A be an observable that can be described by a smooth function S(E, ω). Then, the only one way to avoid that the dispersion of relevant gaps σ_G associated to an observable A diverges in the macroscopic limit is that the matrix-representation of the observable A in the energy basis is banded. More specifically, the density S(ω)ρ(ω) has a standard deviation σ_A that is independent of the system size

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{n\to \infty }{\sigma }_{A}=\mathop{\mathrm{lim}}\limits_{n\to \infty }\displaystyle \frac{{\int }_{\epsilon }^{\infty }{\rm{d}}\omega {(\omega -{\mu }_{A})}^{2}S(\omega ){\rho }_{\epsilon }(\omega )}{{\int }_{\epsilon }^{\infty }{\rm{d}}\omega \,S(\omega ){\rho }_{\epsilon }(\omega )}\lt \infty ,\end{eqnarray} \tag{ 53 }$$

where μ_A is the first moment.

Note that a divergent dispersion of relevant gaps σ_G is expected to imply an equilibration time that tends to zero in the thermodynamic limit. In other words, observables which are not banded in their energy representation are expected to be always equilibrated, since the amount of time that they can be out of equilibrium is negligible.

It is worth mentioning that generic observables have a flat S(ω) and fulfil ${\sigma }_{A}\propto \sqrt{n}$ due to the domination of the density of gaps ρ(ω). Thus, they will have microscopically short equilibration times [20]. Observables with the property (53) turn out to be both rare and physically relevant.

In [39] it is shown for several concrete examples that indeed the matrix elements S(ω) decrease exponentially or super-exponentially with ω from a certain threshold independent of the system size (see also section 4.3.1.2 of [38]). In the following result, we show that this is a property of any local observable, that is, operators that only act non-trivially on a finite region of the system.

Result 5 Let us consider a local Hamiltonian

$$\begin{eqnarray}&&H=\displaystyle \sum _{(x,x^{\prime} )\in { \mathcal E }}{h}_{(x,x^{\prime} )}\end{eqnarray} \tag{ 54 }$$

acting on a Hilbert space ${ \mathcal H }={\bigotimes }_{x\in V}{{ \mathcal H }}_{x}$ with $\dim ({{ \mathcal H }}_{x})=d$, with a locality structure given by a graph with a vertex set V and edge set ${ \mathcal E }$. Then, the matrix elements in the energy eigenbasis of a local operator $\langle {E}_{i}| {A}_{x}| {E}_{j}\rangle$ acting on a site x fulfil the condition

$$\begin{eqnarray}&&| \langle {E}_{i}| {A}_{x}| {E}_{j}\rangle | \leqslant \parallel {A}_{x}\parallel {{\rm{e}}}^{\mathrm{log}\left(\tfrac{{\rm{e}}({E}_{i}-{E}_{j})}{J(1+\alpha )}\right)-c({E}_{i}-{E}_{j})/J},\end{eqnarray} \tag{ 55 }$$

where $c=\mathrm{log}(1+{\alpha }^{-1})$ is the decay rate, $J={\max }_{(x,x^{\prime} )\in { \mathcal E }}{\parallel h\parallel }_{\infty }$ is the strength of the local interactions and α is the lattice animal constant [40] of the graph $(V,{ \mathcal E })$.

The proof of result 5 is presented in appendix C. In particular, note that the lattice animal constant mentioned above is a parameter that captures the connectivity of the underlying graph of the Hamiltonian. For D-dimensional cubic lattices, it can be bounded as $\alpha \leqslant 2\,D{\rm{e}}$ (lemma 2 in [40]).

In the period of finishing this manuscript we have been alerted to the existence of a result very similar to our result 5, due to Arad et al (theorem 2.1 in [23]). Both proofs are similar in spirit and give similar decay rates. For a D-dimensional cubic lattice with interactions in the edges, they obtain a decay rate of $1/(8D)$ while we get $\mathrm{log}(1+{(2{\rm{e}}D)}^{-1})\simeq 1/(2{\rm{e}}D)$. The main difference is that while we bound the number of terms by counting lattice animals, they use a combinatoric argument.

Of course, this behaviour of being banded in the energy basis extends to global operators that can be decomposed into a sum of local terms, as well as for operators that are not local in real space but in momentum space when the Hamiltonian is also local in the momentum representation. Note that, indeed, most observables considered in the literature are of this type.

Let us now consider the relevant scenario of a local quench [33], in which the system is brought out of equilibrium in only a local region of the system. In such a case, the width of the energy density of the initial state is independent of the system size and related to the operator norm of the perturbation applied on the system. Unlike the global quench scenario, now even the observables that are not banded (and are initially out of equilibrium) will take a finite non-negligible time to relax. The equilibration timescale is then governed by whichever energy scale is smallest: the energy fluctuations of the state, or the dispersion σ_A of the observable. Note that our results also allow for having equilibration times that increase with the system size, as long as either σ_A and σ_E shrink with it.

We illustrate our results using the XXZ model in a transverse field and with next-nearest-neighbour coupling. We choose to use this particular model since it is not integrable, and hence does not have an exponential number of degenerate gaps. The hamiltonian is $H=J{\sum }_{i}\,{S}_{i}^{x}{S}_{i+1}^{x}+{S}_{i}^{y}{S}_{i+1}^{y}+{\rm{\Delta }}{S}_{i}^{z}{S}_{i+1}^{z}+{J}_{2}{S}_{i}^{z}{S}_{i+2}^{z}+{h}_{z}{S}_{i}^{z}$. As our equilibrating observable we choose the magnetisation density in the x direction, ${M}^{x}={\sum }_{i}{S}_{i}^{x}/N$ and, as our initial state, the fully magnetised state in the x direction. In figure 4, we show the evolution of the time signal $| {g}_{M}(t){| }^{2}$ of this observable, in the sense of equation (5). The calculations were done using full exact diagonalization of H with Δ = 0.5, J₂ = 1.0 and h_z = 0.2, for various system sizes.

Zoom In Zoom Out Reset image size

We expect g_M(t) to go to zero when the system equilibrates. Indeed, this is what happens initially, for all system sizes, and we can notice that the equilibration time (the time when $| {g}_{M}(t){| }^{2}$ becomes negligible, just before Jt=20) does not depend much on n. Furthermore, we can compare this value with our heuristic estimate for the equilibration time, ${T}_{{\rm{eq}}}\sim \pi /{\sigma }_{G}$, where we use equation (15) to calculate σ_G from the numerically obtained eigenvalues. The results are shown in table 1. We can see that the estimated equilibration times also depend only weakly on n, and are in good agreement with the timescale indicated by figure 4.

Table 1.  Estimated equilibration times ${T}_{{\rm{eq}}}$ for the XXZ model with next-nearest-neighbour coupling and an external magnetic field (Δ = 0.5, J₂ = 1.0 and h_z = 0.2). The gap dispersion σ_G was obtained by explicitly calculating equation (15) from the numerically obtained energy spectrum.

n	${T}_{{\rm{eq}}}\sim \pi /{\sigma }_{G}$
2	21
4	19
6	20
8	22
10	23
12	24

Of course, due to the small size of the simulated systems, the time signals g_M(t) also exhibit strong fluctuations. However, one can already see that, as n increases, the size of these fluctuations tends to decrease, and their onset happens later. Our numerical results therefore seem to corroborate our expectation that the observable M^x does indeed equilibrate in the limit of large n, and that this equilibration does happen at a timescale roughly given by ${T}_{{\rm{eq}}}\sim \pi /{\sigma }_{G}$.

To better illustrate the dephasing mechanism behind the equilibration process, in figure 5 we plot the amplitudes ${v}_{\alpha }{{\rm{e}}}^{{\rm{i}}{G}_{\alpha }t}$ for this same situation, in the case n = 10. Starting from an initial condition where all the amplitudes are in phase (in this case, all real and negative), one can see them rotating at different speeds and becoming more spread out in the complex plane as time goes by, resulting in the decay seen in figure 4. Note that the approximate four-way symmetry exhibited at Jt = 20 (implying g(t) ≃ 0 at this time) is already a symptom of a future recurrence: clearly, after four times this interval, all of the amplitudes will have rotated approximately back to their initial position. Indeed, one can see in figure 4 that a recurrence occurs at around Jt = 80.

Zoom In Zoom Out Reset image size

It is also instructive to compare the exact time evolution of the magnetisation with coarse-grained versions derived according to the procedures described in section 4. In figure 6 we plot again the exact time signal g_M(t) for the chain with n = 12 spins (the dark blue line in figure 4), represented here by the black dotted curve. We also plot two coarse-grained time signals g_M(t), with  = 0.4 (red line) and  = 0.02 (blue line). These curves were obtained by Fourier transforming coarse-grained frequency spectra (such as the one in figure 3) that were numerically calculated according to equation (26) (i.e., we did not simply dampen the exact time signal using equation (30)). It can be seen that all three signals are essentially indistinguishable up to the equilibration time (the slight deviation close to t = 0 is an artifact of our having discarded terms with v_α < 10⁻⁴ when calculating the sum in equation (26)). We can see that, for small (= 0.02), the signals remain indistinguishable up to and including the recurrence time. However, by choosing a value of that is sufficiently large (0.4), the coarse-grained signal faithfully reproduces the exact one during the equilibration phase, but suppresses later recurrences. Note that must still be chosen sufficiently small ( ≪ σ_G ) in order to avoid suppressing the signal even before equilibration has occurred. As discussed in previous sections, under these circumstances we may use the width of the (square-integrable) coarse-grained signal as a measure of the equilibration time of the original signal.

Zoom In Zoom Out Reset image size

Finally, let us remark that we also obtained similar results for other values of Δ and h_z, and also for the XY model. However, in the latter case, the fluctuations do not decrease exponentially with n, but only polynomially, since there are an exponential number of degenerate gaps.

7.1. Reinterpretation of previous results

It is useful now to reinterpret some previous results on equilibration times from our dephasing point of view. For example, Short and Farrelly [15] obtain a rigorous upper bound for the equilibration time of any observable by studying the time-averaged fluctuations ${\langle | g{| }^{2}\rangle }_{T}$ in equation (8) above. They are able to determine a value T₀≳ d_E/ΔE, where ΔE is the range of energies in the system, and d_E the number of different energy levels, such that averages taken over intervals longer than T₀ are negligible. This implies that the equilibration time must be upper bounded by T₀. Unfortunately, for a typical many-body system with n degrees of freedom, ΔE scales only polynomially with n, while d_E scales exponentially with n—and thus so does T₀. In other words, although this upper bound is mathematically sound, it vastly overestimates the actual equilibration time of most observables. This suggests that its derivation must be incomplete, in the sense of missing or disregarding an essential physical ingredient [20].

We now argue that this ingredient is, in a word, dephasing. Roughly speaking, in the course of their derivation, the authors bound ${\langle | g{| }^{2}\rangle }_{T}$ by separately bounding the absolute value of every term in its Fourier expansion, disregarding the interference between different terms due to dephasing. Each of these terms, which rotate according to $\exp [{\rm{i}}({G}_{\alpha }-{G}_{\beta })T]$ does gets individually dephased, due to the time averaging, but only on a time scale t_α β ∼ 1/(G_α − G_β). The bound in [15] hence corresponds to the amount of time needed for the slowest term in the Fourier sum to average out over time. For many-body systems, the gaps G_α, and also the differences between different gaps, can be exponentially small in n, hence the exponentially long upper bound.

Although this bound is therefore much too large to be a reasonable estimate the equilibration time of most observables, it must be stressed that one can always construct a specific observable A which saturates it. This is not entirely unexpected, as the bound itself is observable-independent. In fact, in [20], Goldstein et al construct such an observable, by considering a direct sum of banded matrices in the energy basis, each of which has a bandwidth that is exponentially small in the system size. From the dephasing mindset, it is straightforward to understand what is going on in this example. Since each component of A is banded with an (exponentially) narrow bandwidth, with no coherences between different bands, the dispersion of relevant gaps, σ_G, is exponentially small. In this case, our estimate π/σ_G for the equilibration time becomes exponentially large, since the small difference in angular speeds means that it must take a long time until all points are dephased and isotropically distributed in the complex plane. Analogously, in [20] it is also shown that one can always construct an observable which equilibrates extremely fast, by defining an A which is far from banded, having coherences between vastly different energies. Again, the dephasing picture intuitively explains the reason for such quick equilibration.

We can also understand some of the results obtained in another approach to the problem of relaxation of many-body systems, namely the study of the survival probability given by the quantum fidelity ${ \mathcal F }(t):= | \langle \psi (0)| \psi (t)\rangle {| }^{2}$, often in situations where the initial state $| \psi (0)\rangle$ is generated after a sudden displacement ('quantum quench') that brings the system out of equilibrium (see [19] and references therein for a review). Note that the quantum fidelity is, up to an additive constant, equivalent to the time evolution of the observable $A=| \psi (0)\rangle \langle \psi (0)|$. In this case, then, both the energy fluctuations of the initial state and the bandwidth of the observable written in the energy basis are σ_E. Hence, our considerations within the dephasing picture predict, for local interacting lattice systems, an equilibration time determined by σ_E, in agreement with the results of [19].

Finally, our outlook and conclusions are also compatible, and in some senses complementary, to recent remarkable results by Reimann et al [17, 18]. In these works, the time signal g(t) in equations (5) and (6) above is rewritten as g(t) = c F(t) + ξ(t), where (using our notations)

$$\begin{eqnarray}&&F(t)=\displaystyle \frac{1}{{d}_{T}({d}_{T}-1)}\displaystyle \sum _{\alpha }{{\rm{e}}}^{{\rm{i}}{G}_{\alpha }t},\end{eqnarray} \tag{ 56 }$$

and c is a constant (=g(0)). Note that F(t) depends only on the gap spectrum, and on the dimension of the system, but not on the initial condition, observable or eigenbasis of the Hamiltonian. It is then proven that, if one averages the time signal over certain ensembles of Hamiltonians with fixed spectra (i.e. varying only their eigenvectors) then both the average value and the standard deviation of ξ(t) become extremely small in the thermodynamic limit. In other words, for any fixed initial condition and observable of a quantum system, and also any given fixed energy spectrum, the 'typical' time signal will always be of the form g(t) = cF(t), up to negligible error. Comparing with equation (6) above, we can see that the effect of the ensemble average is to uniformize the different amplitudes v_α, both in amplitude and in phase. As a result, the 'typical' equilibration dynamics described by equation (56) becomes a pure dephasing process, i.e, a sum of uniform-length vectors in the complex plane, all initially pointing along the positive real axis, which then rotate at different speeds.

As we have argued above, dephasing/equilibration should then occur on a timescale of order π/σ_G. Indeed, this can be seen in the examples worked out in these references. For instance, one situation considered in [17] is a system described by a continuous 'microcanonical' density of energy eigenstates, of the form $\rho (E-y)={{c}{\rm{e}}}^{-y/{k}_{B}T}$, where y > 0 and c is a normalisation constant. It is also assumed that the system's initial state has support restricted to a narrow energy band [E − ΔE, E]. In this case F(t) can be calculated exactly; in particular it is shown that, for ΔE ≫ k_BT, it has the form of a Lorentzian function,

$$\begin{eqnarray}&&F(t)=\displaystyle \frac{1}{1+{(\gamma t)}^{2}},\end{eqnarray} \tag{ 57 }$$

where γ = k_B T / ℏ.

Let us now analyse this result from the point of view of our approach. First of all, as we have noted in section 3, restricting the initial state to a finite energy band, and therefore a finite gap spectrum, ensures that, for all observables, dephasing/equilibration must occur in a finite timescale.

We can estimate this timescale from the (continuous) gap density for this model:

$$\begin{eqnarray}&&{\rho }_{G}(\omega )=\int \rho (E^{\prime} )\rho (E^{\prime} +\omega ){\rm{d}}{E}^{\prime} =\displaystyle \frac{\gamma {c}^{2}{{\rm{e}}}^{-{\rm{\Delta }}E/\gamma }}{2}\sinh \left[\displaystyle \frac{{\rm{\Delta }}E-| \omega | }{\gamma }\right]\sim \displaystyle \frac{\gamma {c}^{2}}{2}{{\rm{e}}}^{-| \omega | /\gamma },\end{eqnarray} \tag{ 58 }$$

where the last equality is valid for ΔE ≫ k_BT. Given that here the amplitudes v_α of each gap in the time signal are uniform, this expression is also proportional to the gap spectrum (i.e., the Fourier transform of F(t)), and we can calculate directly the gap dispersion

$$\begin{eqnarray*}&&{\rm{\Delta }}\omega =\displaystyle \frac{\int {\omega }^{2}{\rho }_{G}^{2}(\omega ){\rm{d}}\omega }{\int {\rho }_{G}^{2}(\omega ){\rm{d}}\omega }=\displaystyle \frac{\gamma }{\sqrt{2}}.\end{eqnarray*}$$

Hence, we estimate an equilibration time ${T}_{{\rm{eq}}}\sim \tfrac{\pi \sqrt{2}}{\gamma }=\tfrac{\pi \sqrt{2}{\hslash }}{{k}_{B}T}$. It can be seen from equation (57), that, after t = T_eq, the signal has decreased by about ∼95% from its initial value, showing that this is indeed a fair estimate of the equilibration timescale. Note that for a time signal governed by equation (57), the standard deviation Δt defined in equation (17) equals 1/γ, which is also the 'half-width at half maximum' of the Lorentzian⁶ . Hence in this case the uncertainty relation, equation (18) is not far from being saturated, which is the condition we are relying on for our estimate.

As a final remark, it is also interesting to note that in this case the roles of the time signal and spectrum are exactly reversed from the well-known example of spontaneous decay of an excited atomic or nuclear state: there it is the spectrum that is a Lorentzian, and the time signal that is an exponential decay.

7.2. Equilibration time-scales and level statistics

A consequence of the coarse-graining machinery introduced above is that the fine-grained details of the spectrum do not affect the dynamics of the system up to very long time scales. More formally, the following result states that if two coarse-grained frequency signals are close in the trace norm, then their corresponding time signals must also be point-wise close for any t < O (⁻¹):

(Spectrum and dynamics).

Result 6 Let g⁽¹⁾(t) and g⁽²⁾(t) be two time signals whose -coarse-grained Fourier transforms satisfy

$$\begin{eqnarray}&&{\parallel {\tilde{g}}_{\epsilon }^{(1)}-{\tilde{g}}_{\epsilon }^{(2)}\parallel }_{1}\leqslant {\delta }_{1}\end{eqnarray} \tag{ 59 }$$

for some δ₁ > 0, and let some 0 < δ₂ < 2 set a distinguishably threshold. Then g⁽¹⁾(t) and g⁽²⁾(t) are indistinguishable up to times $t\leqslant \sqrt{{\delta }_{2}}/\epsilon$, in the sense that

$$\begin{eqnarray}&&| {g}^{(1)}(t)-{g}^{(2)}(t)| \lt {\delta }_{1}+{\delta }_{2}.\end{eqnarray} \tag{ 60 }$$

Proof. By using twice the triangular inequality, we have

$$\begin{eqnarray}&&| {g}^{(1)}(t)-{g}^{(2)}(t)| \leqslant | {g}^{(1)}(t)-{g}_{\epsilon }^{(1)}(t)| +| {g}_{\epsilon }^{(1)}(t)-{g}_{\epsilon }^{(2)}(t)| +| {g}_{\epsilon }^{(2)}(t)-{g}^{(2)}(t)| .\end{eqnarray} \tag{ 61 }$$

The first and third terms are analogous and are bounded by equation (30) and the fact that $| g(t)| \leqslant 1$. Concerning the second term, we note that, since ${\tilde{g}}_{\epsilon }(\omega )\in {L}^{2}$, then the uniform continuity statement of functional analysis allows us to write

$$\begin{eqnarray}&&{\parallel {g}_{\epsilon }^{(1)}-{g}_{\epsilon }^{(2)}\parallel }_{\infty }=\mathop{\sup }\limits_{t}| {g}_{\epsilon }^{(1)}(t)-{g}_{\epsilon }^{(2)}(t)| \leqslant {\parallel {\tilde{g}}_{\epsilon }^{(1)}-{\tilde{g}}_{\epsilon }^{(2)}\parallel }_{1}\leqslant {\delta }_{1},\end{eqnarray} \tag{ 62 }$$

and putting everything together, we get

$$\begin{eqnarray}&&| {g}^{(1)}(t)-{g}^{(2)}(t)| \leqslant 2| 1-{{\rm{e}}}^{-\tfrac{1}{2}{\varepsilon }^{2}{t}^{2}}| +{\delta }_{1}.\end{eqnarray} \tag{ 63 }$$

Finally, the bound (60) is respected as long as $1/2{\epsilon }^{2}{t}^{2}\leqslant -\mathrm{log}(1-{\delta }_{2}/2)$. The inequality $x\leqslant -\mathrm{ln}(1-x)$ for 0 < x < 2 implies that $1/2{\epsilon }^{2}{t}^{2}\leqslant {\delta }_{2}/2\leqslant -\mathrm{log}(1-{\delta }_{2}/2)$. Thus, equation (60) is guaranteed to hold for times $t\leqslant \sqrt{{\delta }_{2}}/\epsilon$ .

Result 6 shows that two frequency signals that become very similar once they are -coarse-grained have indistinguishable dynamics up to times ⁻¹. This is particularly relevant in many-body systems where the separation between consecutive energy levels shrinks exponentially in the system size⁷ . One can then consider the possibility of two many-body Hamiltonians with qualitatively different level statistics, for example one where the distribution of gaps E_i + 1 − E_i between consecutive energy levels follows a Poisson distribution, and another with Wigner–Dyson statistics, giving rise to time signals that are nevertheless indistinguishable in practice for time up to and beyond the equilibration time.

In appendix D we present an example of this kind. We have estimated the one-norm distance between two coarse-grained frequency signals evolving according to two Hamiltonians with identical eigenbases but different level statistics (one Poissonian and the other Wigner–Dyson). For this case we find that

$$\begin{eqnarray}&&{\parallel {g}_{\epsilon }^{(1)}-{g}_{\epsilon }^{(2)}\parallel }_{1}\leqslant \displaystyle \frac{C}{\epsilon }\sqrt{\displaystyle \frac{{n}^{3}}{{d}_{\mathrm{eff}}}},\end{eqnarray} \tag{ 64 }$$

where C is a numerical constant, n is the system size, the coarse-graining parameter and d_eff the effective dimension.

Note now that, assuming that d_eff increases exponentially in n (d_eff ∼ exp(cn), for some constant c), then choosing  ∼ exp(−cn/4) in equation (64) and using result 6 implies that the dynamics would not be affected up to times t ∼ exp(cn/4). (Here we are also assuming, as in the discussion following equation (38), ≫ E_i + 1 − E_i, which also decrease exponentially. This requires choosing d_eff to increase sufficiently slowly (c sufficiently small)). If these conditions are met, we obtain time signals originating from two Hamiltonians with different level-statistics but that for all practical purposes display the same dynamical behaviour.

Recall now that it is a well-known conjecture that a Poissonian nearest-neighbour gap distribution is a manifestation of integrability, and Wigner–Dyson statistics are a signature of quantum chaos. Our example seems therefore to show that it is possible for both kinds of Hamiltonian to lead to identical time signals, with identical equilibration times, at least some specific cases.

It is less clear what will happen in a more realistic example in which the two Hamiltonians H₁ and H₂ with different level statistics also have different eigenbases (as is in practice always the case). In these situations the two Hamiltonians are related by a perturbation V which does not commute with the integrable Hamiltonian H₁ and where both H₁ and V have a locality structure of the type in (41). In such a scenario, the coarse-grained energy density of the initial state is not affected by the perturbation since it keeps the Hamiltonian local. Thus, any change in the one-norm distance between frequency signals, and thereby in the dynamical behaviour of the system, must come from a drastic change in the matrix-elements of the observable. A question that arises beyond the scope of this paper is then how integrability, non-integrability, and chaos can be identified in the behaviour of the matrix-elements of the observable in the energy basis.

7.3. The dephasing mindset for quadratic Hamiltonians

In this work we have argued that equilibration in closed quantum systems should be understood as a process of dephasing of complex numbers in the complex plane. From this mechanism, we have heuristically estimated the equilibration time-scale as roughly the inverse of the dispersion of the relevant gaps. We have seen that, under physically relevant circumstances, the equilibration time-scale estimated in this way depends at most weakly on the system size, in agreement with realistic situations. Although our argument does not result in a rigorous bound, we claim that it captures the correct way in which the time-scale depends on the physical properties of the system. In particular, we have seen that the coherences of the observables of interest in the energy basis, $\langle {E}_{i}| A| {E}_{j}\rangle$, play a fundamental role: in order to attain a finite equilibration time for generic initial states, these coherences must become small as E_i − E_j increases.

We have also observed that the size of the system only plays a role in the typical size of the fluctuations, but not in the time of equilibration, and thus small systems fail to equilibrate not because their equilibration time is large, but because their fluctuations are big. We illustrate these results with numerical simulations of spins chains.

Finally, we have applied the dephasing mindset to give an intuitive interpretation to earlier works on equilibration times. Our results satisfactorily reproduce many particular cases of determining equilibration time scales found in the literature. Of course, further work is still required to put our claims on a more rigorous foundation. For example, we conjecture that in a sufficiently wide range of locally interacting n-body systems, the coarse-grained frequency signal ${\tilde{g}}_{\epsilon }(\omega )$ may itself approach a Gaussian in the limit of large n. In this case, the Heisenberg-like uncertainty principle we have been using to heuristically estimate the equilibration timescale would become close to saturated, and would therefore indeed be a bona-fide measure for it.

We would like to thank J Eisert and his co-authors in [22] for the coordination in the publication of this manuscript and their work, and A Molnar and H Wilming for letting us know about the existence of [23] which contains a very similar theorem to our result 5. TRO and DJ are supported by the Brazilian National Institute for Science and Technology of Quantum Information (INCT-IQ). TRO also thanks the National Counsel of Technological and Scientific Development (CNPq). CC, ML and AR thank financial support from Spanish MINECO (QIBEQI FIS2016-80773-P, FISICATEAMO FIS2016-79508-P and Severo Ochoa Grant No. SEV-2015-0522), Fundació Privada Cellex, Generalitat de Catalunya (Grant No. SGR 874, 875, and CERCA Programme) and the European Commission [EQuaM (FP7-ICT-2013-C No.323714), OSYRIS (ERC-2013-AdG No. 339106), QUIC (H2020-FETPROACT-2014 No. 641122), SIQS (FP7-ICT-2011-9 No. 600645)]. AR is also supported by the Beatriu de Pinós fellowship (BP-DGR 2013) and the CELLEX-ICFO-MPQ research fellowship.

A.1. Energy uncertainty of the initial state

A standard way to bring quantum systems out of equilibrium is to initialise the system as the ground state of a local Hamiltonian, and change the Hamiltonian of the system sufficiently fast such that the state is kept unchanged. This procedure is called a quantum quench [33].

It is well known that the ground state of local Hamiltonians exhibit a clustering of correlations [42]. In fact, correlations between observables supported at different lattice points decay algebraically with the distance for critical systems (gapless Hamiltonians) and exponentially for systems off criticality (gapped Hamiltonians).

Consider a local Hamiltonian of the type (41) and a quantum state $| \psi \rangle \in { \mathcal H }$ with either exponentially or algebraically decaying correlations such that

$$\begin{eqnarray}&&| \langle {h}_{u}{h}_{v}\rangle -\langle {h}_{u}\rangle \langle {h}_{v}\rangle | \leqslant \displaystyle \frac{c}{d{(u,v)}^{D+1}},\end{eqnarray} \tag{ A.1 }$$

where c > 0, d(u, v) is the graph distance between the edges u and v, and D is the spatial dimension in which the graph can be embedded. Then, the variance of the energy distribution is upper bounded by ${\mathrm{var}}_{\psi }({H}_{n})\leqslant c^{\prime} n$ with c' > 0.

By using (A.1), it is easy to see that

$$\begin{eqnarray}&&\displaystyle \sum _{v\in E}(\langle {h}_{u}{h}_{v}\rangle -\langle {h}_{u}\rangle \langle {h}_{v}\rangle )\leqslant c^{\prime} \,\forall \ u.\end{eqnarray} \tag{ A.2 }$$

Hence, writing the variance in terms of the local terms of the Hamiltonian, and using the previous bound, one gets

$$\begin{eqnarray}&&{{\rm{var}}}_{\psi }({H}_{n})={\langle {H}^{2}\rangle }_{\psi }-{\langle H\rangle }_{\psi }^{2}=\displaystyle \sum _{u,v\in E}({\langle {h}_{u}{h}_{v}\rangle }_{\psi }-{\langle {h}_{u}\rangle }_{\psi }{\langle {h}_{v}\rangle }_{\psi })\leqslant {cn}.\end{eqnarray} \tag{ A.3 }$$

A.2. Effective dimension

In this section we provide the proof of result 2. The coarse grained Fourier transform of the time signal naturally decomposes

$$\begin{eqnarray}&&{\tilde{g}}_{\epsilon }(\omega )=\displaystyle \sum _{\alpha }v({G}_{\alpha }){h}_{\epsilon }(\omega -{G}_{\alpha })+\displaystyle \sum _{\alpha }\delta {v}_{\alpha }{h}_{\epsilon }(\omega -{G}_{\alpha }).\end{eqnarray} \tag{ B.1 }$$

The second term in the right-hand side is a random variable with variance given by

$$\begin{eqnarray*}&&\delta {v}_{\epsilon }{(\omega )}^{2}:= {\rm{var}}\left(\displaystyle \sum _{\alpha }\delta {v}_{\alpha }{h}_{\epsilon }(\omega -{G}_{\alpha })\right)=\displaystyle \sum _{\alpha }{h}_{\epsilon }^{2}(\omega -{G}_{\alpha }){\rm{var}}(\delta {v}_{\alpha })=\displaystyle \sum _{\alpha }{h}_{\epsilon }^{2}(\omega -{G}_{\alpha }){\gamma }^{2}({G}_{\alpha }).\end{eqnarray*}$$

By means of the density of gaps, the variance of the fluctuation term in equation (B.1) becomes

$$\begin{eqnarray}&&\delta {v}_{\epsilon }{(\omega )}^{2}=\displaystyle \frac{\sqrt{\pi }}{\epsilon }\int {\rm{d}}\omega ^{\prime} {N}_{\epsilon /\sqrt{2}}(\omega -\omega ^{\prime} )\rho (\omega ^{\prime} )\gamma (\omega ^{\prime} ),\end{eqnarray} \tag{ B.2 }$$

where we have taken into account that ${h}_{\epsilon }{(\omega )}^{2}=\sqrt{\pi }{\epsilon }^{-1}{N}_{\epsilon /\sqrt{2}}(\omega )$. Now we use the fact that γ(ω) is almost constant in an interval , i.e. the Lipshitz constant K of γ is such that K ≪ 1. The idea is to approximate δv(ω)² by

$$\begin{eqnarray}&&\delta {v}_{\epsilon }{(\omega )}^{2}\simeq \displaystyle \frac{\sqrt{\pi }}{\epsilon }{\gamma }^{2}(\omega )\int {\rm{d}}\omega ^{\prime} {N}_{\epsilon /\sqrt{2}}(\omega -\omega ^{\prime} )\rho (\omega ^{\prime} )=\displaystyle \frac{\sqrt{\pi }}{\epsilon }{\gamma }^{2}(\omega ){\rho }_{\epsilon /\sqrt{2}}(\omega ).\end{eqnarray} \tag{ B.3 }$$

More precisely, such approximation has 2 steps. The first step is to restrict the domain of integration around ω, where the Gaussian ${N}_{\epsilon /\sqrt{2}}(\omega -\omega ^{\prime} )$ is centred. The error of such step can be upper bounded by

$$\begin{eqnarray}&&\left|\delta {v}_{\epsilon }{(\omega )}^{2}-\displaystyle \frac{\sqrt{\pi }}{\epsilon }{\int }_{\omega -{\rm{\Delta }}}^{\omega +{\rm{\Delta }}}{\rm{d}}\omega ^{\prime} {N}_{\epsilon /\sqrt{2}}(\omega -\omega ^{\prime} )\rho (\omega ^{\prime} )\gamma (\omega ^{\prime} )\right|\leqslant \displaystyle \frac{\sqrt{\pi }}{\epsilon }{\rho }_{\epsilon /\sqrt{2}}(\omega )\alpha {{\rm{e}}}^{-2{{\rm{\Delta }}}^{2}/{\epsilon }^{2}},\end{eqnarray} \tag{ B.4 }$$

where we have used the Chernoff bound of the Gaussian distribution ${\int }_{x}^{\infty }{\rm{d}}{{yN}}_{\sigma }(y)\leqslant 2{{\rm{e}}}^{-{x}^{2}/{\sigma }^{2}}$ and α > 0 is a constant that depends on the density of gaps.

The second step is considering the worst case in the variation of γ(ω) within the interval of integration [ω − Δ, ω + Δ]. To do so, let us consider two real positive functions a(x) and b(x), where b(x) has Lipschitz constant K. Then, by using trivial calculus and the definition of Lipschitz constant, we get

$$\begin{eqnarray*}&&{\int }_{\omega -{\rm{\Delta }}}^{\omega +{\rm{\Delta }}}{\rm{d}}{xa}(x)b(x)\leqslant \left(\mathop{\max }\limits_{x\in [\omega -{\rm{\Delta }},\omega +{\rm{\Delta }}]}b(x)\right){\int }_{\omega -{\rm{\Delta }}}^{\omega +{\rm{\Delta }}}{\rm{d}}{xa}(x)\leqslant (b(\omega )+K{\rm{\Delta }}){\int }_{\omega -{\rm{\Delta }}}^{\omega +{\rm{\Delta }}}{\rm{d}}{xa}(x).\end{eqnarray*}$$

A straight forward application of the former statement leads to

$$\begin{eqnarray*}&&{\int }_{\omega -{\rm{\Delta }}}^{\omega +{\rm{\Delta }}}{\rm{d}}\omega ^{\prime} {\gamma }^{2}(\omega ^{\prime} ){N}_{\epsilon /\sqrt{2}}(\omega -\omega ^{\prime} )\rho (\omega ^{\prime} )-{\gamma }^{2}(\omega ){\int }_{\omega -{\rm{\Delta }}}^{\omega +{\rm{\Delta }}}{\rm{d}}\omega ^{\prime} {N}_{\epsilon /\sqrt{2}}(\omega -\omega ^{\prime} )\rho (\omega ^{\prime} )\leqslant {\rho }_{\epsilon /\sqrt{2}}(\omega )K{\rm{\Delta }}.\end{eqnarray*}$$

By putting everything together, we get

$$\begin{eqnarray}&&\left|\delta {v}_{\epsilon }{(\omega )}^{2}-\displaystyle \frac{\sqrt{\pi }}{\epsilon }{\gamma }^{2}(\omega ){\rho }_{\epsilon /\sqrt{2}}(\omega )\right|\leqslant \displaystyle \frac{\sqrt{\pi }}{\epsilon }{\rho }_{\epsilon /\sqrt{2}}(\omega )(2\alpha {{\rm{e}}}^{-2{{\rm{\Delta }}}^{2}/{\epsilon }^{2}}+K{\rm{\Delta }}).\end{eqnarray} \tag{ B.5 }$$

Finally, an optimisation over Δ gives

$$\begin{eqnarray}&&\left|\delta {v}_{\epsilon }{(\omega )}^{2}-\displaystyle \frac{\sqrt{\pi }}{\epsilon }{\gamma }^{2}(\omega ){\rho }_{\epsilon /\sqrt{2}}(\omega )\right|\leqslant \alpha ^{\prime} {\epsilon }^{-1}{\rho }_{\epsilon /\sqrt{2}}(\omega )K\epsilon ,\end{eqnarray} \tag{ B.6 }$$

where α' > 0 contains all the numerical factors. Equation (B.6) implies a relative error in the of the order K.

The standard deviation δv(ω) sets the order of magnitude for the fluctuations of the random variable. We will see that it is negligible in comparison with the first term in equation (B.1). Such a term can also be written in terms of the density of gaps ρ(ω)

$$\begin{eqnarray}&&{R}_{\epsilon }(\omega ):= \displaystyle \sum _{\alpha }v({G}_{\alpha }){h}_{\epsilon }(\omega -{G}_{\alpha })=\int {\rm{d}}\omega ^{\prime} \rho (\omega ^{\prime} )v(\omega ^{\prime} ){h}_{\epsilon }(\omega -\omega ^{\prime} ).\end{eqnarray} \tag{ B.7 }$$

The function R(ω) can be shown to be close to $\sqrt{2\pi }\ v(\omega ){\rho }_{\epsilon }(\omega )$. More precisely, by following an analogous argument of the one above, it can be proven that

$$\begin{eqnarray}&&| {R}_{\epsilon }(\omega )-\sqrt{2\pi }\ v(\omega ){\rho }_{\epsilon }(\omega )| \leqslant \beta {\rho }_{\epsilon }(\omega )K\epsilon ,\end{eqnarray} \tag{ B.8 }$$

where β > 0 is an order one factor.

We are now ready to bound the difference ${\tilde{g}}_{\epsilon }(\omega )-\sqrt{2\pi }v(\omega ){\rho }_{\epsilon }(\omega )$. By using equation (B.1), definition (B.7) and the triangular inequality, we get

$$\begin{eqnarray}&&| {\tilde{g}}_{\epsilon }(\omega )-\sqrt{2\pi }v(\omega ){\rho }_{\epsilon }(\omega )| \leqslant | {R}_{\epsilon }(\omega )-\sqrt{2\pi }v(\omega ){\rho }_{\epsilon }(\omega )| +\left|\displaystyle \sum _{\alpha }\delta {v}_{\alpha }{h}_{\epsilon }(\omega -{G}_{\alpha })\right|.\end{eqnarray} \tag{ B.9 }$$

Finally, from equations (B.3) and (B.8), the error is bounded by

$$\begin{eqnarray}&&| {\tilde{g}}_{\epsilon }(\omega )-\sqrt{2\pi }v(\omega ){\rho }_{\epsilon }(\omega )| \leqslant {\rho }_{\epsilon }(\omega )\left(\beta K\epsilon +\displaystyle \frac{{\pi }^{1/4}\gamma (\omega )}{\sqrt{\epsilon {\rho }_{\epsilon }(\omega )}}\right).\end{eqnarray} \tag{ B.10 }$$

Here, we have used that ${\rho }_{\epsilon /\sqrt{2}}(\omega )\leqslant {\rho }_{\epsilon /\sqrt{2}}(\omega )(1+K\epsilon )$ where K also upper-bounds the Lipshitz constant of ρ(ω).

In this section we prove result 5. We first introduce the imaginary time evolution of the observable

$$\begin{eqnarray}&&{A}_{x}(\beta ):= {{\rm{e}}}^{\beta H}{A}_{x}{{\rm{e}}}^{-\beta H}.\end{eqnarray} \tag{ C.1 }$$

By Taylor expanding of the exponentials, it can be written as

$$\begin{eqnarray}&&{A}_{x}(\beta )={A}_{x}+\displaystyle \sum _{m=1}^{\infty }\displaystyle \frac{{\beta }^{m}}{m!}[H,\,\ldots \,[H,{A}_{x}]\ldots ],\end{eqnarray} \tag{ C.2 }$$

where we have used that $[H,\,\ldots \,[H,{A}_{x}]\ldots ]={\sum }_{k=0}^{m}\left(\genfrac{}{}{0em}{}{m}{k}\right){H}^{k}{A}_{x}{(-H)}^{m-k}$.

We use now the locality of the Hamiltonian given by its underlying graph $(V,{ \mathcal E })$. In order to do so and for the sake of clarity, let us introduce some notation. We refer to edges $\lambda \in { \mathcal E }$ as letters, to the edge set ${ \mathcal E }$ as an alphabet, and call the sequences of edges as words. For any sub-alphabet $F\subset { \mathcal E }$, we denote by F^l the set of words with letters in F and length l, where the length $| w|$ of a word w is the total number of letters it contains. In turn, the set of words with letters in F and arbitrary length l is defined as ${F}^{* }:= {\bigcup }_{l=0}^{\infty }{F}^{l}$ . The size of a sub-alphabet F is denoted by $| F|$ and corresponds to the number of letters it contains. We call a sub-alphabet F an animal if it forms a set of connected edges, and a word connected or cluster $c\in {{ \mathcal E }}^{* }$ if the set of letters in c is an animal, i.e., connected. That is, clusters are connected sequences of edges where the edges can also occur multiple times, while animals are connected sub-graphs (sets of connected edges without order or repetition). For $w\in {E}^{* }$ and any sub-alphabet $G\subset { \mathcal E }$, we write $G\subset w$ if every letter in G also occurs in w. When a sub-alphabet F contains at least one edge adjacent to a vertex x, we denote it by $F\ni x$.

With the above notation, the order n commutator can be written as

$$\begin{eqnarray}&&[H,\,\ldots \,[H,{A}_{x}]\ldots ]=\displaystyle \sum _{w\in {{ \mathcal E }}^{m}}[{h}_{{w}_{1}},\,\ldots \,[{h}_{{w}_{m}},{A}_{x}]\ldots ],\end{eqnarray} \tag{ C.3 }$$

where w_k is the kth letter of the word w and ${h}_{{w}_{k}}$ the corresponding Hamiltonian term. Note that the only way for the w-commutator $[{h}_{{w}_{1}},\,\ldots \,[{h}_{{w}_{n}},{A}_{x}]\ldots ]$ to be nonzero is if w forms a single cluster with at least one of its letters adjacent to the vertex x. More formally,

$$\begin{eqnarray}&&[H,\,\ldots \,[H,{A}_{x}]\ldots ]=\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{F\subset { \mathcal E }:}{\genfrac{}{}{0em}{}{| F| \leqslant n}{F\ni x}}}\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{w\in {F}^{n}:}{F\subset w}}[{h}_{{w}_{1}},\,\ldots \,[{h}_{{w}_{n}},{A}_{x}]\ldots ],\end{eqnarray} \tag{ C.4 }$$

where the first sum runs over all the animals F that contain at least an edge adjacent to the vertex x and are smaller or equal than n. Then, the imaginary time evolution of the observable can be written as

$$\begin{eqnarray}&&{A}_{x}(\beta )={A}_{x}+\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{F\subset { \mathcal E }:}{F\ni x}}\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{w\in {F}^{* }:}{F\subset w}}\displaystyle \frac{{\beta }^{| w| }}{| w| !}[{h}_{{w}_{1}},\,\ldots \,[{h}_{{w}_{| w| }},{A}_{x}]\ldots ].\end{eqnarray} \tag{ C.5 }$$

We can upper bound the operator norm of _Ax(β) as

$$\begin{eqnarray}&&\parallel {A}_{x}(\beta )\parallel \leqslant \parallel {A}_{x}\parallel \left(1+\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{F\subset { \mathcal E }:}{F\ni x}}\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{w\in {F}^{* }:}{F\subset w}}\displaystyle \frac{{(J\beta )}^{| w| }}{| w| !}\right),\end{eqnarray} \tag{ C.6 }$$

where we have used that for any operator O, $\parallel [{h}_{{w}_{k}},O]\parallel \leqslant \parallel {h}_{{w}_{k}}\parallel \parallel O\parallel$, whose iteration implies $\parallel [{h}_{{w}_{1}},\,\ldots \,[{h}_{{w}_{| w| }},{A}_{x}]\ldots ]\parallel \leqslant \parallel {A}_{x}\parallel {J}^{| w| }$. We recall now lemma 5 from [43] which states that

$$\begin{eqnarray}&&\displaystyle \sum _{w\in {F}^{* }:F\subset w}\displaystyle \frac{| \beta J{| }^{| w| }}{| w| !}={({{\rm{e}}}^{| \beta J| }-1)}^{| F| }\end{eqnarray} \tag{ C.7 }$$

and implies

$$\begin{eqnarray}&&\parallel {A}_{x}(\beta )\parallel \leqslant \parallel {A}_{x}\parallel \left(1+\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{F\subset { \mathcal E }:}{F\ni x}}{({{\rm{e}}}^{| \beta J| }-1)}^{| F| }\right).\end{eqnarray} \tag{ C.8 }$$

By decomposing the sum over all animals which contain x according to their size, we get

$$\begin{eqnarray}&&\parallel {A}_{x}(\beta )\parallel \leqslant \parallel {A}_{x}\parallel \left(1+\displaystyle \sum _{l=1}^{\infty }\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{F\subset { \mathcal E }:}{\genfrac{}{}{0em}{}{| F| =l}{F\ni x}}}{({{\rm{e}}}^{\beta J}-1)}^{| F| }\right).\end{eqnarray} \tag{ C.9 }$$

The number of lattice animals of size l is upper-bounded by means of the so called lattice animal constant α. Given a graph $(V,{ \mathcal E })$ and denoting by a_l the number of lattice animals of size l containing the fixed vertex x, then the animal constant α is the smallest constant satisfying

$$\begin{eqnarray}&&{a}_{l}:= \displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{F\subset { \mathcal E }:}{\genfrac{}{}{0em}{}{| F| =l}{F\ni x}}}1\leqslant {\alpha }^{l}.\end{eqnarray} \tag{ C.10 }$$

Regular lattices have finite animal constants[44]. For example, the animal constant of a D-dimensional cubic lattice can be bounded as $\alpha \leqslant 2D{\rm{e}}$ (lemma 2 in [40]), where ${\rm{e}}$ is Euler's number.

Equation (C.10) allows us to bound (C.9) as

$$\begin{eqnarray}&&\parallel {A}_{x}(\beta )\parallel \leqslant \parallel {A}_{x}\parallel \displaystyle \sum _{l=0}^{\infty }{(\alpha ({{\rm{e}}}^{\beta J}-1))}^{l},\end{eqnarray} \tag{ C.11 }$$

which is a geometric series with common ratio $\alpha ({{\rm{e}}}^{2\beta J}-1)$. In order for the series to converge we require that β is such that $\alpha ({{\rm{e}}}^{\beta J}-1)\lt 1$, for which

$$\begin{eqnarray}&&\parallel {A}_{x}(\beta )\parallel \leqslant \parallel {A}_{x}\parallel \displaystyle \frac{1}{1-\alpha ({{\rm{e}}}^{\beta J}-1)}.\end{eqnarray} \tag{ C.12 }$$

We compute the absolute value of the matrix-element $\langle {E}_{i}| \cdot | {E}_{j}\rangle$ in equation (C.1)

$$\begin{eqnarray}&&{{\rm{e}}}^{\beta ({E}_{i}-{E}_{j})}| \langle {E}_{i}| {A}_{x}| {E}_{j}\rangle | =| \langle {E}_{i}| {A}_{x}(\beta )| {E}_{j}\rangle | \leqslant \parallel {A}_{x}(\beta )\parallel \end{eqnarray} \tag{ C.13 }$$

and, together with equation (C.12), it implies that

$$\begin{eqnarray}&&| \langle {E}_{i}| {A}_{x}| {E}_{j}\rangle | \leqslant \parallel {A}_{x}\parallel \displaystyle \frac{{{\rm{e}}}^{-\beta J({E}_{i}-{E}_{j})/J}}{1-\alpha ({{\rm{e}}}^{\beta J}-1)}\end{eqnarray} \tag{ C.14 }$$

for any β such that $\alpha ({{\rm{e}}}^{\beta J}-1)\lt 1$.

In order to get an explicit bound independent of β, we optimise over β. Given some energy difference E_i − E_j, we look for the β that minimises the upper bound of equation (C.14). To do so, it is useful to rewrite the bound in terms of a new parameter $z={{\rm{e}}}^{\beta J}$ such that $| \langle {E}_{i}| {A}_{x}| {E}_{j}\rangle | \leqslant \parallel {A}_{x}\parallel f(z)$ with

$$\begin{eqnarray}&&f(z):= \displaystyle \frac{1}{\alpha (1+{\alpha }^{-1}-z){z}^{\omega }}\end{eqnarray} \tag{ C.15 }$$

and ω = (E_i − E_j)/J. By means of an optimisation over z, we find the minimum at $z=\tfrac{\omega }{\omega +1}\tfrac{\alpha +1}{\alpha }$ which fulfils the convergence condition 1 < z < 1 + α⁻¹ as long as ω > α. The minimum value of f(z_*) becomes

$$\begin{eqnarray}&&f({z}_{* })=\displaystyle \frac{\omega +1}{\alpha +1}{\left(\displaystyle \frac{\omega }{\omega +1}\right)}^{-\omega }{(1+{\alpha }^{-1})}^{-\omega }\leqslant \displaystyle \frac{\omega +1}{\alpha +1}{\rm{e}}{(1+{\alpha }^{-1})}^{-\omega },\end{eqnarray} \tag{ C.16 }$$

where the inequality comes from the fact that ${(\omega /(\omega +1))}^{-\omega }\leqslant {\rm{e}}$. The matrix element $| \langle {E}_{i}| {A}_{x}| {E}_{j}\rangle |$ is then bounded by

$$\begin{eqnarray}&&| \langle {E}_{i}| {A}_{x}| {E}_{j}\rangle | \leqslant \parallel {A}_{x}\parallel {{\rm{e}}}^{\mathrm{log}\left(\tfrac{{\rm{e}}({E}_{i}-{E}_{j})}{J(1+\alpha )}\right)-c({E}_{i}-{E}_{j})/J},\end{eqnarray} \tag{ C.17 }$$

where $c=\mathrm{log}(1+{\alpha }^{-1})$ is the decay rate.

Two Hamiltonians H⁽¹⁾ and H⁽²⁾ with identical eigenstates but eigenvalues with different level statistics give rise to time signals for equal initial states and observables given by

$$\begin{eqnarray}&&{g}^{(1)}(t)=\displaystyle \sum _{\alpha }{v}_{\alpha }{{\rm{e}}}^{{\rm{i}}{G}_{\alpha }^{(1)}t},\end{eqnarray} \tag{ D.1 }$$

and analogously for g⁽²⁾(t). Here, G_α⁽¹⁾ are the gaps associated the spectrum of H₁, and the relevances v_α are the same for both systems. These time-signals have Fourier transforms

$$\begin{eqnarray}&&{\tilde{g}}^{(1)}(\omega )=\displaystyle \sum _{\alpha }{v}_{\alpha }\delta (\omega -{G}_{\alpha }^{(1)})\end{eqnarray} \tag{ D.2 }$$

that once coarse grained read

$$\begin{eqnarray}&&{\tilde{g}}_{\epsilon }^{(1)}(\omega )=\displaystyle \sum _{\alpha }{v}_{\alpha }{h}_{\epsilon }(\omega -{G}_{\alpha }^{(1)}).\end{eqnarray} \tag{ D.3 }$$

The one-norm distance between the frequency signals ${\tilde{g}}_{\epsilon }^{(1)}$ and ${\tilde{g}}_{\epsilon }^{(2)}$ is defined by

$$\begin{eqnarray}&&\parallel {\tilde{g}}_{\epsilon }^{(1)}-{\tilde{g}}_{\epsilon }^{(2)}{\parallel }_{1}:= \int {\rm{d}}\omega | {g}^{(1)}(\omega )-{g}^{(2)}(\omega )| .\end{eqnarray} \tag{ D.4 }$$

Let us now introduce the function

$$\begin{eqnarray}&&{u}_{\delta }(\omega )=2\epsilon {h}_{\epsilon }(\omega ){h}_{\epsilon }(\delta )\sinh \left(\displaystyle \frac{\delta }{\epsilon }\displaystyle \frac{\omega }{\epsilon }\right).\end{eqnarray} \tag{ D.5 }$$

By plugging (D.3), we can write the difference in frequency signals as

$$\begin{eqnarray}&&{g}^{(1)}(\omega )-{g}^{(2)}(\omega )=\displaystyle \sum _{\alpha }{v}_{\alpha }\ {u}_{\delta {G}_{\alpha }}(\omega -{\bar{G}}_{\alpha }),\end{eqnarray} \tag{ D.6 }$$

where we have used that ${u}_{\delta {G}_{\alpha }}(\omega -{\bar{G}}_{\alpha })={h}_{\epsilon }(\omega -{G}_{\alpha }^{(1)})-{h}_{\epsilon }(\omega -{G}_{\alpha }^{(2)})$ with $\delta {G}_{\alpha }={G}_{\alpha }^{(1)}-{G}_{\alpha }^{(2)}$ and ${\bar{G}}_{\alpha }=({G}_{\alpha }^{(1)}+{G}_{\alpha }^{(2)})/2$.

In order to understand how the variables ${\bar{G}}_{\alpha }$ and δG_α behave, let us introduce the spectrum of H₁ (with Poisson level statistics) by $\{{E}_{k}^{(1)}\}{}_{k=1}^{{d}_{E}}$, and the one of H₂ (with Wigner–Dyson level statistics) by $\{{E}_{k}^{(2)}\}{}_{k=1}^{{d}_{E}}$. The energy levels can be iteratively constructed as

$$\begin{eqnarray}&&{E}_{k+1}={E}_{k}+{s}_{k},\end{eqnarray} \tag{ D.7 }$$

where s_k $\forall \ k$ are i.i.d. random variables sampled from a Poisson distribution in the integrable case and from a Wigner–Dyson distribution in the chaotic case. Both distributions are assumed to have the same mean μ_s and standard deviation $\sqrt{{\mu }_{s}}$ . The gaps read then

$$\begin{eqnarray}&&{G}_{\alpha }={E}_{j}-{E}_{i}=\displaystyle \sum _{k=i}^{j-1}{s}_{k},\end{eqnarray} \tag{ D.8 }$$

where we have assumed i < j. For i ≪ j we can apply the central limit theorem and obtain that the gaps are random variables

$$\begin{eqnarray}&&{G}_{\alpha }\sim N({\mu }_{s}(j-i),\sqrt{{\mu }_{s}(j-i)}),\end{eqnarray} \tag{ D.9 }$$

where N(μ, σ) is the Gaussian distribution with mean μ and standard deviation σ. The difference of gaps generated by different probability distributions

$$\begin{eqnarray}&&\delta {G}_{\alpha }:= {G}_{\alpha }^{(1)}-{G}_{\alpha }^{(2)},\sim N(0,\sqrt{2(j-i){\mu }_{s}}).\end{eqnarray} \tag{ D.10 }$$

Note the standard deviation increasing with the separation of the energy levels. This is indeed the behaviour of δG_α if both s_k⁽¹⁾ and s_k⁽²⁾ are independent random variables. However, in our example the spectrum E_k⁽²⁾ is built from E_k⁽¹⁾ with the single goal of changing its level-statistics. This can be made shifting the energy levels by an amount independent of the separation between them. One way to achieve so is to group the energy levels $\{{E}_{k}^{(1)}\}{}_{k=1}^{d}$ in different sets $\{{E}_{k}^{(1)}\}{}_{k=j\,L}^{(j+1)L-1}$ of consecutive L energy levels, where j labels the different sets. Then, the first and the last energy levels of every set are kept fixed and the other energy levels are shifted according to the other level-statistics. In such a case, the random variable δG_α is bounded by

$$\begin{eqnarray}&&\delta {G}_{\alpha }\leqslant K\sqrt{{\mu }_{S}}\end{eqnarray} \tag{ D.11 }$$

with $K=\sqrt{2L}$ a constant independent of the system size.

One strategy to bound the one point distance between the coarse-grained frequency signals would be to use the triangular inequality as follows

$$\begin{eqnarray}&&| {g}_{\epsilon }^{(1)}(\omega )-{g}_{\epsilon }^{(2)}(\omega )| \leqslant \displaystyle \sum _{\alpha }| {v}_{\alpha }| | {u}_{\delta {G}_{\alpha }}(\omega -{\bar{G}}_{\alpha })| \end{eqnarray} \tag{ D.12 }$$

and the one norm distance

$$\begin{eqnarray}&&{\parallel {g}_{\epsilon }^{(1)}-{g}_{\epsilon }^{(2)}\parallel }_{1}\leqslant \displaystyle \sum _{\alpha }| {v}_{\alpha }| \int {\rm{d}}\omega | {u}_{\delta {G}_{\alpha }}(\omega -{\bar{G}}_{\alpha })| .\end{eqnarray} \tag{ D.13 }$$

In order to estimate the scaling of the $| {v}_{\alpha }|$ in the system size, let us consider that there are d_eff energy levels with a non-negligible occupation in the initial state. Thus, let us note that

$$\begin{eqnarray}&&\displaystyle \sum _{\alpha }| {v}_{\alpha }{| }^{2}\sim {d}_{\mathrm{eff}}^{2}| {v}_{\alpha }{| }^{2}\leqslant \displaystyle \frac{1}{{d}_{\mathrm{eff}}},\end{eqnarray} \tag{ D.14 }$$

and therefore $| {v}_{\alpha }| \sim {d}_{\mathrm{eff}}^{-3/2}$.

By splitting the integrals and a convenient change of variables,

$$\begin{eqnarray}\displaystyle \int {\rm{d}}\omega | {u}_{\sigma }(\omega )| & = & 2\left|{\displaystyle \int }_{-\infty }^{\tfrac{\sigma }{2}}{\rm{d}}\omega {h}_{\epsilon }(\omega )-{\displaystyle \int }_{-\infty }^{-\tfrac{\sigma }{2}}{\rm{d}}\omega {h}_{\epsilon }(\omega )\right|\\ & = & 2\left|{\displaystyle \int }_{-\tfrac{\sigma }{2}}^{\tfrac{\sigma }{2}}{\rm{d}}\omega {h}_{\epsilon }(\omega )\right|\leqslant \displaystyle \frac{2\sigma }{\epsilon }.\end{eqnarray} \tag{ D.15 }$$

Note that for $\delta {G}_{\alpha }\leqslant K\sqrt{{\mu }_{s}}$ with μ_s ∼ n / d_eff, the individual terms of the sum behave as $K\sqrt{n}/(\epsilon {d}_{\mathrm{eff}}^{2})$. However, the sum over α contains d_eff² many terms, and therefore the bound becomes useless. The reason for that is that, as it happened for the time signal, the quantity $| {g}^{(1)}(\omega )-{g}^{(2)}(\omega )|$ is small due to the cancellations between its different contributions ${u}_{\delta {G}_{\alpha }}(\omega )$.

Instead, let us exploit interference by applying the triangular inequality as follows

$$\begin{eqnarray}&&| {g}^{(1)}(\omega )-{g}^{(2)}(\omega )| \leqslant \displaystyle \sum _{k}\left|\displaystyle \sum _{{\bar{G}}_{\alpha }\in {I}_{k}}{v}_{\alpha }{u}_{\delta {G}_{\alpha }}(\omega -{\bar{G}}_{\alpha })\right|,\end{eqnarray} \tag{ D.16 }$$

where ${I}_{k}:= [k\delta \omega ,(k+1)\delta \omega )$.

For given k, the contribution to the one-norm behaves as

$$\begin{eqnarray}&&\int {\rm{d}}\omega \left|\displaystyle \sum _{{\bar{G}}_{\alpha }\in {I}_{k}}{v}_{\alpha }{u}_{\delta {G}_{\alpha }}(\omega -{\bar{G}}_{\alpha })\right|\sim \sqrt{{d}_{k}}{\langle {v}_{\alpha }\rangle }_{k}\int {\rm{d}}\omega | {u}_{c\sqrt{{\mu }_{s}}}(\omega )| \leqslant \displaystyle \frac{K}{\epsilon }\sqrt{\displaystyle \frac{{d}_{k}{\mu }_{s}}{{d}_{\mathrm{eff}}^{3}}},\end{eqnarray} \tag{ D.17 }$$

where ${d}_{k}:= | \{{\bar{G}}_{\alpha }\in {I}_{k}\}|$ is the number of gaps in the interval I_k and we have used $| {v}_{\alpha }| \sim {d}_{\mathrm{eff}}^{-3/2}$ and equation (D.15). Note that in (D.17) and due to interference between different ${u}_{\delta {G}_{\alpha }}(\omega -{\bar{G}}_{\alpha })$ the scaling with the number of gaps within an interval I_k is $\sqrt{{d}_{k}}$ and not d_k. This behaviour has been inferred from numerical simulations and can be understood by means of the central limit theorem.

Finally, let us note that the amount of gaps in an interval I_k scales as d_k ∼ C_k d_eff² where d_eff² is the total number of gaps, n the the system size, and C_k the fraction of gaps in the interval I_k. Thus, adding all the contributions of every interval

$$\begin{eqnarray}&&{\parallel {g}_{\epsilon }^{(1)}-{g}_{\epsilon }^{(2)}\parallel }_{1}\leqslant \displaystyle \frac{C}{\epsilon }\sqrt{\displaystyle \frac{{n}^{3}}{{d}_{\mathrm{eff}}}},\end{eqnarray} \tag{ D.18 }$$

where we have used that the average level spacing μ_s scales as μ_s ∼ K₂n / d_eff with C and K₂ positive constants independent of the system size.