# The internal clock of many-body delocalization

Ferdinand Evers  Institute of Theoretical Physics, University of Regensburg,
D-93053 Regensburg, Germany  Ishita Modak  Department of Physics, Indian
Institute of Technology Bombay, Mumbai 400076, India  Soumya Bera  Department
of Physics, Indian Institute of Technology Bombay, Mumbai 400076, India

###### Abstract

After a decade of many claims to the opposite, there now is a growing
consensus that generic disordered quantum wires, e.g. the XXZ-Heisenberg
chain, do not exhibit many-body localization (MBL) - at least not in a strict
sense within a reasonable window of disorder values $W$. Specifically,
computational studies of short wires exhibit an extremely slow but
unmistakable flow of physical observables with increasing time and system size
(“creep”) that is consistently directed away from (strict) localization. Our
work sheds fresh light on delocalization physics: Strong sample-to-sample
fluctuations indicate the absence of a generic time scale, i.e. of a naive
“clock rate”; however, the concept of an “internal clock” survives, at least
in an ensemble sense. Specifically, we investigate the relaxation of the
imbalance $\mathscr{I}(t)$ and its temporal fluctuations $\mathscr{F}(t)$, the
entanglement and Renyi entropies, $\mathscr{S}_{\mathrm{e}}(t)$ and
$\mathscr{S}_{2}(t)$, in a 1D system of interacting disordered fermions. We
observe that adopting $\mathscr{S}_{\mathrm{e}}(t),\mathscr{S}_{2}(t)$ as a
measure for the internal time per sample reduces the sample-to-sample
fluctuations but does not eliminate them. However, a (nearly) perfect collapse
of the average $\overline{\mathscr{I}}(t)$ and $\overline{\mathscr{F}}(t)$ for
different $W$ is obtained when plotted against
$\overline{\mathscr{S}}_{\mathrm{e}}(t)$ or $\overline{\mathscr{S}}_{2}(t)$,
indicating that the average entropy appropriately models the ensemble-averaged
internal clock. We take the tendency for faster-than-logarithmic growth of
$\overline{\mathscr{S}}_{\mathrm{e}}(t)$ together with smooth dependency on
$W$ of all our observables within the entire simulation window as support for
the cross-over scenario, discouraging an MBL transition within the traditional
parametric window of computational studies.

##  I Introduction

Many-body localization (MBL) is a spatio-temporal phenomenon that is believed
to exist in interacting fermion systems at strong enough disorder [1, 2, 3, 4,
5, 6, 7, 8, 9]. It manifests itself as a strongly reduced (preMBL) or even
fully inhibited tendency (proper MBL) towards thermalization at large enough
disorder strength. The requirement of strong disorder implies that all
remnants of relaxation dynamics in the preMBL-regime are necessarily very
slow, reminiscent of the familiar behavior of glasses. “Slowness” originates
from the fact that long-time behavior typically is dominated by collective
reorganization processes in a highly disordered many-body-energy landscape.
The salient collective events are associated with a broad distribution of time
scales; there is no single characteristic rate that would lend itself as a
measure of time [10, 11].

We briefly elaborate on the issue of time- and length scales by formulating a
scenario: Consider an ensemble of disordered quantum wires of a finite length
$L$. At strong enough disorder, typical many-body states can be thought about
essentially as single (dressed) Slater-determinants constructed from localized
single-particle wavefunctions; this is the jest of the celebrated concept of
local integrals of motion (LIOM) [12, 13, 14, 15, 16]. While many, perhaps
even most, wires may follow the LIOM-paradigm, a fraction ${\mathfrak{q}}(L)$
of the wires will not; they form a thermalizing sub-ensemble that we refer to
as “ergodic bubbles” [17, 18]. A new ensemble of wires with the length $2L$
can be formed by combining samples giving rise to a new fraction
${\mathfrak{q}}(2L)$. The evolution of ${\mathfrak{q}}(2L)$ with system size
has been investigated in RG-studies of toy models [14, 19, 20, 21, 22, 23,
24].

We emphasize two important implications of this rough picture for
observations, numerical and experimental, in thermalizing phases when growing
the system size at ${\mathfrak{q}}(L)\ll 1$: (i) strong finite size effects
occur because thermalization sets in only after a certain length-scale have
been reached long enough to include an ergodic bubble; (ii) the time scale for
relaxation can be exponentially long because it reflects how long it takes a
bubble to thermalize an almost localized and large sample region. We refer to
this slow relaxation process as ‘creep’. Since creep emerges from
destabilizing an interacting but localized (i.e. LIOM-dominated) sample region
via remote thermal bubbles, creep indicates the dominant mechanism for many-
body de-localization and is prevalent in the preMBL-regime.

Signals of creep are exponentially long observation times required in
numerical and experimental simulations as well as an appreciable dependency of
transport-sensitive observables on the simulation volume. We interpret the
gradual increase of the estimate for the disorder strength beyond which the
localization transition was expected to occur as a signature of creep: while
early works favored a value $W{\approx}3.8$ [25, 26], later authors reported
significant finite size effects and favored larger values $W{\approx 4.5-5.5}$
[27, 28, 29].

Creep manifests, in particular, in the spatio-temporal behavior of correlation
functions; in the preMBL-regime, they exhibit a spatial decay much slower than
exponential with tails slowly increasing with time; for a discussion of how
creep appears in various observables see Weiner  _et al._ [30]. It is only
recently that creep in its diverse manifestations and the resulting
implications for the (pre)MBL phenomenology in finite size systems have been
appreciated: A weak tendency towards equilibration has been analyzed in the
relaxation dynamics even at large disorder in the local charge density, the
sublattice imbalance, and the density matrix [31, 30, 32, 33, 34, 35]; the
evolution of the entanglement and the number entropy has been investigated
[36, 33, 37, 38, 39, 40, 41, 42, 43]; attempts have been made to interpret the
shift of intersection points (“critical points”) and the flow of the spectral
function with increasing system sizes [44, 45, 46, 28, 47, 48, 49]. The status
as we see it is summarized in the phase diagram Fig. 1.

Figure 1: Schematic phase diagram of the disordered t-V (isotropic Heisenberg)
model (1)  with $V$=1. Shaded region: thermalizing regime exhibiting
(transient) power-laws, e.g., for the decay of imbalance (or time-
fluctuations) and growth of entanglement and second Renyi entropy.
$W_{c}{\gtrsim 4.5}$: crossover from accelerating to decelerating dynamics as
indicated by the temporal increase/decrease of an effective diffusion exponent
$\beta(t)$ [31, 30]. Beyond $W_{\text{freeze}}{\approx}10$: fractal dimension
of the dynamically active part of the many-body Hilbert space vanishes
(freezing) [50]. This region is largely unexplored in finite size and time
numerical studies and hosts the proper MBL phase, if it exists.

An implication of creep is the absence of characteristic time scales; one
rather expects broad distributions of rates characterizing dynamical
processes. This poses the question of at which time it is meaningful to
compare the dynamical status of two samples that nominally belong to the same
ensemble but differ by the specific disorder realization. We here pursue the
implications of the following hypothesis: The time evolution in both samples
can at least partially be synchronized when introducing an ‘internal sample
time’. We here propose to model the internal time by a form of entropy. While
the model has its limitations when adopted on a per-sample level, it works
incredibly well for the disorder-averaged ‘internal ensemble time’. To
illustrate the power of the concept, we plot in Fig. 2 the average charge
imbalance $\overline{\mathscr{I}}(t)$ over the entanglement entropy
$\overline{\mathscr{S}}_{\mathrm{e}}(t)$: a data collapse is observed in a
wide time window for a range of disorder values situated in the (transient)
sub-diffusive regime (Fig. 1).

On a more intuitive level, two considerations motivate us to consider an
entropic entity ${\mathscr{S}}_{L}(t)$ as a model of the internal clock of a
system with spatial extension $L$: (i) The time evolution of the entanglement
entropy is, in a sense, the fastest and most stable relaxation process: It is
relatively fast, because, unlike the redistribution of, e.g., energy and
particle number it is not restricted by a local conservation law: even in the
proper MBL phase, one expects a logarithmic-in-time growth of the entanglement
entropy, when energy and particle densities have long ceased to relax. The
saturation dynamics of physical observables will thus be represented by a
plateau when plotted against $\mathscr{S}_{\mathrm{e}}$ in the thermodynamic
limit (ii) On a more formal level, the rate of entropy production can be
considered as a generalized driving force for the relaxation of (ensemble-
averaged) physical observables [51].

Figure 2: Evolution of the ensemble-averaged imbalance after quenching a
product (Neél)-state in a disordered wire of length $L{=}16,20,24,26$ at four
moderate disorder strengths $W{=}1.5,2.0,2.5,3.0$. Time is measured in terms
of the average entanglement entropy in units of the Page-value
${\mathscr{S}}_{\text{Page}}{\coloneqq}\ln 2(L/2){-}1/2$.  Inset: Scaling
factor $f_{L}(W)$ as a function of $W$. By definition, rescaling the ordinate
with $f_{L}(W)$ collapses imbalance traces for different disorder $W$. The
factor $f_{L}(W)$ has been extracted from our data performing a standard
scaling analysis. It consists of rescaling the y-axis until the best data
collapse is observed.

Subjecting the synchronization concept to a further test on a second
observable, we apply it to the time fluctuations of the local density,
$\overline{\mathscr{F}}(t)$. Again an excellent scaling collapse is observed,
demonstrating the power of the concept.

Further, we find that
$\overline{\mathscr{F}}\propto\overline{\mathscr{S}}_{\mathrm{e}}^{-\rho}$,
where the exponent depends on the disorder strength, $\rho(W)$; a similar
relationship also holds for typical values of ${\mathscr{F}}$. The
corresponding exponent functions, $\rho_{\text{ave}}(W)$ and
$\rho_{\text{typ}}(W)$, reveal opposing trends with increasing disorder
strength, in particular, they intersect at a disorder value
$W_{{\mathscr{F}}}$. Since $W_{{\mathscr{F}}}$ is sufficiently close to
$W_{c}$ we take this as fresh evidence in favor of the existence of (at least)
two subphases in the thermalizing regime indicated in the phase diagram Fig 1
(‘Subphase’ is meant here in a weak sense indicating that the thermalizing
phase has regions with qualitatively different relaxation behavior that are
connected via a cross-over.).

##  II Theoretical setting

###  II.1 Model and Method

We consider spinless fermions within the $t{-}V$ model

$\displaystyle\mathscr{H}=$
$\displaystyle-\frac{1}{2}\sum_{i=1}^{L-1}c_{i}^{\dagger}c_{i+1}+\text{h.c}+\sum_{i=1}^{L}\epsilon_{i}(n_{i}-1/2)$
$\displaystyle+V\sum_{i=1}^{L-1}(n_{i}-1/2)(n_{i+1}-1/2),$ (1)

where $n_{i}\coloneqq c_{i}^{\dagger}c_{i}$, $i$ denotes the site index, $L$
the system size and $V$ the nearest neighbor interaction strength; $V{=}1$
throughout this work. The on-site potentials, $\epsilon_{i}$, are uncorrelated
and distributed within $[-W,W]$; for comparison,  we present a case of
correlated disorder in appendix A.3. All calculations are done at half-
filling, $N/L{=}1/2$. We study quench dynamics starting from the (product
state) charge density wave state $|\Psi\rangle=|1010\ldots\rangle$.  To
investigate possible initial-state dependencies, a brief investigation of an
alternative initializing state, $|\Psi\rangle=|1100\ldots\rangle$, is reported
in appendix A.1.  Time propagation employs the standard Chebyshev expansion of
the time evolution operator, which reduces the exponentiation of $\mathscr{H}$
to several sparse matrix-vector multiplications [52, 31]. A comparison of the
Chebyshev method with exact time evolution is presented in appendix D.

###  II.2 Observables

In most of this work, we focus on two main observables. For the “internal
clock” we adopt the entanglement entropy

$\displaystyle\mathscr{S}_{\mathrm{e}}\coloneqq-\text{Tr}_{\text{A}}\rho_{\text{A}}\ln\rho_{\text{A}},$
(2)

where $\rho_{\text{A}}$ denotes the density operator of subsystem A as
obtained after integrating out the complement of A, i.e. subsystem B. We use
an equal partition $L_{\text{A}}{=}L_{\text{B}}{=}L/2$. With an eye to the
measurement of internal times, we mention that $\mathscr{S}_{\mathrm{e}}(t)$
behaves similarly to the second Renyi-entropy
${\mathscr{S}}_{2}\coloneqq-\ln\text{Tr}_{\text{A}}\hat{\rho}_{\text{A}}^{2},$
see appendix C. Recently, $\mathscr{S}_{\mathrm{e}}$ [53] as well as the Renyi
entropy ${{\mathscr{S}}_{2}}$ have been measured [54] in ion trap experiments.

Our second observable is the sublattice imbalance

$\displaystyle{\mathscr{I}}(t)\coloneqq\sum_{j}^{L}(-)^{j}\ n_{j}(t)$ (3)

and derived quantities such as the density fluctuations

$\displaystyle{\mathscr{F}}(t)\coloneqq
1/L\sum_{j=1}^{L}\\{[n_{j}(t)-\\{n_{j}(t)\\}_{\Delta t}]^{2}\\}_{\Delta t},$
(4)

where $\\{\ldots\\}_{\Delta t}$ denotes a sliding time window average over a
set of given time traces $n_{j}(t),j=1,\ldots L$, also see [50]; $\Delta
t=12.5$ in our calculations. Observables averaged over different disorder
configurations will be denoted by an overline, e.g.,
$\overline{\mathscr{S}}_{\mathrm{e}}(t),\overline{\mathscr{I}}(t)$.  We also
checked our result with a different observable; the detail is provided in App.
A.2.

##  III Clean case - Results & Discussion

Figs. 3 and 4 display the evolution of $\mathscr{S}_{\mathrm{e}}(t)$ and
${\mathscr{I}}(t)$ in the clean system, i.e., $W{=}0$.

###  III.1 Entanglement entropy - evolution in time and system size

Figure 3: (a) Time evolution of the sublattice imbalance, ${\mathscr{I}}(t)$
and (b) the entanglement entropy, $\mathscr{S}_{\mathrm{e}}(t)$, after a
quench in a clean quantum wire, Eq. (1). $\alpha$ denotes the saturation value
of $\mathscr{S}_{\mathrm{e}}$ in units of the Page value
${\mathscr{S}}_{\text{Page}}$ for the respective system sizes. The inset shows
the infinite-system size extrapolation of the saturation value,
$\mathscr{S}_{\mathrm{e}}^{\infty}(L)$, as read off the main plot. We take the
curvature seen in the data as an indication of corrections to scaling that are
non-analytic in $1/L$, possibly logarithmic [55]. Fitting the trace to Eq. (5)
yields $a_{\text{V}}=0.50(1),a_{\text{S}}=-0.29(4),a=1.8(1).$

We identify three windows of characteristic times the traces
$\mathscr{S}_{\mathrm{e}}(t)$, see Fig. 3b): Let the characteristic velocity
of pair-excitations be $v$; following the kinetic energy given in Eq. (1) one
expects $v\approx 1$. In the intermediate window, $1\lesssim t\lesssim L/2v$,
we have $\mathscr{S}_{\mathrm{e}}(t)=c_{\mathscr{S}}t+\text{subleading
terms}$, reflecting the asymptotic behavior expected for the large-system
limit, $L\to\infty$, of the ballistic system [55]. Based on our data Fig. 3(b)
we obtain $c_{\mathscr{S}}\approx 0.52$. Due to the limited time window and
the unknown subleading terms, reliable error bars are difficult to establish.

Relaxation processes in the short time window occur on rates given by the
(inverse) bandwidth; they reflect the local physics of the clean system. As
seen in Fig. 3a), the imbalance decays on these ultra-short time scales.

At large times, $t\gtrsim L/2v$, the entanglement entropy
$\mathscr{S}_{\mathrm{e}}(t)$ reaches a saturation value,
$\mathscr{S}_{\mathrm{e}}^{\infty}(L)\coloneqq\mathscr{S}_{\mathrm{e}}(L;t\to\infty)$,
that exhibits pronounced system size effects. Its analytical form is well
described by

$\displaystyle\mathscr{S}_{\mathrm{e}}^{\infty}(L)=a_{\text{V}}L/2+a_{\text{S}}\ln(L/2a)+\ldots$
(5)

where $\ldots$ indicate terms that vanish in the thermodynamic limit and $a$
denotes a microscopic scale that accounts for terms $\mathcal{O}(L^{0})$. The
first term incorporates the expectation that the entanglement entropy of the
ballistic system after saturation is an extensive observable, so the leading
behavior is described by a volume law. Based on quantum field theory, one
would expect a prefactor ratio $a_{\text{V}}/c_{\mathscr{S}}=1$, which is
consistent with our numerical estimate, $\approx 96$%, within our
extrapolation uncertainties [55]. An earlier estimate of this ratio, 88%, has
been reported for free fermions (XX model) [56].

The second term in Eq. (5) we read as a manifestation of interface (surface)
effects: we interpret $L^{d-1}$ in the limit $d{\to}1$ as a logarithm. It
accounts for interaction mediated long-range correlations in the one-
dimensional bulk phase that lead to anomalously strong surface effects [55].

For a fully chaotic system, $\mathscr{S}_{\mathrm{e}}^{\infty}(L)$ is expected
to saturate at the Page value,
${\mathscr{S}}_{\text{Page}}\coloneqq\ln(2)L/2-1/2$ 111In Ref. [57] it is
shown that the average entanglement entropy of random pure states to leading
order is of the form  $\displaystyle{\mathscr{S}}_{m,n}=\ln
m-\frac{m^{2}-1}{2mn}+\ldots$ where $n,m$ is the Hilbert space dimensions of
the respective subsystems assumed to be large. For the present case, we have
$n=m=2^{L/2}$ motivating the definition of ${\mathscr{S}}_{\text{Page}}$ given
in the text. . Confirming results of earlier authors [58], we find that clean
systems do not exhaust the Page limit, meaning $\alpha<1$, see Fig. 3b), with
$\alpha(L)\coloneqq\mathscr{S}_{\mathrm{e}}(t\to\infty)/{\mathscr{S}}_{\text{Page}}$
and $\alpha(L\to\infty)=a_{\text{V}}/\ln 2$. For the latter coefficient, we
here report an estimate significantly below unity, $\alpha\approx 0.72$. We
interpret this result, $\alpha{<}1$, as a consequence of the fact that the
clean model exhibits local conservation laws, such as momentum conservation,
that (apparently) can limit the phase-space volume accessible during
relaxation as compared to a fully chaotic situation.

###  III.2 Clock rate and scaling

With an eye on the disordered case and its internal clock, we mention that
entanglement has been introduced as a suitable time unit before already in the
clean case. Calabrese and Cardy [55] observe that the entanglement evolution
$\mathscr{S}_{\mathrm{e}}(L;t)$ corresponding to quenches with different
initializing states collapse onto a master curve after normalization with
respect to the asymptotic value, i.e. when plotting the ratio
$\mathscr{S}_{\mathrm{e}}(L;t)/\mathscr{S}_{\mathrm{e}}^{\infty}(L)$. At
intermediate times, $1\lesssim t\lesssim L/2v$,
$\mathscr{S}_{\mathrm{e}}(t)\propto t$ and therefore the collapse implies
$\dot{\mathscr{S}}_{\mathrm{e}}(t)/\mathscr{S}_{\mathrm{e}}^{\infty}$
independent of time and the initial condition. In this sense, one might say
that the rate of entropy production is controlled by the (saturation)
entanglement $\mathscr{S}_{\mathrm{e}}^{\infty}$.

Taking the idea further, Kim and Huse [58] have postulated a scaling form

$\displaystyle\mathscr{S}_{\mathrm{e}}(t)\approx\mathscr{S}_{\mathrm{e}}^{\infty}(L)\
F(t/\mathscr{S}_{\mathrm{e}}^{\infty}(L))$ (6)

for a situation in which the system size $L$ rather than the initializing
state is varied; $F(x)$ denotes a scaling function. Also here, the saturation
value $\mathscr{S}_{\mathrm{e}}^{\infty}$ sets the rate for entanglement
growth.

In Fig. 4 we test the hypothesis (6) and show the scaling function
corresponding to the data given in Fig. 3. While a reasonable collapse is
observed, corrections to scaling are appreciable that we tentatively assign to
logarithmic corrections, such as displayed in (5). Hence, existing reports on
quenches in clean systems lend support to the main hypothesis of our work,
i.e. that the entanglement entropy can be a good model for an internal clock
in interacting quantum wires.

Figure 4: (Near) scaling collapse of the entanglement entropy after employing
the saturation value,
$\mathscr{S}_{\mathrm{e}}^{\infty}(L){\coloneqq}\mathscr{S}_{\mathrm{e}}(L;t\to\infty)$,
as a measure of time. Alternatively, the inset employs $L/2v$ for units in
order to highlight the finite-size convergence near $t\approx L/2v$ consistent
with expectations based on Ref. [55].

##  IV Disorder - Results & Discussion

###  IV.1 Time evolution of entanglement entropy

Figure 5: Ensemble-averaged entanglement entropy
$\mathscr{S}_{\mathrm{e}}(L,t)$ for system sizes $L{=}16,20,24$ and four
moderate disorder values $W$. (a) Three temporal regimes are identified:
ballistic, (sub-)diffusive, and saturating. (b-c): logarithmic time derivative
$\dot{\mathscr{S}}_{\mathrm{e}}$ reveals a (nearly) power-law asymptotics. The
shaded region indicates the window underlying the fits following Eq. (7).
(Parameters in Tab. 1.) (d) Saturation value is seen in the left panel in
units of the Page limit ${\mathscr{S}}_{\text{Page}}$ over the linear system
size.

We will begin our analysis of disordered systems with the ensemble-averaged
entanglement entropy $\overline{\mathscr{S}}_{\mathrm{e}}(t)$ displayed in
Figs. 5(a) and 6. Similar to the clean case, also here the time traces suggest
introducing three time regimes: In the short time regime, all traces (nearly)
collapse irrespective of disorder $W$ (ballistic period). This regime we will
here leave unattended since disorder effects are (relatively) weak. Our focus
is on the sub-diffusive regime (significant disorder effects, system size
converged) and on the saturation regime (strong effects of disorder and system
sizes).

####  IV.1.1 Sub-diffusive Regime

The sub-diffusive regime will eventually display the asymptotic dynamics of
bulk systems at large enough $L$ and long enough times, i.e. the thermodynamic
limit. Our nomenclature reflects an important result of our investigations
also indicated in Fig. 1: In the parametric window that we have investigated
we consistently observe dynamical signatures, which are expected for the
preMBL-regime; MBL proper is never observed.

Early researchers have characterized the time evolution of observables in the
sub-diffusive regime by power-laws, e.g.
$\overline{\mathscr{S}}_{\mathrm{e}}(t){\propto}t^{1/{z_{\text{ee}}}}$ [59,
60], especially in the range up to $W\simeq 3$ displayed in Fig. 5. However,
as is clearly visible in Fig. 5(a) for the paradigmatic case of entanglement
entropy, the putative time window of power-law dynamics has not yet opened up
in the range of computationally available system sizes. 222To be extracted
from this data is only a lower bound for the dynamical exponent
${z_{\text{ee}}}$, i.e. the slope at the turning point in the (sub-)diffusive
regime. It thus is apparent that finite size effects in disordered wires are
dramatically stronger as compared to the clean case, Fig. 3.

Figure 6: (a) Traces similar to Fig. 5(a) for disorder in the regime of
decelerating dynamics Fig. 1. For $W=4.0,6.0$ system sizes are $L=16,20,24$
and else $L=20,24$. (b) Logarithmic time derivative of
$\overline{\mathscr{S}}_{\mathrm{e}}$ for $W=4.0,6.0$. The plot highlights the
emergence of a power-law increase with effective exponents that quickly grow
with the system size $L$. An increase faster than $\ln t$ rules out MBL within
the respective disorder regime.

The same conclusion also holds at larger disorder in the regime of
decelerating dynamics shown in Fig. 6(a). Notice the the left-hand side (lhs)
curvature in the intermediate time window where the traces are nearly system-
size converged (subdiffusive regime): the crucial observation is that the lhs
curvature becomes more pronounced in a larger window of time with ever growing
system size. As illustrated in Fig. 6(b) the curvature signalizes an emergent
power law with effective exponents that rapidly increase with the growing
system size: a change in system size by 50% increases the effective exponent
at long times by 200% ($W{=}4$) resp. 50% ($W{=}6$). The phase diagram
proposed in Fig. 1 is based on the assumption that the trend observed in Fig.
6 is indicative of the thermodynamic limit and, therefore, does not reverse at
much larger system sizes. While trend reversal towards insulating behavior
(“reentrance”) cannot be rigorously excluded in finite-size studies, we don’t
see any hints pointing toward its existence in numerical work in the MBL
context 333A kind of trend reversal, if only a different one, is observed in
the model of regular-random graphs (RRG) [61]: this reversal is similar to
creep in the sense that the reversed flow in RRG is directed away from the
localized to the ergodic regime. Opposed to this, the reentrance described in
the main text implies a flow reversal towards localizing behavior. . In this
sense, our result Fig. 6 is inconsistent with the expected behavior in the MBL
phase, $\overline{\mathscr{S}}_{\mathrm{e}}(t)\propto\ln t$ [62, 63, 64], and
indicates the absence of many-body localization in the respective window of
disorder values.

We mention the implications of our result Fig. 6 for Refs. [36, 37]. These
authors are quenching random half-filled product states and studying the
number entropies ${\mathscr{S}}_{\mathrm{N}}^{(\alpha)}$ that they relate to
the second Renyi entropy
$\overline{{\mathscr{S}}}_{\mathrm{N}}^{(2)}{\sim}\ln\overline{{\mathscr{S}}}_{2}$
and to the entanglement entropy
$\overline{{\mathscr{S}}}_{\mathrm{N}}^{(1)}{\sim}\ln\overline{\mathscr{S}}_{\mathrm{e}}$
[36, 37]. Fitting their data in the time domain, they conclude a double-
logarithmic growth of the number entropies,
$\overline{{\mathscr{S}}}_{\mathrm{N}}^{(\alpha)}\sim\ln\ln t$. These
statements imply
$\overline{\mathscr{S}}_{\mathrm{e}}(t),\overline{{\mathscr{S}}}_{2}(t)\propto\ln
t$, which taken at face value is inconsistent with our claim of growth of
$\overline{\mathscr{S}}_{\mathrm{e}}(t)$ faster than logarithmic, Fig. 6.

Figure 7: Plot illustrating extremely strong sample-to-sample fluctuations.
Left: Imbalance ${\mathscr{I}}(t)$ over simulation time $t$ for $10$ samples
and four disorder values. The plot illustrates the logarithmically wide
distribution of ${\mathscr{I}}$ at large times. Second column: Logarithmic
derivative (effective exponent $\beta(t)$) of data shown in lhs column. The
logarithmically broad distribution of $\beta$ is illustrated. Third column:
Corresponding entanglement entropy $\mathscr{S}_{\mathrm{e}}(t)$. At large
disorder, at a given time $\mathscr{S}_{\mathrm{e}}$ is logarithmically
distributed. Right: Imbalance as a function of entanglement
${\mathscr{S}}_{\mathrm{e}}(t)$; sample-to-sample fluctuations are reduced if
times with similar values for ${\mathscr{S}}_{\mathrm{e}}(t)$ are compared.

We assign this apparent discrepancy to a significant difference in data
analysis. The $\ln t$ behavior in Refs. [36, 37] has been extracted from time
traces obtained from system sizes up to $L{=}24$. In particular, the fit for
$L{=}24$ was targeting the largest observation times, similar to our data Fig.
6(a). Indeed, in this regime also our data is consistent with a logarithmic
fitting, as is seen by the lack of curvature in the time traces for $L{=}24$
in Fig. 6(a) at $t{\gtrsim}10^{2}$. In this sense, there is no discrepancy
here.

However, as is also seen in Fig. 6(a) the regime of large times is strongly
affected by finite-size effects: with increasing system sizes the traces
quickly develop an ever-growing tendency for lhs-curvature that connects
smoothly to the lhs-curvature already secured in the sub-diffusive time
window. In Fig. 6(b) this flow away from logarithmic towards power-law
dynamics is particularly apparent.

The main result of Refs. [36, 37], i.e. the absence of MBL-proper in the
investigated parameter regime, remains unchallenged here; the reported
$\overline{\mathscr{S}}_{\mathrm{e}}(t)\sim\ln t$ from our point of view
represents a lower bound 444In our analysis we adopt the perspective that
there is no qualitative difference between quenching a Neel-state - as we do -
and a random product state as done in Ref. [36, 37] for the entanglement
dynamics. However, subsequent work noted that averaging over random product
states could make it more difficult to extrapolate the long-time dynamics
[38]..

####  IV.1.2 Saturation regime - disorder enhanced entanglement

Confirming earlier results, we readily observe in Fig. 5 that also in the
disordered situation, $\overline{\mathscr{S}}_{\mathrm{e}}(t)$ does not
exhaust the Page-value, ${\mathscr{S}}_{\text{Page}}$, in the limit of long
observation times [65]. To extract the saturation value,
$\overline{\mathscr{S}}_{\mathrm{e}}(L,W)$, we feed the asymptotic form

$\displaystyle\overline{\mathscr{S}}_{\mathrm{e}}(L,W;t)=\overline{\mathscr{S}}_{\mathrm{e}}^{\infty}(L,W)-c(L,W)/t^{\gamma_{L}(W)}$
(7)

into a three-parameter fit of the data Fig. 5 (a-c). To estimate the amplitude
$c(L,W)$ and the exponent $\gamma_{L}(W)$, we employ the asymptotics of the
time derivative, $\dot{\overline{\mathscr{S}}}_{\mathrm{e}}{}_{W}(L,t)$,
reproduced in Fig. 5(b,c). With this input, the saturation value
$\overline{\mathscr{S}}_{\mathrm{e}}^{\infty}(L,W)$ is fitted adopting (7) to
the data in the main panel, Fig. 5(a); the fitting parameters are reproduced
in the appendix, Tab. 1.

We emphasize that the exponent $\gamma_{L}(W)$ is an effective one. Close
scrutiny of Fig. 5(b,c) reveals a slow time evolution of the slope (exponent)
that is nearly but (perhaps) not fully converged within our window of
observation times. As a consequence, the fitting parameters contain a
systematic error that is not accounted for in the error estimates given in the
appendix, Tab. 1.

The saturation values of the entanglement entropy,
$\overline{\mathscr{S}}_{\mathrm{e}}(L,W)$ are reproduced in Fig. 5(d). A
striking aspect revealed here is disorder enhanced entanglement; there is a
wide regime in the $L,W$-plane, in which the subsystem-entanglement grows with
the disorder. This is always the case in the limit of weak disorder, where it
perhaps is less surprising: weak disorder destroys integrability and therefore
brings the system closer towards chaotic dynamics, so
$\overline{\mathscr{S}}_{\mathrm{e}}(L,W)$ moves upwards towards the Page
limit. This increase has been seen previously in Ref. [66]. Remarkably, the
trend appears to continue even to moderate and strong disorder, $1.5\lesssim
W\lesssim 3.0$, provided the system size is large enough - as seen from the
traces with family parameter $W$ intersecting each other in Fig. 5(d).

###  IV.2 Entanglement entropy as an internal clock

In this subsection, we explore the potential of entanglement entropy as a
model for an internal clock. We begin by pointing out a possible advantage of
using internal clocks. Time traces that display the relaxation of two
different transport-related observables both measured in a given sample tend
to be correlated. For example consider $\mathscr{S}_{\mathrm{e}}(t)$ and
${\mathscr{I}}(t)$: When employing $\mathscr{S}_{\mathrm{e}}(t)$ as a measure
of time on a per-sample basis, sample-to-sample fluctuations should reduce
considerably, because the slow growth of the first will, in general, indicate
strong localization and therefore hints at slow growth of the latter; ergo,
${\mathscr{I}}(\mathscr{S}_{\mathrm{e}})$ fluctuates less than
${\mathscr{I}}(t)$ between different samples and therefore should be easier to
interpret. The reduction of sample-to-sample fluctuations is indeed seen in
Fig. 7, right most column.

We elaborate on the qualitative argument. It can be rephrased employing
concepts underlying the traditional mode-coupling theory [51]. It foresees
that two slow observables in their time evolution can each follow their own
internal clock rate, but only if they are ‘orthogonal’, i.e., sufficiently
decoupled. The generic expectation is that orthogonality is an exception; at
least quantities deriving from the same (conserved) field, e.g. the particle
density, follow the same clock rate.

Adopting $\mathscr{S}_{\mathrm{e}}$ as a model for the internal clock
illustrates a benefit of the concept, i.e. reducing sample-to-sample
fluctuations in Fig. 7 significantly. At the heart of the concept is the
hypothesis that improved models (as compared to $\mathscr{S}_{\mathrm{e}}$)
exist, at least in principle, that further reduce sample-to-sample
fluctuations, nearly eliminating them in generic observables at large enough
$L$. The optimal model for the internal time on a per-sample level still has
to be found. However as we will see in the next paragraphs,
$\overline{\mathscr{S}}_{\mathrm{e}}$ is close to an optimal model for the
disorder-averaged internal clock, i.e. the ‘internal ensemble time’.

####  IV.2.1 Excursion: the arrow of time - progress and regression

A conservative interpretation of $\mathscr{S}_{\mathrm{e}}(t)$ as the internal
system clock requires that the time evolution of the entanglement entropy
exhibits a monotonous growth not only on average but also on a sample-to-
sample basis. This requirement is not obviously met, because the entanglement
entropy is not bound to strict growth in time far enough in non-equilibrium.
As seen in Fig. 7, first and third column, similar to the charge imbalance
[50] ${\mathscr{I}}(t)$, also the time evolution of the entanglement entropy
$\mathscr{S}_{\mathrm{e}}(t)$ exhibits strong sample-to-sample fluctuations.
And indeed, as seen in Fig. 7, the third column at strong disorder $W\gtrsim
4$, $\mathscr{S}_{\mathrm{e}}(t)$ becomes highly non-monotonous. Hence, a
conservative interpretation of $\mathscr{S}_{\mathrm{e}}(t)$ as an internal
system clock will confine itself to moderate disorder $W\lesssim 4$.

We would like to argue, however, that strict monotony is not a necessary
condition for an observable to be a useful measure of time. To this end we
recall that the arrow of time is bound to point in a single direction only in
the macroscopic limit; relaxing finite-size systems may always exhibit
transient periods of entropy shrinkage (‘regression’), so that the condition
of strict monotony can be released, at least in principle.

####  IV.2.2 Collapsing averaged imbalance traces -
$\overline{\mathscr{I}}(t)$ over $\overline{\mathscr{S}}_{\mathrm{e}}(t)$

In the spirit of the preceding paragraphs, we adopt
$\overline{\mathscr{S}}_{\mathrm{e}}(t)$ as a legitimate model for the
internal clock of the ensemble dynamics also at larger disorder values, even
if individual samples display periods of regression. In fact, we will continue
to bravely use $\overline{\mathscr{S}}_{\mathrm{e}}(t)$ also at large disorder
as displayed, e.g., in Fig. 9.

Figure 8: Data from Fig. 2 plotted over the second Renyi-entropy. A good
collapse is achieved here, also, confirming that
$\overline{\mathscr{S}}_{\mathrm{e}}(t)$ and $\overline{{\mathscr{S}}}_{2}(t)$
model the internal clock with comparable quality. A direct comparison of
$\overline{\mathscr{S}}_{\mathrm{e}}(t)$ and $\overline{{\mathscr{S}}}_{2}$ is
given in the appendix, Fig. 16. Parameters: $L{=}16,20,24,26$ at four moderate
disorder strengths $W{=}1.5,2.0,2.5,3.0$.  Inset: prefactor $f_{L}(W)$ as a
function of $W$. The factor $f_{L}(W)$ has been extracted from our data in a
similar process as Fig. 2.

As we have argued, internal clocks are useful to the extent that they reveal
similarities in the time evolution of two different samples that are not
apparent when using the lab-time $t$. A particularly strong case in favor of
internal clocks can be made if similarities can be revealed in the dynamics of
systems - or ensembles - that nominally exhibit strong differences, e.g., in
the disorder strength, so that at least ‘naively’ similarities are not
expected. Such a strong validation of the “internal clock” concept has been
given with Figs. 2 and 8: imbalance traces $\overline{\mathscr{I}}(t)$ taken
from samples of fixed length $L$, and for disorder varying from the moderate
($W{\approx}1.5$) to the beginning of the strong disorder regime
($W{\approx}3$) collapse to a master-curve - within the available window of
observations times.

For collapsing $\overline{\mathscr{I}}$, rescaling of the abscissa is not
required, which reflects the absence of relevant microscopic time scales - at
least within the observation window and outside the short-time regime. The
fact that the data collapse requires a rescaling of the ordinate is not
unexpected, but still merits a remark: Our computations operate in the limit
of infinite temperature, $T^{-1}{=}0$; the corresponding equilibrium density
is spatially homogeneous, so that the equilibrium imbalance vanishes,
${\mathscr{I}}{=}0$. Therefore, this observable or quantities derived thereof
do not lend themselves to compare to the relaxation behaviour seen in
$\overline{\mathscr{I}}(t)$. Due to the apparent lack of an obvious scale
derived from an equilibrium quantity, it is not inconceivable that the scaling
factor of the ordinate, $f_{L}(W)$ (inset of Fig. 2), reflects an initial
condition, here the Néel state, and changes upon other choices.  The factor
$f_{L}(W)$ has been extracted from our data performing a standard scaling
analysis. It consists of rescaling the y-axis until the best data-collapse is
observed.
We mention that at late times when entanglement approaches its saturation
value, the master curve exhibits a large - possibly diverging - slope. The
implication is that entanglement propagates fast and may saturate long before
other physical observables do, which are subject to a local conservation law.

Figure 9: Approximate collapse to a single master curve of the traces for
$\overline{\mathscr{F}}(t)$ for different disorder values
$W{=}1.5,2.0,2.5,3.0,4.0,5.0,6.0$ displayed in Fig. 10 for $L=24$.
$\xi_{\mathrm{sp}}$ denotes the non-interacting localization length (further
details in Ref. [31]).  The $\overline{\mathscr{S}}_{\mathrm{e}}$ traces are
taken from Fig. 5 and Fig. 6. Inset: Prefactor $g(W)$ defined by the scaling
collapse as a function of $1/W^{2}$.

####  IV.2.3 Collapsing of density fluctuations – $\overline{\mathscr{F}}$
over $\overline{\mathscr{S}}_{\mathrm{e}}$

Traditional mode-coupling theory suggests that if
$\overline{\mathscr{S}}_{\mathrm{e}}$ is suitable as an internal clock for
$\overline{\mathscr{I}}$ then it should also be for related observables that
derive from the evolution of density modes. In this spirit, we show in Fig. 9
that also the time evolution of imbalance fluctuations (4) can be collapsed
employing the entanglement evolution, $\overline{\mathscr{S}}_{\mathrm{e}}$,
as an internal clock.  In appendix A.2 we offer preliminary evidence
suggesting that a collapse is also achieved for a third variable, i.e. the
density auto-correlation function.

Since we have already demonstrated the concept of a system-internal clock for
$\overline{\mathscr{I}}(t)$ in the regime of accelerated dynamics, $W\lesssim
3$, the excellent scaling observed for $\overline{\mathscr{F}}$ in this regime
is comforting but not, perhaps, surprising. Interestingly, a reasonably good
data collapse is seen for $\overline{\mathscr{F}}(t)$ even in the regime of
large disorder, $W\approx 8$, which we take as an additional encouragement for
the concept of internal clocks, here proposed.

Incidentally, it is implied here that at least with respect to the observable
density fluctuations (as seen within our observation time) there is no
indication of a phase transition all the way from moderate to strong disorder.

###  IV.3 Intermediate power-laws in average & typical
$\overline{\mathscr{F}}$

The raw data underlying the master curve is displayed in Fig. 10. The data
confirms a trend already observed in our previous work Nandy  _et al._ [50]:
the imbalance fluctuations exhibit a rather benign convergence with the system
size. We, therefore, are confident to identify for each trace an intermediate
regime between short times (plateau region) and longest times (pronounced
system-size dependency) that defines an (effective) power law

$\displaystyle\log\overline{\mathscr{F}}/\overline{\mathscr{F}}_{1}$
$\displaystyle=-(\rho_{\text{ave}}\ \xi_{\text{sp}}\log
2)\log\overline{\mathscr{S}}_{\mathrm{e}}$ (8)

with $\overline{\mathscr{F}}_{1}$ denoting the prefactor; $\xi_{\text{sp}}$ is
the non-interacting localization length extracted in the long time limit of
the second moment of the density-density correlation function at infinite
temperature [31].

Figure 10: Ensemble-averaged temporal fluctuations of the imbalance
${\mathscr{I}}(t)$ as defined in (4). Following our earlier work [50], we have
dressed $\overline{\mathscr{F}}$ with a power $1/\xi_{\text{sp}}\log 2$. For
selected traces, two system sizes,
$L=20~{}(\mathrm{lines}),24~{}(\mathrm{symbols})$, are given to expose finite-
size effects.

The exponent function $\rho_{\text{ave}}(W)$ is extracted fitting the data
Fig. 10 to the form (8); an analogous fitting procedure can also be performed
for the typical traces of ${\mathscr{F}}$ and $\mathscr{S}_{\mathrm{e}}$
yielding $\rho_{\text{typ}}(W)$.

The results for both exponents are displayed in Fig. 11. We first note that
the evolution of the traces with the disorder is rather smooth and slow; there
is no indication of a nearby MBL-transition. 555 One might ask why we expect
$\rho(W)$ to exhibit a signature of the transition into the proper MBL phase:
the definition of the exponent,
$\overline{\mathscr{F}}\propto\overline{\mathscr{S}}_{\mathrm{e}}^{\rho}$, at
least partially reflects the fact that $\overline{\mathscr{F}}(t)$ and
$\overline{\mathscr{S}}_{\mathrm{e}}(t)$ exhibit a power-law growth at long
times (in a large-enough system). At least for the entanglement entropy,
$\overline{\mathscr{S}}_{\mathrm{e}}(t)\propto t^{1/z_{\text{ee}}}$, the
growth becomes slower than any power, $1/z_{\text{ee}}\to 0$, so that
generically the proportionality
$\overline{\mathscr{F}}\propto\overline{\mathscr{S}}_{\mathrm{e}}^{\rho}$
should cease to hold.

Remarkably, the exponent functions display opposing trends with increasing
disorder. The typical value follows the naive expectation, which is that with
the increasing disorder, the interaction-induced damping of temporal
fluctuations is less effective, so $\rho_{\text{typ}}$ has a tendency to
decrease. While the typical exponent $\rho_{\text{typ}}(W)$ decreases, the
average exponent $\rho_{\text{ave}}(W)$ increases with disorder and there is,
indeed, a common intersection point $W_{{\mathscr{F}}}\approx 6$. At weaker
disorder, $W{<}W_{{\mathscr{F}}}$, rare samples exist that exhibit rather long
decoherence times, so in this regime $\rho_{\text{typ}}{>}\rho_{\text{ave}}$.
In the other regime, $W{>}W_{\overline{\mathscr{F}}}$ the situation is
reversed and the average is dominated by rare samples in which damping is
unusually effective. Notice that the crossover value is situated near the
regime in which the system dynamics changes from accelerated to decelerated,
$W_{\overline{\mathscr{F}}}\approx W_{c}$, cf. Fig. 1.

Summarizing, we interpret Fig. 11 in conjunction with Fig. 6 as supplying
fresh support in favor of the phase-diagram Fig. 1 and its main claim: In the
range $W_{c},W_{\overline{\mathscr{F}}}\approx 3-6$ there is a crossover
between two markedly different thermalizing regimes; there is no evidence that
the $t-V$-model exhibits a phase transition at disorder values below
$W\lesssim 10$. We thus advocate the point of view that until recently this
crossover has been widely misinterpreted in numerical work as indicating an
ergodicity-breaking (i.e. MBL-) transition, e.g. in Refs. [67, 68, 69, 6, 25,
70, 71, 26, 28, 72, 73, 29].

Figure 11: (Effective) exponents characterizing the evolution of the imbalance
fluctuations with the entanglement entropy,
${\mathscr{F}}{\propto}\mathscr{S}_{\mathrm{e}}^{\rho}$, for typical values,
$\rho_{\text{typ}}$, and average values $\rho_{\text{ave}}$. The intersection
point indicates the crossover from accelerated to decelerated dynamics, see
Fig. 1.

Recent numerical work claims that the MBL transition if it exists, occurs at a
disorder strength $W^{*}\gtrsim 20$ [34, 35]. We notice that already at
disorder value $W{\approx}3$ the non-interacting localization length
$\xi_{\mathrm{sp}}$ is of the order the lattice spacing $a$, which implies
$\xi_{\mathrm{sp}}(W^{*})\ll a$. In this sense, if the MBL transition occurs
at all, it for sure takes place at an extremely strong disorder. This
immediately prompts the question of its physical significance.

###  IV.4 Sample-to-sample fluctuations and thermalization

####  IV.4.1 Distribution of exponents

The importance of sample-to-sample fluctuations is demonstrated in Fig. 7, lhs
column. The plot shows the distribution of ${\mathscr{I}}$ being
logarithmically broad at large observation times already at moderate disorder,
$W{=}1.5$. In fact, a fraction of samples exhibit traces ${\mathscr{I}}(t)$
that does not indicate a discernible trend towards equilibration,
${\mathscr{I}}(t){\to}0$ in the limit $t\to\infty$, at all. This point is also
illustrated in Fig. 7, second column. It shows an effective exponent

$\displaystyle\beta(t)\coloneqq\frac{\partial\log{\mathscr{I}}(t)}{\partial\log
t}$ (9)

that characterizes the time-evolution of ${\mathscr{I}}(t)$ for an individual
sample; the corresponding (logarithmic) distribution functions
${\mathscr{P}}(\log\beta)$ are given in Fig. 12. As is clearly illustrated by
the data, the exponent distribution is logarithmically wide. Extremely strong
sample-to-sample fluctuations have been observed also by other authors. As an
interesting example, we mention Doggen  _et al._ [27]: these authors use a
machine-learning algorithm in order to distinguish localized from thermal
samples. Within a regime of intermediate disorder values and system sizes the
algorithm identifies the existence of two different sample types: those with
insulating character and manifestly thermalizing others.

Figure 12: Distribution function ${\mathscr{P}}$ of the effective exponents
$\beta(t)$ (Eq. (9)) that characterize the long-term behavior of the imbalance
${\mathscr{I}}(t)$ after a quench in systems of sizes $L=16,20,24$ and for two
disorder values in the (sub-)diffusive regime.

####  IV.4.2 Restoration of self-averaging

The tails of ${\mathscr{P}}$ are seen to be shrinking with increasing system
size in Fig. 12, if only exceedingly slowly. The general expectation is that
at $W<W^{*}$ the imbalance, ${\mathscr{I}}$, in interacting, disordered wires
is self-averaging. The statement implies that in the limit of large systems,
$L\gg 32$, the distribution ${\mathscr{P}}$ acquires zero width, eventually.
In other words, in the thermodynamic limit all samples, except for a small
number of measure zero, are expected to thermalize with a relaxation behavior
that is characterized by the same ($W$-dependent) exponents.

Due to computational limitations in system sizes and observation times, we can
neither confirm nor dispute this expectation. If it is true, then all traces
shown in Fig. 11 will undergo a slow evolution so that they eventually
collapse in the thermodynamic limit. However, also in this case it should not
go unnoticed that for the range of system sizes traditionally studied in
numerics, e.g. $L\lesssim 32$, individual samples hardly ever show the typical
behavior because the exponent distribution ${\mathscr{P}}$ is so wide.

###  IV.5 Relation to previous work: creep and RG

For further discussion, we recall the RG-scenario outlined in the
introduction. It assumes that samples fall into roughly two categories,
thermalizing (ergodic) “bubbles” and non-thermalizing “grains”. Fig. 7, lhs
column supports this picture and allows rough quantitative estimates: the
thermalizing fraction of samples comprises seven or eight out of ten samples
at $W{=}1.5$, so ${\mathfrak{q}}(L{=}24,W{=}1.5){\approx}80$%; this fraction
is rapidly decreasing for a larger disorder, e.g., we have only
${\mathfrak{q}}(L{=}24,W{=}2.5){\approx}50$% at $W{=}2.5$ and
${\mathfrak{q}}(L{=}24,W{=}5.0){\approx}0$.

One can consider growing the sample size, e.g., by merging two samples of
$L=24$, each. The RG assumes that the bubble and grain scenario continues to
hold after merging and then predicts the evolution of ${\mathfrak{q}}(L,W)$
with $L$ by postulating phenomenological growths rules [14, 19, 21, 20, 22,
23] Of interest to us is the general picture of the relaxation dynamics that
take place immediately after merging two samples. Specifically, at moderate to
strong disorder, $W\gtrsim 2.5$, most samples behave like grains since
${\mathfrak{q}}<1/2$. Hence, even when doubling the sample size the most
likely outcome is that a grain is added to a grain, so the relaxation dynamics
after merging is extremely slow and does not support thermalization. Only upon
growing the system size even further a bubble will be added eventually, which
can foster the delocalization process. Since a rather small bubble will need
to equilibrate a large grain, thermalization is especially slow. An
exponentially slow equilibration process that is enhanced by growing the
system size is the hallmark of “creep” [31, 30]. In this sense, our numerical
observations of large sample-to-sample fluctuations seen in Fig. 7, the slow
flow seen in Fig. 12, and creep are all qualitatively consistent with the
simplified RG-picture of thermalization through bubbles and grains. 666 In
recent numerical works, attempts have been made to substantiate the picture
and identify correlated/entangled clusters in disordered samples [71, 74, 75].
Since the data analysis in these papers underlies the assumption of a critical
point near $W\approx 3.8$, it remains to be seen to what extent the
conclusions are affected once creep is taken into account.

The alert reader will have realized that we adopted from the RG-concept its
intuitive phenomenological building blocks, i.e. bubbles and grain, and left
aside an important result of the reported RG-flow, which predicts a critical
point separating a localized from an ergodic phase [22, 20, 23]. We have
neglected this part of the RG studies, because, within the parameter regimes
investigated, we don’t see support for a critical scenario provided by “ab-
initio” simulations of microscopic models. We mention two possible reasons for
the discrepancy - other than reentrance behavior at a larger system size
outside of our observation window. (i) The critical point might be situated at
extremely large disorder values, as suggested by Morningstar  _et al._ [34]
and Sels [35]. (ii) The phenomenological RG-equations are oversimplifications
that miss relevant terms with a delocalizing effect.

##  V Conclusion

Short quantum wires at relatively weak interactions and intermediate disorder
strength can exhibit localization behavior: initial deviations from the
equilibrium density are highly resilient against equilibration; a non-
vanishing fraction of them may not, in fact, equilibrate at all. Understanding
the fate of shorter samples with respect to their relaxation dynamics when
successively growing their length is the central theme of many-body (de-)
localization (MBL). The numerical investigation we have presented in this work
motivates three statements:

(i) We follow how the relaxation dynamics of the sublattice imbalance
$\overline{\mathscr{I}}$ and its fluctuations $\overline{\mathscr{F}}$ evolve
over a large range of disorder values, reaching from below the clean bandwidth
to four times its value. The temporal flow of the ensemble dynamics is
conveniently parameterized by $\overline{\mathscr{S}}_{\mathrm{e}}$, which
acts as a model for an internal (ensemble) clock. The usefulness of the
“internal clock” concept becomes most apparent by demonstrating that for a
broad range of disorder values time traces $\overline{\mathscr{I}}(t)$ and
$\overline{\mathscr{S}}_{\mathrm{e}}(t)$ can be collapsed to a single master
curve - for fixed system size. The collapse works within the entire window of
investigated disorder values, which implies the absence of a localization
transition even when the disorder exceeds the clean bandwidth by a factor of
four.

(ii) We observe extremely strong sample-to-sample fluctuations: For systems
with $L=24$ sites and moderate disorder non-ergodic samples (grains) coexist
with highly ergodic samples (bubbles). Previous work has reported an extremely
slow flow toward equilibration when growing the system size (“creep”) [30]. In
this work we have argued that the creep phenomenon is closely related to the
strong sample-to-sample fluctuations here observed; the relationship has been
established by borrowing ideas from a real-space renormalization group
approach and the avalanche concept [18].

(iii) While the existing evidence points to a thermodynamic limit that
represents an equilibrating thermal phase, the flow towards this limit
suggests the existence of subphases. Specifically, when the disorder strength
reaches about two times the bandwidth, accelerated dynamics gives way to
decelerated dynamics as indicated in Fig. 1. In this work, we present
additional numerical evidence for this scenario, based on the observation that
the evolution of ${\mathscr{F}}(t)$ exhibits two regimes: At weak disorder,
typical fluctuations are damped more strongly than average fluctuations; at
larger disorder, it is the other way round.

##### Outlook.

As an outlook, we express our belief that the neglect of creep in many, if not
most, of the earlier computational studies of MBL, invalidates in parts the
data interpretation that has been offered in these works, presumably often in
crucial ways. At the time being it is too early writing an MBL review from the
perspective of many-body de-localization (MBdL), i.e. “creep”. Nonetheless, a
firm ground has been laid that allows for a more careful analysis of the
physical phenomena to be encountered. Many fascinating ideas have been
expressed in the past, some of them formulated as mathematical theorems, some
of them in terms of toy models, all of them making predictions that merit a
careful computational test. Of course, creep will have to be included as a
hallmark of delocalization physics in the data analysis - and no longer be
ignored.

The most important statement supporting the existence of MBL goes back to
Imbrie [76, 16]. While Imbrie’s proof is believed to be rigorous, it also
relies upon assumptions, e.g., concerning the spectral statistics of the
sample Hamiltonian. It will be interesting to see in future work, whether the
lack of evidence for MBL in the XXZ-Heisenberg model is due to the disorder
being still too weak, due to the proof being not fully complete yet, or due to
a much simpler reason, which is that the theorem does not apply to the XXZ-
chain, because its assumptions are not met.

Another open pressing question concerns the nature of MBdL in correlated
disorder, such as represented, e.g., by the Andre-Aubry potential (AA). In our
previous studies Weiner  _et al._ [30] of charge-density relaxations, we did
not detect a qualitative difference between correlated and fully uncorrelated
randomness  (also see preliminary data in App. A.3); similarly, also here we
have no indication that the dominating physics is related to the effects of
rare regions. Both observations together prompt the expectation that the AA-
model exhibits creep, i.e., MBdL. If true, the interpretation of cold-atom
experiments in terms of the observation of MBL proper is challenged [77, 78,
79].

One more crucial issue that merits closer scrutiny relates to the choice of
the initial state used for time propagation, e.g., a Neél state. In
thermalizing phases, it is tempting to assume that the qualitative dynamics
seen in a quench mostly reflect properties of the phase rather than the
initial state;  confirming preliminary evidence is given in App. A.1. To what
extent this remains true also in marginally thermalizing situations requires
further study.

We conclude with a caveat: In this article, we have adopted a manner of
speaking that is established in the MBL community according to which a phase
is “equilibrating” or “thermal” if the temporal evolution of the entanglement
entropy is at large times faster than $\ln t$ and if simultaneously the
sublattice imbalance and derived quantities in the thermodynamic limit decay
to zero. It should be noted that the traditional meaning of equilibration
implies that all local observables relax towards their equilibrium value,
eventually. Whether the equilibrating/thermalizing phase(s) indicated in Fig.
1 also satisfies such a stronger condition is a matter of ongoing research.

##  VI ACKNOWLEDGMENTS

We would like to thank J. Bardarson, I. Gornyi, A. Mirlin, M. Kiefer-
Emmanouilidis, J. Sirker, and J. Zakrzewski for critical reading of the
manuscript, and for their valuable comments, which improved our manuscript
considerably. SB would like to thank S. Nandy for several discussions, and for
an earlier collaboration on a closely related topic. FE expresses his
gratitude to A. Rosch for an enjoyable set of conversations on the topic. SB
acknowledges support from SERB-DST, India, through Matrics (No.
MTR/2019/000566), and MPG for funding through the Max Planck Partner Group at
IITB. Funding from the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) through EV30/11-2, EV30/12-1, EV30/14-1, EV30/14-2 is gratefully
acknowledged. IM acknowledges financial support from Prime Minister’s Research
Fellows (PMRF) scheme offered by the Ministry of Education, Government of
India.

##  Appendix A Further confirmation of the generality of the internal-clock
concept

In this section we provide additional evidence for the general applicability
of the concept of the (ensemble) internal clock for synchronizing time traces
obtained in different disorder regimes.

###  A.1 Different initial condition: Neel-block states

Figure 13:  Time evolution of the imbalance, Eq. 10, after quenching a Neél
block state, $|\uparrow\uparrow\downarrow\downarrow\ldots\rangle$ for two
system sizes $L=16,20$. The figure is fully analogous to Fig. 2, which is for
quenching a Neél state. In both cases the data collapse is seen to be of a
similar quality, underlining the generality of the internal clock as a means
to synchronize time evolution in different ensembles. (Number of samples 200
($W<4$) and 600 ($W>4$).

Synchronizing by using the evolving entanglement entropy as a measure of time
also works for initial states different from the Neél state. We demonstrate
this by quenching from a Neél-block state
$|\uparrow\uparrow\downarrow\downarrow\ldots\rangle$, which similar to the
Neél state is situated in the middle of the many-body spectrum. For the block
state, we define the imbalance in the following way:

$\overline{\mathscr{I}}(t)=\sum_{j=1}^{L}\ n_{4j-3}(t)+\ n_{4j-2}(t)-\
n_{4j-1}(t)-\ n_{4j}(t)$ (10)

After quenching Neél-block states, we observe an approximate collapse of all
traces for different disorder values, see Fig. 13, similar to the case of the
Neél state, Fig. 2.

###  A.2 Third observable: Correlation function

Figure 14:  Time evolution of the nearest, $\overline{\mathscr{C}}_{1}(t)$
(upper row), and next nearest, $\overline{\mathscr{C}}_{1}(t)$ (lower row),
density correlations for system size $L=20$. Similar to the imbalance, Fig. 2,
and the density fluctuations, Fig. 9, after employing
$\overline{\mathscr{S}}_{\mathrm{e}}(t)$ as an internal clock and rescaling
the ordinate by a scale factor $\tilde{f}_{1,2}(W)$ a (reasonably) good
scaling collapse is achieved, see rhs column. Fitting parameters:
$\tilde{f}_{1}(W)=[1,1,1,1,1,1.07]$ and
$\tilde{f}_{2}(W)=[1,0.967,0.96,0.977,1.05,1.2]$ for $W$ ranging from 2.0 to
5.0. (Number of samples: 200 at $W<4$ and 600 at $W>4$.)

In order to demonstrate that $\overline{\mathscr{S}}_{\mathrm{e}}(t)$
potentially synchronizes a large set of observables, we here give a third
synchronizing example, in addition to the average imbalance
$\overline{\mathscr{I}}(t)$ and density fluctuations
$\overline{\mathscr{F}}(t)$ considered in the main text. Specifically, we
consider the density auto-correlation function defined as

${\mathscr{C}_{l}}(t)=\frac{1}{L-l}\sum_{i=1}^{L-l}C_{i,i+l}(t),$ (11)

with
$C_{i,i+l}(t)=\langle\hat{n}_{i}\hat{n}_{i+l}\rangle-\langle\hat{n}_{i}\rangle\langle\hat{n}_{i+l}\rangle;$
$\langle\ldots\rangle$ represents expectation value in the state
$|\psi(t)\rangle=e^{-iHt/\hbar}|\psi(0)\rangle$. We consider ensemble
averages, which restore the translational invariance:
$\overline{\mathscr{C}}_{l}(t)$ is the ensemble averaged correlator where the
overline stands for the ensemble average.

Fig. 14 shows our results for nearest-neighbour and next nearest-neighbour
correlations, i.e. $\overline{\mathscr{C}}_{1}(t)$ and
$\overline{\mathscr{C}}_{2}(t)$. As one infers from the right-hand side
column, a reasonable collapse is also achieved for this observable when using
$\overline{\mathscr{S}}_{\mathrm{e}}(t)$ as an effective ensemble time.

###  A.3 Quasi-periodic potential

As a sanity check, we here demonstrate that synchronization with
$\overline{\mathscr{S}}_{\mathrm{e}}(t)$ also works in situations of
correlated randomness. To this end we adopt the Aubry-Andre model which is
given by eq.(1) with $\epsilon_{i}=\frac{W}{2}\cos({2\pi\sigma i+\phi})$,
where, $\sigma=\frac{\sqrt{5}-1}{2}$ is the golden ratio and $\phi$ is a
random number chosen uniformly from the range 0 to $2\pi$. As seen in Fig. 15,
the quasi-periodic model shows a collapse of the imbalance trances similar to
the situation with uncorrelated disorder model, Fig. 2. This strengthens the
general applicability of the concept of a synchronizing interal clock.

Figure 15:  Evolution of ensemble-averaged particle imbalance after quenching
Nèel state at system size $L=16,20$ for the case of correlated disorder
(Aubry-Andre model). The figure is fully analogous to Fig. 2, which features
uncorrelated disorder. Inset: Scaling factor $f_{L}(W)$ as a function of $W$.
(Number of disorder realizations: 600).

##  Appendix B Fitted parameters

The list of the parameters obtained from fitting power-law corrections, (7) to
the data given in Fig. 5.

$W$ $L$ | $\overline{\mathscr{S}}_{\mathrm{e}}^{\infty}$ | $\overline{\mathscr{S}}_{\mathrm{e}}^{\infty}/\mathscr{S}_{\mathrm{Page}}$ | $c_{0}$ | $\gamma_{L}(W)$
---|---|---|---|---
0.0 | 16 | 3.568(4) | 0.707 | - | -
| 20 | 4.504(21) | 0.700 | - | -
| 24 | 5.453(12) | 0.697 | - | -
| 26 | 5.929(6) | 0.696 | - | -
1.25 | 16 | 3.90(1) | 0.773 | - | -
| 20 | 5.118(17) | 0.796 | - | -
| 24 | 6.253(20) | 0.800 | - | -
| 26 | 6.848(24) | 0.805 | - | -
1.5 | 16 | 3.871(3) | 0.767 | - | -
| 20 | 5.158(5) | 0.802 | - | -
| 24 | 6.387(8) | 0.817 | 83(10) | 0.852(40)
| 26 | 7.012(8) | 0.824 | 74(4) | 0.796(23)
2.0 | 16 | 3.524(7) | 0.698 | 40(10) | 0.807(48)
| 20 | 5.117(10) | 0.796 | 13.96(30) | 0.445(54)
| 24 | 6.645(30) | 0.850 | 18.66(24) | 0.40(9)
| 26 | 7.476(9) | 0.878 | 19.68(13) | 0.369(3)
2.5 | 16 | 2.928(1) | 0.580 | 19.32(7) | 0.653(1)
| 20 | 4.450(5) | 0.692 | 10.95(10) | 0.362(4)
| 24 | 6.307(5) | 0.807 | 13.13(4) | 0.269(1)
| 26 | 7.372(10) | 0.866 | 14.58(8) | 0.239(2)
3.0 | 16 | 2.371(4) | 0.470 | 10.75(80) | 0.545(8)
| 20 | 3.497(11) | 0.544 | 9.61(41) | 0.356(7)
| 24 | 5.187(10) | 0.663 | 9.96(19) | 0.226(7)
| 26 | 6.70(20) | 0.788 | 10.87(17) | 0.166(10)
Table 1:  Parameters obtained from fitting the expression (7) for
$\overline{\mathscr{S}}_{\mathrm{e}}(t)$ to the data given in Fig. 5a). Fits
of the clean limit, Fig. 3 right column, have also been included.

Details of the procedure are described in the main text; the ratio
$\overline{\mathscr{S}}_{\mathrm{e}}^{\infty}/{\mathscr{S}}_{\mathrm{Page}}$
is plotted in Fig. 5(d). At weak disorder, the saturation values have been
read directly from the original traces. At stronger disorder, $W>3.0$ the
saturation behavior moves out of the time window where the ansatz (7) would
apply.

To get a reasonable estimate of the error in the fitting parameters, we use
the bootstrap method [80, 81]. This involves generating 100 new random sample
traces around the original trace with a variance that is given by the sampling
error. All these traces are then fitted with Eq. 7 to get a spread with $95\%$
confidence interval in the fitting parameters
$\overline{\mathscr{S}}_{\mathrm{e}}^{\infty},\gamma_{L}(W)$.

Figure 16: The time evolution of the (second) Renyi entropy in units of
$\overline{\mathscr{S}}_{\mathrm{e}}$ in a broad range of disorder values, and
system sizes $L=16,20,24$.

##  Appendix C $\mathscr{S}_{2}$ vs $\mathscr{S}_{\mathrm{e}}$

To address the issue we depict in Fig. 16 the ratio of
$\overline{{\mathscr{S}}}_{2}/\overline{\mathscr{S}}_{\mathrm{e}}$. The ratio
deviates from (roughly) a constant at short times and at long times, where
system-size effects enter. We conclude that the entropy model for the internal
clock has the best chance to work in an intermediate time regime where the
ratio adopts a plateau-type behavior. Finite-size effects set in as soon as
$\overline{\mathscr{S}}_{\mathrm{e}}(t)\propto{\mathscr{S}}_{\mathrm{Page}}$
with a prefactor that shrinks with increasing $W$.

As a model for the internal clock, we have adopted the entanglement entropy
$\mathscr{S}_{\mathrm{e}}(t)$. Alternative choices are available, such as the
generalized Renyi-entropies

$\displaystyle{\mathscr{S}}^{(\alpha)}\coloneqq(1-\alpha)^{-1}\ln\text{Tr}_{\mathrm{A}}\hat{\rho}_{\mathrm{A}}^{\alpha}$
(12)

with $\mathscr{S}_{2}={\mathscr{S}}^{(2)}$ and
$\mathscr{S}_{\mathrm{e}}={\mathscr{S}}^{(1)}$. Therefore, a question arises
concerning the equivalence of different choices.

Figure 17: Shows the time evolution of entanglement entropy and imbalance for
both the exact and the Chebyshev method for two randomly chosen disorder
realizations (red and blue) for $L=16$, and $W=1.5$.

##  Appendix D Convergence of time evolution

The validity of the conclusions presented in this work relies upon the
reliability of our simulation data.

To demonstrate the accuracy of the Chebyshev expansion as we have implemented
it, we have performed a per-sample comparison with data from exact
diagonalization (ED). Figure 17 displays the result for two typical, i.e.
randomly chosen samples. As seen there, within our simulation time $t\sim
10^{3}$ the traces are indistinguishable. The main reason for the usage of the
Chebyshev expansion is that it allows treating systems larger than those
affordable with ED, i.e. $L\gtrsim 16$ .

## References

  * Basko  _et al._ [2006] D. M. Basko, I. L. Aleiner, and B. L. Altshuler, Metal insulator transition in a weakly interacting many electron system with localized single particle states, Ann. Phys. 321, 1126 (2006).
  * Gornyi  _et al._ [2005] I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Interacting electrons in disordered wires: Anderson localization and low-$t$ transport, Phys. Rev. Lett. 95, 206603 (2005).
  * Oganesyan and Huse [2007] V. Oganesyan and D. A. Huse, Localization of interacting fermions at high temperature, Phys. Rev. B 75, 155111 (2007).
  * Nandkishore and Huse [2015] R. Nandkishore and D. A. Huse, Many-body localization and thermalization in quantum statistical mechanics, Ann. Rev. Condens. Matter Phys. 6, 15 (2015).
  * Bera  _et al._ [2015] S. Bera, H. Schomerus, F. Heidrich-Meisner, and J. H. Bardarson, Many-body localization characterized from a one-particle perspective, Phys. Rev. Lett. 115, 046603 (2015).
  * Devakul and Singh [2015] T. Devakul and R. R. P. Singh, Early breakdown of area-law entanglement at the many-body delocalization transition, Phys. Rev. Lett. 115, 187201 (2015).
  * Prelovšek  _et al._ [2017] P. Prelovšek, M. Mierzejewski, O. Barišić, and J. Herbrych, Density correlations and transport in models of many-body localization, Ann. Phys. (Berlin) 529, 1600362 (2017).
  * Abanin  _et al._ [2019] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Colloquium: Many-body localization, thermalization, and entanglement, Rev. Mod. Phys. 91, 021001 (2019).
  * Alet and Laflorencie [2018] F. Alet and N. Laflorencie, Many-body localization: An introduction and selected topics, C. R. Phys. (2018).
  * Binder and Young [1986] K. Binder and A. P. Young, Spin glasses: Experimental facts, theoretical concepts, and open questions, Rev. Mod. Phys. 58, 801 (1986).
  * Rieger and Young [1997] H. Rieger and A. P. Young, Quantum spin glasses, in  _Complex Behaviour of Glassy Systems_ , edited by M. Rubí and C. Pérez-Vicente (Springer Berlin Heidelberg, Berlin, Heidelberg, 1997) pp. 256–265.
  * Serbyn  _et al._ [2013a] M. Serbyn, Z. Papić, and D. A. Abanin, Local conservation laws and the structure of the many-body localized states, Phys. Rev. Lett. 111, 127201 (2013a).
  * Huse  _et al._ [2014] D. A. Huse, R. Nandkishore, and V. Oganesyan, Phenomenology of fully many-body-localized systems, Phys. Rev. B 90, 174202 (2014).
  * Vosk  _et al._ [2015] R. Vosk, D. A. Huse, and E. Altman, Theory of the many-body localization transition in one-dimensional systems, Phys. Rev. X 5, 031032 (2015).
  * O’Brien  _et al._ [2016] T. E. O’Brien, D. A. Abanin, G. Vidal, and Z. Papić, Explicit construction of local conserved operators in disordered many-body systems, Phys. Rev. B 94, 144208 (2016).
  * Imbrie  _et al._ [2017] J. Z. Imbrie, V. Ros, and A. Scardicchio, Local integrals of motion in many-body localized systems, Ann. Phys. (Berlin) 529, 1600278 (2017).
  * De Roeck and Huveneers [2017] W. De Roeck and F. m. c. Huveneers, Stability and instability towards delocalization in many-body localization systems, Phys. Rev. B 95, 155129 (2017).
  * Thiery  _et al._ [2018] T. Thiery, F. Huveneers, M. Müller, and W. De Roeck, Many-body delocalization as a quantum avalanche, Phys. Rev. Lett. 121, 140601 (2018).
  * Zhang  _et al._ [2016] L. Zhang, B. Zhao, T. Devakul, and D. A. Huse, Many-body localization phase transition: A simplified strong-randomness approximate renormalization group, Phys. Rev. B 93, 224201 (2016).
  * Dumitrescu  _et al._ [2019] P. T. Dumitrescu, A. Goremykina, S. A. Parameswaran, M. Serbyn, and R. Vasseur, Kosterlitz-thouless scaling at many-body localization phase transitions, Phys. Rev. B 99, 094205 (2019).
  * Dumitrescu  _et al._ [2017] P. T. Dumitrescu, R. Vasseur, and A. C. Potter, Scaling theory of entanglement at the many-body localization transition, Phys. Rev. Lett. 119, 110604 (2017).
  * Goremykina  _et al._ [2019] A. Goremykina, R. Vasseur, and M. Serbyn, Analytically solvable renormalization group for the many-body localization transition, Phys. Rev. Lett. 122, 040601 (2019).
  * Morningstar and Huse [2019] A. Morningstar and D. A. Huse, Renormalization-group study of the many-body localization transition in one dimension, Phys. Rev. B 99, 224205 (2019).
  * Morningstar  _et al._ [2020] A. Morningstar, D. A. Huse, and J. Z. Imbrie, Many-body localization near the critical point, Phys. Rev. B 102, 125134 (2020).
  * 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).
  * Macé  _et al._ [2019] N. Macé, F. Alet, and N. Laflorencie, Multifractal scalings across the many-body localization transition, Phys. Rev. Lett. 123, 180601 (2019).
  * 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).
  * Sierant  _et al._ [2020a] P. Sierant, M. Lewenstein, and J. Zakrzewski, Polynomially filtered exact diagonalization approach to many-body localization, Phys. Rev. Lett. 125, 156601 (2020a).
  * Doggen  _et al._ [2021] E. V. Doggen, I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Many-body localization in large systems: Matrix-product-state approach, Ann. Phys. 435, 168437 (2021).
  * Weiner  _et al._ [2019] F. Weiner, F. Evers, and S. Bera, Slow dynamics and strong finite-size effects in many-body localization with random and quasiperiodic potentials, Phys. Rev. B 100, 104204 (2019).
  * Bera  _et al._ [2017] 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. 118, 196801 (2017).
  * Lezama  _et al._ [2019] T. L. M. Lezama, S. Bera, and J. H. Bardarson, Apparent slow dynamics in the ergodic phase of a driven many-body localized system without extensive conserved quantities, Phys. Rev. B 99, 161106 (2019).
  * Sierant and Zakrzewski [2022] P. Sierant and J. Zakrzewski, Challenges to observation of many-body localization, Phys. Rev. B 105, 224203 (2022).
  * Morningstar  _et al._ [2022] A. Morningstar, L. Colmenarez, V. Khemani, D. J. Luitz, and D. A. Huse, Avalanches and many-body resonances in many-body localized systems, Phys. Rev. B 105, 174205 (2022).
  * Sels [2022] D. Sels, Bath-induced delocalization in interacting disordered spin chains, Phys. Rev. B 106, L020202 (2022).
  * Kiefer-Emmanouilidis  _et al._ [2020] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Evidence for unbounded growth of the number entropy in many-body localized phases, Phys. Rev. Lett. 124, 243601 (2020).
  * Kiefer-Emmanouilidis  _et al._ [2021a] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Slow delocalization of particles in many-body localized phases, Phys. Rev. B 103, 024203 (2021a).
  * Kiefer-Emmanouilidis  _et al._ [2021b] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Unlimited growth of particle fluctuations in many-body localized phases, Annals of Physics 435, 168481 (2021b).
  * Kiefer-Emmanouilidis  _et al._ [2022a] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Particle fluctuations and the failure of simple effective models for many-body localized phases, SciPost Phys. 12, 034 (2022a).
  * Luitz and Lev [2020] D. J. Luitz and Y. B. Lev, Absence of slow particle transport in the many-body localized phase, Phys. Rev. B 102, 100202 (2020).
  * Ghosh and Žnidarič [2022] R. Ghosh and M. Žnidarič, Resonance-induced growth of number entropy in strongly disordered systems, Phys. Rev. B 105, 30 (2022).
  * Kiefer-Emmanouilidis  _et al._ [2022b] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Comment on “Resonance-induced growth of number entropy in strongly disordered systems”, arXiv:2203.06689 (2022b).
  * Ghosh and Žnidarič [2022] R. Ghosh and M. Žnidarič, Response to ‘Comment on Theory of growth of number entropy in disordered systems’, arXiv:2205.05768 (2022).
  * Šuntajs  _et al._ [2020a] J. Šuntajs, J. Bonča, T. Prosen, and L. Vidmar, Quantum chaos challenges many-body localization, Phys. Rev. E 102, 062144 (2020a).
  * Panda  _et al._ [2019] R. K. Panda, A. Scardicchio, M. Schulz, S. R. Taylor, and M. Žnidarič, Can we study the many-body localisation transition?, Epl 128, 10.1209/0295-5075/128/67003 (2019), arXiv:1911.07882 .
  * Sierant  _et al._ [2020b] P. Sierant, D. Delande, and J. Zakrzewski, Thouless time analysis of anderson and many-body localization transitions, Phys. Rev. Lett. 124, 186601 (2020b).
  * Šuntajs  _et al._ [2020b] J. Šuntajs, J. Bonča, T. Prosen, and L. Vidmar, Ergodicity breaking transition in finite disordered spin chains, Phys. Rev. B 102, 064207 (2020b).
  * Sels and Polkovnikov [2021] D. Sels and A. Polkovnikov, Dynamical obstruction to localization in a disordered spin chain, Phys. Rev. E 104, 054105 (2021).
  * Abanin  _et al._ [2021] 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 427, 168415 (2021).
  * Nandy  _et al._ [2021] S. Nandy, F. Evers, and S. Bera, Dephasing in strongly disordered interacting quantum wires, Phys. Rev. B 103, 085105 (2021).
  * Brenig [1989] W. Brenig,  _Statistical Theory of Heat: Nonequilibrium Phenomena_ (Springer-Verlag, 1989).
  * Weiße  _et al._ [2006] A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, The kernel polynomial method, Rev. Mod. Phys. 78, 275 (2006).
  * Lukin  _et al._ [2019] A. Lukin, M. Rispoli, R. Schittko, M. E. Tai, A. M. Kaufman, S. Choi, V. Khemani, J. Léonard, and M. Greiner, Probing entanglement in a many-body–localized system, Science 364, 256 (2019).
  * Brydges  _et al._ [2019] T. Brydges, A. Elben, P. Jurcevic, B. Vermersch, C. Maier, B. P. Lanyon, P. Zoller, R. Blatt, and C. F. Roos, Probing Renyi entanglement entropy via randomized measurements, Science 364, 260 (2019).
  * Calabrese and Cardy [2005] P. Calabrese and J. Cardy, Evolution of entanglement entropy in one-dimensional systems, J. Stat. Mech.: Theory Exp. 2005 (04), P04010.
  * Zhao  _et al._ [2016] Y. Zhao, F. Andraschko, and J. Sirker, Entanglement entropy of disordered quantum chains following a global quench, Phys. Rev. B 93, 205146 (2016).
  * Page [1993a] D. N. Page, Average entropy of a subsystem, Phys. Rev. Lett. 71, 1291 (1993a).
  * Kim and Huse [2013] H. Kim and D. A. Huse, Ballistic spreading of entanglement in a diffusive nonintegrable system, Phys. Rev. Lett. 111, 127205 (2013).
  * 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).
  * Luitz and Lev [2017] D. J. Luitz and Y. B. Lev, The ergodic side of the many-body localization transition, Ann. Phys. (Berlin) 529, 1600350 (2017).
  * Tikhonov  _et al._ [2016] K. S. Tikhonov, A. D. Mirlin, and M. A. Skvortsov, Anderson localization and ergodicity on random regular graphs, Phys. Rev. B 94, 220203 (2016).
  * Žnidarič  _et al._ [2008] M. Žnidarič, T. Prosen, and P. Prelovšek, Many-body localization in the heisenberg xxz magnet in a random field, Phys. Rev. B 77, 064426 (2008).
  * Bardarson  _et al._ [2012] J. H. Bardarson, F. Pollmann, and J. E. Moore, Unbounded growth of entanglement in models of many-body localization, Phys. Rev. Lett. 109, 017202 (2012).
  * Serbyn  _et al._ [2013b] M. Serbyn, Z. Papić, and D. A. Abanin, Universal slow growth of entanglement in interacting strongly disordered systems, Phys. Rev. Lett. 110, 260601 (2013b).
  * Page [1993b] D. N. Page, Average entropy of a subsystem, Phys. Rev. Lett. 71, 1291 (1993b).
  * Singh  _et al._ [2016] R. Singh, J. H. Bardarson, and F. Pollmann, Signatures of the many-body localization transition in the dynamics of entanglement and bipartite fluctuations, New. J. Phys. 18, 023046 (2016).
  * 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).
  * Bar Lev and Reichman [2014] Y. Bar Lev and D. R. Reichman, Dynamics of many-body localization, Phys. Rev. B 89, 220201 (2014).
  * Khemani  _et al._ [2017] V. Khemani, S. P. Lim, D. N. Sheng, and D. A. Huse, Critical properties of the many-body localization transition, Phys. Rev. X 7, 021013 (2017).
  * Herviou  _et al._ [2019] L. Herviou, S. Bera, and J. H. Bardarson, Multiscale entanglement clusters at the many-body localization phase transition, Phys. Rev. B 99, 134205 (2019).
  * Laflorencie  _et al._ [2020] N. Laflorencie, G. Lemarié, and N. Macé, Chain breaking and kosterlitz-thouless scaling at the many-body localization transition in the random-field heisenberg spin chain, Phys. Rev. Research 2, 042033 (2020).
  * Chanda  _et al._ [2020] T. Chanda, P. Sierant, and J. Zakrzewski, Time dynamics with matrix product states: Many-body localization transition of large systems revisited, Phys. Rev. B 101, 035148 (2020).
  * Szołdra  _et al._ [2021] T. Szołdra, P. Sierant, K. Kottmann, M. Lewenstein, and J. Zakrzewski, Detecting ergodic bubbles at the crossover to many-body localization using neural networks, Phys. Rev. B 104, L140202 (2021).
  * Hémery  _et al._ [2022] K. Hémery, F. Pollmann, and A. Smith, Identifying correlation clusters in many-body localized systems, Phys. Rev. B 105, 064202 (2022).
  * Imbrie [2016] J. Z. Imbrie, Diagonalization and many-body localization for a disordered quantum spin chain, Phys. Rev. Lett. 117, 027201 (2016).
  * 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).
  * 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).
  * Rispoli  _et al._ [2019] M. Rispoli, A. Lukin, R. Schittko, S. Kim, M. E. Tai, J. Léonard, and M. Greiner, Quantum critical behaviour at the many-body localization transition, Nature 573, 385 (2019).
  * Young [2012] P. Young, Everything you wanted to know about data analysis and fitting but were afraid to ask, arXiv:1210.3781 (2012).
  * Press  _et al._ [2007] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,  _Numerical Recipes 3rd Edition: The Art of Scientific Computing_ (Cambridge University Press, London, 2007).