# Slow Dynamics and Kohlrausch Relaxation in Isolated Disordered Many-Body
Systems

Asmi Haldar  Max Planck Institute for the Physics of Complex Systems,
Nöthnitzer Straße 38, D-01187, Dresden, Germany

###### Abstract

The Kohlrausch(-Williams-Watts) law of stretched exponential relaxation has
been observed for more than a century and a half in diverse complex classical
systems. Here we show that this law describes relaxation quite generically in
closed (executing Schrödinger dynamics), interacting disordered many-body
systems across a range of system sizes using interaction range and disorder
strength as primary tuning parameters. This we observe for both time-
independent and periodically driven (Floquet) systems. Finite-size analysis
indicates the persistence of this non-thermal relaxation regime in the
thermodynamic limit thus defining a distinct dynamical regime. This regime
exhibits a peak in the time-scale of the perceptible relaxation, upon crossing
over from weak to strong disorder. We provide a simple picture of this
behavior, which naturally accounts for its general occurrence. Formation of
spin-glass – one of the possible mechanisms for stretched relaxation appears
incidental to the occurrence of Kohlrausch law in our context. Finally, we
provide a simple non-Hermitian Hamiltonian formulation for the dynamics of a
single spin embedded in the disordered chain. This provides an analytical
formula that captures not only the Kohlrausch relaxation of the disorder
averaged auto-correlation but also captures the largely diverse dynamics of an
arbitrary target spin in the system. Our work hence also provides a concrete
quantification of the “pre-thermal slowness” in many-body disordered system.

##  I Introduction

Strong Disorder is believed to result in absolute localization in many-
particle systems even in the presence of interactions - a phenomenon known as
many-body localization (MBL) [1, 2, 3, 4, 5, 6]. MBL is believed to persist
even in the presence of periodic drive and can exhibit interesting new phases,
collectively known as Floquet MBL phases [7, 8]. A striking example of Floquet
MBL phase is the discrete time crystalline (DTC) phase where the dynamics of a
periodically driven system breaks the discrete time-translational symmetry of
the time-periodic Hamiltonian (see e.g., [9, 10, 11] and references therein).
In general, the actual stability of an MBL (and hence a DTC) phase in the
thermodynamic limit is a matter of considerable debate [12, 13, 14, 15, 16,
17, 18, 19], though the stability of a DTC over experimentally relevant length
and time-scale has been established already [20, 21, 22, 23, 24].

Here, we leave this debate to the side and address a different aspect of a
Many-Body Disordered (MBD) system – what, if any qualitative difference in the
relaxation behaviour of local correlations occur as a consequence of disorder.
Our results not only provide a positive answer to this question, but also
uncovers a connection between the non-thermal dynamics of an MBD system, and a
broadly occurring yet somewhat enigmatic law that has been observed in the
relaxation of complex systems in the classical world over more than a century
and a half.

The relaxation dynamics of a many-body system may slow down qualitatively due
to disorder and inhomogeneity. In the classical regime, the mechanism of
slowing down are believed to be quite different in different cases, but a
surprisingly general law of relaxation applies under remarkably diverse
circumstances. This is known as the famous Kohlrausch law of stretched
exponential relaxation [25], where

${\rm
Auto~{}Correlation}~{}\sim\exp{\left\\{-\left(\frac{t}{\tau}\right)^{\gamma}\right\\}},$
(1)

where $0<\gamma\leq 1,$ with $\gamma=1$ giving the simple exponential
relaxation. Examples include slow dynamics of disordered classical systems
like structural glasses and super-cooled liquids [26], spin glasses [27] and
visco-elastic media [25]. The origin and ubiquity of the Kohlrausch law has
remained an intriguing open problem of physics and chemistry in spite of
powerful attempts to address it (see, e.g., [28, 29]). Here we show that the
Kohlrausch law also appears in the zero-temperature Schrödinger dynamics of
interacting quantum disordered many-body systems. The Kohlrausch form with its
stretching exponent ($\gamma$) and the time-scale ($\tau$) provides a
quantification of the well-known “slowness” of heating in the so called
“prethermal regime” [30, 31, 32, 33] in the context of localized quantum many-
body systems [22, 34, 24].

Figure 1: Kohlrausch law of stretched exponential relaxation of the disorder-
averaged auto-correlation $\overline{\tilde{C}}$ from least (bottom) to most
localized regime (top) for $L=26$. Representative cuts across the
multidimensional parameter space of the Floquet system (Eq. 2) are shown
(several others taken – all show similar behaviour; data not shown). Dotted
lines are numerical data and the solid lines are the fittings of the data with
a stretched exponential function (Eq. 7). From Left to Right, as a function of
(a) $\alpha$ (keeping $W=0.4\pi$, $\epsilon=0.2$ ), (b) $W$ (keeping
$\alpha=2.6,\epsilon=0.2$)and (c) $\epsilon$ (keeping $W=0.4\pi,\alpha=2.6$),
respectively.

We demonstrate the occurrence of the Kohlrausch law in largely diverse
settings - both with long and short range interactions, periodically driven
(Floquet) Hamiltonians and static ones. For concreteness, we primarily focus
on periodically driven quantum spin chains with power-law interactions, whose
range can be tuned by varying the exponent. A short-range version of the model
exhibits a stable DTC-MBL phase in appropriate parameter regimes in finite
systems [35]. For this system we characterize the phenomenological features of
the slow regime. Then we demonstrate those in other settings. We numerically
uncover distinct early and late dynamical regimes, and in none of those we
find any significant signature of the putative dynamical phase of “critical
DTC” [36, 10] with scale-free relaxation in our one-dimensional geometry.

Finally, we provide an analytical formulation based on a non-hermitian
Hamiltonian approach to capture the dynamics of an arbitrary target spin
embedded in the disordered system. The idea is to describe the localization
behaviour of the system in terms of its action as a bath on an arbitrary spin
embedded in it. While the approach is incapable of capturing multi-spin
correlations, it captures the single-spin dynamics sufficiently accurately,
and hence reproduces Kohlrausch relaxation of the disorder-averaged auto-
correlation function with excellent accuracy.

The plan of the paper is as follows. In Sec II we demonstrate the Kohlrausch
law in a power-law interacting Floquet setting, where the degree of
localization is controlled by tuning various parameters. We also illustrate
the salient phenomenological aspects of the Kohlrausch regime, especially the
occurrence of slowest perceptible dynamics (a “peak” in the relaxation
timescale) between the strongly and weakly localized regimes. In Sec. III we
demonstrate the phenomenology for next nearest-neighbour interacting systems –
both Floquet DTC-type and static ones. In Sec. IV.1 we provide a simple
physical picture explaining the occurrence of the peak. The picture is based
on identifying roughly three different relaxation regimes separated by their
characteristic time-scales. In Sec.IV.2, we resort to a non-hermitian
Hamiltonian description of a single spin in the system. For this we provide an
analytical expression for the dynamics of the single-spin auto-correlation
function, which, upon disorder averaging, exhibits the Kohlrausch relaxation
with remarkable accuracy. In Sec. V we investigate the occurrence of the
putative critical DTC phase in one spatial dimension.

##  II The Kohlrausch Relaxation in Long-Range Interacting driven MBD System

###  II.1 Model and Observables

Figure 2:  The Kohlrausch regime and its phenomenology: (a) Behaviour of the
stretching exponent $\gamma$ as a function of the exponent $\alpha$ of the
power-law interaction (fixing $W=0.4\pi$, $\epsilon=0.2$). In the parameter
regime marked with yellow, $\gamma$ shows systematic deviation from unity as
$L$ is increased indicating stability of Kohlrausch law with increasing $L$ in
that region. (b) Perceptible relaxation time-scale $\tau$ vs $\alpha$ for
various $L.$ A clear peak is observed in the intermediate values of $\alpha$.
(c) The residual auto-correlation, $a$ (see discussion after Eq. 7), as
$t\to\infty$ obtained from evolution up to a finite time $T_{obs}=10^{4}$).
Lower Panel: $\gamma$ vs $L$ for various variations of parameters:(d) $\alpha$
(fixing $W=0.4\pi$, $\epsilon=0.2$) , (e) $W$ (fixing
$\alpha=2.6,\epsilon=0.2$) and (f) $\epsilon$ (fixing $W=0.4\pi$,
$\alpha=2.6$). In each frame, the curves in shades of red shows increase in
$\gamma$ towards unity as $L$ is increased (instability of the Kohlrausch
regime with increasing $L$), and those in shades of blue shows departure of
$\gamma$ from unity with increasing $L,$ indicating stability of the
Kohlrausch regime against increase in $L$.

In this section we consider the following Floquet-Unitary drive Hamiltonian in
$1d$.

$U_{F}=\exp{\left[-i\sum_{ij}\frac{J_{ij}}{r_{ij}^{\alpha}}X_{i}X_{j}\right]}\exp{\left[-i(\pi+\epsilon)\sum_{i}Z_{i}\right]}$
(2)

where $J_{ij}$ is chosen uniformly from $[0,W]$, and $X_{i},Z_{i}$ are Pauli
operators on site $i.$ A cycle consists of the application of this unitary
once on the wave-function. The evolution within a cycle can be viewed as an
evolution over a period $T=1+\phi,$ where $\phi=(\pi+\epsilon)$ under the
following piece-wise local time-periodic Hamiltonian.

$\displaystyle H(t)=$ (3) $\displaystyle H_{z}=\sum_{i}Z_{i},~{}~{}~{}nT\leq
t<nT+\phi$ $\displaystyle
H_{x}=\sum_{ij}\frac{J_{ij}}{r_{ij}^{\alpha}}X_{i}X_{j},$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}nT+\phi\leq
t<(n+1)T.$

where $n$ are positive integers including $0.$

The long-range interactions above can be realizable and often even natural in
experimental setups like [20, 21, 24, 37, 38, 39, 40, 41]. Those setups are
likely candidates for implementing variants of the Hamiltonian given in Eq.
(3).

The numerically exact simulation of the Floquet dynamics of the above form can
be carried out by efficiently employing the Hadamard transform, reaching
system sizes larger than those accessible via exact diagonalization (see, e.g.
[42]). Here we consider the (site-resolved) auto-correlation function defined
for spin $i$ as

$\displaystyle C_{auto}^{i}(nT)$ $\displaystyle=$ $\displaystyle\langle
X_{i}(0)X_{i}(t)\rangle,~{}t=nT$ (4) $\displaystyle C_{auto}(nT)$
$\displaystyle=$ $\displaystyle\frac{1}{L}\sum_{i}C_{auto}^{i},$ (5)

where $X_{i}(t)$ is the operator $X_{i}$ after $n$ complete cycles in the
Heisenberg picture.

For $\epsilon=0,$ $C_{auto}$ is trivially frozen to unity. For small but
finite $\epsilon$ and large enough $W$ the above Hamiltonian has been argued
to exhibit a DTC-MBL order for a short (effective) interaction range, i.e.
large interaction exponent $\alpha$ (see, e.g. [9, 10]). In a DTC, the sign of
the auto-correlator flips every period, and this sub-harmonic oscillation
persists in spite of the fact that the system is interacting and non-
integrable. This happens for generic initial states, since each Floquet
eigenstate is argued to be perturbatively close (in $\epsilon$) to a cat state
formed by a simultaneous eigenstate of $\\{\sigma_{i}^{x}\\}_{i=1,L}$ and its
spin-flipped partner. It is known from numerical studies that regardless of
whether there is an underlying stable MBL-DTC phase, the initial relaxation in
an MBD system can slow down due to the disorders (see, e.g. [43, 44]). To
capture this slowness, we track the dynamics of the time staggered DTC order
after each cycle, given by

$\begin{split}\tilde{C^{i}}(nT)=(-1)^{(n+1)}C^{i}_{auto}(nT)~{}{\rm and}\\\
\tilde{C}(nT)=(-1)^{(n+1)}C_{auto}(nT).\end{split}$ (6)

The Kohlrausch relaxation is manifested in the disorder-averaged values of the
above quantities denoted by $\overline{\tilde{C^{i}}(nT)}$ and
$\overline{\tilde{C}(nT)}$ respectively. The overhead bar denotes the
averaging over disorder realizations as well as over various random $x-$bit-
string initial states.

The Kohlrausch form for MBD Systems: The presence of disorder in our system
necessitates a formally slight but crucial modification of the standard form
of the Kohlrausch law of stretched exponential form (Eq. 1): we need to add a
constant to it (e.g., to take care of a stable localization within a finite
system-size or a finite observation time). This gives us the following form.

$y(t)=a+(1-a)\exp{\left\\{-\left(\frac{t}{\tau}\right)^{\gamma}\right\\}}.$
(7)

The added constant $a$ accounts for the residual (frozen) correlation as can
be estimated from the data measured over a finite interval. This might or
might not be the true residual auto-correlation within our resolution as
$t\to\infty,$ depending on our measurement interval and/or the existence of a
stable MBL phase (see the end of Sub-Sec. II.2.1 for further discussions on
this).

###  II.2 Main Numerical Results

Here we summarize the main numerical results. In our system, the degree of
localization can be tuned by three parameters. It can be increased by
increasing the disorder strength $W,$ while can be decreased either by
increasing $\epsilon$ (which enhances delocalization over the $x-$bit-string
states) or increasing the effective range of interaction by decreasing the
exponent $\alpha.$ Fig. 1, shows various cuts in the parameter space
representing the salient features of the relaxation: the disorder-averaged
real-time dynamics has been shown (dotted lines). For all the cases,
throughout the entire parameter regime, the numerical data is fitted
(continuous lines) very accurately with the Kohlrausch law (Eq. 7): The
largest absolute least square error for fitting is $<10^{-5}.$ This fits much
better than the other standard relaxations forms, namely, exponential and
power-laws with and without additive constants.

In Fig. 1$(a)$, plots for various values of $\alpha$ are shown for a fixed
value of $\epsilon$ and $W$. As $\alpha$ is decreased, the dynamics switches
from a strongly localized DTC-MBL regime for our finite system (where
$\tilde{C}$ is frozen close to its initial value) to the delocalized regime
(where $\tilde{C}$ decays rapidly to zero). In $(b)$ and $(c)$, the parameter
regimes are traversed (from stronger to weaker localization) by increasing
$\epsilon$ and decreasing $W$ respectively, keeping other parameters fixed.

####  II.2.1 Phenomenology in the light of Relaxation Parameters:

The relaxation parameters ($\gamma,\tau$ and $a$) are extracted from the
fitting of the Kohlrausch form (Eq. 7) to $\tilde{C}(t)$ measured up to a
finite observation time $T_{obs}.$ The aim of the fitting is to provide a
simple yet sufficiently accurate (leading order) description of a complicated
relaxation dynamics. The exact dynamics, of course have, several finer details
and also various dynamical regimes. Hence there is apriori no reason to expect
that the entire relaxation should be well approximatable by a single simple
function. The values of the best fitting parameters as well as the suitable
fitting functions can thus depend upon choice of the temporal regime under
analysis (i.e., on $T_{obs}$). In our case, though the parameters do exhibit
such a dependence (see Fig. 8), the functional form, quite interestingly, does
not – Kohlrausch form always provides the most optimal (over the set of usual
physical relaxation forms like exponential, power-law etc with or without an
additive constant etc) and also a very accurate description in all regimes.
Though we have fitted the whole regime with a single set of parameters in most
cases here, the separation between the early time and late time regimes are
actually quite pronounced (the short, early time dynamics and the longer late
time dynamics does not fit equally well with a single set of parameters). This
has been addressed in further details in Sec. V. In the following we discuss
the phenomenology and the significance of the relaxation parameters.

The Stretching Exponent $\gamma$ and the Stability of Kohlrausch Regime:
Though the Kohlrausch’s stretched exponential relaxation does not probably
look qualitatively very different from the thermal exponential relaxation, a
unified understanding of the stretching has still remained elusive and
enigmatic. The efforts in this direction (models with finite density of long-
lived metastable states, hierarchical constraints etc) indicate that a unified
explanation would require elements fundamentally different from those
responsible for exponential relaxation [26]. Hence the fate of the stretching
exponent $\gamma$ in the thermodynamic limit, and identification of those
regime where it is non-trivial in that limit, are issues of fundamental
importance.

In order to investigate the stability of the stretched exponential relaxation
against the increase in system-size $L$, we extracted the stretching exponent
$\gamma$ for various values of $L$ from real-time dynamics data over $10^{4}$
cycles averaged over disordered configurations. The results are shown in Fig.
2(a),(d),(e),(f). The Fig. 2(a) shows that there is a regime of $\alpha$
(shaded with yellow) where the stretching exponent $\gamma$ actually decreases
systematically away from unity with increasing system-size $L.$ This is what
we identify as the “slow-regime” or “Kohlrausch relaxation regime”. Away from
the slow regime, $\gamma$ increases steadily with $L,$ and as is intuitively
believed (see, e.g. [43]), $\gamma$ vs $L$ might have unity as its asymptote,
implying a thermal relaxation behaviour as $L\to\infty$. Between these two
regimes, the trend is non-monotonic and ambiguous in our finite system
simulation, and it is difficult to determine whether the two phases are
separated by a transition or a crossover.

The three frames in the lower panel of Fig. 2 shows the variation of $\gamma$
with $L$ both within (in shades of blue) and outside (in shades of red) the
Kolhrausch regime. In Fig. 2$(d)$ , the results are shown for various values
of $\alpha$ for a fixed value of $W$ and $\epsilon$ (corresponding to those in
$(a)$). In Fig. 2$(e)$ and $(f)$ the parameter space are traversed by varying
$W$ and $\epsilon$ respectively.

The possibility of the existence of an intermediate phase/regime between a
thermal and the non-ergodic phase in a disordered many body system has been
indicated long ago (see, e.g. [45] for a review) and some of the dynamical
consequences of such an anticipated regime have also been considered more
recently [46, 47]. The Kohlrausch regime puts all those observations on a
unified ground, yet it does not necessarily require a stable MBL phase as an
asymptote. Strong enough disorder which gives the local integrals of motions
[6, 5] a life-time $>T_{obs}$ is sufficient.

Figure 3: Kohlrausch relaxation in up to next-nearest neighbour interacting
Floquet system (Eq. 8). (a) Auto-correlator for various values of $\epsilon$,
for a fixed disorder strength, for $L=24.$ Dotted lines are actual data from
the numerics, and solid lines are fitting of the data with a stretched
exponential function defined at Eq. (7). (b) Time scale $\tau$ as a function
of $\epsilon.$ (c) Stretching exponent $\gamma$ as a function of $\epsilon.$
The inset shows the dependence of $\gamma$ on the system size for at various
values of $\epsilon.$ The stable and unstable flows of Kohlrausch exponents
with $L$ are shown respectively in blue ($\epsilon=0.25$) and red
($\epsilon=0.40$). Similar behaviour can be observed as a function of the
disorder strength for fixed values of $\epsilon$ (data not shown).  Figure 4:
Stretched exponential relaxation for quench with a time-independent MBD
Hamiltonian (Eq. 9). The inset shows the variation of the stretching exponent
$\gamma$ with the disordered strength $h.$

Peak in Relaxation Time-scale $\tau$: Next we focus on another striking aspect
of our phenomenology. The relaxation time-scale $\tau$ does not monotonically
increase as one moves from a relatively weakly localized regime to a more
strongly localized one. Instead, $\tau$ exhibits a lofty peak somewhere
between the two said regimes. This is shown in Fig. 2(b). This might appear
counter-intuitive at first, since it is sensible to associate slower dynamics
(larger $\tau$) with a more strongly localized regime. However, for a finite
$T_{obs}$ and a finite fitting resolution, there could be a substantial
fraction of spin which would appear completely frozen, and our finite time
observation will not sense their participation in the relaxation and club it
in the constant $a$ (discussed in further details in the following). An
explanation and interpretation of the peak and its behaviour with the change
in parameters have been given in terms of a simple physical picture based on a
few very reasonable physical assumptions in Sec. IV.1.

Interpretation of the Residual auto-correlation $a$: The parameter $a$
provides a measure of the degree of localization as extractible from the data
for dynamics over a finite $T_{obs}.$ Our finite-time characterization of weak
and strong localization is based on the value of $a$ ($a=0$ implies
delocalization where $a=1$ means perfectly frozen). However, it may be noted
that $a$ is neither the residual auto-correlation after time $T_{obs}$ nor
(necessarily) the remnant correlation observed after infinite time (i.e., for
$T_{obs}\to\infty$). It is the asymptotic auto-correlation estimated from the
data obtained over a finite $T_{obs}.$ If there were only one dynamical regime
that could be captured accurately by a single fitting function with a fixed
set of parameter values, then of course everything would be independent of
$T_{obs},$ and $a$ could be interpreted as the residual correlation after
infinite time. But as behooves a complicated many-body system, an MBD system
exhibits various distinct stages of dynamics in the sense that a single
function with a fixed set of parameter values cannot capture the dynamics over
all time scales. Hence, $a$ just gives an estimate of the part of the
correlation that looks completely frozen over an observation time-scale
$T_{obs}$ within the accuracy with which the Kohlrausch description is
expected to hold. Thus, a non-vanishing $a$ extracted by finite-size analysis
from a finite $T_{obs}$ data does not necessarily imply existence of a stable
MBL.

##  III Short Range Variants: Floquet and Static MBDs

In this section we consider another Floquet-Unitary drive for a strictly short
range system with interactions ranging up to the next nearest neighbour on the
same lattice geometry as before. The evolution operator over a period is given
as follows.

$\displaystyle U_{F}$ $\displaystyle=$
$\displaystyle\exp{\left[-iH_{x}\right]}\exp{\left[-i(\pi+\epsilon)H_{z}\right]}$
$\displaystyle where,~{}H_{z}$ $\displaystyle=$ $\displaystyle\sum_{i}Z_{i},$
$\displaystyle H_{x}$ $\displaystyle=$
$\displaystyle-\sum_{i}J_{i}X_{i}X_{i+1}+J_{2}\sum_{i}X_{i}X_{i+2}$ (8)

$J_{i}$’s are chosen uniformly from $[0.1\pi,0.4\pi].$ A static version of the
model in a transverse field is believed to show a localization-delocalization
transition [48] as the disorder strength is tuned. We choose a parameter
regime where the localization is strong enough, though, like in the earlier
case, is not based on existence of a stable MBL phase.

The evolution shows qualitatively same phenomenology as obtained for the long-
range system, as shown in Fig. 3. Here the degree of localization is
controlled by tuning the imperfection $\epsilon$ ($\epsilon=0$ is trivially
localized). The left frame shows real-time dynamics fitted with the stretched
exponential form in Eq. 7. The middle panel show the peak in the relaxation
time-scale $\tau$ as a function of $\epsilon.$ The right frame shows the value
of the stretching exponent $\gamma$ as a function of $\epsilon$ for various
values of $L.$ A stable regime for Kohlrausch relaxation is clearly observed
(where $\gamma$ decreases with increasing $L.$) The inset shows the variation
of $\gamma$ with $L$ in (blue) and away (red) from the stable regime.

Next, we consider the anisotropic Heisenberg model [49], in a disordered
field, and study the relaxation dynamics of a state under this Hamiltonian:

$H=J\sum_{i}(X_{i}X_{i+1}+Y_{i}Y_{i+1})+J_{z}\sum_{i}Z_{i}Z_{i+1}+\sum_{i}h_{i}Z_{i}$
(9)

$h_{i}$ are chosen uniformly from $[h,-h]$, $J_{z}=0.5J,J=1.$

##  IV A Phenomenological Picture of the Kohlrausch Scenario and A Non-
hermitian Analytical Approach

###  IV.1 A Simple Physical Picture

In this section, we provide a cartoon of a qualitative picture underpinning
the phenomenology entailing Kohlrausch relaxation described in Sec. II. A
remarkable heterogeneity in the microscopic dynamics underpins the simple
Kohlrausch law, as can be gleaned from a representative plot of the real-time
dynamics of single spins chosen randomly from various samples (Fig. 5) for a
set of parameter values for which the Kohlrausch law remains stable with
increasing $L$. Here we provide clarifications of central aspects of the
phenomenon: we explain salient features of the Kohlrausch phenomenology (e.g.,
the occurrence of the peak in $\tau,$ $L-$dependence of $\gamma,$ the behavior
of the peak in the $T{obs}\to\infty$ as well as $L\to\infty$ limit) employing
a simple physical picture.

Figure 5:  Real-time dynamics of single spins sampled completely randomly from
various realizations, for $\alpha=2.5,W=0.4\pi,\epsilon=0.2.$ The plot shows a
high degree of variability in the local dynamics with different types of
behavior appearing with appreciable probabilities.

Figure 6: The physical picture: Top: A sequence of cartoons of three
relaxation scenarios. (a) Strongly localized, dominated by frozen auto-
correlation (blue arrows), (b) Kohlrausch regime, dominated by slowly decaying
auto-correlation (green arrows), and (c) mostly delocalized, dominated by
rapidly relaxing auto-correlation (red arrows).
Bottom: Histogram of single-spin auto-correlations $C_{auto}^{i}$ (Eq. 5 )
sampled over an ensemble of randomly chosen 5000 spins, corresponding to the
respective cartoon above. Three different colors show results after three
successive decades of evolution. The mean single-spin auto-correlation of the
ensemble at three different times is shown by the vertical lines.

The Basic Ingredients: The Slow, The Fast, and The Frozen: Our physical
picture is based on a cartoon (Fig. 6 (a)-(c)). The observation window
$T_{obs}$ defines the reference timescale. Though there is a broad
distribution of relaxation times, the cartoon is built on three variants in
particular: spins that relax fast (red arrows in the cartoon), slow spins –
whose relaxation is tiny but discernible during the process of the observation
(green arrows), and spins whose relaxation (if any) within $T_{obs}$ is beyond
our resolution for detecting any dynamics (blue arrows). The dynamics can be
broadly classified into three regimes (as described below), and its features
can be qualitatively captured in terms of the dynamics of these cartoon spins.
The relative populations of these three variants of spins marks three regimes
of dynamics in our cartoon, as described below.

The Strongly Localized Regime: This is a regime where a large fraction of the
auto-correlation is frozen (Fig. 6(a)). This corresponds to the topmost blue
(almost frozen) curves in Fig. 1. In this regime, within $T_{obs}=10^{4},$ the
relaxation happens only due to residual small fraction $(1-a)$ of fast (red
spins) and slowly decaying (green spins) parts of the auto-correlation This
corresponds, for example, to the large $\alpha$ regime on the right of the
peak of $\tau$ in Fig. 2(b). Here the $L-$dependence of the stretching
exponent does not exhibit any systematic trend within system-sizes accessible
to us, and stability of the Kohlrausch law cannot be inferred.

The Kohlrausch (Slow) Regime and the Occurrence of the Peak in $\tau$: Fig.
6(b) shows the situation as we tune the system more towards the delocalized
regime either by increasing $\epsilon$ or decreasing $\alpha$ and/or $W.$ A
considerable number of frozen (blue) spins are liberated marginally from the
strong constraints (that had kept them frozen) to have restricted dynamics and
thereby got converted into slow (green) spins. In comparison, the growth of
the fraction of the fast (red) spins (presumably, mostly from the already
existing slow spins) is not appreciable yet. The dynamics is hence dominated
by the slow spins. This corresponds to the curves with intermediate relaxation
rates (the pale orange plots in the midway of Fig. 1), and the region around
the peak in $\tau$ in Fig. 2(b). Here the growth rate of the red spins takes
over the growth rate of the green spins as one moves towards a more
delocalized regime, and the peak occurs roughly at the point where the number
of red spins overtakes the number of green (and rarely occurring blue) spins.
This is the regime (marked by vertical yellow stripe in Fig. 2) where the
Kohlrausch law is stable to increase in $L$ -(i.e. $\gamma$ drifts off from
unity as $L$ is increased). The stability is concurrent with the dominance of
the slow spins.

The Delocalized Regime: Here most of the spins are free to relax fast,
resulting in dynamics which is very close to thermal (Fig. 6(c)). This
corresponds to the most rapidly decaying red curves in Fig. 1 and the low
$\alpha$ regime on the left of the $\tau$ peak in Fig. 2 (b). In this regime
the stretching exponent $\gamma$ systematically increases with $L$ approaching
unity, indicating an exponential relaxation for $L\to\infty.$

Numerical Signatures of the Regimes: The frames in the lower panel of Fig. 6
show the number distribution of spins over auto-correlation (i.e., single-spin
auto-correlation versus the fraction of spins having that particular value of
the auto-correlation), over an ensemble of $5000$ spin taken from 5000 various
disorder realizations, for each of which a random $x-$bit-string state (a
simultaneous eigenstates of all $\sigma_{i}s$) is taken as the initial state.
Each frame in the lower panel corresponds to the physical the situation
depicted in the cartoon frame right above it in the upper panel. The plots
illustrate the corresponding cartoons. Fig. 6 (d) represents the most strongly
localized regime, where there is no significant change in the mean of the
distribution (vertical lines representing the ensemble average of the auto-
correlation with colors corresponding to various times) with time, especially
in the later decades. Fig. 6(e) represents the Kohlrausch regime, where the
change in the mean takes place through the later decades, but still a major
fraction of the auto-correlation remains intact. Fig. 6(f) represents the fast
relaxing regime, where rapid decay is observed throughout the entire duration,
and also shifting the mean from 1 to 0. The consistency of the above behavior
with the numerical results for the Kohlrausch phenomenology can be gleaned, at
once by recalling, for example, the correspondences of (d), (e), and (f) with
the three regimes – the one to the right of the yellow band (strongly
localized), the one within the yellow band (the Kohlrausch), and the one to
the left of the yellow band (delocalized) of frames (a), (b), and (c) of Fig.
2.

###  IV.2 Single-Spin Analytics: A Non-Hermitian Random Hamiltonian Approach

Here we resort to a simple formulation suitable for capturing the key elements
behind the Kohlrausch law and its associated phenomenology in an MBD system.
Since the single-spin auto-correlation is the central object of interest here,
the key idea is to focus on the dynamics of a single spin (an arbitrary target
spin), treating the rest of the system as a bath. The effects of the disorder,
interaction, and drive are thus encapsulated in the action of the bath on the
target spin. Here we resort to a minimal description for the above setting,
consisting of a spin evolving under a single-spin non-hermitian Hamiltonian.
The non-hermiticity is essential for modeling aspects like dissipation, energy
gain, and dephasing of the dynamics of the target spin. The stroboscopic
dynamics under the time-periodic Hamiltonian can be viewed as stroboscopic
observation of the dynamics under the time-independent Floquet Hamiltonian
$H_{eff}$ (where the evolution operator over a period $U(T;0)=e^{-iH_{eff}}$),
the non-hermitian Hamiltonian, in principle, is to be obtained by tracing over
all the degrees of freedom except the target spin, evolving under $H_{eff}$
(time-independent). Thus for stroboscopic dynamics, it is sufficient to
consider a time-independent non-hermitian Hamiltonian whose parameters encode
the information of the rest of the system as well as the drive. The problem is
soluble, and we treat the parameters in the solution as free ones to be
extracted from the numerical results. To this end, we devise a simple single-
spin non-hermitian Hamiltonian formulation as follows.

####  IV.2.1 The Formulation

We describe the dynamics of an arbitrary target spin in our system (observed
stroboscopically at $t=nT$) by that governed by a time-independent non-
hermitian Hamiltonian as follows.

$\displaystyle\tilde{H}$ $\displaystyle=$ $\displaystyle H_{h}+H_{nh},~{}{\rm
where}$ $\displaystyle H_{h}$ $\displaystyle=$
$\displaystyle-\Omega\sigma^{x}~{}{\rm and}$ $\displaystyle H_{nh}$
$\displaystyle=$ $\displaystyle-i\left[a_{0}\mathcal{I}_{2\times
2}+a_{1}\sigma^{x}+a_{2}\sigma^{y}+a_{3}\sigma^{z}\right]$ (10)

Here $\sigma^{x/y/z}$ are the Pauli matrices representing the components of
the target spin, $a_{i}$ $(i=1,4)$ are random parameters related to the
parameters of the corresponding random realization of the time-dependent many-
body system in which the target spin is embedded. Note that when focused on a
single spin, $\Omega$ is just an overall scale factor. However, for an
ensemble of spins embedded in a random many-body system, $\Omega$ will be an
overall relative scale of relaxations among the spins, e.g., a larger $\Omega$
will correspond to a faster spin, while a smaller $\Omega$ will correspond to
a slower spin, given all other parameters are same between them. Our
approximation consists of the assumption that $\tilde{H}$ encapsulates the
interaction of the target spins with the environment as well as the drive.
This is a rather drastic approximation, and is suitable only for calculating
single-spin observables. Since we are calculating auto-correlations in
$X$-components of the lattice spins, we represent the $x-$component $X_{i}$ of
a target spin (say, at a site $i;$ see Eq. 3) to be diagonal in the $\sigma$
representation, i.e., $(-1)^{n+1}X_{i}(nT)\leftrightarrow\sigma^{z}(nT).$ Note
that, since we are representing the dynamics of a spin by a time-independent
$\tilde{H},$ it is necessary that we incorporate the factor $(-1)^{n+1}$ in
the correlations to set the equivalence correctly.

We switch to a re-defined set of parameters (after [50]) to link them more
clearly to various physically intuitive (as has been discussed down the line)
quantities as follows.

$a_{0}=\gamma,~{}a_{1}=\gamma\beta,~{}a_{2}=\nu,~{}a_{3}=V.$ (11)

In general, the state of the target spin is expected to be represented by a
density $2\times 2$ matrix $\rho(t)$ which satisfies the following Schrödinger
equation.

$i\partial_{t}\rho(t)=-i[\tilde{H},\rho(t)].$ (12)

Since $\tilde{H}$ is non-hermitian, $Tr[\rho(t)]$ is not necessarily conserved
(here $Tr[.]$ denotes trace). The expectation value of any observable ${\cal
O}$ is hence needed to be normalized as follows.

$\langle{\cal O}\rangle(t)=Tr[\rho(t){\cal O}]/Tr[\rho(t)].$ (13)

The Schrödinger equation Eq. (12) is exactly soluble [50] for any real set of
values of the parameters in $\tilde{H}$. The general solution reads:

$\begin{split}\rho_{mn}&=\frac{e^{-\Gamma
t}}{2({\beta}^{2}+{\gamma}^{2})}\left[A_{mn}\cos{(\omega
t)}+B_{mn}\sin{(\omega t)}\right.\\\ &\left.+C_{mn}\cosh{(\Gamma
t)}+D_{mn}\sinh{(\Gamma t)}\right],\end{split}$ (14)

where

$\Gamma=2\Omega\gamma,~{}\omega=2\Omega\beta.$ (15)

Here $\rho_{mn},$ ($m,n\in\\{1,2\\}$) are the four elements of the density
matrix $\rho(t),$ and $A_{mn},B_{mn},C_{mn}$ and $D_{mn}$ are constants
depending on the initial condition and the Hamiltonian parameters for a given
disorder realization. This solution is general enough to cover various
scenarios including (but not limited to) exponential and polynomial
relaxations with or without oscillations, energy-conserving dynamics, pure
dephasing, etc as special cases depending on the values of the Hamiltonian
parameters even for a given initial state (see, e.g., [50]). This versatility
of the exact solution makes it suitable for capturing the diversity of single-
spin relaxation (as illustrated in Fig. 5). Here we focus on a fully polarized
eigenstate of $\sigma^{z}$ ($X_{i}$) with eigenvalue $+1,$ i.e., all spins up
in the $x-$direction in the chain in Eq. 3. Under this initial condition, the
time-staggered single-spin auto-correlation

$\tilde{C}^{i}_{auto}=(-1)^{n+1}\langle
X_{i}(t)\rangle=\langle\sigma_{z}(t)\rangle,$ (16)

where $t=nT.$

The expression for $\tilde{C}^{i}_{auto}$ (expressed in terms of $\sigma^{z}$)
for the above-mentioned initial condition reads (using Equation 14)

$\langle\sigma^{z}(t)\rangle_{Analytic}=\frac{[(1+\gamma^{2}-\nu-V^{2})\cos{(\omega
t)}+(\beta V)\sin{(\omega t)]\operatorname{sech}{(\Gamma
t)}+(\nu+\beta^{2}+V^{2}-1)+\gamma V\tanh{(\Gamma
t)}}}{[(\nu+\beta^{2}-1)\cos{(\omega t)}+(\beta V)\sin{(\omega
t)}]\operatorname{sech}{(\Gamma t)}+(1+\gamma^{2}-\nu)+\gamma V\tanh{(\Gamma
t)}},$ (17)

The expression contains the parameters $\Omega,\gamma,\beta,\nu$ and $V$, that
are essential for characterizing a $2\times 2$ non-hermitian matrix (Eqs. 10
and 11), and are expected to encode the necessary effects of the drive and the
other spins on the relaxation of the target spin.

Figure 7: Fitting of the exact numerical results for the time staggered
single-spin auto-correlations with the analytical formula from the non-
hermitian formulation. Left: main frame shows some representative fitting of
the numerical data (colored) by the analytical function (solid black lines).
The inset shows the fitting of the numerical disorder averaged relaxations
with the disorder averaged analytical expression corresponding to the case
$\alpha=2.8,W=0.4\pi,\epsilon=0.2,L=26$ (Kohlrausch regime). The fitting is
quite accurate and satisfies Kohlrausch’s law. Middle: frame shows the
probability distributions $P_{{}_{\gamma}}$ of $\gamma.$ Right: frame shows
the probability distributions $P_{{}_{\Omega}}$ of $\Omega.$

It is easy to glean from Eq. (17) that $\gamma$ sets the time-scale of energy
loss/gain by the target spin, $\beta$ provides the characteristic frequency of
the damped oscillation executed by the target spin.

Since energy is pumped into the system, there are generically exponentially
growing terms in the subsystem dynamics, as well as exponentially decaying
terms due to both drive and the loss of energy to the rest of the system. A
steady state appears when those two rapidly changing terms counterbalance each
other. In the above form (Eq. 17), the expression for
$\langle\sigma^{z}(t)\rangle_{Analytic}$ this is in-built, and averts the
usual problem with divergences associated with individual exponentially
increasing terms as $t\to\infty$ and the associated numerical problems at the
large but finite time: the magnitude of both $\operatorname{sech}{(x)}$ and
$\tanh{(x)}$ are bounded between $0$ and $1$.

####  IV.2.2 Auto-Correlation Dynamics and the Kohlrausch law from the
Analytics

We fit the single-spin auto-correlation data obtained from exact numerics (as
shown in Fig. 5) with the analytical expression for
$\langle\sigma^{z}(t)\rangle_{Analytic}$ in Eq. (17), for a large number of
disorder realizations, and extract the values of the parameters
$\Omega,\gamma,\beta,\nu,$ and $V$ in each case. Plugging such a set of
parameter values into the expression of
$\langle\sigma^{z}(t)\rangle_{Analytic},$ we get our analytical expression for
$\tilde{C}_{i}(t)$ with numerical values for its parameters. We average this
expression over various disorder realizations.

The large varieties of single-spin relaxations obtained numerically (Fig. 5)
is captured by our analytic expression obtained from a simple formulation
(Fig. 7, Left frame, main) for a few random samples. The inset shows the fit
of the disorder-averaged numerical result by the disordered averaged
analytical expression. They are both consistent with the Kohlrausch law to
very good accuracy over multiple decades.

The probability distribution $P_{{}_{\gamma}}$ of the relaxation rate $\gamma$
(middle frame) also shows two clear peaks. The population at the peak at
$\gamma=0$ exhibits no relaxation at all (blue spins), while those within the
width of it show slow relaxation (green spin), and the second peak represents
spins with considerable relaxation. The bimodal structure appearing in the
distribution of these two parameters with an appreciable width around the peak
at zero is consistent with the classification of our spins into three
categories, namely, the frozen (blue), the slow (green) and the fast (red)
based on their relaxation rates. The distributions for $\nu,$ $\beta$ and $V$
are shown in Appendix A. Similarly, the distribution of $\Omega$ (Fig. 7,
right frame) exhibits two clear peaks - one at $\Omega\approx 0$, another
around $\Omega\approx 0.06$. The peak at $\Omega=0$ represents the spins
frozen to our resolution (the blue spins). The populations within the width
(suitably defined) of this peak represent the slow (green) spins. The second
peak and those within its width represent spins with faster relaxations (the
red spins in our cartoon).

##  V Early vs Late Stages of Relaxation: A Numerical Search for the
“critical” DTC Phase

Figure 8:  Dependence of the relaxation parameters with the observation time
$T_{obs}$ indicates the existence of two distinct dynamical stages. At the
early stage (roughly corresponding to $T_{obs}<500$), the relaxation
parameters are quite sensitive to the observation time $T_{obs}$ (especially
for the shorter-ranged systems). At the later times (the late stage), the
relaxation parameters tend to saturate and their dependence on $T_{obs}$ is
much weaker or negligible.

The existence of various stages of relaxation is not uncommon in interacting
many-body systems. For example, a system described by a fractional Fokker-
Planck equation can exhibit relaxation that crosses over from stretched
exponential to an inverse power law (power law with negative exponent) [51].
Here we found, quite interestingly, that the stretched exponential form of Eq.
7 continues to fit both the early and late stages the best (compared to other
more common relaxation forms like exponential and power-laws) and also
remarkably accurately. The fitting parameters are distinctly different for the
two stages distinguising one from the other as shown in Fig. 8. One should
however keep in mind, that the time-scale of the emergence of the late regime
is still negligibly small compared to the characteristic time-scale of
relaxation (if any) of the frozen (blue) spins (the latter being infinity to
our resolution).

The phenomenology of relaxation in the early and late stages is distinct as
discussed below. Fig. 8 (top frame) shows, initially the stretching exponent
$\gamma$ is closer to unity (i.e., the relaxation is more thermal-like) for
larger $\alpha$ (shorter interaction range). But between $T_{obs}=100$ and
$200,$ the trend reverses, and $\gamma$ drifts further from unity for larger
$\alpha.$ The distinction between the two stages of relaxation is most
pronounced for the short-ranged systems and gradually levels off as the range
of interactions is increased. This can be roughly understood as follows. In
shorter-ranged MBD systems, for a sufficiently strong disorder, there is
localization (regardless of whether it persists as $t\to\infty$) at the
initial stage of the dynamics (assuming we start with a localized state). The
early stage consists of dynamics within the associated localization length.
This is followed by a much slower relaxation process involving equilibration
between different localized regions. With an increase in the range of the
interaction, the size of the localized patches increase. The relaxation
dynamics hence involve cooperative dynamics of a larger number of spins. This
makes the relaxation slower than the shorter-range system at early times (see,
e.g., [52]). The early stage is thus characterized by faster relaxation
(larger $\gamma$) for the shorter-ranged systems. The late stage is marked, on
the other hand, by the equilibration dynamics between the localized patches
(beyond the localization length), which is much slower in shorter-range
systems due to stronger localization compared to the longer-ranged ones. This
is hence characterized by a smaller value of $\gamma$ (more stretched
relaxation) for shorter-ranged systems.

The persistence of Kohlrausch law in almost all stages of relaxation is also
interesting in view of a prediction of a dynamical phase called “critical” DTC
in long-range interacting Floquet systems with unitaries of the form given in
Eq. 2 in various spatial dimensions[36, 10]. A perturbative expansion of an
effective Hamiltonian in a suitable frame under specific conditions predicts a
scale-free relaxation of a DTC order parameter (the auto-correlation
$\tilde{C}$ we measured) at the early time [36, 10] in a power-law interacting
setting whose one-dimensional version we studied (Eq. 2). The dynamical phase
is called critical DTC, and is predicted to occur when the exponent $\alpha$
of the power-law interaction equals the dimension of the lattice ($\alpha=1$
for our setting).

Focusing on the early-time data (the first two decades) we could not find the
proposed critical DTC phase in spite of a careful scanning of the parameter
space. Instead, we found Kohlrausch law persists across the board for
$\alpha=1,$ as shown in Fig. 9.

Figure 9:  Short-time behaviors and fitting: Showing no sign of power-law
decay. A log-log plot does not show a straight line and fits better (least
squares being always $<10^{-6}$) with a stretched exponential function (We
tried variants of exponential fits, which always gives least squares
$~{}10^{-2}$, data not shown). Top Panel: Main plot is for
$L=24,\epsilon=0.2,\alpha=1.0$ and the inset shows the system size scaling of
the exponent $\gamma$ for various values of $W.$ Bottom Panel: Main plot is
for $L=24,W=0.4\pi,\alpha=1.0.$ The inset shows the absence of power-law
relaxation for a range of values of $\alpha$ around $\alpha=1.$

##  VI Discussions and Outlook

Spin Glass Mechanism Plays No Essential Role in Kohlrausch Stretching of
Relaxation: Glass formation is known to one of the key mechanisms that give
rise to stretched relaxation (see, e.g. [26]). In our spin set-up, formation
of a quantum spin glass (see, e.g., [53] for a review) could hence have been a
candidate behind the Kohlrausch relaxation. However, two essential ingredients
of a spin glass are (a) disorder and (b) frustration [54, 55]. Our results
indicate that Kohlrausch stretching of relaxation is in fact independent of
the latter.

In our setting, frustration is introduced if a fraction of the interactions
$J_{ij}$ are made negative in the long-range case – there is no frustration
when all the bonds are ferromagnetic, while the presence of a sufficient
fraction of anti-ferromagnetic bonds can introduce the frustration necessary
to stabilize a spin-glass in a long-range interacting disorder system [27]. We
find that Kohlrausch stretching is present even in systems with completely
ferromagnetic couplings ($J_{ij}\geq 0$). As a strong example, we analytically
show (see Appendix B) that taking any arbitrary bit-string state in $x$
direction as the initial state, the dynamics with all $J_{ij}>0,$ in the
interaction part $H_{x}$ (Eq. 3), i.e., completely ferromagnetic interactions
are exactly identical for that in which each $J_{ij}$ has flipped sign (highly
frustrated one).

On Plausible Saturation of the Peak with increasing $L$: The fraction of
frozen spins observed within $T_{obs}$ depends on $L.$ All we assume is that
it asymptotes to some non-zero value even as $L\to\infty.$ This only requires
the characteristic time for melting of the Local Integrals Of
Motion/$\ell-$bits (see, e.g., [6]) to be larger than $T_{obs}.$ Though we
have not seen the saturation within the system size we can access numerically,
on physical grounds, the assumption seems quite justified for strong disorder
- the thermalization is shown to be extremely slow even where MBL is argued to
be absent [17], and hence the frozen part of the correlation is expected to
persist even for substantially long $T_{obs}$.

As $L$ is increased, the peak in $\tau$ grows in height (as can be seen from
Fig. 2(b)). However, for a finite $T_{obs},$ the peak height will eventually
saturate beyond some very large value $L\sim L_{max}$ (not accessible to
numerics). This may be understood as follows. We assume that the number
density of spins with effective relaxation time-scale $\tau$ is $n(\tau),$
than the average separation between such spins will be $1/n(\tau).$ Then it is
easy to see that the average height of the peak can grow as $L$ is increased
till some $L_{max}\sim 1/n(\tau_{max})$ where $\tau_{max}$ is the smallest
value of $\tau$ for which $\int_{\tau_{max}}^{\infty}\tau n(\tau)d\tau\approx
0$ (provided such a $\tau_{max}$ exists). Till $L$ reaches $L_{max},$
increasing $L$ might result in an inclusion of much slower spins (much larger
$\tau$) with an appreciable probability. This might lift the average value of
$\tau$ and hence the peak height (as it does in our case). In the following,
we argue the existence of a $\tau_{max}$ in our setting. Since we assumed a
finite $T_{obs}$ and finite resolutions of the auto-correlation measurement,
we cannot resolve a very large $\tau$ from infinity, from the auto-correlation
measurement (for example, for $T_{obs}=10^{4},$ a spin with $\tau=10^{15}$
will be considered as frozen with a resolution that cannot register any change
in correlation taking place over a period of $10^{-11}\times\tau$). The cut-
off value of $\tau$ beyond which the change in auto-correlation cannot be
registered thus serves as a natural $\tau_{max}.$ We consider spins with
$\tau>\tau_{max}$ as completely frozen and hence non-participating in the
measurable correlation decay (their contribution constitutes the atemporal
part $a$ of the auto-correlation). Hence, taking into account only the
remaining spins, we have $n(\tau)=0$ for $\tau>\tau_{max}.$ Since this is
indeed the case ($a\neq 0$) in the Kohlrausch regime where the peak occurs, a
$\tau_{max}$ exists in this regime for which $\int_{\tau_{max}}^{\infty}\tau
n(\tau)d\tau=0,$ and the peak will saturate beyond $L_{max}$ associated to it.

The peak also moves deeper into the more localized regime as $L$ is increased.
That could simply happen because the putative MBL phase seems to move further
away with increasing system size [14, 17, 18] in most cases. However, since we
are focusing on a finite $T_{obs},$ by our assumption that there will always
be a residual auto-correlation even in the thermodynamic limit for
$t<T_{obs},$ the peak should converge to some finite value as $L\to\infty.$
Unfortunately, this saturation is also beyond the system sizes accessible to
our numerics.

The Peak as $T_{obs}\to\infty:$ If there is a stable MBL for a finite value of
the disorder strength, then the above picture in principle can hold even as
$T_{obs}\to\infty.$ This is because there will always be some blue spins, and
it is plausible that as one moves towards a less localized regime from a more
localized one, some of them will be liberated to give rise to the green spins.
Of course, there will be a competition between the infinities: the infinity
associated with the time-scale of stability of the MBL and the infinity
associated with $T_{obs}.$ On the other hand, if there is no stable MBL
anywhere, then as $T_{obs}$ is increased further, more blue spins will
actually qualify as the green spins, and the peak height will increase, and
eventually saturates to the longest delocalization time-scale as
$T_{obs}\to\infty$.

Conclusions:

We showed that the century-old law for relaxation of complex classical systems
(structural glass, super-cooled liquid, etc) – the famous Kohlrausch law –
also appears ubiquitously in the domain of isolated MBD systems – both driven
and time-independent – evolving under pure Schrödinger dynamics.

This helps us quantify the slowness of the well-known prethermal relaxation in
MBD systems.

Finite-size analysis shows systematic stabilization of this law with
increasing system size in certain parameter regimes, allowing us to define a
dynamical phase characterizing the effect of disorder.

The relaxation time $\tau$ of auto-correlation over a finite observation time
$T_{obs}$ shows a peak between the weakly and most strongly localized
(characterized by the frozen part of the auto-correlation as appears from a
finite time observation) regimes.

The phenomenon does not necessarily depend on the existence of a stable MBL
phase in the thermodynamic limit at a finite parameter strength (e.g.,
disorder strength $W$, power $\alpha$ of the power-law interaction, etc) that
controls the degree of localization, though, the existence is assured if there
is a stable MBL phase.

We provide a simple analytical formulation of the problem, targeting a single
spin in the MBD system. The approach is based on random non-hermitian
Hamiltonians. Our analytical formula provides an accurate description of the
single spin dynamics even at the level of each individual target spin and
consequently captures the Kohlrausch law very accurately upon disorder
averaging. The characteristics of the MBD phase are captured through the
behavior of the system as a bath to a random target spin, and are encapsulated
in the distribution of the parameters in the single-spin non-hermitian
Hamiltonian.

Outlook:

The existence of some sort of slow dynamical phase beyond the MBL regime is
probably already in our indirect knowledge: It was known that one sees
logarithmic growth of entanglement in the MBL phase [49], but now it is known
that MBL is unstable in a substantial part of that regime [14, 15, 18]. Our
results are hence expected to trigger a new direction of study toward
understanding disordered material in terms of new dynamical phases.

The diversity of single spin relaxations, especially in the moderately
localized regime where the Kohlrausch law is stable, indicates that the many-
body bath appears to exhibit a broad range of behaviors to various spins –
ranging from sustained oscillatory to over-damped. Determining the
distribution of such bath characteristics might be an interesting way to
characterize an MBD system.

Investigating the Markovianity of the dynamics and information scrambling in
the Kohlrausch phase might reveal interesting facets of relaxation in MBD
systems.

Acknowledgments: The author is deeply grateful to V. Khemani and R. Moessner
for several crucial discussions at various stages of the project starting
right from its inception, as well as for their continuous encouragement and
support. The author also thanks A. Chandran, A. Das, A. Polkovnikov and T.
Prosen for many useful discussions.

##  Appendix A The Probability Distributions of $\beta,$ $\nu$ and $V$

Figure 10:  The Probability Distribution for $\beta,$ $\nu$ and $V$
corresponding to the case $\alpha=2.8,W=0.4\pi,\epsilon=0.2,L=26$ (Kohlrausch
regime). The results are obtained (for this Fig. as well as the later ones) by
averaging over 5000 disorder realizations for each data point, and for each
realization, a random $x$-bit-string state is taken as the initial state.

Here we provide the probability distributions of some of the parameters of the
non-hermitian Hamiltonians obtained from the fitting of the analytical formula
(Eq. 17) to the numerical results obtained for various disorder realizations.

##  Appendix B Ruling out Spin-Glass Mechanism as Essential in Kohlrausch
Stretching:

In order to strike out any role of frustration in the slow dynamics in our
setting, we show that taking any arbitrary bit-string state in $x$ direction
as the initial state, the dynamics with all $J_{ij}>0,$ in the interaction
part $H_{x}$ (Eq. 3), i.e., completely ferromagnetic interactions are exactly
identical for that in which each $J_{ij}$ has flipped sign (highly frustrated
one). Our proof subsumes exact dynamical equivalence between a system with
random ferromagnetic interactions and that with Sherrington-Kirkpatrick (SK)
spin glass interactions [56] in the above sense for each of the ($2^{L}$)
eigenstates of the $\\{\sigma^{x}_{i=1,L}\\}$ as the initial states having all
kinds of energies. For a glassy system, a low energy initial state is
generically expected to get trapped in extremely long-lived metastable states
(surrounded by lofty energy barriers) during relaxation, while in a random
ferromagnet, the dynamics do not slow down due to such long-lived glassy
traps. Here we show that the slowing down happens equally in the case of a
glassy $H_{x}$ (being completely antiferromagnetic, long-ranged, and random
[57]) as well as a purely random ferromagnetic $H_{x},$ ruling out any
essential role of the frustration in $H_{x},$ in the stretching of the
relaxation.

To show this, we consider Floquet evolution of an observable $O_{x}$ diagonal
on $x$-basis under a unitary (over a period) of the form

$\displaystyle U$ $\displaystyle=$ $\displaystyle U_{z}U_{x},~{}{\rm where}$
$\displaystyle U_{x}$ $\displaystyle=$ $\displaystyle
e^{-iH_{x}(\\{\sigma_{i}^{x}\\})},~{}U_{z}=e^{-i(\pi+\epsilon)\sum_{j}\sigma_{j}^{z}},$
(18)

and show that its dynamics exactly maps to that of $O_{x}$ evolved with

$\displaystyle\tilde{U}$ $\displaystyle=$ $\displaystyle
U_{z}\tilde{U_{x}},~{}{\rm where}$ $\displaystyle\tilde{U}_{x}$
$\displaystyle=$ $\displaystyle e^{+iH_{x}}$ (19)

(that is, with a sign-flipped $H_{x}$), provided the initial states
$|\psi(0)\rangle$ (for evolution with $U$) and $|\tilde{\psi}(0)\rangle$ (for
evolution with $\tilde{U}$) satisfy the following condition:

$[a_{\alpha}^{(0)}]^{*}=e^{-i\gamma_{\alpha}^{(0)}}\tilde{a}_{\alpha}^{(0)},$
(20)

for some $\gamma_{\alpha}^{(0)},$ where

$\displaystyle|\psi(0)\rangle$ $\displaystyle=$
$\displaystyle\sum_{\alpha}a_{\alpha}^{(0)}|x_{\alpha}\rangle~{}{\rm and}$
$\displaystyle|\tilde{\psi}(0)\rangle$ $\displaystyle=$
$\displaystyle\sum_{\alpha}\tilde{a}_{\alpha}^{(0)}|x_{\alpha}\rangle,$ (21)

Here $H_{x}$ denotes the interaction part that depends only on
$\sigma_{i}^{x}$s, and $|x_{\alpha}\rangle$ denotes the $x$-basis states. Note
that this condition includes the same initial states (e.g., same random bit-
string $|x_{\alpha}\rangle$) for the dynamics under $U$ and $\tilde{U}$.
Taking $H_{x}=-\sum_{i>j}J_{ij}/r_{ij}^{\alpha}\sigma_{i}^{x}\sigma_{j}^{x}$
with $J_{ij}>0$

We prove this by induction. We assume the states after the $n-$th complete
drive cycles are

$\displaystyle|\psi(n)\rangle$ $\displaystyle=$
$\displaystyle\sum_{\alpha}a_{\alpha}^{(n)}|x_{\alpha}\rangle~{}{\rm and}$
$\displaystyle|\tilde{\psi}(n)\rangle$ $\displaystyle=$
$\displaystyle\sum_{\alpha}\tilde{a}_{\alpha}^{(n)}|x_{\alpha}\rangle,$ (22)

with the coefficients satisfying

$[a_{\alpha}^{(n)}]^{*}=e^{-i\gamma_{\alpha}^{(n)}}\tilde{a}_{\alpha}^{(n)},$
(23)

then one can show (see appendixB ) that after the $(n+1)-$th cycle, the
coefficients (for the expansion of $|\psi(n+1)\rangle$ and
$|\tilde{\psi}(n+1)\rangle$) will satisfy

$[a_{\alpha}^{(n+1)}]^{*}=e^{-i\gamma_{\alpha}^{(n+1)}}\tilde{a}_{\alpha}^{(n+1)},$
(24)

Note that Eq. 24 automatically implies equality of the mod-squares of the
coefficients of the individual $x-$basis states between the states evolved by
$U$ and $\tilde{U}$ respectively at all time, hence the equality of $\langle
O_{x}\rangle(n)$ in both cases.

We denote

$\displaystyle U_{x}|x_{\alpha}\rangle$ $\displaystyle=$ $\displaystyle
e^{-i\phi_{\alpha}}|x_{\alpha}\rangle~{}{\rm and}$ (25) $\displaystyle
U_{z}|x_{\alpha}\rangle$ $\displaystyle=$
$\displaystyle\sum_{\beta}b^{\alpha}_{\beta}|x_{\beta}\rangle.$ (26)

Then we have

$\displaystyle|\psi(n)\rangle$ $\displaystyle=$
$\displaystyle\sum_{\alpha}a_{\alpha}^{(n)}|x_{\alpha}\rangle~{}{\rm and}$
$\displaystyle|\tilde{\psi}(n)\rangle$ $\displaystyle=$
$\displaystyle\sum_{\alpha}\tilde{a}_{\alpha}^{(n)}|x_{\alpha}\rangle$ (27)
$\displaystyle=$
$\displaystyle\sum_{\alpha}e^{i\gamma_{\alpha}^{(n)}}[\tilde{a}_{\alpha}^{(n)}]^{*}|x_{\alpha}\rangle.$

Last step is obtained using Eq. (23). Now applying $U$ on the expression of
$|\psi(n)\rangle$ above, we get

$\displaystyle|\psi(n+1)\rangle$ $\displaystyle=$
$\displaystyle\sum_{beta}a^{(n+1)}_{\beta}|x_{\beta}\rangle,~{}{\rm where}$
(28) $\displaystyle a^{(n+1)}_{\beta}$ $\displaystyle=$
$\displaystyle\sum_{\alpha}a_{\alpha}^{(n)}e^{-i\phi_{\alpha}}b^{\alpha}_{\beta}$
(29)

Similarly, applying $\tilde{U}$ on $|\tilde{\psi}(n)\rangle$ in Eq. 27

$\displaystyle|\tilde{\psi}(n+1)\rangle$ $\displaystyle=$
$\displaystyle\sum_{beta}\tilde{a}^{(n+1)}_{\beta}|x_{\beta}\rangle,~{}{\rm
where}$ (30) $\displaystyle\tilde{a}^{(n+1)}_{\beta}$ $\displaystyle=$
$\displaystyle\sum_{\alpha}[a_{\alpha}^{(n)}]^{*}e^{i\phi_{\alpha}}e^{i\gamma_{\alpha}^{(n)}}b^{\alpha}_{\beta}.$
(31)

Now the strategy is to put the above expressions in the expressions (Eq. 24)
we want to prove and see if we get something correct (that is, an assignment
of $\gamma_{\alpha}$’s that satisfies it) or something absurd (we find no such
assignment possible). Eq. (24) to be true, we must have (using Eq. 31 in it),

$\displaystyle[a_{\beta}^{(n+1)}]^{*}=e^{-i\gamma_{\beta}^{(n+1)}}\sum_{\alpha}[a_{\alpha}^{(n)}]^{*}e*{i\phi_{\alpha}}e^{i\gamma_{\alpha}^{(n)}}b^{\alpha}_{\beta}.$
(32)

Comparing this with the complex conjugate of Eq. 29, we must have

$\sum_{\alpha}[a_{\alpha}^{n}]^{*}e^{i\phi_{\alpha}}[b^{\alpha}_{\beta}]^{*}=e^{-i\gamma_{\beta}^{(n+1)}}\sum_{\alpha}e^{-i\gamma_{\alpha}^{(n)}}[a_{\alpha}^{n}]^{*}e^{i\phi_{\alpha}}b^{\alpha}_{\beta}.$
(33)

Since $\phi_{\alpha}$ should depend on the specifics of $H_{x},$ in order to
have a condition independent of that, we must have term-by-term equality for
the above sum. That is, we must have (canceling out the common factors)

$e^{i(\gamma_{\beta}(n)-\gamma_{\alpha}(n))}=\frac{b^{\alpha}_{\beta}}{[b^{\alpha}_{\beta}]^{*}}.$
(34)

At this stage, it is clear that this will be a feasible condition with
$\gamma_{\alpha}$s independent of $n,$ so we drop $n-$dependence of
$\gamma_{\alpha}$s.

Now,

$b^{\alpha}_{\beta}=\langle x_{\alpha}|U_{z}|x_{\beta}\rangle.$ (35)

Let us denote by $n_{\alpha\beta},$ the number of spins needed to be flipped
to generate $|x_{\alpha}\rangle$ from $|x_{\beta}\rangle$. Then from the
Taylor expansion of $U_{z},$ it can be shown that when $n_{\alpha\beta}$ is
odd, then $b^{\alpha}_{\beta}$ is imaginary hence
$\frac{b^{\alpha}_{\beta}}{[b^{\alpha}_{\beta}]^{*}}=-1$ and when
$n_{\alpha\beta}$ is even, then $b^{\alpha}_{\beta}$ is real, hence
$\frac{b^{\alpha}_{\beta}}{[b^{\alpha}_{\beta}]^{*}}=+1.$ Thus, we can write

$\frac{b^{\alpha}_{\beta}}{[b^{\alpha}_{\beta}]^{*}}=e^{i\pi
n_{\alpha\beta}}=e^{i(n_{\alpha}-n_{\beta})}.$ (36)

The proof of the last step of the above equation (the assignment) follows this
part of the proof. This is equivalent for the present purpose to the
individual assignments of $\gamma_{\alpha}$s as

$\gamma_{\alpha}=n_{\alpha}\pi,$ (37)

where $n_{\alpha}=$ the number of spins required to be flipped to create
$|x_{\alpha}\rangle$ from some pre-assigned reference state $|R\rangle$ (e.g.,
all-up $x-$basis state). This gives the concrete form of the necessary
starting condition for the induction proof (Eq. 23, 24). The condition now
reads: if after $n-$th cycles of evolution by $U$ and $\tilde{U}$ the
coefficients satisfy the condition

$[a_{\alpha}^{(n)}]^{*}=e^{-i\pi n_{\beta}}\tilde{a}_{\alpha}^{(n)},$ (38)

then the same must hold after $(n+1)-$th step, i.e, we must have

$[a_{\alpha}^{(n+1)}]^{*}=e^{-i\pi n_{\beta}}\tilde{a}_{\alpha}^{(n)},$ (39)

This statement can be proved immediately following the steps sketched above.
Now it is easy to see that this satisfied if we take
$|\psi(0)\rangle=|x_{\alpha}\rangle$ (any $x-$basis state). This completes the
proof.

Proof of the Assignment: Here we prove that $e^{i\pi
n_{\alpha\beta}}=e^{i\pi(n_{\alpha}-n_{\beta})}.$

Let $n_{\alpha}=$ the number of spins required to be flipped to create
$|x_{\alpha}\rangle$ from the all-up $x-$basis state, and $n_{\beta}$ is
defined likewise.

Let $n_{c}=$ the number of common spins that are flipped in both
$|x_{\alpha}\rangle$ and $|x_{\beta}\rangle$ compared to the reference state
$|R\rangle$.

Then $\tilde{n}_{\alpha}=n_{\alpha}-n_{c}=$ No. of spins that are flipped in
$|x_{\alpha}\rangle$ (w.r.t $|R\rangle$) but not flipped in
$|x_{\beta}\rangle.$

Similarly, $\tilde{n}_{\beta}=n_{\beta}-n_{c}=$No. of spins that are flipped
in $|x_{\beta}\rangle$ but not flipped in $|x_{\alpha}\rangle.$

Thus $n_{\alpha\beta}=\tilde{n}_{\alpha}+\tilde{n}_{\beta}.$

Now, all we care about is the parity (whether odd/even) of
$P(n_{\alpha\beta}).$ But from above we have

$n_{\alpha\beta}=\tilde{n}_{\alpha}+\tilde{n}_{\beta}=n_{\alpha}+n_{\beta}-2n_{c}.$

Since $2n_{c}$ is always even, we hence have

$\displaystyle P(n_{\alpha\beta})=P(\tilde{n}_{\alpha}+\tilde{n}_{\beta})$
$\displaystyle=$ $\displaystyle P(n_{\alpha}+n_{\beta})$

But it for any integers $n_{\alpha}$ and $n_{\beta},$ it is evident that
$P(n_{\alpha}+n_{\beta})=P(n_{\alpha}-n_{\beta}).$ Thus we finally have

$\displaystyle P(n_{\alpha\beta})$ $\displaystyle=$ $\displaystyle
P(n_{\alpha}-n_{\beta})$ $\displaystyle\implies e^{i\pi n_{\alpha\beta}}$
$\displaystyle=$ $\displaystyle
e^{i\pi(n_{\alpha}-n_{\beta})}~{}{\rm(proved)}.$ (40)

## References

  * Basko  _et al._ [2006] D. Basko, I. Aleiner, and B. Altshuler, Annals of Physics 321, 1126 (2006).
  * Gornyi  _et al._ [2005] I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Phys. Rev. Lett. 95, 206603 (2005).
  * Pal and Huse [2010] A. Pal and D. A. Huse, Phys. Rev. B 82, 174411 (2010).
  * Huse  _et al._ [2014] D. A. Huse, R. Nandkishore, and V. Oganesyan, Phys. Rev. B 90, 174202 (2014).
  * Imbrie  _et al._ [2017] J. Z. Imbrie, V. Ros, and A. Scardicchio, Annalen der Physik 529, 1600278 (2017).
  * Abanin and Papić [2017] D. A. Abanin and Z. Papić , Annalen der Physik 529, 1700169 (2017).
  * Lazarides  _et al._ [2015] A. Lazarides, A. Das, and R. Moessner, Phys. Rev. Lett. 115, 030402 (2015).
  * Ponte  _et al._ [2015] P. Ponte, Z. Papić, F. Huveneers, and D. A. Abanin, Phys. Rev. Lett. 114, 140401 (2015).
  * Khemani  _et al._ [2019] V. Khemani, R. Moessner, and S. L. Sondhi, A brief history of time crystals (2019).
  * Else  _et al._ [2020] D. V. Else, C. Monroe, C. Nayak, and N. Y. Yao, Annual Review of Condensed Matter Physics 11, 467 (2020), https://doi.org/10.1146/annurev-conmatphys-031119-050658 .
  * Haldar and Das [2022] A. Haldar and A. Das, Journal of Physics: Condensed Matter 34, 234001 (2022).
  * 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, Annals of Physics 427, 168415 (2021).
  * Panda  _et al._ [2020] R. K. Panda, A. Scardicchio, M. Schulz, S. R. Taylor, and M. Žnidarič, EPL (Europhysics Letters) 128, 67003 (2020).
  * Šuntajs  _et al._ [2020a] J. Šuntajs, J. Bonča, T. c. v. Prosen, and L. Vidmar, Phys. Rev. E 102, 062144 (2020a).
  * Šuntajs  _et al._ [2020b] J. Šuntajs, J. Bonča, T. c. v. Prosen, and L. Vidmar, Phys. Rev. B 102, 064207 (2020b).
  * Kiefer-Emmanouilidis  _et al._ [2021] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Phys. Rev. B 103, 024203 (2021).
  * Sels and Polkovnikov [2021a] D. Sels and A. Polkovnikov, Thermalization of dilute impurities in one dimensional spin chains (2021a), arXiv:2105.09348 [quant-ph] .
  * Morningstar  _et al._ [2022] A. Morningstar, L. Colmenarez, V. Khemani, D. J. Luitz, and D. A. Huse, Phys. Rev. B 105, 174205 (2022).
  * Jana  _et al._ [2022] A. Jana, V. R. Chandra, and A. Garg, Universal properties of single particle excitations across the many-body localization transition (2022).
  * Choi  _et al._ [2017] S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, C. von Keyserlingk, N. Y. Yao, E. Demler, and M. D. Lukin, Nature 543, 221 (2017).
  * Zhang  _et al._ [2017] J. Zhang, P. W. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I.-D. Potirniche, A. C. Potter, A. Vishwanath, N. Y. Yao, and C. Monroe, Nature 543, 217 (2017).
  * Luitz  _et al._ [2020] D. J. Luitz, R. Moessner, S. L. Sondhi, and V. Khemani, Phys. Rev. X 10, 021046 (2020).
  * Ippoliti  _et al._ [2021] M. Ippoliti, K. Kechedzhi, R. Moessner, S. Sondhi, and V. Khemani, PRX Quantum 2, 030346 (2021).
  * Mi et. al. [2021] X. Mi et. al., Nature 601, 531–536 (2021).
  * Kohlrausch [1847] R. Kohlrausch, Ann. Phys. (Leipzig) 12, 393 (1847).
  * Phillips [1996] J. C. Phillips, Reports on Progress in Physics 59, 1133 (1996).
  * Binder and Young [1986] K. Binder and A. P. Young, Rev. Mod. Phys. 58, 801 (1986).
  * Palmer  _et al._ [1984] R. G. Palmer, D. L. Stein, E. Abrahams, and P. W. Anderson, Phys. Rev. Lett. 53, 958 (1984).
  * Spohn [1989] H. Spohn, Communications in Mathematical Physics 125, 3 (1989).
  * Kuwahara  _et al._ [2016] T. Kuwahara, T. Mori, and K. Saito, Annals of Physics 367, 96 (2016).
  * Mori  _et al._ [2016] T. Mori, T. Kuwahara, and K. Saito, Phys. Rev. Lett. 116, 120401 (2016).
  * Abanin  _et al._ [2017a] D. A. Abanin, W. De Roeck, W. W. Ho, and F. m. c. Huveneers, Phys. Rev. B 95, 014112 (2017a).
  * Abanin  _et al._ [2017b] D. Abanin, W. De Roeck, W. W. Ho, and F. Huveneer, Commun. Math. Phys. 354, 809 (2017b).
  * Else  _et al._ [2017] D. V. Else, B. Bauer, and C. Nayak, Phys. Rev. X 7, 011026 (2017).
  * Khemani  _et al._ [2016] V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. 116, 250401 (2016).
  * Ho  _et al._ [2017] W. W. Ho, S. Choi, M. D. Lukin, and D. A. Abanin, Phys. Rev. Lett. 119, 010602 (2017).
  * Saffman  _et al._ [2010] M. Saffman, T. G. Walker, and K. Mølmer, Rev. Mod. Phys. 82, 2313 (2010).
  * SchauÃŸ  _et al._ [2012] P. SchauÃŸ, M. Cheneau, M. Endres, T. Fukuhara, S. Hild, A. Omran, T. Pohl, C. Gross, S. Kuhr, and I. Bloch, Nature 491, 87 (2012).
  * Childress  _et al._ [2013] L. Childress, M. V. Gurudev Dutt, J. M. Taylora, A. S. Zibrovf, F. Jelenzko, J. Wrachtrup, P. R. Hemmer, and M. D. Lukin, Science 314, 281 (2013).
  * Islam  _et al._ [2013] R. Islam, C. Senkow, W. C. Campbell, S. Korenblit, J. Smith, A. LEE, E. E. Edward, C.-C. J. Wang, J. K. Freericks, and C. Monroe, Science 340, 583 (2013).
  * Hung  _et al._ [2016] C.-L. Hung, A. Gonzalez-Tudela, J. I. Cirac, and H. J. Kimble, Proc. Natl. Acad. Sci. U.S.A 113, E4946 (2016).
  * Prosen [1998] T. c. v. Prosen, Phys. Rev. Lett. 80, 1808 (1998).
  * Lezama  _et al._ [2019] T. L. M. Lezama, S. Bera, and J. H. Bardarson, Phys. Rev. B 99, 161106 (2019).
  * Lezama  _et al._ [2021] T. L. M. Lezama, E. J. Torres-Herrera, F. Pérez-Bernal, Y. Bar Lev, and L. F. Santos, Phys. Rev. B 104, 085117 (2021).
  * Luitz and Lev [2017] D. J. Luitz and Y. B. Lev, Annalen der Physik 529, 1600350 (2017), https://onlinelibrary.wiley.com/doi/pdf/10.1002/andp.201600350 .
  * Bar Lev and Reichman [2014] Y. Bar Lev and D. R. Reichman, Phys. Rev. B 89, 220201 (2014).
  * Sels and Polkovnikov [2021b] D. Sels and A. Polkovnikov, Phys. Rev. E 104, 054105 (2021b).
  * Kjäll  _et al._ [2014] J. A. Kjäll, J. H. Bardarson, and F. Pollmann, Phys. Rev. Lett. 113, 107204 (2014).
  * Bardarson  _et al._ [2012] J. H. Bardarson, F. Pollmann, and J. E. Moore, Phys. Rev. Lett. 109, 017202 (2012).
  * Segri and Zloshchastiev [2013] A. Segri and K. G. Zloshchastiev, International Journal of Modern Physics B 27, 1350163 (2013), https://doi.org/10.1142/S0217979213501634 .
  * Metzler and Klafter [2002] R. Metzler and J. Klafter, Journal of Non-Crystalline Solids 305, 81–87 (2002).
  * Bhakuni  _et al._ [2021] D. S. Bhakuni, L. F. Santos, and Y. B. Lev, Phys. Rev. B 104, L140301 (2021).
  * [53] H. Rieger and A. P. Young, in  _Complex Behaviour of Glassy Systems_ (Springer Berlin Heidelberg) pp. 256–265.
  * Mezard  _et al._ [1986] M. Mezard, G. Parisi, and M. Virasoro,  _Spin Glass Theory and Beyond_ (WORLD SCIENTIFIC, 1986) https://www.worldscientific.com/doi/pdf/10.1142/0271 .
  * Suzuki  _et al._ [2013] S. Suzuki, J. Inoue, and B. K. Chakrabarti,  _Quantum Ising Phases and Transitions in Transverse Ising Models_ , Vol. 862 (Lecture Notes in Physics, Springer, 2013).
  * Sherrington and Kirkpatrick [1975] D. Sherrington and S. Kirkpatrick, Phys. Rev. Lett. 35, 1792 (1975).
  * Chandra  _et al._ [2009] A. K. Chandra, J. ichi Inoue, and B. K. Chakrabarti, Journal of Physics: Conference Series 143, 012013 (2009).