Self-averaging in many-body quantum systems out of equilibrium: Chaotic systems
Mauro Schiulaz E. Jonathan Torres-Herrera Francisco Pérez-Bernal Lea F. Santos









    Department of Physics, Yeshiva University, New York City, New York, 10016, USAInstituto de Física, Benemérita Universidad Autónoma de Puebla,
        Apt. Postal J-48, Puebla, 72570, MexicoDep. CC. Integradas y Centro de Estudios Avanzados en Física,
        Matemáticas y Computación. Fac. CC. Experimentales, Universidad de Huelva, Huelva 21071, & Instituto Carlos I de Física Teórica y Computacional, Universidad de Granada, Granada 18071, SpainDepartment of Physics, Yeshiva University, New York City, New York, 10016, USA
        Despite its importance to experiments, numerical simulations, and the development of theoretical models, self-averaging in many-body quantum systems out of equilibrium remains underinvestigated. Usually, in the chaotic regime, self-averaging is taken for granted. The numerical and analytical results presented here force us to rethink these expectations. They demonstrate that self-averaging properties depend on the quantity and also on the time scale considered. We show analytically that the survival probability in chaotic systems is not self-averaging at any time scale, even when evolved under full random matrices. We also analyze the participation ratio, Rényi entropies, the spin autocorrelation function from experiments with cold atoms, and the connected spin-spin correlation function from experiments with ion traps. We find that self-averaging holds at short times for the quantities that are local in space, while at long times, self-averaging applies for quantities that are local in time. Various behaviors are revealed at intermediate time scales.
        Self-averaging in many-body quantum systems out of equilibrium: Chaotic systems
            Lea F. Santos
            December 30, 2023
        ===============================================================================
        
        
        
        § INTRODUCTION
        
        
        The property of self-averaging is at the heart of studies about disordered systems <cit.> and random matrices <cit.>. It holds when the distribution of the quantity of interest is peaked around its average and its relative variance goes to zero as the system size increases. For sufficiently large systems, the distribution converges to a delta function. This implies that, by increasing the system size, one can reduce the number of samples used in experiments and in statistical analysis. If the system exhibits self-averaging, its physical properties are independent of the specific realization, which allows for the construction of theoretical models to describe finite samples.  Lack of self-averaging, on the other hand, means sample to sample fluctuations even in the thermodynamic limit, so ensemble averages are needed no matter how large the system size is. In this case, scaling analyses become quite challenging.
        
        Absence of self-averaging typically happens near critical points of disordered systems <cit.>. This is the case of one-body <cit.> and many-body <cit.> systems in the vicinity of the metal-insulator transition. Using results from Anderson localization, it has been shown, for example, that the entanglement entropy is not self-averaging <cit.>. Self-averaging has also been the theme of works about spin glass <cit.>, the kinetics of domain growth <cit.>, and diffusion in disordered media <cit.>, often in comparison to ergodicity.
        
        Ergodicity refers to temporal averages <cit.> and has recently received extensive attention in studies about equilibration and thermalization of isolated quantum systems <cit.>. Self-averaging, on the other hand, is associated with averages over disorder realizations and it has got little consideration in the context of quantum systems out of equilibrium, apart from few recent works on driven systems <cit.> and studies about the two-level form factor <cit.>, which is a quantity used to study spectral properties in the time domain.
        
        The present work addresses the mostly uncharted territory of self-averaging during the evolution of interacting many-body quantum systems. The focus is on the chaotic regime, where self-averaging is usually assumed to hold. Our results invalidate these expectations. 
        
        We show analytically and confirm numerically that the survival probability evolving under full random matrices from the Gaussian orthogonal ensemble (GOE) is not self-averaging at any time scale.  We also study the survival probability in a chaotic disordered spin model of experimental interest and verify that it is nowhere self-averaging. This is a consequential result, since this spatially non-local quantity is a very common tool in studies of nonequilibrium quantum dynamics <cit.>. As we show, if self-averaging is assumed and one decreases the number of random realizations used as the system size increases, one misses essential features of the evolution of the survival probability. The same is expected for other non-self-averaging quantities.
        
        Our analysis is extended also to other non-local and local quantities evolved with both the GOE and the disordered spin model. As examples of non-local quantities, we consider the participation ratio and Rényi entropies. They measure the spread of the initial state in the many-body Hilbert space and are connected with the out-of-time ordered correlator <cit.>. As local observables, we investigate the spin autocorrelation function, which is similar to the density imbalance measured in experiments with cold atoms <cit.>, and the connected spin-spin correlation function, which is used in experiments with ion traps <cit.>. The results are rather non-trivial, being dependent on the observable and time scales, although a general picture emerges for short and long times.
        
        The two spatially local quantities are self-averaging at the short times currently accessible experimentally, which is reassuring. After equilibration, when there are only small fluctuations around infinite-time averages, the connected spin-spin correlation function, the participation ratio, and the Rényi entropies are all self-averaging, while the survival probability and the spin autocorrelation function are not. These two latter quantities are autocorrelation functions, being therefore non-local in time. 
        
        This paper is organized as follows. The definition of self-averaging, as well as the  models, initial states, and quantities investigated are presented in Sec. <ref>. Section <ref> shows the entire evolution of the average values of the observables under the GOE and the spin model. The plots display the different time scales involved in the dynamics and they serve as references for the following core sections. Section <ref> provides the analytical and numerical results for the survival probability. The other non-local quantities are studied in Sec. <ref>, and the local experimental observables are examined in Sec. <ref>. An explicit example of the consequences one may face when decreasing the number of random realizations for a quantity that is non-self-averaging is given in Sec. <ref>. Conclusions and future directions are presented in Sec. <ref>.  There are also two appendices. 
        
        
        
        
        
        
        § GENERAL DEFINITIONS
        
        
        In this section, we define the concept of self-averaging, and introduce the Hamiltonians and quantities examined in this work.
        
        
        
         §.§ Self-averaging
        
        
        Self-averaging implies that a single large system is enough to represent the whole statistical ensemble. By analyzing the sample to sample fluctuations, a quantity O is said to be self-averaging when the ratio between its variance σ^2_O and the square of its mean, that is, its relative variance <cit.>,
         R_O(t)=σ^2_O(t)/<O(t)>^2 = <O^2(t)>-<O(t)>^2/<O(t)>^2 ,
        
        goes to zero as the system size L increases. In our studies, ⟨.⟩ includes the averages over both disorder realizations and initial states taken in a narrow energy window around the middle of the spectrum. Notice that our R_O(t) is time-dependent, since we investigate whether the observable is self-averaging not only at equilibrium, but during its entire time evolution.
        
        It is common to distinguish strong self-averaging, when R_O(t) ∼ L^-1, from weak self-averaging, when R_O(t) ∼ L^-ν for  0<ν<1. In this work, we find also more extreme cases, in the sense that the relative variance of O can decrease or increase exponentially in system size. This sort of “super" self-averaging or “super" non-self-averaging behavior occurs at large times, when the dynamics of a chaotic many-body system become analogous to those of full random matrices <cit.>. At such long times, the initial state is spread over the many-body Hilbert space, which is exponentially large in L.
        
        Equation (<ref>) is the standard definition of self-averaging <cit.>. One may, of course, have observables for which <O(t)>^2 goes to zero at long times as the system size increases, which is indeed the case of all quantities studied here, except for the entropies. The question we address is whether the variance goes to zero faster than <O(t)>^2.
        
        
        
        
         §.§ Models and Initial States
        
        
        We study quantum Hamiltonians of the form
        
        H=H_0+ V,
        
        where H_0 is the integrable part of H and V is a strong perturbation that takes the system into the chaotic regime. We denote by |n⟩ the eigenstates of H_0. The eigenstates and eigenvalues of H are  |α⟩ and E_α, respectively.
        
        
        
          §.§.§ GOE model
        
        One of the models that we study is formed by GOE full random matrices of dimension D<cit.>. For this model, H_0 corresponds to the diagonal part of the matrix, while V contains the off-diagonal elements. All entries are real numbers independently drawn from a Gaussian distribution with mean value <<H_ij>>=0 and variance
        <<H_ij^2>>=
        2     i=j 
        
        1     i≠ j
        .
        
        This model is unphysical, because it assumes interactions between all degrees of freedom, but it has the advantage of allowing for analytical calculations. This is possible, because the eigenvalues of the GOE model are highly correlated and the eigenstates are normalized random vectors. This model is also relevant, because it correctly reproduces the spectral correlations and the late time dynamics of realistic models <cit.>.
        
        
        
          §.§.§ Disordered spin model
        
        
        We consider a realistic disordered spin-1/2 chain in the strong chaotic regime. It has local two-body interactions only and its Hamiltonian is given by,
        
        H_0    =    J ∑_k=1^L (h_k S_k^z +  S_k^z S_k+1^z) 
        
        V    =    J∑_k=1^L (S_k^x S_k+1^x + S_k^y S_k+1^y).
        
        In the above, ħ=1, S_k^x,y,z are spin operators on site k, L is the number of spins in the lattice,  J sets the energy scale, and periodic boundary conditions are taken. The Zeeman splittings h_i are independent random numbers uniformly distributed in [-h,h], where h is the disorder strength. The total magnetization in the z-direction is conserved. We work in the largest subspace, namely the one with zero total z-magnetization, which has dimension D=L!/(L/2)!^2. To be in the fully chaotic region, we fix h=0.75. This model has been extensively studied in the context of many-body localization, both theoretically <cit.> and experimentally <cit.>.
        
        While, for the sake of concreteness, our calculations and numerical simulations are done for this model, they can be extended to a large class of different systems <cit.>. The important elements that the model has to satisfy are to be strongly chaotic, in the sense of level statistics comparable to those of random matrix theory, and have interactions that are strictly local and two body only.
        
        
        
        
          §.§.§ Initial state and notation for R(t)
        
        
        The initial state |Ψ (0)> is an eigenstate of H_0, as often considered in experiments with cold atoms and ion traps. It has energy close to the middle of the spectrum,
        
        E_0 = ⟨Ψ(0)|H|Ψ(0)⟩ = ∑_α|c_α^0|^2 E_α∼ 0,
         
        where c_α^0=<α|Ψ (0)>. This is the region of the spectrum where the energy eigenstates are chaotic <cit.>. Such initial states are very far from equilibrium, which results in an extremely fast initial evolution under the full Hamiltonian H.
        
        In our numerical simulations, our averages are performed over 0.01D  initial states with energy closest to the center of the spectrum and over many disorder realizations. The total amount of data considered in each average is 10^4.
        
        For clarity, we refer to the relative variance obtained for GOE matrices as R^GOE_O(t) and as R^spin_O(t) the one obtained for the chaotic spin model.
        
        
        
        
        
         §.§ Quantities
        
        
        We consider both non-local quantities and local experimental observables. 
        
        
        
        
          §.§.§ Survival Probability
         
        
        The survival probability is the squared overlap between the initial state and its time evolved counterpart, 
        
        P_S(t)=|<Ψ (0)|e^-iHt|Ψ (0)>|^2 .
        
        This quantity is non-local in space and also in time, since it compares the state at time t with the state at time t=0. It has been studied in many different contexts, from the decay of unstable nuclei <cit.>, the quantum speed limit <cit.>, and the onset of power-law decays <cit.>, to quench dynamics <cit.>, ground state and excited state quantum phase transitions <cit.>, and multifractality in one-body and many-body systems <cit.>. The survival probability is related to the two-level form factor studied in <cit.>, but this one contains information about the eigenvalues only, while the survival probability contains information about the initial state also, being therefore more appropriate for studies of dynamics. 
        
        The survival probability can be written in the following useful integral representation,
        
        P_S(t)=| ∫ dE e^-iEtρ_0(E) |^2,
        
        where
        ρ_0(E)=∑_α|c_α^0|^2δ(E-E_α)
        
        is the energy distribution of the initial state, known as local density of states (LDOS). The width Γ_0 of the LDOS is related to the number of states |n⟩  that are directly coupled to |Ψ (0)> according to
        Γ_0^2 = ∑_n ≠ 0 |  ⟨ n |H| Ψ(0) ⟩ |^2.
        
        In the above, the sum runs over all states |n>, apart from |Ψ(0)>. For the GOE model, the average over initial states and disorder realizations naturally gives 
        <Γ_0^2>^GOE=D,
        
        while the sparsity of the spin Hamiltonian implies that 
        <Γ_0^2>^spin = J^2 L^2/8 (L-1)∼J^2 /8  L .
        
        This difference has important consequences for the time scales involved in the evolution of the mean value of the observables and in their self-averaging behavior. In what follows, we use the notation Γ = √(<Γ_0^2>).
        
        
        
        
          §.§.§ Inverse Participation Ratio and Rényi Entropies
         
        
        The inverse participation ratio and the Rényi entropies are non-local quantities in space, but they are local in time. In contrast to the survival probability, they compare the state at time t with all states |n⟩, not only with | Ψ(0) ⟩. 
        
        The  inverse participation ratio is defined as
        IPR(t)=∑_n|<n|e^-iHt|Ψ (0)>|^4 .
        
        It quantifies the spread of the initial many-body state in the basis of unperturbed many-body states  |n><cit.>.  At t=0, when the initial state is fully localized in this basis, IPR(0)=1. A state completely delocalized at time t has IPR(t)∼ 1/D.
        
        The second-order Rényi entropy is related to the inverse participation ratio as
        
        S(t)=-ln [ IPR(t)] ,
        
        where no partial trace of degrees of freedom is involved. While the asymptotic value of IPR(t) scales with the inverse of the size of the exponentially large Hilbert space, the maximum value of S(t) is proportional to L. The motivation to study not only IPR(t), but also its logarithm, comes from the knowledge that the logarithm cuts the tails of the distribution, therefore enhancing self-averaging properties. 
        
        The results for S(t) are equivalent to those for the Shannon entropy (or first-order Rényi entropy), which is written as
        
        Sh (t) = - ∑_n |<n|e^-iHt|Ψ (0)>|^2 ln|<n|e^-iHt|Ψ (0)>|^2 .
        
        This entropy is often used in studies of quantum chaos (see <cit.> and references therein).
        
        
        
        
          §.§.§ Spin Autocorrelation Function
         and Connected Correlation Function
         
        
        The spin autocorrelation function and the connected spin-spin correlation function are experimental quantities. They are both local in space, but only the latter is also local in time. 
        
        The spin autocorrelation function  is given by
        
        I(t)=4/L∑_k=1^L <Ψ (0) |S^z_k e^iHt S^z_k e^-iHt|Ψ (0) >.
        
        It measures the average over all sites of the proximity of the orientation of a spin k at time t to its orientation at t=0. By mapping the spin system to hardcore bosons, one finds that this quantity is analogous to the density imbalance between even and odd sites, which is measured in cold atom experiments <cit.>.
        
        The connected spin-spin correlation function is defined as
        
        C(t)   =   4/L∑_k[<Ψ(t)|S_k^zS_k+1^z|Ψ(t)>.
           -   .<Ψ(t)|S_k^z|Ψ(t)><Ψ(t)|S_k+1^z|Ψ(t)>].
        Similar to IPR(t), it quantifies how far |Ψ(t)⟩ is from the classical states |n⟩. This quantity has been measured in experiments with ion traps <cit.>. 
        
        
        
        § DYNAMICS OF MEAN VALUES
        
        
        Before studying in detail the behavior of the relative variance of the quantities above, we briefly explain how the average values, <O(t)>, change with time. We outline the main stages of the evolution and the time scales associated with these steps, so that in the following sections, we can analyze how the fluctuations behave in each of these regions.
        
        In Fig. <ref> (a) and Fig. <ref> (b), we show the evolution of the mean value of the survival probability for the GOE model and for the chaotic spin model, respectively. The entire dynamics is depicted, from the moment the system is quenched out of equilibrium to the moment a new equilibrium is reached, which happens when <P_S(t)> only fluctuates around a finite asymptotic value. This saturation point corresponds to
        < P_S> = <  lim_t→∞ P_S(t) > = <∑_α|c_α^0|^4>.
        
        
        The analytical expression for the evolution of the survival probability was obtained for large Hamiltonian matrices in <cit.>. For both the GOE and the realistic model, it is given by
        
            <P_S(t)>=1- < P_S>/D-1[D b_1^2( Γ t)-b_2(Γ t/μ D)]+ < P_S>.
        
        The first term in the equation above is determined by the shape and energy bounds of the LDOS <cit.>, which depend on the model and initial state. For GOE matrices, the shape is semicircular, up to corrections that are subdominant in 1/D, which gives <cit.>
        b^2_1(Γ t)=𝒥^2_1(2 Γ t)/Γ^2 t^2,
        
        where J_1 is the Bessel function of the first kind. For realistic chaotic many-body systems, the LDOS is Gaussian <cit.> and b_1 involves an early Gaussian decay, e^-Γ^2 t^2, and a later power-law behavior ∝ t^-2, as given in <cit.>. 
        
        The b_1 function controls the initial decay of the survival probability.  For Γ t≪  1, it leads to the universal 1 - Γ^2 t^2 behavior, where 1/Γ is the characteristic time for the depletion of the initial state. Later, oscillations emerge that decay as a power law <cit.>. The power-law behavior continues until the minimal value of <P_S(t)> is reached at t_Th. This time is referred to as Thouless time and it marks the point of the complete spread of the initial state in the many-body Hilbert space, as explained in Ref. <cit.>.
        
        Beyond t_Th, the dynamics become universal and determined by the second term in Eq. (<ref>), which is the two-level form factor,
        
        b_2(t) =
        1-2t + t ln(2 t+1)     t≤1 
        
        t ln( 2 t+12 t -1) -1     t>1
        .
         
        This function is responsible for the dip below < P_S>, which is known as correlation hole and exists only when the eigenvalues are correlated <cit.>. 
        Since the GOE and the realistic chaotic model have similar level statistics, the same equation for b_2(t) is used for both cases. In Eq. (<ref>), μ=2 for the GOE model and μ = √(2 π) for the spin model <cit.>. The b_2 function initially grows linearly and later shows a power-law behavior up to saturation, which happens at the relaxation time t_R.
        
        We therefore have four regions in time that exhibit different behaviors, as indicated in Fig. <ref> (a) and Fig. <ref> (b):
        
          *  The short time region, for t ≪ 1/Γ.
        	
          *  The power-law decay, happening for 1/Γ< t <t_Th. 
        	
          *  The interval for the correlation hole, t_Th< t <t_R. The time t_Th to reach the minimum of the hole is a constant for the GOE model, but grows exponentially with system size for the spin model <cit.>. This exceedingly long time is a consequence of the spatial locality of the initial state and of couplings of the realistic model.
        	
          *  The saturation region, for t> t_R. The relaxation time (or Heisenberg time) is the largest time scale of the system and is given by the inverse of the mean level spacing <cit.>. 
        
        
        Four distinct behaviors, at the same time scales identified for the survival probability, appear also for the spin autocorrelation function, as seen in Fig. <ref> (e). Similarly to what one finds for ⟨ P_S(t) ⟩, the correlation hole is evident for ⟨ I(t) ⟩ as well. Note, however, that we know analytically that the ratio between the saturation point and the minimum value of <P_S(t)> remains constant as the system size increases, but this analysis has not been done for ⟨ I(t) ⟩.
        
        For the quantities that are local in time, we observe two different behaviors before the Thouless time, one for t<1/Γ and another one for 1/Γ< t <t_Th, as evident in the plots for the inverse participation ratio [Fig. <ref> (c)], second-order Rényi entropy [Fig. <ref> (d)], and connected spin-spin correlation function [Fig. <ref> (f)], and also for the Shannon entropy (not shown). However, beyond t_Th, the effects of the correlation hole are minor for these quantities, and we basically see only fluctuations around their infinite-time averages. We can then say that the dynamics of ⟨IPR(t) ⟩, ⟨ S(t) ⟩, ⟨ Sh(t) ⟩, and ⟨ C(t) ⟩ saturate before the Thouless time. 
        
        
        
        
        
        
        
        
        
        § SURVIVAL PROBABILITY
        
        
        Despite being a central quantity in the analysis of systems out of equilibrium, not much is known about the self-averaging properties of the survival probability. Some of the existing works have focused on the two-level form factor, which corresponds to the long-time part of the survival probability. They include numerical studies about the spectral correlations of the hydrogen atom in a magnetic field <cit.> and theoretical arguments <cit.>, that both indicate the lack of self-averaging of the two-level form factor.
        
        Here, we provide an analytical expression for R_P_S^GOE(t) at all times, and show that the survival probability is not self-averaging at any time scale. We confirm numerically that this picture holds for physical chaotic models as well. We start the discussion below with estimates for R_P_S (t) for both models at short and long times, and then proceed with the presentation of the analytical result for R_P_S^GOE(t) and numerical results for both models.
        
        
        
         §.§ Short times
        
        
        For short times, t ≪ 1/Γ, one can expand the survival probability as
        
        P_S(t)=1-Γ_0^2 t^2+ O(t^4).
        
        From this expansion, one finds that the relative variance is given by
         R_P_S=σ^2_Γ^2t^4+ O(t^6),
        
        where σ^2_Γ^2=<Γ_0^4>-<Γ_0^2>^2. For the GOE and the realistic chaotic model, one has 
               (σ^2_Γ^2)^GOE=2 D , 
               (σ^2_Γ^2)^spin = J^4 L^2 (L-2)/64 (L-1)^2∼J^4 /64 L ,
        
        so the relative variance grows linearly with matrix and system size, respectively. The survival probability is non-self-averaging for both the GOE and the spin model.
        
        
        
        
         §.§ Long times
        
        
        Strong evidence for the absence of self-averaging of the survival probability at long times was already hinted at by studies of temporal fluctuations. In chaotic systems after saturation, the dispersion of the temporal fluctuations of P_S(t) is proportional to the value of the infinite-time average, P_S<cit.>, so the relative variance remains constant as L increases. The same result is obtained also for the relative variance of P_S(t) over the ensemble of realizations, because for chaotic systems and time intervals beyond t_R, temporal averages become comparable to ensemble averages.
        
        The survival probability in Eq. (<ref>) can also be written as
        
        P_S(t)=∑_α≠β |c_α^(0)|^2 |c_β^(0)|^2e^-i(E_α-E_β)t  + ∑_α | c_α^(0) |^4. 
        
        For times t>t_R, the ensembles averages cancel out the first term in Eq. (<ref>), so < P_S (t>t_R) > ∼⟨P_S⟩ =  ⟨∑_α | c_α^(0) |^4  ⟩. The eigenvectors of GOE random matrices are statistically equivalent to normalized Gaussian random vectors <cit.>, which gives ⟨P_S⟩∼ 3/D. For physical models, we also have ⟨P_S⟩∝ 1/D for eigenstates away from the borders of the spectrum. 
        
        To obtain the variance σ_P_S^2(t), in addition to ⟨ P_S(t)⟩ ^2, one needs the ensemble average of the square of the survival probability as well,
        <P_S^2(t)> =<∑_α,β,γ,δ e^-i(E_α-E_β+E_γ-E_δ)t|c_α^0|^2 |c_β^0|^2 |c_γ^0|^2 |c_δ^0|^2 >.
        
        At large times, only the terms with α=β and γ=δ, or α=δ and γ=β matter, because the phase factors average to zero.  Since 2 ⟨∑_α | c_α^(0) |^4   ∑_β | c_β^(0) |^4  ⟩∼ 2 ⟨∑_α | c_α^(0) |^4  ⟩⟨∑_β | c_β^(0) |^4  ⟩, one finds that σ^2_P_S∝⟨P_S⟩^2, and as a consequence
         R_P_S  (t>t_R)  = O(1).
        
        At long times, the relative variance is therefore independent of system size. 
        
        
        
        
         §.§ Analytical result for all times
        
        
        The derivation of the analytical expression for the mean value of the squared survival probability evolved under random matrices is lengthy, and is explained in Appendix <ref>. The formula reads
        <P_S^2(t)>    =    g_4(t)+ ⟨P_S⟩ g_3(t)+4 ⟨P_S⟩<P_S(t)-   ⟨P_S⟩> 
        
           +   ⟨P_S⟩^2 <P_S(2t) -   ⟨P_S⟩>+2  ⟨P_S⟩^2  .   
        
        In the above, the functions g_4(t) and g_3(t) are related to the Fourier transforms of the four- and three-point spectral correlation functions, respectively. The agreement between Eq. (<ref>) and the numerical results is perfect, as can be seen in the figure in Appendix <ref>. Combining Eq. (<ref>) and Eq. (<ref>), one obtains R_P_S^GOE(t) analytically. 
        
        
        
        In Fig. <ref> (a), we plot the numerical results for R_P_S^GOE(t) for various matrix sizes.  It is clear that the survival probability is not self-averaging at any time scale. The analytical expression for R_P_S^GOE(t) agrees very well with the numerics for t > 1/Γ, while for very short times, where the Fourier transform of the LDOS dominates the dynamics, finite size corrections are important. In this case, one computes R^GOE_P_S(t) using Eq. (<ref>). This is shown in the inset of Fig. <ref> (a), where the dots are the coefficients for lim_t→0 R_P_S^GOE(t)/t^4 extracted numerically and the dashed line shows the analytical behavior predicted from Eq. (<ref>) and Eq. (<ref>), lim_t→0 R_P_S^GOE(t)/t^4 ∝ D. 
        
        A qualitatively analogous behavior is shown for the spin model in Fig. <ref> (b). At short times, R^spin_P_S(t) grows with system size, and at large times, it is  size independent. The inset confirms the prediction by Eq. (<ref>) and Eq. (<ref>), indicating that  lim_t→0 R_P_S^spin(t)/t^4 is proportional to L. This coefficient grows slower for the spin model than for the GOE model, because of the sparseness of the realistic Hamiltonian. 
        
        At intermediate times, where in Fig. <ref>, < P_S(t) > shows oscillations decaying as a power law, R_P_S(t) also oscillates for both models, as observed in Fig. <ref>. In the case of the physical model, the relative variance in this region reaches values above 1. Since the power-law decay of < P_S(t) > is caused by the bounds of the spectrum of finite systems <cit.>, at such intermediate times,  the state |Ψ(t)> acquires weight on eigenstates closer to the edges of the spectrum, which are not described by random matrix theory. The values of R^spin_P_S(t) above 1 could be a manifestation of correlations between the components of these states. 
        
        Beyond the region of the power-law decay, the relative variance of the survival probability behaves similarly for any time, R_P_S(t>t_Th) ∼ 1. This suggests that the correlation hole, which is clearly manifested in the mean value of the survival probability in Fig. <ref>, does not affect the self-averaging properties of this quantity in any way different from what one sees for t>t_R. 
        
        
        
        
        
        § INVERSE PARTICIPATION RATIO AND RÉNYI ENTROPIES
        
        
        The inverse participation ratio and Rényi entropies are defined in Eq. (<ref>), Eq. (<ref>), and Eq. (<ref>). While their self-averaging properties are similar for t>1/Γ, they differ for short times.
        
        
        
         §.§ Inverse participation ratio
        
        There is a fundamental difference between P_S(t) and IPR(t): the survival probability measures the distance of the evolved state |Ψ(t)> from the initial state (indicated as n=0), while the inverse participation ratio measures the distance of |Ψ(t)> from any unperturbed many-body state |n>. The inverse participation ratio can be seen as a generalization of the squared survival probability, where in addition to the term with n=0, which gives P_S^2(t) itself, it contains also all other terms with n≠ 0. This parallel allows us to intuitively understand the behavior of the relative variance of IPR. At short times, |Ψ(t)> is still very close to |Ψ(0)>, so the term with n=0 dominates the evolution and IPR(t) behaves analogously to P_S^2(t)<cit.>. This means that self-averaging is absent at short times. At large times, on the other hand, the average over all unperturbed many-body states drastically reduces the fluctuations, and the inverse participation ratio becomes very strongly self-averaging, in the sense that R_IPR∝ D^-1 for both the GOE and the spin model. 
        
        For times t ≪ 1/Γ, IPR(t) can be expanded as
        IPR(t)=1-2Γ_0^2t^2+ O(t^4),
        
        which is exactly the same expression one has for P_S^2(t) at leading order in t. This means that 
         R_IPR(t) ∝ R_P_S(t) ∝σ_Γ^2^2t^4+ O(t^6)
        
        at short times, for any model, so the same kind of non-self-averaging behavior found in Sec. <ref> emerges here also.
        
        At large times, we study the ensemble average of 
        IPR(t)   =   ∑_n∑_α,β,γ,δe^-i(E_α-E_β+E_γ-E_δ)t
           ×    c_α^0 c_α^n* c_β^n c_β^0* c_γ^0 c_γ^n* c_δ^n c_δ^0*
        and  of IPR^2(t) using the same arguments employed in Sec. <ref>, namely that the eigenstates of GOE matrices imply that < c_α> ∼ 0, < |c_α|^2 > ∼ 1/D, and < ∑_α ,β e^-i(E_α-E_β)t> ∼ 0 unless α = β. One finds that
         R_IPR^GOE (t > t_R) ∝1/D.
        
        Thus, unlike the survival probability, IPR(t) is actually self-averaging at large times.
        
        
        
        Not all eigenstates of realistic chaotic models are close to normalized Gaussian random vectors, but they are the majority, so we should expect a similar behavior for the realistic spin model as well.
        
        In Fig. <ref>, we plot the numerical data for R_IPR^GOE(t) in panel (a), and for R_IPR^spin(t) in panel (d). In both cases, the lack of self-averaging at short times and the very strong self-averaging at large times are clearly visible. This is shown quantitatively in the other panels. In Fig. <ref> (b) and (e), we plot the short time coefficient lim_t→0 R_IPR(t)/t^4 for the GOE and the spin model, respectively. In Fig. <ref> (c) and (f), we show the value of R_IPR(t) for a long time t ≥ t_R.  In all four cases, the numerical values (circles) agree very well with our analytical estimates (dashed lines) in Eq. (<ref>) and Eq. (<ref>), except for a small noticeable deviation for the spin model with L=18 in Fig. <ref> (f).
        
        At intermediate times, in the region of the power-law decay of P_S(t) and similarly to the behavior of R_P_S(t), the relative variance of the inverse participation ratio oscillates and reaches its largest values, as seen in Fig. <ref> (a) and Fig. <ref> (d). Beyond this region, the curves are pretty much flat and R_IPR is  unaware of the correlation hole.
        
        
        
         §.§ Rényi entropies
        
        
        Similarly to the behavior of the inverse participation ratio, the second-order Rényi entropy is super self-averaging at long times and non-self-averaging at the time scales of the power-law decay of P_S(t) (see figures for R_S(t) in Appendix <ref>). The two quantities differ, however, at short times. This happens, because for t ≪ 1/Γ, 
        
        S(t)=2Γ_0^2t^2+ O(t^4),
        
        so S(t)→0 for t→0, while the inverse participation ratio goes to 1 for t→0 [see Eq. (<ref>)].  Contrary to R_IPR(t), the time dependence of R_S(t) cancels out at the lowest order in t,
         R_S(t)=σ^2_Γ^2/<Γ_0^2>^2+ O(t^2).
        
        Using Eq. (<ref>), Eq. (<ref>), and Eq. (<ref>), one finds that the second-order Rényi entropy is self-averaging at short times for both random matrices and physical models,
         R_S^GOE(t)=2/D+ O(t^2),    R_S^spin(t)=1/L+ O(t^2).
        
        
        The difference in the behavior of IPR(t) and -ln [IPR(t)] is somewhat reminiscent to what happens in  the Anderson model, where the transmission amplitude, which scales multiplicatively with the system size, is not self-averaging, while its logarithm, which scales additively, is self-averaging <cit.>. The fact that  S(t) is self-averaging at short times makes it more appealing for experiments than IPR(t). 
        
        From the point of view of self-averaging properties, whether one uses the second-order Rényi entropy or the Shannon entropy is indifferent. We verified that both R_S(t) and R_Sh(t) exhibit equivalent behaviors.
        
        
        
        § EXPERIMENTAL LOCAL QUANTITIES
        
        
        The two experimental quantities considered, the spin autocorrelation function and the connected spin-spin correlation function, are local in space. Since for random matrices, the notion of locality is meaningless, we study these quantities for the spin model only. 
        
        The main difference between the two observables is that the connected spin-spin correlation function is also local in time, while the spin autocorrelation function is not. In this sense, the spin autocorrelation function is the spatially local counterpart of the survival probability, since both are measured with respect to the state of the system at t=0, and the connected correlation function is the spatially local counterpart of the inverse participation ratio, both involving only the state at t and both dealing also with averages over all unperturbed many-body states.
        
        The differences between these quantities are reflected in their self-averaging properties. At short times, spatial locality ensures that both observables are self-averaging. Their dynamics involve only a finite number of spins, and the spatial averages ensure that the relative variances are reduced as the system size increases. At long times, on the other hand, the spin autocorrelation function, just as the survival probability, is not self-averaging, while the connected correlation function, similarly to the inverse participation ratio, is.
        
        
        
         §.§ Spin autocorrelation function
        
        
        At short times, t ≪ 1/Γ, one can expand the spin autocorrelation function in the following way,
        
        I(t) =1-Γ_0^2t^2+4t^2/L∑_k=1^LS_k^00∑_n≠0| ⟨ n|H|Ψ(0)⟩|^2 S_k^nn+ O(t^4),
        
        where S_k^nn=<n|S_k^z|n>. The third term on the right-hand side of Eq. (<ref>) is zero, unless |n⟩ is directly coupled with |Ψ(0)⟩. These states |n⟩ differ from |Ψ(0)⟩ by two neighboring sites only. Therefore, using the definition of Γ_0^2 in Eq. (<ref>), one finds that 
        
            ∑_k=1^LS_k^00∑_n≠0| ⟨ n|H|Ψ(0)⟩|^2S_k^nn = L-4/4Γ_0^2 ,
        
        which gives 
        
        I(t)=1-4  Γ_0^2  t^2/L+ O(t^4).
        
        From this expansion, we obtain the relative variance, 
         R_I(t) = 16 σ_Γ^2^2 t^4/L^2+ O(t^6) ∝J^4 t^4/L,
        
        where we used that σ_Γ^2^2 ∝ L [see Eq. (<ref>)]. 
        
        The estimates above provide us with two important results about the spin autocorrelation function at short times: (i) Its mean value <I(t)> is independent of system size, since according to Eq. (<ref>), < Γ_0^2 > ∝ L; and (ii) the quantity is strongly self-averaging, since the relative variance decays as 1/L. Both these features are a consequence of the local structure of the observables and of the Hamiltonian, and not of any peculiar property of I(t). As a consequence, we can claim that the local quantities studied here, or any other involving only sums of local operators and evolving under local Hamiltonians, are self-averaging for short times.
        
        
        
        The relative variance R_I(t) is plotted  in Fig. <ref> (a). At short times, we observe the expected power-law behavior ∝ t^4, with a coefficient lim_t→0 R_I(t)/t^4∝ 1/L, as shown in Fig. <ref> (b). The self-averaging behavior persists up to times currently reachable experimentally, so self-averaging can be safely assumed in real experiments. 
        
        At large times, however, an inversion happens, and R_I(t) starts growing with system size, as seen in Fig. <ref> (a) and Fig. <ref> (c). This happens at times of the order of the Thouless time, when ⟨ I(t) ⟩ enters the correlation hole [see Fig. <ref> (e)], and the dynamics crossover from a model-dependent regime at short times, to a universal regime at long times <cit.>.  After this point, self-averaging is lost, and in an even stronger sense than for the survival probability, for which R_P_S(t) for t > t_Th is system size independent. At these very long time scales, the spin autocorrelation function presents positive and negative values very close to zero, which results in large fluctuations.
        
        
        
         §.§ Connected spin-spin correlation function
        
        
        Using a short-time expansion for the connected spin-spin correlation function, one finds that, 
        
        C(t)    =   4 t^2/L∑_k=1^L  ∑_n ≠ 0( S_k^00 S_k+1^00  - 2 S_k^00 S_k+1^nn + S_k^nn S_k+1^nn) ×
              | ⟨ n|H|Ψ(0)⟩|^2 + O(t^4)  = -2 Γ_0^2 t^2/L + O(t^4).
        
        This implies that, as a consequence of space locality and since ⟨Γ_0^2 ⟩∝ L, the mean value of this quantity is system size independent, just like the spin autocorrelation function, and it is also self-averaging, 
         R_C(t) = σ_Γ^2^2 /⟨Γ_0^2 ⟩^2 + O(t^6)
        ∼1/L.
        
        
        The main difference between the connected spin-spin correlation function and the spin autocorrelation at short times is that C(0)=0, while I(0)=1. Therefore, the relative variance R_C(t) tends to a finite value as t→0, while R_I(t)→0. In Fig. <ref>  (d), we plot R_C(t) for the spin model. At short times, the relative variance indeed saturates to a finite value. Figure <ref> (e) confirms that this value decays as 1/L.
        
        At large times, R_C (t) saturates to an asymptotic value that decreases exponentially with L, as shown in Fig. <ref> (f). This is similar to what happens for the inverse participation ratio, and in contrast to the behavior of the spin autocorrelation function. We conjecture that this very strong self-averaging behavior at long times is associated to the fact that no memory of the initial state is encoded in the connected spin-spin correlation function.
        
        
        
        
        
        
        § IMPACT OF LACK OF SELF-AVERAGING
        
        
        Lack of self-averaging, and therefore not being able to decrease the number of random realizations (or samples) as the system size L increases, can make numerical and experimental analysis quite challenging. It is also very dangerous to assume self-averaging, as one often does when studying chaotic systems, and then decrease the number of realizations as L increases for a quantity that is in fact non-self-averaging. The purpose of this section is to provide an example of this case.
        
        
        
        In Fig. <ref>, we show the survival probability for the chaotic spin model (<ref>) averaged over a number of random realizations. This number decreases as L increases in Fig. <ref> (a), but remains constant as L increases in Fig. <ref> (b). Main features of the evolution of the survival probability are missed in Fig. <ref> (a). There is a hint for the correlation hole for L=10, which is entirely lost for L=18. There is an indication of a power-law decay before the correlation hole for L=14, but this is also lost for L=18. As a result, motivated by the curve for the largest system size, one could assume that what one sees for the smaller sizes are just finite size effects, and that after the Gaussian decay, the dynamics simply saturates. Worse still, since Γ∝√(L) and the Gaussian decay e^-Γ^2 t^2 gets faster as L increases, one could be led to conclude that the relaxation time decreases as the chain grows, while in reality the relaxation time for ⟨ P_S (t) ⟩ increases exponentially with L (see Sec. <ref> and Ref. <cit.>).
        
        In Fig. <ref> (b), both the power-law decay and the correlation hole are visible for L=14 and L=18 confirming that these features are indeed present. The results would be even more convincing if instead of keeping the number of realizations fixed, we would increase it.
        
        A critical ongoing debate in the current literature about nonequilibrium quantum dynamics is whether existing conclusions may be incomplete or even incorrect due to finite size effects. Our results show that self-averaging is another important element that needs to be taken into account. For a quantity that is non-self-averaging, one can reach wrong conclusions if one assumes self-averaging and decreases the number of realizations.
         
        
        
        
        
        
        
        § CONCLUSIONS
        
        
        This work analyzes the self-averaging behavior of many-body quantum systems out of equilibrium. The focus is on the chaotic regime, where self-averaging is often taken for granted. By examining different non-local and local quantities in space, we bring forward a rich variety of behaviors and deduce that self-averaging is not an intrinsic consequence of quantum chaos, but depends strongly on the quantity and on the time scale. 
        
        On the bright side, the local quantities studied here and measured in experiments with cold atoms and ion traps, namely the spin autocorrelation function and the connected spin-spin correlation function, are self-averaging for the times that are currently experimentally reachable. The same arguments that we employ for these quantities can be extended to any observable comprising only sums of spatially local operators and evolving under local Hamiltonians, so they should also be self-averaging at short times, a result that is reassuring for experimentalists. 
        
        Numerical studies, on the other hand, where long times and an array of quantities are available, should be cautious. Autocorrelation functions, such as the survival probability and the spin autocorrelation function, are not self-averaging at long times, so one needs large statistics even when pushing towards very large system sizes. In fact, as we showed analytically for full random matrices, the survival probability is not self-averaging at any time scale. Extra care should therefore be taken when dealing with this quantity, which has a central role in studies of nonequilibrium quantum dynamics.
        
        The time evolution of the fluctuations of observables is still uncharted territory. There are multiple interesting directions that the study initiated here could take, from the analysis of non-chaotic models, such as those approaching many-body localization, to time-dependent Hamiltonians and open systems. Our results lay the groundwork for such future analysis.
        
        
        
        
        M.S. and L.F.S. were supported by the NSF Grant No. DMR-1603418. E.J.T.-H. acknowledges funding from VIEP-BUAP (Grant Nos. MEBJ-EXC19-G, LUAGEXC19-G), Mexico. He is also grateful to LNS-BUAP for allowing use of their supercomputing facility. F.P.B. thanks the Consejería de Conocimiento, Investigación y Universidad, Junta de Andalucía and European Regional Development Fund (ERDF), ref. SOMM17/6105/UGR. Additional computer resources supporting this work were provided by the Universidad de Huelva CEAFMC High Performance Computer  located in the Campus Universitario el Carmen and funded by FEDER/MINECO project UNHU-15CE-2848. LFS is supported by the NSF Grant No. DMR-1936006. Part of this work was performed at the Aspen Center for Physics, which is supported by the NSF grant PHY-1607611.
        
        
        
        § ANALYTICAL EXPRESSION FOR <P_S^2(T)>
        
        
        We provide here the analytical derivation for R_P_S^GOE(t). The analytical expression for <P_S(t)> is given in Eq. (<ref>) and was already obtained in <cit.>. It remains to show how we obtain Eq. (<ref>) for <P_S^2(t)>, which is a much more involved derivation.
        
        Before presenting the steps of the derivation, we compare in Fig. <ref>  the expression from Eq. (<ref>) (dashed lines) with the numerical results (solid lines) for GOE matrices of different sizes. The agreement is indeed perfect.
        
        
        
        To obtain the relative variance of the survival probability, we need to compute
        <P_S^2(t)> = Ξ_1 + Ξ_2 + Ξ_3 + Ξ_4,
        
        where
        Ξ_1    =   <∑_α≠γ≠β≠δ
        e^-i(E_α-E_β+E_γ-E_δ)t|c^(0)_α|^2|c^(0)_β|^2|c^(0)_γ|^2|c^(0)_δ|^2> ,
        
        Ξ_2    =   <∑_α≠β≠γe^-i(2E_α-E_β-E_γ)t|c^(0)_α|^2|c^(0)_β|^4|c^(0)_γ|^2> 
           +   <∑_α≠β≠γe^-i(E_α-2E_β+E_γ)t|c^(0)_α|^4|c^(0)_β|^2|c^(0)_γ|^2> ,
        Ξ_3    =   <∑_α≠βe^-2i(E_α-E_β)t|c^(0)_α|^4|c^(0)_β|^4> 
           +    4<∑_α≠βe^-i(E_α-E_β)t|c^(0)_α|^2|c^(0)_β|^2∑_γ|c^(0)_γ|^4>  , 
        Ξ_4    =    2<∑_α,β|c^(0)_α|^4|c^(0)_β|^4> .  
        Equation (<ref>) is obtained by splitting the sum in Eq. (<ref>) into all possible combination of equal indexes α,β,γ,δ. For example, Ξ_1 contains the terms with all indexes different, Ξ_2 the terms with either α=γ or β=δ, and so on. We now compute each one of these terms, starting from Ξ_4 and moving upwards.
        
        
        
         §.§ Term Ξ_4
        
        
        The fourth term in Eq. (<ref>) can be computed straightforwardly, applying the results of Refs. <cit.>. Indeed, up to subleading corrections in D^-1, one finds
        
        2<∑_α,β|c^(0)_α|^4|c^(0)_β|^4>=2<P_S>^2+O(D^-3).
        
        Note that, at large times, this is the only term that does not vanish. This is because it is the only term which does not contain fluctuating phases in time. We can therefore compute the asymptotic value of the relative variance already,
        lim_t→∞ℛ_P_S^GOE(t)   =   2<P_S>^2-<P_S>^2/<P_S>^2+O(D^-1)
           =   1+O(D^-1).
        
        This means that, as explained in Sec. <ref>, the relative variance is system size independent at large times, and self-averaging is not present.
        
        
        
         §.§ Term Ξ_3
        
         Let us first compute the following term of Ξ_3,
         
         4<∑_α≠βe^-i(E_α-E_β)t|c^(0)_α|^2|c^(0)_β|^2∑_γ|c^(0)_γ|^4>.
         
         In order to do this, let us define the time dependent part of the survival probability
         P_S(t)=P_S(t)-P_S=∑_α≠βe^-i(E_α-E_β)t|c^(0)_α|^2|c^(0)_β|^2,
         
         which tends to zero as t→∞. From Eq. (<ref>), one gets for its average value
         <P_S(t)>=1-<P_S>/D-1[D b_1^2( Γ t)-b_2(Γ t/2D)].
        
        As a consequence, one finds, up to subleading corrections in D^-1, that
              4 <∑_α≠βe^-i(E_α-E_β)t|c^(0)_α|^2|c^(0)_β|^2∑_γ|c^(0)_γ|^4>
           =   4<P_S><P_S(t)>+O(D^-2).
        
        
        We now compute the other term of Ξ_3,
        <∑_α≠βe^-2i(E_α-E_β)t|c^(0)_α|^4|c^(0)_β|^4>.
        
        For this, one needs to recall that for GOE matrices, eigenvalues and eigenvectors are statistically independent <cit.>. This means that the averages over the components and the eigenenergies factorize. This fact will be used multiple times in this derivation. For this particular term, it implies that
              <∑_α≠βe^-2i(E_α-E_β)t|c^(0)_α|^4|c^(0)_β|^4> 
           =   <∑_α≠β|c^(0)_α|^4|c^(0)_β|^4><e^-2i(E_α-E_β)t>
           =   <P_S(2t)><P_S>^2+O(D^-3)/1-<P_S>
           ∼   <P_S(2t)>(<P_S>^2+O(D^-2)).
        
         §.§ Term Ξ_2
        
        
        We now need to compute Ξ_2,
        <∑_α≠β≠γe^-i(2E_α-E_β-E_γ)t|c^(0)_α|^2|c^(0)_β|^4|c^(0)_γ|^2>   +   
        <∑_α≠β≠γe^-i(E_α-2E_β+E_γ)t|c^(0)_α|^4|c^(0)_β|^2|c^(0)_γ|^2>.
        
        Using again the statistical independence of eigenvalues and eigenvectors, we get for the components of the initial state,
        <∑_α≠β≠γ|c^(0)_α|^2|c^(0)_β|^4|c^(0)_γ|^2>=<P_S>+O(D^-2).
        
        
        It remains to obtain the following function,
        
        g_3(t)=<e^-i(2E_α-E_β-E_γ)t+e^-i(E_α-2E_β+E_γ)t>.
        
        Following Refs. <cit.>, we write this average as
        
        g_3(t)   =   (D-3)!/D!∫ dE e^-iEt∫ dE_α dE_β dE_γ R_3(E_α,E_β,E_γ)
           ×   [δ(E-(2E_α-E_β-E_γ))+δ(E-(E_α-2E_β+E_γ))],
        where R_3(E_α,E_β,E_γ) is the three-point spectral correlation function <cit.>. It can be written as
        
        R_3(E_α,E_β,E_γ)   =   R_1(E_α)R_1(E_β)R_1(E_γ)-R_1(E_α)T_2(E_β,E_γ)
           -   R_1(E_β)T_2(E_α,E_γ)-R_1(E_γ)T_2(E_α,E_β)
           +   T_3(E_α,E_β,E_γ).
        
        In the above,
        
        R_1(E)=1/π√(2D-E^2)
        is the density of states, while T_2 and T_3 are the two- and three-point cluster functions, respectively. T_3 represents the connected part of the three-point correlation function, while all other terms of R_3 represent all possible disconnected contributions.
        
        We first compute the term which depends on the density of states only,
        (D-3)!/D!∫ dE e^-iEt   ×   ∫ dE_α dE_β dE_γ R_1(E_α)R_1(E_β)R_1(E_γ)[δ(E-(2E_α-E_β-E_γ))+δ(E-(E_α-2E_β+E_γ))]
           =   (D-3)!/D!∫ dE_α e^-2iE_α tR_1(E_α)∫ dE_β e^-iE_β tR_1(E_β)∫ dE_γ e^-iE_γ tR_1(E_γ)+c.c.
           =   2(D-3)!/D!(2D)^3b_1(2Γ t)b_1^2(Γ t).
        
        
        We now consider the terms containing both R_1 and T_2. As an example, let us compute
        (D-3)!/D!∫ dE e^-iEt   ×   ∫ dE_α dE_β dE_γ R_1(E_α)T_2(E_β,E_γ)[δ(E-(2E_α-E_β-E_γ))+δ(E-(E_α-2E_β+E_γ))]
           =   (D-3)!/D!∫ dE_α R_1(E_α)e^-2iE_α t∫ dE_β dE_γ T_2(E_β,E_γ)e^-i(-E_β-E_γ)t
           +   (D-3)!/D!∫ dE_α R_1(E_α)e^-iE_α t∫ dE_β dE_γ T_2(E_β,E_γ)e^-i(-2E_β+E_γ)t.
        From Ref. <cit.>, we know that, for any n≥2, integrals of the form
        ∫ dE_1dE_2⋯ dE_nT_n(E_1,E_2,_⋯ E_n)e^-i∑_j=1^n k_j E_j
        are non-vanishing if and only if ∑_j=1^n k_j=0. Since this is not the case in the integrals of Eq. (<ref>), all terms in g_3 containing both R_1 and T_2 vanish.
        
        Finally, we consider the term containing the T_3 connected correlation function,
        (D-3)!/D!∫ dE e^-iEt   ×   ∫ dE_α dE_β dE_γ T_3(E_α,E_β,E_γ)[δ(E-(2E_α-E_β-E_γ))+δ(E-(E_α-2E_β+E_γ))]
           =   (D-3)!/D!∫ dE_α dE_β dE_γ T_3(E_α,E_β,E_γ)[e^-i(2E_α-E_β-E_γ)t+e^-i(E_α-2E_β+E_γ)t].
        
        We perform the change of variables ξ_α=E_α/R_1(0), and similarly for E_β and E_γ. This corresponds to rescale the energies by their mean level spacing. Calling Y_3(ξ_α,ξ_β,ξ_γ)=R_1^3(0)T_3(E_α,E_β,E_γ), this integral becomes
        (D-3)!/D!∫ dξ_α dξ_β dξ_γ Y_3(ξ_α,ξ_β,ξ_γ)   ×   
        [e^-i(2ξ_α-ξ_β-ξ_γ)R_1(0)t+e^i(ξ_α-2ξ_β+ξ_γ)R_1(0)t] ,
        
        which can be computed using the following formula, found in Ref. <cit.>,
              ∫ dξ_1dξ_2⋯ dξ_nY_n(ξ_1,ξ_2,_⋯,ξ_n)e^-i∑_j=1^nξ_jτ_j
           =   δ(τ_1+τ_2+_⋯+τ_n)
        ∫_-∞^∞ dτ[  ∑_Pf(τ)f(τ+τ_P(1))_⋯.
        
               f(τ+τ_P(1)+ _⋯+  . τ_P(n-1))   ]_(0).
        In this formula, P labels all permutations of the indexes 1,2,3, … n, and P(i) is the permuted counterpart of i, according to P. f(τ) is a matrix valued function, which reads
        
        f(τ)=[     f_2(τ)   τ f_2(τ); f_2(τ)-1/τ     f_2(τ) ],
          
        f_2(τ)=
        1    |τ|<1/2 
        
        0    |τ|≥1/2
        .
        
        With the notation [·]_(0), we mean the following: Taken a generic 2×2 matrix A, we call
        [A]_(0)=A_11+A_22/2,
        
        where A_11 and A_22 are the diagonal entries of A. In particular, for Eq. (<ref>), we have τ_1=2R_1(0), τ_2=τ_3=-R_1(0) for the first term, and τ_1=τ_3=R_1(0), τ_2=-2R_1(0) for the second. Thus, we obtain for Eq. (<ref>), 
              (D-3)!/D!∫ dξ_α dξ_β dξ_γ Y_3(ξ_α,ξ_β,ξ_γ) 
              ×[e^-i(ξ_α-ξ_β-ξ_γ)R_1(0)t+e^i(ξ_α-2ξ_β+ξ_γ)R_1(0)t] 
           =   (D-3)!/D!Db_3(Γ t/2D),
        
        with
        
        b_3(τ)=
        4[1-4τ+3τlog(2τ+1)]     0<τ≤1/2
        
        8[τ-1+3/2τlog(4-τ/2+τ)]    τ>1/2.
        
        
        Putting all these terms together, we get
        
        g_3(t)=(D-3)!/D![16D^3b_1(2Γ t)b_1^2(Γ t)-Db_3(Γ t/2D)].
        
        
         §.§ Term Ξ_1
        
        
        Finally, we compute the term Ξ_1,
        <∑_α≠γ≠β≠δe^-i(E_α-E_β+E_γ-E_δ)t|c^(0)_α|^2|c^(0)_β|^2|c^(0)_γ|^2|c^(0)_δ|^2>.
        
        Using again the independence of eigenvalues and eigenvectors, we factorize the average, and compute the part depending on the eigenvectors first,
        <∑_α≠γ≠β≠δ|c^(0)_α|^2|c^(0)_β|^2|c^(0)_γ|^2|c^(0)_δ|^2>   =   
        
        1-6<P_S>+O(D^-2)   ∼   1.
        
        
        We now compute the average over the energy levels,
        
        g_4(t)   =   <e^-i(E_α-E_β+E_γ-E_δ)>
           =   (D-4)!/D!∫ dE e^-iEt
           ×   ∫ dE_α dE_β dE_γ dE_δ R_4(E_α,E_β,E_γ,E_δ)
           ×   δ(E-(E_α-E_β+E_γ-E_δ)).
        R_4(E_α,E_β,E_γ,E_δ) is the four-point correlation function, which can be written as <cit.>
        R_4(E_α,E_β,E_γ,E_δ)   =   R_1(E_α)R_1(E_β)R_1(E_γ)R_1(E_δ)-R_1(E_α)R_1(E_β)T_2(E_γ,E_δ)-R_1(E_α)R_1(E_γ)T_2(E_β,E_δ)
           -   R_1(E_α)R_1(E_δ)T_2(E_β,E_γ)-R_1(E_β)R_1(E_γ)T_2(E_α,E_δ)-R_1(E_β)R_1(E_δ)T_2(E_α,E_γ)
           -   R_1(E_γ)R_1(E_δ)T_2(E_α,E_β)+T_2(E_α,E_β)T_2(E_γ,E_δ)+T_2(E_α,E_γ)T_2(E_β,E_δ)
           +   T_2(E_1,E_4)T_2(E_2,E_3)+R_1(E_1)T_3(E_2,E_3,E_4)+R_1(E_2)T_3(E_1,E_3,E_4)
           +   R_1(E_3)T_3(E_1,E_2,E_4)+R_1(E_4)T_3(E_1,E_2,E_3)-T_4(E_α,E_β,E_γ,E_δ).
        
        Once again, this amounts to consider the connected correlation function T_4(E_α,E_β,E_γ,E_δ) and all disconnected components. We now study the resulting integrals, one by one.
        
        We first consider the term containing R_1 only:
        (D-4)!/D!   ×   ∫ dE_α dE_β dE_γ dE_δ R_1(E_α)R_1(E_β)R_1(E_γ)R_1(E_δ)e^-i(E_α-E_β+E_γ-E_δ)t
           =   (D-4)!/D!|∫ dE R_1(E)e^-iEt|^4=(D-4)!/D!(2D)^4b_1^4(Γ t).
        
        
        Next, we consider terms containing both R_1 and T_2, such as
        
            (D-4)!/D!∫ dE_α dE_β dE_γ dE_δ R_1(E_α) R_1(E_β) T_2(E_γ,E_δ) e^-i(E_α-E_β+E_γ-E_δ)t.
        
        This term can be rewritten as
        (D-4)!/D!   ×   |∫ dE R_1(E)e^-iEt|^2
           ×   ∫ dE_γ dE_δ T_2(E_γ,E_δ)e^-i(E_α-E_β)t
           =   (D-4)!/D!(2D)^2b_1^2(Γ t)Db_2(Γ t/2D).
        Note that this integral is non zero only if the two energies in the T_2 function come with opposite signs in the corresponding exponential. For this reason, this term comes with a combinatorial factor equal to 4.
        
        Next are the terms containing T_2 only,
        (D-4)!/D!∫ dE_α dE_β dE_γ dE_δ T_2(E_α,E_β)T_2(E_γ,E_δ)e^-i(E_α-E_β+E_γ-E_δ)t   =   
        (D-4)!/D!(∫ dE_α dE_β T_2(E_α,E_β)e^-i(E_α-E_β)t)^2   =   (D-4)!/D!D^2b_2^2(Γ t).
        
        This term comes with a combinatorial factor equal to 2.
        
        Then come the terms with R_1 and T_3, such as
        (D-4)!/D!   ×   ∫ dE_α dE_β dE_γ dE_δ
           ×   R_1(E_α)T_3(E_β,E_γ,E_δ)e^-i(E_α-E_β+E_γ-E_δ)t.
        However, as explained for the g_3(t) function, all the integrals of this form vanish. 
        
        Finally, comes the term containing the T_4 function:
        (D-4)!/D!   ×   ∫ dE_α dE_β dE_γ dE_δ T_4(E_α,E_β,E_γ,E_δ)
           ×   e^-i(E_α-E_β+E_γ-E_δ)t=(D-4)!/D!Db_4(Γ t/2D).
        We compute the function b_4 with a procedure analogous to the one used for b_3. The result is
        
        b_4(τ)=4τ^3-20τ^2-4τ+3/2τ+1+6τlog(2τ+1)    τ≤1/2
        
        22τ^3+8τ^2-2τ-3/2τ+1+3τ[2log(2τ+1)-3log(4τ-1)+log2]    1/2<τ<1 
        
        26τ^2-1/4τ^2-1+3τlog(2τ+1/2τ-1)
        .
        
        So the g_4 function reads
        
        g_4(t)=(D-4)!/D![16D^4b_1^4(Γ t)-16D^3b_1^2(Γ t)b_2(Γ t/2D)+2D^2b_2^2(Γ t/2D)-Db_4(Γ t/2D)].
        
        Combining all the terms together, one recovers Eq. (<ref>) from the main text.
        
        
        
        § PLOTS FOR THE RELATIVE VARIANCE OF THE SECOND-ORDER RÉNYI ENTROPY
        
        
        For completeness, we present in Fig. <ref> (a) the relative variance of the second-order Rényi entropy as a function of time. The behavior for the relative variance of the Shannon entropy is very similar (not shown). 
        
        Contrary to the inverse participation ratio, the second-order Rényi entropy is self-averaging at short times. Indeed, as seen in Fig. <ref> (b), the relative variance follows very well the prediction in Eq. (<ref>) that  R_S^spin(t) ∝ 1/L, as we find also for the spatially local quantities discussed in Sec. <ref>. 
        
        At long times, the second-order Rényi entropy behaves similarly to the inverse participation ratio and the connected spin-spin correlation function, being super self-averaging, since R_S^spin(t) ∝ 1/D, as shown in Fig. <ref> (c).  The behavior of R_S^spin(t) for short- and long-times is therefore similar to that of R_C^spin(t).
        
        
        
        Between the two extremes of short and long times, in the region of  the power-law behavior of ⟨ P_S (t) ⟩, the second-order Rényi entropy is not self-averaging. This means that the relative variance of S (same for Sh) exhibits two crossing points as time increases from zero, a feature that contrasts with those of the other quantities studied in this paper, P_S, IPR, I, and C.
        
        We fitted R_S^spin (t>t_R) in Fig. <ref> (c) and R_IPR^spin (t>t_R) in Fig. <ref> (f) with 1/D for making an analogy with the results for the GOE model, but the exponent ν in the best fit 1/D^ν is 1± 0.2. In Fig. <ref> (f) for the connected spin-spin correlation function, where we do not make a comparison with the GOE model, we actually show the best fit. Independently of the exact value of the exponent ν, the relative variance at long times for any of these three quantities decreases exponentially fast with L.
        
        
        
        
        
        
        
        
        
        
        
        84
        fxundefined [1]
         ifx#1
        fnum [1]
          #1firstoftwo
         secondoftwo
         
        fx [1]
          #1firstoftwo
         secondoftwo
         
        
        
        
        
        
        noop [0]secondoftwo
        
        ref[1]@startlink#1@href
        href[1]#1@endlink
        anitize@url [0] `
        12 `$12
          `&12 `#12 `1̂2 `_12 `%12
        startlink[1]
        endlink[0]
        
        rl [1]href #1
        
        
        
        
        
        
        
        
        
        
        
        
        @bib@innerbibempty
        
         [Lifshitz et al.(1988)Lifshitz, Gredeskul, and Pastur]LifshitzBook
          authorauthorI. M. Lifshitz, authorS. A. Gredeskul,  and authorL. A. Pastur, @noop titleIntroduction to the
          Theory of Disordered Systems (publisherWiley, addressNew York, year1988)NoStop
         [Mehta(1991)]MehtaBook
          authorauthorM. L. Mehta, @noop titleRandom Matrices  (publisherAcademic Press, addressBoston,year1991)NoStop
         [Wiseman and Domany(1995)]Wiseman1995
          authorauthorS. Wiseman and authorE. Domany, titletitleLack of
          self-averaging in critical disordered systems, 
          10.1103/PhysRevE.52.3469journaljournalPhys.
          Rev. E volume52, pages3469
          (year1995)NoStop
         [Aharony and Harris(1996)]Aharony1996
          authorauthorA. Aharony and authorA. B. Harris, titletitleAbsence of
          self-averaging and universal fluctuations in random systems near critical
          points,  10.1103/PhysRevLett.77.3700journaljournalPhys. Rev. Lett. volume77, pages3700 (year1996)NoStop
         [Wiseman and Domany(1998)]Wiseman1998
          authorauthorS. Wiseman and authorE. Domany, titletitlefinite-size
          scaling and lack of self-averaging in critical disordered systems,  10.1103/PhysRevLett.81.22journaljournalPhys. Rev. Lett. volume81, pages22 (year1998)NoStop
         [Castellani and Cavagna(2005)]Castellani2005
          authorauthorT. Castellani and authorA. Cavagna, titletitleSpin-glass
          theory for pedestrians,  10.1088/1742-5468/2005/05/p05012journaljournalJ. Stat. Mech. Th. Exp.volume2005, pagesP05012 (year2005)NoStop
         [Malakis and Fytas(2006)]Malakis2006
          authorauthorA. Malakis and authorN. G. Fytas, titletitleLack of
          self-averaging of the specific heat in the three-dimensional random-field
          Ising model,  10.1103/PhysRevE.73.016109journaljournalPhys. Rev. E volume73, pages016109 (year2006)NoStop
         [Roy and Bhattacharjee(2006)]Roy2006
          authorauthorS. Roy and authorS. M. Bhattacharjee, titletitleIs
          small-world network disordered? 
          https://doi.org/10.1016/j.physleta.2005.10.105journaljournalPhys. Lett. A volume352,pages13 (year2006)NoStop
         [Monthus(2006)]Monthus2006
          authorauthorC. Monthus, titletitleRandom
          Walks and Polymers in the Presence of Quenched Disorder,  10.1007/s11005-006-0122-2journaljournalLett. Math. Phys. volume78, pages207 (year2006)NoStop
         [Efrat and Schwartz(2014)]Efrat2014
          authorauthorA. Efrat and authorM. Schwartz, titletitleLack of
          self-averaging in random systems - Liability or asset?  https://doi.org/10.1016/j.physa.2014.06.071journaljournalPhys. A Stat. Mech. Appl. volume414, pages137 (year2014)NoStop
         [Müller and Delande()]MullerARXIV
          authorauthorC. A. Müller and authorD. Delande, @noop titleDisorder and
          interference: localization phenomena, noteLes Houches 2009 -
          Session XCI: Ultracold Gases and Quantum Information, C. Miniatura, L.-C.
          Kwek, M. Ducloy, B. Gremaud, B.-G. Englert, L.F. Cugliandolo, A. Ekert, eds.
          (Oxford University Press, Oxford 2011); arXiv:1005.0915NoStop
         [Serbyn et al.(2017)Serbyn,
          Papi ćć, and Abanin]Serbyn2017
          authorauthorM. Serbyn, authorZ. Papi ćć,  and authorD. A. Abanin, titletitleThouless energy and multifractality across the
          many-body localization transition, 
          10.1103/PhysRevB.96.104201journaljournalPhys.
          Rev. B volume96, pages104201
          (year2017)NoStop
         [Pastur and Slavin(2014)]Pastur2014PRL
          authorauthorL. Pastur and authorV. Slavin, titletitleArea Law
          Scaling for the Entropy of Disordered Quasifree Fermions,  10.1103/PhysRevLett.113.150404journaljournalPhys. Rev. Lett. volume113, pages150404 (year2014)NoStop
         [Wreszinski and Bolina(2004)]Wreszinski2004
          authorauthorW.F. Wreszinski and authorO. Bolina, titletitleA self-averaging
          “order parameter” for the Sherrington-Kirkpatrick spin glass model, 10.1023/B:JOSS.0000041743.24497.63journaljournalJ. Stat. Phys. volume116,pages1389 (year2004)NoStop
         [Milchev et al.(1986)Milchev, Binder, and Heermann]Milchev1986
          authorauthorA. Milchev, authorK. Binder,   and authorD. W. Heermann,titletitleFluctuations and lack of
          self-averaging in the kinetics of domain growth, 
          10.1007/BF01726202journaljournalZ. Phys. B
          Condensed Matter volume63, pages521
          (year1986)NoStop
         [Bouchaud and Georges(1990)]Bouchaud1990
          authorauthorJ.-P. Bouchaud and authorA. Georges, titletitleAnomalous
          diffusion in disordered media: Statistical mechanisms, models and physical
          applications, 
          https://doi.org/10.1016/0370-1573(90)90099-NjournaljournalPhys. Rep. volume195, pages127 (year1990)NoStop
         [Akimoto et al.(2016)Akimoto, Barkai, and Saito]AkimotoPRL2016
          authorauthorT. Akimoto, authorE. Barkai,   and authorK. Saito, titletitleUniversal Fluctuations of
          Single-Particle Diffusivity in a Quenched Environment,  10.1103/PhysRevLett.117.180602journaljournalPhys. Rev. Lett. volume117, pages180602 (year2016)NoStop
         [Russian et al.(2017)Russian, Dentz, and Gouze]Russian2017
          authorauthorA. Russian, authorM. Dentz,   and authorP. Gouze, titletitleSelf-averaging and weak ergodicity
          breaking of diffusion in heterogeneous media, 
          10.1103/PhysRevE.96.022156journaljournalPhys.
          Rev. E volume96, pages022156
          (year2017)NoStop
         [Akimoto et al.(2018)Akimoto, Barkai, and Saito]AkimotoPRE2018
          authorauthorT. Akimoto, authorE. Barkai,   and authorK. Saito, titletitleNon-self-averaging behaviors and
          ergodicity in quenched trap models with finite system sizes,  10.1103/PhysRevE.97.052143journaljournalPhys. Rev. E volume97, pages052143 (year2018)NoStop
         [Mithun et al.(2018)Mithun,
          Kati, Danieli, and Flach]Mithun2018
          authorauthorThudiyangalMithun, authorYagmurKati, authorCarlo Danieli,  and authorSergej Flach, titletitleWeakly nonergodic
          dynamics in the Gross-Pitaevskii lattice, 
          10.1103/PhysRevLett.120.184101journaljournalPhys. Rev. Lett. volume120, pages184101 (year2018)NoStop
         [Deutsch(1991)]Deutsch1991
          authorauthorJ. M. Deutsch, titletitleQuantum
          statistical mechanics in a closed system, @noop journaljournalPhys. Rev. A volume43, pages2046 (year1991)NoStop
         [Srednicki()]SrednickiARXIV
          authorauthorM. Srednicki, @noop titleDoes Quantum
          Chaos Explain Quantum Statistical Mechanics? noteArXiv:9410046NoStop
         [Rigol et al.(2008)Rigol,
          Dunjko, and Olshanii]Rigol2008
          authorauthorM. Rigol, authorV. Dunjko,   and authorM. Olshanii,titletitleThermalization and its
          mechanism for generic isolated quantum systems, http://dx.doi.org/10.1038/nature06838journaljournalNature volume452, pages854 EP – (year2008)NoStop
         [Alessio et al.(2016)Alessio, Kafri, Polkovnikov, andRigol]Dalessio2016
          authorauthorL. D' Alessio, authorY. Kafri,
          authorA. Polkovnikov,  andauthorM. Rigol, titletitleFrom quantum chaos and eigenstate
          thermalization to statistical mechanics and thermodynamics,  10.1080/00018732.2016.1198134journaljournalAdv. Phys. volume65, pages239–362 (year2016)NoStop
         [Reimann(2008)]Reimann2008
          authorauthorPeter Reimann, titletitleFoundation of
          statistical mechanics under experimentally realistic conditions,@noop journaljournalPhys. Rev. Lett.volume101, pages190403 (year2008)NoStop
         [Short(2011)]Short2011
          authorauthorA. J. Short, titletitleEquilibration of
          quantum systems and subsystems, @noop journaljournalNew J. Phys. volume13,pages053009 (year2011)NoStop
         [Short and Farrelly(2012)]Short2012
          authorauthorA. J. Short and authorT. C. Farrelly, titletitleQuantum
          equilibration in finite time, @noop journaljournalNew J. Phys. volume14,pages013063 (year2012)NoStop
         [Zangara et al.(2013)Zangara, Dente, Torres-Herrera,
          Pastawski, Iucci, and Santos]Zangara2013
          authorauthorPablo R.Zangara, authorAxel D.Dente, authorE. J. Torres-Herrera, authorHoracio M.Pastawski, authorA. Iucci,  and authorLea F. Santos, titletitleTime
          fluctuations in isolated quantum systems of interacting particles,@noop journaljournalPhys. Rev. Evolume88, pages032913 (year2013)NoStop
         [Nation and Porras(2018)]Nation2018
          authorauthorCharlie Nation and authorDiego Porras, titletitleOff-diagonal
          observable elements from random matrix theory: distributions, fluctuations,
          and eigenstate thermalization,  10.1088/1367-2630/aae28fjournaljournalNew J. Phys. volume20, pages103003 (year2018)NoStop
         [Łobejko et al.(2018)Łobejko, Dajka, and Łuczka]Lobejko2018
          authorauthorM. Łobejko, authorJ. Dajka,   and authorJ. Łuczka,titletitleSelf-averaging of random
          quantum dynamics,  10.1103/PhysRevA.98.022111journaljournalPhys. Rev. A volume98, pages022111 (year2018)NoStop
         [Mukherjee(2018)]Mukherjee2018
          authorauthorB. Mukherjee, titletitleFloquet
          topological transition by unpolarized light, 
          10.1103/PhysRevB.98.235112journaljournalPhys.
          Rev. B volume98, pages235112
          (year2018)NoStop
         [Leviandier et al.(1986)Leviandier, Lombardi, Jost, andPique]Leviandier1986
          authorauthorLuc Leviandier, authorMaurice Lombardi, authorRémi Jost,
           and authorJean Paul Pique,titletitleFourier transform: A tool to
          measure statistical level properties in very complex spectra,  10.1103/PhysRevLett.56.2449journaljournalPhys. Rev. Lett. volume56, pages2449–2452 (year1986)NoStop
         [Argaman et al.(1993)Argaman, Dittes, Doron, Keating, Kitaev, Sieber, andSmilansky]Argaman1993b
          authorauthorN. Argaman, authorF.-M. Dittes,
          authorE. Doron, authorJ. P. Keating, authorA. Yu. Kitaev, authorM. Sieber,  and authorU. Smilansky, titletitleCorrelations in the actions of periodic orbits derived
          from quantum chaos,  10.1103/PhysRevLett.71.4326journaljournalPhys. Rev. Lett. volume71, pages4326–4329 (year1993)NoStop
         [Eckhardt and Main(1995)]Eckhardt1995
          authorauthorB. Eckhardt and authorJ. Main, titletitleSemiclassical
          Form Factor of Matrix Element Fluctuations, 
          10.1103/PhysRevLett.75.2300journaljournalPhys. Rev. Lett. volume75, pages2300–2303 (year1995)NoStop
         [Prange(1997)]Prange1997
          authorauthorR. E. Prange, titletitleThe spectral
          form factor is not self-averaging, 
          10.1103/PhysRevLett.78.2280journaljournalPhys. Rev. Lett. volume78, pages2280 (year1997)NoStop
         [Braun and Haake(2015)]Braun2015
          authorauthorPetr Braun and authorFritz Haake, titletitleSelf-averaging
          characteristics of spectral fluctuations, 
          10.1088/1751-8113/48/13/135101journaljournalJ.Phys. A volume48, pages135101
          (year2015)NoStop
         [Fonda et al.(1978)Fonda,
          Ghirardi, and Rimini]Fonda1978
          authorauthorL. Fonda, authorG. C. Ghirardi,
           and authorA. Rimini,titletitleDecay theory of unstable
          quantum systems, @noop journaljournalRep. Prog. Phys. volume41, pages587 (year1978)NoStop
         [Bhattacharyya(1983)]Bhattacharyya1983
          authorauthorK. Bhattacharyya, titletitleQuantum
          decay and the Mandelstam-Tamm-energy inequality, @noop journaljournalJ. Phys. A volume16, pages2993 (year1983)NoStop
         [Ufink(1993)]Ufink1993
          authorauthorJ. Ufink, titletitleThe rate of
          evolution of a quantum state,  10.1119/1.17368journaljournalAm. J. Phys. volume61, pages935–936 (year1993)NoStop
         [Wilkie and Brumer(1991)]Wilkie1991
          authorauthorJoshua Wilkie and authorPaul Brumer, titletitleTime-dependent
          manifestations of quantum chaos, 
          10.1103/PhysRevLett.67.1185journaljournalPhys. Rev. Lett. volume67, pages1185–1188 (year1991)NoStop
         [Alhassid and Levine(1992)]Alhassid1992
          authorauthorY. Alhassid and authorR. D. Levine, titletitleSpectral
          autocorrelation function in the statistical theory of energy levels, 10.1103/PhysRevA.46.4650journaljournalPhys. Rev. A volume46, pages4650–4653 (year1992)NoStop
         [Heller()]Heller1991
          authorauthorE. J. Heller, @noop titleWavepacket dynamics
          and quantum chaology, noteLes Houches, Session LII, Quantum
          and Quantum Physics, ed. M.-J. Giannoni, A. Voros, and J. Zinn Justin
          (1989)NoStop
        Khalfin1958  L. A. Khalfin,  “Contribution to
          the decay theory of a quasi-stationary state”, J. Exptl. Theoret. Phys.  33, 1371 (1957) [Sov. Phys. JETP  6, 1053 (1958)].
         [Muga et al.(2009)Muga,
          Ruschhaupt, and del Campo]MugaBook
          authorauthorJ. G. Muga, authorA. Ruschhaupt,   and authorA. del Campo,@noop titleTime in Quantum Mechanics, vol.
          2 (publisherSpringer, addressLondon,year2009)NoStop
         [del Campo(2016)]Campo2016NJP
          authorauthorA. del
          Campo, titletitleExact quantum
          decay of an interacting many-particle system: the
          Calogero–Sutherland model, 
          10.1088/1367-2630/18/1/015014journaljournalN.
          J. Phys. volume18, pages015014
          (year2016)NoStop
         [Távora et al.(2016)Távora, Torres-Herrera, and Santos]Tavora2016
          authorauthorMarco Távora, authorE. J. Torres-Herrera,  and authorLea F.Santos, titletitleInevitable power-law behavior of isolated many-body quantum systems and how
          it anticipates thermalization, 
          10.1103/PhysRevA.94.041603journaljournalPhys.
          Rev. A volume94, pages041603
          (year2016)NoStop
         [Távora et al.(2017)Távora, Torres-Herrera, and Santos]Tavora2017
          authorauthorMarco Távora, authorE. J. Torres-Herrera,  and authorLea F.Santos, titletitlePower-law decay exponents: A dynamical criterion for predicting
          thermalization,  10.1103/PhysRevA.95.013604journaljournalPhys. Rev. A volume95, pages013604 (year2017)NoStop
         [Torres-Herrera and Santos(2014)]Torres2014PRA
          authorauthorE. J. Torres-Herrera and authorLea F.Santos, titletitleQuench dynamics of isolated many-body quantum systems, 
          10.1103/PhysRevA.89.043620journaljournalPhys.
          Rev. A volume89, pages043620
          (year2014)NoStop
         [Torres-Herrera et al.(2014)Torres-Herrera, Vyas, and Santos]Torres2014NJP
          authorauthorE. J. Torres-Herrera, authorM Vyas,
           and authorLea F. Santos,titletitleGeneral features of the
          relaxation dynamics of interacting quantum systems, @noop journaljournalNew J. Phys. volume16, pages063010 (year2014)NoStop
         [Torres-Herrera et al.(2015)Torres-Herrera, Kollmar, and Santos]TorresKollmar2015
          authorauthorE. J. Torres-Herrera, authorD. Kollmar,  and authorL. F. Santos, titletitleRelaxation and
          thermalization of isolated many-body quantum systems, @noop journaljournalPhys. Scr. T volume165, pages014018 (year2015)NoStop
         [Mazza et al.(2016)Mazza,
          Stéphan, Canovi, Alba,
          Brockmann, and Haque]Mazza2016
          authorauthorPaolo P Mazza, authorJean-Marie Stéphan, authorElena Canovi, authorVincenzo Alba,
          authorMichael Brockmann,  andauthorMasudul Haque,titletitleOverlap distributions for
          quantum quenches in the anisotropic Heisenberg chain, 
          10.1088/1742-5468/2016/01/013104journaljournalJ. Stat. Mech. volume2016, pages013104 (year2016)NoStop
         [Volya and Zelevinsky()]VolyaARXIV
          authorauthorA. Volya and authorV. Zelevinsky, @noop titleTime-dependent
          relaxation of observables in complex quantum systems, noteArXiv:1905.11918NoStop
         [Heyl et al.(2013)Heyl,
          Polkovnikov, and Kehrein]Heyl2013
          authorauthorM. Heyl, authorA. Polkovnikov,   and authorS. Kehrein,titletitleDynamical quantum phase
          transitions in the transverse-field ising model, 
          10.1103/PhysRevLett.110.135704journaljournalPhys. Rev. Lett. volume110, pages135704 (year2013)NoStop
         [Santos and Pérez-Bernal(2015)]Santos2015
          authorauthorL. F. Santos and authorF. Pérez-Bernal, titletitleStructure of eigenstates and quench dynamics at an excited-state quantum
          phase transition,  10.1103/PhysRevA.92.050101journaljournalPhys. Rev. A volume92, pages050101 (year2015)NoStop
         [Ketzmerick et al.(1992)Ketzmerick, Petschel, and Geisel]Ketzmerick1992
          authorauthorR. Ketzmerick, authorG. Petschel,  and authorT. Geisel, titletitleSlow decay of
          temporal correlations in quantum systems with Cantor spectra,  10.1103/PhysRevLett.69.695journaljournalPhys. Rev. Lett. volume69, pages695–698 (year1992)NoStop
         [Torres-Herrera and Santos(2015)]Torres2015
          authorauthorE. J. Torres-Herrera and authorLea F.Santos, titletitleDynamics at the many-body localization transition, 
          10.1103/PhysRevB.92.014208journaljournalPhys.
          Rev. B volume92, pages014208
          (year2015)NoStop
         [Bera et al.(2018)Bera,
          De Tomasi, Khaymovich, and Scardicchio]Bera2018
          authorauthorSoumya Bera, authorGiuseppe De Tomasi, authorIvan M. Khaymovich,  and authorAntonelloScardicchio, titletitleReturn probability for the Anderson model on the random regular graph, 10.1103/PhysRevB.98.134205journaljournalPhys. Rev. B volume98, pages134205 (year2018)NoStop
         [Lerma-Hernández et al.(2018)Lerma-Hernández, Chávez-Carlos, Bastarrachea-Magnani, Santos, and Hirsch]Lerma2018
          authorauthorSergio Lerma-Hernández, authorJorgeChávez-Carlos, authorMiguel A Bastarrachea-Magnani, authorLea F Santos,  and authorJorge G Hirsch, titletitleAnalytical description of the survival
          probability of coherent states in regular regimes, 
          10.1088/1751-8121/aae2c3journaljournalJ.
          Phys. A volume51, pages475302
          (year2018)NoStop
         [Reimann(2016)]Reimann2016
          authorauthorPeter Reimann, titletitleTypical fast
          thermalization processes in closed many-body systems, 
          http://dx.doi.org/10.1038/ncomms10821journaljournalNat. Comm. volume7, pages10821 (year2016)NoStop
         [Reimann(2019)]Reimann2019
          authorauthorPeter Reimann, titletitleTransportless
          equilibration in isolated many-body quantum systems, 
          10.1088/1367-2630/ab1a63journaljournalNew J.
          Phys. volume21, pages053014
          (year2019)NoStop
         [Fan et al.(2017)Fan,
          Zhang, Shen, and Zhai]Fan2017
          authorauthorRuihua Fan, authorPengfei Zhang,
          authorHuitao Shen,  andauthorHui Zhai, titletitleOut-of-time-order correlation for
          many-body localization, 
          https://doi.org/10.1016/j.scib.2017.04.011journaljournalScience Bulletin volume62, pages707 – 711 (year2017)NoStop
         [Borgonovi et al.(2019a)Borgonovi, Izrailev, and Santos]Borgonovi2019b
          authorauthorFausto Borgonovi, authorFelix M. Izrailev,  and authorLea F. Santos, titletitleTimescales in
          the quench dynamics of many-body quantum systems: Participation ratio versus
          out-of-time ordered correlator, 
          10.1103/PhysRevE.99.052143journaljournalPhys.
          Rev. E volume99, pages052143
          (year2019a)NoStop
         [Schreiber et al.(2015)Schreiber, Hodgman, Bordia, Lüschen, Fischer, Vosk, Altman, Schneider, and Bloch]Schreiber2015
          authorauthorM. Schreiber, authorS. S. Hodgman, authorPr. Bordia,
          authorH. P. Lüschen,
          authorM. H. Fischer, authorR. Vosk, authorE. Altman, authorU. Schneider,  and authorI. Bloch, titletitleObservation of many-body localization of interacting fermions in a
          quasirandom optical lattice,  10.1126/science.aaa7432journaljournalScience volume349, pages842–845 (year2015)NoStop
         [Richerme et al.(2014)Richerme, Gong, Lee, Senko,
          Smith, Foss-Feig, Michalakis,
          Gorshkov, and Monroe]Richerme2014
          authorauthorP. Richerme, authorZ.-X. Gong,
          authorA. Lee, authorCr. Senko, authorJ. Smith, authorM. Foss-Feig, authorS. Michalakis, authorA. V.Gorshkov,  and authorC. Monroe, titletitleNon-local propagation of correlations in quantum systems with long-range
          interactions, @noop journaljournalNature volume511, pages198–201
          (year2014)NoStop
         [Schiulaz et al.(2019)Schiulaz, Torres-Herrera, and Santos]Schiulaz2019
          authorauthorMauro Schiulaz, authorE. Jonathan Torres-Herrera,  and authorLea F.Santos, titletitleThouless and relaxation time scales in many-body quantum systems,  10.1103/PhysRevB.99.174313journaljournalPhys. Rev. B volume99, pages174313 (year2019)NoStop
         [Santos et al.(2004)Santos,
          Rigolin, and Escobar]SantosEscobar2004
          authorauthorL. F. Santos, authorG. Rigolin,   and authorC. O. Escobar,titletitleEntanglement versus chaos in
          disordered spin systems, @noop journaljournalPhys. Rev. A volume69, pages042304 (year2004)NoStop
         [Dukesz et al.(2009)Dukesz,
          Zilbergerts, and Santos]Dukesz2009
          authorauthorF. Dukesz, authorM. Zilbergerts,
           and authorL. F. Santos,titletitleInterplay between
          interaction and (un)correlated disorder in one-dimensional many-particle
          systems: delocalization and global entanglement, 
          10.1088/1367-2630/11/4/043026journaljournalNew
          J. Phys. volume11, pages043026
          (year2009)NoStop
         [Pal and Huse(2010)]Pal2010
          authorauthorA. Pal and authorD. A. Huse,titletitleMany-body localization phase
          transition,  10.1103/PhysRevB.82.174411journaljournalPhys. Rev. B volume82, pages174411 (year2010)NoStop
         [Luitz and Lev(2017)]Luitz2017
          authorauthorD. Luitz and authorY. Bar Lev, titletitleThe ergodic side of
          the many-body localization transition, 
          10.1002/andp.201600350journaljournalAnn.
          Phys.(Berlin) volume529, pages1600350 (year2017)NoStop
         [Santos et al.(2012)Santos,
          Borgonovi, and Izrailev]Santos2012PRE
          authorauthorL. F. Santos, authorF. Borgonovi,   and authorF. M. Izrailev,titletitleOnset of chaos and
          relaxation in isolated systems of interacting spins: Energy shell
          approach,  10.1103/PhysRevE.85.036209journaljournalPhys. Rev. E volume85, pages036209 (year2012)NoStop
         [Chan et al.(2018)Chan,
          De Luca, and Chalker]Chan2018PRX
          authorauthorA. Chan, authorA. De Luca,   and authorJ. T. Chalker,titletitleSolution of a Minimal
          Model for Many-Body Quantum Chaos, 
          10.1103/PhysRevX.8.041019journaljournalPhys.
          Rev. X volume8, pages041019
          (year2018)NoStop
         [Bertini et al.(2018)Bertini, Kos, and Prosen]Bertini2018
          authorauthorB. Bertini, authorP. Kos,  andauthorT. Prosen, titletitleExact Spectral Form Factor in a
          Minimal Model of Many-Body Quantum Chaos, 
          10.1103/PhysRevLett.121.264101journaljournalPhys. Rev. Lett. volume121, pages264101 (year2018)NoStop
         [Suntajš et al.()Suntajš, Bonča, Prosen, andVidmar]SuntajsARXIV
          authorauthorJ. Suntajš, authorJ. Bonča, authorT. Prosen,   and authorL. Vidmar,@noop titleQuantum chaos challenges
          many-body localization, noteArXiv:1905.06345NoStop
         [Friedman et al.()Friedman,
          Chan, Luca, and Chalker]FriedmanARXIV
          authorauthorAaron J.Friedman, authorAmosChan, authorAndrea DeLuca,  and authorJ. T.Chalker, @noop titleSpectral statistics and many-body quantum chaos with conserved charge,noteArXiv:1906.07736NoStop
         [Borgonovi et al.(2019b)Borgonovi, Izrailev, and Santos]Borgonovi2019
          authorauthorFausto Borgonovi, authorFelix M. Izrailev,  and authorLea F. Santos, titletitleExponentially
          fast dynamics of chaotic many-body systems, 
          10.1103/PhysRevE.99.010101journaljournalPhys.
          Rev. E volume99, pages010101
          (year2019b)NoStop
        Torres2016 E. J. Torres-Herrera, J. Karp, M. Távora, and L. F. Santos,
        Entropy  18, 359 (2016).  
         [Torres-Herrera and Santos(2017)]Torres2017Philo
          authorauthorE. J. Torres-Herrera and authorLea F.Santos, titletitleDynamical manifestations of quantum chaos: correlation hole and bulge, 10.1098/rsta.2016.0434journaljournalPhilos. Trans. Royal Soc. A volume375,pages20160434 (year2017)NoStop
         [Arranz et al.(2019)Arranz,
          Benito, and Borondo]Arranz2019
          authorauthorF. J. Arranz, authorR. M. Benito,
           and authorF. Borondo,titletitleShannon entropy at avoided
          crossings in the quantum transition from order to chaos, 
          10.1103/PhysRevE.99.062209journaljournalPhys.
          Rev. E volume99, pages062209
          (year2019)NoStop
         [Torres-Herrera et al.(2018)Torres-Herrera, García-García, andSantos]Torres2018
          authorauthorE. J. Torres-Herrera, authorAntonio M.García-García,  and authorLea F. Santos, titletitleGeneric dynamical features of quenched
          interacting quantum systems: Survival probability, density imbalance, and
          out-of-time-ordered correlator, 
          10.1103/PhysRevB.97.060303journaljournalPhys.
          Rev. B volume97, pages060303
          (year2018)NoStop
         [Santos and Torres-Herrera(2017)]SantosTorres2017AIP
          authorauthorL. F. Santos and authorE. J. Torres-Herrera, titletitleAnalytical expressions for the evolution of many-body quantum systems
          quenched far from equilibrium,  10.1063/1.5016140journaljournalAIP Conference Proceedingsvolume1912, pages020015 (year2017)NoStop
         [Zelevinsky et al.(1996)Zelevinsky, Brown, Frazier, andHoroi]ZelevinskyRep1996
          authorauthorV. Zelevinsky, authorB. A. Brown, authorN. Frazier,   and authorM. Horoi, titletitleThe nuclear shell model as a testing
          ground for many-body quantum chaos, @noop journaljournalPhys. Rep. volume276,pages85–176 (year1996)NoStop
         [Flambaum and Izrailev(2001)]Flambaum2001a
          authorauthorV. V. Flambaum and authorF. M. Izrailev, titletitleUnconventional
          decay law for excited states in closed many-body systems, @noop journaljournalPhys. Rev. E volume64, pages026124 (year2001)NoStop
         [Kota and Sahu(2001)]Kota2001PRE
          authorauthorV. K. B.Kota and authorR. Sahu, titletitleStructure of wave
          functions in (1+2)-body random matrix ensembles, 
          10.1103/PhysRevE.64.016219journaljournalPhys.
          Rev. E volume64, pages016219
          (year2001)NoStop
         [foo()]footOrtho
          @noop noteOrthogonality actually causes deviations
          that are subleading in 1/D, so they can be neglected in our
          calculations <cit.>NoStop
         [Cotler et al.(2017)Cotler,
          Hunter-Jones, Liu, and Yoshida]Cotler2017GUE
          authorauthorJ. Cotler, authorN. Hunter-Jones, authorJ. Liu,   and authorB. Yoshida,titletitleChaos, complexity, and
          random matrices,  10.1007/JHEP11(2017)048journaljournalJ. High Energy Phys. volume2017, pages48 (year2017)NoStop
