# Phenomenology of many-body localization in bond-disordered spin chains

Adith Sai Aramthottil  Szkoła Doktorska Nauk Ścisłych i Przyrodniczych,
Uniwersytet Jagielloński, Łojasiewicza 11, PL-30-348 Kraków, Poland  Instytut
Fizyki Teoretycznej, Uniwersytet Jagielloński, Łojasiewicza 11, PL-30-348
Kraków, Poland  Piotr Sierant  ICFO-Institut de Ciències Fotòniques, The
Barcelona Institute of Science and Technology, Av. Carl Friedrich Gauss 3,
08860 Castelldefels (Barcelona), Spain  Maciej Lewenstein  ICFO-Institut de
Ciències Fotòniques, The Barcelona Institute of Science and Technology, Av.
Carl Friedrich Gauss 3, 08860 Castelldefels (Barcelona), Spain  ICREA, Passeig
Lluis Companys 23, 08010 Barcelona, Spain  Jakub Zakrzewski  Instytut Fizyki
Teoretycznej, Uniwersytet Jagielloński, Łojasiewicza 11, PL-30-348 Kraków,
Poland  Mark Kac Complex Systems Research Center, Uniwersytet Jagielloński,
PL-30-348 Kraków, Poland

###### Abstract

Many-body localization (MBL) hinders the thermalization of quantum many-body
systems in the presence of strong disorder. In this work, we study the MBL
regime in bond-disordered spin-1/2 XXZ spin chain, finding the multimodal
distribution of entanglement entropy in eigenstates, sub-Poissonian level
statistics, and revealing a relation between operators and initial states
required for examining the breakdown of thermalization in the time evolution
of the system. We employ a real space renormalization group scheme to identify
these phenomenological features of the MBL regime that extend beyond the
standard picture of local integrals of motion relevant for systems with
disorder coupled to on-site operators. Our results pave the way for
experimental probing of MBL in bond-disordered spin chains.

#### Introduction.

The Eigenstate thermalization hypothesis (ETH) [1, 2, 3, 4] conjectures that
quantum many-body systems thermalize and reach equilibrium state determined
solely by a few macroscopic quantities, irrespective of the details of the
initial state. Many-body localization (MBL) [5, 6, 7, 8, 9, 10, 11] has been
proposed as a mechanism to avoid thermalization for a large class of strongly
disordered systems, including spin [12, 13, 14], bosonic [15, 16, 17, 18],
fermionic [19, 20, 21, 22] and Floquet systems [23, 24, 25, 26, 27, 28, 29,
30, 31, 32].

Figure 1:  The RSRG-X scheme approximates the eigenstates and eigenvalues of
the bond-disorder XXZ spin chain, capturing the phenomenology of the MBL
regime at large disorder strength $W$. (a) RSRG-X identifies and decimates the
bonds between spins approximating eigenvalues ($\epsilon$) for a single
realization of Hamiltonian (1) with $M=8$ spins, see panel (b). (c) Properties
captured by the RSRG-X procedure: (i) nontrivial entanglement entropy
distribution, (ii) sub-Poissonian mean gap ratio, (iii) lack of thermalization
of $\hat{\mathcal{O}}_{n}$ at large $W$.

Explorations of MBL have been primarily driven by numerical studies restricted
by the computationally accessible system sizes [33, 34] and timescales [35,
36, 37]. Consequently, it is unclear whether the observed MBL regime remains
stable in the asymptotic limit of infinite system size and infinite time [38,
39, 40, 41], motivating studies of persistent slow dynamics and prevalent
system size drifts of ergodicity breaking indicators in strongly disordered
interacting systems [42, 43, 44, 45, 46, 47, 48, 49, 50], as reviewed in [11].

The MBL regime is characterized, on experimentally relevant timescales and
system sizes, by the slow growth of entanglement in dynamics [51, 52, 53, 54],
power-law decay of local correlation functions [55, 56], Poissonian level
statistics [57, 58, 59, 60, 61], and area-law entanglement of eigenstates [62,
63]. Local integrals of motions (LIOMs) [64, 65] provide a unifying framework
that captures the phenomenology of MBL [62, 64, 10] in the vast class of
systems [12, 13, 14, 19, 20, 15, 16, 17, 18, 21, 22, 23, 24, 25, 26, 27, 28,
29, 30, 31, 32]. The LIOMs in the MBL regime can be constructed as long-time
averages of local operators [66], with perturbative approaches [67], with flow
equation methods [68, 69, 70, 71] or direct numerical optimization [72, 73,
74, 75]. In the paradigmatic disordered XXZ spin chain, LIOMs,
$\hat{\tau}_{i}^{z}$, are the $z$-component of the spin operator
$\hat{s}^{z}_{i}$ transformed by a quasilocal unitary [67, 76, 10]. This
paradigm applies whenever the disorder is coupled to on-site operators, which
become exact LIOMs in the strong disorder limit.

This work focuses on many-body systems in which the disorder is not coupled to
on-site operators. We argue that the standard picture of LIOMs does not
capture the phenomenology of MBL in such systems. Instead, the real space
renormalization group for excited states (RSRG-X) [77] provides a
phenomenological description of the system’s behavior at large disorder
strengths. As a concrete example, we consider a random-bond XXZ spin-1/2 chain
relevant for experiments with ultracold Rydberg atoms [78, 79, 80, 81]. For
the considered model, the RSRG-X predicts sub-Poissonian level statistics,
multimodal structure of the entanglement entropy distribution in eigenstates
at large disorder strengths, and allows for identification of operators that
may be conserved in the time evolution of experimentally relevant initial
states and hence serve as order parameters for the MBL regime, see Fig. 1. We
show that the non-trivial memory effects of the RSRG-X procedure prevent us
from finding a global set of LIOMs for the considered system even in the
infinite disorder limit. We discuss how to extend our procedure to obtain
experimentally verifiable predictions for other types of many-body systems.

#### The model.

We consider the random-bond XXZ spin chain

$\hat{\mathcal{H}}=\sum_{i=1}^{M-1}J_{i,i+1}(\hat{s}_{i}^{x}\hat{s}_{i+1}^{x}+\hat{s}_{i}^{y}\hat{s}_{i+1}^{y}+\Delta_{i,i+1}\hat{s}_{i}^{z}\hat{s}_{i+1}^{z})$
(1)

with open boundary conditions, where, $\hat{s}^{\beta}_{i}$ are the spin-$1/2$
operators acting on the $i$ site with $\beta\in\\{x,y,z\\}$, $J_{i,i+1}$ are
couplings and $M$ is the number of spins. We focus exclusively on the zero-
magnetization sector and avoid the special case of the non-abelian SU(2)
symmetry by imposing $|\Delta_{i,i+1}|\neq 1$ [82]. The remaining symmetry of
the model is $\mathbb{Z}_{2}$ spin inversion symmetry $\mathcal{C}$. The
Hamiltonian (1) is relevant for experiments with Rydberg atoms confined in an
array of microtraps [83, 84]. Investigations of these models have demonstrated
a weak breaking of the ETH in the presence of interactions [81] and glassy
dynamics of selected operators [79, 80]. With that motivation in mind, we
assume that the couplings in (1) take values according to a power-law
$J_{i,i+1}=|x_{i}-x_{i+1}|^{-\alpha}$, where $x_{i}$ denotes the position of
the $i$-th spin within the length $L$ occupied by the system, and
$\Delta=-0.73$, $\alpha=6$ relevant for a Rydberg state of Rubidium 87. Fixing
the Rydberg blockade [85] diameter as $D_{rb}$ results in an excluded density
$\rho=D_{rb}M/L$. With decreasing density $\rho$, there are more possible
configurations present, leading to an increase of the effective disorder
strength $W=1/\rho-1$ in the model [86] and thus to the MBL regime. In S-II,
we describe our algorithm to generate particle configurations $\\{x_{i}\\}$
avoiding the issue of the Rényi parking constant [87].

#### Decimation Scheme.

To understand the MBL phenomenology of (1), we employ the RSRG-X scheme [77].
The scheme proceeds iteratively by identifying the largest energy gap,
diagonalizing the corresponding local operator, setting the state of the
relevant spins into one of the obtained eigenstates, and perturbatively
modifying the nearby effective local operators and the associated energy gaps
[88, 77]. The energy gaps encountered at each step of this procedure, which we
refer to as a decimation, are progressively smaller. Hence, the chosen
eigenstates fix the time dynamics at longer times, eventually producing an
approximation of the eigenstate of (1) when the last spin is fixed. If the
distribution of energy gaps is steadily broadening in the subsequent
decimation steps, the scheme is increasingly well controlled and may become
asymptotically exact [89, 90, 91, 92]. Here, relying on the broadness of the
initial gap distribution, the RSRG-X is employed to obtain insights about
relatively small systems relevant for numerical and experimental studies.

The RSRG-X proceeds as follows. The Hamiltonian (1) can be written as a sum of
2-site local operators of the form

$\hat{\mathcal{H}}_{i,i+1}=J_{i,i+1}\biggl{(}\frac{(\hat{s}_{i}^{+}\hat{s}_{i+1}^{-}+h.c)}{2}+\Delta_{i,i+1}\hat{s}_{i}^{z}\hat{s}_{i+1}^{z}\biggl{)},$
(2)

with eigenbasis $|Z_{\pm}\rangle,|\pm\rangle$ where,
$|Z_{+/-}\rangle\equiv|\uparrow\uparrow\rangle$/$|\downarrow\downarrow\rangle$,
and
$|\pm\rangle\equiv\frac{1}{\sqrt{2}}(|\uparrow\downarrow\rangle\pm|\downarrow\uparrow\rangle)$.
The relevant energy gap associated with each local operator is taken as the
maximum of its energy gaps, so that any of the eigenstates can be chosen
freely. This is a valid approximation if the energy scales are controlled by
the couplings $J_{i,i+1}$. Let us assume that the largest relevant energy gap
is associated with the local operator $\hat{\mathcal{H}}_{4,5}$. In that case,
we select the bond between spins $i=3$, $i=4$ and fix the state of the two
spins as $|Z_{\pm}\rangle,|\pm\rangle$. We now describe how the neighboring
operators, i.e., $\hat{\mathcal{H}}_{2,3}$ and $\hat{\mathcal{H}}_{4,5}$ are
perturbatively modified for each of the choices.

For states $|\pm\rangle$, there are no contributions to the nearest operators
within the first-order perturbation. The second-order perturbation results in
the formation of a bond between spins $i=2$, $i=5$, described by
$\hat{\mathcal{H}}_{2,5}$ of the form as in (2) with couplings modified to
$\tilde{J}_{2,5}=\frac{J_{2,3}J_{4,5}}{(1\mp\Delta_{3,4})J_{3,4}}$,
$\tilde{\Delta}_{2,5}=\frac{\mp\Delta_{2,3}\Delta_{4,5}(1\mp\Delta_{3,4})}{2}$.
The choice $|Z_{\pm}\rangle$, within the first-order perturbation modifies the
operators $\hat{\mathcal{H}}_{2,3}$ and $\hat{\mathcal{H}}_{4,5}$ to take a
chemical potential form [93]

$\hat{\mathcal{H}}_{2,3}^{\prime}=\pm\frac{1}{2}(\Delta_{2,3}J_{2,3}\hat{s}_{2}^{z}),\hat{\mathcal{H}}_{4,5}^{\prime}=\pm\frac{1}{2}(\Delta_{4,5}J_{4,5}\hat{s}_{5}^{z}),$
(3)

while the second order term leads to a bond between $i=2,5$ with
$\tilde{J}_{2,5}=-\frac{J_{2,3}J_{4,5}}{2J_{3,4}}\left(\frac{1}{(1-\Delta_{3,4})}+\frac{1}{(1+\Delta_{3,4})}\right)$
and $\tilde{\Delta}_{2,5}=0$. The first-order modification of
$\hat{\mathcal{H}}^{\prime}$ is typically much stronger than the second-order
bond and the neighboring operators admit the chemical potential form in the
following steps of the scheme. For instance, a decimation of the bond between
sites $i=4$, $i=5$, would result in fixing the eigenstate at site $5$ to
$|\uparrow\rangle/|\downarrow\rangle$ and modifying the next operator as
$\hat{\mathcal{H}}_{5,6}^{\prime}=\frac{1}{2}(\Delta_{5,6}J_{5,6}\hat{s}_{6}^{z})/-\frac{1}{2}(\Delta_{5,6}J_{5,6}\hat{s}_{6}^{z})$.
This may lead to a process of growth of domains of spins fixed in states
$|\uparrow\rangle/|\downarrow\rangle$ which stops if the next operator to be
modified already has the form $\hat{\mathcal{H}}^{\prime}$. The spin polarized
domains $|\uparrow\rangle/|\downarrow\rangle$ may extend over already
decimates pair of sites with states fixed as $|\pm\rangle$.

Figure 2:  (a) Distribution ${P}_{\mathcal{S}}$ of half-chain entanglement
entropy ${\mathcal{S}}$ for $M=20$ and varying disorder strength $W$. (b) The
mean gap ratio obtained by RSRG-X ($\bar{r}_{RG}$, blue) is compared with ED
results ($\bar{r}_{ED}$, orange). The inset shows
$\delta_{\bar{r}}=\bar{r}_{ED}-\bar{r}_{RG}$. The system size is $M=10$, with
$25000$ disorder realizations and all the eigenvalues are considered. (c) The
rescaled energies $\epsilon$ obtained with ED are compared with RSRG-X result
for a single disorder realization, (d) is a zoomed image of the same. All the
results for the $\mathcal{C}=+1$ symmetry sector.

The generation of $\hat{\mathcal{H}^{\prime}}$ operators through the first
order perturbation suffices in all scenarios except when degenerate states
$|K_{+}\rangle\equiv|\uparrow\rangle|Z_{-}\rangle|\uparrow\rangle$,
$|K_{-}\rangle\equiv|\downarrow\rangle|Z_{+}\rangle|\downarrow\rangle$ are
formed. Then, the second order perturbation needs to be taken into account
leading to operators of the form $\hat{\mathcal{H}}$ with modified parameters,
see S-III.

A generic approximation of an eigenstate stemming from RSRG-X is associated
with arbitrary choices of $|Z_{\pm}\rangle$, and $|\pm\rangle$ throughout the
decimation scheme. To qualitatively understand the properties of these
eigenstates, we start by assuming that only states $|\pm\rangle$ have been
chosen in the decimation scheme. This leads to formation of bonds
$\hat{\mathcal{H}}_{l_{i},r_{i}}$ between arbitrarily distant spins, resulting
in states $|\pm\rangle_{l_{i},r_{i}}$. The spatial distribution of the
correlations of states within this subsystem is analogous to a random singlet
phase found for the ground state of (1) [88, 89], or for the excited states in
the absence of interactions ($\Delta=0$) [94, 95]. In contrast, fixing of the
eigenstate as $|Z_{\pm}\rangle$ in the decimation procedure leads to growth of
a domain which breaks the spin inversion symmetry
$|\uparrow/\downarrow\rangle_{l_{m}}\cdots|\uparrow/\downarrow\rangle_{l_{2}}|Z_{+/-}\rangle_{l_{1},r_{1}}|\uparrow/\downarrow\rangle_{r_{2}}\cdots|\uparrow/\downarrow\rangle_{r_{m^{\prime}}}$
(except for the configurations involving $|K_{+}\rangle$ and $|K_{-}\rangle$).
When the energy scales are not controlled by the couplings $J$, i.e., when
$\Delta\gg 1$, such spin inversion symmetry breaking domains extend through
the entire system [93]. In contrast, in our case, a generic eigenstate is
formed by domains with spins polarized along the $z$ direction and subsystems
in which the state is a product of triplet/singlet states $|\pm\rangle$ 111
With probability $\approx 0.02$, the illustrated subsystem of states
$|\pm\rangle$ contains states of the form
$\frac{1}{\sqrt{2}}(|K_{+}\rangle\pm|K_{-}\rangle)$ .

#### Eigenstates and eigenvalues.

To demonstrate that the eigenstates of the random-bond XXZ spin chain (1) at
large $W$ indeed feature the structure predicted by the RSRG-X, we analyze the
half-chain entanglement entropy
$\mathcal{S}=-\mathrm{tr}(\rho\log_{2}(\rho))$, where
$\rho=\mathrm{tr_{B}}(\ket{\psi}\bra{\psi})$ is the reduced density obtained
by tracing out the spins in half of the system and $\ket{\psi}$ is the
considered eigenstate. The distribution of $\mathcal{S}$, calculated by exact
diagonalization (ED) and shown in Fig. 2(a), exhibits peaks at integer values.
The peaks correspond to different numbers of the bonds in triplet/singlet
states $\ket{\pm}$ cut by the bipartition of the system, confirming the
intuitive picture predicted by the RSRG-X. The peaks get broadened at smaller
disorder strengths $W$ due to the contribution of the higher order terms not
taken into the account by our decimation scheme. The observed distribution of
$\mathcal{S}$ is starkly different from the XXZ spin-1/2 chain which, in the
large $W$ limit, collapses to a single peak at $\mathcal{S}=0$ [97].

The RSRG-X also approximates the eigenvalues $\epsilon_{k}$ of (1). We compare
the RSRG-X prediction with the ED results by calculating the gap ratio
$r_{k}\equiv\frac{min\\{\delta_{k},\delta_{k+1}\\}}{max\\{\delta_{k},\delta_{k+1}\\}}$
(here $\delta_{k}\equiv\epsilon_{k+1}-\epsilon_{k}$) averaged over the entire
spectrum of the system and disorder realizations, shown in Fig. 2(b). The
RSRG-X qualitatively reproduces the behavior of $\bar{r}$, and, in particular,
predicts the sub-Poissonian value $\bar{r}<0.386=\bar{r}_{PS}$ at larger $W$,
where $\bar{r}_{PS}\approx 0.386$ [98] is the value for Poissonian level
statistics found in the MBL regime of models with a complete set of LIOMs
[10]. The sub-Poissonian level statistics arises due to multiple
quasidegeneracies in the spectrum, qualitatively captured by the RSRG-X
procedure. A closer inspection reveals that the $\bar{r}$ from the RSRG-X is
systematically underestimating the ED value, which we attribute to higher
order processes neglected by the RSRG-X. In S-III, we show that this
discrepancy is smaller when the distribution of $J$’s is wider, i.e., given by
a power-law with $\alpha>6$.

#### Memory effects along the RSRG-X branch.

Each step $n$ of the RSRG-X is associated with choice of the states
$|Z_{\pm}\rangle$, $|\pm\rangle$ or, in case of the single-site decimations,
with $\ket{\uparrow/\downarrow}$. While fixing $\ket{\pm}$ results only in
decimation of the involved spins and renormalization of the parameters of the
bond linking their neighbors, the selection of $\ket{Z_{\pm}}$ affects non-
trivially the subsequent steps of the scheme. Each choice of $\ket{Z_{\pm}}$
introduces a possibility of the growth of the domain of spins polarized along
$z$ direction. This may arise whenever the interactions are present in the
system, i.e., when $|\Delta_{i,i+1}|>0$. The triplet/singlet states
$\ket{\pm}$ can not extend over the already fixed spin polarized domains.
Consequently, various branches of RSRG-X lead to eigenstates corresponding to
different partitions of the system into spin polarized and triplet/singlet
domains. Therefore, in contrast to the models in which the disorder is coupled
to the on-site operators, we cannot find a complete set of (quasi-)local
operators that would commute with the Hamiltonian (1).

Figure 3:  (a) The decimation scheme for one configuration of spins at $W=3$
and $M=16$; connected spins represent decimation steps with $\ket{\pm}$ (in
red) or $\ket{Z_{\pm}}$ (in green) fixed, and the numbers denote the number
$n$ of the step; the blue box denotes an ergodic bubble (b) Time evolution of
the average values $\mathcal{O}_{n}$ of operators corresponding to the RSRG-X
branch shown in (a). (c) Time evolution of $\mathcal{O}_{n}(t)$ for the step
$n=5$; the bold (dashed) line corresponds to $W=3$ ($W=7$).

#### Probing the MBL regime with time dynamics.

Time evolution of the average value $s(t)$ of $\hat{s}^{z}_{i}$ operator is
experimentally realizable way of studying the breakdown of thermalization in
the standard disordered XXZ spin chain. For an initial state taken as the
eigenstate $\ket{\psi_{Z}}$ of $\hat{s}^{z}_{i}$ operators, $s(t)$ decays to
zero in the ergodic regime, while, deep in the MBL regime, for large $W$,
$s(t)$ remains close to unity even at long times [55, 56, 50]. This procedure,
however, would not work for the bond-disordered XXZ spin chain. Due to the
non-trivial relations between decimations along various branches of RSRG-X and
the presence of eigenstates at any energy density with triplet/singlet
domains, the average $s(t)$ may dephase to zero irrespective of the disorder
strength $W$.

To circumvent this issue, we propose to set as an initial state,
$\ket{\psi_{\mathrm{in}}}$, the approximation of eigenstate obtained with
RSRG-X and to measure average values of operators $\hat{\mathcal{O}}_{n}$
related to the RSRG-X branch leading to $\ket{\psi_{\mathrm{in}}}$. For
decimations involving $|\pm\rangle$ we choose
$\hat{\mathcal{O}}_{n}=\hat{T}_{i,j}\equiv(\hat{s}_{i}^{+}\hat{s}_{j}^{-}+h.c.)$,
while for $|Z_{\pm}\rangle$, we select
$\hat{\mathcal{O}}_{n}=\hat{P}_{i,j}\equiv 4\hat{s}_{i}^{z}\hat{s}_{j}^{z}$
222In rare situation of a decimation with
$\frac{1}{\sqrt{2}}(|K_{+}\rangle\pm|K_{-}\rangle)$, we consider as operators
a combination of two $\hat{P}_{i,j}$.. For single-site decimations, we
consider the operator $\hat{\mathcal{O}}_{n}=2\hat{s}_{i}^{z}$. The obtained
average values $\mathcal{O}_{n}$ decay to zero in the ergodic regime, while
remain close to their starting value as presented in Fig. 1(c)(iii). In Fig.
3(a), we show a particular position configuration and decimation scheme along
a selected RSRG-X branch. The initial state has a simple tensor product
structure
$\ket{\psi_{\mathrm{in}}}=|Z_{+}\rangle^{1}\otimes|-\rangle^{2}\otimes|-\rangle^{3}\otimes|+\rangle^{4}\otimes|+\rangle^{5}\otimes|-\rangle^{6}\otimes|-\rangle^{7}\otimes|Z_{-}\rangle^{8}$,
where the superscripts represent the decimation step $n$. In Fig. 3(b), we
present the time evolution of the operators $\hat{\mathcal{O}}_{n}$ and find
that the expectation values are preserved except for $\hat{\mathcal{O}}_{3}$
and $\hat{\mathcal{O}}_{5}$ which reveal oscillations. The oscillations are
due to a local narrow distribution of $J$s, represented by a box in Fig. 3(a).
This is an example of an ergodic bubble [100, 101, 102] where local disorder
is weaker than its average over the whole system. Fig. 3(c) details the time
oscillations at the ergodic bubble, whose amplitude and frequency decrease
with increase of $W$.

To further demonstrate the validity of the procedure to probe the breakdown of
ergodicity in time dynamics, we consider a long-time averaged expectation
value $|\bar{\mathcal{O}}_{n}|=|\int_{0}^{T}\mathcal{O}_{n}(t)dt|$. Disorder
averaged $\langle|\bar{\mathcal{O}}_{n}|\rangle$ serves as ergodicity breaking
indicator. As shown in Fig. 4(a), $\langle|\bar{\mathcal{O}}_{n}|\rangle$
decays to zero in the ergodic regime, where the mean gap ratio $\bar{r}$ (for
detailed finit-size scaling analysis, see S-II) admits the GOE value
$\bar{r}_{GOE}\approx 0.53$ characteristic for ergodic systems. In contrast,
in the MBL regime, where $\bar{r}$ approaches (or is smaller than) the
Poissonian value $\bar{r}_{PS}$, the $\langle|\bar{\mathcal{O}}_{n}|\rangle$
is close to unity, showing clearly the lack of thermalization in the system.

Figure 4: (a) Long-time (up to $t=5000$) averages
$\langle|\bar{\mathcal{O}}|\rangle$ of operators found in the RSRG-X procedure
(see text), averaged over $10$ random RSRG-X branches and over $50$ disorder
realizations. (b) Mean gap ratio $\bar{r}$ as a function of disorder strength
$W$ for number of spins $M$ with $\mathcal{C}=+1$

.

#### Generalization.

Let us now consider a generic Hamiltonian of the form
$H=\sum_{[\mu]}h_{[\mu]}\otimes_{i}\hat{s}^{\mu_{i}}$, where each $h_{[\mu]}$
is drawn from a certain, sufficiently broad, distribution [103], and
$\otimes_{i}\hat{s}^{\mu_{i}}$ is a local operator. In general, the disorder
is not coupled to the on-site operators, hence the standard LIOMs are
insufficient to grasp the phenomenological properties of the system at large
$W$. Instead, in the absence of any special symmetry that would lead to
resonances, see e.g. [82], the RSRG-X scheme would suffice to identify the
properties of the MBL regime in the system. Decimating the largest energy gap,
a second-order perturbation theory would lead to virtual connections between
the neigboring subsystems, giving rise to non-local operators conserved at
large $W$. The first-order perturbation theory would fix the states of the
nearby subsystems and could suppress the long-range correlations.

#### Conclusion.

We studied the MBL regime arising at large disorder strengths in the bond-
disordered spin-1/2 XXZ chain. The phenomenology of the MBL regime in the
considered system differs significantly from models with disorder coupled to
on-site operators, widely studied in the context of MBL [10, 11].
Consequently, the standard picture of LIOMs does not describe the MBL regime
in the bond-disordered XXZ spin chain, which is characterized by multimodal
distribution of entanglement entropy in eigenstates, sub-Poissonian level
statistics, non-trivial interplay of operators and the initial state, which
has to be accounted for when probing the breakdown of ergodicity with time
dynamics. The RSRG-X scheme provides a qualitative understanding of the MBL
regime by approximating the highly excited eigenstates as products of domains
of spins polarized along the $z$ direction and regions in which spins are
paired into triplet/singlet states. The interplay between the two types of
domains leads to non-trivial memory effects in the RSRG-X scheme, disallowing
the identification of a global set of LIOMs for the model. We argued that the
predictions of RSRG-X become quantitatively valid at strong disorder for
system sizes accessible to exact diagonalization, provided that the initial
distribution of couplings in the model is sufficiently broad.

The bond-disordered XXZ spin-1/2 chain is relevant for Rydberg atoms
platforms, which provide a setting to experimentally test our results. The
idea of applying the RSRG-X scheme to extract phenomenological properties of
MBL regime can be straightforwardly extended to other systems in which the
disorder is not coupled to on-site operators. Another question involves
examination of the other phenomenological features of the MBL regime in bond-
disordered spin chains. For instance, the participation entropies that
characterize the spread of eigenstates in the eigenbasis $\hat{s}_{i}^{z}$,
see e.g. [104], would not be vanishing even in the $W\to\infty$ limit due to
the presence of triplet/singlet domains in the eigenstates. These questions
are left for further studies.

###### Acknowledgements.

The work of A.S.A. has been realized within the Opus grant
2019/35/B/ST2/00034, financed by National Science Centre (Poland). P.S.
acknowledges support from: ERC AdG NOQIA; MICIN/AEI
(PGC2018-0910.13039/501100011033, CEX2019-000910-S/10.13039/501100011033, Plan
National FIDEUA PID2019-106901GB-I00, FPI; MICIIN with funding from European
Union NextGenerationEU (PRTR-C17.I1): QUANTERA MAQS PCI2019-111828-2);
MCIN/AEI/ 10.13039/501100011033 and by the “European Union NextGeneration
EU/PRTR" QUANTERA DYNAMITE PCI2022-132919 within the QuantERA II Programme
that has received funding from the European Union’s Horizon 2020 research and
innovation programme under Grant Agreement No 101017733Proyectos de I+D+I
“Retos Colaboración” QUSPIN RTC2019-007196-7); Fundació Cellex; Fundació Mir-
Puig; Generalitat de Catalunya (European Social Fund FEDER and CERCA program,
AGAUR Grant No. 2021 SGR 01452, QuantumCAT U16-011424, co-funded by ERDF
Operational Program of Catalonia 2014-2020); Barcelona Supercomputing Center
MareNostrum (FI-2024-1-0043); EU (PASQuanS2.1, 101113690); EU Horizon 2020
FET-OPEN OPTOlogic (Grant No 899794); EU Horizon Europe Program (Grant
Agreement 101080086 — NeQST), ICFO Internal “QuantumGaudi” project; European
Union’s Horizon 2020 research and innovation program under the Marie-
Skłodowska-Curie grant agreement No 101029393 (STREDCH) and No 847648 (“La
Caixa” Junior Leaders fellowships ID100010434: LCF/BQ/PI19/11690013,
LCF/BQ/PI20/11760031, LCF/BQ/PR20/11770012, LCF/BQ/PR21/11840013). E.P. is
supported by “Ayuda (PRE2021-098926) financiada por MCIN/AEI/
10.13039/501100011033 y por el FSE+". The work of J.Z. was funded by the
National Science Centre, Poland, project 2021/03/Y/ST2/00186 within the
QuantERA II Programme that has received funding from the European Union
Horizon 2020 research and innovation programme under Grant agreement No
101017733. A partial support by the Strategic Programme Excellence Initiative
at Jagiellonian University is acknowledged. Views and opinions expressed in
this work are, however, those of the authors only and do not necessarily
reflect those of the European Union, European Climate, Infrastructure and
Environment Executive Agency (CINEA), nor any other granting authority.
Neither the European Union nor any granting authority can be held responsible
for them. No part of this work was written with the help of AI.

## References

  * 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).
  * Srednicki [1999] M. Srednicki, The approach to thermal equilibrium in quantized chaotic systems, Journal of Physics A: Mathematical and General 32, 1163 (1999).
  * D'Alessio  _et al._ [2016] L. D'Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics, Advances in Physics 65, 239 (2016).
  * 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).
  * Pal and Huse [2010] A. Pal and D. A. Huse, Many-body localization phase transition, Phys. Rev. B 82, 174411 (2010).
  * Nandkishore and Huse [2015] R. Nandkishore and D. A. Huse, Many-Body Localization and Thermalization in Quantum Statistical Mechanics, Annual Review of Condensed Matter Physics 6, 15 (2015).
  * Alet and Laflorencie [2018] F. Alet and N. Laflorencie, Many-body localization: An introduction and selected topics, Comptes Rendus Physique 19, 498 (2018).
  * 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).
  * Sierant  _et al._ [2024] P. Sierant, M. Lewenstein, A. Scardicchio, L. Vidmar, and J. Zakrzewski, Many-body localization in the age of classical computing (2024), arXiv:2403.07111 [cond-mat.dis-nn] .
  * Santos  _et al._ [2004] L. F. Santos, G. Rigolin, and C. O. Escobar, Entanglement versus chaos in disordered spin chains, Phys. Rev. A 69, 042304 (2004).
  * Oganesyan and Huse [2007] V. Oganesyan and D. A. Huse, Localization of interacting fermions at high temperature, Phys. Rev. B 75, 155111 (2007).
  * 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).
  * Sierant  _et al._ [2017] P. Sierant, D. Delande, and J. Zakrzewski, Many-body localization due to random interactions, Phys. Rev. A 95, 021601 (2017).
  * Sierant and Zakrzewski [2018] P. Sierant and J. Zakrzewski, Many-body localization of bosons in optical lattices, New Journal of Physics 20, 043032 (2018).
  * Orell  _et al._ [2019] T. Orell, A. A. Michailidis, M. Serbyn, and M. Silveri, Probing the many-body localization phase transition with superconducting circuits, Phys. Rev. B 100, 134504 (2019).
  * Hopjan and Heidrich-Meisner [2020] M. Hopjan and F. Heidrich-Meisner, Many-body localization from a one-particle perspective in the disordered one-dimensional Bose-Hubbard model, Phys. Rev. A 101, 063617 (2020).
  * Mondaini and Rigol [2015] R. Mondaini and M. Rigol, Many-body localization and thermalization in disordered Hubbard chains, Phys. Rev. A 92, 041601 (2015).
  * Prelovšek  _et al._ [2016] P. Prelovšek, O. S. Barišić, and M. Žnidarič, Absence of full many-body localization in the disordered Hubbard chain, Phys. Rev. B 94, 241104 (2016).
  * Zakrzewski and Delande [2018] J. Zakrzewski and D. Delande, Spin-charge separation and many-body localization, Phys. Rev. B 98, 014203 (2018).
  * Kozarzewski  _et al._ [2018] M. Kozarzewski, P. Prelovšek, and M. Mierzejewski, Spin subdiffusion in the disordered Hubbard chain, Phys. Rev. Lett. 120, 246602 (2018).
  * Lazarides  _et al._ [2015] A. Lazarides, A. Das, and R. Moessner, Fate of many-body localization under periodic driving, Phys. Rev. Lett. 115, 030402 (2015).
  * Ponte  _et al._ [2015a] P. Ponte, Z. Papić, F. m. c. Huveneers, and D. A. Abanin, Many-body localization in periodically driven systems, Phys. Rev. Lett. 114, 140401 (2015a).
  * Ponte  _et al._ [2015b] P. Ponte, A. Chandran, Z. Papić, and D. A. Abanin, Periodically driven ergodic and many-body localized quantum systems, Annals of Physics 353, 196 (2015b).
  * Abanin  _et al._ [2016] D. A. Abanin, W. De Roeck, and F. Huveneers, Theory of many-body localization in periodically driven systems, Annals of Physics 372, 1 (2016).
  * Zhang  _et al._ [2016] L. Zhang, V. Khemani, and D. A. Huse, A floquet model for the many-body localization transition, Phys. Rev. B 94, 224202 (2016).
  * Bairey  _et al._ [2017] E. Bairey, G. Refael, and N. H. Lindner, Driving induced many-body localization, Phys. Rev. B 96, 020201 (2017).
  * Sahay  _et al._ [2021] R. Sahay, F. Machado, B. Ye, C. R. Laumann, and N. Y. Yao, Emergent ergodicity at the transition between many-body localized phases, Phys. Rev. Lett. 126, 100604 (2021).
  * Garratt and Chalker [2021] S. J. Garratt and J. T. Chalker, Many-body delocalization as symmetry breaking, Phys. Rev. Lett. 127, 026802 (2021).
  * Sonner  _et al._ [2021] M. Sonner, M. Serbyn, Z. Papić, and D. A. Abanin, Thouless energy across the many-body localization transition in floquet systems, Phys. Rev. B 104, L081112 (2021).
  * Sierant  _et al._ [2023] P. Sierant, M. Lewenstein, A. Scardicchio, and J. Zakrzewski, Stability of many-body localization in floquet systems, Physical Review B 107, 115132 (2023).
  * Pietracaprina  _et al._ [2018] F. Pietracaprina, N. Macé’, D. J. Luitz, and F. Alet, Shift-invert diagonalization of large many-body localizing spin chains, SciPost Phys. 5, 45 (2018).
  * Sierant  _et al._ [2020a] P. Sierant, M. Lewenstein, and J. Zakrzewski, Polynomially filtered exact diagonalization approach to many-body localization, Phys. Rev. Lett. 125, 156601 (2020a).
  * Bera  _et al._ [2017] S. Bera, G. De Tomasi, F. Weiner, and F. Evers, Density propagator for many-body localization: Finite-size effects, transient subdiffusion, and exponential decay, Phys. Rev. Lett. 118, 196801 (2017).
  * Doggen  _et al._ [2018] E. V. H. Doggen, F. Schindler, K. S. Tikhonov, A. D. Mirlin, T. Neupert, D. G. Polyakov, and I. V. Gornyi, Many-body localization and delocalization in large quantum chains, Phys. Rev. B 98, 174202 (2018).
  * Chanda  _et al._ [2020a] T. Chanda, P. Sierant, and J. Zakrzewski, Time dynamics with matrix product states: Many-body localization transition of large systems revisited, Phys. Rev. B 101, 035148 (2020a).
  * Šuntajs  _et al._ [2020] J. Šuntajs, J. Bonča, T. Prosen, and L. Vidmar, Ergodicity breaking transition in finite disordered spin chains, Phys. Rev. B 102, 064207 (2020).
  * Sierant  _et al._ [2020b] P. Sierant, D. Delande, and J. Zakrzewski, Thouless time analysis of anderson and many-body localization transitions, Phys. Rev. Lett. 124, 186601 (2020b).
  * Abanin  _et al._ [2021] D. Abanin, J. Bardarson, G. D. 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).
  * Panda  _et al._ [2020] R. K. Panda, A. Scardicchio, M. Schulz, S. R. Taylor, and M. Znidaric, Can we study the many-body localisation transition?, Europhysics Letters 128, 67003 (2020).
  * Šuntajs  _et al._ [2020] J. Šuntajs, J. Bonča, T. Prosen, and L. Vidmar, Quantum chaos challenges many-body localization, Phys. Rev. E 102, 062144 (2020).
  * Kiefer-Emmanouilidis  _et al._ [2020] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Evidence for unbounded growth of the number entropy in many-body localized phases, Phys. Rev. Lett. 124, 243601 (2020).
  * Kiefer-Emmanouilidis  _et al._ [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).
  * 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).
  * Sels and Polkovnikov [2021] D. Sels and A. Polkovnikov, Thermalization of dilute impurities in one dimensional spin chains (2021), 2105.09348 [quant-ph] .
  * Morningstar  _et al._ [2022] A. Morningstar, L. Colmenarez, V. Khemani, D. J. Luitz, and D. A. Huse, Avalanches and many-body resonances in many-body localized systems, Phys. Rev. B 105, 174205 (2022).
  * Sels [2022] D. Sels, Bath-induced delocalization in interacting disordered spin chains, Phys. Rev. B 106, L020202 (2022).
  * Krajewski  _et al._ [2022] B. Krajewski, L. Vidmar, J. Bonča, and M. Mierzejewski, Restoring ergodicity in a strongly disordered interacting chain, Phys. Rev. Lett. 129, 260601 (2022).
  * Sierant and Zakrzewski [2022] P. Sierant and J. Zakrzewski, Challenges to observation of many-body localization, Phys. Rev. B 105, 224203 (2022).
  * Ž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).
  * Nanduri  _et al._ [2014] A. Nanduri, H. Kim, and D. A. Huse, Entanglement spreading in a many-body localized system, Phys. Rev. B 90, 064201 (2014).
  * 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).
  * Serbyn  _et al._ [2014] M. Serbyn, Z. Papić, and D. A. Abanin, Quantum quenches in the many-body localized phase, Phys. Rev. B 90, 174302 (2014).
  * 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).
  * Serbyn and Moore [2016] M. Serbyn and J. E. Moore, Spectral statistics across the many-body localization transition, Phys. Rev. B 93, 041424 (2016).
  * Bertrand and García-García [2016] C. L. Bertrand and A. M. García-García, Anomalous thouless energy and critical statistics on the metallic side of the many-body localization transition, Phys. Rev. B 94, 144201 (2016).
  * Sierant and Zakrzewski [2019] P. Sierant and J. Zakrzewski, Level statistics across the many-body localization transition, Phys. Rev. B 99, 104205 (2019).
  * Buijsman  _et al._ [2019] W. Buijsman, V. Cheianov, and V. Gritsev, Random matrix ensemble for the level statistics of many-body localization, Phys. Rev. Lett. 122, 180601 (2019).
  * Sierant and Zakrzewski [2020] P. Sierant and J. Zakrzewski, Model of level statistics for disordered interacting quantum many-body systems, Phys. Rev. B 101, 104201 (2020).
  * Serbyn  _et al._ [2013a] M. Serbyn, Z. Papić, and D. A. Abanin, Universal slow growth of entanglement in interacting strongly disordered systems, Phys. Rev. Lett. 110, 260601 (2013a).
  * Bauer and Nayak [2013] B. Bauer and C. Nayak, Area laws in a many-body localized state and its implications for topological order, Journal of Statistical Mechanics: Theory and Experiment 2013, P09005 (2013).
  * Serbyn  _et al._ [2013b] 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 (2013b).
  * 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).
  * 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).
  * Ros  _et al._ [2015] V. Ros, M. Mueller, and A. Scardicchio, Integrals of motion in the many-body localized phase, Nuclear Physics B 891, 420 (2015).
  * Rademaker and Ortuño [2016] L. Rademaker and M. Ortuño, Explicit local integrals of motion for the many-body localized state, Phys. Rev. Lett. 116, 010404 (2016).
  * Monthus [2016] C. Monthus, Many-body-localization transition: strong multifractality spectrum for matrix elements of local operators, Journal of Statistical Mechanics: Theory and Experiment 2016, 073301 (2016).
  * Thomson and Schiró [2018] S. J. Thomson and M. Schiró, Time evolution of many-body localized systems with the flow equation approach, Phys. Rev. B 97, 060201 (2018).
  * Thomson and Schirò [2023] S. J. Thomson and M. Schirò, Local integrals of motion in quasiperiodic many-body localized systems, SciPost Phys. 14, 125 (2023).
  * O’Brien  _et al._ [2016] T. E. O’Brien, D. A. Abanin, G. Vidal, and Z. Papić, Explicit construction of local conserved operators in disordered many-body systems, Phys. Rev. B 94, 144208 (2016).
  * Peng  _et al._ [2019] P. Peng, Z. Li, H. Yan, K. X. Wei, and P. Cappellaro, Comparing many-body localization lengths via nonperturbative construction of local integrals of motion, Phys. Rev. B 100, 214203 (2019).
  * Leipner-Johns and Wortis [2019] B. Leipner-Johns and R. Wortis, Charge- and spin-specific local integrals of motion in a disordered Hubbard model, Phys. Rev. B 100, 125132 (2019).
  * Adami  _et al._ [2022] S. Adami, M. Amini, and M. Soltani, Structural properties of local integrals of motion across the many-body localization transition via a fast and efficient method for their construction, Phys. Rev. B 106, 054202 (2022).
  * Imbrie [2016] J. Z. Imbrie, Diagonalization and many-body localization for a disordered quantum spin chain, Phys. Rev. Lett. 117, 027201 (2016).
  * Pekker  _et al._ [2014] D. Pekker, G. Refael, E. Altman, E. Demler, and V. Oganesyan, Hilbert-glass transition: New universality of temperature-tuned many-body dynamical quantum criticality, Phys. Rev. X 4, 011052 (2014).
  * Whitlock  _et al._ [2019] S. Whitlock, H. Wildhagen, H. Weimer, and M. Weidemüller, Diffusive to nonergodic dipolar transport in a dissipative atomic medium, Phys. Rev. Lett. 123, 213606 (2019).
  * Signoles  _et al._ [2021] A. Signoles, T. Franz, R. Ferracini Alves, M. Gärttner, S. Whitlock, G. Zürn, and M. Weidemüller, Glassy dynamics in a disordered heisenberg quantum spin system, Phys. Rev. X 11, 011011 (2021).
  * Franz  _et al._ [2022a] T. Franz, S. Geier, C. Hainaut, N. Thaicharoen, A. Braemer, M. Gärttner, G. Zürn, and M. Weidemüller, Observation of universal relaxation dynamics in disordered quantum spin systems, arXiv preprint arXiv:2209.08080 (2022a).
  * Franz  _et al._ [2022b] T. Franz, S. Geier, C. Hainaut, A. Signoles, N. Thaicharoen, A. Tebben, A. Salzinger, A. Braemer, M. Gärttner, G. Zürn, _et al._ , Absence of thermalization in an interacting system of thousands of quantum spins, arXiv preprint arXiv:2207.14216 (2022b).
  * Protopopov  _et al._ [2020] I. V. Protopopov, R. K. Panda, T. Parolini, A. Scardicchio, E. Demler, and D. A. Abanin, Non-abelian symmetries and disorder: A broad nonergodic regime and anomalous thermalization, Phys. Rev. X 10, 011025 (2020).
  * Whitlock  _et al._ [2017] S. Whitlock, A. W. Glaetzle, and P. Hannaford, Simulating quantum spin models using rydberg-excited atomic ensembles in magnetic microtrap arrays, Journal of Physics B: Atomic, Molecular and Optical Physics 50, 074001 (2017).
  * Robert-de Saint-Vincent  _et al._ [2013] M. Robert-de Saint-Vincent, C. S. Hofmann, H. Schempp, G. Günter, S. Whitlock, and M. Weidemüller, Spontaneous avalanche ionization of a strongly blockaded rydberg gas, Phys. Rev. Lett. 110, 045004 (2013).
  * Lukin  _et al._ [2001] M. D. Lukin, M. Fleischhauer, R. Cote, L. M. Duan, D. Jaksch, J. I. Cirac, and P. Zoller, Dipole blockade and quantum information processing in mesoscopic atomic ensembles, Phys. Rev. Lett. 87, 037901 (2001).
  * Braemer  _et al._ [2022] A. Braemer, T. Franz, M. Weidemüller, and M. Gärttner, Pair localization in dipolar systems with tunable positional disorder, Phys. Rev. B 106, 134212 (2022).
  * Rényi [1958] A. Rényi, Egy egydimenziós véletlen térkitöltési problémáról, Publ. Math. Inst. Hung. Acad. Sci 3, 109 (1958).
  * Dasgupta and Ma [1980] C. Dasgupta and S.-k. Ma, Low-temperature properties of the random heisenberg antiferromagnetic chain, Phys. Rev. B 22, 1305 (1980).
  * Fisher [1994] D. S. Fisher, Random antiferromagnetic quantum spin chains, Physical review b 50, 3799 (1994).
  * Vosk and Altman [2013] R. Vosk and E. Altman, Many-body localization in one dimension as a dynamical renormalization group fixed point, Phys. Rev. Lett. 110, 067204 (2013).
  * Turkeshi  _et al._ [2020] X. Turkeshi, P. Ruggiero, and P. Calabrese, Negativity spectrum in the random singlet phase, Phys. Rev. B 101, 064207 (2020).
  * Ruggiero and Turkeshi [2022] P. Ruggiero and X. Turkeshi, Quantum information spreading in random spin chains, Phys. Rev. B 106, 134205 (2022).
  * Vasseur  _et al._ [2016] R. Vasseur, A. J. Friedman, S. A. Parameswaran, and A. C. Potter, Particle-hole symmetry, many-body localization, and topological edge modes, Phys. Rev. B 93, 134207 (2016).
  * Huang and Moore [2014] Y. Huang and J. E. Moore, Excited-state entanglement and thermal mutual information in random spin chains, Phys. Rev. B 90, 220202 (2014).
  * Pouranvari and Yang [2015] M. Pouranvari and K. Yang, Entanglement spectrum and entangled modes of highly excited states in random $xx$ spin chains, Phys. Rev. B 92, 245134 (2015).
  * Note [1] With probability $\approx 0.02$, the illustrated subsystem of states $|\pm\rangle$ contains states of the form $\frac{1}{\sqrt{2}}(|K_{+}\rangle\pm|K_{-}\rangle)$.
  * Luitz [2016] D. J. Luitz, Long tail distributions near the many-body localization transition, Phys. Rev. B 93, 134201 (2016).
  * 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).
  * Note [2] In rare situation of a decimation with $\frac{1}{\sqrt{2}}(|K_{+}\rangle\pm|K_{-}\rangle)$, we consider as operators a combination of two $\hat{P}_{i,j}$.
  * Herviou  _et al._ [2019] L. Herviou, S. Bera, and J. H. Bardarson, Multiscale entanglement clusters at the many-body localization phase transition, Phys. Rev. B 99, 134205 (2019).
  * Szołdra  _et al._ [2021] T. Szołdra, P. Sierant, K. Kottmann, M. Lewenstein, and J. Zakrzewski, Detecting ergodic bubbles at the crossover to many-body localization using neural networks, Phys. Rev. B 104, L140202 (2021).
  * Szołdra  _et al._ [2024] T. Szołdra, P. Sierant, M. Lewenstein, and J. Zakrzewski, Catching thermal avalanches in the disordered xxz model, Phys. Rev. B 109, 134202 (2024).
  * You  _et al._ [2016] Y.-Z. You, X.-L. Qi, and C. Xu, Entanglement holographic mapping of many-body localized system by spectrum bifurcation renormalization group, Phys. Rev. B 93, 104205 (2016).
  * Macé  _et al._ [2019] N. Macé, F. Alet, and N. Laflorencie, Multifractal scalings across the many-body localization transition, Phys. Rev. Lett. 123, 180601 (2019).
  * Vidmar and Rigol [2017] L. Vidmar and M. Rigol, Entanglement entropy of eigenstates of quantum chaotic hamiltonians, Phys. Rev. Lett. 119, 220603 (2017).
  * Vosk  _et al._ [2015] R. Vosk, D. A. Huse, and E. Altman, Theory of the many-body localization transition in one-dimensional systems, Phys. Rev. X 5, 031032 (2015).
  * Potter  _et al._ [2015] A. C. Potter, R. Vasseur, and S. A. Parameswaran, Universal properties of many-body delocalization transitions, Phys. Rev. X 5, 031033 (2015).
  * Goremykina  _et al._ [2019] A. Goremykina, R. Vasseur, and M. Serbyn, Analytically solvable renormalization group for the many-body localization transition, Phys. Rev. Lett. 122, 040601 (2019).
  * Dumitrescu  _et al._ [2019] P. T. Dumitrescu, A. Goremykina, S. A. Parameswaran, M. Serbyn, and R. Vasseur, Kosterlitz-thouless scaling at many-body localization phase transitions, Phys. Rev. B 99, 094205 (2019).
  * Morningstar and Huse [2019] A. Morningstar and D. A. Huse, Renormalization-group study of the many-body localization transition in one dimension, Phys. Rev. B 99, 224205 (2019).
  * Morningstar  _et al._ [2020] A. Morningstar, D. A. Huse, and J. Z. Imbrie, Many-body localization near the critical point, Phys. Rev. B 102, 125134 (2020).
  * 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).
  * Yousefjani and Bayat [2023] R. Yousefjani and A. Bayat, Mobility edge in long-range interacting many-body localized systems, Phys. Rev. B 107, 045108 (2023).
  * Šuntajs  _et al._ [2021] J. Šuntajs, T. Prosen, and L. Vidmar, Spectral properties of three-dimensional anderson model, Annals of Physics 435, 168469 (2021).
  * Šuntajs and Vidmar [2022] J. Šuntajs and L. Vidmar, Ergodicity breaking transition in zero dimensions, Physical Review Letters 129, 060602 (2022).
  * Pawlik  _et al._ [2023] K. Pawlik, P. Sierant, L. Vidmar, and J. Zakrzewski, Many-body mobility edge in quantum sun models (2023), arXiv:2308.01073 [cond-mat.dis-nn] .
  * Haegeman  _et al._ [2011] J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pižorn, H. Verschelde, and F. Verstraete, Time-dependent variational principle for quantum lattices, Phys. Rev. Lett. 107, 070601 (2011).
  * Koffel  _et al._ [2012] T. Koffel, M. Lewenstein, and L. Tagliacozzo, Entanglement entropy for the long-range ising chain in a transverse field, Phys. Rev. Lett. 109, 267203 (2012).
  * Haegeman  _et al._ [2016] J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Unifying time evolution and optimization with matrix product states, Phys. Rev. B 94, 165116 (2016).
  * Paeckel  _et al._ [2019] S. Paeckel, T. Köhler, A. Swoboda, S. R. Manmana, U. Schollwöck, and C. Hubig, Time-evolution methods for matrix-product states, Annals of Physics 411, 167998 (2019).
  * Goto and Danshita [2019] S. Goto and I. Danshita, Performance of the time-dependent variational principle for matrix product states in the long-time evolution of a pure state, Phys. Rev. B 99, 054307 (2019).
  * Chanda  _et al._ [2020b] T. Chanda, P. Sierant, and J. Zakrzewski, Many-body localization transition in large quantum spin chains: The mobility edge, Phys. Rev. Research 2, 032045 (2020b).
  * [123] https://www.scipy.org.

Supplemental Material: Phenomenology of many-body localization in bond-
disordered spin chains

##  S-I Brief Overview

We divide the supplementary material into three broad sections. In the section
S-II, unlike in the letter, we consider a long-range random-bond XXZ model to
make finite-system studies. This section includes details for generating
position configurations, definitions for observables, finite-size scaling of
the observables, and comparing properties with the thermodynamic limit. The
following section S-III gives additional steps for the real space
renormalization group for excited states (RSRG-X) with the calculation of
trivial energy corrections, comparison of the mean gap ratio with the varying
power law exponents, the illustration of thermalization, and the
quantification of the spread of operators associated with the RSRG-X branch.
In the final section S-IV, we provide details of the numerical methods used in
the paper.

##  S-II Properties of long-range random-bond XXZ Hamiltonian

###  S-II.1 Long-range random-bond XXZ model

We consider here the random-bond XXZ Hamiltonian with long-range couplings
with periodic boundary conditions, similar to [86] of the form:

$\hat{\mathcal{H}}_{lr}=\frac{1}{2}\sum_{i\neq
j}J_{ij}(\hat{S}_{x}^{i}\hat{S}_{x}^{j}+\hat{S}_{y}^{i}\hat{S}_{y}^{j}+\Delta\hat{S}_{z}^{i}\hat{S}_{z}^{j}).$
(S.1)

With the same parameters as in the letter, i.e., with $\Delta=-0.73$ and a
power-law decay for both the interaction and hopping strength given as
$J_{ij}=\frac{1}{|x_{i}-x_{j}|^{\alpha}}$, with $\alpha=6$. Again with a
Rydberg blockade diameter $D_{rb}$, for a system with $M$ number of spins and
within the length $L$, defining the spin density $\rho=D_{rb}M/L$. We will
exclusively focus on the zero-magnetization sector, and in this section, we
will only consider the spin inversion symmetry sector $\mathcal{C}=+1$.

###  S-II.2 Generation of position configuration

The disorder can enter long-range random-bond XXZ Hamiltonian due to the
different position configurations that can emerge with decreasing density
($\rho<1$). Here, we describe a method to generate position configurations,
$\\{x_{i}\\}$, while avoiding the problem of Rényi parking constant [87].
Initially, shrink the Rydberg-blockaded hard spheres into points, giving a new
length for the system $\tilde{L}=L-D_{rb}M$. A set of independently drawn
positions is taken from the length $\tilde{L}$, $\\{\tilde{x}_{i}\\}$ to which
$D_{rb}$ padding is added between each $\tilde{x}_{i}$, thus giving the
original length $L$. Henceforth, we shall assume that the diameter of the
Rydberg blockade $D_{rb}$ is unity; thus, the system size $L$ will exclusively
control the density.

Fluctuations of the parameters in space characterize disorders in spin
systems. In $\hat{\mathcal{H}}_{lr}$, the disorder occurs as fluctuations in
the parameter $J_{i,j}$, with both the fluctuation and the mean value around
which the fluctuation occurs dependent on the density. Thus, defining disorder
strength is natural as a fluctuation around the effective $J$ term as

$W=\frac{\sigma_{J^{i}_{eff}}}{\langle J^{i}_{eff}\rangle}$ (S.2)

wherein, $J_{eff}^{i}$ is $\sum_{j}J_{i,j}$ and $\sigma$ represents standard
deviation. The disorder strength increases monotonically with decreasing
density, as shown in Fig. S-1(a). Henceforth, we shall refer to a particular
position configuration as a disorder realization.

Figure S-1: (a) The disorder strength defined as in (S.2) $W$ varying with
inverse density $1/\rho$. With a number of spins $M=16$ and $25000$ different
position configurations. (b) The density of states (DOS)
$\mathcal{P}(\epsilon)$ for the rescaled energy $\epsilon$ for different
system sizes $M\in[12,22]$ for a disorder strength in the transition
($W=0.615$). The inset shows the location of the $\epsilon$ at the maximum
density of states, varying with the number of spins. The dashed orange line
corresponds to the $\epsilon$ value $0.333$.

###  S-II.3 Observables

To depict the crossover between the ergodic and MBL regimes, we consider the
mean gap ratio $\bar{r}$[13] as in the letter and the average rescaled half-
chain entanglement entropy (EE) $S$. The half-chain EE is defined as the von
Neumann entropy of the reduced density matrix $\rho_{M/2}$:

$S_{M/2}=-Tr[\rho_{M/2}\ln\rho_{M/2}]$ (S.3)

where $\rho_{M/2}=Tr_{1,2,\cdots,M/2}|\Psi\rangle\langle\Psi|$ is obtained by
tracing out half of the system. To minimize the dependence of the system size
for EE, we rescale it by the corresponding random matrix theory value
${S}_{RMT}=(M/2)\ln(2)+[1/2+\ln(2)]/2-1/2$ as $S_{\Psi}=S_{M/2}/S_{RMT}$
[105]. Finally, similar to that of $\bar{r}$, we average the rescaled half-
chain EE over the eigenstates corresponding to a particular disorder
realization within an eigenvalue range and then over different disorder
realizations, $S$.

For each disorder strength $W$ and the number of spins $M$, the observables
are calculated by averaging over eigenstates near the rescaled energy

$\epsilon=\frac{E-E_{min}}{E_{max}-E_{min}}=0.333$ (S.4)

for a given disorder realization and subsequently over different disorder
realizations. Note that we do not consider the typical rescaled energy of
$\epsilon=0.5$ to average the observables here, since the density of states
(DOS) peaks around $\epsilon\approx 0.333$ for higher number of spins ($M\geq
14$) for the disorder strength of interest to us, as evident in Fig. S-1(b),
this is a result of the smaller energy gap between $|Z_{\pm}\rangle$ and
$|-\rangle$. The number of eigenvalues and positional configurations used in
this study for different $M$s and the corresponding Hilbert space dimension
$\mathcal{N}$ considered are enumerated in Table 1. The mean gap ratio and the
average rescaled half-chain EE dependence on disorder strength and number of
spins are shown in Fig. S-2.

Table 1: Table showing the approximate number of states considered and the minimum number of disorder realizations used for the statistical observables. $M$ | $\mathcal{N}$ | No. eigenstates | No. realizations
---|---|---|---
$10$ | $126$ | $3$ | $2\cdot 10^{5}$
$12$ | $462$ | $5$ | $10^{6}$
$14$ | $1716$ | $17$ | $10^{5}$
$16$ | $6435$ | $64$ | $8\cdot 10^{3}$
$18$ | $24310$ | $243$ | $2\cdot 10^{3}$
$20$ | $92378$ | $10^{3}$ | $10^{3}$
$22$ | $352716$ | $10^{3}$ | $700$
Figure S-2: The dependence of (a) the mean gap ratio $\bar{r}$, (b) the
rescaled average half-chain EE ${S}$, with varying strengths $W$ for different
number of spins $M\in[10,22]$.

###  S-II.4 Finite-size scaling analysis

In the case of a phase transition, the correlation length $\xi$ diverges near
the critical disorder strength and acts as the dominant length scale. As a
result, the normalized observable, $X$, takes a functional form

$X=\mathcal{G}(L/\xi),$ (S.5)

where $\mathcal{G}(.)$ is a continuous function.

The two prominent proposals for how the correlation length $\xi$ diverges in
the MBL transition is either as a standard second-order phase transition where
$\xi$ diverges as a power-law fashion $\xi_{0}=\frac{1}{|W-W^{*}|^{\nu}}$,
$\nu$ is the critical exponent and $W^{*}$ is the critical disorder strength
based on the early renormalization group approaches [106, 107], or assuming a
Berezinskii-Kosterlitz-Thouless (BKT) like divergence:
$\xi_{B}=\exp\left\\{\frac{b_{\pm}}{\sqrt{|W-W^{*}|}}\right\\}$, where
$b_{\pm}$ are parameters on the two sides of the transition, based on the more
recent real-space renormalization group approaches [108, 109, 110, 111] with
the leading mechanism as the avalanche scenario of delocalization in the MBL
phase.

To best approximate the correlation length parameters, we follow the cost
function minimization method introduced by Suntajs et al. [38, 42] for the XXZ
model in the presence of random disorder. This method has been used for other
MBL models [112, 113]; recent studies have also shown that it predicts the
critical point for the 3-D Anderson [114] and the zero-dimensional quantum sun
model [115, 116] accurately.

Table 2:  Cost function $\mathcal{C}_{X}$ comparison for finite-size scaling
with number of spins dependent drifts in the critical disorder strength.

| $\mathcal{C}_{\bar{r}}[\xi_{0}]$ | $\mathcal{C}_{\bar{r}}[\xi_{B}(b_{+}=b_{-})]$ | $\mathcal{C}_{\bar{r}}[\xi_{B}(b_{+}\neq b_{-})]$
---|---|---|---
$W_{0}$ | $1.1358$ | $1.566$ | $1.226$
$W_{0}+W_{1}M$ | $0.312$ | $0.226$ | $0.222$
$W_{0}+W_{1}\ln(M)$ | $0.304$ | 0.183 | 0.183
$W_{0}+W_{1}/{M}$ | $0.686$ | $0.462$ | $0.462$
$W_{0}+{W_{1}}/{[\ln(M)]}$ | $0.336$ | ${0.199}$ | $0.199$
$W_{0}+W_{1}M^{\gamma}$ | $0.304$ | 0.140 | 0.140
$W(M)$ | $0.282$ | 0.106 | 0.106

| $\mathcal{C}_{S}[\xi_{0}]$ | $\mathcal{C}_{S}[\xi_{B}(b_{+}=b_{-})]$ | $\mathcal{C}_{S}[\xi_{B}(b_{+}\neq b_{-})]$
---|---|---|---
$W_{0}$ | $1.117$ | $1.338$ | $1.211$
$W_{0}+W_{1}M$ | $0.122$ | $0.0.279$ | $0.168$
$W_{0}+W_{1}\ln(M)$ | $0.102$ | $0.131$ | $0.098$
$W_{0}+W_{1}/{M}$ | $0.397$ | $0.542$ | $0.499$
$W_{0}+{W_{1}}/{[\ln(M)]}$ | $0.109$ | 0.0504 | 0.0504
$W_{0}+W_{1}M^{\gamma}$ | $0.102$ | 0.0507 | 0.0507
$W(M)$ | $0.086$ | 0.030 | 0.030

Suppose that the observables of interest are continuous monotonic functions of
$sgn[W-W^{*}]L/\xi$. The cost function for a quantity $X=\\{X_{j}\\}$ that
consists of $N$ different $W$ and $L$ is defined as

$\mathcal{C}_{X}=\frac{\sum_{j=1}^{N-1}|X_{j+1}-X_{j}|}{\max\\{X_{j}\\}-\min\\{X_{j}\\}}$
(S.6)

where $X_{j}$’s are sorted according to non-decreasing values of
$sgn[W-W^{*}]L/\xi$. For an ideal collapse with $X$ being a monotonic function
of $sgn[W-W^{*}]L/\xi$, we must have
$\sum_{j}|X_{j+1}-X_{j}|=\max\\{X_{j}\\}-\min\\{X_{j}\\}$, and thus
$\mathcal{C}_{X}=0$. Thus, the best collapse corresponds to the global minima
of $\mathcal{C}_{X}$ for different correlation lengths and functional forms of
$W^{*}(L)$. The final section details the numerical parameters for the cost
function minimization performed.

To perform an unbiased finite size scaling analysis for
$\hat{\mathcal{H}}_{lr}$, we consider observables of system sizes
$M\in[12,22]$ within disorder strengths $W\in[0.3,0.725]$, minimizing the cost
function for the correlation lengths $\xi_{0},\xi_{B}$ with different
functional forms of the critical disorder, $W^{*}$ on the number of spins
($M$) as shown in Table 2.

With a fixed critical disorder strength $W^{*}=W_{0}$, the power law
correlation length furnishes a smaller value $\mathcal{C}_{X}$ for both the
mean gap ratio, $\bar{r}$, and the average rescaled half-chain EE, $S$.
Finding critical disorder strengths as $W^{*}_{r}=0.584$ and $W^{*}_{S}=0.607$
and critical exponents $\nu_{\bar{r}}=0.661$ and $\nu_{S}=0.662$.

While assuming different functional forms for the critical disorder strength,
the best-minimized cost function value is a BKT-like correlation with
sublinear scaling, Table 2; a functional form $W^{*}=W_{0}+W_{1}L^{\gamma}$
has the best cost function value for both $\bar{r}$ and $S$. The minimized
values for the drift $(W_{0},\gamma)$ are quite different for
$\bar{r}\rightarrow(1.36,-0.32)$ and $S\rightarrow(1.00,-0.59)$. With the
nonuniversal parameter, $b$, for $\bar{r}$ as $0.804$ and for $S$ as $1.09$.
There appears to be no improvement in $\mathcal{C}_{X}$ when considering
different nonuniversal parameters $b_{+}\neq b_{-}$, as the observables
considered appear mostly in the localized regime. Considering a generic
function for $W^{*}=W(M)$, where each $W(M)$ is chosen individually for $M$,
we find the functional form $W^{*}=W_{0}+W_{1}L^{\gamma}$ to approximate
$W(M)$ as shown in Fig. S-3.

These results, though, do not indicate a strong interpretation of $W_{0}$ in
the functional form $W^{*}=W_{0}+W_{1}L^{\gamma}$, as the critical disorder
strength at the thermodynamic limit as in the 3-D Anderson model [114], this
is because the values found $W^{*}$ for the considered number of spins $M$ are
pretty far from $W_{0}$. However, $2\gamma_{\bar{r}}\sim\gamma_{S}$ indicates
that $S$ might more evidently show properties of the thermodynamic limit in
the $M$ considered than in $\bar{r}$.

Figure S-3: Finite-size scaling for mean gap-ratio ($\bar{r}$) and rescaled
averaged half-chain EE ($S$), with correlation length $\xi_{B}$ and the
critical disorder strength with a drift of functional form
$W^{*}=W_{0}+W_{1}M^{\gamma}$ is shown in (a) and (b). While (c) and (d) show
the dependency of free disorder strength $W^{*}(M)$, on $M$ in comparison with
the functional form $W^{*}=W_{0}+W_{1}M^{\gamma}$.

###  S-II.5 How well estimated are the finite-spin observables compared to the
thermodynamic limit?

To address the question of how well estimated the finite-spin observables are
to properties in the thermodynamic limit, we define two disorder strengths as
in [45, 32]

  1. 1.

$W_{T}^{X}(M)$: the disorder strength at which, for a given number of spins
$M$, the quantity $X$ deviates from a percentage of its ergodic value compared
to the localized non-ergodic value.

  2. 2.

$W_{\times}^{X}(M)$: the disorder strength at which the curves $X(W)$ cross
for systems with the number of spins $M-1$ and $M+1$.

The disorder strength $W_{T}^{X}(M)$ can be considered as the boundary of the
ergodic regime. In contrast, the disorder strength $W_{\times}^{X}(M)$ when
$M\rightarrow\infty$ approximates the critical disorder strength assuming MBL
persists in the thermodynamic limit, $W_{\infty}$. Depending on the type of
transition, the regime between $W_{T}^{X}(M)$ and $W_{\times}^{X}(M)$ either
vanishes for $L\rightarrow\infty$ or tends to a constant.

For the random disordered XXZ model, the disorder strength $W_{T}^{X}(M)$
increases monotonously with the number of spins, $W_{T}^{X}(M)\sim M$ and
$W_{\times}^{X}(M)\sim W_{\infty}-a_{1}/L$ [34]. However, these scalings are
incompatible with the existence of MBL in the large length limit when
$W_{T}^{X}(M)$ exceeds $W_{\times}^{X}(M)$ at $M_{0}^{XXZ}\approx 50$. Thus,
approaching the number of spins $M_{0}^{XXZ}$ indicates the presence or
absence of MBL and can be extended to compare how well-approximated
observables are with the thermodynamic limit.

In case of $\hat{\mathcal{H}}_{lr}$, the results of which are shown in Fig.
S-4, with $W_{\times}^{X}(M)$ we consider a functional form
$W_{\times}^{X}(M)=W_{\infty}+a_{1}/M+a_{2}/M^{2}$ from which $W_{\infty}$ is
extracted for $\bar{r}$ and $S$, respectively, as
$W_{\infty}^{\bar{r}}=0.951\pm 0.157$ and $W_{\infty}^{S}=0.764\pm 0.134$. A
substantial percentage difference is considered from the ergodic value to
estimate $W_{T}^{X}(M)$ because, for $M\leq 14$, $X$ does not saturate to the
ergodic value even with small disorder strengths, $8.5\%$ for $\bar{r}$ and
$20\%$ for $S$. Similar to the Floquet model considered in [32] $M\geq 20$,
$W_{T}^{X}(M)$ deviates from a linear function, so we make a linear fit with
$M\geq 20$ to make a comparison with the XXZ model. From the curves extracted
from $W_{T}^{X}(M)$ and $W_{\times}^{X}(M)$ we find that the value $M$ at
which the two curves cross as $M_{0}^{\bar{r}}\in[34.5,43.1]$ and
$M_{0}^{S}\in[31.9,38.1]$ is less than the XXZ model with random disorder [34]
hinting at better approximating thermodynamic properties, but higher than the
Floquet model in [32].

Figure S-4: (a) The crossing point disorder strength $W_{\times}^{\bar{r}}$
and the near ergodic disorder strength $W_{T}^{\bar{r}}$ varying with $M$ for
$\bar{r}$. The bold blue line denotes a fit on $W_{\times}^{\bar{r}}$ with a
functional form $W_{\infty}+a_{1}/M+a_{2}/M^{2}$ the grayed area denotes the
associated errors, while the dashed red line denotes a linear fit for
$W_{T}^{\bar{r}}$. In (b), the same results as in (a) are shown for $S$.

##  S-III Extension of RSRG-X results

###  S-III.1 Random-bond XXZ model

For this section, we define the random-bond XXZ model with nearest-neighbor
interaction as in the letter again for ease of reading.

$\hat{\mathcal{H}}=\sum_{i=1}^{M-1}J_{i,i+1}(\hat{s}_{i}^{x}\hat{s}_{i+1}^{x}+\hat{s}_{i}^{y}\hat{s}_{i+1}^{y}+\Delta\hat{s}_{i}^{z}\hat{s}_{i+1}^{z}).$
(S.7)

The parameters are the same as in (S.7). The disorder strength is defined
$W\equiv 1/\rho-1$ as in the letter.

###  S-III.2 Second order perturbation for $|Z_{\pm}\rangle$ states.

Fixing spins to a state $|Z_{\pm}\rangle$ at a decimation leads to effective
nearby operators modified by a first-order perturbation to the chemical
potential form ${\hat{\mathcal{H}}^{\prime}}_{i,j}$ thereby typically leading
to a process of growing domain spins. There is an exception when the
generation of domain spins leads to degenerate states of the form
$|K_{+}\rangle\equiv|\uparrow\rangle|Z_{-}\rangle|\uparrow\rangle$,
$|K_{-}\rangle\equiv|\downarrow\rangle|Z_{+}\rangle|\downarrow\rangle$. The
states $|K_{\pm}\rangle$ could be coupled through second-order perturbation if
the resulting energy gap is greater than the effective neighbouring operators
${\hat{\mathcal{H}}^{\prime}}_{i,j}$ to $|K_{\pm}\rangle$. Eventually, this
coupling would result in operators of the form $\hat{\mathcal{H}}_{i,j}$.

Figure S-5: (a) Show the absolute value of energy difference between the
approximate states of $|\phi_{+}\rangle$ and $|\phi_{-}\rangle$, given as
$|E_{+}-E_{-}|$ represented with bold line, while dashed lines represents
$\tilde{J}_{1,6}^{+}$ from (S.9). (b) The probability of spread of operators,
$p_{l}$, rescaled with $p_{1}$, where $l$ represents the distance between
sites. The rescaled values shown in the above figure were calculated for $10$
different decimations $\gamma$s, with $10^{4}$ positional configurations for
$M=100$ with $\mathcal{C}=+1$ .

Consider the strongest bond is between the spins $i=3$ and $i=4$. Fixing the
spins with the choice $|Z_{\pm}\rangle$ results in operators of the form
$\hat{\mathcal{H}}^{\prime}_{4,5}$ , $\hat{\mathcal{H}}^{\prime}_{2,3}$ these
operators in further decimation steps can give states $|K_{\pm}\rangle$ and
generation of operators $\hat{\mathcal{H}}^{\prime}_{5,6}$ and
$\hat{\mathcal{H}}^{\prime}_{1,2}$. In such a scenario, in addition to the
operators $\hat{\mathcal{H}}^{\prime}_{5,6}$,
$\hat{\mathcal{H}}^{\prime}_{1,2}$, a second-order perturbation couples the
states $|K_{\pm}\rangle$. To elucidate this coupling better, assume that the
spins $i=3$ , $i=4$ are instead fixed with the states
$\frac{1}{\sqrt{2}}(|Z_{+}\rangle\pm|Z_{-}\rangle)$ which would result in an
effective operator of the form

$\displaystyle\hat{\mathcal{H}}^{\prime\prime}_{2,5}$ $\displaystyle=$
$\displaystyle\frac{J^{\prime}}{2}\sum_{\Theta=\pm
1}\biggl{[}\frac{\pm\Theta}{(1+\Theta\Delta_{3,4})}(\hat{s}^{+}_{2}\hat{s}^{+}_{5}+h.c.)$
(S.8)
$\displaystyle-\frac{1}{(1+\Theta\Delta_{3,4})}(\hat{s}_{2}^{+}\hat{s}_{5}^{-}+h.c.)\biggl{]}$

where, $J^{\prime}$ is defined as
$J^{\prime}\equiv\frac{J_{2,3}J_{4,5}}{2J_{3,4}}$. The second term in (S.8) is
just the operator considered with second-order perturbation with choice of
states $|Z_{\pm}\rangle$. The first term, which vanishes when
$|\Delta|\rightarrow 0$, results in states
$\frac{1}{\sqrt{2}}(|Z_{+}\rangle\pm|Z_{-}\rangle)$ on spins $i=2$, $i=5$.
Combined with restrictions due to $\hat{\mathcal{H}}^{\prime}_{4,5}$ and
$\hat{\mathcal{H}}^{\prime}_{2,3}$ and local $U(1)$ symmetry, leads to states
of the form $\frac{1}{\sqrt{2}}(|K_{+}\rangle\pm|K_{-}\rangle)$ . If the
energy gap resulting from the effective neighbour operators
$\hat{\mathcal{H}}^{\prime}_{5,6}$, $\hat{\mathcal{H}}^{\prime}_{1,2}$, are
smaller than $\hat{\mathcal{H}}^{\prime\prime}_{2,5}$ we would make a
decimation fixing to states
$\frac{1}{\sqrt{2}}(|K_{+}\rangle\pm|K_{-}\rangle)$, which would result in
modifying the nearby operators via second-order perturbation to take the form
$\hat{\mathcal{H}}_{1,6}$ with parameters

$\tilde{J}_{1,6}^{\pm}=\sum_{\Theta=\pm
1}\frac{J_{1,2}J_{5,6}}{2J^{\prime}\frac{1}{1\mp\Theta\Delta_{3,4}}+\Theta(J_{2,3}\Delta_{2,3}+J_{4,5}\Delta_{4,5})}$
(S.9)

and

$\tilde{J}_{1,6}^{\pm}\tilde{\Delta}_{1,6}^{\pm}=\frac{1}{4}\left(\frac{J_{1,2}\Delta_{1,2}J_{5,6}\Delta_{5,6}}{\pm
J^{\prime}(\frac{1}{1+\Delta_{3,4}}-\frac{1}{1-\Delta_{3,4}})}\right).$ (S.10)

To check the validity of (S.9) in Fig. S-5(a) we consider a particular
disorder realization that satisfies the above condition for a system with
$M=6$ and make a comparison between $\tilde{J}^{+}_{1,6}$ and the absolute
energy difference $|E_{+}-E_{-}|$ between two states with maximum overlap to
the states
$|\phi_{\pm}\rangle\equiv\frac{1}{\sqrt{2}}(|K_{+}\rangle+|K_{-}\rangle)|\pm\rangle_{1,6}$.
We find that with varying $\Delta$, the absolute energy difference is well
approximated by $\tilde{J}^{+}_{1,6}$ except when $|\Delta|\rightarrow 1$.

###  S-III.3 Trivial energy corrections.

At each decimation step, one encounters trivial energy corrections $\delta E$,
which are scalar values that modify the energy. These corrections are
essential for the precise calculation of energy; in Table 3, we show the
trivial energy correction associated with each decimation step. Where,
$J^{\prime}$ is again defined as
$J^{\prime}\equiv\frac{J_{2,3}J_{4,5}}{2J_{3,4}}$.

Table 3: Trivial energy $\delta E$ corrections with second order perturbation for each choice of state. State | $\delta E$
---|---
$|\pm\rangle$ | $\mp\frac{J_{2,3}^{2}+J_{4,5}^{2}}{4J_{3,4}}\biggl{(}\frac{\Delta_{3,4}^{2}}{4}+\frac{1}{1\pm\Delta_{3,4}}\biggl{)}$
$\frac{1}{\sqrt{2}}(|Z_{+}\rangle\pm|Z_{-}\rangle)$ | $\frac{\left(J_{2,3}^{2}+J_{4,5}^{2}\right)}{8J_{3,4}}\left(\frac{1}{(1+\Delta_{3,4})}-\frac{1}{(1-\Delta_{3,4})}\right)$
$|K_{\pm}\rangle$ | $\frac{J_{2,3}^{2}+J_{4,5}^{2}}{-4\left(J_{2,3}\Delta_{2,3}+J_{4,5}\Delta_{4,5}\right)+8J^{\prime}\left(\frac{1}{1\pm\Delta_{2,3}}\right)}$
|
$-\frac{J_{2,3}^{2}+J_{4,5}^{2}}{4\left(J_{2,3}\Delta_{2,3}+J_{4,5}\Delta_{4,5}\right)+8J^{\prime}\left(\frac{1}{1\mp\Delta_{2,3}}\right)}$
|
$\pm\frac{J_{2,3}^{2}\Delta_{2,3}^{2}+J_{4,5}^{2}\Delta_{4,5}^{2}}{8J^{\prime}\left(\frac{1}{1+\Delta_{2,3}}-\frac{1}{1-\Delta_{2,3}}\right)}$

###  S-III.4 Mean gap ratio with broader distribution

On considering exponents for the power law controlling $J$’s, $\alpha$, we
find that RSRG-X well approximates the mean gap ratio $\bar{r}$ compared to
the exact diagonalization results, Fig. S-6 (a). More significantly, the
distribution of the gap ratio averaged over a disorder realization
$P_{\delta_{r_{s}}}$ decreases with $\alpha$ strongly Fig. S-6 (b).

Figure S-6: (a) The mean gap ratio ($\bar{r}$) obtained by RSRG-X and exact
diagonalization (ED) varying with the exponent $\alpha$. The inset shows the
relative difference $\delta_{\bar{r}}=\bar{r}_{ED}-\bar{r}_{RG}$. (b)
Distribution of relative differences in the gap ratios averaged over a
disorder realization ($P_{\delta{r_{s}}}$) with the exponent $\alpha$. The
inset shows the variation of the standard deviations of $P_{\delta{r_{s}}}$
with respect to $\alpha$, given as $\sigma_{\delta_{r_{s}}}$. The system size
considered is $M=10$ in the $\mathcal{C}=+1$ sector with fixed disorder
strength $W=8.0$, with $25000$ disorder realizations and considering all the
eigenvalues

###  S-III.5 Illutration of thermalization

Figure S-7: The deviation of the absolute value of long-time averaged
operators $|\bar{O_{i}}|$ at each site with varying densities. (a) $M=30$, (b)
$M=16$ for the same positional configuration and set of operators as in Fig.
(3) in the letter. In both cases averages are taken up to time $t=5000$.

In [66], infinite time averages of local operators are taken, giving the
correlation length $\zeta$ associated with the local integral of motions
(LIOMs). Here, we look for the absolute value of infinite time averages of the
operators corresponding to an RSRG-X branch
$|\bar{\mathcal{O}}_{n}|=|\int_{0}^{T}\mathcal{O}_{n}(t)dt|$, an unchanging
value will indicate that RSRG-X branch is stable. A variation from the
stability that does not spread through the system, which typically arises from
ergodic bubbles, indicates that $\hat{\mathcal{O}}_{n}$ are not the best form
of local operators in the given RSRG-X branch. However, the system as a whole
is still localized. While the operators’ instability throughout the system
would indicate thermalization. To visualize these operators spatially, in Fig.
S-7 we consider that the values $|\bar{\mathcal{O}}_{n}|$ are distributed
equally in space,
$|\bar{\mathcal{O}}_{i}|=\frac{1}{2}|(\bar{\mathcal{O}}_{n}^{i,j})|$. Here
$i,j$ denotes indices associated of spins. Fig. S-7(a) considers a position
configuration with $M=30$ for a strong disorder with a particular disorder
realization and RSRG-X branch. In contrast, in Fig. S-7(b), we consider a
smaller system size with the same position configuration and operators as in
Fig. (3) of the letter to explore weaker disorder regimes. In Fig. S-7(a), we
find that at strong disorder, most of the operators are preserved, while with
decreasing disorder strength, new ergodic bubbles pop up as seen in the region
between sites $12-15$ but stay strictly local in most cases except between
sites $3-6$ at disorder strength $W=1$. While in Fig. S-7(b), we find as we
approach the region of weak disorder $W\sim 0.1$, the ergodic bubble between
sites $3-6$ starts to thermalize the nearby region, eventually thermalizing
operators at all sites indicating ergodicity.

###  S-III.6 Localization length

To contrast the spread of operators associated with an RSRG-X branch
$\hat{\mathcal{O}}_{n,k}^{i_{a},i_{b}}$ with the correlation length of LIOMs
with on-site disodered systems [64, 66, 67] we look for the spread of
operators on the spins. Here, $i_{a}$ denotes the spin indices and $k$ denotes
the associated eigen index, by defining the distance between the spins of
these operators as $l=|i_{a}-i_{b}|$, for a single site operator $l=0$, we
define a probability of the spread as

$p_{l}=\frac{1}{\eta}\sum_{n,k}||\hat{\mathcal{O}}_{n,k}^{l}||$ (S.11)

Where $||\hat{\mathcal{O}}_{n,k}^{l}||$ denotes the Frobenius norm of the
operator and $\eta$ is a normalization constant such that $\sum_{l}p_{l}=1$.
If the probability distributions are exponentially distributed,
$p_{l}\propto\exp{(-l/\zeta_{0})}$, where $\zeta_{0}$ is the length of the
localization at asymptotic disorder strength, we would identify our system as
localized.

In Fig. S-5(b), using RSRG-X and assuming the Frobenius norm associated with
the operator of an RSRG-X branch as unity, we compare $p_{l}$ for a finite
$\Delta$ with that of $\Delta=0$. For $\Delta=-0.73$, we find that $p_{l}$
decreases as a stretched exponential with the distance between sites. In
contrast, it decreases as a power law in the non-interacting case. Thus, we
conclude that the interaction case is quasilocal in the asysmptotic disorder
strength. The other extreme of relevance is $|\Delta|\gg 1$ which will have
$p_{l\neq 0}=0$.

##  S-IV  Description of numerics

###  S-IV.1 Diagonalization

Statistical observables, namely the mean gap ratio, $\bar{r}$, and the
rescaled average half-chain EE, $S$, required the extraction of eigenenergies
and eigenstates of the Hamiltonian in Eq. S.1 at high energies
($\epsilon=0.333$). Note that while extracting high energy states in the main
letter as in Fig. 2(a) and Fig. 4(b) we consider again states around
$\epsilon=0.333$ with half the disorder realization as in Table 1. For a
system with a small number of spins $M\leq 16$, we use the standard exact
diagonalization (ED) method for dense matrices, while for $M\geq 18$ the
recently developed state-of-the-art method for sparse matrices polynomially
filtered exact diagonalization (POLFED) [34] was used.

###  S-IV.2 Time evolution

For the time evolution of a system with a small number of spins, i.e., $M\leq
16$, we employed exact diagonalization to make a time evolution. In contrast,
for larger particle numbers, we employ the matrix product state (MPS) based
time-dependent variational principle (TDVP) method [117, 118, 119].
Specifically, we use the ‘hybrid’ scheme of TDVP [120, 121, 37, 122] with step
size $\delta t=0.1$, where the two-site variant of TDVP is initially used to
dynamically grow the bond dimension of the MPS until the dimension of the
largest MPS bond reaches a predetermined value $\chi_{max}$. For subsequent
dynamics, a one-site variant of TDVP is employed. In our calculations, we set
the maximum value the bond dimension can grow to as $\chi_{max}=1024$, while
the bond dimensions do not reach the value $\chi_{max}$ for the maximum time
up to which we make the time evolution $t=5000$.

###  S-IV.3 Parameters for cost-function minimization

To minimize the cost function, we employ the differential evolution method
implemented in SciPy [123]. We run a series of independent differential
evolution algorithms with different random seeds for each $\mathcal{C}_{X}$
shown in 2. We employ $1000$ independent realizations to find the optimal
minimization. In each realization, we allow up to $10^{4}$ iterations with a
relative tolerance of convergence $10^{-4}$ and a population size of $250$.

###  S-IV.4 Calculation of errors

Figure S-8: An example of extracting the crossing point $W_{\times}$ between
systems with two different numbers of spins $M$ for both (a) the mean gap
ratio $\bar{r}$ and (b) the average rescaled half-chain EE $S$. The black
dashed line denotes the identified crossing point, while the grayed part
signifies the error margins.

After averaging over different disorder realizations, the statistical errors
associated with the mean value of observables $\bar{X}$ are calculated by
initially averaging over the disorder realization, $X_{s}$, then identifying
the uncertainty [32]

$\sigma_{X}=\frac{(\langle(X_{s}-\bar{X})^{2}\rangle)^{1/2}}{N_{dis}^{1/2}}$
(S.12)

here $N_{dis}$ is the number of disorder realizations.

To calculate the errors associated with the disorder strengths $W_{T}^{X}(M)$
and $W_{\times}^{X}(M)$ and the errors associated with fitting them to a
functional form, we consider the random sampling of errors that are assumed to
be Gaussian distributed. For example, to calculate the crossing point disorder
strength $W_{\times}^{X}(M)$ and the associated errors. Resample from a
Gaussian distribution with the standard deviation given as $\sigma_{X}$ around
a mean value $\bar{X}$ and fit the associated curves of the samples of system
sizes $M-1$ and $M+1$ in the transition regime with a $5$th order polynomial
and extract the point where they cross. Re-sampling over a $10000$ times, we
find a distribution of intercepts, the mean value of which gives the
associated crossing point, while the standard deviation gives the associated
error as shown in Fig. S-8. A similar procedure calculates the error
associated with $W_{T}^{X}(M)$ with a cubic spline instead of a polynomial
fit. While the errors associated with the fits, shown as the gray area in S-4,
are extracted by resampling and fitting with the functional form
$W_{\times}^{X}(M)=W_{\infty}+a_{1}/M+a_{2}/M^{2}$.


# A note on powerful numbers in short intervals

Tsz Ho Chan

###### Abstract

In this note, we are interested in obtaining uniform upper bounds for the
number of powerful numbers in short intervals $(x,x+y]$. We obtain
unconditional upper bounds $O(\frac{y}{\log y})$ and $O(y^{11/12})$ for all
powerful numbers and $y^{1/2}$-smooth powerful numbers respectively.
Conditional on the $abc$-conjecture, we prove the bound
$O(\frac{y}{\log^{1+\epsilon}y})$ for squarefull numbers and the bound
$O(y^{(2+\epsilon)/k})$ for $k$-full numbers when $k\geq 3$. They are related
to Roth’s theorem on arithmetic progressions and the conjecture on non-
existence of three consecutive squarefull numbers.

##  1 Introduction and Main Result

A number $n$ is squarefull if its prime factorization
$n=p_{1}^{a_{1}}p_{2}^{a_{2}}\cdots p_{r}^{a_{r}}$ satisfies $a_{i}\geq 2$ for
all $1\leq i\leq r$. Similarly, a number $n$ is $k$-full if the exponents
$a_{i}\geq k$ for all $1\leq i\leq r$. For example, $72=2^{3}\cdot 3^{2}$ is
squarefull and $243=3^{5}$ is $5$-full. Let $Q_{k}(x)$ denote the number of
$k$-full numbers which are less than or equal to $x$. It is known that

$Q_{k}(x)=\prod_{p}\Bigl{(}1+\sum_{m=k+1}^{2k-1}\frac{1}{p^{m/h}}\Bigr{)}x^{1/k}+O(x^{1/(k+1)})$
(1)

where the product is over all primes (see [4] or [1] for example.) There are
also interests in counting the number of $k$-full numbers in short intervals
$(x,x+y]$ with $y=o(x)$. For moderate size $y$, there are some asymptotic
results: For example, Trifonov [8] and Liu [6] obtained,

$Q_{2}\bigl{(}x+x^{1/2+\theta}\bigr{)}-Q_{2}(x)\sim\frac{\zeta(3/2)}{2\zeta(3)}x^{\theta}\;\;\text{
for }\;\;19/154=0.12337...<\theta<1/2,$ (2)

and

$Q_{3}\bigl{(}x+x^{2/3+\theta}\bigr{)}-Q_{3}(x)\sim\frac{\zeta(4/3)}{3\zeta(4)}x^{\theta}\;\;\text{
for }\;\;5/42=0.11904...<\theta<1/3$ (3)

respectively.

What happens when $y$ is very small such as $y\ll x^{1/2}$ or even $y\ll\log
x$? For such short intervals, one can only expect suitable upper bounds rather
than asymptotic formulas. Thus, in this note, we are interested in finding
uniform upper bounds for $Q_{k}(x+y)-Q_{k}(x)$ with $1\leq y\leq x$ that are
independent of $x$. By comparing $k$-full numbers with perfect $k$-th powers,
we suspect the following to be true.

###### Conjecture 1

Given an integer $k\geq 2$ and a real number $x\geq 1$. There exists some
constant $C_{k}\geq 1$ such that

$Q_{k}(x+y)-Q_{k}(x)\leq C_{k}y^{1/k}$

uniformly over $1\leq y\leq x$.

We are far from proving this at the moment. The current best upper bound

$Q_{k}(x+y)-Q_{k}(x)\ll\frac{y\log\log(y+2)}{\log(y+2)}$ (4)

was obtained by De Koninck, Luca and Shparlinski [3]. In this note, we improve
(4) slightly.

###### Theorem 1

Given an integer $k\geq 2$ and a real number $x\geq 1$. We have

$Q_{k}(x+y)-Q_{k}(x)\ll\frac{y}{\log(y+1)}$ (5)

uniformly over $1\leq y\leq x$.

In fact, we shall prove the following more general result concerning
squarefull numbers in arithmetic progressions over short intervals which gives
Theorem 1 immediately as $k$-full numbers are included in squarefull numbers.

###### Theorem 2

Given real numbers $x\geq 1$ and $0<\alpha<1$. Suppose $q>0$ and $r$ are two
integers with $\text{gcd}(r,q)=1$. We have

$\mathop{\mathop{\sum_{x<n\leq x+y}}_{n\text{ squarefull}}}_{n\equiv
r\;(\bmod{q})}1\ll_{\alpha}\frac{y}{\phi(q)\log(y+1)}$

uniformly over $1\leq y\leq x$ and $1\leq q\leq y^{1-\alpha}$.

Using similar technique, one can obtain some power savings over (5) for smooth
$k$-full numbers in short intervals.

###### Theorem 3

Given an integer $k\geq 2$ and a real number $x\geq 1$. We have

$\mathop{\mathop{\sum_{x<n\leq x+y}}_{n\text{ $k$-full}}}_{p^{+}(n)\leq
y^{1/2}}1\leq\mathop{\mathop{\sum_{x<n\leq x+y}}_{n\text{
squarefull}}}_{p^{+}(n)\leq y^{1/2}}1\ll y^{11/12}.$ (6)

uniformly over $1\leq y\leq x$. Here $p^{+}(n)$ stands for the largest prime
factor of $n$. Remark: One may increase the exponent $1/2$ up to $1$ and
obtain similar power saving upper bound.

The bound (6) lends evidence towards Conjecture 1 and shows that the
difficulty lies with non-smooth $k$-full numbers. Another piece of evidence
comes from the famous $abc$-conjecture. It was proved in [3] that, given any
$\epsilon>0$, the interval

$(x,x+x^{1-(2+\epsilon)/k}]$ (7)

contains at most one $k$-full number for sufficiently large $x$ under the
$abc$-conjecture. From this, one has

###### Theorem 4

Given an integer $k\geq 3$ and real numbers $\epsilon>0$ and $x\geq 1$. We
have

$Q_{k}(x+y)-Q_{k}(x)\ll_{\epsilon}y^{(2+\epsilon)/k}$ (8)

uniformly over $1\leq y\leq x$.

Unfortunately, the proof in [3] concerning (7) was a little inaccurate as the
$a$, $b$, $c$ in its application of the $abc$-conjecture might not be
relatively prime. So, we shall give a complete proof of Theorem 4 later.
Furthermore, (7) gives us nothing when $k=2$. To remedy this, we shall prove
the following conditional result which improves on (5) slightly by a small
power of logarithm.

###### Theorem 5

The $abc$-conjecture implies that, for some absolute constant $c>0$,

$Q_{2}(x+y)-Q_{2}(x)\ll\frac{y}{\log^{1+c}(y+1)}$

uniformly over $1\leq y\leq x$.

Its proof relies on the following recent breakthrough result of Bloom and
Sisask [2] on density of integer sequence without three term arithmetic
progressions.

######  Theorem 6 (Bloom-Sisask)

Let $N\geq 2$ and $A\subset\\{1,2,...,N\\}$ be a set with no non-trivial
three-term arithmetic progressions, i.e. solutions to $x+y=2z$ with $x\neq y$.
Then

$|A|\ll\frac{N}{(\log N)^{1+c}},$

where $c>0$ is an absolute constant.

This paper is organized as follows. First, we will prove Theorems 2 and 3
using Brun-Titchmarsh inequality and ideas from Shiu’s generalization [7].
Then we will prove Theorem 4 using the $abc$-conjecture. Finally, we will
prove Theorem 5 by establishing non-existence of three term arithmetic
progressions for squarefull numbers in short intervals.

Notation. We use $|A|$ to denote the number of elements in a finite set $A$
and $\lfloor x\rfloor$ to denote the greatest integer less than or equal to
$x$. We let $p_{-}(n)$ and $p^{+}(n)$ be the smallest and the largest prime
factor of $n$ respectively. The symbols $f(x)=O(g(x))$, $f(x)\ll g(x)$ and
$g(x)\gg f(x)$ are equivalent to $|f(x)|\leq Cg(x)$ for some constant $C>0$.
Also, $f(x)=O_{\lambda}(g(x))$, $f(x)\ll_{\lambda}g(x)$ or
$g(x)\gg_{\lambda}f(x)$ mean that the implicit constant may depend on
$\lambda$. Furthermore, $f(x)=o(g(x))$ means
$\lim_{x\rightarrow\infty}f(x)/g(x)=0$ and $f(x)\sim g(x)$ means
$\lim_{x\rightarrow\infty}f(x)/g(x)=1$. Finally, the summation symbol
$\mathop{\sum\nolimits^{\prime}}$ signifies that a sum is over squarefull
numbers only.

##  2 Some Preparations

###### Lemma 1

For any $X\geq 1$,

$\mathop{\sum\nolimits^{\prime}}_{X<n\leq X^{2}}\frac{1}{n}\ll X^{-1/2}.$

Proof: From (1), we have $Q_{2}(X)\ll X^{1/2}$. By partial summation, the
above sum is

$=\int_{X}^{X^{2}}\frac{1}{u}dQ(u)=\frac{Q(X^{2})}{X^{2}}-\frac{Q(X)}{X}+\int_{X}^{X^{2}}\frac{Q(u)}{u^{2}}du\ll\frac{1}{X^{1/2}}+\int_{X}^{X^{2}}\frac{1}{u^{3/2}}du\ll\frac{1}{X^{1/2}}.$

######  Lemma 2 (Brun-Titchmarsh inequality)

Let $q\geq 1$ and $r$ be integers satisfying $\text{gcd}(r,q)=1$. Suppose
$q<y\leq x$ and $z\geq 2$. Then

$\mathop{\mathop{\sum_{x<n\leq x+y}}_{n\equiv
r\;(\bmod{q})}}_{p_{-}(n)>z}1\ll\frac{y}{\phi(q)\log z}+z^{2}.$

Remark: The above bound is still true when $y\leq q$ or $y<1$ since there is
at most one term in the sum.

Proof: This follows from Selberg upper bound sieve method (see [5], page 104
for example).

Finally, let us recall the $abc$-conjecture. For any nonzero integer $m$, let

$\kappa(m):=\prod_{p|m}p$

be the kernel of $m$.

######  Conjecture 2 ($abc$-conjecture)

For any $\epsilon>0$, there exists a constant $C_{\epsilon}>0$ such that, for
any integers $a,b,c$ with $a+b=c$ and $\text{gcd}(a,b)=1$, the bound

$\max\\{|a|,|b|,|c|\\}\leq C_{\epsilon}\kappa(abc)^{1+\epsilon}$

holds.

##  3 Proof of Theorem 2

Our proof is inspired by Shiu [7] on Brun-Titchmarsh theorem for
multiplicative functions. We may assume that $y\geq 2^{2/\alpha}$ for the
theorem is clearly true when $1\leq y<2^{2/\alpha}$ by choosing a large enough
implicit constant. Recall that $1\leq q\leq y^{1-\alpha}$ for some $\alpha>0$.
Let $z=y^{\alpha/2}\geq 2$. Any squarefull number $n$ in $[x,x+y]$ can be
factored as

$n=\underbrace{p_{1}^{a_{1}}\cdots
p_{j}^{a_{j}}}_{b_{n}}\underbrace{p_{j+1}^{a_{j+1}}\cdots
p_{s}^{a_{s}}}_{d_{n}}\;\;\text{ with }\;\;p_{1}<p_{2}<\cdots<p_{s}$

where $j$ is the greatest index such that $p_{1}^{a_{1}}\cdots
p_{j}^{a_{j}}\leq z$. Hence, $b_{n}\leq z<b_{n}p_{j+1}^{a_{j+1}}$. Note that
$j$ may be $0$ (the product is an empty product) if $p_{1}^{a_{1}}>z$. In this
case, $b_{n}=1$ and $d_{n}=n$. Also, since $n\equiv r\pmod{q}$ with
$\text{gcd}(r,q)=1$, we must have $\text{gcd}(b_{n},q)=1=\text{gcd}(d_{n},q)$.

Case 1: $b_{n}>z^{1/2}$. As $q\leq y^{1-\alpha}$ and $z=y^{\alpha/2}$, the
number of such squarefull numbers is bounded by

$\mathop{\mathop{\sum\nolimits^{\prime}}_{z^{1/2}<b\leq
z}}_{\text{gcd}(b,q)=1}\mathop{\mathop{\sum_{x<n\leq x+y}}_{b|n}}_{n\equiv
r\;(\bmod{q})}1\leq\mathop{\sum\nolimits^{\prime}}_{z^{1/2}<b\leq
z}\Bigl{(}\frac{y/b}{q}+1\Bigr{)}\ll\frac{y}{qz^{1/4}}+z^{1/2}\ll_{\alpha}\frac{y}{\phi(q)\log
y}$ (9)

by (1) and Lemma 1.

Case 2: $b_{n}\leq z^{1/2}$ and $p_{-}(d_{n})\leq z^{1/2}$. Then $p_{j+1}\leq
z^{1/2}$ and $p_{j+1}^{a_{j+1}}>z^{1/2}$ which implies
$p_{j+1}^{-a_{j+1}}\leq\min(z^{-1/2},p_{j+1}^{-2})$ as $a_{j+1}\geq 2$. Hence,
the sum

$\sum_{p_{j+1}\leq z^{1/2}}\frac{1}{p_{j+1}^{a_{j+1}}}\leq\sum_{p_{j+1}\leq
z^{1/4}}z^{-1/2}+\sum_{z^{1/4}<p_{j+1}\leq
z^{1/2}}\frac{1}{p_{j+1}^{2}}\ll\frac{1}{z^{1/4}}.$

Therefore, by replacing $p_{j+1}^{a_{j+1}}$ with a generic $p^{a}$, the number
of squarefull numbers in this case is bounded by

$\mathop{\sum_{p\leq
z^{1/2}}}_{\text{gcd}(p,q)=1}\mathop{\mathop{\sum_{x<n\leq
x+y}}_{p^{a}|n}}_{n\equiv r\;(\bmod{q})}1\leq\sum_{p\leq
z^{1/2}}\Bigl{(}\frac{y/p^{a}}{q}+1\Bigr{)}\ll\frac{y}{qz^{1/4}}+z^{1/2}\ll_{\alpha}\frac{y}{\phi(q)\log
y}$ (10)

since $q\leq y^{1-\alpha}$ and $z=y^{\alpha/2}$.

Case 3: $b_{n}\leq z^{1/2}$ and $p_{-}(d_{n})>z^{1/2}$. As $q\leq
y^{1-\alpha}$ and $z=y^{\alpha/2}$, the number of such squarefull numbers is
bounded by

$\mathop{\mathop{\sum\nolimits^{\prime}}_{b\leq
z^{1/2}}}_{\text{gcd}(b,q)=1}\mathop{\mathop{\sum_{x/b<n/b\leq(x+y)/b}}_{p_{-}(n/b)>z^{1/2}}}_{(n/b)\equiv
r\overline{b}\;(\bmod q)}1\ll\mathop{\sum\nolimits^{\prime}}_{b\leq
z}\Bigl{(}\frac{y/b}{\phi(q)\log z}+z\Bigr{)}\ll\frac{y}{\phi(q)\log
z}+z^{3/2}\ll_{\alpha}\frac{y}{\phi(q)\log y}$ (11)

by (1), Lemma 2 and the convergence of the sum of reciprocal of squarefull
numbers (which follows from Lemma 1 for instance). Here $\overline{b}$ denotes
the multiplicative inverse of $b\;(\bmod{q})$, i.e. $b\overline{b}\equiv
1\pmod{q}$.

Combining (9), (10) and (11), we have Theorem 2.

##  4 Proof of Theorem 3

It is very similar to the proof of Theorem 2. So, we just highlight the
necessary adjustments. We set $q=1$ and $z=y^{1/3}$. The arguments for Case 1
and Case 2 are exactly the same as (9) and (10) and we get the bound

$\frac{y}{z^{1/4}}+z^{1/2}\ll y^{11/12}.$

It remains to deal with Case 3 where $b_{n}\leq z^{1/2}$ and
$z^{1/2}<p_{-}(d_{n})\leq y^{1/2}$ as the squarefull numbers are assumed to be
$y^{1/2}$-smooth. Thus, with $p:=p_{-}(d_{n})$ and $d_{n}:=p^{2}d$, the number
of squarefull numbers in this case is bounded by

$\displaystyle\mathop{\sum\nolimits^{\prime}}_{b\leq
z^{1/2}}\sum_{z^{1/2}<p\leq
y^{1/2}}\mathop{\sum_{x/b<p^{2}d\leq(x+y)/b}}_{p_{-}(d)>z^{1/2}}1=$
$\displaystyle\mathop{\sum\nolimits^{\prime}}_{b\leq
z^{1/2}}\sum_{z^{1/2}<p\leq
y^{1/2}}\mathop{\sum_{\frac{x}{bp^{2}}<d\leq\frac{x+y}{bp^{2}}}}_{p_{-}(d)>z^{1/2}}1$
$\displaystyle\ll$ $\displaystyle\mathop{\sum\nolimits^{\prime}}_{b\leq
z^{1/2}}\sum_{z^{1/2}<p\leq y^{1/2}}\Bigl{(}\frac{y/(bp^{2})}{\log
z}+z\Bigr{)}\ll\frac{y}{z^{1/2}\log z}+\frac{z^{5/4}y^{1/2}}{\log y}\ll
y^{11/12}$

by (1), Lemma 2 and the convergence of the sum of reciprocal of squarefull
numbers. The above bounds together yield Theorem 3.

##  5 Proof of Theorem 4

Given integer $k\geq 3$ and small real number $\epsilon>0$. We claim that, for
some sufficiently small constant $c_{\epsilon}>0$, the interval

$(x,x+c_{\epsilon}x^{1-(2+\epsilon)/k}]$

contains at most one $k$-full number for all sufficiently large $x\geq C$ (in
terms of $\epsilon$ and $c_{\epsilon}$). We shall abbreviate
$c_{\epsilon}x^{1-(2+\epsilon)/k}$ as $z$.

Suppose the contrary that the interval $(x,x+z]$ contains two $k$-full numbers
$b<c$. Then $c=a+b$ for some integer $0<a\leq z$. Let $d=\text{gcd}(a,b)$.
Then the integers $a/d$, $b/d$ and $c/d$ are pairwise relatively prime. Note
that $\kappa(n)\leq n^{1/k}$ for any $k$-full numbers. Apply the
$abc$-conjecture to the equation $a/d+b/d=c/d$, we get

$\displaystyle\frac{x}{d}<\frac{c}{d}$
$\displaystyle\ll_{\epsilon}\Bigl{(}\kappa\Bigl{(}\frac{a}{d}\frac{b}{d}\frac{c}{d}\Bigr{)}\Bigr{)}^{1+\epsilon/2}=\Bigl{(}\kappa\Bigl{(}\frac{a}{d}\frac{b}{d}\Bigr{)}\kappa\Bigl{(}\frac{c}{d}\Bigr{)}\Bigr{)}^{1+\epsilon/2}$
$\displaystyle\leq\Bigl{(}\kappa\Bigl{(}\frac{ab}{d}\Bigr{)}\kappa(c)\Bigr{)}^{1+\epsilon/2}\leq\Bigl{(}\frac{a^{1/k}b^{1/k}}{d^{1/k}}c^{1/k}\Bigr{)}^{1+\epsilon/2}\leq\frac{z^{(1+\epsilon/2)/k}(2x)^{2(1+\epsilon/2)/k}}{d^{(1+\epsilon/2)/k}}$

as $ab/d=\text{lcm}(a,b)$ is still $k$-full, $0<a\leq z$ and $b,c\leq 2x$.
Since $0<d\leq a\leq z$, the above implies

$x^{1-(2+\epsilon)/k}\ll_{\epsilon}z$

which contradicts with $z=c_{\epsilon}x^{1-(2+\epsilon)/k}$ when
$c_{\epsilon}>0$ is small enough. Thus, the claim is true.

Clearly, Theorem 4 is true for $1\leq y\leq C$ by picking an appropriate
implicit constant. Now, for $C\leq y\leq x$, the above claim implies that the
interval

$(x,x+c_{\epsilon}y^{1-(2+\epsilon)/k}]$

contains at most one $k$-full number. By dividing the interval $(x,x+y]$ into
subintervals of length $c_{\epsilon}y^{1-(2+\epsilon)/k}$, we obtain

$Q_{k}(x+y)-Q_{k}(x)\ll\frac{y}{c_{\epsilon}y^{1-(2+\epsilon)/k}}\cdot
1\ll_{\epsilon}y^{\frac{2+\epsilon}{k}}$

which gives Theorem 4.

##  6 Proof of Theorem 5

First, we suppose $y\leq x^{0.2}$. We claim that there is no non-trivial three
term arithmetic progression among the squarefull numbers in the interval
$(x,x+y]$ under the abc-conjecture. Suppose the contrary. Then we have three
squarefull numbers
$x<a_{1}^{2}b_{1}^{3}<a_{2}^{2}b_{2}^{3}<a_{3}^{2}b_{3}^{2}\leq x+y$ such that

$a_{1}^{2}b_{1}^{3}=a_{2}^{2}b_{2}^{3}-d\;\;\text{ and
}\;\;a_{3}^{2}b_{3}^{3}=a_{2}^{2}b_{2}^{3}+d$

for some positive integer $d$ with $2d\leq y$. Multiplying the above two
equations, we get

$a_{1}^{2}a_{3}^{2}b_{1}^{3}b_{3}^{3}=a_{2}^{4}b_{2}^{6}-d^{2}\;\;\text{ or
}\;\;a_{1}^{2}a_{3}^{2}b_{1}^{3}b_{3}^{3}+d^{2}=a_{2}^{4}b_{2}^{6}.$

Say $D^{2}=\text{gcd}(a_{2}^{4}b_{2}^{6},d^{2})$ as the numbers are perfect
squares. Then, the three integers

$\frac{a_{1}^{2}a_{3}^{2}b_{1}^{3}b_{3}^{3}}{D^{2}},\;\;\frac{d^{2}}{D^{2}},\;\;\frac{a_{2}^{4}b_{2}^{6}}{D^{2}}$

are pairwise relatively prime and we have the equation

$\frac{a_{1}^{2}a_{3}^{2}b_{1}^{3}b_{3}^{3}}{D^{2}}+\frac{d^{2}}{D^{2}}=\frac{a_{2}^{4}b_{2}^{6}}{D^{2}}.$

Now, by the abc-conjecture, we have

$\displaystyle\frac{x^{2}}{D^{2}}\leq\frac{a_{2}^{4}b_{2}^{6}}{D^{2}}$
$\displaystyle\ll_{\epsilon}\kappa\Bigl{(}\frac{a_{1}^{2}a_{3}^{2}b_{1}^{3}b_{3}^{3}}{D^{2}}\frac{d^{2}}{D^{2}}\frac{a_{2}^{4}b_{2}^{6}}{D^{2}}\Bigr{)}^{1+\epsilon}$
$\displaystyle\leq\kappa\bigl{(}a_{1}^{2}a_{3}^{2}b_{1}^{3}b_{3}^{3}\bigr{)}^{1+\epsilon}\kappa\Bigl{(}\frac{d^{2}}{D^{2}}\Bigr{)}^{1+\epsilon}\kappa\bigl{(}a_{2}^{4}b_{2}^{6}\bigr{)}^{1+\epsilon}$
$\displaystyle\leq(a_{1}b_{1}a_{2}b_{2}a_{3}b_{3})^{1+\epsilon}\Bigl{(}\frac{d}{D}\Bigr{)}^{1+\epsilon}\ll
x^{3/2+3\epsilon/2}\frac{y^{1+\epsilon}}{D^{1+\epsilon}}.$

Since $1\leq D\leq d\leq y$, this implies
$x^{1/2-3\epsilon/2}\ll_{\epsilon}D^{1-\epsilon}y^{1+\epsilon}\ll y^{2}\leq
x^{0.4}$ which is a contradiction for small enough $\epsilon$, say
$\epsilon=0.01$, and sufficiently large $x>C$ (in terms of the implicit
constant).

Clearly, the theorem is true for $1\leq y\leq C$ by picking an appropriate
implicit constant. So, we may assume $y>C$. Since arithmetic progressions are
invariant under translation, we may shift the interval $(x,x+y]$ to $(0,y]$.
Therefore, by Theorem 6, we have

$Q_{2}(x+y)-Q_{2}(x)\ll\frac{y}{\log^{1+c}y}$

which gives the theorem.

Now, if $y>x^{0.2}$, one can simply divide the interval $(x,x+y]$ into
subintervals of length $x^{0.2}$:

$(x,x+x^{0.2}]\cup(x+x^{0.2},x+2x^{0.2}]\cup\cdots\cup\Bigl{(}x+\Bigl{\lfloor}\frac{y}{x^{0.2}}\Bigr{\rfloor}x^{0.2},x+\Bigl{(}\Bigl{\lfloor}\frac{y}{x^{0.2}}\Bigr{\rfloor}+1\Bigr{)}x^{0.2}\Bigr{]}$

Then, over each interval $(x+ix^{0.2},x+(i+1)x^{0.2}]$, we have the bound

$Q_{2}(x+(i+1)x^{0.2})-Q_{2}(x+ix^{0.2})\ll\frac{x^{0.2}}{\log^{1+c}x}.$

Summing over $\lfloor\frac{y}{x^{0.2}}\rfloor+1$ of these intervals, we have

$Q_{2}(x+y)-Q_{2}(x)\ll\frac{y}{x^{0.2}}\cdot\frac{x^{0.2}}{\log^{1+c}x}\ll\frac{y}{\log^{1+c}y}$

which gives the theorem as well.

## References

  * [1] P.T. Bateman and E. Grosswald, On a theorem of Erdős and Szekeres, Illinois J. Math. 2 (1958), 88–98.
  * [2] T.F. Bloom and O. Sisask, Breaking the logarithmic barrier in Roth’s theorem on arithmetic progressions, preprint, arXiv:2007.03528.
  * [3] J.M. De Koninck, F. Luca and I.E. Shparlinski, Powerful numbers in short intervals. Bull. Austral. Math. Soc. 71 (2005), no. 1, 11–16.
  * [4] P. Erdős and G. Szekeres, Über die Anzahl der Abelschen Gruppen gegebener Ordnung und über ein verwandtes zahlentheoretisches Problem, Acta Univ. Szeged, 7 (1934-1935), 95–102.
  * [5] H. Halberstam and H.E. Richert, Sieve methods, London Mathematical Society Monographs, No. 4. Academic Press [Harcourt Brace Jovanovich, Publishers], London-New York, 1974. xiv+364 pp.
  * [6] H.Q. Liu, The number of cubefull numbers in an interval, Funct. Approx. Comment. Math. 43 (2010), part 2, 105–107.
  * [7] P. Shiu, A Brun-Titchmarsh theorem for multiplicative functions. J. Reine Angew. Math. 313 (1980), 161–170.
  * [8] O. Trifonov, Lattice points close to a smooth curve and squarefull numbers in short intervals, J. London Math. Soc. (2) 65 (2002), 303–319.

Mathematics Department
Kennesaw State University
Marietta, GA 30060
<EMAIL_ADDRESS>