Typical fast thermalization processes in closed many-body systems

The lack of knowledge about the detailed many-particle motion on the microscopic scale is a key issue in any theoretical description of a macroscopic experiment. For systems at or close to thermal equilibrium, statistical mechanics provides a very successful general framework to cope with this problem. However, far from equilibrium, only very few quantitative and comparably universal results are known. Here a quantum mechanical prediction of this type is derived and verified against various experimental and numerical data from the literature. It quantitatively describes the entire temporal relaxation towards thermal equilibrium for a large class (in a mathematically precisely defined sense) of closed many-body systems, whose initial state may be arbitrarily far from equilibrium.

I n a macroscopic object, which is spatially confined and unperturbed by the rest of the world, every single atom exhibits an essentially unpredictable, chaotic motion ad infinitum, yet the system as a whole seems to approach in a predictable and often relatively simple manner some steady equilibrium state. Paradigmatic examples are compound systems, parts of which are initially hotter than others, or a simple gas in a box, streaming through a little hole into an empty second box. While such equilibration and thermalization phenomena are omnipresent in daily life and extensively observed in experiments, they entail some very challenging fundamental questions: why are the macroscopic phenomena reproducible though the microscopic details are irreproducible in any real experiment? How can the irreversible tendency towards macroscopic equilibrium be reconciled with the basic laws of physics, implying a perpetual and essentially reversible motion on the microscopic level?
Here we will further extend these findings in two essential respects: instead of upper bounds for some suitably defined characteristic timescale, as in refs 49-51, the entire temporal relaxation will be approximated in the form of an equality. As an even more decisive generalization of ref. 49-51, we will admit largely arbitrary observables. Finally, and actually for the first time within the realm of the above-mentioned analytical approaches 1,8-12,38-51 , we will compare our predictions with various experimental as well as numerical data from the literature. In fact, most of those data have not been quantitatively explained by any other analytical theory before. Adopting a 'typicality approach' similar in spirit to random matrix theory [9][10][11][12] , our result covers the vast majority (in a suitably defined mathematical sense) of initial conditions, observables and system Hamiltonians.
On the other hand, many commonly considered observables and initial conditions actually seem to be rather special in that they are close to or governed by a hidden conserved quantity and therefore thermalize 'untypically slowly'.

Results
Set-up. Employing textbook quantum mechanics, we consider time-independent Hamiltonians H with eigenvalues E n and eigenvectors |ni on a Hilbert space H of large (but finite) dimensionality D ) 1. As usual, system states (pure or mixed) are described by density operators r : H ! H and observables by Hermitian operators A : H ! H with matrix elements r mn : ¼ hm|r|ni and A mn : ¼ hm|A|ni, respectively. Expectation values are given by hAi r : ¼ Tr{rA} and the time evolution by The main examples are closed many-body systems with a macroscopically well-defined energy, that is, all relevant eigenvalues E 1 , ..., E D are contained in some microcanonical energy window [E À DE, E], where DE is small on the macroscopic but large on the microscopic scale. For systems with f ) 1 degrees of freedom, D is then exponentially large in f (refs 9,41). Accordingly, the relevant Hilbert space H is spanned by the eigenvectors n j i f g D n¼1 and is sometimes also named energy shell or active Hilbert space, see, for example, refs 8-12 and Supplementary Note 2 for more details.
Analytical results. Our main players are the three Hermitian operators H (Hamiltonian), A (observable) and r(0) (initial state), each with its own eigenvalues (spectrum) and eigenvectors (basis of H). In the following, the three spectra will be considered as arbitrary but fixed, while the eigenbases will be randomly varied relatively to each other. More precisely, all unitary transformations U : H ! H between the eigenbases of H and A are considered as equally likely (Haar distributed [8][9][10][11] ), while the basis of r(0) relatively to that of A is arbitrary but fixed. (Equivalently, we could let 'rotate' H relatively to r(0) while keeping A fixed relatively to r(0).) In particular, the initial expectation value hAi r(0) can be chosen arbitrary but then remains fixed (U independent). It is only for times t40 that the randomness of the unitary U also randomizes (via H) the further temporal evolution of r(t) and thus of hAi r(t) .
The basic idea behind this randomization of U is akin to random matrix theory 9-12 , namely, to derive an approximation for hAi r(t) , which applies to the overwhelming majority of all those randomly sampled U's, hence it typically should apply also to the particular (non-random) U of the actual system of interest. A more detailed justification of this 'typicality approach' will be provided in the section 'Typicality of thermalization'.
Since A mn refers to the basis of H, these matrix elements depend on U, and likewise for r mn (0) (the explicit formulae are provided in the section 'Basic matrices'). Indicating averages over U by the symbol [?] U and exploiting that all basis transformations U are equally likely, it follows for symmetry reasons that [r nn (0)A nn ] U must be independent of n. Likewise, [r mn (0)A nm ] U must be independent of m and n for all man. We thus can conclude that for any n and that for any man for arbitrary n. Likewise, equation (3) yields for arbitrary man. Upon separately averaging in equation (1), the summands with m ¼ n and those with man over U, and then exploiting equations (4) and (5) one readily finds that FðtÞ : where f(t) is the Fourier transform of the spectral density from ref. 46 (see also refs 51-53) The following results can be derived in principle along similar lines (symmetry arguments being one key ingredient), but since the actual details are quite tedious, they are postponed to Methods. As a first result, one obtains where r mc : ¼ I/D is the microcanonical density operator and I the identity on H. As a second result, one finds for the statistical fluctuations the estimate for arbitrary t, where D A is the range of A, that is, the difference between the largest and smallest eigenvalues of A. Since averaging over U and integrating over t are commuting operations, equation (11) implies that 1 for arbitrary t 2 4t 1 . Considering t in equation (11) as arbitrary but fixed, equation (10) and D ) 1 imply (obviously or by exploiting Chebyshev's inequality 1,9,40,43,45 ) that hAi r(t) is practically indistinguishable from the average in equation (6) for the vast majority of all unitaries U. Indeed, the fraction (normalized Haar measure) of exceptional U's is unimaginably small for typical macroscopic systems with, say, fE10 23 degrees of freedom, since D in equation (11) is exponentially large in f (see below equation (1)). Likewise, considering an arbitrary but fixed time interval [t 1 ,t 2 ] in equation (12), it follows for all but a tiny fraction of U's that the time average over x 2 (t) on the left-hand side of equation (12) must be unimaginably small, and hence also the integrand x 2 (t) itself must be exceedingly small for the overwhelming majority of all tA[t 1 ,t 2 ]. Accordingly, hAi r(t) must remain extremely close to equation (6) simultaneously for all those tA[t 1 ,t 2 ].
Due to equation (9) and D ) 1, we furthermore can safely approximate A h i r av in equation (6) by A h i r mc . Altogether, we thus can conclude that in very good approximation for the vast majority of unitaries U and times t. As detailed in Methods, the neglected corrections in equation (13) consist of a systematic (U independent) part, which is bounded in modulus by D A /(D 2 À 1) for all t, and a random (U dependent) part (namely x(t)), whose typical order of magnitude is D A ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Tr r 2 ð0Þ f g=D p (for most U and t, cf. equations (11) and (12)), that is, x(t) is dominating by far (note that 1 ! Tr r 2 ð0Þ f g! Tr r 2 mc È É ¼ 1=D). Moreover, the correlations of x(t) decay on timescales comparable to those governing F(t).
These are our main formal results. In the rest of the paper, we discuss their physical content.
Basic properties of F(t). Equation (8) implies that f(0) ¼ 1, f( À t) ¼ f*(t) and |f(t)|r1. With equation (7) and D ) 1, it follows that in very good approximation and thus Fð0Þ ¼ 1; 0 FðtÞ 1; Fð À tÞ ¼ FðtÞ: Indicating averages over all tZ0 by an overbar, one can infer from equations (8) and (14) that where k labels the eigenspaces of H with mutually different eigenvalues and d k denotes their dimensions. Since P k d k ¼ D, we thus obtain FðtÞ max k ðd k =DÞ. Excluding extremely large multiplicities (degeneracies) of energy eigenvalues, it follows that the time average FðtÞ is negligibly small, and hence 1,9,40,43,45 that F(t) itself must be negligibly small for the overwhelming majority of all sufficiently large t, symbolically indicated as Fðt ! 1Þ*0: ð16Þ Note that there still exist arbitrarily large exceptional t's owing to the quasi-periodicity of f(t) implied by equation (8). We also emphasize that our main result (13) itself admits arbitrary degeneracies of H.
As an example, we focus on the microcanonical set-up introduced below equation (1) and on not too large times, so that equation (8) is well approximated by where r(x) represents the (smoothened and normalized) density of energy levels E n in the vicinity of the reference energy x. If the level density is constant throughout the energy window [E À DE, E], we thus obtain with equation (14) Next, we recall Boltzmann's entropy formula S(x) ¼ k B ln (O(x)), where O(x) counts the number of E n 's below x and k B is Boltzmann's constant. Hence, O 0 (x) must be proportional to the level density r(x) from above. Furthermore, T: ¼ 1/S 0 (E) is the usual microcanonical temperature of a system with energy E at thermal equilibrium. A straightforward expansion then yields the approximation rðE À yÞ ¼ ce À y=k B T for yZ0, where c is fixed via The omitted higher-order terms are safely negligible for all yZ0 and systems with f ) 1 degrees of freedom, see also ref. 53. With equations (14) and (17), one thus finds where a :¼ e À DE=kBT . For DE ( k B T, one recovers equation (18), and for DE ) k B T, one obtains Typicality of thermalization. Equations (13) and (16) imply thermalization in the sense that the expectation value hAi r(t) becomes (for most U) practically indistinguishable from the microcanonical average A h i r mc for the overwhelming majority of all sufficiently large t. Exceptional t's are, for instance, due to quantum revivals, which, in turn, are apparently closely related to the quasi-periodicities of F(t).
The usual time inversion invariance on the fundamental, microscopic level 7 is maintained by equation (13) due to equation (15). Surprisingly, and in accordance with the second law of thermodynamics, the latter symmetry persists even if it is broken in the microscopic quantum dynamics, for example, by an external magnetic field.
By propagating r(0) backward in time (with respect to one particular U) and taking the result as new initial state, one may easily tailor 41 examples of the very rare U's and t's, which notably deviate from the typical behaviour (13). Equivalently, one may back-propagate A instead of r(0) (Heisenberg picture).
Note that S and T were introduced below equation (18) not in the sense of associating some entropy and temperature to the non-equilibrium states r(t) but rather as a convenient levelcounting tool. However, we now can identify them a posteriori with the pertinent entropy and temperature after thermalization.
The randomization via U (see the section 'Analytical results') can be viewed in two ways: either one considers r(0), A and the spectrum of H as arbitrary but fixed, while the eigenbasis of H is sampled from a uniform distribution (Haar measure). Or one considers H and the spectra of r(0) and A as arbitrary but fixed and randomizes the eigenvectors of A and r(0). In doing so, a key point is that the relative orientation of the eigenbases of r(0) and A can be chosen arbitrarily but then is kept fixed. Indeed, it is well known 12,49 that for 'most' such orientations the expectation values hAi r(0) and A h i r mc are practically indistinguishable, that is, an initial hAi r(0) far from equilibrium requires a careful fine-tuning of r(0) relatively to A.
In reality, there is usually nothing random in the actual physical systems one has in mind. Hence, results such as equation (13), which (approximately) apply to the overwhelming majority of unitaries U, should be physically interpreted according to the common lore of random matrix theory 9,10,12 , namely, as to apply practically for sure to a concrete system under consideration, unless there are particular reasons to the contrary.
Such reasons arise, for instance, when A is known to be a conserved quantity, implying a common eigenbasis of A and H, that is, the basis transformations U must indeed be very special. Furthermore, this non-typicality is structurally stable against sufficiently small perturbations of A and/or H so that the eigenvectors remain 'almost aligned' (each eigenvector of A mainly overlaps with one or a few eigenvectors of H), and hence A remains 'almost conserved' (almost commuting with H).
Analogous non-typical U's are expected when r(0) is known to be (almost) conserved (commuting with H).
Further well-known exceptions are integrable systems, for which thermalization in the above sense may be absent for certain r(0) and A 4,32 (but not for others 22 ), systems exhibiting many-body localization 34,36

or trivial cases with non-interacting subsystems (Supplementary Note 2).
Our present focus is different: taking thermalization for granted is the temporal relaxation well approximated by equation (13)?
Typical fast relaxation and prethermalization. Equation (20) is governed by the Boltzmann time t B : ¼ h/k B T, amounting to t B E10 À 13 s at room temperature. Equation (19) gives rise to comparably short timescales, unless the temperature is exceedingly low or the energy window DE is unusually small. Such relaxation times are much shorter than commonly observed in real systems 46,[49][50][51] . Moreover, the temporal decay is typically non-exponential (see, for example, equations (18)-(20)), again in contrast to the usual findings.
This seems to imply that typical experiments correspond to non-typical unitaries U. Plausible explanations are as follows: to begin with, the above-predicted typical relaxation times are so short that they simply could not be observed in most experiments. Second (or as a consequence), the usual initial conditions and/or observables are indeed quite 'special' with respect to the prominent role of almost conserved quantities (see previous section), in particular, 'local descendants' of globally conserved quantities such as energy, charge, particle numbers and so on: examples are the amount of energy, charge and so on within some subdomain of the total system or, more generally, local densities, whose content within a given volume can only change via transport currents through the boundaries of that volume. As a consequence, the global relaxation process becomes 'unusually slow' if the densities between macroscopically separated places need to equilibrate (small surface-to-volume ratio) or if there exists a natural 'bottleneck' for their exchange (weakly interacting subsystems).
Put differently, our present theory is meant to describe the very rapid relaxation towards local equilibrium, but not any subsequent global equilibration. Only if there exists a clear-cut timescale separation between these two relaxation steps (or if there is no second step at all), can we hope to quantitatively capture the first step by our results. Conversely, the timescale separation usually admits some Markovian approximation for the second step, yielding an exponential decay, whose timescale still depends on many details of the system. Natural further generalizations include the closely related concepts of hindered equilibrium, quasi-equilibrium (metastability) and, above all, prethermalization 29,54,55 , referring, for example, to a fast partial thermalization within a certain subset of modes, (quasi-)particles or other generalized degrees of freedom. (Like in ref. 54, we do not adopt here the additional requirement 55 that the almost conserved quantities originate from a weak perturbation of an integrable system.) In short, our working hypothesis is that the theory (13) describes the temporal relaxation of hAi r(t) for any given pair (r(0), A) unless one of them is exceptionally close to or in some other way slowed down by an (almost) conserved quantity.
Comparison with experimental results. We focus on experiments in closed many-body systems in accordance with the above general requirements. In comparing them with our theory (13), we furthermore assume that the (pre-)thermalized system occupies a microcanonical energy window with some (effective) ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms10821 temperature T and DE ) k B T, so that equation (20) applies. Finally, the asymptotic values hAi r(0) and A h i r mc in equation (13) are either obvious or will be estimated from the measurements, hence no further knowledge about the often quite involved details of the experimental observables will be needed! Figure 1 demonstrates the very good agreement of the theory with the rapid initial prethermalization of a coherently split Bose gas, observed by the Schmiedmayer group in ref. 29.
In Fig. 2, the theory is compared with the pump-probe experiment by the Bigot group from ref. 56. The finite widths of the pump and the probe laser pulses are roughly accounted for by convoluting equation (13) with a Gaussian of 35 fs FWHM (full width at half maximum). In ref. 56, the FWHM of the pump pulse is estimated as 20 fs and the combined FWHM for both pulses as 22 fs, implying a FWHM of 9 fs for the probe pulse. The latter value seem quite optimistic to us. A second 'excuse' for our slightly larger FWHM value of 35 fs is that the tails of the experimental pulse shape may be considerably broader than those of a Gaussian with the same FWHM (see, for example, Fig. 2c in the Supplementary Material of ref. 57). Finally, the convolution of equation (13) with a Gaussian represents a rather poor 'effective description' in the first place: our entire theoretical approach becomes strictly speaking invalid when the duration of the perturbation becomes comparable to the thermalization time.
A similar comparison with the pump-probe experiments from ref. 58 is presented in Fig. 3. As before, we adopted a slightly larger FWHM of 100 fs than the estimate of 76 fs in ref. 58. Due to the above-mentioned fundamental limitations of our theory for such rather large FWHM values, the temperatures adopted in Fig. 3 should still be considered as quite crude estimates. Apart from that, Fig. 3 nicely confirms the predicted temperature dependence from equation (20).
We close with three remarks: first, refs 56,58 also implicitly confirm our prediction that the essential temporal relaxation (encapsulated by F(t) in equation (13)) is generically the same for different observables. Second, similar pump-probe experiments abound in the literature, but usually the pulse widths are too large for our purposes. Third, the temporal relaxation in Figs 1-3 has also been investigated numerically, but closed analytical results have not been available before 29,58 .
Comparison with numerical results. Figure 4 illustrates the very good agreement of our theory with Rigol's numerical findings from ref. 32, both for an integrable and an non-integrable example. A similar agreement is found for all other parameters and also for an analogous hardcore boson model examined in refs 31,32. On the other hand, a second observable considered in ref. 32, deriving from the momentum distribution function, exhibits in all cases a significantly slower and also qualitatively different temporal relaxation. According to the discussion in the section 'Typical fast relaxation and prethermalization', it is quite plausible that the latter observable is indeed 'non-typical' in view of the fact that it represents a conserved quantity for fermions In Fig. 5, we compare our theory with the simulations of a different one-dimensional electron model from ref. 59. In doing so, the pertinent temperature T has been estimated as follows: the textbook Sommerfeld expansion for N electrons in a one- , where E is their total energy, E 0 ¼ (1/3)NE F the ground-state energy, E F ¼ (p' N/gL) 2 /2m the Fermi energy, L the box length, m the electron mass and g: ¼ 2s þ 1 ¼ 2 (s ¼ 1/2 for electrons). Assuming that the pulse acts solely on the small well implies N ¼ 16, LC15 nm (ref. 59), and E À E 0 C0.045 eV (see Fig. 8a in ref. 59). Altogether, we thus obtain TC170 K.
The remnant 'fluctuations' of the numerical data in Figs 4 and 5 can be readily explained as finite particle number effects (see Fig. 4 in ref. 32 and Fig. 10 in ref. 59), and their temporal correlations are as predicted below equation (13). The seemingly rather strong fluctuations in Fig. 5 are a fallacy since the systematic changes themselves are very small.
Next, we turn to the numerical findings for a qubit in contact with a spin bath by the Trauzettel group from ref. 60. The agreement with our theory in Fig. 6 is as good as it possibly can be for such a rather small dimensionality of D ¼ 2 7 . Indeed, the remaining differences nicely confirm the predictions below equation (13), regarding both their typical order of magnitude   (13) and (20) with T ¼ 5 nK. The pertinent effective temperature has also been roughly estimated in ref. 29 (see Fig. 2b therein) and is still compatible with our present fit T ¼ 5 nK. As discussed at the end of the section 'Typical fast relaxation and prethermalization', the depicted prethermalization is followed by a much slower, global thermalization 29 , which is omitted in the present figure.
Our final example is Bartsch and Gemmer's random matrix model from ref. 17. Referring to the notation and definitions in the caption of Fig. 7, one readily sees that the considered observable A is a conserved quantity for the unperturbed Hamiltonian (l ¼ 0). In agreement with our discussion in the section 'Typical fast relaxation and prethermalization', A is therefore still 'almost conserved' for small l and indeed exhibits a slow, exponential decay towards A h i r mc ¼ 0 (see Fig. 1a in ref. 17). Upon increasing l, one recovers the much faster, non-exponential decay of our present theory (see Fig. 1b in ref. 17). Unfortunately, the l-value 1.77 Â 10 À 3 from Fig. 1b of ref. 17 is still somewhat too small and the eigenvalues E 1 , ..., E 6,000 are not any more available (I asked the authors). Therefore, we repeated the numerics from ref. 17 on our own for l ¼ 7 Â 10 À 3 . The resulting agreement with equation (13) in Fig. 7 is very good, and the temporal correlations of the deviations as well as their typical order of magnitude D A ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Trfr 2 ð0Þg=D p ¼ 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1=6; 000 p ' 0:03 are as predicted below equation (13).
We close with two remarks: first, there is no fit parameter in any of the above examples apart from hAi r(0) in Fig. 4 (13) and (20) with temperatures as indicated and convoluted with a Gaussian of 100 fs FWHM (see also main text). The conversion of a given fluence into a temperature change of the electron gas is not obvious. In particular, the estimates provided in ref. 58 seem not very reliable to us: first of all, Fig. 6 in ref. 58 indicates a temperature of ca. 250 K at four different time points B200 fs before the pump pulse, while the actual temperature of the unperturbed system is known to be 130 K. Second, the temperature error bars in Fig. 6b Supplementary  Fig.3b of ref 66 should coincide, while their actual agreement is only moderately better than for the 'bare' curves in Supplementary Fig. 3a. For all these reasons, we used the temperature as a fit parameter in the present figure.
A quantum quench generates an initial pure state out of equilibrium, whose energy corresponds to that of a canonical ensemble with temperature T ¼ 2. As detailed in ref. 32, the considered observable dN k (t) is a dimensionless descendant of the density-density structure factor. The depicted data are from Fig. 1g,j of ref. 32. Lines: theoretical predictions (13) and (20) with T ¼ 2.  (13) and (20), exploiting the estimate T ¼ 170 K from the main text, and neglecting the finite temporal width (20 fs) of the pulse. As in Figs 1-3, we are actually dealing with a prethermalization process within the smaller well. The subsequent global thermalization is much slower due to the high barrier between the wells. Considering that A h i r mc is the only remaining fit parameter in the theory from equations (13) and (20), the agreement with the simulations is remarkably good. In particular, the two very different observables CT(t) and DP(t) are indeed governed by the same F(t), as predicted by equations (13) and (20). model in Fig. 4, one may question whether the considered system exhibits thermalization in the first place, as is tacitly assumed in equation (13). In Supplementary Note 2, we argue that equation (13) indeed is expected to still remain valid in such cases if A h i r mc is replaced by the pertinent non-thermal long-time asymptotics (which, in turn, is estimated from the numerical data in Fig. 4).

Discussion
Our main result (13) implies thermalization in the sense that a generic non-equilibrium system with a macroscopically welldefined energy becomes practically indistinguishable from the corresponding microcanonical ensemble for the overwhelming majority of all sufficiently late times. Apart from the concrete initial and long-time expectation values (that is, hAi r(0) and A h i r mc in equation (13)), the temporal relaxation (that is, F(t) in equation (13)) depends only on the spectrum of the Hamiltonian within the pertinent interval of non-negligibly populated energy eigenstates, but not on any further details of the initial condition or the observable. This represents one of the rare instances of a general quantitative statement about systems far from equilibrium.
The theory agrees very well with a wide variety of experimental and numerical results from the literature (though none of them was originally conceived for the purpose of such a comparison). We are in fact not aware of any other quantitative analytical explanation of those data comparable to ours. Indeed, the usual paradigm to identify and then analytically quantify the main physical mechanisms seems almost hopeless here. In a sense, our present approach thus amounts to a different paradigm: there is no need of any further 'explanations', since the observed behaviour is expected with overwhelming likelihood from the very beginning, that is, unless there are special a priori reasons to the contrary.
Similarly as in refs 46,49-51, generic thermalization is found to happen extremely quickly (unless the system's energy or temperature is exceedingly low). Moreover, the temporal decay is typically non-exponential. A main prediction of our theory is that these features should in fact be very common (at least in the form of prethermalization), but often they are unmeasurably fast or they have simply not been looked for so far. Conversely, most of the usually considered observables and initial conditions are actually quite 'special', namely, exceptionally slow, 'almost conserved' quantities. A better understanding of those principally untypical but practically very common thermalization processes remains an open problem [49][50][51] .

Methods
Basic matrices. According to the section 'Analytical results', the unitary U represents the basis transformation between the eigenvectors |ni (n ¼ 1, ..., D) of the Hamiltonian H and those of the observable A. Denoting the eigenvalues of A by l n and the eigenvectors by c v j i (n ¼ 1, ..., D), the matrix elements of U are thus U nn : ¼ n j c v h i . Accordingly, the matrix elements of r(0) in the basis of H are related to those in the basis of A via where r mn :¼ c u rð0Þ j jc v h i . Similarly, the matrix elements of A satisfy and hence As announced below equation (3), we work (without loss of generality) in a reference frame (or reference basis of H) so that only H (and thus |ni) depends on U, while A and r(0) (and thus c v j i) are independent of U. Hence, r mn and l x on the right-hand side of equations (21)-(23) are independent of U.
Derivation of equation (9). As a simple first exercise, let us average equation (23) over all uniformly (Haar) distributed unitaries U, as specified in the section 'Analytical results'. Since the factors r mn l x on the right-hand side are independent of U, we are left with averages over the U matrix elements. Such averages have been evaluated repeatedly and often independently of each other in the literature, see, for example, refs 5,61-63, a key ingredient being symmetry arguments due to the invariance of the Haar measure under arbitrary unitary transformations. Particularly convenient for our present purposes is the formalism adopted by Brouwer and Beenakker, see ref. 63, and further references therein. The general structure of such averages is provided by equation (2.2) in ref. 63, reading Quoting verbatim from ref. 63, 'the summation is over all permutations P and P 0 of the numbers 1, ..., n. The coefficients V P,P 0 depend only on the cycle structure of the permutation P À 1 P 0 . Recall that each permutation of 1, ..., n has a unique factorization in disjoint cyclic permutations ('cycles') of lengths c 1 , ...., c k (where n ¼ P k j¼1 c j ). The statement that V P,P 0 depends only on the cycle structure of P À 1 P 0 means that V P,P 0 depends only on the lengths c 1 , ..., c k of the cycles in the factorization of P À 1 P 0 . One may therefore write V c1 ; :::; ck instead of V P,P 0 .' The  (13), (14) and (8). Due to the above-mentioned initial condition and the quite small dimension D ¼ 2 7 , the approximations (18)- (20) are not very well satisfied by the actual energy eigenvalues E 1 , ... , E 128 (kindly provided by the authors of ref. 60). Hence, we have evaluated F(t) in equation (13) directly via equations (14) and (8).  (13), (14) and (8). Similarly as in Fig. 6, the numerically obtained energies E 1 , ..., E 6,000 were found to satisfy equations (18)-(20) not very well, hence we have directly evaluated equations (8) and (14).
explicit numerical values of all V c1; :::; ck with nr5 are provided by the columns 'CUE' of Tables II and IV in ref. 63. Further remarks: the labels m and n in equation (24) have nothing to do with those in equation (23). Equation (24) equals zero unless m ¼ n. Every label a j must have a 'partner', that is, its value must coincide with one of the a j 's and vice versa, since otherwise the product over the Kronecker delta's d aja PðjÞ in equation (24) would be zero for all P's. Note that some a j 's may assume the same value, but then an equal number of a j 's also must assume that value. Likewise, every b j needs a 'partner' among the b j 's and vice versa.
Adopting the abbreviation and the renamings a 1 : ¼ m, a 2 : ¼ n, b 1 : ¼ m, b 2 : ¼ x and b 3 : ¼ n, equation (23) yields The connection with equation (24) is established via the identifications a 1 : ¼ a 1 , a 2 : ¼ a 2 , b 1 : ¼ b 2 and b 2 : ¼ b 3 . Therefore, if b 1 ab 2 , then the only potential 'partner' of b 1 is b 2 , and only if their values coincide, that is, b 3 ¼ b 1 , the corresponding summands may be non-zero. The same conclusion can be drawn if b 1 ¼ b 2 . We thus can rewrite equation (26) with equation (24) as There are two permutations of the numbers 1 and 2, namely, the identity and one, which exchanges 1 and 2. Denoting them as P 1 and P 2 , respectively, and observing that b j ¼ b P2ðjÞ , equation (27) can be rewritten as For l ¼ 1, the two Kronecker delta's in equation (29) both require that b 1 ¼ b 2 and hence The last equality can be verified by evaluating the trace in the eigenbasis of A, see above equation (21). In the same way, one finds that In the last equation, we exploited that Tr{r(0)} ¼ 1 and r mc : ¼ I/D, see below equation (9). Observing that the two Kronecker delta's in equation (28) equal one if k ¼ 1 or if k ¼ 2 and a 1 ¼ a 2 , the overall result is where, as usual, hAi r(0) : ¼ Tr{r(0)A} and A h i r mc :¼ Tr r mc A f g. Finally, the coefficients V Pk;Pl are evaluated as explained below equation (24): if k ¼ l, then P À 1 l P k ¼ P 1 factorizes in two cycles of lengths c 1 ¼ c 2 ¼ 1, that is, V Pk;Pl ¼ V c1;c2 ¼ V 1;1 . Likewise, if kal, then P À 1 l P k ¼ P 2 consists of one cycle with c 1 ¼ 2, that is, V Pk;Pl ¼ V 2 . Referring to columns 'CUE' and rows 'n ¼ 2' of Tables II and IV in ref. 63 yields V 1,1 ¼ 1/(D 2 À 1) and Returning to the original labels m and n in equation (25), we thus can rewrite equation (32) as As a consequence, we can infer from equations (4) and (25) that A h i r av ¼ DX nn and with equation (33) that Hence, one readily recovers equation (9). A relation remarkably similar to our present equation (9), albeit in a quite different physical context, has been previously obtained also in ref. 64 (see equation (2) therein).
Derivation of equation (11). Without any doubt, there are much faster ways to obtain equations (33) or (34). The advantage of our present way is that it can be readily adopted without any conceptual differences (albeit the actual calculations become more lengthy) to more demanding cases like see equation (10).
Derivation of equation (13). While the essential steps in deriving equation (13) have been outlined already in the main text, we still have to provide the details of the statements below equation (13): our first observation is that R 1 (t) in equation (36) amounts to the systematic (U independent) part of the omitted corrections in equation (13), and equation (41) to the bound announced below equation (13).
By means of a straightforward (but again very tedious) generalization of the calculations from the preceding subsection, one finds that where C(t, s) has the following six properties: first, C(t, s) ¼ C(s, t) ¼ C( À t, À s) for all t, s. Second, |C(t, s)|r9 for all t,s. Third, C(t, 0) ¼ 0 for all t. Fourth, C(t, s)*0 for |t À s|-N, cf. equation (16). Fifth, Cðt; sÞFðt À sÞhðA À A h i r mc Þ 2 i r mc for t, s-N. Sixth, given s, the behaviour of C(t, s) as a function of t is roughly comparable to that of F(t À s) for most t.
Though we did not explicitly evaluate the last term in equation (53), closer inspection of its general structure shows that it can be bounded in modulus by cD 2 A =D 2 for some c, which is independent of t, s, D, A, r(0) and H. Moreover, there is no indication of any fundamental structural differences in comparison with the leading and next-to-leading order terms, which we did evaluate. In other words, the last term in equation (53) is expected to satisfy properties analogous to those mentioned below equation (53). Recalling that the purity Tr{r 2 (0)} satisfies the usual bounds 1 ! Tr r 2 ð0Þ f g!Tr r 2 mc È É ¼ 1=D, we thus recover the properties of x(t) announced below equation (13).