Equilibration of Quantum Gases

Finding equilibration times is a major unsolved problem in physics with few analytical results. Here we look at equilibration times for quantum gases of bosons and fermions in the regime of negligibly weak interactions, a setting which not only includes paradigmatic systems such as gases confined to boxes, but also Luttinger liquids and the free superfluid Hubbard model. To do this, we focus on two classes of measurements: (i) coarse-grained observables, such as the number of particles in a region of space, and (ii) few-mode measurements, such as phase correlators and correlation functions. We show that, in this setting, equilibration occurs quite generally despite the fact that the particles are not interacting. Furthermore, for coarse-grained measurements the timescale is generally at most polynomial in the number of particles N, which is much faster than previous general upper bounds, which were exponential in N. For local measurements on lattice systems, the timescale is typically linear in the number of lattice sites. In fact, for one dimensional lattices, the scaling is generally linear in the length of the lattice, which is optimal. Additionally, we look at a few specific examples, one of which consists of N fermions initially confined on one side of a partition in a box. The partition is removed and the fermions equilibrate extremely quickly in time O(1/N).


I. INTRODUCTION
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.
FIG. 1: Microscopically, a cup of tea is never in equilibrium: the molecules are constantly moving around, but we cannot measure this. What we do measure is the temperature, according to which a hot cup of tea can reach equilibrium (room temperature). This highlights an important point about equilibration, which is that it only occurs when we account for physical restrictions on what we can actually measure. * Electronic address: farreltc@tcd.ie In [3][4][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][7][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][14][15][16][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 ul-tracold 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.

II. EQUILIBRATION
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 23 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] D M (ρ, σ) = 1 2 max where M denotes the set of measurements we can do, and {M i } denotes a POVM measurement, with the positive operators M i satisfying i M i = 1 1. POVM (Positive operator valued measure) 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 ρ(t) by ρ . If D M (ρ(t), ρ ) is small most of the time, then for all practical purposes ρ(t) is indistinguishable from its time average ρ 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 where M 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 M . If equilibration is to occur, we require that most of the time ∆ M (t) is smaller than some , with chosen so that M 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 tr[ρ(t)M ] with extremely high probability. Therefore, equilibration truly occurs.

III. GASES OF BOSONS AND FERMIONS
The key step in getting estimates of the equilibration time for N particle systems is equation (6) below, which will allow us to equate ∆ M (t) to the distinguishability for a single particle.
First, it will be useful to introduce some notation. Let H be a single-particle Hilbert space, and let |i denote an orthonormal basis. Then we can define creation operators a † i , acting on a fermionic Hilbert space, that create fermions corresponding to these states. Equivalently, we may say a † i creates a fermion in mode i. The fermionic Hilbert space is spanned by states with varying numbers of creation operators acting on |0 , the empty state. To avoid confusion, any state vectors written as kets are in the single-particle Hilbert space H, with the exception of |0 , which represents the empty state in a fermionic (or bosonic) system.
The creation operator that creates a particle corresponding to the single-particle state |ψ = i c i |i is a † (|ψ ) = i c i a † i . Suppose we have a single-particle Hamiltonian with discrete spectrum, where E labels the energies. There is a corresponding fermionic Hamiltonian, given by For any single-particle state |ψ , we also have 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 ρ(t) = U (t)ρ(0)U † (t) of N non-interacting bosons or fermions and a measurement operator counting the number of particles in some orthog- Then, there exist orthonormal single-particle states |ψ α (t) , evolving via the corresponding single-particle Hamiltonian, such that where n α are occupation numbers adding up to N , P = i |φ i φ i | and ψ α (t) = |ψ α (t) ψ α (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 †n1 1 ...a †n k n k |0 , where a † α = a † (|ψ α ), and |ψ α are some orthonormal singleparticle states. In this case, the states |ψ α mentioned in the theorem are already given.

Proof. Expand
where c i,α are complex numbers. Then where n α is the number of particles in mode α. Next, we use Therefore, where P = i |φ i φ i |. To incorporate the dependence on time, we use a α (t) = U (t)a α U † (t) and The end result is Notice that linearity of the time average, together with equation (6) implies A. Coarse-grained measurements We can apply this to ∆ M (t), noting that for the applications we are interested in M = 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 N , where is some very small constant. For equilibration to occur, we need ∆ M (t) to be small compared to most of the time. We get where σ(t) = 1 N j |ψ j (t) ψ 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 1 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 (tr 2 ) 1/2 , are bounded above by √ N , which is proved in appendix B. In fact, a large class of fermion systems have time-averaged fluctuations bounded above by √ N , as seen in appendix B. For large numbers of particles, (comparable to 10 23 , for example) √ N is small compared to our experimental precision N , and the fluctuations are not practically observable. Even for dilute gases with O(10 4 ) particles, √ N ∼ 100, so the fluctuations are of the order of 1% of the total particle number, which is still quite small.

B. Few-mode measurements
We are not just restricted to coarse-grained measurements. We can also discuss measurements involving a few modes. These could be single-site densities or correlation functions in the setting of lattice models. Or they could be phase correlators tr[ρ(t)a † i a j ], which are typically inferred from time-of-flight measurements [25].
We will return to this in section VI B, where we will see that for a large class of lattice systems any measurement on a small number of modes (small compared to the lattice size) will equilibrate. And the timescale will be relatively fast.

IV. EXAMPLE I: PARTICLES IN A BOX.
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. Above is a plot of the resulting single-particle distinguishability for the one dimensional example with N = 10 or 100 fermions and any number of bosons. Time is measured in units of the recurrence time, though for the initial state here there is another recurrence at half the recurrence time.
In general, for fermions the equilibration time is O(1/N ). For bosons, the system does not equilibrate, as can be seen from the figure. These plots were generated using equation (20).
Using equation (14), we can replace this N particle problem by a single-particle one, so which is plotted in figure 2. Here σ(t) is a state of a free particle in a box, σ is its time average and P is the projector onto the left hand side of the box.

A. 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 σ(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 where n > 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 where again k > 0 is an integer. The initial state of the single-particle system is with matrix elements σ nm = n|σ(0)|m . Similarly, the matrix elements of the projector onto the left half of the box are P mn = m|P |n . 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 rec = 2π ν since the energy levels are E n = νn 2 , where n is an integer greater than zero, and ν = π 2 2mL 2 . This is because all phases of density matrix elements in the energy basis e −i(En−Em)t are 1 at t = 2π ν . As in [26], this means that we need only study the system over the interval [0, T rec ]. In fact, with the particular initial state below, a recurrence actually occurs at T rec /2.
Evaluating the distinguishability at time t, we get where we used the fact that σ 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 for n = 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 D P (σ(t), σ ) = O(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 Here a is a small constant that we choose such that the distinguishability at t = T eq is small: which is also derived in appendix C. Interestingly, the timescale decreases with increasing particle number.

B. Bosons
The situation for N bosons is simpler. As they are initially at zero temperature, all N bosons are in the ground state. The corresponding initial single-particle state σ(0) is just the lowest energy state for a particle trapped on the left of the partition. This does not depend on N . By looking at the plot of the distinguishability in figure 2, it is clear that this system does not equilibrate because the distinguishability is large for most times.
So the behaviour of N bosons is very different from the fermion case. This is because of the exclusion principle: in the fermion case, the fermions had to occupy different energy levels and so the corresponding single-particle state was spread out over many energy levels. This is not the case for bosons at zero temperature.

V. EXAMPLE II: BOSONS AFTER A QUENCH
For our second example, suppose we have N bosons at zero temperature in a one dimensional harmonic trap with frequency ω 0 . We will consider what happens after two different quenches.  (14), we see that where ψ(t) = |ψ(t) ψ(t)| is a pure state of a singleparticle 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 infinitetime average of D P (ψ(t), ψ ) is numerically shown to scale like where it was assumed that the width of the initial wavefunction is small compared to the length of the box, meaning σ = 1/ √ 2mω 0 L. So for sufficiently narrow potentials (or sufficiently large boxes), equilibration occurs. Furthermore, the timescale for equilibration is shown to be [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.

B. 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 behaviour was observed [28]. Such behaviour 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 ω 0 and post-quench frequency ω was close to one: ω 0 /ω 1.3. Here we see equilibration when ω 0 /ω 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 where the P is the projector onto the region [−l, l]. The distinguishability is So we need only see if tr[P ψ(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] K(x, y, t) = which leads to the expression As ψ(x) is a Gaussian wavefunction, the y 1 and y 2 integrals are straightforward, leading to where with γ = ω 0 /ω. Next we use the approximation for the error function [30] erf where the maximum error for any x is around 0.00012, and b 0.147. The result is that 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; γ = ω 0 /ω, which is the ratio of trap strengths before and after the quench; and mω 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ω 0 ). As we can see from figure 3, as γ becomes bigger and bigger, tr[P ψ(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, tr[P ψ(t)] tends to zero. This holds for all times, except when ωt = nπ, with n ∈ Z.
In [28], oscillatory behaviour was seen at γ = 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 tr[P ψ(t)] to reach p 1. Using equation (34) and log(1 − p 2 ) −p 2 , we get  [28]. We see oscillatory behaviour for γ approximately between 4 and 10. We see that, as γ becomes larger, tr[P ψ(t)] is small most of the time, and the system equilibrates. These plots were generated using equation (34).
Since p is small, this requires α(t)l 2 to be small. Using the earlier choice l 2 = 4/(mω 0 ), we get For large γ, this is satisfied at where we assumed that t was small compared to 1/ω, and used sin(x) x for small x.

VI. EQUILIBRATION IN GENERAL
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][7][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: 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 A , and let σ(t) be a state unitarily evolving via a Hamiltonian H. Denote the infinite-time average of σ(t) by σ . Assuming that we can make the density of states approximation, meaning we replace E by dE n(E), where n(E) is the density of states, we get where · T denotes the time average over [0, T ], and we have constants c 1 = e √ π/2 and c 2 = 4 √ π. Also, n max = max E n(E), and the effective dimension of the state σ(t) is defined by 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 → ∞, equilibration is guaranteed if The effective dimension d eff measures how spread out over energy levels the initial state is. If d 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: The main task now is to use this single-particle equilibration result to find timescales for N particle systems.

A. 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 singleparticle picture via equation (14), we want to bound where the third line follows from concavity of the square root. The last line follows from the result for singleparticle equilibration from the last section, namely equation (38). To see whether equilibration occurs at all, let T → ∞ to get the infinite-time average. And so we must estimate 1 d eff . Defining N E to be the operator that counts the number of particles in energy level E, we get where n E = tr[ρN E ]. Getting the second line used equation (6) from theorem 1. So we see that 1 d 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 σ(t) is an equal mixture of N orthogonal states. In that case, 1/d eff ≤ n d N 2 E n E ≤ n d /N , where n d is the degeneracy of the most degenerate energy level. As a result, 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 σ(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 ∝ V ∝ 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 max ∼ d 3 for a system whose energy levels go like E n ∝ 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 behaviour. For example, when E n ∝ n 2 , one gets n max ∼ 1.
Putting this all together, if equilibration occurs, the timescale is typically for some positive integer k. This is far better than the bounds of [6][7][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][16][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][16][17] for equilibration of local Hamiltonian systems. Nevertheless, [15][16][17] do say something about equilibration timescales of fully interacting models, which is very interesting.

B. 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 † (|φ )a(|φ ), where |φ is a singleparticle state localized on at most l sites (which need not be near each other). Then we have where d is the dimension of the corresponding singleparticle Hilbert space, and we have constants c 2 = 4 √ π and c 1 = (e √ π/2). Also, n max = 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 tr[ρ(t)a † x a y ], where x and 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 where we can express tr[ρ(t)a † x a y ] in terms of single-mode densities. And so via the triangle inequality, we can upper bound the time average of |tr[ρ(t)a † x a y ] − tr[ρ(t)a † x a y ]| using corollary 1.
Interestingly, these results apply to a vast range of initial states ρ(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.
(48) where d is the dimension of the corresponding singleparticle Hilbert space, and we have constants c 2 = 4 √ π and c 1 = (e √ π/2). Also, n max = 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 corollary 2, the upper bound for the equilibration timescale scales like T eq ∝ n max . So it remains to estimate n max .
In appendix E 1, we show that for these lattice models, we can effectively take n max ∝ V , where V is the number of lattice sites. Therefore, we get In particular, for one dimensional systems, we get T eq ∝ 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 eq ∝ 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].

VII. DISCUSSION AND OUTLOOK.
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].

Acknowledgments
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 1. Take a state ρ(t) = U (t)ρ(0)U † (t) of N non-interacting bosons or fermions and a measurement operator counting the number of particles in some modes Then there exist orthonormal single-particle states |ψ α (t) , evolving via the corresponding single-particle Hamiltonian, such that where n α are occupation numbers adding up to N , P = There exists a complete set of independent annihilation operators a α = a(|ψ α ), such that To see this, take any complete set of independent annihilation operators d α . The matrix tr ρ d † α d β = C αβ is Hermitian. So there exists a unitary U αβ , such that in terms of tr ρ(0) a † α a β is diagonal. So it makes sense to work with the a α , which determine the states |ψ α (0) = |ψ α mentioned in the statement of the theorem via a α = a(|ψ α ). Next, expand where c i,α are complex numbers. Then where n α is the number of particles in mode α. Next, we use where P = i |φ i φ i |. For a state like a † 1 ...a † N |0 , the first N occupation numbers n α are one, so we get To account for the dependence on time in the formula, we use a(|ψ α (t) ) =a α (t) = U (t)a α U † (t) and The end result is

Appendix B: Fluctuations
It will be useful to prove the following lemma.
Lemma 3. Let a † α = a † (|ψ α ) be a complete set of creation operators with {a α , a † β } = δ αβ . And suppose we have a state a † 1 ...a † N |0 with corresponding density operator ρ. And take a measurement operator counting the number of particles in some modes Each term can be written as which holds for bosons or fermions. To make sense of this, we write where c i,α are complex numbers. It will turn out below that only terms with α ∈ {1, ..., N } = V contribute. Now we use the identity For fermions, we use the identity Then the second term in equation (B9) becomes which is positive. Therefore, Putting everything together gives For bosons, the result is similar, but the identity in equation (B8) is replaced by Because of this, following similar steps to those used to get equation (B9), we get The extra term arising in our expression for tr ρM 2 is where P = i |φ i φ i | and we used rank(Q) = N .
A corollary of this is that for pure states of bosons or fermions of the form a † 1 ...a † N |0 , the fluctuations satisfy 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 σ 2 M = O(N ). To say something more general about the fluctuations in fermion systems, we can also prove the following result.
Theorem 4. 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 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 M , the time-average observable in the Heisenberg picture, is proportional to the total number operator.
Proof. Key to this result is the following inequality 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 As ρ is the infinite-time average of ρ(t), it follows that where p E is a normalized probability distribution and ω 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 . Now, the support of ω E is spanned by the states with E i1 + ... + E i N = E. Therefore, given two different configurations with energy E, {i 1 , ..., i N } = {j 1 , ..., j N }, because single-particle energy levels are non-degenerate. This implies where i is short for Similarly, it is a consequence of (single-particle) non-degenerate energy gaps that Now we can use lemma 3 to upper bound Next, we use tr[ωM ] = m independent of the state ω for fixed total particle number N , when ω is a time-average state, a special case of which is ω = P ( i). So Finally, using the fact that the expectation value of M is bounded above by N on the N particle subspace leads to and therefore

(C1)
And from equation (18), we have Evaluating the inner products, we get where (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 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 (σ(t), σ ) becomes Furthermore, using sin(rπ/2) = (−1) (r−1)/2 , which holds for odd r, we get for n = 2k. Next, we can find a bound for the equilibration time. First, we make the substitution n = 2k + l, noting that l > −2k since n > 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 ∈ S, where S contains all odd integers from −K to K. To quantify the resulting error, we use The sum of all terms with l / ∈ S can be bounded above So neglecting terms corresponding to l / ∈ S results in an upper bound for D P (σ(t), σ ) of which follows from the triangle inequality and | cos(x)| ≤ 1. Next, as f (2k + l, 2k) is awkward to work with, we use where we used the triangle inequality and the fact that 1/(4k + l) < 1/(2k), since l > −2k. This allows us to upper bound D P (σ(t), σ ) by 2 N l∈S N k=1 cos (4kl + l 2 )νt 1 π 2 l 2 + µ, where µ = 16 To get this, we used the triangle inequality and the inequality R r=1 1/r < ln(R) + 1. Using the triangle inequality again, we get In the third line we used the identity [33] N −1 k=0 cos(αk + φ) = sin[N α/2] cos[(N − 1)α/2 + φ] sin(α/2) .

(C15)
Let us look at each term in the sum in the last line of equation (C14) separately. They have period π 2lν , so we need only focus on this interval to find the time average of D P (σ(t), σ ).
When 2lνt is close to 0 or π, the sin[2lνt] in the denominator in the last line in equation (C14) is small. So for t such that 2lνt ∈ [0, 1 N a ) or 2lνt ∈ (π − 1 N a , π], where a is a small constant we will choose at the end, we bound When 2lνt ∈ [ 1 N a , π 2 ], we can use the inequality sin(x) ≥ 2x/π for all x ∈ [0, π/2] to get To find the time average of D P (σ(t), σ ), we use the fact that the average over [ π 4lν , π 2lν ] is the same as that over [0, π 4lν ] by symmetry. So we need only average each term over [0, π 4lν ]. The result is Plugging this into equation (C14), we get where we used l∈S 1/l 2 ≤ 2 ∞ l=1 1/l 2 = π 2 /3. If we choose K = N , then we see that 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.
In figure 2 we saw that D P (σ(t), σ ) 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 = 1 2N aν into the bound in equation (C17), gives the bound Here we choose a to be a small constant such that the distinguishability is small at t = 1 2N aν . Then the equilibration time is bounded above by

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 where J is the set of three-component vectors with components in {1, ..., J max }, so we have N = J 3 max . And let |ψ j be the energy eigenstate for a particle in the left half of a box labelled by j analogous to |ψ k 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 where P = P x ⊗ 1 1 y ⊗ 1 1 z is the projector onto the left half of the box. Here σ x (t) is the reduced state of the system on H x , the Hilbert space corresponding to the x degrees of freedom. We have where |ψ j 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 eq = O(1/J max ) = O(1/N 1/3 ).
Appendix D: Single-particle equilibration 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 5. Take a finite dimensional system evolving via a time independent Hamiltonian in the state σ(t). For any operator A and time T > 0 where c 1 = e √ π 2 and G α = E i − E j denote the non-zero energy gaps, so E i = E j . Also, where P E is the projector onto the energy eigenspace corresponding to energy E.
Proof. First, we will take σ(t) to be pure, extending the result to mixed states at the end. Because σ(t) = |ψ(t) ψ(t)| is pure, we can choose an eigenbasis of H where |ψ(t) only overlaps with a single eigenstate |n for each energy level. This means that degenerate energy levels will not cause any problems. The state at time t is where c n = n|ψ(0) . The time-average state is σ = n |c n | 2 |n n|, and the effective dimension is given by 1/d eff = n |c n | 4 = tr σ 2 .
Using the notation A ij = i|A|j , we have To make our expressions more concise, we denote nonzero energy gaps by G β = E i − E j , with β = (i, j), where i = j. We also define the vector and the Hermitian matrix As tr A † B defines an inner product for operators, we may use the Cauchy-Schwartz inequality. We can also use the inequality tr[P Q] ≤ P tr[Q], which holds for P and Q positive operators. Then we get Next, we can use the inequality for matrix norms [34] M ≤ |||M ||| ∞ |||M ||| 1 where |||M ||| 1 and |||M ||| ∞ are the column and row matrix norms respectively. The second line holds because M is hermitian, implying |||M ||| ∞ = |||M ||| 1 . Our next task is to deal with M αβ = e i(Gα−G β )t 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) ≥ 1/T on the interval [0, T ].
We will choose the weight to be a Gaussian. Then for any positive g(t), we have Here, e = exp (1). Employing this weighted averaging, the matrix elements of M are So we get where c 1 = e √ π/2. From equation (D10), we get Finally, we arrive at 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 H S . Then we can define a pure state |ψ(0) on H S ⊗H A , with dim(H S ) = dim(H A ), with the property that the reduced state on the first system is tr A [|ψ(0) ψ(0)|] = σ(0). We recover the original evolution σ(t) of the first system by evolving |ψ(t) under the joint Hamiltonian H ⊗ 1 1. Crucially, the expectation value of any operator A on the state σ(t) is the same as the expectation value of A ⊗ 1 1 on the purified state on the joint system. Also, A = A ⊗ 1 1 . Finally, the effective dimension of the purified system is equal to the effective dimension of the original system, which can be seen from from tr[P E σ(0)] = tr[P E ⊗ 1 1|ψ(0) ψ(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 α = G β , and evaluate the sum by making the density of states approximation. We make the replacement E = dE n(E), where n(E) is the density of states, so that where n max = max E n(E) and d is the dimension of the particle's state space. We also used dE n(E) = d to get the last line.
We define D G to be the maximum number of gaps G α of the same size. In other words, So finally we get It is good to point out that the density of states approximation misses degenerate energy gaps because they are measure zero.

Appendix E: Free lattice models
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 † (|φ )a(|φ ), which counts the number of particles in the state |φ . By applying equation (6) The trick now is to switch to the Heisenberg picture.
Because of this, we can think of |φ(−t) φ(−t)| = σ(−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 Π = N j=1 n j ψ j ≤ 1 1. (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.
Applying equation (38) in theorem 2, we get The first task is to estimate 1/d eff . Because the measurement operator M = a † (|φ(0) )a(|φ(0) ) is local, the state |φ(0) φ(0)| is localized. But for free lattice models, the energy eigenstates are spread out over the whole lattice: they often have the form |φ( p) | p , where | p is the momentum state with momentum p, and |φ( p) is the state of the extra degree of freedom (spin or particle type, for example). In that case, given an energy eigenstate |E we get | E|φ(0) | 2 ≤ l/V , where the factor of l appears because |φ(0) is spread over at most l lattice sites. This implies 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 where c 3 = c 1 n 2 d s 2 . The last task is to estimate n max = max E n(E).

Density of states for lattice models
For lattice models, estimating max E n(E) causes problems because n(E) often diverges. Fortunately, we can truncate to a slightly smaller energy range, such that 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 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) = cos(p), where p = 2πk/L, with k ∈ {1, ..., L}. (In making the density of states approximation, we assume that L is large but finite.) The density of states is n(E) = L π dp dE = L π At E = ±1 this is infinite. However, we can truncate the state, neglecting all terms with E ∈ [−1, − cos(p 0 )) ∪ (cos(p 0 ), 1], where we choose a fixed p 0 = 2πk 0 /L > 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 0 to be the projector onto the subspace corresponding to the omitted energy range, we get |φ − (1 1 − P 0 )|φ 2 2 = 1 − φ|(1 1 − P 0 )|φ where we used | E|φ(0) | 2 ≤ l/V from the previous section to get the last line. And N 0 is the number of states with energy in the excluded set. Next, we can use that N 0 = 4k 0 = 2p 0 L/π to get So the error in approximating |φ by a state restricted to the smaller energy range can be made small by choosing p 0 to be small. Furthermore, we have that max E n(E) = L π 1 |sin(p 0 )| .
Crucially, p 0 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 = 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 max E n(E) ∝ V , where V is the number of lattice sites.