# Late-time universal distribution functions of observables in one-dimensional
many-body quantum systems

I. Vallejo-Fabila  E. Jonathan Torres-Herrera  Instituto de Física, Benemérita
Universidad Autónoma de Puebla, Puebla, 72570, México

###### Abstract

We study the probability distribution function of the long-time values of
observables being time-evolved by Hamiltonians modeling clean and disordered
one-dimensional chains of many spin-1/2 particles. In particular, we analyze
the return probability and its version for a completely extended initial
state, the so-called spectral form factor. We complement our analysis with the
spin autocorrelation and connected spin-spin correlation functions, both of
interest in experiments with quantum simulators. We show that the distribution
function has a universal shape provided the central limit theorem holds.
Explicitly, the shape is exponential for the return probability and spectral
form factor, meanwhile it is Gaussian for the few-body observables. We also
discuss implications over the so-called many-body localization. Remarkably,
our approach requires only a single sample of the dynamics and small system
sizes, which could be quite advantageous when dealing specially with
disordered systems.

##  I Introduction

Probability distribution functions (pdfs) are of fundamental interest in
different branches of science. In physics, in particular, pdfs are important
not only for the purpose of computing mean values or expectation values of
relevant physical quantities but also by themselves because they could provide
fine information about fluctuations around the mean values of properties of a
given system. For a comprehensive review on fluctuations in physical systems
see, for instance, Ref. [1].

In the realm of quantum mechanics, the characterization of fluctuations in
time domain has been useful to understand reversibility in generic chaotic and
integrable systems [2], to establish a Gaussian scenario for equilibration in
noninteracting models [3, 4], where fluctuations can decay as the square root
of the system size. Although later is was shown in Ref. [5] that for many-body
systems the fluctuations of generic observables decay exponentially with
system size in both, integrable and chaotic systems.

Late-time pdfs are also important in the context of thermalization and many-
body localization (MBL). Just some examples are: Late-time exponential
probability distribution functions experimentally observed for the squared-
contrast of a one-dimensional Bose gas and related to prethermalization and
thermal equilibrium in Ref. [6]. Gaussian momentum distribution functions for
noninteracting spinless fermions and hard-core bosons in one-dimensional (1D)
lattices were observed in Refs. [7, 8]. Gaussian distributions are
experimentally obtained for the equilibrium value of a temporal
autocorrelation function, number entropy and Hamming distance in Ref. [9]. In
numerical experiments employing the time dependence variational principle for
matrix product states and machine learning, late-time pdfs where analyzed for
the spin imbalance, entanglement entropy and so-called Schmidt gap [10]. Pdfs
were useful for the determination of the universality class of MBL under
speckle disorder [11]. Long-tails and bimodal probability distribution
functions of the entanglement entropy near and at the MBL transition,
respectively, were previously analyzed in Refs. [12, 13].

From our point of view a clear determination of the conditions under which the
formerly found, experimentally and theoretically, pdfs of observables in
quantum systems is missing. Our aim in this work is twofold, one is to
generalize the Gaussian scenario established in Refs. [3, 4] and the other one
is to set, in a systematic way, the basic conditions that lead to universal
late-time probability distribution functions of physical quantities and
observables in 1D interacting quantum systems. We will show that the central
limit theorem (CLT) is behind the universal late-time pdfs, and will also
analyze the conditions under which this almost ubiquitous and undoubtedly
celebrated theorem holds in the time evolutions generated by paradigmatic
models of spin-1/2 particles, with or without interactions, clean or
disordered, chaotic or integrable. In particular, we conjecture that a
sufficient condition to have universal probability distribution functions is
that of linearly uncorrelated identically distributed (u.i.d.) variables.
Projections of the initial state in the energy eigenstates, as well as
eigenvalues are the main characters playing a relevant role in the fulfillment
of the CLT. Of course, when dealing with observables, matrix elements in
energy representation could be also relevant for the fulfillment of the CLT.

Our analysis will employ different probes of the dynamics, namely, return
probability, spectral form factor, spin autocorrelation function and connected
spin-spin correlation function. We will put main emphasis in the first one,
the return probability, this because it has been prominently studied since the
early years of quantum mechanics. Let us give another nonexhaustive but larger
list of examples: In the context of the time-energy uncertainty relation [14,
15], nonexponential late-time decays of a quasistationary state [16], quantum
Zeno’s paradox [17], quantum speed limits [18, 19, 20, 21], quantum energy
flow [22], nonperiodic substitution potentials [23], connection with the time
operator through a generalized Weyl relation [24], cosmology [25], quantum
walks and complex networks [26, 27, 28, 29], dynamics of a thermofield double
state [30, 31], and matter-radiation interaction models like Dicke [32] and
Bose-Hubbard [33]. The return probability has been also measured in
experiments to observe time-resolved level repulsion of chaotic systems [34],
fluorescence experiments with molecules [35, 36], ultracold atoms in magneto-
optic traps [37], atom chip [38], and prethermalization in Floquet systems
[39].

The rest of the manuscript is organized as follows. Models and quantities are
introduced in Secs. II and III, respectively. In Sec. IV we provide an
overview of the dynamics from the time it is initiated up to equilibrium. In
Sec. V we present a justification of the universal distributions based on
theoretical grounds. Results for the addressed systems are shown and discussed
in Sec. VI. Conclusions are finally presented in Sec. VII.

##  II Models

We consider a system of spin-$1/2$ particles that can interact with nearest-
neighbors in a 1D lattice. Particles can also be subjected to an on-site
random potential. The general Hamiltonian describing such a system is

$\mathcal{H}=\sum_{k}(S_{k}^{x}S_{k+1}^{x}+S_{k}^{y}S_{k+1}^{y}+\Delta
S_{k}^{z}S_{k+1}^{z})+\sum_{k=1}^{L}h_{k}S_{k}^{z}.$ (1)

With spin-$1/2$ operators $S_{k}^{x,y,z}$ acting on the spin located at the
$k$-th site. We set $\hbar=1$. The anisotropy parameter $\Delta$ and the on-
site potential amplitudes $h_{k}$ are tuned to obtain different instances of
model (1). Setting $\Delta=0$ together with $h_{k}=0$ leads to a
noninteracting model, known as $XX$ model which is solved by applying a
Jordan-Wigner transformation [40, 41]. With $\Delta=1$ and $h_{k}=0$ we obtain
the isotropic $XXZ$ model with Ising-like interactions in $z$ direction, this
model in also solvable but via the celebrated Bethe ansatz [42, 43, 44, 45,
46]. We arrive to a paradigmatic model for the MBL transition by setting again
$\Delta=1$ but with $h_{k}$ as uncorrelated random numbers uniformly
distributed in the interval $(-h,h)$, with $h$ as the disorder strength. A
critical point $h_{c}$ has been identified, $3.75\lesssim h_{c}$, but still
without consensus about its precise location in the $h$-axis, see for instance
Refs. [47, 48, 49, 50, 51, 52, 10, 53]. For the $XX$ and $XXZ$ models we
consider open boundaries conditions, meanwhile for the disordered model we use
periodic boundary conditions, $\hat{S}_{L+1}=\hat{S}_{1}$, implying that in
the first summation of Hamiltonian (1) we have the index $k$ running from $1$
to $L-1$ for the two former models and from $1$ to $L$ for the latter.

Hamiltonian (1) conserves the total spin in the $z$-direction,
${\cal{S}}^{z}=\sum_{k}S_{k}^{z}$, it only moves an excitation (spin up)
through the chain due to the $XX$ term which is also known as the flip-flop
term. This allows to focus our study on the largest subspace with
${\cal{S}}^{z}=0$ and dimension ${\cal{N}}=L!/(L/2)!^{2}$.

##  III Quantities and observables

In this section we describe the time-dependent quantities and observables
employed in our study. We also present the level spacing distribution
traditionally used as a probe of repulsion between adjacent energy levels, a
fingerprint of quantum chaos.

###  III.1 Return probability and spectral form factor

The return probability, RP, also known as survival or nondecay probability is
a dynamical quantity defined as the probability to find an initial state
$\ket{\Psi(0)}$ at a future time $t$. It is given by

$\text{RP}(t)=|S(t)|^{2}=\left|\displaystyle{\sum_{\alpha=1}^{\cal{N}}|c_{\alpha}^{0}|^{2}e^{-iE_{\alpha}t}}\right|^{2},$
(2)

where $S(t)=\bra{\Psi(0)}{\cal{U}}(t)\ket{\Psi(0)}$ is the return amplitude
and ${\cal{U}}(t)=\exp(-i{\mathcal{H}}t)$ the unitary time evolution operator.
For the Hamiltonian generating the dynamics we have the eigenvalue equation
$\mathcal{H}\ket{\psi_{\alpha}}=E_{\alpha}\ket{\psi_{\alpha}}$,
$\alpha=1,\,2,\,\dots,\,{\cal{N}}$. Meanwhile
$c_{\alpha}^{0}=\braket{\psi_{\alpha}}{\Psi(0)}$ are the projections of the
initial state into the energy eigenstates.

When the initial state $\ket{\Psi(0)}$ is completely extended in the energy
eigenbasis, its components goes as $1/\sqrt{\cal{N}}$ and from Eq. (2) we
recover the so-called spectral form factor,

$K(t)=\frac{1}{{\cal{N}}^{2}}\left|\sum_{\alpha=1}^{\cal{N}}e^{-iE_{\alpha}t}\right|^{2},$
(3)

where only energy eigenvalues are involved.

###  III.2 Spin autocorrelation function

How close the spin configuration at time $t$ is to the one at $t=0$ can be
measured through so-called spin autocorrelation function, given by

$I(t)=\frac{4}{L}\sum_{i=1}^{L}\braket{\Psi(0)}{\hat{S}_{i}^{z}e^{i\hat{\mathcal{H}}t}\hat{S}_{i}^{z}e^{-i\hat{\mathcal{H}}t}}{\Psi(0)}.$
(4)

For particular initial states, like a quantum Néel-like state, $I(t)$ is
equivalent to the density imbalance measured in experimental platforms
studying ultracold atoms [54, 55].

###  III.3 Connected spin-spin correlation function

Another quantity also of interest in quantum simulators is the connected spin-
spin correlation function, defined through

$\begin{array}[]{ll}C(t)&=\displaystyle{\frac{4}{L}\sum_{k=1}^{L-1}\left[\bra{\Psi(t)}S_{k}^{z}S_{k+1}^{z}\ket{\Psi(t)}\right.}\\\
&\displaystyle{\left.\hskip
28.45274pt-\bra{\Psi(t)}S_{k}^{z}\ket{\Psi(t)}\bra{\Psi(t)}S_{k+1}^{z}\ket{\Psi(t)}\right]}.\end{array}$
(5)

It quantifies time-dependent correlations between neighboring spins in the
chain. It has been measured, for instance, in experiments with trapped ions
studying systems with long-range interactions [56].

###  III.4 Level spacings

We present in this subsection the theoretical predictions for the distribution
$P(s)$ of spacings between adjacent energy levels,
$s_{\alpha}=(E_{\alpha+1}-E_{\alpha})/\delta E$, with $\delta E$ being the
mean level spacing of the system. For integrable systems the level spacings
typically follow an exponential distribution, they obey Poisson-like
statistics [57],

$P_{\text{P}}(s)=\exp(-s).$ (6)

No level repulsion is apparent from Eq. (6) since $P(s\to 0)\to 1$. Meanwhile
spacings between energy levels of time-reversal invariant systems, just like
the ones we deal with, that are chaotic follow statistics from the Gaussian
orthogonal ensemble (GOE) of random matrix theory (RMT) [58],

$P_{\text{GOE}}(s)=\frac{\pi}{2}s\exp\left(-\dfrac{\pi}{4}s^{2}\right),$ (7)

and level repulsion is manifested as indicated by the limit case $P(s\to 0)\to
s$. We stress that $P(s)$ will be used just to signal the regime of the
Hamiltonian, being integrable or chaotic (ergodic, thermalizing), and also to
expose possible degeneracies of the energy spectrum, like in the $XX$ model.

##  IV Dynamics overview

We focus on the late-time behavior of the dynamical quantities and observables
introduced lines above, that is, beyond the so-called Heisenberg time,
$t_{\text{H}}=2\pi/\delta E$. However, it is instructive to illustrate their
time evolution from $t=0$ up to saturation. Figure 1 depicts the full-time
evolution generated by the disordered model with fixed number of sites $L=16$
for $\text{RP}(t)$ [Fig. 1(a)], $K(t)$ [Fig. 1(b)], $I(t)$ [Fig. 1(c)], and
$C(t)$ [Fig. 1(d)]. Three different disorder strengths are considered, $h=0.5$
(red), $h=3.75$ (yellow), and $h=6.0$ (blue). Dark-colored curves are averages
over samples obtained from $130$ disorder realizations and $78$ initial states
with energy closest to zero,
$E_{0}=\bra{\Psi(0)}{\mathcal{H}}\ket{\Psi(0)}\approx 0$. Then a total of
$\sim 10^{3}$ samples were considered for the average. Light-colored curves
are results from five individual samples, i.e., each one representing the
dynamics for a single initial state with energy closest to zero and only one
disorder realization. Note that individual samples for $h=3.75$ are omitted,
this is done only to avoid overcrowding of figure. For $K(t)$ only disorder
realizations were considered. The individual samples present huge sample-to-
sample fluctuations [59] that could hide or obscure some relevant aspects of
the dynamics, as the ones shown by the averaged quantities [60, 61, 62].

Figure 1: Full-time dynamics of the quantities and observables treated in our
study. Time evolution is generated by the disordered model in the ergodic
phase with $h=0.5$ (red), $h=3.75$ (orange), and $h=6.0$ (blue). Return
probability (a), spectral form factor (b), spin autocorrelation function (c),
and connected spin-spin correlation function (d). Light-colored curves are
results from five samples of the dynamics, meanwhile dark-colored curves
correspond to an average over $\sim 10^{4}$ samples originated from different
initial states and random realizations. System size is $L=16$.

In fact, from the averages shown in Fig. 1 a clear power-law decay is seen for
the return probability and apparently also for the imbalance [63, 64],
although for this last quantity dissipative dynamics includes a stretched
exponential decay [65, 66]. The power-law decay is followed by the correlation
hole [67, 68, 69] which is visible for both quantities but also for the
spectral form factor. For $C(t)$ the correlation hole is not apparent but it
is present, a deep close up is needed to observe it [70]. We note that also in
Ref. [70] it was shown how these features depend on both, system size and kind
of observable.

Figure 2 shows the whole dynamics for the clean $XXZ$ and $XX$ models, but
considering only the return probability [Figs. 2(a) and 2(c)] and spectral
form factor [Figs. 2(b) and 2(d)]. In Figs. 2(a) and 2(c) the dark-blue curves
are averages over $200$ initial states together with an additional moving
average to reduce fluctuations. Light-blue curves are dynamics for five
individual initial states with energies closest to zero. In Figs. 2(b) and
2(d) the dark curves are moving averages over a single sample. While the
light-colored curves are the rough data. Although interesting transient
behavior is observed in Figs. 1 and 2, it is worth to emphasize that our main
analysis of the distribution functions will be carried out for times beyond
the Heisenberg time, $t_{\text{H}}$. However, we eventually will provide
examples of distribution functions for different timescales. The vertical
dashed line in all panels of both figures marks the point where
$t/t_{\text{H}}=1$.

Figure 2: Full-time dynamics of $\text{RP}(t)$ and $K(t)$ generated by the
$XXZ$ (a, b) and $XX$ (c, d) models. The light-colored curves in panels (a, c)
are for individual samples of $\text{RP}(t)$ obtained for five different
initial states, meanwhile the dark-colored ones correspond to an average over
$200$ initial states, together with an additional moving average performed to
further smooth the curves. For $K(t)$ in panels (b, d) the light-colored
curves represent the dynamics from a single initial state and the dark-colored
curves correspond to moving averages over the same sample. System size is
$L=16$.

##  V Central Limit Theorem and distribution functions

Equation (2) for the return probability can be rewritten as

$\text{RP}(t)=\left|\sum_{\alpha=1}^{\cal{N}}X_{\alpha}(t)-i\sum_{\alpha=1}^{\cal{N}}Y_{\alpha}(t)\right|^{2},$
(8)

where $X_{\alpha}(t)=|c^{0}_{\alpha}|^{2}\cos(E_{\alpha}t)$ and
$Y_{\alpha}(t)=|c^{0}_{\alpha}|^{2}\sin(E_{\alpha}t)$. At this point we remind
that the classical central limit theorem roughly states that if
$X_{\alpha}(t)$ or $Y_{\alpha}(t)$ are independent and identically distributed
(i.i.d.) random variables with finite second moment, then the distribution of
their sum is well approximated by a Gaussian shape. Thus, if we assume that
$X_{\alpha}(t)$ and $Y_{\alpha}(t)$ behave as i.i.d. in a time window
$[t_{i},t_{f}]$ with initial time $t_{i}$ and final time $t_{f}$, then Eq. (8)
can be considered as the squared norm of a complex random variable with
Gaussian real and imaginary parts. In such a case the distribution of
$\text{RP}(t)$ should be exponential [71],

$P(\text{RP})=\lambda e^{-\lambda\text{RP}},\hskip 14.22636pt\text{RP}\geq
0,\hskip 14.22636pt\lambda>0.$ (9)

Considering the terms appearing in Eq. (8), it is clear that if both, initial
state components $c^{0}_{\alpha}$ and energy eigenvalues $E_{\alpha}$ behave
simultaneously as i.i.d., then the CLT should hold. It is important to mention
that in this work (see Appendix A), we show numerical evidence that a less
restrictive condition for the fulfillment of the CLT could be that of
uncorrelated random variables. Assuming that $X_{\alpha}(t)$ and
$Y_{\alpha}(t)$ remain each one identically distributed, three possible
scenarios where the CLT does not hold are when initial state components are
autocorrelated, energy eigenvalues are autocorrelated or both of them are
correlated, each of these cases mean that $X_{\alpha}(t)$ and $Y_{\alpha}(t)$
are not independent. Note, however, that even for weakly enough correlated
variables the CLT could remain true [72].

In the case of $K(t)$ the exponential distribution should be achieved when
correlations between energy eigenvalues are absent, just see the expression
defining $K(t)$ in Eq. (3). The scenario for the expectation value of generic
observables is more involved. From its time evolution

$\left\langle{\cal{O}}(t)\right\rangle=\sum_{\alpha,\beta}c_{\alpha}^{0}c_{\beta}^{0}{\cal{O}}_{\alpha\beta}e^{-i(E_{\alpha}-E_{\beta})t},$
(10)

with
${\cal{O}}_{\alpha\beta}=\left\langle\psi_{\alpha}|{\cal{O}}|\psi_{\beta}\right\rangle$,
we see that matrix elements in energy representation of the observable are
also included. The simplest case arises with the assumption that all terms
inside the summation in Eq. (10) behave as i.i.d. random variables, then by
the CLT we should expect the distribution to be Gaussian. We note that having
a Gaussian distribution for the late-time values of a generic observable could
be an indication that the matrix elements of the observable in energy
representation are u.i.d. by themselves. However, if the distribution is not
Gaussian, then we cannot ensure that the matrix elements are correlated, this
because the presence of energy eigenvalues and eigenstates inside the
summation that could also be correlated. Although interesting, we leave the
analysis of those fine details for a future work.

At this point it is fair to mention that Aurich and Steiner conjetured in Ref.
[73] that given an extended initial state, the distribution function of a
normalized version of the return amplitude $S(t)$ for a chaotic quantum system
should be universally described by Rayleigh’s law,

$P(S)=\frac{\pi}{2}S\exp\left(-\frac{\pi}{4}S^{2}\right),\hskip
28.45274ptS\geq 0.$ (11)

This conjecture was established in the context of 2D and 3D quantum billiards.
It is straightforward to show that the square root of an exponential
distribution, like the one for the return probability, results in a Rayleigh
distribution just as Eq. (11), this is consistent with Eq. (9) which was
derived using only basics from probability theory. Also fair is to note that
in the context of RMT, Kunz arrived to the rigorous result that the pdf of
both, $K(t)$ and $\text{RP}(t)$ are exponential, independently of the initial
state in the later case [74, 75].

##  VI Results

Now we present and describe our results on the pdfs, starting with the return
probability $\text{RP}(t)$, to which we devote majority of this paper. Results
about level statistics are also presented, but of course they are no new in
the existing literature. As explained before, we use them only as a reference.

###  VI.1 Return probability and spectral form factor

Figure 3 depicts $P(s)$ and $P(\text{RP})$ for the disordered model and fixed
number of spins $L=18$. The upper panels $P(s)$ show the behavior already
known in the vicinity of the MBL transition [76]. For small disorder strength,
$h=0.5$ [Fig. 3 (a)], $P(s)$ has the shape described by Eq. (7). At
intermediate disorder strength, $h=3.75$, the shape of $P(s)$ is neither GOE-
like nor Poisson-like, although closer to the Poisson shape than to GOE. Once
at the localized phase, $h=6.0$, $P(s)$ gets the shape predicted by Eq. (6).

Figure 3: Level spacing distribution $P(s)$ (upper panels) and distribution of
the late-time values of $\text{RP}(t)$ (lower panels) for the disordered
model. The time intervals are $t/t_{\text{H}}\in(11.38,14.93)$,
$t/t_{\text{H}}\in(29.25,38.38)$, and $t/t_{\text{H}}\in(33.34,43.74)$ for the
disorder strengths $h=0.5$ (a), $h=3.75$ (b) and $h=6.0$ (c), respectively.
Solid lines in the upper panels depict the theoretical values for $P(s)$, Eq.
(6) [color red in panels (b) and (c)] and Eq. (7) [color blue in panels (a)
and (b)]. Solid curves in the lower panels correspond to Eq. (9) with mean 1.
System size is $L=18$.

Distributions of $\text{RP}(t)$ with $t/t_{\text{H}}>1$ are shown in the lower
panels of Fig. 3, specific time intervals are indicated in the captions. An
initial state with energy $E_{0}$ closest to zero is chosen and a single
disorder realization is considered. Note that in the whole work we have
normalized the values of $\text{RP}(t)$ with their standard deviation in the
respective time interval. Figure 3 (d) confirms the anticipated in Sec. V, for
$h=0.5$ the distribution of RP (histogram in color blue) conforms with the
exponential distribution given by Eq. (9) with $\lambda=1$. Then we could
assume that $c_{\alpha}^{0}$, $E_{\alpha}$, $X_{\alpha}(t)$ and
$Y_{\alpha}(t)$ behave as uncorrelated random variables and therefore the
universal exponential pdf is achieved (see discussion in Appendix A, where
evidence is given in the sense that the assumption is actually a fact).The
increase of disorder strength $h$ is accompanied by the lack of correspondence
between the exponential and the distributions of RP, which is clearly seen in
Figs. 3(e) and 3(f) for $h=3.75$ and $h=6.0$, respectively. This is a relevant
observation because the distribution $P(\text{RP})$ could detect the MBL
transition, and remarkably from a single sample of the dynamics. However, we
emphasize that detection of the MBL transition is not the one of the aims of
this work, which instead is the establishment of a Gaussian scenario for the
equilibrium state of interacting quantum systems.

The shape of $P(\text{RP})$ different from the exponential is a signature of
autocorrelations, be in the energy eigenvalues or between initial state
components. For this last statement, we can apply the negation of the CLT, if
the distribution of the sums involved in Eq. (8) does not tend to a Gaussian,
then the random variables are not statistical independent, they are
correlated. Correlations between eigenstates components have been shown to
exist around the MBL transition; see, for instance, Refs. [63, 77] and more
recently Ref. [78]. In fact, the slow dynamics produced by those correlations
have been observed experimentally [79].

Figure 4: Distribution of the late-time values of $K(t)$ for the same time
window as in Fig. 3 for the disordered model with $h=0.5$ (a), $h=3.75$ (b),
and $h=6.0$ (c). Solid curves are given by $P(K)=\exp(-K)$. System size is
$L=18$.

However, up to our knowledge, information about autocorrelations of energy
levels does not exist. To fill out this gap, we recur to the spectral form
factor and analyze its distribution for $t/t_{\text{H}}>1$. Figure 4 is
conclusive, it shows that for any disorder strength in the MBL model the
distribution of $K(t)$ has an exponential shape. This exponential shape can be
explained analogously to the case of $\text{RP}(t)$, that is, in terms of the
CLT. As suggested by the results depicted Fig. A.1.1 (a) of Appendix A, the
energy eigenvalues remain uncorrelated for any disorder strength. Then the
sole cause for the unconform shapes of $P(\text{RP})$ in Figs. 3(e) and 3(f)
should be the autocorrelations between initial state components
$c_{\alpha}^{0}$. To reinforce our statement about autocorrelations between
initial state components, we present in Fig. 5 the distribution $P(\text{RP})$
departing from a random initial state with tuned degree of autocorrelations
and evolved by the disordered model in the chaotic region, with fixed $h=0.5$.
There the components $c_{\alpha}^{0}$ of the initial state are random numbers
from a uniform distribution and with autocorrelation degree $q$, generated by
the recipe given in Ref. [80] and recently used to show the effects of
correlated disorder in the region around the MBL transition in Ref. [81].
Figure 5 shows that as the degree of autocorrelations moves from a regime of
strong autocorrelations (a), passing trough intermediate (b) and finally
arriving to weak autocorrelations (c), the exponential distribution of
$P(\text{RP})$ is eventually achieved. The autocorrelations of initial state
components are confirmed through Fig A.1.1 of Appendix A, there it is shown
that while, as already said, energy eigenvalues remain uncorrelated for any
disorder strength, the initial state components get autocorrelated once the
disorder strength is strong enough, say $h=3.75$. Thus apparently inducing
autocorrelations in $X_{\alpha}(t)$ and $Y_{\alpha}(t)$, as suggested also by
the results shown in Figs. A.1.1(c) and A.1.1(d), respectively.

Figure 5: Distribution of the late-time values of $\text{RP}(t)$ under the
disordered model in the ergodic phase with $h=0.5$. The considered time window
is $t/t_{\text{H}}\in(11.38,14.93)$. The initial state is a random one with
correlated components $c_{\alpha}^{0}$. Distributions for different degrees of
correlation are depicted: strong $q=0.5$ (a), intermediate $q=0.33$ (b) and
weak $q=0.2$ (c). Red solid curve in panel (c) corresponds to
$P(\text{RP})=\exp(-\text{RP})$. System size is $L=18$.

####  VI.1.1 XX and XXZ models

As already discussed, autocorrelations between energy levels could block the
fulfillment of the CLT. Those correlations can be caused, for instance, by
degeneracies in the energy spectrum. A suitable model to show this is the $XX$
model [$\Delta=0$ and $h=0$ in Eq. (1)], which energy spectrum contains a high
amount of degeneracies [5]. Figure 6 (a) showing $P(s)$ for the $XX$ model
confirms the degeneracies with the Shnirelman’s peak as a witness [82].
Consequently, the distribution not only of $\text{RP}(t)$ but also of $K(t)$
should not be exponential, which is confirmed by Figs. 6(c) and 6(e) for
$P(\text{RP})$ and $P(K)$, respectively.

We complement our analysis with the study of the time evolution generated by
the $XXZ$ model [$\Delta=1$ and $h=0$ in Eq. (1)]. The Ising interactions
present in the integrable $XXZ$ model break some symmetries of the $XX$ model
and thus removes degeneracies in the energy spectrum. This is reflected into
$P(s)$, Fig. 6(b), which displays Poisson-like statistics. The distribution of
the return probability, $P(\text{RP})$, is exponential meaning that energy
eigenvalues and initial state components behave as uncorrelated random
variables.

Figure 6: Distribution of spacings $P(s)$ (upper panels), distribution of the
late-time values of $\text{RP}(t)$ (middle panels) and $K(t)$ (lower panels).
Left column contains results for the $XX$ model, meanwhile in the right
columns the results are for the $XXZ$ model. The considered time windows are
$t/t_{\text{H}}\in(9.092,11.928)$ for the $XX$ model and
$t/t_{\text{H}}\in(9.859,12.935)$ for the $XXZ$ model. $L=18$.

Figure 6 for $P(K)$ confirms absence of autocorrelations between eigenvalues.
To confirm the last same for initial state components, we employed the so-
called inverse participation ratio,
$\text{IPR}=\sum_{\alpha=1}^{\cal{N}}|\braket{\psi_{\alpha}}{\Psi(0)}|^{4}\propto{\cal{N}}^{-D_{2}}$,
with $0\leq D_{2}\leq 1$ known as the correlation dimension that measures the
degree of correlations between eigenstate components [83]. Two extreme cases
$D_{2}=0$ and $D_{2}=1$ mean correlations and absence of correlations,
respectively. Our scaling analysis of IPR with dimension $\cal{N}$ (not shown)
for the initial state $\ket{\Psi(0)}$ used in Fig. 6 (d) lead us to
$D_{2}\approx 0.75$, this value suggests weak correlations between initial
state components $c_{\alpha}^{0}$, at least weak enough to have the CLT still
valid.

Our results for the integrable $XXZ$ model reaffirm the idea that the
universal late-time exponential distribution is achieved once the conditions
for the CLT fulfillment are provided, that is, simultaneous absence or weak
autocorrelations in the initial state components and in the energy spectrum of
the Hamiltonian generating the time evolution.

####  VI.1.2 Timescales

Up to now, the analysis has been carried out using timescales greater than the
Heisenberg time $t_{\text{H}}$, an interesting question is if the universal
distribution could appear for smaller timescales. Figure 7 displays the
distribution of RP for the disordered model in the ergodic phase with $h=0.5$
and for different time intervals, at very short times [Fig. 7(a)], power-law
decay [Fig. 7(b)], around the correlation hole [Fig. 7(c)], and saturation
[Fig. 7(d)]. The reader is referred to Fig. 3(a) and discussion in Sec. IV to
remind all those timescales.

Figure 7: Distribution $P(\text{RP})$ of the return probability for the
disordered model in the ergodic phase with $h=0.5$ in different time windows.
Initial decay, $t/t_{\text{H}}\in(0.000034,0.000089)$ (a), power-law decay,
$t/t_{\text{H}}\in(0.00034,0.00089)$ (b), around the time where the minimum of
the correlation hole (Thouless time) is located,
$t/t_{\text{H}}\in(0.0091,0.0228)$ (c), and saturation,
$t/t_{\text{H}}\in(11.38,14.93)$ (d). Solid curves are given by
$P(\text{RP})=\exp(-\text{RP})$. System size is $L=18$.

Once the initial state begins the transition to other basis states, as time
increases, the distribution $P(\text{RP})$ moves gradually to an exponential
one. At short times the return probability is very similar in values and a
peaked distribution is obtained, Fig. 7(a). The exponential shape is
completely achieved when the time is large enough like in Fig. 7(d); however,
note that for times when the correlation hole appears (Thouless time) [Fig.
7(c)] the distribution of RP also looks like an exponential but not quite good
as in Fig. 7(d). Then we infer that to observe the universal exponential
distribution one should wait for times beyond the Heisenberg time,
$t_{\text{H}}$. Now we have the opportunity to mention that the exponential
distribution of $\text{RP}(t)$ in a time interval coincides with the
exponential observed in Ref. [62] also at large times, but there fixing a
single time $t_{0}$ and computing the pdf over different disorder realizations
and initial states. This is consistent with the fact that to study ergodicity
one should do it once at equilibrium.

###  VI.2 Spin autocorrelation function and connected spin-spin correlation
function

In this section we focus on the analysis of the pdf of two dynamical
observables of great interest in quantum simulations with cold atoms and ion
traps, the spin autocorrelation function $I(t)$ and the connected spin-spin
correlation function $C(t)$ given by Eqs. (5) and (4), respectively. We see in
the upper panels of Fig. 8 the pdf for the spin autocorrelation function,
$P[I(t/t_{\text{H}}>1)]$, meanwhile in the lower panels of Fig. 8 we have the
pdf for the connected spin-spin correlation function,
$P[C(t/t_{\text{H}}>1)]$. Both quantities are evolved by the disordered model
with different disorder strengths. As predicted in Sec. V, when $h=0.5$ the
late-time pdf are Gaussian for both observables [Figs. 8(a) and 8(d)].

Figure 8: Late-time probability distribution functions corresponding to the
spin autocorrelation function (upper panels) and the connected spin-spin
correlation function (lower panels) evolved under the disordered model. The
considered time windows are the same as in Fig. 3. Disorder strengths are
$h=0.5$ (a, d), $h=3.75$ (b, e), and $h=6.0$ (c, f). Solid curves are Gaussian
functions, $P(x)=(\sqrt{2\pi})^{-1/2}e^{-(x-\mu)^{2}/2}$, with $\mu$ the mean
of the corresponding observable. System size is $L=18$.

Once the system is outside the ergodic phase, the pdf for both observables do
not conform with the Gaussian shape. Comparing Figs. 8(b) and 8(e), for
$h=3.75$, we see that $P(I)$ is more sensitive to the change in disorder than
$P(C)$. Finally, when the disorder strength is $h=6.0$, the lack of
correspondence between the pdf of both observables and the Gaussian is even
more evident. Again, the shape of the pdf could be used as a probe of the MBL
transition. Interestingly, our results in the ergodic side of the MBL are in
line with results for noninteracting systems in Refs. [3, 4].

##  VII Conclusions

We analyzed the late-time probability distribution functions of dynamical
quantities of interest for experimental platforms. We argued that the pdfs are
universal in the sense that the unique requirement is the CLT to hold. This
needs uncorrelated initial state components and energy spectrum. For chaotic
quantum systems these conditions are fulfilled, but also for integrable
systems like the $XXZ$ model. Strong degeneracies in the energy spectrum, like
in the $XX$ model, prevent the CLT to hold, as well as correlated initial
state components, like in the disordered model for intermediate and strong
disorder strengths. For generic observables, an additional condition should be
fulfilled, uncorrelated elements in energy representation. This happens for
the two quantities studied in this work, the spin autocorrelation function and
the connected spin-spin correlation function. No averages are needed, our
analysis was based on a single sample of the dynamics. Details about
autocorrelations between initial state components, energy levels and elements
of observables deserve a deeper analysis, but this is leaved for a future
work. We expect our results to motivate further studies of late-time pdfs in
several quantum systems.

###### Acknowledgements.

I.V.-F. and E.J.T.-H. are grateful to LNS-BUAP for their supercomputing
facility. E.J.T.-H. is also grateful for financial support from VIEP-BUAP,
Project No. 00270-2022 and CONAHCYT under Project Ciencia de Frontera No.
CF-2023-I-1748.

##  Appendix A Autocorrelations

To test the lack or presence of autocorrelations between eigenvalues
$E_{\alpha}$, initial state components $c_{\alpha}^{0}$ and the variables
$X_{\alpha}(t)$ and $Y_{\alpha}(t)$ defined in the context of Eq. (8), we
employ the so-called autocorrelation function, given by

$R_{x}(k)=\frac{1}{({\cal{N}}-k)\sigma^{2}}\sum_{\alpha=1}^{{\cal{N}}-k}(x_{\alpha}-\mu)(x_{\alpha+k}-\mu).$
(12)

This function measures the linear correlation between the sequence of values
$x_{\alpha}$ and a copy of it at lag $k$, that is, $x_{\alpha+k}$. In Eq.
(12), $\mu=\left\langle x_{\alpha}\right\rangle$ and $\sigma^{2}=\left\langle
x_{\alpha}^{2}\right\rangle-\left\langle x_{\alpha}\right\rangle^{2}$ are the
mean and variance, respectively, of $x_{\alpha}$. Having $R_{x}(k)=0$ for all
lag $k$ means absence of autocorrelations. Meanwhile, $R_{x}(k)=1$ for some
$k$ indicates maximum level of autocorrelations, which by definition is always
the case for $k=0$. In the following we avoid $k=0$ and restrict ourselves to
$k=1,\,\dots,\,{\cal{N}}/2$, this because of in our case $k$ larger than
${\cal{N}}/2$ do not provide additional information.

Our results for the autocorrelation function $R_{x}(k)$ are shown in Fig.
A.1.1. We consider the MBL model with $h=0.5$ (red), $h=3.75$ (orange,
squares), and $h=6.0$ (blue, circles). The system size if fixed to $L=18$.
Note that we normalize $k$ with ${\cal{N}}/2$. All panels of Fig. A.1.1
include a pair of horizontal dashed lines delimiting a standard confidence
interval of 95% computed as $\pm 1.96/\sqrt{\cal{N}}$, thus $R_{x}(k)$ staying
inside this interval for all $k$ suggests lack of autocorrelations.

Figure A.1.1(a) confirms what was already anticipated as a naive assumption in
the discussion around Fig. 4 for the late-time distribution function of the
spectral form factor $K(t)$, energy eigenvalues remain uncorrelated for any
value of the disorder strength $h$. This is witnessed by the behavior of
$R_{E_{\alpha}}(k)$, which is inside of the confidence interval for all lag
$k$.

Figure A.1.1(b) for $R_{c_{\alpha}^{0}}(k)$ is quite interesting, for $h=0.5$
its value for any lag $k$ lies inside the confidence interval, suggesting that
the initial state components $c_{\alpha}^{0}$ are not autocorrelated. In
contrast, since $R_{c_{\alpha}^{0}}(k)$ lies notably outside of the confidence
interval for some lags, the image for $h=3.75$ and $h=6.0$ is that of
autocorrelated components. In more detail, we observe that for $h=6.0$ some
values of $R_{c_{\alpha}^{0}}(k)$ are larger than for $h=3.75$.

Figure A.1: Autocorrelation function $R$ for $E_{\alpha}$ (a),
$c_{\alpha}^{0}$ (b), $X_{\alpha}(t)$ (c), and $Y_{\alpha}(t)$ (d). Results
are for the MBL model with $h=0.50$ (red), $h=3.75$ (orange) and $h=6.00$
(blue). For panels (c, d), we consider $t/t_{H}\approx 4.553$ ($h=0.5$),
$t/t_{H}\approx 11.701$ ($h=3.75$), and $t/t_{H}\approx 13.334$ ($h=6.0$).
Dashed lines determine a confidence interval indicated in the main text.
System size is $L=18$.

Apparently the autocorrelations in $c_{\alpha}^{0}$ are the origin of
autocorrelations in $X_{\alpha}(t)$ and $Y_{\alpha}(t)$, as shown in Figs.
A.1.1(c) and A.1.1(d). The behavior of $R_{X_{\alpha}}(k)$ and
$R_{Y_{\alpha}}(k)$ is similar to the one just described for
$R_{c_{\alpha}^{0}}(k)$, having values larger than the confidence interval
only for $h=3.75$ and $h=6.0$, while for $h=0.5$ the values are inside the
confidence interval for any lag $k$. These observations are fully consistent
with Figs. 3 and 5 for the return probability and the discussion around them.
They confirm the conditions for the fulfillment of the CLT in the context of
Eq. (8).

Our results with smaller system sizes $L=8-16$, not shown here, are consistent
with the ones for $L=18$. Of course, the smaller the system size the wider the
confidence interval of 95%. We leave a deeper analysis on autocorrelations for
a future work.

Let us make a final comment. Certainly, it is known that eigenvalues of the
MBL model in the thermalizing regime, say with $h=0.5$, show level repulsion
just as eigenvalues of full random matrices from a Gaussian orthogonal
ensemble [76]. This is what in the context of RMT is known as a correlated
energy spectrum. Then, our claims based on the behavior of $R_{E_{\alpha}}(k)$
depicted in Fig. A.1.1 could appear as intriguing and surprising. However, one
should remind that the autocorrelation function $R_{x}(k)$ measures only
linear autocorrelations, and the lack of those kind of autocorrelations
appears as enough to achieve fulfillment of the CLT and consequently universal
probability distribution functions.

##  Appendix B System-size dependence

Our results in the main text where carried out for a fixed system size,
$L=18$, which certainly is not that big, but large enough for standard exact
diagonalization. Here we explore a much smaller system size, say $L=6$. Figure
B.1.1 shows the pdf of the return probability $\text{RP}(t)$ and the spin
autocorrelation function, $I(t)$. The observed images are similar to the lower
panels of Fig. 3 for $\text{RP}(t)$ and upper panels of Fig. 8 for $I(t)$.
Exponential [Fig. B.1.1 (a)] and Gaussian [Fig. B.1.1 (d)] pdfs for
$\text{RP}(t)$ and $I(t)$ when the system is in the ergodic phase with
$h=0.5$. For stronger disorder strength, $h=6.0$, the analysis leads to pdfs
that do not conform with exponential [Fig. B.1.1 (c)] or Gaussian [Fig. B.1.1
(f)]. The picture for $h=3.75$ goes in line with the one for $h=6.0$.

Figure B.1: Probability distribution function of the return probability (upper
panels) and imbalance (lower panels) for the MBL model with $L=6$. $h=0.5$ in
panels (a) and (d), $h=3.75$ in panels (b) and (e), $h=6.0$ in panels (c) and
(f). Time intervals as in Figs. 3 and 8. Solid curves in the upper panels are
Eq. (9) with $\lambda=1$, while in the lower panels they are Gaussian
fittings.

The considered system size, $L=6$, is the smallest one we can take into
account to avoid what apparently are finite size effects. Our results for
$L=4$ (not shown here) with half-filling display pdfs that do not conform for
any disorder strength with the Gaussian scenario that we have established for
interacting quantum systems. Something comparable is expected to happen for
other models like the clean $XXZ$ model.

##  Appendix C Disorder realization dependence

We devote this Appendix to the analysis of the dependence on the disorder
realization of our results for the MBL model. It is justified since one reader
could think, with certain reason, that independent random disorder
realizations will result in different pdfs. Certainly, it is known that the
time evolution of physical quantum models depends on the initial state, then
in the case of the MBL model, for each disorder realization, if the initial
state is not fixed by hand, then a different initial state could be picked up
and consequently a different pdf could be obtained for the quantities and
observables that we are addressing in this work. As stated in the main text,
to study the pdf we choose one initial state with energy closest to zero, this
is our fixed criteria.

Figure C.1: Pdf of the return probability for the MBL model from two different
disorder realizations (left and right). $h=0.5$ in panels (a, d), $h=3.75$ in
panels (b, e), $h=6.0$ in panels (c, f). Solid curves are the exponential pdf
given by Eq. (9) with $\lambda=1$. $L=16$.

Figure C.1.1 shows the pdf of the return probability from two different random
realizations of the MBL model and three different disorder strengths. There is
seen that the behavior is in general very similar for both realizations (left
and right). Specifically, the pdf for $h=0.5$ is the universal exponential
distribution, meanwhile for $h=3.75$ and $h=6.0$ the pdf does not conform with
the exponential, although for this two later values the pdf is different for
each disorder realization. This can be explained because in the ergodic phase
and close to the middle of the energy spectrum, the structure of basis states
(one of which is our initial state) is pretty similar and stable,
independently of the disorder realization. In contrast, for $h$ around $3.75$
and larger, the distribution of basis states in energy is disperse and the
chosen initial state at the middle of the spectrum will have a different
structure depending of the disorder realization, but its components will
always contain a degree of correlations strong enough to avoid the fullfilment
of the central limit theorem.

## References

  * Clusel and Bertin [2008] M. Clusel and E. Bertin, Global fluctuations in physical systems: A subtle interplay between sum and extreme value statistics, Int. J. Mod. Phys. B 22, 3311 (2008).
  * Peres [1984] A. Peres, Stability of quantum motion in chaotic and regular systems, Phys. Rev. A 30, 1610 (1984).
  * Venuti and Zanardi [2010] L. C. Venuti and P. Zanardi, Universality in the equilibration of quantum systems after a small quench, Phys. Rev. A 81, 032113 (2010).
  * Venuti and Zanardi [2013] L. C. Venuti and P. Zanardi, Gaussian equilibration, Phys. Rev. E 87, 012106 (2013).
  * Zangara  _et al._ [2013] P. R. Zangara, A. D. Dente, E. J. Torres-Herrera, H. M. Pastawski, A. Iucci, and L. F. Santos, Time fluctuations in isolated quantum systems of interacting particles, Phys. Rev. E 88, 032913 (2013).
  * Gring  _et al._ [2012] M. Gring, M. Kuhnert, T. Langen, T. Kitagawa, B. Rauer, M. Schreitl, I. Mazets, D. A. Smith, E. Demler, and J. Schmiedmayer, Relaxation and prethermalization in an isolated quantum system, Science 337, 1318 (2012).
  * Gramsch and Rigol [2012] C. Gramsch and M. Rigol, Quenches in a quasidisordered integrable lattice system: Dynamics and statistical description of observables after relaxation, Phys. Rev. A 86, 053615 (2012).
  * He  _et al._ [2013] K. He, L. F. Santos, T. M. Wright, and M. Rigol, Single-particle and many-body analyses of a quasiperiodic integrable system after a quench, Phys. Rev. A 87, 063637 (2013).
  * Gong  _et al._ [2021] M. Gong, G. D. de Moraes Neto, C. Zha, Y. Wu, H. Rong, Y. Ye, S. Li, Q. Zhu, S. Wang, Y. Zhao, _et al._ , Experimental characterization of the quantum many-body localization transition, Phys. Rev. Res. 3, 033043 (2021).
  * Doggen  _et al._ [2018] E. V. H. Doggen, F. Schindler, K. S. Tikhonov, A. D. Mirlin, T. Neupert, D. G. Polyakov, and I. V. Gornyi, Many-body localization and delocalization in large quantum chains, Phys. Rev. B 98, 174202 (2018).
  * Maksymov  _et al._ [2020] A. Maksymov, P. Sierant, and J. Zakrzewski, Many-body localization in a one-dimensional optical lattice with speckle disorder, Phys. Rev. B 102, 134205 (2020).
  * Luitz [2016] D. J. Luitz, Long tail distributions near the many-body localization transition, Phys. Rev. B 93, 134201 (2016).
  * Yu  _et al._ [2016] X. Yu, D. J. Luitz, and B. K. Clark, Bimodal entanglement entropy distribution in the many-body localization transition, Phys. Rev. B 94, 184202 (2016).
  * L. Mandelstam [1945] I. T. L. Mandelstam, The uncertainty relation between energy and time in nonrelativistic quantum mechanics, J. Phys. USSR 9, 249–254 (1945).
  * Bhattacharyya [1983] K. Bhattacharyya, Quantum decay and the Mandelstam-Tamm-energy inequality, J. Phys. A 16, 2993 (1983).
  * Khalfin [1958] L. A. Khalfin, Contribution to the decay theory of a quasi-stationary state, Sov. Phys. JETP 6, 1053 (1958).
  * Misra and Sudarshan [1977] B. Misra and E. G. Sudarshan, The Zeno’s paradox in quantum theory, J. Math. Phys. 18, 756 (1977).
  * Pfeifer [1993] P. Pfeifer, How fast can a quantum state change with time?, Phys. Rev. Lett. 70, 3365 (1993).
  * Ufink [1993] J. Ufink, The rate of evolution of a quantum state, Am. J. Phys. 61, 935 (1993).
  * Margolus and Levitin [1998] N. Margolus and L. B. Levitin, The maximum speed of dynamical evolution, Phys. D: Nonlinear Phenom. 120, 188 (1998).
  * Giovannetti  _et al._ [2003] V. Giovannetti, S. Lloyd, and L. Maccone, Quantum limits to dynamical evolution, Phys. Rev. A 67, 052109 (2003).
  * Schofield  _et al._ [1995] S. A. Schofield, P. G. Wolynes, and R. E. Wyatt, Computational study of many-dimensional quantum energy flow: From action diffusion to localization, Phys. Rev. Lett. 74, 3720 (1995).
  * De Oliveira and Pellegrino [1999] C. R. De Oliveira and G. Q. Pellegrino, Quantum return probability for substitution potentials, J. Phys. A Math. Gen. 32, L285 (1999).
  * Miyamoto [2001] M. Miyamoto, A generalized Weyl relation approach to the time operator and its connection to the survival probability, J. Math. Phys. 42, 1038 (2001).
  * Krauss and Dent [2008] L. M. Krauss and J. Dent, Late time behavior of false vacuum decay: Possible implications for cosmology and metastable inflating states, Phys. Rev. Lett. 100, 171301 (2008).
  * Xu and Liu [2008] X. Xu and F. Liu, Continuous-time quantum walks on Erdös–Rényi networks, Phys. Lett. A 372, 6727 (2008).
  * Mülken and Blumen [2011] O. Mülken and A. Blumen, Continuous-time quantum walks: Models for coherent transport on complex networks, Phys. Rep. 502, 37 (2011).
  * Ampadu [2012] C. Ampadu, Return probability of the Fibonacci quantum walk, Commun. Theor. Phys. 58, 220 (2012).
  * Riascos and Mateos [2015] A. Riascos and J. L. Mateos, Fractional quantum mechanics on networks: Long-range dynamics and quantum transport, Phys. Rev. E 92, 052814 (2015).
  * del Campo  _et al._ [2017] A. del Campo, J. Molina-Vilaplana, and J. Sonner, Scrambling the spectral form factor: Unitarity constraints and exact results, Phys. Rev. D 95, 126008 (2017).
  * del Campo  _et al._ [2018] A. del Campo, J. Molina-Vilaplana, L. F. Santos, and J. Sonner, Decay of a thermofield-double state in chaotic quantum systems, Eur. Phys. J. Spec. Top. 227, 247 (2018).
  * Lerma-Hernández  _et al._ [2019] S. Lerma-Hernández, D. Villaseñor, M. Bastarrachea-Magnani, E. Torres-Herrera, L. F. Santos, and J. Hirsch, Dynamical signatures of quantum chaos and relaxation timescales in a spin-boson system, Phys. Rev. E 100, 012218 (2019).
  * de la Cruz  _et al._ [2020] J. de la Cruz, S. Lerma-Hernández, and J. G. Hirsch, Quantum chaos in a system with high degree of symmetries, Phys. Rev. E 102, 032208 (2020).
  * Leviandier  _et al._ [1986] L. Leviandier, M. Lombardi, R. Jost, and J. P. Pique, Fourier transform: A tool to measure statistical level properties in very complex spectra, Phys. Rev. Lett. 56, 2449 (1986).
  * Gruebele and Wolynes [2004] M. Gruebele and P. G. Wolynes, Vibrational energy flow and chemical reactions, Acc. Chem. Res. 37, 261 (2004).
  * Rothe  _et al._ [2006] C. Rothe, S. Hintschich, and A. Monkman, Violation of the exponential-decay law at long times, Phys. Rev. Lett. 96, 163601 (2006).
  * Wilkinson  _et al._ [1997] S. R. Wilkinson, C. F. Bharucha, M. C. Fischer, K. W. Madison, P. R. Morrow, Q. Niu, B. Sundaram, and M. G. Raizen, Experimental evidence for non-exponential decay in quantum tunnelling, Nature 387, 575 (1997).
  * Gherardini  _et al._ [2017] S. Gherardini, C. Lovecchio, M. M. Müller, P. Lombardi, F. Caruso, and F. S. Cataliotti, Ergodicity in randomly perturbed quantum systems, Quantum Sci. Technol. 2, 015007 (2017).
  * Singh  _et al._ [2019] K. Singh, C. J. Fujiwara, Z. A. Geiger, E. Q. Simmons, M. Lipatov, A. Cao, P. Dotti, S. V. Rajagopal, R. Senaratne, T. Shimasaki, _et al._ , Quantifying and controlling prethermal nonergodicity in interacting Floquet matter, Phys. Rev. X 9, 041021 (2019).
  * Jordan and Wigner [1928] P. Jordan and E. P. Wigner, Über das Paulische äquivalenzverbot, Z. Phys. 47, 631 (1928).
  * Lieb  _et al._ [1961] E. Lieb, T. Schultz, and D. Mattis, Two soluble models of an antiferromagnetic chain, Ann. Phys. 16, 407 (1961).
  * Bethe [1931] H. Bethe, Zur theorie der metalle, Z. Phys. 71, 205 (1931).
  * Hulthén [1938] L. Hulthén,  _Über das austauschproblem eines kristalles_ , Ph.D. thesis, Almqvist & Wiksell (1938).
  * Des Cloizeaux and Pearson [1962] J. Des Cloizeaux and J. Pearson, Spin-wave spectrum of the antiferromagnetic linear chain, Phys. Rev. 128, 2131 (1962).
  * Yang and Yang [1966] C.-N. Yang and C.-P. Yang, One-dimensional chain of anisotropic spin-spin interactions. I. Proof of Bethe’s hypothesis for ground state in a finite system, Phys. Rev. 150, 321 (1966).
  * Kirillov and Reshetikhin [1987] A. N. Kirillov and N. Y. Reshetikhin, Exact solution of the integrable $XXZ$ Heisenberg model with arbitrary spin. I. The ground state and the excitation spectrum, J. Phys. A Math. Gen. 20, 1565 (1987).
  * Oganesyan and Huse [2007] V. Oganesyan and D. A. Huse, Localization of interacting fermions at high temperature, Phys. Rev. B 75, 155111 (2007).
  * Pal and Huse [2010] A. Pal and D. A. Huse, Many-body localization phase transition, Phys. Rev. B 82, 174411 (2010).
  * Berkelbach and Reichman [2010] T. C. Berkelbach and D. R. Reichman, Conductivity of disordered quantum lattice models at infinite temperature: Many-body localization, Phys. Rev. B 81, 224429 (2010).
  * Kjäll  _et al._ [2014] J. A. Kjäll, J. H. Bardarson, and F. Pollmann, Many-body localization in a disordered quantum Ising chain, Phys. Rev. Lett. 113, 107204 (2014).
  * Luitz  _et al._ [2015] D. J. Luitz, N. Laflorencie, and F. Alet, Many-body localization edge in the random-field Heisenberg chain, Phys. Rev. B 91, 081103 (2015).
  * Devakul and Singh [2015] T. Devakul and R. R. Singh, Early breakdown of area-law entanglement at the many-body delocalization transition, Phys. Rev. Lett. 115, 187201 (2015).
  * Sierant  _et al._ [2020] P. Sierant, D. Delande, and J. Zakrzewski, Thouless time analysis of Anderson and many-body localization transitions, Phys. Rev. Lett. 124, 186601 (2020).
  * Schreiber  _et al._ [2015] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, 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 349, 842 (2015).
  * Bordia  _et al._ [2017] P. Bordia, H. Lüschen, U. Schneider, M. Knap, and I. Bloch, Periodically driving a many-body localized quantum system, Nat. Phys. 13, 460 (2017).
  * Richerme  _et al._ [2014] 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 511, 198 (2014).
  * Berry and Tabor [1977] M. V. Berry and M. Tabor, Level clustering in the regular spectrum, Proc. R. Soc. Lond. A. Math. Phys. Sc. 356, 375 (1977).
  * Bohigas  _et al._ [1984] O. Bohigas, M.-J. Giannoni, and C. Schmit, Characterization of chaotic quantum spectra and universality of level fluctuation laws, Phys. Rev. lett. 52, 1 (1984).
  * Prange [1997] R. E. Prange, The spectral form factor is not self-averaging, Phys. Rev. Lett. 78, 2280 (1997).
  * Schiulaz  _et al._ [2020] M. Schiulaz, E. J. Torres-Herrera, F. Pérez-Bernal, and L. F. Santos, Self-averaging in many-body quantum systems out of equilibrium: Chaotic systems, Phys. Rev. B 101, 174312 (2020).
  * Torres-Herrera  _et al._ [2020a] E. J. Torres-Herrera, G. De Tomasi, M. Schiulaz, F. Pérez-Bernal, and L. F. Santos, Self-averaging in many-body quantum systems out of equilibrium: Approach to the localized phase, Phys. Rev. B 102, 094310 (2020a).
  * Torres-Herrera  _et al._ [2020b] E. J. Torres-Herrera, I. Vallejo-Fabila, A. J. Martínez-Mendoza, and L. F. Santos, Self-averaging in many-body quantum systems out of equilibrium: Time dependence of distributions, Phys. Rev. E 102, 062126 (2020b).
  * Torres-Herrera and Santos [2015] E. J. Torres-Herrera and L. F. Santos, Dynamics at the many-body localization transition, Phys. Rev. B 92, 014208 (2015).
  * Luitz  _et al._ [2016] D. J. Luitz, N. Laflorencie, and F. Alet, Extended slow dynamical regime close to the many-body localization transition, Phys. Rev. B 93, 060201 (2016).
  * Fischer  _et al._ [2016] M. H. Fischer, M. Maksymenko, and E. Altman, Dynamics of a many-body-localized system coupled to a bath, Phys. Rev. Lett. 116, 160401 (2016).
  * Levi  _et al._ [2016] E. Levi, M. Heyl, I. Lesanovsky, and J. P. Garrahan, Robustness of many-body localization in the presence of dissipation, Phys. Rev. Lett. 116, 237203 (2016).
  * Torres-Herrera and Santos [2017a] E. J. Torres-Herrera and L. F. Santos, Dynamical manifestations of quantum chaos: Correlation hole and bulge, Philos. Trans. Royal Soc. A 375, 20160434 (2017a).
  * Torres-Herrera  _et al._ [2018] E. J. Torres-Herrera, A. M. García-García, and L. F. Santos, Generic dynamical features of quenched interacting quantum systems: Survival probability, density imbalance, and out-of-time-ordered correlator, Phys. Rev. B 97, 060303(R) (2018).
  * Torres-Herrera and Santos [2019] E. J. Torres-Herrera and L. F. Santos, Signatures of chaos and thermalization in the dynamics of many-body quantum systems, Eur. Phys. J. Spec. Top. 227, 1897 (2019).
  * Lezama  _et al._ [2021] T. L. Lezama, E. J. Torres-Herrera, F. Pérez-Bernal, Y. B. Lev, and L. F. Santos, Equilibration time in many-body quantum systems, Phys. Rev. B 104, 085117 (2021).
  * Siddiqui [1962] M. M. Siddiqui, Some problems connected with Rayleigh distributions, J. Res. Natl. Bur. Stand. (U. S.) 66, 167 (1962).
  * Hilhorst [2009] H. Hilhorst, Central limit theorem for correlated variables: Some critical remarks, Brazilian J. Phys. 39, 371 (2009).
  * Aurich and Steiner [1999] R. Aurich and F. Steiner, Temporal quantum chaos, Int. J. Mod. Phys. B 13, 2361 (1999).
  * Kunz [1999] H. Kunz, The probability distribution of the spectral form factor in random matrix theory, J. Phys. A Math. Gen. 32, 2171 (1999).
  * Kunz [2002] H. Kunz, Quantum dynamics and random matrix theory, Int. J. Mod. Phys. B 16, 2003 (2002).
  * Serbyn and Moore [2016] M. Serbyn and J. E. Moore, Spectral statistics across the many-body localization transition, Phys. Rev. B 93, 041424 (2016).
  * Torres-Herrera and Santos [2017b] E. J. Torres-Herrera and L. F. Santos, Extended nonergodic states in disordered many-body quantum systems, Ann. Phys. (Berlin) 529, 1600284 (2017b).
  * Tikhonov and Mirlin [2021] K. S. Tikhonov and A. D. Mirlin, Eigenstate correlations around the many-body localization transition, Phys. Rev. B 103, 064204 (2021).
  * Lüschen  _et al._ [2017] H. P. Lüschen, P. Bordia, S. Scherg, F. Alet, E. Altman, U. Schneider, and I. Bloch, Observation of slow dynamics near the many-body localization transition in one-dimensional quasiperiodic systems, Phys. Rev. Lett. 119, 260401 (2017).
  * Willemain and Desautels [1993] R. Willemain and A. Desautels, A method to generate autocorrelated uniform random numbers, J. Stat. Comput. Simul. 45, 23 (1993).
  * Vallejo-Fabila and Torres-Herrera [2022] I. Vallejo-Fabila and E. J. Torres-Herrera, Effects of autocorrelated disorder on the dynamics in the vicinity of the many-body localization transition, Phys. Rev. B 106, L220201 (2022).
  * Chirikov and Shepelyansky [1995] B. Chirikov and D. Shepelyansky, Shnirelman peak in level spacing statistics, Phys. Rev. Lett. 74, 518 (1995).
  * Wegner [1980] F. Wegner, Inverse participation ratio in $2+\varepsilon$ dimensions, Z. Phys. B Con. Mat. 36, 209 (1980).