# Many-body localization in the infinite-interaction limit and
the discontinuous eigenstate phase transition

Chun Chen<EMAIL_ADDRESS>School of Physics and Astronomy, Key
Laboratory of Artificial Structures and Quantum Control (Ministry of
Education), Shenyang National Laboratory for Materials Science, Shanghai Jiao
Tong University, Shanghai 200240, China  Max-Planck-Institut für Physik
komplexer Systeme, Nöthnitzer Straße 38, 01187 Dresden, Germany  Yan Chen
<EMAIL_ADDRESS>Department of Physics and State Key Laboratory of
Surface Physics, Fudan University, Shanghai 200433, China  Collaborative
Innovation Center of Advanced Microstructures, Nanjing University, Nanjing
210093, China  Xiaoqun Wang<EMAIL_ADDRESS>School of Physics and
Astronomy, Key Laboratory of Artificial Structures and Quantum Control
(Ministry of Education), Shenyang National Laboratory for Materials Science,
Shanghai Jiao Tong University, Shanghai 200240, China  Tsung-Dao Lee
Institute, Shanghai Jiao Tong University, Shanghai 200240, China  Beijing
Computational Science Research Center, Beijing 100084, China

###### Abstract

Can localization persist when interaction grows infinitely stronger than
randomness? If so, is it many-body Anderson localization? How about the
associated localization transition in the infinite-interaction limit? To
tackle these questions, we study many-body localization (MBL) in a spin-chain
model mimicking the Rydberg-blockade quantum simulator with both infinite-
strength projection and moderate quasiperiodic modulation. Employing exact
diagonalization, Krylov-typicality technique, and time-evolving block
decimation, we identify evidence for a constrained MBL phase stabilized by a
pure quasirandom transverse field. Remarkably, the constrained MBL transition
may embody a discontinuous eigenstate phase transition, whose discontinuity
nature significantly suppresses the finite-size drifts that plague most
numerical studies of conventional MBL transition. Through quantum dynamics, we
find that rotating the modulated field from parallel toward perpendicular to
the projection axis induces an eigenstate transition between the diagonal and
constrained MBL phases. Intriguingly, the entanglement-entropy growth in
constrained MBL follows a double-log form, whereas it changes to a power law
in approaching the diagonal limit. By unveiling the significance of confined
nonlocal effects in integrals of motion of constrained MBL, we show that this
newfound insulating state is not a many-body Anderson insulator. Our
predictions can be tested in Rydberg experiments.

INTRODUCTION

The framework of many-body localization lays its foundation on noninteracting
Anderson insulator [1] to address the quest of ergodicity breaking [2, 3] and
instability toward delocalization and eigenstate thermalization [4, 5] under
weak many-body interactions in low spatial dimensions [6, 7].

This short-range weak-interaction picture forms the backbone of conventional
MBL. Nonetheless, it also raises a question of whether there can arise many-
body non-Anderson localization in circumstances where interaction strengths
are not weak but infinitely strong. See [8, 9, 10, 11, 12, 13, 14, 15, 16] for
different considerations on uniform or random-interaction systems.

Phenomenologically, isolated many-body Anderson insulators may be describable
by the emergent extensive set of local integrals of motion (LIOMs or
$\ell$-bits) [17, 18, 19, 20], at least in one dimension ($1$D) [21]. Then, is
it conceivable that localization persists but owing to restriction or
frustration, the LIOM-based picture breaks down? It is known that finite
interaction activates more resonance channels for dephasing, so it is expected
to suppress localization. In this regard, a better and affirmative route to
achieving the unconventional MBL might be associated with the presence of
restriction or frustration, for instance, in disordered Rydberg-blockade
chains [22, 23], where, as a consequence of strong, coherent dipole-dipole van
der Waals repulsions, two nearest-neighbouring Rydberg atoms cannot be
simultaneously excited. This energy constrained dynamics is modelled by a
projection action of infinite strength.

Specifically, would there be a singular boundary separating different phases
of MBL due to abrupt distortion rather than a progressive dressing of
$\ell$-bits?

MBL phase and MBL transition are two interdependent topics central to
ergodicity breaking in statistical mechanics. Recently, there is a debate
questioning MBL as a viable state in the thermodynamic limit. The issue stems
from the strong finite-size drift of the MBL-thermal phase boundary seen in
nearly all numerical scaling analyses of small chains [24, 25, 26, 27, 28, 29,
30, 31, 32, 33, 34, 35, 36, 37, 38]. Because the critical disorder strength
keeps shifting toward infinity under the increase of system’s size, it was
inferred that no MBL transition occurs within these models such that MBL might
be a finite-size crossover phenomenon that ultimately gives way to the normal
process of thermalization.

How about the nature of MBL transition in the presence of infinite
interparticle interaction?

Counterintuitively, we find through solving a concrete lattice model that the
eigenstate transition between MBL and thermal regimes may contrastingly be
discontinuous when (off-diagonal) constrained limit is taken, a feature
probably enabled by the infinite interaction that considerably reduces the
adverse effects of finite-size drifts at transition points, thereby
strengthening the stabilization of unconventional constrained MBL state and,
as a byproduct, the robustness of diagonal MBL phase.

RESULTS

The minimal model

The aforementioned physics might be visible in disordered and locally
constrained chain models [23]. The simplest of such category takes the
following archetypal form,

$H_{\mathrm{qp}}=\sum_{i}\left(g_{i}\widetilde{X}_{i}+h_{i}\widetilde{Z}_{i}\right),$
(1)

where $\widetilde{X}_{i},\widetilde{Z}_{i}$ are projected Pauli matrices,
$\widetilde{X}_{i}\coloneqq P\sigma^{x}_{i}P$ and $\widetilde{Z}_{i}\coloneqq
P\sigma^{z}_{i}P$. The global operator
$P\coloneqq\prod_{i}(\frac{3+\sigma^{z}_{i}+\sigma^{z}_{i+1}-\sigma^{z}_{i}\sigma^{z}_{i+1}}{4})$
prohibits motifs of $\downarrow\downarrow$-configuration over any adjacent
sites, hence rendering the Hilbert space of model (1) locally constrained.

In Ref. [23], we showed that a random version of model (1) by quenched
disorder exhibits tentative signatures of a constrained MBL (cMBL) phase;
nevertheless, as being in proximity to the nearby transition, the Griffiths
effect therein proliferates, which impedes an identification and thus a direct
investigation of this unconventional nonergodic state of matter. In current
work, we improve our prior construction by conceiving an experiment-pertinent
quasiperiodic constrained model with open and periodic boundary conditions
(BCs), i.e., choosing [39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52,
53, 54, 55]

$g_{i}=g_{x}+W_{x}\cos(2\pi\xi i+\phi_{x}),\ \ h_{i}=W_{z}\cos(2\pi\xi
i+\phi_{z}),$ (2)

where the wavenumber $\xi=\sqrt{2}$ is irrational, $i=1,\ldots,L$, and
$\phi_{x},\phi_{z}\in[-\pi,\pi)$ are different sample-dependent random overall
phase shifts.

Since Hamiltonian (1) is real, time-reversal symmetry ${\sf T}\coloneqq K$ is
preserved, giving rise to the Gaussian orthogonal ensemble (GOE) in phases
obeying the eigenstate thermalization hypothesis (ETH) [56]. Additionally,
when $W_{z}=0$ there is a particle-hole symmetry ${\sf
P}\coloneqq\prod_{i}\sigma^{z}_{i}$ that anticommutes with $H_{\mathrm{qp}}$.
To our knowledge, no discrete Abelian symmetry is present in Hamiltonian (1),
so the possibility of a localization-protected spontaneous symmetry breaking
[57] is excluded.

To manifest the fundamental interplay between finite tunable randomness and
infinite interparticle interaction, we introduce hardcore boson operators
$b^{\dagger},b$ on each site to describe the local pseudospin-$\frac{1}{2}$
system that emulates the Rydberg lattice gas with ground state
$|g\rangle=|\\!\uparrow\rangle$ and Rydberg excitation state
$|r\rangle=|\\!\downarrow\rangle$. In terms of hardcore bosons, the Pauli spin
matrices can be couched as follows,

$\displaystyle b^{\dagger}+b=|r\rangle\langle g|+|g\rangle\langle
r|=|\\!\downarrow\rangle\langle\uparrow\\!|+|\\!\uparrow\rangle\langle\downarrow\\!|=\sigma^{x},$
(3) $\displaystyle b^{\dagger}b=n=|r\rangle\langle
r|=|\\!\downarrow\rangle\langle\downarrow\\!|=\frac{(1-\sigma^{z})}{2},$ (4)

where $n=0,1$ is the local occupation number of boson. Armed with the above
expressions, Hamiltonian (1) can then be exactly mapped onto an array of
neutral atoms in the Rydberg blockade regime,

$\displaystyle H_{\textrm{qp}}$ $\displaystyle=H_{x}+H_{z}+H_{V},$ (5)
$\displaystyle H_{x}$ $\displaystyle=\sum_{i}g_{i}(b^{\dagger}_{i}+b_{i}),$
(6) $\displaystyle H_{z}$ $\displaystyle=\sum_{i}h_{i}(1-2n_{i}),$ (7)
$\displaystyle H_{V}$ $\displaystyle=\sum_{i}V_{1}n_{i}n_{i+1},\ \ \
V_{1}=\infty.$ (8)

Here $g_{i},h_{i}$ are proportional to onsite Rabi frequency and frequency
detuning, respectively; the repulsive van der Waals interaction in Eq. (8) is
truncated to retain only the nearest-neighbour interaction whose strength
$V_{1}$ is lifted to infinity, producing a blockade radius of $a<R_{b}<2a$.
Clearly, $H_{x}$ breaks system’s particle-number conservation, so the total
energy is the only conserved quantity of the model.

Instructively, using spin operators in Eqs. (3) and (4), Hamiltonian (5) can
also be recast into the standard mixed-field Ising model, for which Imbrie
[21] proved in a mathematically quasiexact way the existence of many-body-
generalized Anderson insulator under conditions of limited level attraction,
weak interaction strengths, and sufficiently strong disorders. In this regard,
although we focus on the infinitely interacting version of such a particular
model, the gained results bear the originality and significance to stimulate
the research about many-body non-Anderson localization. Further, without
$H_{V}$, $H_{0}=H_{x}+H_{z}$ is a free Hamiltonian describing decoupled spins,
each undergoing an independent Larmor precession about the local random
fields.

Therefore, the constrained Rydberg atomic chain we consider consists of two
pieces: a randomized but noninteracting term $H_{0}$ and a nearest-neighbour
density-density interacting term $H_{V}$ featured by an infinite repulsion.

Such a compact form with a single “spin-like” sector and the reduction of
onsite Hilbert-space dimension from the usual value of $2$ to the golden ratio
$\phi=1.618\ldots$ prompt us to regard the bare bones model (5) [or Eq. (1)]
as the fundamental building block for studying the more generic constrained
quantum systems, such as the $t$-$J$ model.

Figure 1:  Spectral diagnostics of the cMBL-cETH transition based on the
scaling analysis of $[r]$ using OBCs. The maximal chain length is pushed to
$L=26$. Under constrained limit $W_{z}=0$, the model’s $[r]$-value forms a
systematic trend in approaching $r_{\textrm{Poi}}$ when
$\frac{g_{x}}{W_{x}}\in(0.5,1]$, demonstrating the realization of cMBL. Panels
(a) and (b) show the proposed respective discontinuity for the left and right
boundary of the transition. Panel (c) schematizes the ideal phase diagram
along the (off-diagonal) constrained line in the thermodynamic limit where, as
denoted by black dots, the chain at the left discontinuous transition point
$g_{x}=0.5$ is thermalized whereas it remains fully localized at the right
discontinuous transition point $g_{x}=1$. This cMBL phase survives to finite
$\frac{W_{z}}{W_{x}}\approx 0.5$, hence forming a dome separated from both the
constrained thermal phase at leading $g_{x}$ and the dMBL state at dominant
$W_{z}$.

Moreover, in light of the following commutation relations,

$[H_{x},H_{V}]\neq 0,\ \ \ [H_{z},H_{V}]=0,$ (9)

the constrained Hamiltonian (5) may accommodate two distinct physical
extremes. (i) When $|W_{z}|\gg|g_{x}|,|W_{x}|$, the system approaches the
diagonal limit during which the role of infinite interaction is effectively
minimized and the resulting diagonal MBL (dMBL) state represents a variant of
many-body Anderson insulator with enhanced robustness [58]. To be pedantic,
throughout this paper, we define the diagonal limit, an analog of Anderson
limit, as specified by $h_{i}\neq 0$ and $g_{i}=0$ in (1); while for dMBL, the
analog of many-body Anderson localization, $|g_{i}|$, although perturbatively
smaller than $|h_{i}|$, is not identically zero. (ii) In comparison, once
$|g_{x}|,|W_{x}|\gg|W_{z}|$, the system enters the off-diagonal constrained
limit—the true “infinite-interaction limit” quoted in the paper’s title—where
mutual impacts from modest randomness and infinite interaction are
contrastingly maximized. Particularly, their constructive interplay gives rise
to the sought infinite-interaction-facilitated MBL state which is different
from the “infinite-randomness-controlled” many-body Anderson localization
stabilized in the opposite limit of weak interaction. Naïvely, no apparent
duality would directly link these two.

It is worth stressing that the kinetic constraint was realized in Rydberg-
blockade chain [22] and the quasiperiodic modulation played a vital role in
experiments [59, 60, 61] to achieve the first signature of MBL in
unconstrained systems. Accordingly, the actual value of model (1) resides
right in its high experimental relevance.

Throughout this paper, $W_{x}=1$ sets the energy scale, i.e., the system is
quasirandom at least along $x$ direction.

Discontinuous cMBL-cETH transition: Spectral analyses

The configuration averaged level-spacing ratio $[r]$ is the unique single-
value quantity routinely adopted to characterize the dynamical states of
matter. One defining feature of the robust localization is the vanishing
repulsion between contiguous gaps and the resulting Poisson distribution of

$r_{n}\coloneqq\frac{\min\\{\delta_{n},\delta_{n-1}\\}}{\max\\{\delta_{n},\delta_{n-1}\\}}$
(10)

with mean $[r]=r_{\textrm{Poi}}\approx 0.386$ where $\delta_{n}\coloneqq
E_{n}-E_{n-1}$ assuming $\\{E_{n}\\}$ an ascending list of eigenvalues [62,
63, 31, 32, 64].

Figures 1(a),(b) show the finite-size evolutions of $[r]$ as a function of
$g_{x}$ along the $W_{z}=0$ axis. Via optimization of the ED algorithm
targeting only the eigenvalues, we obtain the full eigenspectra of the chain
for $1000$ independent quasirandom samples up to system size $L=26$ and the
corresponding Hilbert-space dimension for such a single sample is $317811$
under OBCs. (Parenthetically, the maximal chain length examined by a similar
work [30] is also $L=26$, as they used PBCs, the corresponding Hilbert-space
dimension increases to $271443$.) Within $0.5<g_{x}\leqslant 1$, we find that
$[r]$ steadily converges to $r_{\textrm{Poi}}$ under the successive increase
of $L$, verifying the stabilization of a desired cMBL phase.

Figure 2:  $[r]$ data at the left [column (a)] and right [column (b)]
transition points. Dashed lines in (a1),(b1) mark the lengths exceeding which
the $[r]$ values for $g_{x}=0.5,1.1$ begin to grow. By contrast, (a2),(b2)
show the monotonic decrease of $[r]$ at $g_{x}=0.6,1$ toward
$r_{\textrm{Poi}}$. All data are taken from Fig. 1 with $W_{z}=0$. (a3),(b3)
plot the $[r]$-difference ratios near the two transition points to assess the
sharpness of the discontinuity for the underlying transition.

The peculiarity of the quasiperiodic Rydberg chain may be rooted in the
discontinuous eigenstate phase transition between cMBL and its nearby
constrained ETH (cETH) phase. To pin down the potential discontinuity, we
perform a finite-size scaling analysis of $[r]$ by pushing the chain length to
$L=26$ and simultaneously selecting a compatibly high resolution $0.1W_{x}$
when tuning $g_{x}$. Figure 1(a) targets the left discontinuous transition
point, from which one observes that for $0\leqslant g_{x}\leqslant 0.5$, there
exists a turning point of $L$ beyond which $[r]$ starts to increase
continually toward GOE. More precisely, we find the following correspondence
between $g_{x}$ and its turning position of the length:

$\begin{bmatrix}g_{x}\\\ L\end{bmatrix}=\begin{bmatrix}0\\\
13\end{bmatrix};\begin{bmatrix}0.1\\\ 15\end{bmatrix};\begin{bmatrix}0.2\\\
17\end{bmatrix};\begin{bmatrix}0.3\\\ 19\end{bmatrix};\begin{bmatrix}0.4\\\
21\end{bmatrix};\begin{bmatrix}0.5\\\ 20\end{bmatrix},$ (11)

namely, the turning point of $L$ increases along with $g_{x}$ up to
$g_{x}=0.5$. Surprisingly, this trend terminates abruptly once $g_{x}>0.5$.
For instance, even after moving upward to $L=26$, we find no signature of such
a turning length for $g_{x}=0.6$ and its $[r]$-value keeps rolling down toward
$r_{\textrm{Poi}}$. This continual decrease of $[r]$ becomes more transparent
for $g_{x}=0.7,0.8,0.9,1$, which suggests the realization of cMBL within
$g_{x}\in(0.5,1]$. Therefore, the sharp distinction of the two trends induced
by an incremental change of $g_{x}$ indicates that the assumed discontinuity
of the transition occurs at $g_{x}=0.5$.

The discontinuity between the two contrasting trends showcases more vigorously
at the right discontinuous transition point. As presented by Fig. 1(b),
although working on the small chain of $L=26$, the $[r]$-value for $g_{x}=1$
is already extremely close to the ideal value of the Poissonian distribution;
while, in a striking comparison, the $[r]$-value for $g_{x}=1.2$ appears to
shoot up toward GOE at the very similar length scale. See also Figs. 3(a),(b).
Crucially, the equality of $[r]$ at $g_{x}=0.9,1$ and the finite jump of $[r]$
at $g_{x}=1.1$ comprise a vivid definition of the discontinuity.

Based on the insights gained from extrapolating the scaling analysis of $[r]$,
we draw in Fig. 1(c) the schematic phase diagram of the quasirandom Rydberg
chain in the thermodynamic limit. Here we exclusively focus on the constrained
limit by fixing $W_{z}=0$. After taking random averages, the sign of $W_{x}$
makes no difference. Together with the unitary transformation rotated by
$\sigma^{z}_{i}$, the sign of $g_{x}$ does not matter either. One can thus
take $W_{x}=1$ and consider $\frac{g_{x}}{W_{x}}\in[0,+\infty)$ without loss
of generality. There are several features about the phase diagram. (i) There
only exist two eigenstate phases along the constrained line $W_{z}=0$, the
cMBL phase and the cETH phase. (ii) The cMBL phase occupies a finite interval
$g_{x}\in(0.5,1]$ and the cMBL-cETH transition is likely discontinuous whose
transition points locate at the phase boundaries $g_{x}=\frac{1}{2}$ and
$g_{x}=1$. (iii) The system appears thermalized at $g_{x}=\frac{1}{2}$; while
on $g_{x}=1$, the chain maintains its full localization character. (iv) At the
special point $g_{x}=0$, where the relative randomness strength is infinite
$W_{x}/g_{x}=\infty$, the chain is well within the cETH phase.

The robustness of cMBL implies the existence of dMBL in the thermodynamic
limit, because adding another source of randomness can only enhance
localization. In this sense, the constraint-induced delocalization reported by
Ref. [30] does not lead to contradictions but rather highlights the importance
of cMBL as the state to foster dMBL. It also becomes clear about the necessity
to include both terms of $g_{x}$ and $W_{x}$ on an equal footing to achieve
the stabilization of localization in a general constrained setting.
Nonetheless, this type of randomness encapsulated by the term $g_{i}$ is not
the topic of [30].

Transition points refined.—We examine the postulated discontinuity of $[r]$ at
the transition points a bit further in Fig. 2, where by zooming in the left
[panel (a1)] and right [panel (b1)] transition zones, we highlight the
opposite scaling trends of $[r]$ for the two $g_{x}$’s that are close in
magnitude. In accordance with the emergent integrability of cMBL, Figs.
2(a2),(b2) show up to $L=26$ the decrease and the convergence of $[r]$ onto
the Poisson value at $g_{x}=0.6,1$. The degree of discontinuity of the
transition may be quantified in terms of the $[r]$-difference ratio defined by
comparing the $[r]$-values at three adjacent $g_{x}$’s, viz., with a
particular $L$, for ordered $g_{x}=g_{1};g_{2};g_{3}$,

$\displaystyle[r]$ -difference ratio
$\displaystyle\coloneqq\frac{\textrm{max}\\!\left\\{\big{|}[r](g_{1})-[r](g_{2})\big{|},\big{|}[r](g_{2})-[r](g_{3})\big{|}\right\\}}{\textrm{min}\\!\left\\{\big{|}[r](g_{1})-[r](g_{2})\big{|},\big{|}[r](g_{2})-[r](g_{3})\big{|}\right\\}}.$
(12)

For continuous transitions, the $[r]$-difference ratio is on the order of $1$;
whereas, if discontinuity arises, then it is predicted to diverge right at the
transition once the thermodynamic limit is taken. By choosing four
$g_{x}$-tuples involving the two transition points, we illustrate via Figs.
2(a3),(b3) the growth of $[r]$-difference ratio above unity under the increase
of $L$. Within the system sizes we probe, the degree of discontinuity on the
right transition point appears stronger than that of the left one. This is
attributed to the different eigenstate phases realized at these two transition
points.

Figure 3:  (a): $[r]$-values of the chain for the two parameter sets
$(g_{x}=0,W_{z}=0)$ and $(g_{x}=1.2,W_{z}=0)$. Both grow promptly toward the
value of GOE with $L$, suggesting the system at both parametric points obeys
ETH in the thermodynamic limit. (b): The probability distribution of the
level-spacing ratio $P(r)$ in cMBL (ensured by $g_{x}=1,W_{z}=0$), which
closely traces the prediction of Poisson statistics [solid line,
$P_{\textrm{Poi.}}(r)=2/(1+r)^{2}$] under the increase of $L$, indicating the
realization of a full MBL. By selecting $g_{x}=0,W_{z}=0$, the $P(r)$
distribution switches to follow the Wigner surmise [dashed line,
$P_{\textrm{GOE}}(r)=(27/4)(r+r^{2})/(1+r+r^{2})^{5/2}$], signalling the
thermalization sets in. Figure 4:  Static diagnostics of the continuous cMBL-
dMBL transition using the variance of $S_{\textrm{vN}}$ under OBCs. (a),(b):
Finite-size scaling analyses show that at fixed $\frac{g_{x}}{W_{x}}=1$, $[r]$
and $[S_{\textrm{vN}}]/S_{\textrm{P}}$ stay around $r_{\textrm{Poi}}$ and $0$
under the adjustment of $W_{z}$ from cMBL toward dMBL. No discontinuity is
observed. (c1): The intrasample deviation of the entanglement entropy across
the entirety of the eigenspectrum,
$[\Delta_{E}[S_{\textrm{vN}}]_{c}]_{s}/S_{\textrm{P}}$. The accompanying cut-
to-cut and sample-to-sample parsings of this entropy deviation are plotted by
(c2) and (c3), respectively. The peak forming at $\frac{W_{z}}{W_{x}}=1$ in
(c3) points toward a transition between cMBL (shaded) and dMBL.

Purely random point.—On phase diagram Fig. 1(c), the point
$(\frac{g_{x}}{W_{x}}=0,\frac{W_{z}}{W_{x}}=0)$ bears multiple physical
meanings. (i) It represents the purely random limit. Fig. 3(a) evidences that
in the presence of infinite interaction, the purely random Rydberg chain
likely reaches full thermalization under the constrained limit with no
signature of localization. This contrasts to the full localization at
$(\frac{g_{x}}{W_{x}}=1,\frac{W_{z}}{W_{x}}=0)$ as is displayed by Fig. 3(b).
(ii) An implication of (i) pertains to the probable connection between the
existence of cMBL and the discontinuity of the cMBL-cETH transition. This is
because if assume the cMBL transition is continuous, then the inevitable
finite-size drifts would dominate and consequently cast cMBL into doubt. (iii)
When $W_{z}=0$, by implementing the unitary symmetry involving
$\sigma^{z}_{i}$ to alter the sign of $g_{i}$ combined with an energy-scale
redefinition, one can exactly map
$(\frac{g_{x}}{W_{x}}=0,\frac{W_{z}}{W_{x}}=0)$ onto
$(\frac{g_{x}}{W_{x}}=1,\frac{W_{z}}{W_{x}}=0)$; however, this mapping is
valid if the randomness is quench disorder and obeys the uniform box
distribution [23]. For quasiperiodic randomness, such a formal equivalence
breaks down, but our numerical data hint that
$(\frac{g_{x}}{W_{x}}=0,\frac{W_{z}}{W_{x}}=0)$ might still be equivalent to a
point infinitesimally close to but different from
$(\frac{g_{x}}{W_{x}}=1,\frac{W_{z}}{W_{x}}=0)$. Symbolically, it reads
$(\frac{g_{x}}{W_{x}}=1+0^{+},\frac{W_{z}}{W_{x}}=0)$. This reasoning implies
that the two disjoint cETH regions on the left and right side of the cMBL
phase are physically equivalent, and the robustness of cMBL may thus be
cemented by the significant degree of discontinuity of the transition at
$g_{x}=1$. In parallel, one can hypothesize that the phase diagram for the
quench disordered Rydberg chain is identical to Fig. 1(c) except that the
black dot at $g_{x}=1$ moves from Poisson to GOE. This conjecture overlooks
the Griffiths rare-region effects and the ensuing avalanche-driven
delocalization [65] in disordered models which can be detrimental to cMBL
phase and transition alike.

Continuous cMBL-dMBL transition: Entanglement variances

The bipartite entanglement entropy $[S_{\textrm{vN}}]$ is another useful proxy
for analyzing MBL as is $[r]$ [31, 32, 66]. For each eigenstate
$|\psi_{n}\rangle$, the half-chain von Neumann entropy is defined by
$S_{\textrm{vN}}\coloneqq-{\sf Tr}_{R}\left[\rho_{R}\log_{2}\rho_{R}\right]$
where $\rho_{R}\coloneqq{\sf Tr}_{L}[|\psi_{n}\rangle\langle\psi_{n}|]$ is the
reduced density matrix of the right half chain. Using ED, we compute
$S_{\textrm{vN}}$ for the entire eigenspectrum of a given sample and obtain
$[S_{\textrm{vN}}]$ after averaging over all available eigenstates for more
than $1000$ independent quasirandom realizations.

Figures 4(a),(b) show the respective evolutions of $[r]$ and
$[S_{\textrm{vN}}]/S_{\textrm{P}}$ as a function of $W_{z}$ along the
$g_{x}=1$ axis. Upon increasing $L$, both $[r]$ and
$[S_{\textrm{vN}}]/S_{\textrm{P}}$ converge to a flattening curve centring
around $r_{\textrm{Poi}}$ and $0$, respectively, demonstrating the system
retains full localization as $W_{z}$ varies from cMBL to dMBL. Because no
discontinuity is found in $[r]$ and $[S_{\textrm{vN}}]/S_{\textrm{P}}$, the
change of the phase structure is expected to be continuous. Here
$[S_{\textrm{vN}}]$ is normalized with respect to the Page value of the
thermal entropy [67] whose estimate suitable for constrained circumstance
equals $S_{\textrm{P}}\approx\log_{2}(F_{L/2+2})-1/(2\ln 2)$ where the half-
chain Hilbert-space dimension equals the Fibonacci number $F_{L/2+2}$ [23].

Differing in entanglement patterns, continuous transition between
unconstrained MBL and ETH phases can be probed via the standard deviation of
$S_{\textrm{vN}}$. As first demonstrated by Refs. [68, 69], the explicit
dependence of $S_{\textrm{vN}}$ on eigenstate wavefunction (“$E$”), randomized
sample (“$s$”), and partition cut (“$c$”) brings about three measures for the
quantity. In the following, we adopt this strategy for the continuous cMBL-
dMBL transition. We borrow the convention of Ref. [23], viz., using
$[\ldots]_{E/s/c}$ to denote the respective averages over the eigenspectrum
entirety, all samples, and all cuts with the unspecified subscripts holding
fixed. Similar definitions carry over to the standard deviations
$\Delta_{E/s/c}(\ldots).$

In accord with the continual changes of $[r]$ and
$[S_{\textrm{vN}}]/S_{\textrm{P}}$ in Figs. 4(a),(b), Fig. 4(c1) illustrates
the subvolume scaling law of
$[\Delta_{E}[S_{\textrm{vN}}]_{c}]_{s}/S_{\textrm{P}}$, the state-to-state
intrasample deviation of $S_{\textrm{vN}}$, upon raising $L$ and its overall
smooth lineshape as a function of $W_{z}$. Likewise, the cut-to-cut
entanglement-entropy deviation as given by Fig. 4(c2) exhibits the
qualitatively consistent tendency for the same evolution between the two MBLs.

By comparison, albeit being less prominent for quasiperiodic arrangements, the
sample-to-sample entanglement-deviation curve presented by Fig. 4(c3) gives
the indication of an emergent peak around $\frac{W_{z}}{W_{x}}=1$, hinting
that the rise of $W_{z}$ at fixed $g_{x}=1$ may drive a continuous transition
from cMBL toward dMBL. Intuitively, although both MBLs are dominantly
constituted by area-law entangled eigenstates, on finite-length chains, the
rates of how they approach the area scaling law may differ in magnitude and
form, which potentially allows for the entropy variance across different
samples near the phase boundary. Nonetheless, in view of the fact that the
absolute value of the deviation is not pronounced and the shape of the curve
keeps flattening, it is possible that cMBL and dMBL might vaguely be
distinguishable from pure static diagnostics.

Eigenstate transition from entanglement growth

Alternatively, the qualitative difference between cMBL and dMBL can be
demonstrated from the angle of real-time evolution of entanglement. Notably,
we find an eigenstate transition between these two dynamical regimes in the
numerical quantum quench experiments.

We use two quantities, the bipartite entanglement entropy and the quantum
Fisher information (QFI). The initial state is randomly selected from the
complete basis of nonentangled product states of $\sigma^{z}_{i}$-spins that
respects the local constraint. For each $L$, we generate more than $1000$
random pairs of $(\phi_{x},\phi_{z})$ for the Hamiltonian, and for each
quasiperiodic arrangement, we let the chain evolve and calculate
$S_{\textrm{vN}}$, QFI by ED and TEBD [70] before averaging (see methods
section).

Figure 5(a) compiles time evolutions of $[S_{\textrm{vN}}]$ along the cut
$g_{x}=1$ with ascending $W_{z}$ at $L=20$ in a log-log format. The salient
feature there is the qualitative functional change in the time-evolution
profiles. This eigenstate transition is elaborated by Figs. 5(c) and (e) where
we focus on the entanglement growth deep inside cMBL and dMBL, respectively.
For concreteness, after a transient period $t\lesssim 1$ of the initial
development, $[S_{\textrm{vN}}]$ in dMBL grows steadily as a power law of $t$
[with an exponent $(\approx 0.1)$] within the next prolonged window (up to
$t\approx 10^{14}$ at $L=20$) but its saturated value is far less than the
thermal entropy $S_{T}\approx\log_{2}(F_{2+L/2})-1/(2\ln 2)-0.06$ [23]. In
stark comparison, the growth of $[S_{\textrm{vN}}]$ in cMBL as displayed by
Fig. 5(c) follows a different functional form: within $10^{2}\lesssim
t\lesssim 10^{9}$ at $L=20$, the double-log function fits the entropy data
reasonably well. Moreover, the equilibrated $[S_{\textrm{vN}}]$ reaches a
subthermal value in cMBL and obeys a volume scaling law.

Figure 5:  cMBL-dMBL transition in dynamics with OBCs and fixed $g_{x}=1$. The
maximum chain length in ED is $L=20$. The top row summarizes functional
changes of the growth of $[S_{\textrm{vN}}]$ and $[f_{Q}]$ as a function of
$W_{z}$. Fits in the middle row suggest that for cMBL at $W_{z}=0$, the
entanglement (QFI) growth follows a double (triple) logarithmic form. The
bottom row targets the dynamics of dMBL at large $W_{z}$: consistent with the
logarithmic rise of $[f_{Q}]$, $[S_{\textrm{vN}}]$ grows as a power law of $t$
in dMBL. The four insets of (c)-(f) present the corresponding TEBD results of
$L=28$. Figure 6:  (a1),(a2): A replot of $[S_{\textrm{vN}}]$ in cMBL
(stabilized by $\frac{g_{x}}{W_{x}}=1,\frac{W_{z}}{W_{x}}=0$) on an open chain
of $L=18$ but now fitted by two different types of functions: $\log\log(t)$
and $\log(t)$. The lower panel is the same as the upper one but in a semi-log
format. Apparently, $\log\log(t)$ gives the better fit. (b1),(b2): The
$[S_{\textrm{vN}}]$ data of a longer open chain $L=22$ with the same model
parameters. In each panel, the time interval that matches the double-log fit
is marked by two dashed lines. Noticeably, this time window broadens from
$t\in(100,10^{7.5})$ for $L=18$ to $t\in(360,10^{9})$ for $L=22$. (c),(d):
Temporal growth of $[f_{Q}]$ at two representative points inside the cMBL dome
using longer chains. Both can be captured by the triple-log functions of time.
(c) corresponds to $W_{z}=0.5$, while (d) repeats Fig. 5(d) for comparison.

Experimentally, the QFI, which sets the lower bound of entanglement, was
measured in trapped-ion chain [71] to witness entanglement growth under the
interplay between MBL and long-range interactions. Following [71], we start
from Néel states in even chains,
$|\psi(t\\!=\\!0)\rangle\\!=\\!|\\!\\!\downarrow\uparrow\\!\ldots\\!\downarrow\uparrow\rangle$,
characterized by a staggered $\mathbb{Z}_{2}$ spin-imbalance operator,
${{I}}\coloneqq\frac{1}{L}\sum^{L}_{i=1}(-1)^{i}\sigma^{z}_{i}$, then the
associated QFI density reduces to the connected correlation function of
${{I}}$,
$f_{Q}(t)=4L(\langle\psi(t)|{{I}}^{2}|\psi(t)\rangle-\langle\psi(t)|{{I}}|\psi(t)\rangle^{2})$,
which links multipartite entanglement to the fluctuations encoded in
measurable quantum correlators. Figure 5(b) is a semi-log plot of averaged
$[f_{Q}]$ along the line $g_{x}=1$ with different $W_{z}$ color-coded the same
way as in Fig. 5(a). Likewise, the notable change in functional form of
$[f_{Q}]$ echoes the same cMBL-dMBL eigenstate transition. Specifically, Fig.
5(d) shows that the long-time growth of $[f_{Q}]$ in cMBL matches a triple-log
form, which reinforces the double-log fit of $[S_{\textrm{vN}}]$ in (c).
Parallel relation between $[S_{\textrm{vN}}]$ and $[f_{Q}]$ carries over to
the dMBL phase where the power-law growth of $[S_{\textrm{vN}}]$ in (e)
transforms into a logarithmic growth of $[f_{Q}]$ in (f).

Table 1 recaps the cMBL-dMBL distinction in the fundamental dynamical aspects
of entanglement and its witness.

To supplement the ED simulation in main panels, we employ TEBD and matrix-
product-operator techniques to verify the cMBL-dMBL transition in larger
system sizes $L=28$. A $4$th-order Suzuki-Trotter decomposition is implemented
at a maximal time step of unity (the appropriate time step is set by the
inverse mean gap and for large $W_{z}$, we find that smaller time step of
about $0.1$ is needed, but as the entanglement growth deep inside dMBL is very
slow, TEBD maintains its efficiency at this small time step), and the
truncation error per step is controlled lower than $10^{-6}$ by selecting a
large bond-dimension threshold $1500$. We check that for $t\lessapprox 10^{3}$
the accumulated total truncation error in typical quasirandom samples is well
below $10^{-5}$. The corresponding results and the fits are consistently
presented in the insets of Fig. 5. However, due to entanglement accumulation,
matrix-product-state algorithms of this type retain effectiveness within
limited time scales $(t\lessapprox 10^{3})$.

Table 1: Hierarchies of dynamic characteristics encompassing constrained, unconstrained, and diagonal MBL phases. | $[S_{\textrm{vN}}]$ |  $[$Quantum Fisher Info.$]$
---|---|---
cMBL | $\log\log\left(t\right)$ | $\log\log\log\left(t\right)$
uMBL | $\log\left(t\right)$ | $\log\log\left(t\right)$
dMBL | $t^{\alpha}$ | $\log\left(t\right)$

Double-log entanglement growth.—As refinement, Figs. 6(a1),(a2) manifest that
for cMBL, the double-log fitting function $\log\log(t)$ matches the
$[S_{\textrm{vN}}]$ data significantly better than the single-log fitting
function $\log(t)$, the hallmark of the unconstrained MBL [72, 73, 74, 75].
Likewise, Figs. 6(b1),(b2) show the entanglement data for larger chain’s
length $L=22$ with the model parameters,
$\frac{g_{x}}{W_{x}}=1,\frac{W_{z}}{W_{x}}=0$, intact. As indicated by pairs
of dashed lines there, it is estimated that the duration that best traces the
double-log fit increases progressively from $t\in(100,10^{7.5})\ [\Delta
t\approx 3.162\times 10^{7}]$ in system $L=18$ to $t\in(360,10^{9})\ [\Delta
t\approx 10^{9}]$ in system $L=22$, thus supporting the double-log
entanglement buildup in cMBL.

Triple-log rise of QFI.—Complementarily, Figs. 6(c),(d) reproduce the temporal
evolutions of $[f_{Q}]$ for two chosen parameters $W_{z}=0,0.5$ on longer
chains. The fits on top of data from various lengths corroborate the
speculation that this triple-log growth of $[f_{Q}]$ may comprise a feature
shared within the dome of cMBL. Similar to Figs. 6(a),(b), here we resort to
ED and raise the maximum open-chain size to $L=22$.

Eigenstate transition from transport

Additionally, there are marked differences between cMBL and dMBL, as reflected
through the chain’s relaxation from the prepared Néel state and the spread of
initialized local energy inhomogeneity. In accordance with the time evolution
of $[S_{\textrm{vN}}]$ and $[f_{Q}]$, the decay of
${{I}}(t)\coloneqq\langle\psi(t)|{{I}}|\psi(t)\rangle$ is examined in Fig.
7(a). Apart from a quick suppression during $t\lessapprox 1$, both MBLs relax
to a steady state with finite magnetization. They thus retain remnants of the
initial spin configuration in contrast to the thermal phase where $[{{I}}(t)]$
vanishes irrevocably. Notice that under the increase of $W_{z}$, the frozen
moment $[{{I}}_{\infty}]$ at infinite $t$ develops monotonously from $\sim
0.5$ in cMBL up to $\sim 0.9$ in dMBL; before equilibration, the intermediate
oscillation of $[{{I}}(t)]$ is also damped more severely in dMBL than in (off-
diagonal) cMBL.

Following [76], the energy transport of the constrained model is scrutinized
by monitoring the spread of a local energy inhomogeneity initialized on the
central site of an odd chain at infinite temperature, i.e., the system’s
initial density matrix assumes
$\rho(t\\!=\\!0)=\frac{1}{\dim\\!\mathcal{H}}(\mathbb{1}+\varepsilon\widetilde{X}_{\frac{L+1}{2}})$,
where $\dim\\!\mathcal{H}$ the dimension of projected Hilbert space and
$\varepsilon$ the disturbance of energy on site $i_{c}\coloneqq(L+1)/2$. The
quantity measuring the effective distance $\varepsilon$ travels is
${{R}}(t)\coloneqq\frac{1}{{\sf
Tr}\left[\widetilde{\rho}(t)H_{\mathrm{qp}}\right]}\sum^{L}_{i=1}\left\\{\left|i-i_{c}\right|{\sf
Tr}\left[\widetilde{\rho}(t)H_{i}\right]\right\\}$, where $H_{i}\coloneqq
g_{i}\widetilde{X}_{i}+h_{i}\widetilde{Z}_{i}$ and the time-independent
background is subtracted via inserting
$\widetilde{\rho}(t\\!=\\!0)\\!\coloneqq\\!\frac{1}{\dim\\!\mathcal{H}}\varepsilon\widetilde{X}_{\frac{L+1}{2}}$.
As per ETH, $\varepsilon$ is eventually smeared uniformly over the chain by
unitary time evolution and in that circumstance
$[{{R}}(t\\!=\\!\infty)]\approx L/4$. Figure 7(b) contrasts the behaviour of
$[{{R}}(t)]$ between cMBL and dMBL. Concretely, for dMBL, $[{{R}}]$ stays
vanishingly small, thereby $\varepsilon$ remains confined to $i_{c}$ and shows
no diffusion toward infinite $t$. In comparison, as the consequence of a fast
expansion within $t\lessapprox 100$, largely due to contributions from nearest
and next-nearest neighbours, $\varepsilon$ spreads over a finite range of the
chain in cMBL. Here, however, the saturated value $[{{R}}_{\infty}]$ after an
oscillatory relaxation remains subthermal.

Integrals of motion and dynamical order parameters

Key distinction between cMBL and dMBL can be further resolved from studying
the long-time limit of the spatial distribution of the energy-inhomogeneity
propagation. We utilize three quantities to access this information
complementarily.

(i) For each quasirandom realization, we parse the definition of $R(t)$ as per
the site index, $\varepsilon_{i}(t)\coloneqq\frac{{\sf
Tr}\left[\widetilde{\rho}(t)H_{i}\right]}{{\sf
Tr}\left[\widetilde{\rho}(t)H_{\mathrm{qp}}\right]}$, which measures in
percentage the extra energy on position $i$ with respect to the total
conserved perturbation $\varepsilon$. Observing that $\varepsilon_{i}$
approaches a constant $\varepsilon_{i,\infty}$ at infinite $t$, one might
implement the trick [77],
$\lim\limits_{T\rightarrow\infty}\frac{1}{T}\int^{T}_{0}\\!O(t)dt\approx\sum\limits_{n}\langle
n|O|n\rangle|n\rangle\langle n|$, to extract its value with the aid of
randomness,

$\varepsilon_{i,\infty}\coloneqq\varepsilon_{i}(t\\!\rightarrow\\!\infty)\approx\frac{\sum\limits_{n}\langle
n|\widetilde{X}_{\frac{L+1}{2}}|n\rangle\langle
n|H_{i}|n\rangle}{\sum\limits_{n}E_{n}\langle
n|\widetilde{X}_{\frac{L+1}{2}}|n\rangle},$ (13)

where $\\{|n\rangle\\}$ comprises an eigenbasis satisfying
$H_{\mathrm{qp}}|n\rangle=E_{n}|n\rangle$. Evidently, the profile of
$\\{\varepsilon_{i,\infty}\\}$ bears information pertaining to the local
structure of integrals of motion (IOMs).

(ii) The summation of $\varepsilon_{i,\infty}$ weighted by the separation
returns the equilibrated value of the effective traveling distance,
$R_{\infty}=\sum_{i=1}^{L}(|i-i_{c}|\cdot\varepsilon_{i,\infty})$.

(iii) Viewing that the contribution from $i_{c}$, i.e., the return
probability, is missing from $R_{\infty}$, one can define
$\varepsilon_{i_{c},\infty}$ as the residual energy density on the release
place,
$\varepsilon_{\mathrm{res}}\coloneqq\varepsilon_{\frac{L+1}{2},\infty}$.

All three quantities defined above are used to distinguish ETH and MBL. Here
we show that they are also the dynamical “order parameters” to differentiate
between the cMBL and dMBL regimes and identify the transition point therein.

Figure 7:  cMBL-dMBL transition in transport with PBCs and fixed
$g_{x}=1,L=21$. (a),(b): Time evolution of the $\mathbb{Z}_{2}$
antiferromagnetic imbalance $[{{I}}(t)]$ and the energy spread $[{{R}}(t)]$ as
a function of $W_{z}$. (c) exemplifies the peak-dip-hump lineshape of
$[\varepsilon_{i,\infty}]$ in cMBL for two irrational wavenumbers $\xi$. One
set of time-profiles of $[\varepsilon_{i_{c}\pm 2,3}(t)]$ that characterizes
the nonmonotonicity of the dip-hump structure is given by (d). (e) shows the
lineshape of $[\varepsilon_{i,\infty}]$ in dMBL; the exponential decay can be
seen from the semi-log inset wherein the cMBL data (pink dots) are overlaid
for comparison. (f): The changes in dynamic “order parameters” $[R_{\infty}]$
and $[\varepsilon_{\mathrm{res}}]$ under the tuning of $W_{z}$ signal the
transition between cMBL and dMBL. Light to solid colours in (d),(f) correspond
to $L=17,19,21$.

LIOMs and positive definiteness of dMBL

Despite the central status of LIOMs in disorder-induced MBL [17, 18, 19, 20],
LIOMs in unconstrained aperiodic MBL systems receive attention only recently
[47, 51]. Ref. [47] constructed LIOMs of MBL as time-averaged local operators
for interacting fermions subject to aperiodic potentials. They found that in
this circumstance $\ell$-bits remain localized even at the vicinity of the
quasiperiodic MBL transition. Likewise, upon continuous unitary transforms,
Ref. [51] computed the real-space support of LIOM in quasirandomness-induced
MBL and revealed that the effective interactions between LIOMs exhibit
features inherited from the underlying aperiodic potential. Interestingly,
both works pointed to the weaker finite-size effects in aperiodic modulations
than in truly disordered arrangements. Exploiting the instability of LIOMs,
they also found that the associated MBL transition may occur at a higher
critical quasirandom strength than previously estimated [40].

Before proceeding to numerics, let’s gain some understanding on dMBL within
the LIOM framework. The first step forward is to introduce
$\breve{Z}_{i}\coloneqq\mathcal{P}_{i+1}\widetilde{Z}_{i}\mathcal{P}_{i-1}$
where $\mathcal{P}_{i}\coloneqq\frac{1}{2}(\mathbb{1}+\sigma^{z}_{i})$ as the
building blocks of constrained $\ell$-bits. The convenience of $\breve{Z}_{i}$
stems from ${\sf Tr}\breve{Z}_{i}=0$, which contrasts to ${\sf
Tr}\widetilde{Z}_{i}>0$, thereby $\breve{Z}_{i}$ behaves like a spin free of
restrictions. Following [23], it can be proved that as long as $W_{z}\gg
g_{x}+W_{x}$, the set of tensor-product operators
$\mathcal{I}_{L}\coloneqq\\{\mathcal{Z}_{i_{1}}\otimes\cdots\otimes\mathcal{Z}_{i_{k}}\\}$
fulfilling $1\leqslant i_{1}\leqslant i_{2}\leqslant\cdots i_{k}\leqslant L,\
i_{a+1}\neq i_{a},\ 1\leqslant k\leqslant\frac{L+1}{2}$ may be constructed as
a complete, mutually commuting, and linearly-independent basis to express any
nontrivial operators that commute with $H_{\textrm{qp}}$; in terms of
quasilocal unitaries, $\mathcal{Z}_{i_{a}}\approx
U\breve{Z}_{i_{a}}U^{\dagger}$. This is because the set of states
$\\{|\mathcal{Z}_{i_{1}}\mathcal{Z}_{i_{2}}\cdots\mathcal{Z}_{i_{k}}\rangle\\}$
derived from $\mathcal{I}_{L}$ reproduces faithfully the effective eigenbasis
of projected Hilbert space for dMBL. Accordingly, the IOM in Eq. (13) is
recast into

$\mbox{dMBL:}\ \ \sum\limits_{n}\langle
n|\widetilde{X}_{i}|n\rangle|n\rangle\langle
n|\approx\sum^{\frac{L-1}{2}}_{m=0}\sum_{r}V^{[i]}_{r,m}\widehat{\mathcal{O}}^{[i]}_{r,m},$
(14)

where $\widehat{\mathcal{O}}^{[i]}_{r,m}$ denotes the element of
$\mathcal{I}_{L}$ possessing the support on site $i$ (i.e., contains
$\mathcal{Z}_{i}$) and whose furthest boundary from $i$ is of distance $m$.
The nonidentical members comprising this specified subset are labelled by $r$.
Besides the finite support of $\mathcal{Z}_{i}$, the key property that
promotes $\sum_{n}\langle n|\widetilde{X}_{i}|n\rangle|n\rangle\langle n|$ to
the LIOM of dMBL is the locality condition of its real coefficients, i.e.,
$V^{[i]}_{r,m}\sim e^{-m/\zeta}$. In terms of LIOM representation, the
universal Hamiltonian governing the dynamics of dMBL may assume the following
form,
$H^{\textrm{dMBL}}_{\textrm{qp}}=\sum_{i}\widetilde{h}_{i}\mathcal{Z}_{i}+\sum_{k}\sum_{i_{1}\ldots
i_{k}}J_{i_{1}\ldots
i_{k}}\mathcal{Z}_{i_{1}}\mathcal{Z}_{i_{2}}\cdots\mathcal{Z}_{i_{k}}$, where
from Figs. 5(e),(f), it is feasible to infer $J_{i_{1}\ldots
i_{k}}\sim|i_{k}-i_{1}|^{-1/\alpha}\cdot\phi^{-|i_{k}-i_{1}|}$, which decays
as an exponentially-suppressed power law of LIOMs’ separation. Here, $\alpha$
is the same exponent in Table 1 and $\phi$ is the golden ratio.

Being the trace of product of two IOMs, one consequence of Eq. (14) is the
positive definiteness of the averaged $[\varepsilon_{i,\infty}]$ featured by
an exponential decay in space. Figure 7(e) illustrates that this is the case
even when $W_{z}\approx g_{x}+W_{x}$.

Peak, dip, hump in cMBL

Now we are in the position to highlight the peak-dip-hump structure and the
occurrence of negativity in $[\varepsilon_{i,\infty}]$ [see Figs. 7(c),(d)] as
the peculiarities of cMBL that distinguish it from both dMBL and unconstrained
MBL (uMBL) by the presence of pronounced nonlocal correlations. The
unambiguous negativity of $[\varepsilon_{i_{c}\pm 2}]$ in Fig. 7(c) and the
nonmonotonicity of $[\varepsilon_{i_{c}\pm 2,3}]$ in Fig. 7(d) point to the
insufficiency of Eq. (14) when addressing the cMBL phase from the dMBL side.
Especially, they highlight the dynamical consequence that in cMBL the
correlation between the centre site (where the initial energy inhomogeneity
locates) and its third nearest neighbours might be stronger than that for its
second nearest neighbours, because phenomenologically the net energy current
flowing into the second nearest neighbouring sites could appear noticeably
less than that flows out. To remedy the inconsistency, we propose as a
scenario that the missing pieces may come from the terms in $\mathcal{I}_{L}$
that are nonlocal with respect to $i$, viz., their support on $i$ vanishes,
hence, for cMBL, $\sum_{n}\langle n|\widetilde{X}_{i}|n\rangle|n\rangle\langle
n|\approx\sum^{\frac{L-1}{2}}_{m=0}\sum_{r,\overline{r}}(V^{[i]}_{r,m}\widehat{\mathcal{O}}^{[i]}_{r,m}+V^{\overline{[i]}}_{\overline{r},m}\widehat{\mathcal{O}}^{\overline{[i]}}_{\overline{r},m})$.
The superscript ${\overline{[i]}}$ signifies the absence of $\mathcal{Z}_{i}$
in the associated expansion. Under the successive decrease of $W_{z}$, it is
anticipated that the weights $V^{\overline{[i]}}_{\overline{r},m}$ of small
$m$ grow significantly such that a finite-size core centred at $i$ forms
wherein nonlocal correlated contributions, albeit confined, become
predominant. On the contrary, for those $m$ beyond the core, the importance of
$V^{\overline{[i]}}_{\overline{r},m}$ diminishes sharply so that the rapid
decay tail and the overall signatures of localization are well maintained.

Alternatively, the core formation may be monitored by $[R_{\infty}]$ and
$[\varepsilon_{\textrm{res}}]$. Figure 7(f) illustrates that the duo
constitutes the desired “order parameters” from quantum dynamics that take
values zero and unity in dMBL and saturate to nontrivial plateaus in cMBL. The
critical $W_{z}$ of the transition is hence estimated to be $\sim 0.55$ at
$g_{x}=1$. Furthermore, from Fig. 7(c), the core where substantial nonlocal
effects take place spans roughly $5$ to $7$ lattice sites which, as per the
saturated value of $[R_{\infty}]$ in Fig. 7(f), is comparable to a thermal
segment of approximately $3$ lattice-spacing long.

This embedded thermal-like core in IOMs plays a crucial role in yielding the
novel Lieb-Robinson bound for cMBL. More relevant mathematical justifications
are in Ref. [58].

cMBL-dMBL transition and return probability

Figure 8:  Finite-size scaling analysis for cMBL-dMBL transition at fixed
$g_{x}=1$. The respective data collapse of random-averaged dynamic order
parameters $[\varepsilon_{\textrm{res}}]$ [panel (a)] and $[R_{\infty}]$
[panel (b)] yields consistently a critical $W^{c}_{z}\sim 0.55W_{x}$. Here,
PBCs are used; for $L=23,25$, the results are extracted from the Krylov-
typicality approach; other sizes are solved by ED. Figure 9:  (a): Time
evolution of return probability as a function of $W_{z}$ at fixed $g_{x}=1$
obtained by applying Krylov-typicality method to periodic chains. The light to
solid colours give the results of $L=21,23,25$. (b): A sample calculation on
small chain illustrates that Krylov-typicality approximation reproduces the
results of ED.

In view of the importance of core formation in driving the change between cMBL
and dMBL, we perform a finite-size scaling analysis on “dynamic order
parameters” $[\varepsilon_{\textrm{res}}]$ and $[R_{\infty}]$ to locate more
precisely the critical value of $W_{z}$ that triggers this eigenstate
transition. Following Ref. [68], the target quantity
$Q=[\varepsilon_{\textrm{res}}]$ (or $[R_{\infty}]$) at fixed $g_{x}$ might
assume a standard scaling form as follows,
$\frac{Q(L,W_{z})}{L^{\sigma_{Q}}}\approx
f_{Q}\\!\left((W_{z}-W^{c}_{z,Q})L^{\alpha_{Q}}\right)$, where $f_{Q}$ is some
unknown function associated to $Q$ and similarly $\sigma_{Q},\alpha_{Q}$ are
the corresponding scaling exponents. The critical strength $W^{c}_{z,Q}$ of
the transition can then be extracted from the proper data collapse of $Q$.

Figure 8 depicts scaling profiles of $[\varepsilon_{\textrm{res}}]$ and
$[R_{\infty}]$ at fixed $g_{x}=1$ for different system sizes ranging from
$L=11$ to $25$ using PBCs. As shown by panel (a), the scaling collapse of the
data for $[\varepsilon_{\textrm{res}}]$ yields the following set of
parameters, $W^{c}_{z,\varepsilon}\approx 0.55,\
\sigma_{\varepsilon}\approx-0.01,\ \alpha_{\varepsilon}\approx 0.2$.

Although $[\varepsilon_{\textrm{res}}]$ and $[R_{\infty}]$ carry compatible
information, as $[\varepsilon_{i_{c}}(t)]$ stands for return probability, an
autocorrelation function easier to measure in experiments, we extend its
evaluation to longer chains $L=23,25$ by resort to Krylov-typicality technique
(see methods section). Figure 9(a) presents the obtained
$[\varepsilon_{\frac{L+1}{2}}(t)]$ for a range of $W_{z}$ close to the
transition of $g_{x}=1$. It is noticeable that due to localization, all
evolution curves relax to their constant lineshapes whose saturation values
strengthen with $W_{z}$ and form two individual plateaus around $\sim 0.45$
for cMBL and $\sim 1$ for dMBL, respectively. This convergent trend allows for
an estimate of $[\varepsilon_{\textrm{res}}]$ by averaging the return
probability over a later period $t\in[1000,2000]$, which produces the data
points of $L=23,25$ in Fig. 8(a). The reliability of such a procedure is
justified by Fig. 9(b) where we check the correctness of Krylov-typicality
approximation and the appropriateness of the chosen time window via a
benchmark test against the exact results of $L=21$.

Because energy diffusion is ceased in localized regions (in other words, the
initial energy imbalance is restricted to the chain centre), adding more sites
on the chain ends generates small finite-size flows of
$[\varepsilon_{\textrm{res}}]$ seen in Fig. 8(a) even after raising the length
limit to $L=25$. Likewise, as plotted by Fig. 8(b), the estimates of scaling
parameters obtained from the data collapse of $[R_{\infty}]$ read
$W^{c}_{z,R}\approx 0.55,\ \sigma_{R}\approx 0.06,\ \alpha_{R}\approx 0.1$.

DISCUSSION

To conclude, we find a cMBL regime and a probably discontinuous cMBL-cETH
transition in the quasirandom Rydberg-blockade chain. The orthogonality
between the field strength and the projection direction renders cMBL and its
discontinuous eigenstate transition fundamentally different from dMBL, uMBL,
and the continuous uMBL-uETH transition. Particularly, the entanglement
entropy in cMBL grows as a double-log function of time, as opposed to the
power-law growth in dMBL and the single-log growth in uMBL.

The presumed discontinuity of the cMBL-cETH transition is evidenced
numerically. A future analytical elucidation of its underpinnings may
potentially improve the existing theoretical framework for MBL transition in a
substantial way.

Even though LIOMs capture the phenomenology of dMBL, the cMBL-dMBL transition
triggered by the rotation of the field orientation accentuates the importance
of nonlocal components in the IOMs of cMBL, which, together with the double-
log entanglement growth, raises doubts about how to define the meaningful
LIOMs and the universal Hamiltonian suitable for cMBL.

The continual investigations on these open questions promise to further our
understanding of unconventional MBL beyond the current scope.

METHODS

Quantum dynamics computations

In this work, to cope with the many-body nonequilibrium problem subject to
intertwining complexities from constraint and randomness, three numerical
approaches, ED, Krylov-typicality, and TEBD, are employed.

ED.—For small chains, we resort to the ED method to access the long-time
limit, where quadruple precision is implemented for achieving the time scale
up to $t\approx 10^{29}$. Within full diagonalization, the infinite-time limit
is resolvable by invoking the diagonal approximation. Further, rather than
removing the $\downarrow\downarrow$-motifs from unconstrained Hilbert space,
we construct the projected spin-$\frac{1}{2}$ basis as a selected set of
binary numbers by fulfilling the constraint rule using combinatorial
reasoning, which is more efficient for larger system sizes.

TEBD.—One alternative to evaluate the time evolution of longer quantum spin
chains, albeit with the limitation of much shorter time scales, is the TEBD
algorithm [70], which is built upon parametrization of a quantum wavefunction
in terms of matrix-product states (MPS) [78],

$|\psi\rangle=\sum_{\sigma_{1},\ldots,\sigma_{L}}A^{[1]}A^{[2]}\cdots
A^{[L]}|\sigma_{1},\ldots,\sigma_{L}\rangle,$ (15)

where $A^{[i]}$ stands for a three-leg tensor at site $i$ carrying one
physical bond $\sigma_{i}=1,2$ for a local spin-$\frac{1}{2}$ system and two
virtual legs of dimension $\chi^{i}_{L}$ and $\chi^{i}_{R}$. TEBD relies on
the low amount of entanglement generation and the Suzuki-Trotter decomposition
of the time evolution operator. Concretely, at $4$th order, this unitary can
be approximated in a symmetric format [78],

$e^{-iH\tau}=U(\tau_{1})U(\tau_{2})U(\tau_{3})U(\tau_{2})U(\tau_{1})+\mathcal{O}(\tau^{5}),$
(16)

where

$\displaystyle U(\tau_{i})$
$\displaystyle=e^{-iH_{\textrm{odd}}\tau_{i}/2}e^{-iH_{\textrm{even}}\tau_{i}}e^{-iH_{\textrm{odd}}\tau_{i}/2},$
(17) $\displaystyle\tau_{1}=\tau_{2}$
$\displaystyle=\frac{\tau}{4-\sqrt[3]{4}},\ \ \ \ \
\tau_{3}=\tau-2(\tau_{1}+\tau_{2}),$ (18)

and we assume that the total inspected Hamiltonian $H$ comprises a sum of two-
site operators that can be divided into the respective $H_{\textrm{even}}$ and
$H_{\textrm{odd}}$ parts living across the even and odd bonds. Evidently,
starting from an arbitrary product state in the projective spin basis, to a
good approximation, the repeated application of the unitary time evolution
will not generate components that violate the constraint.

The calculation of QFI entails the evaluation of $I^{2}$, which is easily
computed within TEBD via recasting $I$ as a matrix-product operator (MPO)
[78]: analogous to the MPS representation in Eq. (15), a generic operator $O$
is rewritten as

$O=\sum_{\begin{subarray}{c}\sigma_{1},\ldots,\sigma_{L}\\\
\sigma^{\prime}_{1},\ldots,\sigma^{\prime}_{L}\end{subarray}}W^{[1]}W^{[2]}\cdots
W^{[L]}|\sigma_{1},\ldots,\sigma_{L}\rangle\langle\sigma^{\prime}_{1},\ldots,\sigma^{\prime}_{L}|,$
(19)

where $W^{[i]}$ is a four-leg tensor on site $i$ equipped with two physical
bonds $\sigma_{i},\sigma^{\prime}_{i}$ and two virtual bonds of dimension
$D\times D$. For spin-imbalance operator $I$, $W$’s are simply given by

$W^{[i]}=\begin{bmatrix}\mathds{1}&f_{i}\sigma^{z}_{i}\\\\[5.0pt]
0&\mathds{1}\end{bmatrix},$ (20)

where $f_{i}=(-1)^{i}/L$, and $\mathds{1},\sigma^{z}_{i}$ are $2\times 2$
Pauli matrices, therefore $D=2$ in this case. Then, $I^{2}$ consists of a
stacking of two identical layers of $W$’s tensors, whose expectation value at
$t$ is obtained by executing optimal contractions.

Krylov-typicality.—For brevity, let us recap the main steps of Krylov-space
technique. By definition, starting from an arbitrary normalized wave-vector
$|\phi_{0}\rangle$, the associated linearly-independent Krylov subspace
$\mathcal{K}$ is generated by consecutively applying the Hamiltonian $H$ onto
$|\phi_{0}\rangle$ $(m-1)$ times,

$\mathcal{K}\coloneqq\left\\{|\phi_{0}\rangle,H|\phi_{0}\rangle,H^{2}|\phi_{0}\rangle,\cdots,H^{m-1}|\phi_{0}\rangle\right\\}.$
(21)

Through clever recombination, an equivalent but more convenient reformulation
of $\mathcal{K}$ exists, which is mutually orthonormal and called the Lanczos
basis derived from $|\phi_{0}\rangle=|v_{0}\rangle$,

$\mathcal{K}\sim\mathcal{L}\coloneqq\left\\{|v_{0}\rangle,|v_{1}\rangle,|v_{2}\rangle,\cdots,|v_{m-1}\rangle\right\\},$
(22)

where, to remedy the loss of orthogonality, the procedure of reorthogonality
is always assumed. The advantage of $\mathcal{L}$ lies in the fact that for
most practical calculations, it suffices to choose the Lanczos dimension
$m\approx 50$ to $100$, which is orders of magnitude smaller than the full
Hilbert-space dimension $\mathcal{D}$.

The above rationale can be recapitulated in terms of the following basis
transformation,

$F^{\dagger}HF=H_{\textrm{lanc}},$ (23)

where the full Hamiltonian $H$ is written in the original physical basis,
while the heavily reduced Hamiltonian $H_{\textrm{lanc}}$ is recast in the
Lanczos basis specific for the neighbourhood of
$|v_{0}\rangle=|\phi_{0}\rangle$. It is easy to prove that for Hermitian
operator $H$, $H_{\textrm{lanc}}$ is a tridiagonal matrix. Moreover, stacking
the Lanczos states yields the fundamental transformation matrices
$F,F^{\dagger}$,

$\displaystyle F=\left[|v_{0}\rangle\ |v_{1}\rangle\ |v_{2}\rangle\ \cdots\
|v_{m-1}\rangle\right]_{\mathcal{D}\times m},$ (24) $\displaystyle
F^{\dagger}=\begin{bmatrix}\langle v_{0}|\\\ \langle v_{1}|\\\ \vdots\\\
\langle v_{m-1}|\end{bmatrix}_{m\times\mathcal{D}}.$ (25)

The essence of Lanczos approximation can then be encapsulated in terms of the
following single relation,

$FF^{\dagger}\approx\mathds{1}_{m\times m},$ (26)

which is exact iff $m$ equals the Hilbert-space dimension $\mathcal{D}$.

Armed with these preparations, we are ready to derive the formula for the
real-time propagation of a normalized vector
$|\psi(t)\rangle=|\phi_{0}\rangle=|v_{0}\rangle$ under the unitary evolution
of the Hamiltonian $H$ up to a small time decimation $\delta$, i.e.,

$\displaystyle|\psi(t+\delta)\rangle\approx
FV_{\textrm{lanc}}e^{-i\frac{\delta}{\hbar}D_{\textrm{lanc}}}V^{\dagger}_{\textrm{lanc}}F^{\dagger}|\psi(t)\rangle,$
(27)

where the tridiagonal Lanczos matrix $H_{\textrm{lanc}}$ is diagonalized by
the unitary matrix $V_{\textrm{lanc}}$,
$H_{\textrm{lanc}}=V_{\textrm{lanc}}D_{\textrm{lanc}}V^{\dagger}_{\textrm{lanc}}$.
Symbolically, one writes the combined vector $F^{\dagger}|\psi(t)\rangle$ in
the explicit form,

$F^{\dagger}|\psi(t)\rangle=\begin{bmatrix}\frac{\langle\psi(t)|\psi(t)\rangle}{\sqrt{\langle\psi(t)|\psi(t)\rangle}}\\\
0\\\ \vdots\\\ 0\end{bmatrix}_{m\times 1},$ (28)

where the zeros are resultant from the orthogonality between different Lanczos
basis states.

Notice that both $|\psi(t)\rangle$ and $|\psi(t+\delta)\rangle$ are
$\mathcal{D}\times 1$ vectors in the computational basis, therefore the
evaluation of the half-chain von Neumann entanglement entropy at $t+\delta$ is
proceeded in the usual way once $|\psi(t+\delta)\rangle$ is available.

For an initial energy inhomogeneity at the central site of a Rydberg chain,
its return probability under the unitary time evolution of the quasiperiodic
Hamiltonian $H_{\textrm{qp}}$ is

$\varepsilon_{\frac{L+1}{2}}(t)\coloneqq\frac{1}{{\sf
Tr}\left[\widetilde{\rho}(t)H_{\textrm{qp}}\right]}{\sf
Tr}\left[\widetilde{\rho}(t)H_{\frac{L+1}{2}}\right].$ (29)

Here we assume the infinite temperature. The initial density matrix
$\widetilde{\rho}(t\\!=\\!0)$ specifies the spatial distribution of the energy
disturbance at $t=0$, whose subsequent dynamics can be couched in the
Heisenberg representation as $(\hbar=1)$,

$\displaystyle\widetilde{\rho}(t)\coloneqq
e^{iH_{\textrm{qp}}t}\cdot\widetilde{\rho}(t=0)\cdot e^{-iH_{\textrm{qp}}t},$
(30)
$\displaystyle\widetilde{\rho}(t=0)\coloneqq\frac{\varepsilon}{\mathcal{D}}\widetilde{X}_{\frac{L+1}{2}}.$
(31)

Note that $H_{\textrm{qp}}$ is discretized in Eq. (29), i.e.,

$\displaystyle H_{\textrm{qp}}=\sum^{L}_{i=1}H_{i},\ \ \ \ \ \mbox{where}\ \ \
\ \ H_{i}=g_{i}\widetilde{X}_{i}+h_{i}\widetilde{Z}_{i}.$ (32)

There exists one extra complication in the computation of
$\varepsilon_{\frac{L+1}{2}}(t)$. By full diagonalization, the trace over the
entire Hilbert space is accomplishable through summing up the contributions of
all eigenvectors with the equal weight,

${\sf Tr}[\ldots]=\sum^{\mathcal{D}}_{n=1}\langle E_{n}|\ldots|E_{n}\rangle.$
(33)

This approach is impractical once $L\geqslant 23$. For instance, when $L=23$
under PBCs, $\mathcal{D}=64079$ and the required ram to compute a single
random realization is over $190$GB.

To make progress, we invoke the trick of dynamical typicality [79, 80, 81] to
approximately evaluate the trace of an operator over the gigantic Hilbert
space. The key idea is to replace Eq. (33) by a single scalar product using a
pure state $|\Psi_{\textrm{gaus}}\rangle$,

$\frac{1}{\mathcal{D}}{\sf
Tr}[\ldots]=\frac{1}{\mathcal{D}}\sum^{\mathcal{D}}_{n=1}\langle
E_{n}|\ldots|E_{n}\rangle\approx\langle\Psi_{\textrm{gaus}}|\ldots|\Psi_{\textrm{gaus}}\rangle,$
(34)

where $|\Psi_{\textrm{gaus}}\rangle$ is written in the computational basis
with dimension $\mathcal{D}$ whose entries are Gaussian random numbers with
zero means. The Gaussian distribution of the complex components of
$|\Psi_{\textrm{gaus}}\rangle$ guarantees that $|\Psi_{\textrm{gaus}}\rangle$
is drawn uniformly on the hypersphere of the full Hilbert space (i.e.,
according to the Haar measure) such that the corresponding probability
distribution is invariant under all unitary transformations within the Hilbert
space. As per formal theory of typicality [79, 80],
$|\Psi_{\textrm{gaus}}\rangle$ could be an effective representative of the
underlying statistical ensemble.

Practically, $|\Psi_{\textrm{gaus}}\rangle$ is constructed in a simple manner,

$|\Psi_{\textrm{gaus}}\rangle=\frac{1}{\mathcal{N}}\sum^{\mathcal{D}}_{a=1}\left(r_{a}+is_{a}\right)|a\rangle,$
(35)

where $|a\rangle$ enumerates the physical computational basis, and
$r_{a},s_{a}$ are real, independent Gaussian random numbers with mean zero and
variance one,
$\mathcal{N}=\sqrt{\sum^{\mathcal{D}}_{a=1}\left[r^{2}_{a}+s^{2}_{a}\right]}$.
As $\mathcal{D}$ is at the order of $10^{7}$ (for example,
$\mathcal{D}=1346269$ for $L=29$ under OBCs), for averages over $2000$
independent quasirandom realizations, a good random number generator with long
period $(\sim 10^{18})$ might be needed.

In real computation, the accuracy of dynamical typicality can be improved by
using multiple pure Gaussian states in a single evaluation of Eq. (34). To
maximize the overall landscape for the random realizations of
$H_{\textrm{qp}}$ and $|\Psi_{\textrm{gaus}}\rangle$, it is economic to invoke
one independent $|\Psi_{\textrm{gaus}}\rangle$ for each different
quasiperiodic $H_{\textrm{qp}}$, and then perform the average over $2000$ such
joint samples. Accordingly, under the assumption of the validity of dynamical
typicality, the random-averaged return probability might be approximated as
follows,

$[\varepsilon_{\frac{L+1}{2}}(t)]\approx\frac{1}{\mathcal{R}}\sum^{\mathcal{R}}_{q=1}\langle\Psi^{q}_{\textrm{gaus}}|\left[\widetilde{X}^{q}_{\frac{L+1}{2}}(t)\widetilde{X}_{\frac{L+1}{2}}\right]|\Psi^{q}_{\textrm{gaus}}\rangle,$
(36)

where $q$ denotes the involvement of the $q$th random sample and the total
number for the joint random samples is over $2000$, $\mathcal{R}\geqslant
2000$. Being a scalar product of pure state, each summand in Eq. (36) is
evaluable by the Lanczos method.

For all results presented in this paper, we perform the random sample
calculations over at least $1000$ independent quasiperiodic configurations of
the model parameters $\phi_{x},\phi_{z}$ and when evaluating entanglement
growth, the initial product states are additionally selected from the
constrained spin basis in another randomized manner [73]. The corresponding
statistical uncertainties are estimated from the normal variance of varied
averaged quantities as per [62].

DATA AVAILABILITY

The data set that supports the findings of the present study can be available
from the corresponding authors via email upon reasonable request.

ACKNOWLEDGEMENTS

The discussion with M. Heyl was acknowledged. This work is supported by the
SKP of China (Grant Nos. 2016YFA0300504 and 2017YFA0304204) and the NSFC Grant
No. 11625416.

AUTHOR CONTRIBUTIONS

All authors contributed equally to this work.

COMPETING INTERESTS

The authors declare no competing financial or non-financial interests.

## References

  * Anderson [1958] P. W. Anderson, Absence of diffusion in certain random lattices, Phys. Rev. 109, 1492 (1958).
  * Basko  _et al._ [2006] D. Basko, I. Aleiner, and B. Altshuler, Metal- insulator transition in a weakly interacting many-electron system with localized single-particle states, Annals of Physics 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).
  * Deutsch [1991] J. M. Deutsch, Quantum statistical mechanics in a closed system, Phys. Rev. A 43, 2046 (1991).
  * Srednicki [1994] M. Srednicki, Chaos and quantum thermalization, Phys. Rev. E 50, 888 (1994).
  * 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).
  * Nandkishore and Huse [2015] R. Nandkishore and D. A. Huse, Many-body localization and thermalization in quantum statistical mechanics, Annu. Rev. Condens. Matter Phys. 6, 15 (2015).
  * Turner  _et al._ [2018] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Weak ergodicity breaking from quantum many-body scars, Nat. Phys. 14, 745 (2018).
  * van Horssen  _et al._ [2015] M. van Horssen, E. Levi, and J. P. Garrahan, Dynamics of many-body localization in a translation-invariant quantum glass model, Phys. Rev. B 92, 100305 (2015).
  * [10] J. M. Hickey, S. Genway, and J. P. Garrahan, Signatures of many-body localisation in a system without disorder and the relation to a glass transition, J. Stat. Mech. (2016), 054047.
  * Pancotti  _et al._ [2020] N. Pancotti, G. Giudice, J. I. Cirac, J. P. Garrahan, and M. C. Bañuls, Quantum east model: Localization, nonthermal eigenstates, and slow dynamics, Phys. Rev. X 10, 021051 (2020).
  * Smith  _et al._ [2017] A. Smith, J. Knolle, D. L. Kovrizhin, and R. Moessner, Disorder-free localization, Phys. Rev. Lett. 118, 266601 (2017).
  * Brenes  _et al._ [2018] M. Brenes, M. Dalmonte, M. Heyl, and A. Scardicchio, Many-body localization dynamics from gauge invariance, Phys. Rev. Lett. 120, 030601 (2018).
  * Sierant  _et al._ [2017] P. Sierant, D. Delande, and J. Zakrzewski, Many-body localization due to random interactions, Phys. Rev. A 95, 021601 (2017).
  * Li  _et al._ [2017] X. Li, D.-L. Deng, Y.-L. Wu, and S. Das Sarma, Statistical bubble localization with random interactions, Phys. Rev. B 95, 020201 (2017).
  * Bar Lev  _et al._ [2016] Y. Bar Lev, D. R. Reichman, and Y. Sagi, Many-body localization in system with a completely delocalized single-particle spectrum, Phys. Rev. B 94, 201116 (2016).
  * Ros  _et al._ [2015] V. Ros, M. Müller, and A. Scardicchio, Integrals of motion in the many-body localized phase, Nuclear Physics B 891, 420 (2015).
  * Imbrie  _et al._ [2017] J. Z. Imbrie, V. Ros, and A. Scardicchio, Local integrals of motion in many-body localized systems, Annalen der Physik 529, 1600278 (2017).
  * 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).
  * Imbrie [2016] J. Z. Imbrie, Diagonalization and many-body localization for a disordered quantum spin chain, Phys. Rev. Lett. 117, 027201 (2016).
  * Bernien  _et al._ [2017] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Probing many-body dynamics on a $51$-atom quantum simulator, Nature (London) 551, 579 (2017).
  * Chen  _et al._ [2018] C. Chen, F. Burnell, and A. Chandran, How does a locally constrained quantum system localize?, Phys. Rev. Lett. 121, 085701 (2018).
  * Š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).
  * Š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).
  * Panda  _et al._ [2020] R. K. Panda, A. Scardicchio, M. Schulz, S. R. Taylor, and M. Žnidarič, Can we study the many-body localisation transition?, EPL (Europhysics Letters) 128, 67003 (2020).
  * Sierant  _et al._ [2020a] P. Sierant, D. Delande, and J. Zakrzewski, Thouless time analysis of anderson and many-body localization transitions, Phys. Rev. Lett. 124, 186601 (2020a).
  * Sierant  _et al._ [2020b] P. Sierant, M. Lewenstein, and J. Zakrzewski, Polynomially filtered exact diagonalization approach to many-body localization, Phys. Rev. Lett. 125, 156601 (2020b).
  * 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).
  * Sierant  _et al._ [2021] P. Sierant, E. G. Lazo, M. Dalmonte, A. Scardicchio, and J. Zakrzewski, Constraint-induced delocalization, Phys. Rev. Lett. 127, 126603 (2021).
  * Mondaini and Rigol [2015] R. Mondaini and M. Rigol, Many-body localization and thermalization in disordered hubbard chains, Phys. Rev. A 92, 041601 (2015).
  * 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).
  * 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).
  * 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).
  * Kiefer-Emmanouilidis  _et al._ [2021] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Slow delocalization of particles in many-body localized phases, Phys. Rev. B 103, 024203 (2021).
  * Sels and Polkovnikov [2021] D. Sels and A. Polkovnikov, Dynamical obstruction to localization in a disordered spin chain, Phys. Rev. E 104, 054105 (2021).
  * Vidmar  _et al._ [2021] L. Vidmar, B. Krajewski, J. Bonča, and M. Mierzejewski, Phenomenology of spectral functions in disordered spin chains at infinite temperature, Phys. Rev. Lett. 127, 230603 (2021).
  * 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).
  * Iyer  _et al._ [2013] S. Iyer, V. Oganesyan, G. Refael, and D. A. Huse, Many-body localization in a quasiperiodic system, Phys. Rev. B 87, 134202 (2013).
  * Khemani  _et al._ [2017a] V. Khemani, D. N. Sheng, and D. A. Huse, Two universality classes for the many-body localization transition, Phys. Rev. Lett. 119, 075702 (2017a).
  * Lee  _et al._ [2017] M. Lee, T. R. Look, S. P. Lim, and D. N. Sheng, Many-body localization in spin chain systems with quasiperiodic fields, Phys. Rev. B 96, 075146 (2017).
  * Nag and Garg [2017] S. Nag and A. Garg, Many-body mobility edges in a one-dimensional system of interacting fermions, Phys. Rev. B 96, 060203 (2017).
  * Dutta  _et al._ [2018] A. Dutta, S. Mukerjee, and K. Sengupta, Many-body localized phase of bosonic dipoles in a tilted optical lattice, Phys. Rev. B 98, 144205 (2018).
  * Zhang and Yao [2018] S.-X. Zhang and H. Yao, Universal properties of many-body localization transitions in quasiperiodic systems, Phys. Rev. Lett. 121, 206601 (2018).
  * Agrawal  _et al._ [2020] U. Agrawal, S. Gopalakrishnan, and R. Vasseur, Universality and quantum criticality in quasiperiodic spin chains, Nat. Commun. 11, 2225 (2020).
  * Macé  _et al._ [2019] N. Macé, N. Laflorencie, and F. Alet, Many-body localization in a quasiperiodic Fibonacci chain, SciPost Phys. 6, 50 (2019).
  * Singh  _et al._ [2021] H. Singh, B. Ware, R. Vasseur, and S. Gopalakrishnan, Local integrals of motion and the quasiperiodic many-body localization transition, Phys. Rev. B 103, L220201 (2021).
  * Duthie  _et al._ [2021] A. Duthie, S. Roy, and D. E. Logan, Localization in quasiperiodic chains: A theory based on convergence of local propagators, Phys. Rev. B 104, 064201 (2021).
  * Sierant and Zakrzewski [2019] P. Sierant and J. Zakrzewski, Level statistics across the many-body localization transition, Phys. Rev. B 99, 104205 (2019).
  * Aramthottil  _et al._ [2021] A. S. Aramthottil, T. Chanda, P. Sierant, and J. Zakrzewski, Finite-size scaling analysis of the many-body localization transition in quasiperiodic spin chains, Phys. Rev. B 104, 214201 (2021).
  * Thomson and Schiró [2021] S. J. Thomson and M. Schiró, Local integrals of motion in quasiperiodic many-body localized systems, arXiv:2110.02906 (2021).
  * Žnidarič and Ljubotina [2018] M. Žnidarič and M. Ljubotina, Interaction instability of localization in quasiperiodic systems, Proc. Natl. Acad. Sci. USA 115, 4595 (2018).
  * Szabó and Schneider [2020] A. Szabó and U. Schneider, Mixed spectra and partially extended states in a two-dimensional quasiperiodic model, Phys. Rev. B 101, 014205 (2020).
  * Doggen and Mirlin [2019] E. V. H. Doggen and A. D. Mirlin, Many-body delocalization dynamics in long Aubry-André quasiperiodic chains, Phys. Rev. B 100, 104203 (2019).
  * Sierant and Zakrzewski [2022] P. Sierant and J. Zakrzewski, Challenges to observation of many-body localization, Phys. Rev. B 105, 224203 (2022).
  * Khemani  _et al._ [2019] V. Khemani, C. R. Laumann, and A. Chandran, Signatures of integrability in the dynamics of Rydberg-blockaded chains, Phys. Rev. B 99, 161101 (2019).
  * Huse  _et al._ [2013] D. A. Huse, R. Nandkishore, V. Oganesyan, A. Pal, and S. L. Sondhi, Localization-protected quantum order, Phys. Rev. B 88, 014206 (2013).
  * Chen  _et al._ [2020] C. Chen, X. Wang, and Y. Chen, Lieb-Robinson bound for constrained many-body localization, arXiv:2011.11363 (2020).
  * Schreiber  _et al._ [2015] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch, Observation of many-body localization of interacting fermions in a quasirandom optical lattice, Science 349, 842 (2015).
  * Bordia  _et al._ [2016] P. Bordia, H. P. Lüschen, S. S. Hodgman, M. Schreiber, I. Bloch, and U. Schneider, Coupling identical one-dimensional many-body localized systems, Phys. Rev. Lett. 116, 140401 (2016).
  * 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).
  * Oganesyan and Huse [2007] V. Oganesyan and D. A. Huse, Localization of interacting fermions at high temperature, Phys. Rev. B 75, 155111 (2007).
  * Atas  _et al._ [2013] Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, Distribution of the ratio of consecutive level spacings in random matrix ensembles, Phys. Rev. Lett. 110, 084101 (2013).
  * Giraud  _et al._ [2022] O. Giraud, N. Macé, E. Vernier, and F. Alet, Probing symmetries of quantum many-body systems through gap ratio statistics, Phys. Rev. X 12, 011006 (2022).
  * De Roeck and Huveneers [2017] W. De Roeck and F. Huveneers, Stability and instability towards delocalization in many-body localization systems, Phys. Rev. B 95, 155129 (2017).
  * Yu  _et al._ [2016] X. Yu, D. J. Luitz, and B. K. Clark, Bimodal entanglement entropy distribution in the many-body localization transition, Phys. Rev. B 94, 184202 (2016).
  * Page [1993] D. N. Page, Average entropy of a subsystem, Phys. Rev. Lett. 71, 1291 (1993).
  * Kjäll  _et al._ [2014] J. A. Kjäll, J. H. Bardarson, and F. Pollmann, Many-body localization in a disordered quantum Ising chain, Phys. Rev. Lett. 113, 107204 (2014).
  * Khemani  _et al._ [2017b] 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 (2017b).
  * Vidal [2004] G. Vidal, Efficient simulation of one-dimensional quantum many-body systems, Phys. Rev. Lett. 93, 040502 (2004).
  * Smith  _et al._ [2016] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, Many-body localization in a quantum simulator with programmable random disorder, Nat. Phys. 12, 907 (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).
  * 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).
  * 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).
  * Chandran  _et al._ [2015] A. Chandran, I. H. Kim, G. Vidal, and D. A. Abanin, Constructing local integrals of motion in the many-body localized phase, Phys. Rev. B 91, 085425 (2015).
  * Schollwöck [2011] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Annals of Physics 326, 96 (2011).
  * Popescu  _et al._ [2006] S. Popescu, A. J. Short, and A. Winter, Entanglement and the foundations of statistical mechanics, Nat. Phys. 2, 754 (2006).
  * Goldstein  _et al._ [2006] S. Goldstein, J. L. Lebowitz, R. Tumulka, and N. Zanghì, Canonical typicality, Phys. Rev. Lett. 96, 050403 (2006).
  * Steinigeweg  _et al._ [2016] R. Steinigeweg, J. Herbrych, F. Pollmann, and W. Brenig, Typicality approach to the optical conductivity in thermal and many-body localized phases, Phys. Rev. B 94, 180401 (2016).