Equilibration time in many-body quantum systems
Talía L. M. Lezama, E. Jonathan Torres-Herrera, Francisco Pérez-Bernal, Yevgeny Bar Lev, Lea F. Santos
TThermalization time in many-body quantum systems
Tal´ıa L. M. Lezama, E. Jonathan Torres-Herrera, Francisco P´erez-Bernal, Yevgeny Bar Lev, and Lea F. Santos Department of Physics, Ben-Gurion University of the Negev, Beer-Sheva 84105, Israel Instituto de F´ısica, Benem´erita Universidad Aut´onoma de Puebla, Apt. Postal J-48, Puebla, 72570, Mexico Dep. CC. Integradas y Centro de Estudios Avanzados en F´ısica,Matem´aticas y Computaci´on. Fac. CC. Experimentales, Universidad de Huelva, Huelva 21071,& Instituto Carlos I de F´ısica Te´orica y Computacional, Universidad de Granada, Granada 18071, Spain Department of Physics, Yeshiva University, New York, New York 10016, USA
Isolated chaotic many-body quantum systems can reach thermal equilibrium, but it is not yet clearhow long they take to do so. To answer this question, we use exact numerical methods and analyzethe entire evolution, from perturbation to thermalization, of a paradigmatic disordered many-bodyquantum system in the chaotic regime. We investigate how the thermalization time depends on thesystem size and observables. We show that if dynamical manifestations of spectral correlations inthe form of the correlation hole (“ramp”) are taken into account, the time for thermalization scalesexponentially with system size, while if they are neglected, the scaling is better described by a powerlaw with system size, though with an exponent larger than expected for diffusive transport.
I. INTRODUCTION
One major question in studies of nonequilibrium dynamics of isolated many-body quantum systems is how longit takes for an experimentally relevant observable to reach equilibrium, that is, to reach a point beyond which itsexpectation value simply fluctuates around its infinite-time average and the size of these fluctuations decreases as thesystem size increases [1–6]. Equilibration in this sense does not necessarily imply thermalization.The variety of approaches taken to address this question and the lack of agreement among the existing results areworrisome. Several analyses are based on assumptions about the spectrum, observables, and initial conditions, andoften provide bounds for the equilibration time. Some suggest that this time should decrease with system size [7, 8],others that it should depend weakly on it [9], and others yet that it should increase with it [3, 4, 10–16], possiblyexponentially [7, 17]. Studies aligned with transport behavior [18–31], on the other hand, expect the equilibrationtime to increase as a power law with system size.Confronted by so many options, it is worth to step back and try first to identify a general scenario. For this purpose,we focus specifically on the time for thermalization, examining many-body quantum systems that are in the chaoticregime. In this case, initial states far from equilibrium and with energy expectation values away from the edges ofthe many-body spectrum lead to thermal equilibrium [32–34]. This allows us to avoid the particularities of integrablemodels and specific initial states.The largest possible timescale in quantum systems is given by the inverse of their mean-level spacing, the so-calledHeisenberg time, which grows linearly with the dimension of the Hilbert space, and thus exponentially with systemsize. The Heisenberg time is the absolute upper bound for the thermalization time, but do experimental observablestake this long to thermalize? This is the main question addressed by this work. The answer is yes [35] if the dynamicalmanifestations of spectral correlations, known as correlation hole [35–42] and sometimes referred to as “ramp” [43–45],are taken into account. However, as we show here, these manifestations can be neglected for some observables, leadinginstead to a power-law scaling of the thermalization time with system size, a result that is in better agreement withstudies of transport behavior [31, 34].We investigate the evolution of the survival (return) probability, which is the probability of finding the system inits initial state later in time and may be accessible to experiments with cold atoms [46, 47]; the inverse participationratio, which quantifies the spreading in time of the initial state over the many-body Hilbert space and whose loga-rithm gives the participation second-order R´enyi entropy; the spin autocorrelation function, which is related to theimbalance measured in experiments with cold atoms [48]; and the connected spin-spin correlation function measuredin experiments with ion traps [49]. Our analysis focuses on the disordered spin-1/2 Heisenberg chain in the chaoticregime and we briefly discuss the dependence of the thermalization time on disorder strength.The article is organized as follows. Section II introduces the disordered model, initial conditions, and observables. InSec. III, we show the timescales for the correlation hole and for the thermalization after it. In Sec. IV, the correlationhole is neglected and a new definition for the thermalization time is proposed. Its dependence on the disorder strengthis also provided. Conclusions are presented in Sec. V. a r X i v : . [ c ond - m a t . d i s - nn ] F e b II. MODEL, INITIAL STATES, AND OBSERVABLESA. Model
The disordered spin-1/2 Heisenberg chain that we consider is a representative model of disordered interacting one-dimensional systems and has been extensively used in studies of many-body localization [50–53]. It is given by thefollowing Hamiltonian, ˆ H = L (cid:88) i =1 J ˆS i · ˆ S i +1 + h i ˆ S zi , (1)where ˆS i ≡ (cid:16) ˆ S xi , ˆ S yi , ˆ S zi (cid:17) are spin-1/2 operators, L is the system size, periodic boundary conditions are assumed,we set J = 1 and h i to be independent and uniformly distributed random variables in [ − W, W ], with W being theonsite disorder strength. The system conserves the total magnetization in the z -direction, ˆ S z tot = (cid:80) Li ˆ S zi , and exhibitsa transition to the many-body localized phase at a critical disorder strength of W c . The value of W c is still underdebate [54–63], some papers estimate that 3 < W c < W c > S z tot = 0 subspace that has the Hilbert space dimension D = (cid:0) LL/ (cid:1) , and consider finite systemsaway from the critical region, 0 . ≤ W (cid:46)
1, although a short discussion for values of disorder closer to the criticalregion 1 < W <
B. Initial states and thermalization
We use initial states | Ψ(0) (cid:105) that are product states in the S z -basis, | φ n (cid:105) , such as | ↑↓↓ . . . ↑(cid:105) , which can beexperimentally realized. We choose initial states with an energy expectation value (cid:104) Ψ(0) | ˆ H | Ψ(0) (cid:105) far from the edgesof the spectrum to guarantee that they fall in the chaotic region of the many-body Hamiltonian and that a givenfew-body observable ˆ O reaches thermal equilibrium. Its infinite-time average expressed in terms of the many-bodyeigenstates ˆ H | α (cid:105) = E α | α (cid:105) is given by ¯ O = (cid:88) α | C α | (cid:104) α | ˆ O | α (cid:105) , (2)where C α ≡ (cid:104) α | Ψ(0) (cid:105) .We use either exact diagonalization or Krylov-space methods to evolve the system in time. All results are averagedover 10 samples composed of 0 . D initial states and 10 / (0 . D ) disorder realizations. The average over samplesis denoted by the angle brackets (cid:104)·(cid:105) . We have also verified that fixing a single initial state and using 10 disorderrealizations leads to equivalent results, though the numerical procedure is less efficient. All statistical errors areaccessed using a bootstrap procedure. C. Observables
We study the time-evolution of two non-local (in real space) observables, the survival probability and the inverseparticipation ratio, and two local quantities, the spin autocorrelation function and the connected spin-spin correlationfunction.
1. Non-local Quantities
The survival probability is defined as P S ( t ) = |(cid:104) Ψ(0) | Ψ( t ) (cid:105)| . (3)Taking the mean gives, (cid:104) P S ( t ) (cid:105) = (cid:104) (cid:88) α (cid:54) = β | C α | | C β | e − i ( E α − E β ) t (cid:105) + (cid:104) (cid:88) α | C α | (cid:105) , (4)which is related to the spectral form factor, (cid:104) (cid:80) α (cid:54) = β e − i ( E α − E β ) t (cid:105) . While (cid:104) P S ( t ) (cid:105) is a true dynamical quantity, whichdepends on the initial state through the components C α , the spectral form factor involves only the eigenvalues of theHamiltonian and is used to study the manifestation of level statistics in the time domain [64]. Filter functions aresometimes added to it [45, 65], but they do not come from the quench dynamics, as the coefficients C α in (cid:104) P S ( t ) (cid:105) .The second non-local quantity that we consider is the inverse participation ratio,I PR ( t ) = (cid:88) n |(cid:104) φ n | Ψ( t ) (cid:105)| , (5)which measures the spreading in time of the initial state | Ψ(0) (cid:105) over the many-body Hilbert space. Its logarithm, − ln I PR ( t ), corresponds to the second-order R´enyi entropy.
2. Local Observables
The spin autocorrelation function measures the proximity of the spin configuration in the z -direction at a time t tothe initial spin configuration, I ( t ) = 4 L L (cid:88) i =1 (cid:104) Ψ | ˆ S zi (0) ˆ S zi ( t ) | Ψ (cid:105) , (6)where ˆ S zi ( t ) = e iHt ˆ S zi (0) e − iHt . For the N´eel initial state, | ↑↓↑↓↑↓ . . . (cid:105) , it reduces to the density imbalance betweeneven and odd sites that is measured in experiments with cold atoms [48].The connected spin-spin correlation function is defined as C ( t ) = 4 L (cid:88) i =1 (cid:104) (cid:104) Ψ( t ) | ˆ S zi ˆ S zi +1 | Ψ( t ) (cid:105) − (cid:104) Ψ( t ) | ˆ S zi | Ψ( t ) (cid:105) (cid:104) Ψ( t ) | ˆ S zi +1 | Ψ( t ) (cid:105) (cid:105) (7)and has been measured in experiments with ion traps [49]. III. THERMALIZATION AFTER THE CORRELATION HOLE
A semi-analytical expression for the entire evolution of the average of the survival probability was derived forrealistic chaotic many-body quantum systems with short-range interactions [35, 42], as the one described by Eq. (1),and it is given by (cid:104) P S ( t ) (cid:105) = 1 − (cid:10) P S (cid:11) ( D − (cid:20) Db (Γ t ) − b (cid:18) Γ t √ πD (cid:19)(cid:21) + (cid:10) P S (cid:11) , (8)where Γ = (cid:104) Ψ(0) | ˆ H | Ψ(0) (cid:105) − (cid:104)
Ψ(0) | ˆ H | Ψ(0) (cid:105) (9)is the squared width of the energy distribution of the initial state, (cid:10) P S (cid:11) is the mean of the infinite-time average of P S ( t ), b (Γ t ) = e − Γ t N (cid:12)(cid:12)(cid:12)(cid:12) erf (cid:18) E max + it Γ √ (cid:19) − erf (cid:18) E min + it Γ √ (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) , (10) N is a normalization constant, erf is the error function, E max and E min are the largest and smallest eigenvalues of ˆ H ,respectively, and b ( t ) = − t + t ln(2 t + 1) t ≤ t ln (cid:18) t + 12 t − (cid:19) − t > , (11)is the two-level form factor. t − − − − h P S i (a) L t − − − − h I P R i W = 0 . (b) t − − − h I i (c) t − h C i (d)
10 12 14 16 18 L . . . . . κ (e) W . . . L κ × − (f)
10 12 14 16 18 L . . . . . κ (g)
10 12 14 16 18 L . . . . . κ × − (h) FIG. 1.
Upper panels : Time evolution of the mean (a) survival probability (cid:104) P S ( t ) (cid:105) , (b) inverse participation ratio (cid:104) I PR ( t ) (cid:105) ,(c) spin autocorrelation function (cid:104)I ( t ) (cid:105) , and (d) connected spin-spin correlation function (cid:104) C ( t ) (cid:105) for different system sizes, asindicated in panel (a), and for disorder strength W = 0 .
5, as shown in panel (b). The horizontal dotted-dashed lines mark theasymptotic values and the stars indicate the crossing time t ∗ . In Fig. 1 (a), the numerical data overlap with the curves forthe semi-analytical expression in Eq. (8). Lower panels : Relative depth of the correlation hole, κ , as a function of L , for (e) (cid:104) P S ( t ) (cid:105) , (f) (cid:104) I PR ( t ) (cid:105) , (g) (cid:104)I ( t ) (cid:105) , and (h) (cid:104) C ( t ) (cid:105) for three values of the disorder strength in the chaotic regime, as indicated inFig. 1 (e). The horizontal dashed line in panel (e) corresponds to the value κ = 1 / The decay of the survival probability is controlled by b (Γ t ), which at the time t m ∝ D / / Γ meets the b functionat the minimum value of the survival probability, ∼ /D . After this point, the evolution is described entirely by the b function, which brings (cid:104) P S ( t ) (cid:105) from its minimum up to the saturation value at ∼ /D . Saturation happens at theHeisenberg time t H ∝ D/ Γ.The time interval governed by the b function, where (cid:104) P S ( t ) (cid:105) < (cid:10) P S (cid:11) , is known as correlation hole [35–42] or“ramp” [43–45]. It is a dynamical manifestation of spectral correlations, and, therefore, a dynamical signature ofquantum chaos. The point in time where the ramp starts, t m , has been referred to as Thouless time [35, 65]. Itcoincides with the time where the inverse participation ratio reaches its minimum value [35], indicating that t m is thetime that it takes for the initial state to maximally spread over the many-body Hilbert space and acquire weight overthe unperturbed many-body states | φ n (cid:105) in its microcanonical energy shell given by the width Γ.The evolution of the mean survival probability for the spin model (1) is shown in Fig. 1 (a) for different system sizes.There is an excellent agreement between the numerical results and expression (8) when W = 0 .
5, which correspondsto the deep chaotic regime. The correlation hole is evident in all the curves and it does not fade away as L increases.This means that the complete equilibration of this quantity takes place only after the hole ends at the Heisenbergtime t H . Since this time is exponentially long in the system size L , we use exact diagonalization to resolve the entiredynamics, which limits the accessible systems sizes to L = 18.A correlation hole is also visible for the spin autocorrelation function, as depicted in Fig. 1 (c), suggesting that forsufficiently small system sizes, where it develops for times that are not exceedingly long and reaches minimum valuesthat are not too small, the hole might be experimentally detected.In contrast to the survival probability and the spin autocorrelation function, the effects of the spectral correlationsin the evolution of (cid:104) I PR ( t ) (cid:105) [Fig. 1 (b)] and of (cid:104) C ( t ) (cid:105) [Fig. 1 (d)] are minor and the correlation hole is hardly discernible.Furthermore, for (cid:104) I PR ( t ) (cid:105) , (cid:104) C ( t ) (cid:105) , and also (cid:104)I ( t ) (cid:105) , the hole gets smaller as the system size increases, which justifiesneglecting it in the definition of the thermalization time for these quantities.The relative depth of the correlation hole for a given observable ˆ O can be measured by the ratio [41, 66, 67] κ = (cid:104) O (cid:105) − (cid:104) O (cid:105) min (cid:104) O (cid:105) , (12)where (cid:104) O (cid:105) min stands for the value of the observable at the minimum of the correlation hole. For the survival probability, κ converges to a constant, while it decreases with system size for the other quantities.When the survival probability is evolved using full random matrices taken from a Gaussian orthogonal ensemble(GOE), it is known analytically that (cid:10) P S (cid:11) ∼ /D and (cid:104) P S (cid:105) min ∼ /D [35, 39], which yields κ = 1 /
3. In Fig. 1 (e),we show κ for the survival probability of the spin model as a function of system size for different disorder strengths − − − h I i (a) l og ( t ∗ ) (b) W . .
75 1 . .
25 1 . . . . . . . γ (c) − t − h C i L = 18 (d) . . . . L )123 l og ( t ∗ ) (e) . . . . . . . W . . . . . . γ (f) FIG. 2.
Upper [Lower] panels : Analysis of the crossing time for the spin autocorrelation function (cid:104)I ( t ) (cid:105) [connected spin-spin autocorrelation function (cid:104) C ( t ) (cid:105) ]. (a) [(d)] Evolution of the mean spin autocorrelation function [connected spin-spinautocorrelation function] for L = 18. Each curve corresponds to a different disorder strength W , as indicated. The horizontaldotted-dashed lines represent the asymptotic value of the observable and the crossing time t ∗ is marked with a star. (b) [(e)]Scaling of the crossing time t ∗ with L . Symbols are the data and solid lines give power-law fits to t ∗ ∝ L γ . Error bars indicatethe standard deviation of the bootstrap procedure. (c) [(f)] Exponent γ extracted from the power-law fit to the data points in(b) [(e)] as a function of W . Error bars indicate the standard deviation on the fitted exponent. in the chaotic regime. The relative depth clearly converges to κ = 1 /
3, which is indicated with a horizontal dashedline. The correlation hole is therefore a robust property of the survival probability.Contrary to the survival probability, the relative depth κ for (cid:104) I PR ( t ) (cid:105) [Fig. 1 (f)], (cid:104)I ( t ) (cid:105) [Fig. 1 (g)], and (cid:104) C ( t ) (cid:105) [Fig. 1 (h)] decays with L . In the case of the inverse participation ratio and the connected spin-spin correlationfunction, the decay is exponential, while for the spin autocorrelation function, the results are more subtle and makeapparent the danger of finite-size effects. While for W = 0 . κ decreases monotonically with L , for W = 0 .
75 and W = 1, κ increases for small values of L and the decay becomes clear only for L >
IV. THERMALIZATION BEFORE THE CORRELATION HOLE
Ignoring the correlation hole, we define the thermalization time as the point where (cid:104) ˆ O ( t ) (cid:105) first crosses the infinite-time average (cid:104) O (cid:105) . We denote this time by t ∗ , which is clearly much smaller than the Heisenberg time, t ∗ (cid:28) t H . Thesecrossing points are marked with stars in Figs. 1 (a)-(d). A. Weak disorder region: . ≤ W ≤ . ≤ W ≤ . ≤ W ≤ For sufficiently weak disorder, W = 0 . t ∗ on system size. At long times, disregarding the correlation hole, the decay of (cid:104) P S ( t ) (cid:105) is givenby [69, 70] (cid:10) P S ( t (cid:29) Γ − ) (cid:11) decay (cid:39) e − E / Γ + e − E / Γ π N Γ t . (13)Using that E min , E max and Γ are extensive, namely proportional to the size of the system, and that the survivalprobability saturates at (cid:10) P S (cid:11) ∼ /D , we find that t ∗ ∝ exp(0 . L ), which agrees very well with the crossing timeobtained numerically for L = 10 , , , , , ,
22. Knowing the saturation value of the survival probability, wecan evolve | Ψ( t ) (cid:105) up to a vicinity of t ∗ only, which is a major saving compared to the evolution up to t H . This allowsus to use Krylov-space methods for the time evolution of the survival probability and access to system sizes L > t ∗ with L is indeed better than apower-law scaling. . . . . . W l og ( t ∗ ) (a) h P S ( t ∗ ) ih I PR ( t ∗ ) i hI ( t ∗ ) ih C ( t ∗ ) i . . . . | W − W c | )3 . . . . l og ( t ∗ ) (b) FIG. 3. Dependence of the crossing time t ∗ on (a) [(b)] the disorder strength W [the difference | W − W c | ; with W c = 3 .
7] forthe survival probability (cid:104) P S ( t ) (cid:105) , inverse participation ratio (cid:104) I PR ( t ) (cid:105) , spin autocorrelation function (cid:104)I ( t ) (cid:105) , and the connectedspin-spin correlation function (cid:104) C ( t ) (cid:105) for system size L = 16. The symbols representing each quantity are shown in panel (a).The dashed [dotted-dashed] lines in (a) [(b)] correspond to exponential [power-law] fits obtained with the points for W ≥ . While an analytical expression is not available for the time-dependence of the inverse participation ratio, in theweak disorder regime, its saturation value is known analytically, (cid:10) I PR (cid:11) ∼ /D [71], so we can also obtain t ∗ for systemsizes L >
18. To do that, we apply a Savitzky-Golay filter to smooth the curves for (cid:104) I PR ( t ) (cid:105) and then extract thecrossing time [72]. For this quantity, we find that a power-law scaling of t ∗ with L actually works better than anexponential scaling. This makes us suspect that the exponential dependence of the crossing time with system sizefound for the survival probability may be related to the prevalence of the correlation hole, a conjecture that is furtherreinforced by the next results.For the two local quantities, we do not have analytical results for the saturation values, hence our analysis isrestricted to five system sizes. In Fig. 2 (a) and Fig. 2 (d), we fix the system size at L = 18 and mark the crossingtime for the curves of (cid:104)I ( t ) (cid:105) and (cid:104) C ( t ) (cid:105) obtained for different disorder strengths. In Fig. 2 (b) and Fig. 2 (e), we presentthe scaling analysis for both quantities and those values of W . We find that, similarly to the inverse participationratio, the dependence of t ∗ with L is better fitted with a power law than an exponential, that is t ∗ ∝ L γ with γ > . (14)The exact value of the exponent γ varies with the disorder strength, as shown in Fig. 2 (c) and Fig. 2 (f). For disorderstrength 0 . ≤ W ≤
1, we find that on average γ ≈ . ± .
3. The value of the exponent is somewhat close to thatargued in [28, 73], but larger than the value which typically appears in studies of transport behavior [31, 34].
B. Near critical region: < W < W c < W < W c < W < W c Figure 2 shows the crossing time t ∗ of the local quantities computed for disorder strengths slightly larger than thecoupling parameter, W (cid:38) J . The power-law fit is still better than the exponential one, although γ [Eq. (14)] for W = 1 . W > t ∗ depends on the disorder strength for afixed L .It was shown that the time for the minimum of the correlation hole, t m , for the survival probability and for thespin autocorrelation function grows exponentially as the disorder strength increases [35], which was later confirmed interms of a spectral analysis [65]. It is thus relevant to examine the dependence of t ∗ on W , which we do in Fig. 3. Thebehavior of the crossing time t ∗ is not conclusive as can be seen from Fig. 3. Considering disorder strengths W ∈ [1 , t ∗ is best fitted with a stretched exponential, but if we consider disorder strengths closer to the critical point, W ≥ .
75, we see that either an exponential dependence t ∗ ∝ exp( W/W (cid:48) ) or a critical form t ∗ ∝ | W − W c | − β with anexponent β ≈ . ± .
08 describe the data reasonably well for all the quantities, except for the inverse participationratio, which exhibits strong fluctuations in its corresponding crossing times. The quality of the fits varies depending onthe quantity: the survival probability is better described by an exponential form [Fig. 3 (a)], while the local quantitiesseem to have critical scaling [Fig. 3 (b)]. We did not see a qualitative change in our conclusions when varying W c from 3.5 to 6. V. CONCLUSIONS
We investigated how the thermalization time of the disordered spin-1/2 Heisenberg chain depends on the systemsize L for four different observables. If the correlation hole is taken into account when defining the thermalizationtime, then it coincides with the Heisenberg time and thus grows exponentially with L . This is what happens forthe survival probability, where the correlation hole persists in the thermodynamic limit. However, for the inverseparticipation ratio, spin autocorrelation function, and connected spin-spin correlation function, the correlation holefades away as L increases, which justifies neglecting it. In this case, we defined the thermalization time as the pointwhere the evolution of the observables first crosses their infinite-time averages. The dependence of this crossing timeon system size is best described by a power law.Chaotic systems with static Hamiltonians conserve at least the total energy, have diffusive energy transport, andtheir thermalization time, namely the time it takes for a non-uniformity in the energy density to spread across thesystem, is expected to be bounded from below by L . For disordered chaotic systems transport is expected to besubdiffusive [18, 19, 24, 29, 58], and their thermalization time to be bounded from below by L γ with γ >
2, thoughfor very weak disorder, γ is expected to approach the lower bound, γ ≈
2, where transport is indistinguishable fromdiffusion. Interestingly, even for the lowest disorder that we study, W ≈ .
5, the thermalization time scales as L γ with γ >
3, so it is parametrically larger than the time it takes to make the energy density homogeneous. We leavefor future studies the question of how the system is away from thermal equilibrium during the times L < t < L .We have provided a brief analysis of the dependence of the crossing time t ∗ on the disorder strength away fromthe critical point. Although it is hard to discern between an exponentially growing t ∗ with W and a critical behavior t ∗ ∝ | W − W c | − β , the survival probability seems to be better described by the former, while the local quantitiesappear to show the critical scaling with an exponent β ≈ . ACKNOWLEDGMENTS
This research was supported by a grant from the United States-Israel Binational Foundation (BSF, Grant No.2019644), Jerusalem, Israel, and the United States National Science Foundation (NSF, Grant No. DMR-1936006).T.L.M.L. acknowledges funding from the Kreitman fellowship. E.J.T.-H. is grateful to LNS-BUAP for their super-computing facility. Computing resources supporting this work were provided by the CEAFMC and Universidad deHuelva High Performance Computer (HPC@UHU) located in the Campus Universitario el Carmen and funded byFEDER/MINECO project UNHU-15CE-2848. F.P.B. thanks Spanish National Research, Development, and Innova-tion plan (RDI plan) under the project PID2019-104002GB-C21 and the Consejer´ıa de Conocimiento, Investigaci´ony Universidad, Junta de Andaluc´ıa and European Regional Development Fund (ERDF), refs. SOMM17/6105/UGRand UHU-1262561. [1] A. Peres, Stability of quantum motion in chaotic and regular systems, Phys. Rev. A , 1610 (1984).[2] H. Tasaki, From quantum dynamics to the canonical distribution: General picture and a rigorous example, Phys. Rev.Lett. , 1373 (1998).[3] P. Reimann, Foundation of statistical mechanics under experimentally realistic conditions, Phys. Rev. Lett. , 190403(2008).[4] A. J. Short, Equilibration of quantum systems and subsystems, New J. Phys. , 053009 (2011).[5] P. R. Zangara, A. D. Dente, E. J. Torres-Herrera, H. M. Pastawski, A. Iucci, and L. F. Santos, Time fluctuations in isolatedquantum systems of interacting particles, Phys. Rev. E , 032913 (2013). [6] T. Kiendl and F. Marquardt, Many-particle dephasing after a quench, Phys. Rev. Lett. , 130601 (2017).[7] S. Goldstein, T. Hara, and H. Tasaki, Time scales in the approach to equilibrium of macroscopic quantum systems, Phys.Rev. Lett. , 140401 (2013).[8] S. Goldstein, T. Hara, and H. Tasaki, Extremely quick thermalization in a macroscopic quantum system for a typicalnonequilibrium subspace, New J. Phys. , 045002 (2015).[9] T. R. de Oliveira, C. Charalambous, D. Jonathan, M. Lewenstein, and A. Riera, Equilibration time scales in closedmany-body quantum systems, New J. Phys. , 033032 (2018).[10] P. Reimann and M. Kastner, Equilibration of isolated macroscopic quantum systems, New J. Phys. , 043020 (2012).[11] P. Reimann, Typical fast thermalization processes in closed many-body systems, Nat. Commun. , 10821 (2016).[12] A. J. Short and T. C. Farrelly, Quantum equilibration in finite time, New J. Phys. , 013063 (2012).[13] T. Monnai, Generic evaluation of relaxation time for quantum many-body systems: Analysis of the system size dependence,J. Phys. Soc. Jpn. , 044006 (2013).[14] D. Hetterich, M. Fuchs, and B. Trauzettel, Equilibration in closed quantum systems: Application to spin qubits, Phys.Rev. B , 155314 (2015).[15] C. Gogolin and J. Eisert, Equilibration, thermalisation, and the emergence of statistical mechanics in closed quantumsystems, Rep. Prog. Phys. , 056001 (2016).[16] L. P. Garc´ıa-Pintos, N. Linden, A. S. L. Malabarba, A. J. Short, and A. Winter, Equilibration time scales of physicallyrelevant observables, Phys. Rev. X , 031027 (2017).[17] A. S. L. Malabarba, L. P. Garc´ıa-Pintos, N. Linden, T. C. Farrelly, and A. J. Short, Quantum systems equilibrate rapidlyfor most observables, Phys. Rev. E , 012121 (2014).[18] Y. Bar Lev, G. Cohen, and D. R. Reichman, Absence of diffusion in an interacting system of spinless fermions on aone-dimensional disordered lattice, Phys. Rev. Lett. , 100601 (2015).[19] K. Agarwal, S. Gopalakrishnan, M. Knap, M. M¨uller, and E. Demler, Anomalous diffusion and griffiths effects near themany-body localization transition, Phys. Rev. Lett. , 160401 (2015).[20] A. C. Potter, R. Vasseur, and S. A. Parameswaran, Universal properties of many-body delocalization transitions, Phys.Rev. X , 031033 (2015).[21] R. Vosk, D. A. Huse, and E. Altman, Theory of the many-body localization transition in one-dimensional systems, Phys.Rev. X , 031032 (2015).[22] M. Serbyn, Z. Papi´c, and D. A. Abanin, Criterion for many-body localization-delocalization phase transition, Phys. Rev.X , 041047 (2015).[23] M. ˇZnidariˇc, A. Scardicchio, and V. K. Varma, Diffusive and subdiffusive spin transport in the ergodic phase of a many-bodylocalizable system, Phys. Rev. Lett. , 040601 (2016).[24] D. J. Luitz, N. Laflorencie, and F. Alet, Extended slow dynamical regime close to the many-body localization transition,Phys. Rev. B , 060201 (2016).[25] D. J. Luitz and Y. Bar Lev, Anomalous thermalization in ergodic systems, Phys. Rev. Lett. , 170404 (2016).[26] S. Bera, G. De Tomasi, F. Weiner, and F. Evers, Density propagator for many-body localization: Finite-size effects,transient subdiffusion, and exponential decay, Phys. Rev. Lett. , 196801 (2017).[27] M. Serbyn, Z. Papi´c, and D. A. Abanin, Thouless energy and multifractality across the many-body localization transition,Phys. Rev. B , 104201 (2017).[28] A. Dymarsky, Bound on eigenstate thermalization from transport, (2018), arXiv:1804.08626 [cond-mat.stat-mech].[29] D. J. Luitz and Y. Bar Lev, The ergodic side of the many-body localization transition, Ann. Phys. , 1600350 (2017).[30] B. Bertini, F. Heidrich-Meisner, C. Karrasch, T. Prosen, R. Steinigeweg, and M. Znidaric, Finite-temperature transportin one-dimensional quantum lattice models, (2020), arXiv:2003.03334 [cond-mat.str-el].[31] S. Gopalakrishnan and S. Parameswaran, Dynamics and transport at the threshold of many-body localization, Phys. Rep. , 1 (2020).[32] V. Zelevinsky, B. A. Brown, N. Frazier, and M. Horoi, The nuclear shell model as a testing ground for many-body quantumchaos, Phys. Rep. , 85 (1996).[33] F. Borgonovi, F. Izrailev, L. Santos, and V. Zelevinsky, Quantum chaos and thermalization in isolated systems of interactingparticles, Phys. Rep. , 1 (2016).[34] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, From quantum chaos and eigenstate thermalization to statisticalmechanics and thermodynamics, Adv. Phys. , 239 (2016).[35] M. Schiulaz, E. J. Torres-Herrera, and L. F. Santos, Thouless and relaxation time scales in many-body quantum systems,Phys. Rev. B , 174313 (2019).[36] L. Leviandier, M. Lombardi, R. Jost, and J. P. Pique, Fourier transform: A tool to measure statistical level properties invery complex spectra, Phys. Rev. Lett. , 2449 (1986).[37] T. Guhr and H. Weidenmuller, Correlations in anticrossing spectra and scattering theory. analytical aspects, Chem. Phys. , 21 (1990).[38] J. Wilkie and P. Brumer, Time-dependent manifestations of quantum chaos, Phys. Rev. Lett. , 1185 (1991).[39] Y. Alhassid and R. D. Levine, Spectral autocorrelation function in the statistical theory of energy levels, Phys. Rev. A ,4650 (1992).[40] T. Gorin and T. H. Seligman, Signatures of the correlation hole in total and partial cross sections, Phys. Rev. E , 026214(2002).[41] E. Torres-Herrera and L. F. Santos, Dynamical manifestations of quantum chaos: correlation hole and bulge, Philos. Trans.R. Soc. London, Ser. A , 20160434 (2017). [42] E. J. Torres-Herrera, A. M. Garc´ıa-Garc´ıa, and L. F. Santos, Generic dynamical features of quenched interacting quantumsystems: Survival probability, density imbalance, and out-of-time-ordered correlator, Phys. Rev. B , 060303 (2018).[43] J. Cotler, N. Hunter-Jones, J. Liu, and B. Yoshida, Chaos, complexity, and random matrices, J. High Energy Phys. (11), 48.[44] H. Gharibyan, M. Hanada, S. H. Shenker, and M. Tezuka, Onset of random matrix behavior in scrambling systems, J.High Energy Phys. (7), 124.[45] M. Winer and B. Swingle, Hydrodynamic theory of the connected spectral form factor (2020), arXiv:2012.01436.[46] S. Gherardini, C. Lovecchio, M. M. M¨uller, P. Lombardi, F. Caruso, and F. S. Cataliotti, Ergodicity in randomly perturbedquantum systems, Quantum Sci. Technol. , 015007 (2017).[47] K. Singh, C. J. Fujiwara, Z. A. Geiger, E. Q. Simmons, M. Lipatov, A. Cao, P. Dotti, S. V. Rajagopal, R. Senaratne,T. Shimasaki, M. Heyl, A. Eckardt, and D. M. Weld, Quantifying and controlling prethermal nonergodicity in interactingfloquet matter, Phys. Rev. X , 041021 (2019).[48] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. L¨uschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch,Observation of many-body localization of interacting fermions in a quasirandom optical lattice, Science , 842 (2015).[49] P. Richerme, Z.-X. Gong, A. Lee, C. Senko, J. Smith, M. Foss-Feig, S. Michalakis, A. V. Gorshkov, and C. Monroe,Non-local propagation of correlations in quantum systems with long-range interactions, Nature , 198 (2014).[50] L. F. Santos, G. Rigolin, and C. O. Escobar, Entanglement versus chaos in disordered spin chains, Phys. Rev. A , 042304(2004).[51] R. Nandkishore and D. A. Huse, Many-body localization and thermalization in quantum statistical mechanics, Annu. Rev.Condens. Matter Phys. , 15 (2015).[52] F. Alet and N. Laflorencie, Many-body localization: An introduction and selected topics, C. R. Phys. , 498 (2018).[53] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Colloquium: Many-body localization, thermalization, and entanglement,Rev. Mod. Phys. , 021001 (2019).[54] V. Oganesyan and D. A. Huse, Localization of interacting fermions at high temperature, Phys. Rev. B , 155111 (2007).[55] A. Pal and D. A. Huse, Many-body localization phase transition, Phys. Rev. B , 174411 (2010).[56] T. C. Berkelbach and D. R. Reichman, Conductivity of disordered quantum lattice models at infinite temperature: Many-body localization, Phys. Rev. B , 224429 (2010).[57] J. A. Kj¨all, J. H. Bardarson, and F. Pollmann, Many-Body Localization in a Disordered Quantum Ising Chain, Phys. Rev.Lett. , 107204 (2014).[58] Y. Bar Lev and D. R. Reichman, Dynamics of many-body localization, Phys. Rev. B , 220201 (2014).[59] D. J. Luitz, N. Laflorencie, and F. Alet, Many-body localization edge in the random-field Heisenberg chain, Phys. Rev. B , 081103(R) (2015).[60] T. Devakul and R. R. P. Singh, Early breakdown of area-law entanglement at the many-body delocalization transition,Phys. Rev. Lett. , 187201 (2015).[61] E. V. H. Doggen, F. Schindler, K. S. Tikhonov, A. D. Mirlin, T. Neupert, D. G. Polyakov, and I. V. Gornyi, Many-bodylocalization and delocalization in large quantum chains, Phys. Rev. B , 174202 (2018).[62] L. Herviou, S. Bera, and J. H. Bardarson, Multiscale entanglement clusters at the many-body localization phase transition,Phys. Rev. B , 134205 (2019).[63] D. Abanin, J. Bardarson, G. De Tomasi, S. Gopalakrishnan, V. Khemani, S. Parameswaran, F. Pollmann, A. Potter,M. Serbyn, and R. Vasseur, Distinguishing localization from chaos: Challenges in finite-size systems, Annals of Physics ,168415 (2021).[64] M. L. Mehta, Random Matrices (Academic Press, Boston, USA, 1991).[65] J. ˇSuntajs, J. Bonˇca, T. Prosen, and L. Vidmar, Quantum chaos challenges many-body localization (2020),arXiv:1905.06345 [cond-mat.str-el].[66] E. J. Torres-Herrera and L. F. Santos, Extended nonergodic states in disordered many-body quantum systems, Annalender Physik , 1600284 (2017).[67] E. J. Torres-Herrera and L. F. Santos, Signatures of chaos and thermalization in the dynamics of many-body quantumsystems, Eur. Phys. J. Spec. Top. , 1897 (2019).[68] Since for W = 0 the system is integrable, the disorder cannot be too small.[69] M. T´avora, E. J. Torres-Herrera, and L. F. Santos, Inevitable power-law behavior of isolated many-body quantum systemsand how it anticipates thermalization, Phys. Rev. A , 041603 R (2016).[70] M. T´avora, E. J. Torres-Herrera, and L. F. Santos, Power-law decay exponents: A dynamical criterion for predictingthermalization, Phys. Rev. A , 013604 (2017).[71] E. J. Torres-Herrera, I. Vallejo-Fabila, A. J. Mart´ınez-Mendoza, and L. F. Santos, Self-averaging in many-body quantumsystems out of equilibrium: Time dependence of distributions, Phys. Rev. E , 062126 (2020).[72] The statistical errors of this procedure were accounted for using a bootstrapping procedure with 1000 bootstrap samples,each composed of N = 9000 random samples taken from the full random data set of over 10 elements.[73] J. Richter, A. Dymarsky, R. Steinigeweg, and J. Gemmer, Eigenstate thermalization hypothesis beyond standard indicators:Emergence of random-matrix behavior at small frequencies, Phys. Rev. E102