Equilibration of quantum gases

Over the past few decades, there has been a major push to understand statistical physics by applying tools from quantum information. One particularly pressing problem is understanding equilibration. From everyday experience, we know it to be universal, as anything from a hot cup of tea to a spinning top will relax to a steady state eventually. See figure 1. However, our understanding of why equilibration occurs and how long it takes remains incomplete. Progress has been dramatically helped by recent advances in experiments [1, 2], where mesoscopic quantum systems can now be controlled extremely well, providing better and better playgrounds to probe properties of many-body systems.

Zoom In Zoom Out Reset image size

In [3–5] it was proved that quantum systems will generally equilibrate with very weak assumptions on the Hamiltonian (which ensure, for one thing, that the system is not a collection of non-interacting subsystems). But very little is known about the timescale. This is crucial: if a system equilibrates but the timescale is the age of the Universe, we will never actually observe it equilibrating in a lab. Unfortunately, the best general upper bounds on the timescale [6–8] are far too large for even mesoscopic systems. This is a consequence of the generality of the results. Indeed, models were constructed effectively saturating these timescale bounds [8, 9].

Imposing physical constraints on Hamiltonians and measurements has led to more realistic timescales in specific cases. In fact, one of the earliest equilibration results was equilibration timescales for bosons evolving via the Hubbard Hamiltonian in the absence of interactions [10, 11]. More recently, equilibration timescales for small subsystems interacting with a large thermal bath were found [12]. Along different lines, equilibration timescales were bounded by averaging over Hamiltonians, measurements or initial states [8, 13–17]. For a review, see [18].

Here we will look at N particle systems in the regime of negligible interactions to see when equilibration occurs. Experimentally, such situations appear often: Luttinger liquids [19] are one example.

We look at two classes of measurement, which are natural for macroscopic and mesoscopic systems. The first are coarse-grained measurements. These include the number of particles in some spatial region, the magnetization of fermions on a lattice, or the number of particles with different values of momentum. The last of these arose in experiments with trapped Bose gases [20], which, in the limit of strong point-like interactions, behave like free fermions. The second type of measurements we consider are few-mode measurements. Such measurements are crucial in many settings, and include correlation functions and phase correlators, which are important in ultracold atom experiments.

First we will look at some examples and then we will show that equilibration of N particle systems in this setting occurs quite generally and appears to be much faster than what general timescale bounds suggest.

Because there are recurrences for quantum systems with discrete spectra [21, 22], the naive definition of equilibration as simply relaxation to a steady state is not sufficient. Instead, we say a system equilibrates if it evolves towards a fixed state and stays close to it for most times. To define what it means for two states to be close, we need a definition of distance between states. For this to be realistic, we need to consider what measurements we can actually do. For example, if we can do any measurement we want on a quantum system, then the distance between two states is best quantified by the trace distance, which allows us to calculate the maximum probability of distinguishing two states by doing a measurement [5, 23].

In reality, for systems beyond a few qubits, there will be restrictions on the measurements we can do; for 10²³ particles, clearly we are restricted to very coarse measurements. With this in mind, a useful measure of distance is given by the distinguishability between states ρ and σ, which is defined to be [5]

$$\begin{eqnarray}&&{D}_{{ \mathcal M }}(\rho ,\sigma )=\displaystyle \frac{1}{2}\underset{\{{M}_{i}\}\in { \mathcal M }}{\mathrm{max}}\displaystyle \sum _{i}| \mathrm{tr}[\rho {M}_{i}]-\mathrm{tr}[\sigma {M}_{i}]| ,\end{eqnarray} \tag{ 1 }$$

where ${ \mathcal M }$ denotes the set of measurements we can do, and $\{{M}_{i}\}$ denotes a positive operator valued measure (POVM)measurement, with the positive operators M_i satisfying ${\sum }_{i}{M}_{i}={\mathbb{1}}$. POVM measurements are more general than projective measurements. This description may be necessary in situations where the measurement is not repeatable, for example. Nevertheless, a POVM measurement is equivalent to a projective measurement on the system together with an ancilla [23].

We denote the infinite-time average of $\rho (t)$ by $\langle \rho \rangle$. If ${D}_{{ \mathcal M }}(\rho (t),\langle \rho \rangle )$ is small most of the time, then for all practical purposes $\rho (t)$ is indistinguishable from its time average $\langle \rho \rangle$ most of the time. In that case, equilibration has occurred.

Another notion of equilibration is equilibration of expectation values [3]. This works as follows. Suppose we have the observable M and we look at the quantity

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{M}(t)=\displaystyle \frac{| \mathrm{tr}[\rho (t)M]-\mathrm{tr}[\langle \rho \rangle M]| }{\parallel M\parallel },\end{eqnarray} \tag{ 2 }$$

where $\parallel M\parallel$ is the operator norm of M. This quantity tells us how close the expectation value of M at time t is to its time average, with the scale set by $\parallel M\parallel$.

If equilibration is to occur, we require that most of the time ${{\rm{\Delta }}}_{M}(t)$ is smaller than some , with chosen so that $\epsilon \parallel M\parallel$ is smaller than our experimental resolution.

There is an important caveat here. Even if expectation values equilibrate, we do not measure expectation values; we measure POVM outcomes. In the examples we consider where equilibration of expectation values occurs, the fluctuations in measurement results are unobservably small. This means that the measured value of M is experimentally indistinguishable from $\mathrm{tr}[\rho (t)M]$ with extremely high probability. Therefore, equilibration truly occurs.

The key step in getting estimates of the equilibration time for N particle systems is equation (6) below, which will allow us to equate ${{\rm{\Delta }}}_{M}(t)$ to the distinguishability for a single particle.

First, it will be useful to introduce some notation. Let ${ \mathcal H }$ be a single-particle Hilbert space, and let $| i\rangle$ denote an orthonormal basis. Then we can define creation operators ${a}_{i}^{\dagger }$, acting on a fermionic Hilbert space, that create fermions corresponding to these states. Equivalently, we may say ${a}_{i}^{\dagger }$ creates a fermion in mode i. The fermionic Hilbert space is spanned by states with varying numbers of creation operators acting on $| 0\rangle$, the empty state. To avoid confusion, any state vectors written as kets are in the single-particle Hilbert space ${ \mathcal H }$, with the exception of $| 0\rangle$, which represents the empty state in a fermionic (or bosonic) system.

The creation operator that creates a particle corresponding to the single-particle state $| \psi \rangle ={\sum }_{i}{c}_{i}| i\rangle$ is ${a}^{\dagger }(| \psi \rangle )={\sum }_{i}{c}_{i}{a}_{i}^{\dagger }$. Suppose we have a single-particle Hamiltonian with discrete spectrum

$$\begin{eqnarray}&&H=\displaystyle \sum _{E}E\;| E\rangle \langle E| ,\end{eqnarray} \tag{ 3 }$$

where E labels the energies. There is a corresponding fermionic Hamiltonian, given by

$$\begin{eqnarray}&&{H}_{{\rm{f}}}=\displaystyle \sum _{E}E\;{a}^{\dagger }(| E\rangle )a(| E\rangle ).\end{eqnarray} \tag{ 4 }$$

For any single-particle state $| \psi \rangle$, we also have

$$\begin{eqnarray}&&{{\rm{e}}}^{-{\rm{i}}{H}_{{\rm{f}}}t}{a}^{\dagger }(| \psi \rangle ){{\rm{e}}}^{{\rm{i}}{H}_{{\rm{f}}}t}={a}^{\dagger }({{\rm{e}}}^{-{\rm{i}}{Ht}}| \psi \rangle )={a}^{\dagger }(| \psi (t)\rangle ).\end{eqnarray} \tag{ 5 }$$

The situation for bosons is similar. The only difference is that, while fermionic creation and annihilation operators obey the canonical anti-commutation relations, bosonic creation and annihilation operators obey the canonical commutation relations.

This is the basic idea behind second quantization, which allows one to take a single-particle system and upgrade it to a multi particle system [24]. Our goal here is to go in the opposite direction and to study equilibration of many-particle systems by moving to the single-particle picture. Let us now give a useful simplification for free bosons or fermions.

Theorem 1. Take a state $\rho (t)=U(t)\rho (0){U}^{\dagger }(t)$ of $N$ non-interacting bosons or fermions and a measurement operator counting the number of particles in some orthogonal modes $M={\sum }_{i}{b}_{i}^{\dagger }{b}_{i}$, where ${b}_{i}=a(| {\phi }_{i}\rangle )$. Then, there exist orthonormal single-particle states $| {\psi }_{\alpha }(t)\rangle$, evolving via the corresponding single-particle Hamiltonian, such that

$$\begin{eqnarray}&&\mathrm{tr}[\rho (t)M]=\displaystyle \sum _{\alpha }{n}_{\alpha }\mathrm{tr}[{\psi }_{\alpha }(t)P],\end{eqnarray} \tag{ 6 }$$

where ${n}_{\alpha }$ are occupation numbers adding up to N, $P={\sum }_{i}| {\phi }_{i}\rangle \langle {\phi }_{i}|$ and ${\psi }_{\alpha }(t)=| {\psi }_{\alpha }(t)\rangle \langle {\psi }_{\alpha }(t)|$.

This is proved for any N particle state in appendix A. Here we will just prove it for the simpler case of an initial state with N bosons or fermions ${a}_{1}^{\dagger {n}_{1}}...{a}_{{n}_{k}}^{\dagger {n}_{k}}| 0\rangle$, where ${a}_{\alpha }^{\dagger }={a}^{\dagger }(| {\psi }_{\alpha }\rangle )$, and $| {\psi }_{\alpha }\rangle$ are some orthonormal single-particle states. In this case, the states $| {\psi }_{\alpha }\rangle$ mentioned in the theorem are already given.

Proof. Expand

$$\begin{eqnarray}&&{a}^{\dagger }(| {\phi }_{i}\rangle )={b}_{i}^{\dagger }=\displaystyle \sum _{\alpha }{c}_{i,\alpha }{a}_{\alpha }^{\dagger },\end{eqnarray} \tag{ 7 }$$

where ${c}_{i,\alpha }$ are complex numbers. Then

$$\begin{eqnarray}\mathrm{tr}[\rho (0){b}_{i}^{\dagger }{b}_{i}] & = & \displaystyle \sum _{\alpha }| {c}_{i,\alpha }{| }^{2}\mathrm{tr}[\rho (0){a}_{\alpha }^{\dagger }{a}_{\alpha }]\\ & = & \displaystyle \sum _{\alpha }{n}_{\alpha }| {c}_{i,\alpha }{| }^{2},\end{eqnarray} \tag{ 8 }$$

where ${n}_{\alpha }$ is the number of particles in mode α. Next, we use ${c}_{i,\alpha }=\{{a}_{\alpha },{b}_{i}^{\dagger }\}=\langle {\psi }_{\alpha }| {\phi }_{i}\rangle$ for fermions, or ${c}_{i,\alpha }=[{a}_{\alpha },{b}_{i}^{\dagger }]=\langle {\psi }_{\alpha }| {\phi }_{i}\rangle$ for bosons, to get

$$\begin{eqnarray}&&\mathrm{tr}[\rho (0){b}_{i}^{\dagger }{b}_{i}]=\displaystyle \sum _{\alpha }{n}_{\alpha }\langle {\psi }_{\alpha }| {\phi }_{i}\rangle \langle {\phi }_{i}| {\psi }_{\alpha }\rangle .\end{eqnarray} \tag{ 9 }$$

Therefore

$$\begin{eqnarray}\mathrm{tr}[\rho (0)M] & = & \displaystyle \sum _{\alpha }{n}_{\alpha }\langle {\psi }_{\alpha }| P| {\psi }_{\alpha }\rangle \\ & = & \displaystyle \sum _{\alpha }{n}_{\alpha }\mathrm{tr}[{\psi }_{\alpha }P],\end{eqnarray} \tag{ 10 }$$

where $P={\sum }_{i}| {\phi }_{i}\rangle \langle {\phi }_{i}|$. To incorporate the dependence on time, we use ${a}_{\alpha }(t)=U(t){a}_{\alpha }{U}^{\dagger }(t)$ and

$$\begin{eqnarray}&&\mathrm{tr}[\rho (0){a}_{\alpha }^{\dagger }{a}_{\beta }]=\mathrm{tr}[\rho (t){a}_{\alpha }^{\dagger }(t){a}_{\beta }(t)].\end{eqnarray} \tag{ 11 }$$

The end result is

$$\begin{eqnarray}&&\mathrm{tr}[\rho (t)M]=\displaystyle \sum _{\alpha }{n}_{\alpha }\mathrm{tr}[{\psi }_{\alpha }(t)P].\end{eqnarray} \tag{ 12 }$$

□

Notice that linearity of the time average, together with equation (6) implies

$$\begin{eqnarray}&&\mathrm{tr}[\langle \rho \rangle M]=\displaystyle \sum _{\alpha }{n}_{\alpha }\mathrm{tr}[\langle {\psi }_{\alpha }\rangle P].\end{eqnarray} \tag{ 13 }$$

3.1. Coarse-grained measurements

We can apply this to ${{\rm{\Delta }}}_{M}(t)$, noting that for the applications we are interested in $\parallel M\parallel =N$ when restricted to the N particle subspace. This occurs, for example, when we are measuring the particle number in a region of space. Put another way, we take the experimental accuracy of our measurements to be at best $\epsilon N$, where $\epsilon$ is some very small constant. For equilibration to occur, we need ${{\rm{\Delta }}}_{M}(t)$ to be small compared to most of the time. We get

$$\begin{eqnarray}\displaystyle \frac{| \mathrm{tr}[\rho (t)M]-\mathrm{tr}[\langle \rho \rangle M]| }{\parallel M\parallel } & = & | \mathrm{tr}[\sigma (t)P]-\mathrm{tr}[\langle \sigma \rangle P]| \\ & = & {D}_{P}(\sigma (t),\langle \sigma \rangle ),\end{eqnarray} \tag{ 14 }$$

where $\sigma (t)=\tfrac{1}{N}{\sum }_{j}| {\psi }_{j}(t)\rangle \langle {\psi }_{j}(t)|$ is a single-particle state. In words, the N particle problem has been replaced by a single-particle problem in terms of the distinguishability given a single measurement with projectors P and ${\mathbb{1}}-P$.

Now recall that equilibration of expectation values does not necessarily imply that equilibration will be observed. For the examples we look at, the fluctuations in the observed value of M, given by ${(\mathrm{tr}[\rho (t){M}^{2}]-\mathrm{tr}{[\rho (t)M]}^{2})}^{1/2}$, are bounded above by $\sqrt{N}$, which is proved in appendix B. In fact, a large class of fermion systems have time-averaged fluctuations bounded above by $\sqrt{N}$, as seen in appendix B. For large numbers of particles, (comparable to 10²³, for example) $\sqrt{N}$ is small compared to our experimental precision $\epsilon N$, and the fluctuations are not practically observable. Even for dilute gases with $O({10}^{4})$ particles, $\sqrt{N}\sim 100$, so the fluctuations are of the order of 1% of the total particle number, which is still quite small.

3.2. Few-mode measurements

Suppose we have a one-dimensional box with a partition at the halfway point (this can be extended to a three-dimensional example as shown in appendix C). On the left of the partition we have N fermions or bosons at zero temperature. We open the partition at t = 0, and the observable we focus on is M, which counts the particles in the left half of the box.

Using equation (14), we can replace this N particle problem by a single-particle one, so

$$\begin{eqnarray}&&\left|\displaystyle \frac{\mathrm{tr}[\rho (t)M]-\mathrm{tr}[\langle \rho \rangle M]}{N}\right|={D}_{P}(\sigma (t),\langle \sigma \rangle ),\end{eqnarray} \tag{ 15 }$$

which is plotted in figure 2. Here $\sigma (t)$ is a state of a free particle in a box, $\langle \sigma \rangle$ is its time average and P is the projector onto the left-hand side of the box.

Zoom In Zoom Out Reset image size

4.1. Fermions

First, let us look at the case where the particles are fermions. The initial state of the N fermion system has all fermions in the left half of the box at temperature zero. This means that the initial single-particle state $\sigma (0)$ is an equal mixture of the lowest N energy levels of a particle trapped in the left half of a box. This can be seen from equation (6).

The energy eigenstates for a particle in a box are given by

$$\begin{eqnarray}&&| n\rangle ={\int }_{0}^{L}{\rm{d}}x\sqrt{\displaystyle \frac{2}{L}}\mathrm{sin}\left(\displaystyle \frac{\pi {nx}}{L}\right)| x\rangle ,\end{eqnarray} \tag{ 16 }$$

where $n\gt 0$ is an integer and L is the length of the box. Similarly, the energy eigenstates for a particle trapped in the left half of the box are given by

$$\begin{eqnarray}&&| {\psi }_{k}\rangle ={\int }_{0}^{L/2}{\rm{d}}x\sqrt{\displaystyle \frac{4}{L}}\mathrm{sin}\left(\displaystyle \frac{2\pi {kx}}{L}\right)| x\rangle ,\end{eqnarray} \tag{ 17 }$$

where again $k\gt 0$ is an integer.

The initial state of the single-particle system is

$$\begin{eqnarray}&&\sigma (0)=\displaystyle \frac{1}{N}\displaystyle \sum _{k=1}^{N}| {\psi }_{k}\rangle \langle {\psi }_{k}| ,\end{eqnarray} \tag{ 18 }$$

with matrix elements ${\sigma }_{{nm}}=\langle n| \sigma (0)| m\rangle$. Similarly, the matrix elements of the projector onto the left half of the box are ${P}_{{mn}}=\langle m| P| n\rangle$.

Let us look at the distinguishability to see if the system equilibrates. In figure 2, the distinguishability as a function of time is plotted. From the plots we can see that, as the number of fermions increases, the average distinguishability gets smaller. Notice that particles in a box have an exact recurrence time of ${T}_{\mathrm{rec}}=2\pi /\nu$ since the energy levels are ${E}_{n}=\nu {n}^{2}$, where n is an integer greater than zero, and $\nu ={\pi }^{2}/2{{mL}}^{2}$. This is because all phases of density matrix elements in the energy basis ${{\rm{e}}}^{-{\rm{i}}({E}_{n}-{E}_{m})t}$ are 1 at $t=2\pi /\nu$. As in [26], this means that we need only study the system over the interval $[0,{T}_{\mathrm{rec}}]$. In fact, with the particular initial state below, a recurrence actually occurs at ${T}_{\mathrm{rec}}/2$.

Evaluating the distinguishability at time t, we get

$$\begin{eqnarray}{D}_{P}(\sigma (t),\langle \sigma \rangle ) & = & \left|\displaystyle \sum _{n\ne m}{{\rm{e}}}^{-{\rm{i}}({n}^{2}-{m}^{2})\nu t}{\sigma }_{{nm}}{P}_{{mn}}\right|\\ & = & \left|\displaystyle \sum _{n\ne m}\mathrm{cos}[({n}^{2}-{m}^{2})\nu t]{\sigma }_{{nm}}{P}_{{mn}}\right|,\end{eqnarray} \tag{ 19 }$$

where we used the fact that ${\sigma }_{{nm}}$ and P_mn are symmetric under swapping n and m because all vectors and operators here are real.

In appendix C, we see that the distinguishability can be written as

$$\begin{eqnarray}&&{D}_{P}(\sigma (t),\langle \sigma \rangle )=\displaystyle \frac{2}{N}\left|\displaystyle \sum _{n\ \mathrm{odd}}\displaystyle \sum _{k=1}^{N}\mathrm{cos}[({n}^{2}-4{k}^{2})\nu t]f{(n,2k)}^{2}\right|,\end{eqnarray} \tag{ 20 }$$

where

$$\begin{eqnarray}&&f{(n,2k)}^{2}=\displaystyle \frac{4}{{\pi }^{2}}\displaystyle \frac{4{k}^{2}}{{({n}^{2}-4{k}^{2})}^{2}},\end{eqnarray} \tag{ 21 }$$

for $n\ne 2k$.

For the system to equilibrate, we need it to spend most of its time indistinguishable from its time-average state. We see in appendix C, that the time-average distinguishability satisfies $\langle {D}_{P}(\sigma (t),\langle \sigma \rangle )\rangle =O(\mathrm{ln}{(N)}^{2}/N)$. Therefore, most of the time the system state is indistinguishable from its time average, provided N is large.

We can also say something about the timescale. We see in appendix C, that the timescale for equilibration is at most

$$\begin{eqnarray}&&{T}_{\mathrm{eq}}=\displaystyle \frac{1}{{Na}\nu }=\displaystyle \frac{2{{mL}}^{2}}{{Na}{\pi }^{2}}.\end{eqnarray} \tag{ 22 }$$

Here a is a small constant that we choose such that the distinguishability at $t={T}_{\mathrm{eq}}$ is small:

$$\begin{eqnarray}&&{D}_{P}(\sigma ({T}_{\mathrm{eq}}),\langle \sigma \rangle )\leqslant \displaystyle \frac{\pi a}{3}+O\left(\displaystyle \frac{\mathrm{log}{(N)}^{2}}{N}\right),\end{eqnarray} \tag{ 23 }$$

which is also derived in appendix C. Interestingly, the timescale decreases with increasing particle number.

4.2. Bosons

For our second example, suppose we have N bosons at zero temperature in a one-dimensional harmonic trap with frequency ${\omega }_{0}$. We will consider what happens after two different quenches.

5.1. Quench to a square well potential

Suppose the Hamiltonian changes suddenly so that the bosons are then confined in a deep square well potential, which we can idealize as a box corresponding to the interval $\left[-\tfrac{L}{2},\tfrac{L}{2}\right]$. Let the measurement operator M count the number of bosons in the central region of the box $\left[-\tfrac{L}{4},\tfrac{L}{4}\right]$. Applying equation (14), we see that

$$\begin{eqnarray}&&\left|\displaystyle \frac{\mathrm{tr}[\rho (t)M]-\mathrm{tr}[\langle \rho \rangle M]}{N}\right|={D}_{P}(\psi (t),\langle \psi \rangle ),\end{eqnarray} \tag{ 24 }$$

where $\psi (t)=| \psi (t)\rangle \langle \psi (t)|$ is a pure state of a single-particle and P is the projector onto the central region of the box.

The equilibration timescale has already been estimated for this single-particle system in [26]. First, the infinite-time average of ${D}_{P}(\psi (t),\langle \psi \rangle )$ is numerically shown to scale like

$$\begin{eqnarray}&&\langle {D}_{P}(\psi (t),\langle \psi \rangle )\rangle \sim {\left(\displaystyle \frac{1}{L}\sqrt{\displaystyle \frac{8\pi }{m{\omega }_{0}}}\right)}^{0.79},\end{eqnarray} \tag{ 25 }$$

where it was assumed that the width of the initial wavefunction is small compared to the length of the box, meaning $\sigma =1/\sqrt{2m{\omega }_{0}}\ll L$. So for sufficiently narrow potentials (or sufficiently large boxes), equilibration occurs. Furthermore, the timescale for equilibration is shown to be [26]

$$\begin{eqnarray}&&{T}_{\mathrm{eq}}=\displaystyle \frac{L}{\pi }\sqrt{\displaystyle \frac{m}{2{\omega }_{0}}}.\end{eqnarray} \tag{ 26 }$$

It would be interesting to observe this experimentally. In fact, it may be feasible to create square-well potentials in practice: in [27] a three-dimensional cylindrical potential was created to trap a Bose–Einstein condensate, so creating potentials with sharply defined walls may be possible.

5.2. Quench to a weaker harmonic trap

Recent experiments have followed the dynamics of Bose gases after a different quench to that of the previous section. By quickly changing the strength of a harmonic trap, oscillatory behavior was observed [28]. Such behavior occurred in both the strongly and weakly interacting regimes. For our purposes, the latter of these regimes is relevant and corresponds to an ideal Bose gas in one-dimension. In [28] the ratio of initial trap frequency ${\omega }_{0}$ and post-quench frequency ω was close to one: ${\omega }_{0}/\omega \simeq 1.3$. Here we see equilibration when ${\omega }_{0}/\omega$ is much larger than one.

For our observable, let us take the number of bosons in the spatial region $[-l,l]$. Again, using equation (14), we can replace this N particle problem by a single-particle one, so

$$\begin{eqnarray}&&\left|\displaystyle \frac{\mathrm{tr}[\rho (t)M]-\mathrm{tr}[\langle \rho \rangle M]}{N}\right|={D}_{P}(\psi (t),\langle \psi \rangle ),\end{eqnarray} \tag{ 27 }$$

where the P is the projector onto the region $[-l,l]$. The distinguishability is

$$\begin{eqnarray}&&{D}_{\psi }(\psi (t),\langle \psi \rangle )=| \mathrm{tr}[P\psi (t)]-\mathrm{tr}[P\langle \psi \rangle ]| .\end{eqnarray} \tag{ 28 }$$

So we need only see if $\mathrm{tr}[P\psi (t)]$ spends most of its time close to its time average.

The problem is simplified by using the propagator for a harmonic oscillator with frequency ω, given by [29]

$$\begin{eqnarray}&&K(x,y,t)=\sqrt{\displaystyle \frac{m\omega }{2\pi {\rm{i}}\mathrm{sin}(\omega t)}}\mathrm{exp}\left[-\displaystyle \frac{m\omega (({x}^{2}+{y}^{2})\mathrm{cos}(\omega t)-2{xy})}{2{\rm{i}}\mathrm{sin}(\omega t)}\right],\end{eqnarray} \tag{ 29 }$$

which leads to the expression

$$\begin{eqnarray}&&\mathrm{tr}[P\psi (t)]=\displaystyle \int {\rm{d}}{y}_{1}{\rm{d}}{y}_{2}{\displaystyle \int }_{-l}^{l}{\rm{d}}x\;{\psi }^{* }({y}_{1}){K}^{* }(x,{y}_{1},t)K(x,{y}_{2},t)\psi ({y}_{2}).\end{eqnarray} \tag{ 30 }$$

As $\psi (x)$ is a Gaussian wavefunction, the y₁ and y₂ integrals are straightforward, leading to

$$\begin{eqnarray}&&\mathrm{tr}[P\psi (t)]=\displaystyle \frac{1}{\sqrt{\pi }}{\int }_{-l\sqrt{\alpha (t)}}^{l\sqrt{\alpha (t)}}{\rm{d}}x\;{{\rm{e}}}^{-{x}^{2}},\end{eqnarray} \tag{ 31 }$$

where

$$\begin{eqnarray}&&\alpha (t)=\displaystyle \frac{m{\omega }_{0}}{{\gamma }^{2}{\mathrm{sin}}^{2}(\omega t)+{\mathrm{cos}}^{2}(\omega t)},\end{eqnarray} \tag{ 32 }$$

with $\gamma ={\omega }_{0}/\omega$. Next we use the approximation for the error function [30]

$$\begin{eqnarray}\mathrm{erf}(x) & = & \displaystyle \frac{2}{\sqrt{\pi }}{\displaystyle \int }_{0}^{x}{\rm{d}}t\;{{\rm{e}}}^{-{t}^{2}}\\ & \simeq & \mathrm{sgn}(x)\sqrt{1-\mathrm{exp}\left[-{x}^{2}\displaystyle \frac{\tfrac{4}{\pi }+{{bx}}^{2}}{1+{{bx}}^{2}}\right]},\end{eqnarray} \tag{ 33 }$$

where the maximum error for any x is around 0.00012, and $b\simeq 0.147$. The result is that

$$\begin{eqnarray}&&\mathrm{tr}[P\psi (t)]\simeq \sqrt{1-\mathrm{exp}\left[-\alpha (t){l}^{2}\displaystyle \frac{\tfrac{4}{\pi }+\alpha (t)b{l}^{2}}{1+\alpha (t)b{l}^{2}}\right]}.\end{eqnarray} \tag{ 34 }$$

Notice that there are four independent parameters that matter: l, which controls the width of the interval the measurement looks at; ω, which is the frequency of the trap after the quench; $\gamma ={\omega }_{0}/\omega$, which is the ratio of trap strengths before and after the quench; and $m{\omega }_{0}$, which determines the width of the initial state. A natural starting point is to choose l so that the initial state is almost entirely contained in $[-l,l]$, so we can fix ${l}^{2}=4/(m{\omega }_{0})$.

As we can see from figure 3, as γ becomes bigger and bigger, $\mathrm{tr}[P\psi (t)]$ spends most of its time close to zero. So for very large γ, equilibration occurs. In fact, we can see directly from equations (32) and (34) that, as γ tends to infinity, $\mathrm{tr}[P\psi (t)]$ tends to zero. This holds for all times, except when $\omega t=n\pi$, with $n\in {\mathbb{Z}}$.

Zoom In Zoom Out Reset image size

In [28], oscillatory behavior was seen at $\gamma =1.3$. Here, this value of γ does not lead to any significant departure from the initial state, as seen in figure 3. The reason for this difference is that in [28] the initial states were at non-zero temperature. Here, we are initially at zero temperature, and we see oscillations at higher values of γ.

To estimate the equilibration time when equilibration does occur, we estimate how long it takes for $\mathrm{tr}[P\psi (t)]$ to reach $p\ll 1$. Using equation (34) and $\mathrm{log}(1-{p}^{2})\simeq -{p}^{2}$, we get

$$\begin{eqnarray}&&\alpha (t){l}^{2}\displaystyle \frac{\tfrac{4}{\pi }+\alpha (t)b\;{l}^{2}}{1+\alpha (t)b\;{l}^{2}}\simeq {p}^{2}.\end{eqnarray} \tag{ 35 }$$

Since p is small, this requires $\alpha (t){l}^{2}$ to be small. Using the earlier choice ${l}^{2}=4/(m{\omega }_{0})$, we get

$$\begin{eqnarray}&&\alpha (t){l}^{2}\displaystyle \frac{4}{\pi }=\displaystyle \frac{16/\pi }{{\gamma }^{2}{\mathrm{sin}}^{2}(\omega t)+{\mathrm{cos}}^{2}(\omega t)}\simeq {p}^{2}.\end{eqnarray} \tag{ 36 }$$

For large γ, this is satisfied at

$$\begin{eqnarray}&&{T}_{\mathrm{eq}}\simeq \displaystyle \frac{4}{\sqrt{\pi }\gamma \omega p}=\displaystyle \frac{4}{\sqrt{\pi }{\omega }_{0}p},\end{eqnarray} \tag{ 37 }$$

where we assumed that t was small compared to $1/\omega$, and used $\mathrm{sin}(x)\simeq x$ for small x.

The examples we looked at were encouraging, but a pressing question is whether one can say anything more general. The answer is actually yes: here we will see general estimates for the equilibration timescale of gases with negligible interactions. The starting point is again to replace the N particle problem by a single-particle problem. Then we can use a single-particle equilibration result, which builds on previous results [6–8].

First, let us state the single-particle equilibration result. We will need to take account of degenerate energy gaps. These occur when two different energy gaps are equal: ${E}_{i}-{E}_{j}={E}_{k}-{E}_{l}$, when ${E}_{i}\ne {E}_{k}$ and ${E}_{i}\ne {E}_{j}$. The degeneracy of the most degenerate gap is denoted by D_G. If all gaps are different, D_G  =  1. For a particle in a box there are some degenerate energy gaps, though the addition of an inhomogeneous potential V(x) would generally change this. The harmonic oscillator also has many degenerate energy gaps. Typically, however, these are very special cases, and we expect most Hamiltonians would have few degenerate energy gaps.

Theorem 2. Suppose we have a single-particle system with a d-dimensional state space. Let $A$ be an operator with operator norm $\parallel A\parallel$, and let $\sigma (t)$ be a state unitarily evolving via a Hamiltonian $H$. Denote the infinite-time average of $\sigma (t)$ by $\langle \sigma \rangle$. Assuming that we can make the density of states approximation, meaning we replace ${\sum }_{E}$ by $\int {\rm{d}}E\;n(E)$, where $n(E)$ is the density of states, we get

$$\begin{eqnarray}&&\displaystyle \frac{{\langle | \mathrm{tr}[\sigma (t)A]-\mathrm{tr}[\langle \sigma \rangle A]{| }^{2}\rangle }_{T}}{{\parallel A\parallel }^{2}}\leqslant \displaystyle \frac{{c}_{1}}{{d}_{\mathrm{eff}}}\left[{D}_{G}+{c}_{2}\displaystyle \frac{{n}_{\mathrm{max}}d}{T}\right],\end{eqnarray} \tag{ 38 }$$

where ${\langle \cdot \rangle }_{T}$ denotes the time average over $[0,T]$, and we have constants ${c}_{1}=e\sqrt{\pi }/2$ and ${c}_{2}=4\sqrt{\pi }$. Also, ${n}_{\mathrm{max}}={\mathrm{max}}_{E}\;n(E)$, and the effective dimension of the state $\sigma (t)$ is defined by

$$\begin{eqnarray}&&\displaystyle \frac{1}{{d}_{\mathrm{eff}}}=\displaystyle \sum _{E}{(\mathrm{tr}[\sigma (0){P}_{E}])}^{2},\end{eqnarray} \tag{ 39 }$$

where ${P}_{E}$ is the projector onto the energy eigenspace corresponding to energy $E$.

Equilibration of the expectation value of A occurs provided the right-hand side of equation (38) is sufficiently small. As $T\to \infty$, equilibration is guaranteed if ${c}_{1}{D}_{G}/{d}_{\mathrm{eff}}\ll 1$. The effective dimension ${d}_{\mathrm{eff}}$ measures how spread out over energy levels the initial state is. If ${d}_{\mathrm{eff}}$ is very large we expect equilibration to occur.

But we can also estimate the timescale: the equilibration timescale can be bounded above by the smallest T such that the right-hand side of equation (38) is small. In other words, when equilibration occurs, we get an upper bound for the timescale:

$$\begin{eqnarray}&&{T}_{\mathrm{eq}}\propto \displaystyle \frac{{n}_{\mathrm{max}}d}{{d}_{\mathrm{eff}}}.\end{eqnarray} \tag{ 40 }$$

The main task now is to use this single-particle equilibration result to find timescales for N particle systems.

6.1. Coarse-grained measurements

Let us start with coarse-grained measurements. We will see that equilibration of coarse-grained observables generally occurs much quicker than what we would expect based on previous timescale bounds from [6, 8].

By mapping an N particle problem to the single-particle picture via equation (14), we want to bound

$$\begin{eqnarray}{\left\langle \left|\displaystyle \frac{\mathrm{tr}[\rho (t)M]-\mathrm{tr}[\langle \rho \rangle M]}{N}\right|\right\rangle }_{T} & = & {\langle {D}_{P}(\sigma (t),\langle \sigma \rangle )\rangle }_{T}\\ & = & {\langle | \mathrm{tr}[\sigma (t)P]-\mathrm{tr}[\langle \sigma \rangle P]| \rangle }_{T}\\ & \leqslant & \sqrt{{\langle {(\mathrm{tr}[\sigma (t)P]-\mathrm{tr}[\langle \sigma \rangle P])}^{2}\rangle }_{T}}\\ & \leqslant & \sqrt{\displaystyle \frac{{c}_{1}}{{d}_{\mathrm{eff}}}\left[{D}_{G}+{c}_{2}\displaystyle \frac{{n}_{\mathrm{max}}d}{T}\right]},\end{eqnarray} \tag{ 41 }$$

where the third line follows from concavity of the square root. The last line follows from the result for single-particle equilibration from the last section, namely equation (38).

To see whether equilibration occurs at all, let $T\to \infty$ to get the infinite-time average. And so we must estimate $\tfrac{1}{{d}_{\mathrm{eff}}}$. Defining N_E to be the operator that counts the number of particles in energy level E, we get

$$\begin{eqnarray}\displaystyle \frac{1}{{d}_{\mathrm{eff}}} & = & \displaystyle \sum _{E}{(\mathrm{tr}[\sigma {P}_{E}])}^{2}\\ & = & \displaystyle \frac{1}{{N}^{2}}\displaystyle \sum _{E}{(\mathrm{tr}[\rho {N}_{E}])}^{2}\\ & = & \displaystyle \sum _{E}{\left(\displaystyle \frac{{n}_{E}}{N}\right)}^{2},\end{eqnarray} \tag{ 42 }$$

where ${n}_{E}=\mathrm{tr}[\rho {N}_{E}]$. Getting the second line used equation (6) from theorem 1. So we see that $\tfrac{1}{{d}_{\mathrm{eff}}}$ is extremely small if the state is spread over many energy levels.

Consider N fermions or the special case of N bosons in orthogonal modes. Then the resulting single-particle density operator $\sigma (t)$ is an equal mixture of N orthogonal states. In that case, $1/{d}_{\mathrm{eff}}\leqslant \tfrac{{n}_{{\rm{d}}}}{{N}^{2}}{\sum }_{E}{n}_{E}\leqslant {n}_{{\rm{d}}}/N$, where n_d is the degeneracy of the most degenerate energy level. As a result

$$\begin{eqnarray}&&\left\langle \left|\displaystyle \frac{\mathrm{tr}[\rho (t)M]-\mathrm{tr}[\langle \rho \rangle M]}{N}\right|\right\rangle \leqslant \sqrt{\displaystyle \frac{{c}_{1}{n}_{{\rm{d}}}{D}_{G}}{N}}.\end{eqnarray} \tag{ 43 }$$

So the bottom line is that for the coarse measurements considered here, equilibration occurs very generally despite the fact that these are non-interacting bosons or fermions.

We can also say something substantial about the equilibration timescale.

We can always restrict our attention to d energy levels of the corresponding single-particle system, which may require an energy cutoff. And suppose d is bounded above by a polynomial in N. This depends on the state $\sigma (t)$ and so ultimately on the state of each of the N bosons or fermions. For example, for the calculations with fermions equilibrating in appendix C, we effectively took a cutoff with $d=O(N)$. In fact, for the bosonic examples, d was independent of N. For lattice systems this is particularly natural if there is a constant density of particles, then $d\propto V\propto N$, where V is the number of lattice sites.

Next, we estimate n_max, which is often polynomial in d, and hence N. For example, ${n}_{\mathrm{max}}\sim {d}^{3}$ for a system whose energy levels go like ${E}_{n}\propto 1/{n}^{2}$, similar to the energies for bound states in a Coulomb potential. Notice that this is a system we would expect to have very many small gaps. Conversely, when the energy level spacings grow with d we would expect better behavior. For example, when ${E}_{n}\propto {n}^{2}$, one gets ${n}_{\mathrm{max}}\sim 1$.

Putting this all together, if equilibration occurs, the timescale is typically

$$\begin{eqnarray}&&{T}_{\mathrm{eq}}\in O({N}^{k})\end{eqnarray} \tag{ 44 }$$

for some positive integer k. This is far better than the bounds of [6–8], which were exponential in N for physical systems. Of course, how n_max scales with d and how d scales with N depend on the system in question, but neither of the requirements above appear unnaturally restrictive.

It is also interesting that each of [15–17] found equilibration timescales that were polynomial (or faster) in the number of particles. In contrast to the setting considered here, these results involved averaging over Hamiltonians with respect to the global unitary Haar measure. Because of this, it is not clear how to interpret the implications of [15–17] for equilibration of local Hamiltonian systems. Nevertheless, [15–17] do say something about equilibration timescales of fully interacting models, which is very interesting.

6.2. Local equilibration

We can also look at equilibration of non-interacting lattice models. This would include the free superfluid regime of the Bose–Hubbard model, for example. We consider local few-mode measurements, where few means that the number of modes is small compared to the number of lattice sites. This setting includes all measurements in some small region of the lattice or correlation functions over long distances. It also includes phase correlators, which are important in ultracold atom systems.

We will state the single-mode result first. This relies on the Hamiltonian being some form of local (not necessarily nearest-neighbour) hopping Hamiltonian: the tight-binding model is one example.

To make the formulas easier to read, we will assume that the maximum energy level degeneracy n_d and the number of modes per site are both one. In the proofs of these results in appendix E we allow other values of these quantities.

Corollary 1. Take a free lattice model, and assume we can make the density of states approximation, as in theorem 2. Let $M={a}^{\dagger }(| \phi \rangle )a(| \phi \rangle )$, where $| \phi \rangle$ is a single-particle state localized on at most $l$ sites (which need not be near each other). Then we have

$$\begin{eqnarray}&&{\langle | \mathrm{tr}[\rho (t)M]-\mathrm{tr}[\langle \rho \rangle M]| \rangle }_{T}\leqslant l\sqrt{{c}_{1}\left[\displaystyle \frac{{D}_{G}}{d}+{c}_{2}\displaystyle \frac{{n}_{\mathrm{max}}}{T}\right]},\end{eqnarray} \tag{ 45 }$$

where $d$ is the dimension of the corresponding single-particle Hilbert space, and we have constants ${c}_{2}=4\sqrt{\pi }$ and ${c}_{1}=(e\sqrt{\pi }/2)$. Also, ${n}_{\mathrm{max}}={\mathrm{max}}_{E}\;n(E)$, where $n(E)$ is the density of states.

For bosons, we needed to assume that the initial state has at most one boson in each mode. Otherwise, the same result holds, but with an extra factor on the right-hand side given by the maximum number of bosons in a given mode.

This is proved in appendix E. Again, we see equilibration provided the right-hand side of equation (45) is small. We will estimate the equilibration timescale below corollary 2. First, let us discuss some consequences of this result.

A simple consequence is that phase correlators equilibrate. Phase correlators are expectation values like $\mathrm{tr}[\rho (t){a}_{\vec{x}}^{\dagger }{a}_{\vec{y}}]$, where $\vec{x}$ and $\vec{y}$ denote lattice sites. (There may be several modes at each lattice site, but for simplicity of notation, we have assumed that there is just one.) Using

$$\begin{eqnarray}&&{a}_{\vec{x}}^{\dagger }{a}_{\vec{y}}=\displaystyle \frac{1}{2}({d}_{1}^{\dagger }{d}_{1}-{d}_{2}^{\dagger }{d}_{2}-{\rm{i}}{d}_{3}^{\dagger }{d}_{3}+{\rm{i}}{d}_{4}^{\dagger }{d}_{4}),\end{eqnarray} \tag{ 46 }$$

where

$$\begin{eqnarray}{d}_{1} & = & \tfrac{1}{\sqrt{2}}({a}_{\vec{x}}+{a}_{\vec{y}})\\ {d}_{2} & = & \tfrac{1}{\sqrt{2}}({a}_{\vec{x}}-{a}_{\vec{y}})\\ {d}_{3} & = & \tfrac{1}{\sqrt{2}}({a}_{\vec{x}}+{\rm{i}}{a}_{\vec{y}})\\ {d}_{4} & = & \tfrac{1}{\sqrt{2}}({a}_{\vec{x}}-{\rm{i}}{a}_{\vec{y}}),\end{eqnarray} \tag{ 47 }$$

we can express $\mathrm{tr}[\rho (t){a}_{\vec{x}}^{\dagger }{a}_{\vec{y}}]$ in terms of single-mode densities. And so via the triangle inequality, we can upper bound the time average of $| \mathrm{tr}[\rho (t){a}_{\vec{x}}^{\dagger }{a}_{\vec{y}}]-\mathrm{tr}[\rho (t){a}_{\vec{x}}^{\dagger }{a}_{\vec{y}}]|$ using corollary 1.

Interestingly, these results apply to a vast range of initial states $\rho (0)$. This means that one could perform a huge variety of quenches to a free lattice system, and the equilibration results here and timescale bounds (which we will discuss below) apply.

Before discussing timescales, we can build on corollary 1 further, getting the corollary below, which is proved in appendix E.2. We only prove the fermionic result, as the bosonic result is essentially the same.

Corollary 2. Take a free lattice model, and let $M$ be an operator on $l$ sites. Suppose the initial state $\rho (0)$ is Gaussian and satisfies $[\rho (0),N]=0$, where $N$ is the total number operator. (This is still quite general, though it rules out BCS states, for example.) Then we get

$$\begin{eqnarray}&&{\langle | \mathrm{tr}[\rho (t)M]-\mathrm{tr}[\langle \rho \rangle M]| \rangle }_{T}\leqslant {2}^{l+2}{{ml}}^{2}\sqrt{{c}_{1}\left[\displaystyle \frac{{D}_{G}}{d}+{c}_{2}\displaystyle \frac{{n}_{\mathrm{max}}}{T}\right]},\end{eqnarray} \tag{ 48 }$$

where d is the dimension of the corresponding single-particle Hilbert space, and we have constants ${c}_{2}=4\sqrt{\pi }$ and ${c}_{1}=(e\sqrt{\pi }/2)$. Also, ${n}_{\mathrm{max}}={\mathrm{max}}_{E}\;n(E)$, where $n(E)$ is the density of states. Finally, $m$ is the maximum coefficient of $M$ when $M$ is expanded in an operator basis of Majorana fermion operators.

Typically m will be order one, which is the case for correlation functions, for example. Therefore, as long as the number of lattice sites that the measurement acts on is quite small, equilibration will also occur for free lattice systems.

Furthermore, we can use these results to upper bound the equilibration timescale. From corollary 1 and 2, the upper bound for the equilibration timescale scales like ${T}_{\mathrm{eq}}\propto {n}_{\mathrm{max}}$. So it remains to estimate n_max. In appendix E.1, we show that for these lattice models, we can effectively take ${n}_{\mathrm{max}}\propto V$, where V is the number of lattice sites. Therefore, we get

$$\begin{eqnarray}&&{T}_{\mathrm{eq}}\propto V.\end{eqnarray} \tag{ 49 }$$

In particular, for one-dimensional systems, we get ${T}_{\mathrm{eq}}\propto L$, where L is the length of the system.

The scaling with system size is quite significant. Previous bounds [6, 8] were exponential in the system size, whereas here we get something linear. Furthermore, in the one-dimensional case, the scaling is optimal. This can be seen from Lieb-Robinson bounds [31], which imply that the time it takes for information to propagate appreciably from one region to another increases linearly with the distance between the regions. So in one-dimension ${T}_{\mathrm{eq}}\propto L$ is the best we can hope for. The only possibility for better scaling is if one restricts the set of initial states under consideration. A good example of such results for special states appeared in [10].

Finding the timescale involved in equilibration is an important problem in physics, especially in light of recent advances in experiments with mesoscopic quantum systems [1, 2]. The timescale results here required us to restrict our attention to a subclass of measurements, which are physically sensible for macroscopic or mesoscopic systems. We focused on the regime of negligible interactions, which includes Luttinger liquids and the Hubbard model in the free superfluid regime. First, we found example equilibration timescale bounds for gases of bosons and fermions. We also saw that equilibration occurs quite generally in this setting of very weak interactions and is very fast compared to the best known general bounds on the equilibration time.

From here the outlook is promising: a natural next step is to extend these results to quasi-free systems, where the Hamiltonian is quadratic in terms of creation and annihilation operators but does not conserve particle number. Such models arise in the theory of superconductivity. Other options are to extend the results to interacting models via perturbation theory or to look at equilibration in terms of fermionic or bosonic generating functions [32].

The author is very grateful to Tobias Osborne, Mathis Friesdorf, Jens Eisert, Marcel Goihl, David Reeb and Shane Dooley for helpful discussions and comments. The author also thanks the anonymous referees for useful suggestions and comments. This work was supported by the ERC grants QFTCMPS and SIQS, and by the cluster of excellence EXC201 Quantum Engineering and Space–Time Research. The publication of this article was funded by the Open Access Fund of the Leibniz Universität Hannover.

Theorem A1. Take a state $\rho (t)=U(t)\rho (0){U}^{\dagger }(t)$ of N non-interacting bosons or fermions and a measurement operator counting the number of particles in some modes $M={\sum }_{i}{b}_{i}^{\dagger }{b}_{i}$, where ${b}_{i}=a(| {\phi }_{i}\rangle )$. Then there exist orthonormal single-particle states $| {\psi }_{\alpha }(t)\rangle$, evolving via the corresponding single-particle Hamiltonian, such that

$$\begin{eqnarray}&&\mathrm{tr}[\rho (t)M]=\displaystyle \sum _{\alpha }{n}_{\alpha }\mathrm{tr}[{\psi }_{\alpha }(t)P],\end{eqnarray} \tag{ A1 }$$

where ${n}_{\alpha }$ are occupation numbers adding up to $N$, $P={\sum }_{i}| {\phi }_{i}\rangle \langle {\phi }_{i}|$ and ${\psi }_{\alpha }(t)=| {\psi }_{\alpha }(t)\rangle \langle {\psi }_{\alpha }(t)|$.

Proof. There exists a complete set of independent annihilation operators ${a}_{\alpha }=a(| {\psi }_{\alpha }\rangle )$, such that

$$\begin{eqnarray}&&\mathrm{tr}[\rho (0){a}_{\alpha }^{\dagger }{a}_{\beta }]=0\ \mathrm{if}\ \alpha \ne \beta .\end{eqnarray} \tag{ A2 }$$

To see this, take any complete set of independent annihilation operators ${d}_{\alpha }$. The matrix $\mathrm{tr}[\rho \;{d}_{\alpha }^{\dagger }{d}_{\beta }]={C}_{\alpha \beta }$ is Hermitian. So there exists a unitary ${U}_{\alpha \beta }$, such that in terms of

$$\begin{eqnarray}&&{a}_{\alpha }^{\dagger }=\displaystyle \sum _{\beta }{U}_{\alpha \beta }{d}_{\beta }^{\dagger }.\end{eqnarray} \tag{ A3 }$$

$\mathrm{tr}[\rho (0){a}_{\alpha }^{\dagger }{a}_{\beta }]$ is diagonal. So it makes sense to work with the ${a}_{\alpha }$, which determine the states $| {\psi }_{\alpha }(0)\rangle =| {\psi }_{\alpha }\rangle$ mentioned in the statement of the theorem via ${a}_{\alpha }=a(| {\psi }_{\alpha }\rangle )$. Next, expand

$$\begin{eqnarray}&&{a}^{\dagger }(| {\phi }_{i}\rangle )={b}_{i}^{\dagger }=\displaystyle \sum _{\alpha }{c}_{i,\alpha }{a}_{\alpha }^{\dagger },\end{eqnarray} \tag{ A4 }$$

where ${c}_{i,\alpha }$ are complex numbers. Then

$$\begin{eqnarray}\mathrm{tr}[\rho (0){b}_{i}^{\dagger }{b}_{i}] & = & \displaystyle \sum _{\alpha }| {c}_{i,\alpha }{| }^{2}\mathrm{tr}[\rho (0){a}_{\alpha }^{\dagger }{a}_{\alpha }]\\ & = & \displaystyle \sum _{\alpha }{n}_{\alpha }| {c}_{i,\alpha }{| }^{2},\end{eqnarray} \tag{ A5 }$$

where ${n}_{\alpha }$ is the number of particles in mode α. Next, we use ${c}_{i,\alpha }=\{{a}_{\alpha },{b}_{i}^{\dagger }\}=\langle {\psi }_{\alpha }| {\phi }_{i}\rangle$ for fermions, or ${c}_{i,\alpha }=[{a}_{\alpha },{b}_{i}^{\dagger }]=\langle {\psi }_{\alpha }| {\phi }_{i}\rangle$ for bosons, to get

$$\begin{eqnarray}&&\mathrm{tr}[\rho (0){b}_{i}^{\dagger }{b}_{i}]=\displaystyle \sum _{\alpha }{n}_{\alpha }\langle {\psi }_{\alpha }| {\phi }_{i}\rangle \langle {\phi }_{i}| {\psi }_{\alpha }\rangle .\end{eqnarray} \tag{ A6 }$$

Therefore

$$\begin{eqnarray}\mathrm{tr}[\rho (0)M] & = & \displaystyle \sum _{\alpha }{n}_{\alpha }\langle {\psi }_{\alpha }| P| {\psi }_{\alpha }\rangle \\ & = & \displaystyle \sum _{\alpha }{n}_{\alpha }\mathrm{tr}[{\psi }_{\alpha }P],\end{eqnarray} \tag{ A7 }$$

where $P={\sum }_{i}| {\phi }_{i}\rangle \langle {\phi }_{i}|$. For a state like ${a}_{1}^{\dagger }...{a}_{N}^{\dagger }| 0\rangle$, the first N occupation numbers ${n}_{\alpha }$ are one, so we get

$$\begin{eqnarray}&&\mathrm{tr}[\rho (0)M]=\displaystyle \sum _{\alpha =1}^{N}\mathrm{tr}[{\psi }_{\alpha }P].\end{eqnarray} \tag{ A8 }$$

To account for the dependence on time in the formula, we use $a(| {\psi }_{\alpha }(t)\rangle )={a}_{\alpha }(t)=U(t){a}_{\alpha }{U}^{\dagger }(t)$ and

$$\begin{eqnarray}&&\mathrm{tr}[\rho (0){a}_{\alpha }^{\dagger }{a}_{\beta }]=\mathrm{tr}[\rho (t){a}_{\alpha }^{\dagger }(t){a}_{\beta }(t)].\end{eqnarray} \tag{ A9 }$$

The end result is

$$\begin{eqnarray}&&\mathrm{tr}[\rho (t)M]=\displaystyle \sum _{\alpha }{n}_{\alpha }\mathrm{tr}[{\psi }_{\alpha }(t)P].\end{eqnarray} \tag{ A10 }$$

□

It will be useful to prove the following lemma.

Lemma B1. Let ${a}_{\alpha }^{\dagger }={a}^{\dagger }(| {\psi }_{\alpha }\rangle )$ be a complete set of creation operators with $\{{a}_{\alpha },{a}_{\beta }^{\dagger }\}={\delta }_{\alpha \beta }$. And suppose we have a state ${a}_{1}^{\dagger }...{a}_{N}^{\dagger }| 0\rangle$ with corresponding density operator $\rho$. And take a measurement operator counting the number of particles in some modes $M={\sum }_{i}{b}_{i}^{\dagger }{b}_{i}$, where ${b}_{i}=a(| {\phi }_{i}\rangle )$. Then, we have that

$$\begin{eqnarray}&&\mathrm{tr}[\rho {M}^{2}]\leqslant \mathrm{tr}{[\rho M]}^{2}+\mathrm{tr}[\rho M]\end{eqnarray} \tag{ B1 }$$

for fermions. And for bosons, we have

$$\begin{eqnarray}&&\mathrm{tr}[\rho {M}^{2}]\leqslant \mathrm{tr}{[\rho M]}^{2}+\mathrm{tr}[\rho M]+N.\end{eqnarray} \tag{ B2 }$$

Proof. First

$$\begin{eqnarray}&&\mathrm{tr}[\rho {M}^{2}]=\displaystyle \sum _{i,j}\mathrm{tr}[\rho ({b}_{i}^{\dagger }{b}_{i})({b}_{j}^{\dagger }{b}_{j})].\end{eqnarray} \tag{ B3 }$$

Each term can be written as

$$\begin{eqnarray}&&\mathrm{tr}[\rho ({b}_{i}^{\dagger }{b}_{i})({b}_{j}^{\dagger }{b}_{j})]=\mathrm{tr}[\rho \;{b}_{j}^{\dagger }{b}_{i}^{\dagger }{b}_{i}{b}_{j}]+{\delta }_{{ij}}\mathrm{tr}[\rho \;{b}_{i}^{\dagger }{b}_{i}],\end{eqnarray} \tag{ B4 }$$

which holds for bosons or fermions. To make sense of this, we write

$$\begin{eqnarray}&&{b}_{i}^{\dagger }=\displaystyle \sum _{\alpha }{c}_{i,\alpha }{a}_{\alpha }^{\dagger },\end{eqnarray} \tag{ B5 }$$

where ${c}_{i,\alpha }$ are complex numbers. It will turn out below that only terms with $\alpha \in \{1,\ldots ,N\}={ \mathcal V }$ contribute. Now we use the identity

$$\begin{eqnarray}&&\langle 0| {a}_{N}...{a}_{1}({a}_{\alpha }^{\dagger }{a}_{\beta }){a}_{1}^{\dagger }...{a}_{N}^{\dagger }| 0\rangle ={\delta }_{\alpha \beta }.\end{eqnarray} \tag{ B6 }$$

to get

$$\begin{eqnarray}&&\mathrm{tr}[\rho \;{b}_{i}^{\dagger }{b}_{i}]=\displaystyle \sum _{\alpha ,\beta \in { \mathcal V }}{c}_{i,\alpha }{c}_{i,\beta }^{* }{\delta }_{\alpha \beta }=\displaystyle \sum _{\alpha \in { \mathcal V }}| {c}_{i,\alpha }{| }^{2}.\end{eqnarray} \tag{ B7 }$$

For fermions, we use the identity

$$\begin{eqnarray}&&\langle 0| {a}_{N}...{a}_{1}({a}_{\alpha }^{\dagger }{a}_{\beta }^{\dagger }{a}_{\gamma }{a}_{\varepsilon }){a}_{1}^{\dagger }...{a}_{N}^{\dagger }| 0\rangle =[{\delta }_{\alpha \varepsilon }{\delta }_{\beta \gamma }-{\delta }_{\alpha \gamma }{\delta }_{\beta \varepsilon }]\end{eqnarray} \tag{ B8 }$$

to get

$$\begin{eqnarray}\mathrm{tr}[\rho {b}_{j}^{\dagger }{b}_{i}^{\dagger }{b}_{i}{b}_{j}] & = & \displaystyle \sum _{\alpha ,\beta ,\gamma ,\varepsilon \in { \mathcal V }}{c}_{j,\alpha }{c}_{i,\beta }{c}_{i,\gamma }^{* }{c}_{j,\varepsilon }^{* }[{\delta }_{\alpha \varepsilon }{\delta }_{\beta \gamma }-{\delta }_{\alpha \gamma }{\delta }_{\beta \varepsilon }]\\ & = & \displaystyle \sum _{\alpha ,\beta \in { \mathcal V }}[| {c}_{j,\alpha }{| }^{2}| {c}_{i,\beta }{| }^{2}-({c}_{j,\alpha }{c}_{i,\alpha }^{* })({c}_{i,\beta }{c}_{j,\beta }^{* })].\end{eqnarray} \tag{ B9 }$$

Now, using ${c}_{j,\alpha }=\{{b}_{j}^{\dagger },{a}_{\alpha }\}=\langle {\psi }_{\alpha }| {\phi }_{j}\rangle$, we get

$$\begin{eqnarray}\displaystyle \sum _{\alpha \in { \mathcal V }}{c}_{j,\alpha }{c}_{i,\alpha }^{* } & = & \langle {\phi }_{i}| Q| {\phi }_{j}\rangle ,\\ \mathrm{where}\ Q & = & \displaystyle \sum _{\alpha \in { \mathcal V }}| {\psi }_{\alpha }\rangle \langle {\psi }_{\alpha }| .\end{eqnarray} \tag{ B10 }$$

Then the second term in equation (B9) becomes

$$\begin{eqnarray}&&\displaystyle \sum _{\alpha ,\beta \in { \mathcal V }}({c}_{j,\alpha }{c}_{i,\alpha }^{* })({c}_{i,\beta }{c}_{j,\beta }^{* })=| \langle {\phi }_{i}| Q| {\phi }_{j}\rangle {| }^{2},\end{eqnarray} \tag{ B11 }$$

which is positive. Therefore

$$\begin{eqnarray}\mathrm{tr}[\rho \;{b}_{j}^{\dagger }{b}_{i}^{\dagger }{b}_{i}{b}_{j}] & \leqslant & \displaystyle \sum _{\alpha ,\beta \in { \mathcal V }}| {c}_{i,\alpha }{| }^{2}| {c}_{j,\beta }{| }^{2}\\ & = & \mathrm{tr}[\rho \;{b}_{i}^{\dagger }{b}_{i}]\mathrm{tr}[\rho \;{b}_{j}^{\dagger }{b}_{j}].\end{eqnarray} \tag{ B12 }$$

Putting everything together gives

$$\begin{eqnarray}&&\mathrm{tr}[\rho {M}^{2}]\leqslant \mathrm{tr}{[\rho M]}^{2}+\mathrm{tr}[\rho M].\end{eqnarray} \tag{ B13 }$$

For bosons, the result is similar, but the identity in equation (B8) is replaced by

$$\begin{eqnarray}&&\;\langle 0| {a}_{N}...{a}_{1}({a}_{\alpha }^{\dagger }{a}_{\beta }^{\dagger }{a}_{\gamma }{a}_{\varepsilon }){a}_{1}^{\dagger }...{a}_{N}^{\dagger }| 0\rangle =[{\delta }_{\alpha \varepsilon }{\delta }_{\beta \gamma }+{\delta }_{\alpha \gamma }{\delta }_{\beta \varepsilon }](1-{\delta }_{\alpha \beta }).\end{eqnarray} \tag{ B14 }$$

Because of this, following similar steps to those used to get equation (B9), we get

$$\begin{eqnarray}&&\mathrm{tr}[\rho \;{b}_{j}^{\dagger }{b}_{i}^{\dagger }{b}_{i}{b}_{j}]\leqslant \displaystyle \sum _{\alpha ,\beta \in { \mathcal V }}| {c}_{j,\alpha }{| }^{2}| {c}_{i,\beta }{| }^{2}+| \langle {\phi }_{i}| Q| {\phi }_{j}\rangle {| }^{2}.\end{eqnarray} \tag{ B15 }$$

The extra term arising in our expression for $\mathrm{tr}[\rho {M}^{2}]$ is

$$\begin{eqnarray}\displaystyle \sum _{{ij}}| \langle {\phi }_{i}| Q| {\phi }_{j}\rangle {| }^{2} & = & \displaystyle \sum _{{ij}}\langle {\phi }_{i}| Q| {\phi }_{j}\rangle \langle {\phi }_{j}| Q| {\phi }_{i}\rangle \\ & = & \mathrm{tr}[{PQPQ}]\leqslant N,\end{eqnarray} \tag{ B16 }$$

where $P={\sum }_{i}| {\phi }_{i}\rangle \langle {\phi }_{i}|$ and we used $\mathrm{rank}(Q)=N$.□

A corollary of this is that for pure states of bosons or fermions of the form ${a}_{1}^{\dagger }...{a}_{N}^{\dagger }| 0\rangle$, the fluctuations satisfy

$$\begin{eqnarray}&&{\sigma }_{M}^{2}=\mathrm{tr}[\rho {M}^{2}]-\mathrm{tr}{[\rho M]}^{2}=O(N),\end{eqnarray} \tag{ B17 }$$

where N is the number of fermions or bosons in the system. Furthermore, one can show that for N bosons in the same mode, one also gets ${\sigma }_{M}^{2}=O(N)$.

To say something more general about the fluctuations in fermion systems, we can also prove the following result.

Theorem B1. Given a non-interacting $N$ fermion system with corresponding single-particle Hamiltonian that has no degenerate energy levels and no degenerate energy gaps, then the time-average fluctuations satisfy

$$\begin{eqnarray}&&\langle {\sigma }_{M}(\rho (t))\rangle \leqslant \sqrt{N},\end{eqnarray} \tag{ B18 }$$

when the expectation value of $M$ on any infinite-time average state of $N$ particles is independent of the initial state. Examples where this is true include the gases discussed in the examples in the main text and systems where $M$ counts the number of particles in a spatial region, provided the Hamiltonian is such that $\langle M\rangle$, the time-average observable in the Heisenberg picture, is proportional to the total number operator.

Proof. Key to this result is the following inequality

$$\begin{eqnarray}\langle {\sigma }_{M}(\rho (t))\rangle & \leqslant & \sqrt{\langle \mathrm{tr}[\rho (t){M}^{2}]-\mathrm{tr}{[\rho (t)M]}^{2}\rangle }\\ & = & \sqrt{\mathrm{tr}[\langle \rho \rangle {M}^{2}]-\langle \mathrm{tr}{[\rho (t)M]}^{2}\rangle }\\ & \leqslant & \sqrt{\mathrm{tr}[\langle \rho \rangle {M}^{2}]-\mathrm{tr}{[\langle \rho \rangle M]}^{2}}\\ & = & {\sigma }_{M}(\langle \rho \rangle ),\end{eqnarray} \tag{ B19 }$$

where we used concavity of the square root in the first line and convexity of $f(x)={x}^{2}$ to get to the second last line.

So it remains to calculate

$$\begin{eqnarray}&&{\sigma }_{M}^{2}(\langle \rho \rangle )=\mathrm{tr}[\langle \rho \rangle {M}^{2}]-\mathrm{tr}{[\langle \rho \rangle M]}^{2}.\end{eqnarray} \tag{ B20 }$$

As $\langle \rho \rangle$ is the infinite-time average of $\rho (t)$, it follows that

$$\begin{eqnarray}&&\langle \rho \rangle =\displaystyle \sum _{E}{p}_{E}\;{\omega }_{E},\end{eqnarray} \tag{ B21 }$$

where p_E is a normalized probability distribution and ${\omega }_{E}$ is a state with support only on the energy eigenspace corresponding to energy E. Here E labels the energy of the N fermion system. So E is a sum of N single-particle energies. Let us write the creation operator that creates a fermion with energy E_i as ${a}_{i}^{\dagger }$. Now, the support of ${\omega }_{E}$ is spanned by the states

$$\begin{eqnarray}&&{a}_{{i}_{1}}^{\dagger }...{a}_{{i}_{N}}^{\dagger }| 0\rangle ,\end{eqnarray} \tag{ B22 }$$

with ${E}_{{i}_{1}}\;+\;\cdots \;+\;{E}_{{i}_{N}}=E$. Therefore, given two different configurations with energy E, $\{{i}_{1},\ldots ,{i}_{N}\}\ne \{{j}_{1},\ldots ,\;{j}_{N}\}$

$$\begin{eqnarray}&&\langle 0| {a}_{{i}_{N}}...{a}_{{i}_{1}}({a}_{\alpha }^{\dagger }{a}_{\beta }){a}_{{j}_{1}}^{\dagger }...{a}_{{j}_{N}}^{\dagger }| 0\rangle =0\end{eqnarray} \tag{ B23 }$$

because single-particle energy levels are non-degenerate. This implies

$$\begin{eqnarray}&&\mathrm{tr}[{\omega }_{E}\;{a}_{\alpha }^{\dagger }{a}_{\beta }]=\displaystyle \sum _{\vec{i}\in {C}_{E}}q(\vec{i})\mathrm{tr}[P(\vec{i}){a}_{\alpha }^{\dagger }{a}_{\beta }],\end{eqnarray} \tag{ B24 }$$

where $\vec{i}$ is short for $\{{i}_{1},\ldots ,{i}_{N}\}$, C_E denotes all $\{{i}_{1},\ldots ,{i}_{N}\}$ such that ${E}_{{i}_{1}}\;+\;\cdots \;+\;{E}_{{i}_{N}}=E$, $q(\vec{i})$ is a normalized probability distribution such that ${\sum }_{\vec{i}\in {C}_{E}}q(\vec{i})=1$ and $P(\vec{i})$ is the projector onto ${a}_{{i}_{1}}^{\dagger }...{a}_{{i}_{N}}^{\dagger }| 0\rangle$. Similarly, it is a consequence of (single-particle) non-degenerate energy gaps that

$$\begin{eqnarray}&&\mathrm{tr}[{\omega }_{E}\;{a}_{\alpha }^{\dagger }{a}_{\beta }^{\dagger }{a}_{\gamma }{a}_{\epsilon }]=\displaystyle \sum _{\vec{i}\in {C}_{E}}q(\vec{i})\mathrm{tr}[P(\vec{i}){a}_{\alpha }^{\dagger }{a}_{\beta }^{\dagger }{a}_{\gamma }{a}_{\epsilon }].\end{eqnarray} \tag{ B25 }$$

Now we can use lemma B1 to upper bound

$$\begin{eqnarray}\mathrm{tr}[\langle \rho \rangle {M}^{2}] & = & \displaystyle \sum _{E}\displaystyle \sum _{\vec{i}\in {C}_{E}}{p}_{E}\;q(\vec{i})\mathrm{tr}[P(\vec{i}){M}^{2}]\\ & \leqslant & \displaystyle \sum _{E}\displaystyle \sum _{\vec{i}\in {C}_{E}}{p}_{E}\;q(\vec{i})(\mathrm{tr}{[P(\vec{i})M]}^{2}+\mathrm{tr}[P(\vec{i})M]).\end{eqnarray} \tag{ B26 }$$

Next, we use $\mathrm{tr}[\omega M]=m$ independent of the state ω for fixed total particle number N, when ω is a time-average state, a special case of which is $\omega =P(\vec{i})$. So

$$\begin{eqnarray}&&\mathrm{tr}[\langle \rho \rangle {M}^{2}]-\mathrm{tr}{[\langle \rho \rangle M]}^{2}\leqslant {m}^{2}+m-{m}^{2}=m.\end{eqnarray} \tag{ B27 }$$

Finally, using the fact that the expectation value of M is bounded above by N on the N particle subspace leads to

$$\begin{eqnarray}&&{\sigma }_{M}{(\langle \rho \rangle )}^{2}\leqslant N\end{eqnarray} \tag{ B28 }$$

and therefore

$$\begin{eqnarray}&&\langle {\sigma }_{M}(\rho (t))\rangle \leqslant \sqrt{N}.\end{eqnarray} \tag{ B29 }$$

□

Equation (19) gives the distinguishability at time t

$$\begin{eqnarray}&&{D}_{P}(\sigma (t),\langle \sigma \rangle )=\left|\displaystyle \sum _{n\ne m}\mathrm{cos}[({n}^{2}-{m}^{2})\nu t]{\sigma }_{{nm}}{P}_{{mn}}\right|.\end{eqnarray} \tag{ C1 }$$

And from equation (18), we have

$$\begin{eqnarray}&&{\sigma }_{{nm}}{P}_{{mn}}=\displaystyle \frac{1}{N}\displaystyle \sum _{k=1}^{N}\langle n| {\psi }_{k}\rangle \langle {\psi }_{k}| m\rangle \langle m| P| n\rangle .\end{eqnarray} \tag{ C2 }$$

Evaluating the inner products, we get

$$\begin{eqnarray}\langle m| P| n\rangle & = & f(m,n)\\ \langle n| {\psi }_{k}\rangle & = & \sqrt{2}f(n,2k),\end{eqnarray} \tag{ C3 }$$

where

$$\begin{eqnarray}&&f(n,m)=\left\{\begin{array}{l}\tfrac{1}{\pi }\left[\tfrac{\mathrm{sin}\left[(n-m)\tfrac{\pi }{2}\right]}{n-m}-\tfrac{\mathrm{sin}\left[(n+m)\tfrac{\pi }{2}\right]}{n+m}\right]\;\mathrm{if}\ n\ne m,\\ \tfrac{1}{2}\ \mathrm{if}\ n=m.\end{array}\right.\end{eqnarray} \tag{ C4 }$$

Notice that, if both x and y are even or odd, unless x = y, then $f(x,y)=0$. The net result of this is

$$\begin{eqnarray}&&{\sigma }_{{nm}}{P}_{{mn}}=\displaystyle \frac{1}{N}\displaystyle \sum _{k=1}^{N}f{(n,2k)}^{2}\end{eqnarray} \tag{ C5 }$$

if $m=2k$ and n is odd and similarly if $n=2k$ and m is odd. All other terms are zero. Then subbing this into equation (C1), the distinguishability ${D}_{P}(\sigma (t),\langle \sigma \rangle )$ becomes

$$\begin{eqnarray}&&\displaystyle \frac{2}{N}\left|\displaystyle \sum _{n\ \mathrm{odd}}\displaystyle \sum _{k=1}^{N}\mathrm{cos}[({n}^{2}-4{k}^{2})\nu t]f{(n,2k)}^{2}\right|.\end{eqnarray} \tag{ C6 }$$

Furthermore, using $\mathrm{sin}(r\pi /2)={(-1)}^{(r-1)/2}$, which holds for odd r, we get

$$\begin{eqnarray}&&f{(n,2k)}^{2}=\displaystyle \frac{4}{{\pi }^{2}}\displaystyle \frac{4{k}^{2}}{{({n}^{2}-4{k}^{2})}^{2}},\end{eqnarray} \tag{ C7 }$$

for $n\ne 2k$.

Next, we can find a bound for the equilibration time. First, we make the substitution $n=2k+l$, noting that $l\gt -2k$ since $n\gt 0$, and l must be odd. It follows that $f{(2k+l,2k)}^{2}$ is peaked around small values of l, so we can focus on terms with $l\in { \mathcal S }$, where ${ \mathcal S }$ contains all odd integers from $-K$ to K. To quantify the resulting error, we use

$$\begin{eqnarray}&&f{(2k+l,2k)}^{2}=\displaystyle \frac{4}{{\pi }^{2}}\displaystyle \frac{4{k}^{2}}{{l}^{2}{(4k+l)}^{2}}\leqslant \displaystyle \frac{4}{{\pi }^{2}{l}^{2}}.\end{eqnarray} \tag{ C8 }$$

The sum of all terms with $l\notin { \mathcal S }$ can be bounded above by

$$\begin{eqnarray}\displaystyle \sum _{\mathrm{odd}\ l\notin { \mathcal S }}f{(2k+l,2k)}^{2} & \leqslant & \displaystyle \sum _{\mathrm{odd}\ l\notin { \mathcal S }}\displaystyle \frac{4}{{\pi }^{2}{l}^{2}}\\ & \leqslant & \displaystyle \sum _{l=-2k+1}^{-K-2}\displaystyle \frac{4}{{\pi }^{2}{l}^{2}}+\displaystyle \sum _{l=K+2}^{\infty }\displaystyle \frac{4}{{\pi }^{2}{l}^{2}}\\ & \leqslant & \displaystyle \frac{8}{{\pi }^{2}}\displaystyle \sum _{l=K+2}^{\infty }\displaystyle \frac{1}{{l}^{2}}\leqslant \displaystyle \frac{8}{{\pi }^{2}}\displaystyle \sum _{l=K+2}^{\infty }\left(\displaystyle \frac{1}{l-1}-\displaystyle \frac{1}{l}\right)\\ & = & \displaystyle \frac{8}{{\pi }^{2}}\displaystyle \frac{1}{K+1}.\end{eqnarray} \tag{ C9 }$$

So neglecting terms corresponding to $l\notin { \mathcal S }$ results in an upper bound for ${D}_{P}(\sigma (t),\langle \sigma \rangle )$ of

$$\begin{eqnarray}&&\displaystyle \frac{2}{N}\left|\displaystyle \sum _{l\in { \mathcal S }}\displaystyle \sum _{k=1}^{N}\mathrm{cos}[(4{kl}+{l}^{2})\nu t]\displaystyle \frac{4}{{\pi }^{2}}\displaystyle \frac{4{k}^{2}}{{l}^{2}{(4k+l)}^{2}}\right|+\displaystyle \frac{16}{{\pi }^{2}}\displaystyle \frac{1}{K+1},\end{eqnarray} \tag{ C10 }$$

which follows from the triangle inequality and $| \mathrm{cos}(x)| \leqslant 1$. Next, as $f(2k+l,2k)$ is awkward to work with, we use

$$\begin{eqnarray}\left|\displaystyle \frac{1}{{\pi }^{2}{l}^{2}}-f{(2k+l,2k)}^{2}\right| & = & \displaystyle \frac{1}{{\pi }^{2}| l| }\left|\displaystyle \frac{8k+l}{{(4k+l)}^{2}}\right|\\ & \leqslant & \displaystyle \frac{3}{2{\pi }^{2}}\displaystyle \frac{1}{k| l| },\end{eqnarray} \tag{ C11 }$$

where we used the triangle inequality and the fact that $1/(4k+l)\lt 1/(2k)$, since $l\gt -2k$.

This allows us to upper bound ${D}_{P}(\sigma (t),\langle \sigma \rangle )$ by

$$\begin{eqnarray}&&\displaystyle \frac{2}{N}\left|\displaystyle \sum _{l\in { \mathcal S }}\displaystyle \sum _{k=1}^{N}\mathrm{cos}[(4{kl}+{l}^{2})\nu t]\displaystyle \frac{1}{{\pi }^{2}{l}^{2}}\right|+\mu ,\end{eqnarray} \tag{ C12 }$$

where

$$\begin{eqnarray}&&\mu =\displaystyle \frac{16}{{\pi }^{2}}\displaystyle \frac{1}{K+1}+\displaystyle \frac{6}{{\pi }^{2}}\displaystyle \frac{[\mathrm{ln}(N)+1][\mathrm{ln}(K)+1]}{N}.\end{eqnarray} \tag{ C13 }$$

To get this, we used the triangle inequality and the inequality ${\sum }_{r=1}^{R}1/r\lt \mathrm{ln}(R)+1$.

Using the triangle inequality again, we get

$$\begin{eqnarray}{D}_{P}(\sigma (t),\langle \sigma \rangle ) & \leqslant & \displaystyle \frac{2}{N}\displaystyle \sum _{l\in { \mathcal S }}\displaystyle \frac{1}{{\pi }^{2}{l}^{2}}\left|\displaystyle \sum _{k=1}^{N}\mathrm{cos}[(4{kl}+{l}^{2})\nu t]\right|+\mu \\ & = & 2\displaystyle \sum _{l\in { \mathcal S }}\left|\displaystyle \frac{\mathrm{sin}[2{Nl}\nu t]\mathrm{cos}[(2(N+1)+l)l\nu t]}{N{\pi }^{2}{l}^{2}\mathrm{sin}[2l\nu t]}\right|+\mu \\ & \leqslant & 2\displaystyle \sum _{l\in { \mathcal S }}\left|\displaystyle \frac{\mathrm{sin}[2{Nl}\nu t]}{N{\pi }^{2}{l}^{2}\mathrm{sin}[2l\nu t]}\right|+\mu .\end{eqnarray} \tag{ C14 }$$

In the third line we used the identity [33]

$$\begin{eqnarray}&&\displaystyle \sum _{k=0}^{N-1}\mathrm{cos}(\alpha k+\phi )=\displaystyle \frac{\mathrm{sin}[N\alpha /2]\mathrm{cos}[(N-1)\alpha /2+\phi ]}{\mathrm{sin}(\alpha /2)}.\end{eqnarray} \tag{ C15 }$$

Let us look at each term in the sum in the last line of equation (C14) separately. They have period $\tfrac{\pi }{2l\nu }$, so we need only focus on this interval to find the time average of ${D}_{P}(\sigma (t),\langle \sigma \rangle )$.

When $2l\nu t$ is close to 0 or π, the $\mathrm{sin}[2l\nu t]$ in the denominator in the last line in equation (C14) is small. So for t such that $2l\nu t\in \left[0,\tfrac{1}{{Na}}\right)$ or $2l\nu t\in \left(\pi -\tfrac{1}{{Na}},\pi \right]$, where a is a small constant we will choose at the end, we bound

$$\begin{eqnarray}&&\left|\displaystyle \frac{\mathrm{sin}[2{Nl}\nu t]}{N{\pi }^{2}{l}^{2}\mathrm{sin}[2l\nu t]}\right|\leqslant \displaystyle \frac{1}{{\pi }^{2}{l}^{2}}.\end{eqnarray} \tag{ C16 }$$

When $2l\nu t\in \left[\tfrac{1}{{Na}},\tfrac{\pi }{2}\right]$, we can use the inequality $\mathrm{sin}(x)\geqslant 2x/\pi$ for all $x\in [0,\pi /2]$ to get

$$\begin{eqnarray}&&\left|\displaystyle \frac{\mathrm{sin}[2{Nl}\nu t]}{N{\pi }^{2}{l}^{2}\mathrm{sin}[2l\nu t]}\right|\leqslant \displaystyle \frac{1}{N\pi | l{| }^{3}4\nu t}.\end{eqnarray} \tag{ C17 }$$

To find the time average of ${D}_{P}(\sigma (t),\langle \sigma \rangle )$, we use the fact that the average over $\left[\tfrac{\pi }{4l\nu },\tfrac{\pi }{2l\nu }\right]$ is the same as that over $\left[0,\tfrac{\pi }{4l\nu }\right]$ by symmetry. So we need only average each term over $\left[0,\tfrac{\pi }{4l\nu }\right]$. The result is

$$\begin{eqnarray}\left\langle \left|\displaystyle \frac{\mathrm{sin}[2{Nl}\nu t]}{N{\pi }^{2}{l}^{2}\mathrm{sin}[2l\nu t]}\right|\right\rangle & \leqslant & \displaystyle \frac{1}{{\pi }^{2}{l}^{2}}\displaystyle \frac{2}{\pi {Na}}+\displaystyle \frac{4l\nu }{\pi }\left({\displaystyle \int }_{1/(2{Nal}\nu )}^{\pi /(4l\nu )}\displaystyle \frac{{\rm{d}}t}{N\pi | l{| }^{3}4\nu t}\right)\\ & = & \displaystyle \frac{2}{{\pi }^{3}{l}^{2}\;{Na}}+\displaystyle \frac{\mathrm{ln}(\pi {Na}/2)}{N{\pi }^{2}{l}^{2}}.\end{eqnarray} \tag{ C18 }$$

Plugging this into equation (C14), we get

$$\begin{eqnarray}&&\langle {D}_{P}(\sigma (t),\langle \sigma \rangle )\rangle \leqslant \displaystyle \frac{4}{3\pi {Na}}+\displaystyle \frac{2\mathrm{ln}(\pi {Na}/2)}{3N}+\mu ,\end{eqnarray} \tag{ C19 }$$

where we used ${\sum }_{l\in { \mathcal S }}1/{l}^{2}\leqslant 2{\sum }_{l=1}^{\infty }1/{l}^{2}={\pi }^{2}/3$. If we choose K = N, then we see that

$$\begin{eqnarray}&&\langle {D}_{P}(\sigma (t),\langle \sigma \rangle )\rangle =O\left(\displaystyle \frac{\mathrm{ln}{(N)}^{2}}{N}\right).\end{eqnarray} \tag{ C20 }$$

So for large N, this is extremely small and equilibration occurs. In fact, as figure 4 shows, the time-average distinguishability decays faster with N than the bound here.

Zoom In Zoom Out Reset image size

In figure 2 we saw that ${D}_{P}(\sigma (t),\langle \sigma \rangle )$ becomes small and then stays small for most times. In order to find the equilibration time, we can upper bound the time it takes for the distinguishability to become small. Plugging $t=\tfrac{1}{2{Na}\nu }$ into the bound in equation (C17), gives the bound

$$\begin{eqnarray}&&{D}_{P}(\sigma (t),\langle \sigma \rangle )\leqslant \displaystyle \sum _{l\in { \mathcal S }}\displaystyle \frac{a}{\pi | l{| }^{3}}+\mu \leqslant \displaystyle \frac{\pi a}{3}+\mu .\end{eqnarray} \tag{ C21 }$$

Here we choose a to be a small constant such that the distinguishability is small at $t=\tfrac{1}{2{Na}\nu }$. Then the equilibration time is bounded above by

$$\begin{eqnarray}&&{T}_{\mathrm{eq}}=\displaystyle \frac{1}{2{Na}\nu }.\end{eqnarray} \tag{ C22 }$$

C.1. Three-dimensions

We can extend this to a three-dimensional example in a way similar to the extension to three-dimensions for a particle in a box in [26]. Suppose the initial state of the N fermion system is

$$\begin{eqnarray}&&\displaystyle \prod _{\vec{j}\in J}{a}^{\dagger }(| {\psi }_{\vec{j}}\rangle )| 0\rangle ,\end{eqnarray} \tag{ C23 }$$

where J is the set of three-component vectors with components in $\{1,\ldots ,{J}_{\mathrm{max}}\}$, so we have $N={J}_{\mathrm{max}}^{3}$. And let $| {\psi }_{\vec{j}}\rangle$ be the energy eigenstate for a particle in the left half of a box labelled by $\vec{j}$ analogous to $| {\psi }_{k}\rangle$ in one-dimension in equation (17). Suppose the observable we are considering is that which counts the number of particles in the left half of the box.

After mapping to the single-particle picture, the distinguishability becomes

$$\begin{eqnarray}&&{D}_{P}(\sigma (t),\langle \sigma \rangle )={D}_{{P}_{x}}({\sigma }_{x}(t),\langle {\sigma }_{x}\rangle ),\end{eqnarray} \tag{ C24 }$$

where $P={P}_{x}\otimes {{\mathbb{1}}}_{y}\otimes {{\mathbb{1}}}_{z}$ is the projector onto the left half of the box. Here ${\sigma }_{x}(t)$ is the reduced state of the system on ${{ \mathcal H }}_{x}$, the Hilbert space corresponding to the x degrees of freedom. We have

$$\begin{eqnarray}&&{\sigma }_{x}(0)=\displaystyle \frac{1}{{J}_{\mathrm{max}}}\displaystyle \sum _{j=1}^{{J}_{\mathrm{max}}}| {\psi }_{j}\rangle \langle {\psi }_{j}| ,\end{eqnarray} \tag{ C25 }$$

where $| {\psi }_{j}\rangle$ is the jth energy eigenstate of a particle trapped in the left-hand side of a one-dimensional box. So the problem is now equivalent to the one-dimensional example. Therefore, the equilibration timescale is at most ${T}_{\mathrm{eq}}=O(1/{J}_{\mathrm{max}})=O(1/{N}^{1/3})$.

Here we will derive a useful formula that shows when equilibration occurs. The proof is very similar to one in [8], mainly differing by using a different weight for the time average.

Lemma D1. Take a finite dimensional system evolving via a time independent Hamiltonian in the state $\sigma (t)$. For any operator A and time $T\gt 0$

$$\begin{eqnarray}&&\displaystyle \frac{{\langle | \mathrm{tr}[\sigma (t)A]-\mathrm{tr}[\langle \sigma \rangle A]{| }^{2}\rangle }_{T}}{{\parallel A\parallel }^{2}}\leqslant \displaystyle \frac{{c}_{1}}{{d}_{\mathrm{eff}}}\underset{\beta }{\mathrm{max}}\displaystyle \sum _{\alpha }{{\rm{e}}}^{-{({G}_{\alpha }-{G}_{\beta })}^{2}\tfrac{{T}^{2}}{16}},\end{eqnarray} \tag{ D1 }$$

where ${c}_{1}=e\tfrac{\sqrt{\pi }}{2}$ and ${G}_{\alpha }={E}_{i}-{E}_{j}$ denote the non-zero energy gaps, so ${E}_{i}\ne {E}_{j}$. Also

$$\begin{eqnarray}&&\displaystyle \frac{1}{{d}_{\mathrm{eff}}}=\displaystyle \sum _{E}{(\mathrm{tr}[\sigma (0){P}_{E}])}^{2},\end{eqnarray} \tag{ D2 }$$

where ${P}_{E}$ is the projector onto the energy eigenspace corresponding to energy $E$.

Proof. First, we will take $\sigma (t)$ to be pure, extending the result to mixed states at the end. Because $\sigma (t)=| \psi (t)\rangle \langle \psi (t)|$ is pure, we can choose an eigenbasis of H, where $| \psi (t)\rangle$ only overlaps with a single eigenstate $| n\rangle$ for each energy level. This means that degenerate energy levels will not cause any problems. The state at time t is

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

where ${c}_{n}=\langle n| \psi (0)\rangle$. The time-average state is $\langle \sigma \rangle ={\sum }_{n}| {c}_{n}{| }^{2}| n\rangle \langle n|$, and the effective dimension is given by $1/{d}_{\mathrm{eff}}={\sum }_{n}| {c}_{n}{| }^{4}=\mathrm{tr}[{\langle \sigma \rangle }^{2}]$.

Using the notation ${A}_{{ij}}=\langle i| A| j\rangle$, we have

$$\begin{eqnarray}&&\langle | \mathrm{tr}[\sigma (t)A]-\mathrm{tr}[\langle \sigma \rangle A]{| }^{2}{\rangle }_{T}={\left\langle {\left|\displaystyle \sum _{i\ne j}({c}_{j}^{* }{A}_{{ji}}{c}_{i}){{\rm{e}}}^{-{\rm{i}}({E}_{i}-{E}_{j})t}\right|}^{2}\right\rangle }_{T},\end{eqnarray} \tag{ D4 }$$

$$\begin{eqnarray}&&\qquad \leqslant \displaystyle \sum _{\begin{array}{c}i\ne j\\ k\ne l\end{array}}\;({c}_{j}^{* }{A}_{{ji}}{c}_{i}){({c}_{l}^{* }{A}_{{lk}}{c}_{k})}^{* }{\langle {{\rm{e}}}^{{\rm{i}}[({E}_{k}-{E}_{l})-({E}_{i}-{E}_{j})]t}\rangle }_{T}.\end{eqnarray} \tag{ D5 }$$

To make our expressions more concise, we denote non-zero energy gaps by ${G}_{\beta }={E}_{i}-{E}_{j}$, with $\beta =(i,j)$, where $i\ne j$. We also define the vector

$$\begin{eqnarray}&&{v}_{\beta }={v}_{(i,j)}={c}_{j}^{* }{A}_{{ji}}{c}_{i}\end{eqnarray} \tag{ D6 }$$

and the Hermitian matrix

$$\begin{eqnarray}&&{M}_{\alpha \beta }={\langle {{\rm{e}}}^{{\rm{i}}({G}_{\alpha }-{G}_{\beta })t}\rangle }_{T}.\end{eqnarray} \tag{ D7 }$$

Equation (D4) becomes

$$\begin{eqnarray}{\langle | \mathrm{tr}[\sigma (t)A]-\mathrm{tr}[\langle \sigma \rangle A]{| }^{2}\rangle }_{T} & = & \displaystyle \sum _{\alpha ,\beta }{v}_{\alpha }^{* }{M}_{\alpha \beta }{v}_{\beta }\\ & \leqslant & \parallel M\parallel \displaystyle \sum _{\alpha }| {v}_{\alpha }{| }^{2}\\ & = & \parallel M\parallel \displaystyle \sum _{i\ne j}| {c}_{i}{| }^{2}| {c}_{j}{| }^{2}| {A}_{{ji}}{| }^{2}\\ & \leqslant & \parallel M\parallel \displaystyle \sum _{i,j}| {c}_{i}{| }^{2}| {c}_{j}{| }^{2}| {A}_{{ji}}{| }^{2}\\ & = & \parallel M\parallel \mathrm{tr}[A\langle \sigma \rangle {A}^{\dagger }\langle \sigma \rangle ].\end{eqnarray} \tag{ D8 }$$

As $\mathrm{tr}[{A}^{\dagger }B]$ defines an inner product for operators, we may use the Cauchy–Schwartz inequality. We can also use the inequality $\mathrm{tr}[{PQ}]\leqslant \parallel P\parallel \mathrm{tr}[Q]$, which holds for P and Q positive operators. Then we get

$$\begin{eqnarray}{\langle | \mathrm{tr}[\sigma (t)A]-\mathrm{tr}[\langle \sigma \rangle A]{| }^{2}\rangle }_{T} & \leqslant & \parallel M\parallel \sqrt{\mathrm{tr}[{A}^{\dagger }A\;{\langle \sigma \rangle }^{2}]\mathrm{tr}[{{AA}}^{\dagger }{\langle \sigma \rangle }^{2}]}\\ & \leqslant & \parallel M\parallel {\parallel A\parallel }^{2}\mathrm{tr}[{\langle \sigma \rangle }^{2}]\\ & = & \displaystyle \frac{\parallel M\parallel {\parallel A\parallel }^{2}}{{d}_{\mathrm{eff}}}.\end{eqnarray} \tag{ D9 }$$

Next, we can use the inequality for matrix norms [34]

$$\begin{eqnarray}\parallel M\parallel & \leqslant & \sqrt{{| | | M| | | }_{\infty }{| | | M| | | }_{1}}\\ & = & \underset{\beta }{\mathrm{max}}\displaystyle \sum _{\alpha }| {M}_{\alpha \beta }| ,\end{eqnarray} \tag{ D10 }$$

where ${| | | M| | | }_{1}$ and ${| | | M| | | }_{\infty }$ are the column and row matrix norms respectively. The second line holds because M is hermitian, implying $\parallel | {M}_{\infty }| \parallel =\parallel | {M}_{1}| \parallel$.

Our next task is to deal with ${M}_{\alpha \beta }={\langle {{\rm{e}}}^{{\rm{i}}({G}_{\alpha }-{G}_{\beta })t}\rangle }_{T}$. This can be simplified by replacing the original time average over the interval $[0,T]$ by a differently weighted average [8, 35]. This works because the quantity we are averaging (see equation (D4)) is positive, and because the new weight f(t) satisfies $f(t)\geqslant 1/T$ on the interval $[0,T]$.

We will choose the weight to be a Gaussian. Then for any positive g(t), we have

$$\begin{eqnarray}\displaystyle \frac{1}{T}{\displaystyle \int }_{0}^{T}{\rm{d}}t\;g(t) & \leqslant & \displaystyle \frac{e}{T}{\displaystyle \int }_{0}^{T}{\rm{d}}t\;g(t){{\rm{e}}}^{-4{(t-T/2)}^{2}/{T}^{2}}\\ & \leqslant & \displaystyle \frac{e}{T}{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}t\;g(t){{\rm{e}}}^{-4{(t-T/2)}^{2}/{T}^{2}}.\end{eqnarray} \tag{ D11 }$$

Here, $e=\mathrm{exp}(1)$.

Employing this weighted averaging, the matrix elements of M are

$$\begin{eqnarray}&&{M}_{\alpha \beta }=\displaystyle \frac{e}{T}{\int }_{-\infty }^{\infty }{\rm{d}}t\;{{\rm{e}}}^{-4{(t-T/2)}^{2}/{T}^{2}}\;{{\rm{e}}}^{{\rm{i}}({G}_{\alpha }-{G}_{\beta })t}.\end{eqnarray} \tag{ D12 }$$

So we get

$$\begin{eqnarray}&&| {M}_{\alpha \beta }| ={c}_{1}\;{{\rm{e}}}^{-{({G}_{\alpha }-{G}_{\beta })}^{2}{T}^{2}/16},\end{eqnarray} \tag{ D13 }$$

where ${c}_{1}=e\sqrt{\pi }/2$. From equation (D10), we get

$$\begin{eqnarray}&&\parallel M\parallel \leqslant {c}_{1}\underset{\beta }{\mathrm{max}}\displaystyle \sum _{\alpha }{{\rm{e}}}^{-{({G}_{\alpha }-{G}_{\beta })}^{2}{T}^{2}/16}.\end{eqnarray} \tag{ D14 }$$

Finally, we arrive at

$$\begin{eqnarray}&&\displaystyle \frac{{\langle | \mathrm{tr}[\sigma (t)A]-\mathrm{tr}[\langle \sigma \rangle A]{| }^{2}\rangle }_{T}}{{\parallel A\parallel }^{2}}\leqslant \displaystyle \frac{{c}_{1}}{{d}_{\mathrm{eff}}}\underset{\beta }{\mathrm{max}}\displaystyle \sum _{\alpha }{{\rm{e}}}^{-{({G}_{\alpha }-{G}_{\beta })}^{2}\tfrac{{T}^{2}}{16}}.\end{eqnarray} \tag{ D15 }$$

The final step of the proof is to extend the result to mixed states by doing a purification [23], as in [5]. Denote the system's Hilbert space by ${{ \mathcal H }}_{S}$. Then we can define a pure state $| \psi (0)\rangle$ on ${{ \mathcal H }}_{S}\otimes {{ \mathcal H }}_{A}$, with $\mathrm{dim}({{ \mathcal H }}_{S})=\mathrm{dim}({{ \mathcal H }}_{A})$, with the property that the reduced state on the first system is ${\mathrm{tr}}_{A}[| \psi (0)\rangle \langle \psi (0)| ]=\sigma (0)$. We recover the original evolution $\sigma (t)$ of the first system by evolving $| \psi (t)\rangle$ under the joint Hamiltonian $H\otimes {\mathbb{1}}$. Crucially, the expectation value of any operator A on the state $\sigma (t)$ is the same as the expectation value of $A\otimes {\mathbb{1}}$ on the purified state on the joint system. Also, $\parallel A\parallel =\parallel A\otimes {\mathbb{1}}\parallel$. Finally, the effective dimension of the purified system is equal to the effective dimension of the original system, which can be seen from from $\mathrm{tr}[{P}_{E}\sigma (0)]={\rm{tr}}[{P}_{E}\otimes {\mathbb{1}}| \psi (0)\rangle \langle \psi (0)| ]$.□

Our remaining task is to simplify things in terms of the density of states. In the sum over α, we can separate out the time dependent term, which has ${G}_{\alpha }\ne {G}_{\beta }$, and evaluate the sum by making the density of states approximation. We make the replacement ${\sum }_{E}\;=\;\int {\rm{d}}E\;n(E)$, where n(E) is the density of states, so that

$$\begin{eqnarray}\underset{\beta }{\mathrm{max}}\displaystyle \sum _{\overset{\alpha }{{G}_{\alpha }\ne {G}_{\beta }}}{{\rm{e}}}^{-{({G}_{\alpha }-{G}_{\beta })}^{2}\tfrac{{T}^{2}}{16}} & = & \underset{\beta }{\mathrm{max}}\displaystyle \sum _{\begin{array}{c}{E}^{\prime }\\ {E}^{\prime }\ne E\\ E-{E}^{\prime }\ne {G}_{\beta }\end{array}}\displaystyle \sum _{E}{{\rm{e}}}^{-{(E-{E}^{\prime }-{G}_{\beta })}^{2}\tfrac{{T}^{2}}{16}}\\ & = & \underset{\beta }{\mathrm{max}}{\int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}{\rm{d}}{E}^{\prime }\;n({E}^{\prime }){\int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}{\rm{d}}E\;n(E){{\rm{e}}}^{-{(E-{E}^{\prime }-{G}_{\beta })}^{2}\tfrac{{T}^{2}}{16}}\\ & \leqslant & {n}_{\mathrm{max}}\underset{\beta }{\mathrm{max}}{\int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}{\rm{d}}{E}^{\prime }\;n({E}^{\prime }){\int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}{\rm{d}}E\;{{\rm{e}}}^{-{(E-{E}^{\prime }-{G}_{\beta })}^{2}\tfrac{{T}^{2}}{16}}\\ & \leqslant & {n}_{\mathrm{max}}{\int }_{{E}_{\mathrm{min}}}^{{E}_{\mathrm{max}}}{\rm{d}}{E}^{\prime }\;n({E}^{\prime })\displaystyle \frac{4\sqrt{\pi }}{T}\\ & = & {n}_{\mathrm{max}}\displaystyle \frac{4\sqrt{\pi }d}{T},\end{eqnarray} \tag{ D16 }$$

where ${n}_{\mathrm{max}}={\mathrm{max}}_{E}\;n(E)$ and d is the dimension of the particle's state space. We also used $\int {\rm{d}}E\;n(E)=d$ to get the last line.

We define D_G to be the maximum number of gaps ${G}_{\alpha }$ of the same size. In other words

$$\begin{eqnarray}&&{D}_{G}=\underset{\beta }{\mathrm{max}}\displaystyle \sum _{\overset{\alpha }{{G}_{\alpha }={G}_{\beta }}}1.\end{eqnarray} \tag{ D17 }$$

So finally we get

$$\begin{eqnarray}&&\displaystyle \frac{{\langle | \mathrm{tr}[\sigma (t)A]-\mathrm{tr}[\langle \sigma \rangle A]{| }^{2}\rangle }_{T}}{{\parallel A\parallel }^{2}}\leqslant \displaystyle \frac{{c}_{1}}{{d}_{\mathrm{eff}}}\left[{D}_{G}+{n}_{\mathrm{max}}\displaystyle \frac{4\sqrt{\pi }d}{T}\right].\end{eqnarray} \tag{ D18 }$$

It is good to point out that the density of states approximation misses degenerate energy gaps because they are measure zero.

We want to see when equilibration occurs for free lattice models, and also estimate the timescale. A key step is to prove equilibration with respect to single-mode measurements. A consequence is corollary 2, which implies equilibration occurs for any measurement on K local modes, provided K is relatively small, and the initial state is Gaussian and commutes with the total number operator.

Consider the observable $M={a}^{\dagger }(| \phi \rangle )a(| \phi \rangle )$, which counts the number of particles in the state $| \phi \rangle$. By applying equation (6), we have

$$\begin{eqnarray}&&\mathrm{tr}[\rho (t)M]=\displaystyle \sum _{j=1}^{N}{n}_{j}\mathrm{tr}[{\psi }_{j}(t)| \phi \rangle \langle \phi | ].\end{eqnarray} \tag{ E1 }$$

The trick now is to switch to the Heisenberg picture. Then

$$\begin{eqnarray}&&\mathrm{tr}[\rho (t)M]=\displaystyle \sum _{j=1}^{N}{n}_{j}\mathrm{tr}[{\psi }_{j}| \phi (-t)\rangle \langle \phi (-t)| ].\end{eqnarray} \tag{ E2 }$$

Because of this, we can think of $| \phi (-t)\rangle \langle \phi (-t)| =\sigma (-t)$ as the state of a particle.

For any fermionic state, or a bosonic state with at most one boson in N orthogonal modes, we have

$$\begin{eqnarray}&&{\rm{\Pi }}=\displaystyle \sum _{j=1}^{N}{n}_{j}{\psi }_{j}\leqslant {\mathbb{1}}.\end{eqnarray} \tag{ E3 }$$

So we can think of this as a measurement (POVM) operator. For a system of bosons with more than one boson in each mode, the result can be extended simply by factoring out the maximum number of bosons in a mode, but we will not include this in the following formulas. Then we have

$$\begin{eqnarray}| \mathrm{tr}[\rho (t)M]-\mathrm{tr}[\langle \rho \rangle M]| & = & | \mathrm{tr}[{\rm{\Pi }}\sigma (-t)]-\mathrm{tr}[{\rm{\Pi }}\langle \sigma \rangle ]| \\ & = & {D}_{{\rm{\Pi }}}(\sigma (-t),\langle \sigma \rangle ).\end{eqnarray} \tag{ E4 }$$

Applying equation (38) in theorem 2, we get

$$\begin{eqnarray}&&{\langle {D}_{{\rm{\Pi }}}(\sigma (-t),\langle \sigma \rangle )\rangle }_{T}\leqslant \sqrt{\displaystyle \frac{{c}_{1}}{{d}_{\mathrm{eff}}}\left[{D}_{G}+{c}_{2}\displaystyle \frac{{n}_{\mathrm{max}}d}{T}\right]}.\end{eqnarray} \tag{ E5 }$$

The first task is to estimate $1/{d}_{\mathrm{eff}}$. Because the measurement operator $M={a}^{\dagger }(| \phi (0)\rangle )a(| \phi (0)\rangle )$ is local, the state $| \phi (0)\rangle \langle \phi (0)|$ is localized. But for free lattice models, the energy eigenstates are spread out over the whole lattice: they often have the form $| \phi (\vec{p})\rangle | \vec{p}\rangle$, where $| \vec{p}\rangle$ is the momentum state with momentum $\vec{p}$, and $| \phi (\vec{p})\rangle$ is the state of the extra degree of freedom (spin or particle type, for example). In that case, given an energy eigenstate $| E\rangle$ we get $| \langle E| \phi (0)\rangle {| }^{2}\leqslant l/V$, where the factor of l appears because $| \phi (0)\rangle$ is spread over at most l lattice sites. This implies

$$\begin{eqnarray}&&\displaystyle \frac{1}{{d}_{\mathrm{eff}}}=\displaystyle \sum _{E}{(\mathrm{tr}[\sigma (0){P}_{E}])}^{2}\leqslant \displaystyle \frac{{n}_{{\rm{d}}}^{2}{l}^{2}d}{{V}^{2}},\end{eqnarray} \tag{ E6 }$$

where d is the dimension of the Hilbert space and n_d is the maximum degeneracy of the energy levels. We also have that $d={sV}$, where s is the number of orthogonal states at each site. For a spin 1/2 particle, we would have s = 2. Plugging this into equation (E5), we get

$$\begin{eqnarray}&&{\langle {D}_{{\rm{\Pi }}}(\sigma (-t),\langle \sigma \rangle )\rangle }_{T}\leqslant l\sqrt{\displaystyle \frac{{c}_{3}}{d}\left[{D}_{G}+{c}_{2}\displaystyle \frac{{n}_{\mathrm{max}}d}{T}\right]},\end{eqnarray} \tag{ E7 }$$

where ${c}_{3}={c}_{1}{n}_{{\rm{d}}}^{2}{s}^{2}$.

The last task is to estimate ${n}_{\mathrm{max}}={\mathrm{max}}_{E}\;n(E)$.

E.1. Density of states for lattice models

For lattice models, estimating $\;{\mathrm{max}}_{E}n(E)$ causes problems because n(E) often diverges. Fortunately, we can truncate to a slightly smaller energy range, such that $\;{\mathrm{max}}_{E}n(E)$ is bounded, and the error caused by the truncation is small.

Let us take a nearest-neighbour hopping model on a line as an example. The corresponding single-particle Hamiltonian is

$$\begin{eqnarray}&&H=\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{L}(| i\rangle \langle i+1| +| i+1\rangle \langle i| ),\end{eqnarray} \tag{ E8 }$$

where we are assuming translational invariance, so the site at $L+1$ is identified with site 1. Switching to momentum space diagonalizes this, and we get the dispersion relation $E(p)=\mathrm{cos}(p)$, where $p=2\pi k/L$, with $k\in \{1,\ldots ,L\}$. (In making the density of states approximation, we assume that L is large but finite.) The density of states is

$$\begin{eqnarray}n(E) & = & \displaystyle \frac{L}{\pi }\displaystyle \frac{{\rm{d}}p}{{\rm{d}}E}=\displaystyle \frac{L}{\pi }\displaystyle \frac{1}{| \mathrm{sin}(\mathrm{arccos}(E))| }\\ & = & \displaystyle \frac{L}{\pi }\displaystyle \frac{1}{\sqrt{1-{E}^{2}}}.\end{eqnarray} \tag{ E9 }$$

At $E=\pm 1$ this is infinite. However, we can truncate the state, neglecting all terms with $E\in [-1,-\mathrm{cos}({p}_{0}))\cup (\mathrm{cos}({p}_{0}),1]$, where we choose a fixed ${p}_{0}=2\pi {k}_{0}/L\gt 0$ to be small. This leads to a constant error in approximation of the state as there is a constant fraction of energy eigenstates with energy in this range. Defining P₀ to be the projector onto the subspace corresponding to the omitted energy range, we get

$$\begin{eqnarray}{(\parallel | \phi \rangle -({\mathbb{1}}-{P}_{0})| \phi \rangle {\parallel }_{2})}^{2} & = & 1-\langle \phi | ({\mathbb{1}}-{P}_{0})| \phi \rangle \\ & = & \langle \phi | {P}_{0}| \phi \rangle \\ & \leqslant & {N}_{0}\displaystyle \frac{l}{L},\end{eqnarray} \tag{ E10 }$$

where we used $| \langle E| \phi (0)\rangle {| }^{2}\leqslant l/V$ from the previous section to get the last line. And N₀ is the number of states with energy in the excluded set. Next, we can use that ${N}_{0}=4{k}_{0}=2{p}_{0}L/\pi$ to get

$$\begin{eqnarray}&&{(\parallel | \phi \rangle -({\mathbb{1}}-{P}_{0})| \phi \rangle {\parallel }_{2})}^{2}\leqslant \displaystyle \frac{2{p}_{0}l}{\pi }.\end{eqnarray} \tag{ E11 }$$

So the error in approximating $| \phi \rangle$ by a state restricted to the smaller energy range can be made small by choosing p₀ to be small. Furthermore, we have that

$$\begin{eqnarray}&&\;\underset{E}{\mathrm{max}}\;n(E)=\displaystyle \frac{L}{\pi }\displaystyle \frac{1}{| \mathrm{sin}({p}_{0})| }.\end{eqnarray} \tag{ E12 }$$

Crucially, p₀ is fixed and does not depend on the number of sites.

The same ideas apply to other dispersion relations (for example, those arising from translationally invariant Hamiltonians, possibly in higher spatial dimensions). See figure 5. The basic idea is to exclude regions where the density of states diverges, which corresponds to points where the dispersion relation is flat. The key point is that the density of states is a fixed function, and the fraction of energy eigenstates corresponding to regions where it is large remains constant. (Notice that the trivial Hamiltonian $H=\mathrm{constant}$ has a completely flat dispersion relation, so that this trick will not work in that case.) But generally for free lattice models we expect $\;{\mathrm{max}}_{E}\;n(E)\propto V$, where V is the number of lattice sites.

Zoom In Zoom Out Reset image size

E.2. From single-mode to multi-mode measurements

It is possible to relate K mode measurements to single-mode measurements if the state is Gaussian. We will only prove this here for fermionic Gaussian states, as the bosonic analogue is similar.

Corollary E1. Suppose $\rho =\rho (0)$ is Gaussian and satisfies $[\rho ,N]=0$, where $N$ is the total number operator. (This is still quite general, but it rules out BCS states.) Let $M$ be a measurement operator acting on $K$ modes on $l$ sites, where the modes are local, but not necessarily near each other. Then we have

$$\begin{eqnarray}&&{\langle | \mathrm{tr}[\rho M(t)]-\mathrm{tr}[\rho \langle M\rangle ]| \rangle }_{T}\leqslant {2}^{K+2}{mlK}\sqrt{{c}_{3}\left[\displaystyle \frac{{D}_{G}}{d}+{c}_{2}\displaystyle \frac{{n}_{\mathrm{max}}}{T}\right]},\end{eqnarray} \tag{ E13 }$$

with $m=\mathrm{max}| {m}_{{r}_{1},\ldots ,{r}_{2K}}|$, where ${m}_{{r}_{1},\ldots ,{r}_{2K}}$ are the coefficients of M when expanded in a fermionic operator basis on the $K$ modes. See equation (E14). And ${c}_{2}=4\sqrt{\pi }$ and ${c}_{3}={c}_{1}{n}_{{\rm{d}}}^{2}{s}^{2}$, where ${c}_{1}=e\sqrt{\pi }/2$.

Proof. Define ${c}_{\alpha }$, with $\alpha \in \{1,\ldots ,2K\}$, to be $2K$ Majorana operators generating the algebra for the K modes that M acts on. We can choose these ${c}_{\alpha }$ to make the covariance matrix simple, as we will see below. We will work in the Heisenberg picture here, but suppress the time dependence and just write ${c}_{\alpha }$ instead of ${c}_{\alpha }(t)$. We can expand M in terms of these Majorana operators:

$$\begin{eqnarray}&&M=\displaystyle \sum _{{r}_{1},\ldots ,{r}_{2K}=0,1}{m}_{{r}_{1},\ldots ,{r}_{2K}}{c}_{1}^{{r}_{1}}...{c}_{2K}^{{r}_{2K}}.\end{eqnarray} \tag{ E14 }$$

Using the triangle inequality, we get

$$\begin{eqnarray}&&| \mathrm{tr}[\rho M(t)]-\mathrm{tr}[\rho \langle M\rangle ]| \leqslant \displaystyle \sum _{{r}_{1},\ldots ,{r}_{2K}=0,1}| {m}_{{r}_{1},\ldots ,{r}_{2K}}| | \mathrm{tr}[\rho {c}_{1}^{{r}_{1}}...{c}_{2K}^{{r}_{2K}}]-\mathrm{tr}[\rho \langle {c}_{1}^{{r}_{1}}...{c}_{2K}^{{r}_{2K}}\rangle ]| .\end{eqnarray} \tag{ E15 }$$

We choose the ${c}_{\alpha }$ such that the covariance matrix, defined by

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{ij}}=\displaystyle \frac{{\rm{i}}}{2}\mathrm{tr}[\rho [{c}_{i},{c}_{j}]],\end{eqnarray} \tag{ E16 }$$

is in block diagonal form [36]

$$\begin{eqnarray}{\rm{\Gamma }}=\underset{n=1}{\overset{K}{\displaystyle \bigoplus }}\left(\begin{array}{ll}0 & {\lambda }_{n}\\ -{\lambda }_{n} & 0\end{array}\right),\end{eqnarray} \tag{ E17 }$$

with ${\lambda }_{n}=\tfrac{i}{2}\mathrm{tr}[\rho [{c}_{2n-1},{c}_{2n}]]$. Let R be the bit string $({r}_{1},\ldots ,{r}_{2K})$, and define $2w={\sum }_{i}{r}_{i}$. Because the state is Gaussian, we have [36]

$$\begin{eqnarray}&&\mathrm{tr}[\rho {c}_{1}^{{r}_{1}}...{c}_{2K}^{{r}_{2K}}]={{\rm{i}}}^{w}\mathrm{pf}({{\rm{\Gamma }}}^{R}),\end{eqnarray} \tag{ E18 }$$

where $\mathrm{pf}$ is the Pfaffian, and ${{\rm{\Gamma }}}^{R}$ is a submatrix of Γ restricted to the rows and columns labeled by i corresponding to r_i = 1. The Pfaffian is zero if any row or column is zero, so $\mathrm{pf}({{\rm{\Gamma }}}^{R})=0$, unless the string R only contains consecutive pairs of ones and zeros, such as $(1,1,0,0,1,1,\ldots )$.

The Pfaffian has two useful properties. The first is that $\mathrm{pf}\;\left(\begin{array}{ll}0 & \lambda \\ -\lambda & 0\end{array}\right)=\lambda$, and the second is that $\mathrm{pf}({A}_{1}\oplus \cdots \;\oplus \;{A}_{n})=\mathrm{pf}({A}_{1})\;...\;\mathrm{pf}({A}_{n})$. For any string R giving a non zero Pfaffian, we get

$$\begin{eqnarray}&&\mathrm{tr}[\rho {c}_{1}^{{r}_{1}}...{c}_{2K}^{{r}_{2K}}]={{\rm{i}}}^{w}\displaystyle \prod _{n=1}^{K}{\lambda }_{n}^{{r}_{2n}}.\end{eqnarray} \tag{ E19 }$$

And so

$$\begin{eqnarray}| \mathrm{tr}[\rho {c}_{1}^{{r}_{1}}...{c}_{2K}^{{r}_{2K}}]-\mathrm{tr}[\rho \langle {c}_{1}^{{r}_{1}}...{c}_{2K}^{{r}_{2K}}\rangle ]| & \leqslant & \left|\displaystyle \prod _{n=1}^{K}{\lambda }_{n}^{{r}_{n}}-\langle \displaystyle \prod _{n=1}^{K}{\lambda }_{n}^{{r}_{n}}\rangle \right|\\ & \leqslant & \left|\displaystyle \prod _{n=1}^{K}{\lambda }_{n}^{{r}_{n}}-\displaystyle \prod _{n=1}^{K}\langle {\lambda }_{n}^{{r}_{n}}\rangle \right|+\left\langle \left|\displaystyle \prod _{n=1}^{K}{\lambda }_{n}^{{r}_{n}}-\displaystyle \prod _{n=1}^{K}\langle {\lambda }_{n}^{{r}_{n}}\rangle \right|\right\rangle \\ & \leqslant & \displaystyle \sum _{n=1}^{K}| {\lambda }_{n}^{{r}_{n}}-\langle {\lambda }_{n}^{{r}_{n}}\rangle | +\displaystyle \sum _{n=1}^{K}\langle | {\lambda }_{n}^{{r}_{n}}-\langle {\lambda }_{n}^{{r}_{n}}\rangle | \rangle ,\end{eqnarray} \tag{ E20 }$$

where the third line follows from the triangle inequality and $| \langle f(t)\rangle | \leqslant \langle | f(t)| \rangle$. The last line follows from the triangle inequality, and repeated applications of $| {xy}-\langle x\rangle \langle y\rangle | \leqslant | x-\langle x\rangle | +| y-\langle y\rangle |$, which uses $| x|$, $| y| \leqslant 1$. So we need to focus on

$$\begin{eqnarray}&&| {\lambda }_{n}-\langle {\lambda }_{n}\rangle | =\left|\displaystyle \frac{{\rm{i}}}{2}\mathrm{tr}[\rho [{c}_{{\rm{i}}},{c}_{j}]]-\displaystyle \frac{{\rm{i}}}{2}\mathrm{tr}[\rho \langle [{c}_{i},{c}_{j}]\rangle ]\right|,\end{eqnarray} \tag{ E21 }$$

where $i=2n-1$ and $j=2n$.

For each pair c_i and c_j, we can always define creation operators ${d}_{i}^{\dagger }$ and ${d}_{j}^{\dagger }$, with

$$\begin{eqnarray}{c}_{i} & = & {\rm{i}}({d}_{i}-{d}_{i}^{\dagger })\\ {c}_{j} & = & {d}_{j}+{d}_{j}^{\dagger }.\end{eqnarray} \tag{ E22 }$$

We are not assuming that these creation operators correspond to orthogonal modes. Using $\mathrm{tr}[\rho {d}_{i}{d}_{j}]=0$, which follows because $[\rho ,N]=0$, we have

$$\begin{eqnarray}\displaystyle \frac{{\rm{i}}}{2}\mathrm{tr}[\rho [{c}_{i},{c}_{j}]] & = & \mathrm{tr}[\rho ({d}_{i}^{\dagger }{d}_{j}-{d}_{i}{d}_{j}^{\dagger })]\\ & = & \mathrm{tr}[\rho ({d}_{i}^{\dagger }{d}_{j}+{d}_{j}^{\dagger }{d}_{i})]-\kappa ,\end{eqnarray} \tag{ E23 }$$

where $\kappa =\{{d}_{i},{d}_{j}^{\dagger }\}$ is a constant. We can rewrite this as

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{i}}}{2}\mathrm{tr}[\rho [{c}_{i},{c}_{j}]]=\mathrm{tr}[\rho {b}_{i}^{\dagger }{b}_{i}]-\mathrm{tr}[\rho {b}_{j}^{\dagger }{b}_{j}]-\kappa ,\end{eqnarray} \tag{ E24 }$$

where ${b}_{i}=\tfrac{1}{\sqrt{2}}({d}_{i}+{d}_{j})$ and ${b}_{j}=\tfrac{1}{\sqrt{2}}({d}_{i}-{d}_{j})$. This allows us to apply corollary 1 to get

$$\begin{eqnarray}&&{\left\langle \left|\displaystyle \frac{{\rm{i}}}{2}\mathrm{tr}[\rho [{c}_{i},{c}_{j}]]-\displaystyle \frac{{\rm{i}}}{2}\mathrm{tr}[\rho \langle [{c}_{i},{c}_{j}]\rangle ]\right|\right\rangle }_{T}\leqslant 2l\sqrt{{c}_{3}\left[\displaystyle \frac{{D}_{G}}{d}+{c}_{2}\displaystyle \frac{{n}_{\mathrm{max}}}{T}\right]},\end{eqnarray} \tag{ E25 }$$

where l is the number of sites the measurement acts on. Putting this all together, we get

$$\begin{eqnarray}&&| \mathrm{tr}[\rho M(t)]-\mathrm{tr}[\rho \langle M\rangle ]| \leqslant 4{m}^{\prime }{lK}\sqrt{{c}_{3}\left[\displaystyle \frac{{D}_{G}}{d}+{c}_{2}\displaystyle \frac{{n}_{\mathrm{max}}}{T}\right]},\end{eqnarray} \tag{ E26 }$$

where

$$\begin{eqnarray}&&{m}^{\prime }=\displaystyle \sum _{{r}_{2},{r}_{4},\ldots ,{r}_{2K}=0,1}| {m}_{{r}_{2},{r}_{2},{r}_{4},{r}_{4},\ldots ,{r}_{2K}}| ,\end{eqnarray} \tag{ E27 }$$

where the sum only counts terms labelled by consecutive pairs of ones and zeros. Then, for simplicity, we can use ${m}^{\prime }\leqslant m{2}^{K}$, where $m=\mathrm{max}| {m}_{{r}_{1},\ldots ,{r}_{2K}}|$, since there are ${2}^{K}$ terms that contribute. It may also simplify things to use $K\leqslant {ls}$. This gives

$$\begin{eqnarray}&&| \mathrm{tr}[\rho M(t)]-\mathrm{tr}[\rho \langle M\rangle ]| \leqslant {2}^{{ls}+2}{{msl}}^{2}\sqrt{{c}_{3}\left[\displaystyle \frac{{D}_{G}}{d}+{c}_{2}\displaystyle \frac{{n}_{\mathrm{max}}}{T}\right]}.\end{eqnarray} \tag{ E28 }$$

□