# Many-Body Mobility Edge in Quantum Sun models

Konrad Pawlik<EMAIL_ADDRESS>Szkoła Doktorska Nauk Ścisłych
i Przyrodniczych, Uniwersytet Jagielloński, ulica Stanisława Łojasiewicza 11,
PL-30-348 Kraków, Poland  Instytut Fizyki Teoretycznej, Wydział Fizyki,
Astronomii i Informatyki Stosowanej, Uniwersytet Jagielloński, ulica
Stanisława Ł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
Lev Vidmar  Department of Theoretical Physics, J. Stefan Institute, SI-1000
Ljubljana, Slovenia  Department of Physics, Faculty of Mathematics and
Physics, University of Ljubljana, SI-1000 Ljubljana, Slovenia  Jakub
Zakrzewski<EMAIL_ADDRESS>Instytut Fizyki Teoretycznej, Wydział
Fizyki, Astronomii i Informatyki Stosowanej, Uniwersytet Jagielloński, ulica
Stanisława Łojasiewicza 11, PL-30-348 Kraków, Poland  Mark Kac Complex Systems
Research Center, Jagiellonian University in Krakow, PL-30-348 Kraków, Poland.

###### Abstract

The quantum sun model is an interacting model that exhibits sharp signatures
of ergodicity breaking phase transition. Here, we show that the model exhibits
a many-body mobility edge. We provide analytical arguments for its existence,
complemented by the state-of-the-art numerical simulations analysing gap
ratios, Thouless times as well as entanglement entropy of eigenstates. We also
introduce the quantum sun model with particle number conservation, and we
argue that it shares many similarities with his unrestricted predecessor.

### Introduction.

In recent years, there has been a rapid change of the perspective on Many body
localization (MBL) [1, 2, 3, 4] proposed as a robust mechanism leading to
breakdown of quantum thermalization in isolated many-body systems [5] in the
presence of disorder and interactions. Intensive theoretical studies (for
reviews see e.g. [6, 7, 8]) supplemented by experiments [9, 10, 11, 12, 13,
14] yielded a picture of MBL as of a dynamical phase of matter in which the
local integrals of motion [15, 16, 17, 18] preserve information about the
initial state, the transport is suppressed [6, 19], and the entanglement
spreads slowly [20, 21, 22, 23, 24]. Disordered Heisenberg spin-1/2 chain
served as one of the paradigmatic models in those investigations of MBL [25,
3, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44,
45, 46, 47, 48]. Hence, it came as a surprise when [49] pointed out that
numerical results for this model could be interpreted in terms of a slow-down
of quantum thermalization with increasing disorder strength rather than in
terms of a transition to an MBL phase, challenging the very existence of MBL
as of a dynamical phase of matter. Subsequent works [50, 51] showed that the
existence of MBL is not directly ruled out by the arguments of [49], arguing
that a convincing numerical demonstration of MBL [52] is beyond reach of
present day supercomputers [53, 54]. Those results sparkled a debate about the
status of MBL [54, 55, 56, 57, 58, 59], during which the role of slow, but
apparently persistent approach towards thermal equilibrium was emphasized [60,
61, 62, 63, 64, 65, 66].

Figure 1:  (a) Many-body density of states (DoS) as a function of energy
density $\varepsilon$. By increasing the system size $L$, the DoS gets
exponentially suppressed away from the middle of the spectrum. (b)
Visualization of the model described by (4) and (5). Particles inside the sun
are described by GOE matrix $R$, while the outside spin $j$ is perturbed by
the magnetic field $h_{j}$ and coupled to one of the spins of the sun with
distance, $u_{j}$, dependent coupling. (c) The many-body mobility edge is
manifested as a variation of the critical coupling $\alpha$ with
$\varepsilon$.

Recently, it was shown [67] that a clear numerical demonstration of ergodicity
breaking phase transition (EBPT) is possible for the ”quantum sun” model
(QSM), a model that was originally investigated in a context of quantum
avalanches [68, 69] seeded by rare ergodic regions in spin chains [70, 71, 72,
73, 74]. The QSM consists of a small delocalized core, described by a random
matrix from Gaussian orthogonal ensemble (GOE) [75] coupled to a chain of
spins by interaction which decays exponentially with the distance, c.f. Fig.
1(b). The position of EBPT can be approximated [69, 67, 76] due to a simple
analytic hybridization condition [68] which studies the scaling properties of
the interaction matrix element relative to the level spacing that is the
inverse of the density of states (DoS).

In this Letter, we argue that the hybridization condition does not only enable
a sharp prediction of the EBPT point in QSM, but it also implies an existence
of a mobility edge separating the spectrum of ergodic and localized
eigenstates. Mobility edges exist for noninteracting systems [77, 78, 79, 80]
and they were postulated also for MBL [81, 82, 83, 84, 85, 86, 87], but at the
same time there are theoretical arguments [88] against their existence in
interacting systems. Therefore, deciding whether many-body mobility edges are
stable, or whether they are only a transient feature that disappears in the
thermodynamic limit, is one of the key problems in non-equilibrium physics of
disordered many-body systems. For the QSM studied here, the key insight is
that the DoS is exponentially suppressed away from the middle of the spectrum,
as sketched in Fig. 1. We demonstrate emergence of the many-body mobility edge
in the original QSM, which does not conserve the total magnetization
$\hat{S}_{\rm tot}^{z}$. We then introduce another variant of the QSM that
conserves $\hat{S}_{\rm tot}^{z}$ and we argue that both models exhibits very
similar properties, verifying that the stability of EBPT in QSM is not
affected by the presence of additional constants of motion.

### Analytical arguments for the mobility edge.

The condition for the EBPT in the QSM of size $L$ is based on the
hybridization condition [68, 67]

${\cal G}\propto\alpha^{L}\sqrt{\rho}={\rm
const.}\;\;\;\rightarrow\;\;\;\alpha_{c}{\propto}\rho^{-\frac{1}{2L}}\;,$ (1)

where $\alpha$ is the amplitude of interaction matrix element and $\rho$ is
the DoS. Usually, one considers the middle of the spectrum ($E\approx 0$) for
which the DoS for spin-1/2 particles scales as $\rho(E=0)\propto 2^{L}$,
giving rise to the transition point $\alpha_{c}=1/\sqrt{2}$. In contrast to
the extended disordered spin chains, the presence of EBPT at
$\alpha\approx\alpha_{c}=1/\sqrt{2}$ in QSM can be demonstrated numerically in
a convincing fashion, as shown in [67] with a few indicators of ergodicity
breaking.

In many-body systems with local interactions, the DoS is well approximated by
a Gaussian shape with possible model dependent deviations [89, 90, 91, 92]

$\rho(E)=2^{L}\frac{1}{\sqrt{2\pi\sigma_{E}^{2}}}e^{-\frac{E^{2}}{2\sigma_{E}^{2}}},$
(2)

where $\sigma_{E}$ is the standard deviation of energy, which has a leading
scaling $\sigma_{E}=b\sqrt{L}$. Expressing the energy $E$ in terms of the
rescaled energy
$\varepsilon=(E-E_{\mathrm{min}})/(E_{\mathrm{max}}-E_{\mathrm{min}})$, where
$E_{\mathrm{min}}$ ($E_{\mathrm{max}}$) is the ground state (highest excited
state) energy, and using the system size dependence of $\sigma_{E}$ and
extensivity of energy, $E_{\mathrm{max}}-E_{\mathrm{min}}=aL$, we find the
leading exponential scaling in (2) is
$(2e^{-a^{2}(\varepsilon-1/2)^{2}/(2b^{2})})^{L}$ with constants $a$ and $b$
that are dependent on the details of the model. Using this scaling in the
criterion (1), we find the position of the transition between ergodic and
localized phases at rescaled energy $\varepsilon$ as

$\alpha_{c}(\varepsilon)=\frac{1}{\sqrt{2}}e^{\frac{a^{2}(\varepsilon-1/2)^{2}}{4b^{2}}}\;.$
(3)

This formula reproduces the analytic result
$\alpha_{c}(\varepsilon)=\frac{1}{\sqrt{2}}$ at $\varepsilon=1/2$ and predicts
the shape of the many-body mobility edge for any $\varepsilon\in[0,1]$. The
critical assumption of the formula (3) is that the interaction matrix element
scales with distance $L$ as $\alpha^{L}$. The criterion (3) arises due to a
competition of DoS with the interaction matrix element. In the following, we
start by investigating the EBPT in a QSM with U(1) symmetry, and then we
verify the assumptions behind (3) by investigating numerically the many-body
mobility edges in QSMs with and without the U(1) symmetry.

### The Quantum Sun models.

The QSM considered in [67] and depicted in Fig. 1(b), consists of a set of
${1/2}$-spins divided into two subsystems: $N$ spins in the delocalized core,
i.e., in the “sun”, enjoy all-to-all interactions, while $L$ other spins,
forming rays of the sun, do not interact among themselves. Each spin $j$ is
coupled only to one randomly chosen spin $n_{j}$ of the sun. The j-th outer
spin is positioned at distance $u_{j}$ from the sun, with $u_{j}$ chosen with
a uniform probability from the interval $[j-\zeta,j+\zeta]$, whereas its
coupling to the spin from the sun is set as ${g_{0}\alpha^{u_{j}}}$. The outer
spins are also exposed to a random magnetic field chosen with a uniform
distribution ${h_{j}\in[1-W,1+W]}$. Internal degrees of freedom of the sun are
described by a ${2^{N}\times 2^{N}}$ GOE matrix $R=\frac{\beta}{2}(B+B^{T})$,
where the matrix elements $B_{ij}$ are sampled independently from the standard
normal distribution. We follow the values of parameters chosen in [67] having:
${N=3,\zeta=0.2,g_{0}=1,\beta=0.3,W=0.5}$. The resulting Hamiltonian of the
original QSM is:

$\hat{H}=R\otimes\mathds{1}+\sum\limits_{j=1}^{L}g_{0}\alpha^{u_{j}}\hat{S}^{x}_{n_{j}}\hat{S}^{x}_{j}+\sum\limits_{j=1}^{L}h_{j}\hat{S}^{z}_{j}.$
(4)

We introduce the U(1) symmetric QSM that adds in a minimal way the conserved
quantity - the projection of total spin on the $z$ axis:
${\hat{S}^{z}_{\mathrm{tot}}=\sum_{j=1}^{L+N}\hat{S}_{j}^{z}}$. To this end
two modifications of the original Hamiltonian, (4) are necessary. First,
$\hat{S}^{z}_{\mathrm{tot}}$ conservation needs to be imposed on the couplings
of the outer spins to the sun. To accomplish this, each
$\hat{S}^{x}\hat{S}^{x}$ term is complemented with a $\hat{S}^{y}\hat{S}^{y}$
term with equal coupling, i.e.:
${\hat{S}^{x}_{n_{j}}\hat{S}^{x}_{j}\to\hat{S}^{x}_{n_{j}}\hat{S}^{x}_{j}+\hat{S}^{y}_{n_{j}}\hat{S}^{y}_{j}}$.
Secondly, the sun, $R$, should also commute with $\hat{S}^{z}_{\mathrm{tot}}$.
This condition makes the matrix representing $R$ block diagonal, with blocks
corresponding to eigenvalues of $\hat{S}^{z}_{\mathrm{tot}}$ and the block
ranks equal to the multiplicity of these eigenvalues. The conservation is
imposed by requiring that any matrix element outside those blocks vanishes,
leading to $R_{\mathrm{cons}}$ matrix. The Hamiltonian of the model with
conserved $\hat{S}^{z}_{\mathrm{tot}}$ reads explicitly:

$\hat{H}_{\mathrm{cons}}=R_{\mathrm{cons}}\otimes\mathds{1}+\sum\limits_{j=1}^{L}g_{0}\alpha^{u_{j}}(\hat{S}^{x}_{n_{j}}\hat{S}^{x}_{j}+\hat{S}^{y}_{n_{j}}\hat{S}^{y}_{j})+\sum\limits_{j=1}^{L}h_{j}\hat{S}^{z}_{j},$
(5)

manifesting the U(1) symmetry. In the following, our analysis of the above
model is restricted to the largest symmetry sector, i.e.,
${S^{z}_{\mathrm{tot}}=0}$ for ${N+L}$ even and ${S^{z}_{\mathrm{tot}}=-1/2}$
for ${N+L}$ odd.

### Ergodicy breaking phase transition.

To analyze the transition between ergodic and localized phases, we calculate
eigenvalues $E_{i}$ using full exact diagonalization for the U(1) symmetric
QSM (5). Results are obtained using $10000,3000$ realizations of disorder for
${L=7-12,13-14}$, respectively. Additionally, $3000$ realizations were taken
for $L=15-16$ using POLFED algorithm [54]. Note that the total length of the
chain is $N+L$ reaching $19$ in our calculations. From the spectrum we compute
the gap ratios ${r_{i}=\min({s_{i},s_{i+1}})/\max({s_{i},s_{i+1}})}$ [3, 93,
94], where ${s_{i}=E_{i+1}-E_{i}}$ is the level spacing. The gap ratios are
averaged over ${n_{\mathrm{ev}}=}{\min(\dim\mathcal{H}/10,500)}$ eigenvalues
from the middle of the spectrum, with $\dim\mathcal{H}$ being the dimension of
the Hilbert space, and finally the average over all available disorder
realizations is performed to obtain $\langle r\rangle$. For ergodic systems,
the value of this parameter is close to the GOE value $\langle r\rangle\approx
0.53$, while for localized systems $\langle r\rangle\approx 0.39$, as in the
Poisson ensemble [93].

Our results are presented in Fig. 2. Conclusions of [67] hold also for our
model (5): there exists $\alpha_{c}\approx 0.76$ at which $\langle r\rangle$
is practically independent of the system size, and all $\langle r\rangle$ vs
$\alpha$ curves intersect at this point. To identify the transition point
$\alpha_{c}$ more precisely, we investigate $\langle r\rangle$ as a function
of $L$ for fixed values of $\alpha$ close to $\alpha=0.76$, as shown in the
inset of Fig. 2. We observe a clear monotonous dependence of $\langle
r\rangle$ on $L$ for $\alpha=0.755$ and $\alpha=0.77$, which indicates that
the system is, respectively, in the localized and ergodic regime at those
values of $\alpha$. Data for $L\leq 14$ is inconclusive for $\alpha=0.765$.
Hence, using the POLFED algorithm [54] we additionally calculate $500$
energies in the middle of the spectrum for ${L=15-19,20}$, using ${1000,250}$
disorder realizations, respectively, showing that the point $\alpha=0.765$ is
in the ergodic phase. This allows us to narrow down the possible position of
the critical point to ${\alpha_{c}\in(0.755,0.765)}$.

Figure 2: Mean gap ratios $\langle r\rangle$ in the middle of the spectrum,
for different couplings $\alpha$ and number of outside spins $L$ for
Hamiltonian (5), cross practically in a single $L$ independent point at
$\alpha_{c}\approx 0.76$. The inset shows the small variation of $\langle
r\rangle$ with system size $L$ around critical $\alpha$ value. Black dotted
lines indicate the values expected from ergodic (GOE) and integrable (Poisson)
systems.

The position of the critical point can be studied more quantitatively by
extracting from Fig. 2 the system size dependent intersection point,
$\alpha_{c}(L)$, between $\langle r\rangle$ curves for system sizes $L-1$ and
$L+1$. Such an analysis, see [95], reveals a slight shift of the intersection
point with system size $L$, which is well described by
$\alpha_{c}(L)=\alpha_{c}(\infty)+a/L$ with $\alpha_{c}(\infty)=0.749(7)$.
This result is in agreement with the previously obtained bounds on
$\alpha_{c}$. Similar mild system size dependence of the intersection point is
also present in the original QSM [67]. Even though system size dependence of
the critical point is found, it is very subtle, in contrast to extended 1D
systems studied in the context of MBL, discussed in the introduction.
Moreover, we have verified that the same conclusions about ergodicity breaking
in the QSM (5) are obtained from analysis of the Thouless time extracted from
the spectral form factor [75, 49, 50, 67], see [95].

### Numerical test of the mobility edge.

The results obtained up till now clearly indicate the EBPT in the
thermodynamic limit at ${\alpha\approx 0.76}$ for the states in the middle of
the spectrum of (5). We now numerically test whether the QSMs exhibit a
mobility edge, i.e. whether the transition point depends on the energy. The
QSMs, both in its original version (4) and in our, the $\hat{S}_{z}$
preserving model (5), seem to be ideal candidates to study this problem due to
the sharpness of the transition.

Having the extensive set of eigenvalues for different system sizes for (5), we
diagonalized the QSM (4) for $L=6-12$ and considered at least $3000$ disorder
realizations. Then, we followed the procedure similar to the one described in
[83]. We determined the scaled eigenenergies,
${\varepsilon=(E-E_{\textrm{min}})/(E_{\textrm{max}}-E_{\textrm{min}})}$,
where $E_{\textrm{min}}$ ($E_{\textrm{max}}$) are the smallest (largest)
eigenvalue for a particular diagonalization. We obtain mean gap ratios within
the energy windows centered at specified values of
${\varepsilon\in\\{0.1,0.15,\ldots,0.9\\}}$ instead of the center of the
spectrum $\varepsilon=0.5$. For both models (4) and (5) we diagonalized two
more system sizes using POLFED [54] and at least $3000$ disorder realizations.

Figure 3: Phase diagram for the QSM. (a) U(1)-preserving QSM (5). (b) The
original QSM (4). Background colors show the value of $\langle r\rangle$ for
biggest fully diagonalized system size considered for each model, namely
$L=14$ for (a) and $L=12$ for the (b) (total system size is $L+3$ taking into
account the GOE core). Solid lines show position of system size dependent
boundary between phases $\alpha^{\varepsilon}_{c}(L)$ defined in the main
text. Dashed white lines give $L\rightarrow\infty$ extrapolation of the border
between localized (left, $\alpha<0.75$) and ergodic regimes (right), while the
band around it corresponds to the $2\sigma$ confidence interval of the
extrapolation, see [95] for details.

Figure 4: Entanglement entropy of eigenstates in U(1) symmetric QSM (5). (a)
Average entropy, $S$, for the bipartition described in the main text, averaged
over the middle of the spectrum (around $\varepsilon=0.5$), for different
system sizes $L$. Note the linear increase in the ergodic regime, flattening
close to the transition at $\alpha\approx 0.76$ and even decreasing for
smaller $\alpha$. (b) The rescaled entanglement entropies, $S_{r}$ (see text)
for different system sizes cross at EBPT point for $\varepsilon=0.5$. Inset
presents the similar crossing for energies around $\varepsilon=0.15$ for
$\alpha\approx 0.85$ indicating also the mobility edge.

The data obtained for different system sizes and at different $\varepsilon$
are presented in [95]. Here, we concentrate on the thermodynamic limit only.
We follow the same procedure as for $\varepsilon=0.5$ above, namely, we
extract the crossing point $\alpha_{c}^{\varepsilon}(L)$ between the average
gap ratios for systems of size $L-1$ and $L+1$, at given $\varepsilon$. More
consistent results were obtained when only a fraction of the energy window
size was considered, namely $0.25n_{\mathrm{ev}}$ for the original (4) and
$0.5n_{\mathrm{ev}}$ for the U(1) symmetric (5) models. Results presented in
Fig. 3 show that the further from the middle of the spectrum, the bigger the
value of $\alpha_{c}^{\varepsilon}(L)$. Moreover, with increasing system size,
the critical point for all values of $\varepsilon$ noticeably shifts towards
smaller values of $\alpha$. To extrapolate the results to the thermodynamic
limit in the QSM (4), we use a two parameter formula
$\alpha_{c}(\varepsilon)=A_{0}e^{\frac{(\varepsilon-1/2)^{2}}{A_{1}}}$
inspired by (3). We find that the two-parameter formula matches the numerical
data perfectly at any system size. This indicates that the mobility edge in
that version of the model is caused solely by the variation of the density of
states with $\varepsilon$. Finally, we extract the fit parameters for each
system size and extrapolate them to the thermodynamic limit leading to the
mobility edge shape presented in Fig. 3(b). The same reasoning, after slight
modification of the assumed density of states [95], applies to the U(1)
symmetric QSM (5), with the result shown in Fig. 3(a).

### Entanglement entropy of eigenstates.

EBPT should manifest itself also in the properties of eigenvectors. We
concentrate here on the entanglement entropy of eigenstates obtained by
splitting the system in two parts $A$ and $B$ and finding
$S(|\psi\rangle)=-\operatorname{Tr}_{A}\rho_{A}\ln\rho_{A}$ where
$\rho_{A}=\operatorname{Tr}_{B}|\psi\rangle\langle\psi|$. A small subtlety
arises as the parts cannot be equal, one of them contains the GOE sun. We
consider the splitting in the middle of the system, with the part $A$
consisting of the GOE sun and $\lfloor(N+L)/2\rfloor-N$ spins closest to it,
and the rest of the spins in the part $B$, here $\lfloor\cdot\rfloor$ denotes
the floor function. The dependence of average value of entropy, $S$, for the
middle part of the spectrum, as a function of the system size is shown in Fig.
4 for the $U(1)$ preserving model. In the extraction of the average value
around $\varepsilon$, states lying in the rescaled energy interval
$[\varepsilon-0.025,\varepsilon+0.025]$ were considered only. Note the visible
even-odd oscillations resulting from taking different $S^{z}_{\rm tot}$
sectors for even and odd $L$. We observe that $S$ for $\alpha>\alpha_{c}$
grows almost linearly with the system size, following a volume law. In
contrast, $S(L)$ becomes practically flat below $\alpha_{c}$ (area law
behavior) even showing a decrease with the system size. The latter effect is
due to the fact that with increasing $L$, the splitting between $A$ and $B$
performed in the middle of the system moves further from the GOE sun, which is
the source of entanglement growth in the system.

Scaling the entanglement entropy with the Page value $S_{P}$ [96, 97] i.e. the
typical value for a random vector from the Hilbert space, we obtained the
rescaled entanglement entropy, $S_{r}=S/S_{P}$ c.f. Fig. 4(b). For $S_{r}$
corresponding to the middle of the spectrum the curves for different system
sizes cross around $\alpha_{c}$ in a good correspondence with gap ratio
statistics. The inset shows a similar crossing for states around
$\varepsilon\approx 0.15$ which is significantly shifted to larger values of
$\alpha$ agreeing with the data of Fig. 3, providing another evidence for the
existence of the mobility edge. Similar data is obtained for the original QSM
[95].

### Conclusions.

The aims of the present work are twofold. First, we have shown that it is
possible to construct the QSM with a conserved total spin projection on the
$z$ axis. Translating the model from spin to particle language (via e.g.
Jordan-Wigner transformation) results in the model with particle number
conservation, more prone to experimental realizations. We have shown that the
resulting U(1) symmetric QSM shares with its predecessor all the essential
features resulting in the clear, unambiguous transition between ergodic and
localized phases. In particular, the introduced model allowed us to
demonstrate that the U(1) symmetry affects the EBPT transition in QSM in a
minor way. The second aim of this work is to provide the evidence for the
existence of the many-body mobility edge in the thermodynamic limit. This is
indeed confirmed for both considered QSMs by the extrapolation of our
numerical results, aided by analytic formulas derived from the hybridization
condition for QSM.

Over the last decades, the mobility edges in non-interacting quantum systems
have been a focus of theoretical studies [98, 99, 100, 77, 78, 101, 80] and
experimental investigations [102, 103, 104, 79, 13]. Our work, employs the
lack of finite size drifts in the quantum sun models, to show an unambiguous
evidence for many-body mobility-edge in a strongly interacting system. The
discussed mechanism for appearance of the stable many-body mobility edge
relies on the exponential suppression of many-body density of states away from
the middle of the spectrum. A question of whether this mechanism could lead to
stable many-body mobility edges in other thermodynamically large disordered
many-body systems remains an outstanding challenge for future investigations.

### Note added.

Since the original submission of this Letter a related work appeared [105]
which revealed the similarity between the original QSM and the ultrametric
ensemble. It finds multifractality of eigenstates at the transition
(consistent with the analysis in [76]) and Fock space localization of
eigenstates beyond EBPT. On the other hand, it has been argued [106] that
additional insight into the MBL transition may be obtained from density-
density correlation functions. It will be also a promising direction for
future studies of the QSM.

###### Acknowledgements.

We gratefully acknowledge discussions with W. De Roeck and A.S. Aramthottil.
We also acknowledge Poland’s high-performance computing infrastructure PLGrid
(HPC Centers: ACK Cyfronet AGH) for providing computer facilities and support
within computational grant no. PLG/2023/016370. This research was funded by
National Science Centre (Poland) under the OPUS call within the WEAVE
programme 2021/43/I/ST3/01142 (K.P. and J.Z.). For the purpose of Open Access,
the authors applied a CC-BY public copyright licence to any Author Accepted
Manuscript (AAM) version arising from this submission. A partial support by
the Strategic Programme Excellence Initiative at Jagiellonian University is
noted. 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 (RES-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). L.V.
acknowledges support by the Slovenian Research and Innovation Agency (ARIS),
Research core funding Grants No. P1-0044, N1-0273 and J1-50005. 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 artificial intelligence software.

## References

  * 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).
  * 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).
  * Oganesyan and Huse [2007] V. Oganesyan and D. A. Huse, Localization of interacting fermions at high temperature, Phys. Rev. B 75, 155111 (2007).
  * Pal and Huse [2010] A. Pal and D. A. Huse, Many-body localization phase transition, Phys. Rev. B 82, 174411 (2010).
  * 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), https://doi.org/10.1080/00018732.2016.1198134 .
  * 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).
  * 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).
  * Smith  _et al._ [2016] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, Many-body localization in a quantum simulator with programmable random disorder, Nature Physics 12, 907 (2016).
  * Lüschen  _et al._ [2017] H. P. Lüschen, P. Bordia, S. Scherg, F. Alet, E. Altman, U. Schneider, and I. Bloch, Observation of slow dynamics near the many-body localization transition in one-dimensional quasiperiodic systems, Phys. Rev. Lett. 119, 260401 (2017).
  * 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).
  * Kohlert  _et al._ [2019] T. Kohlert, S. Scherg, X. Li, H. P. Lüschen, S. Das Sarma, I. Bloch, and M. Aidelsburger, Observation of many-body localization in a one-dimensional system with a single-particle mobility edge, Phys. Rev. Lett. 122, 170403 (2019).
  * Léonard  _et al._ [2023] J. Léonard, S. Kim, M. Rispoli, A. Lukin, R. Schittko, J. Kwan, E. Demler, D. Sels, and M. Greiner, Probing the onset of quantum avalanches in a many-body localized system, Nature Physics 19, 481 (2023).
  * Serbyn  _et al._ [2013a] M. Serbyn, Z. Papić, and D. A. Abanin, Local conservation laws and the structure of the many-body localized states, Phys. Rev. Lett. 111, 127201 (2013a).
  * Huse  _et al._ [2014] D. A. Huse, R. Nandkishore, and V. Oganesyan, Phenomenology of fully many-body-localized systems, Phys. Rev. B 90, 174202 (2014).
  * 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).
  * 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).
  * Žnidarič  _et al._ [2016] M. Žnidarič, A. Scardicchio, and V. K. Varma, Diffusive and subdiffusive spin transport in the ergodic phase of a many-body localizable system, Phys. Rev. Lett. 117, 040601 (2016).
  * Chiara  _et al._ [2006] G. D. Chiara, S. Montangero, P. Calabrese, and R. Fazio, Entanglement entropy dynamics of heisenberg chains, Journal of Statistical Mechanics: Theory and Experiment 2006, P03001 (2006).
  * Žnidarič  _et al._ [2008] M. Žnidarič, T. Prosen, and P. Prelovšek, Many-body localization in the Heisenberg XXZ magnet in a random field, Phys. Rev. B 77, 064426 (2008).
  * Bardarson  _et al._ [2012] J. H. Bardarson, F. Pollmann, and J. E. Moore, Unbounded growth of entanglement in models of many-body localization, Phys. Rev. Lett. 109, 017202 (2012).
  * Serbyn  _et al._ [2013b] M. Serbyn, Z. Papić, and D. A. Abanin, Universal slow growth of entanglement in interacting strongly disordered systems, Phys. Rev. Lett. 110, 260601 (2013b).
  * Iemini  _et al._ [2016] F. Iemini, A. Russomanno, D. Rossini, A. Scardicchio, and R. Fazio, Signatures of many-body localization in the dynamics of two-site entanglement, Physical Review B 94, 214206 (2016).
  * 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).
  * Berkelbach and Reichman [2010] T. C. Berkelbach and D. R. Reichman, Conductivity of disordered quantum lattice models at infinite temperature: Many-body localization, Phys. Rev. B 81, 224429 (2010).
  * Serbyn  _et al._ [2015] M. Serbyn, Z. Papić, and D. A. Abanin, Criterion for many-body localization-delocalization phase transition, Phys. Rev. X 5, 041047 (2015).
  * Agarwal  _et al._ [2015] K. Agarwal, S. Gopalakrishnan, M. Knap, M. Müller, and E. Demler, Anomalous diffusion and griffiths effects near the many-body localization transition, Phys. Rev. Lett. 114, 160401 (2015).
  * Bera  _et al._ [2015] S. Bera, H. Schomerus, F. Heidrich-Meisner, and J. H. Bardarson, Many-body localization characterized from a one-particle perspective, Phys. Rev. Lett. 115, 046603 (2015).
  * Luitz  _et al._ [2016] D. J. Luitz, N. Laflorencie, and F. Alet, Extended slow dynamical regime close to the many-body localization transition, Phys. Rev. B 93, 060201 (2016).
  * Luitz and Bar Lev [2016] D. J. Luitz and Y. Bar Lev, Anomalous thermalization in ergodic systems, Phys. Rev. Lett. 117, 170404 (2016).
  * 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).
  * Colmenarez  _et al._ [2019] L. Colmenarez, P. A. McClarty, M. Haque, and D. J. Luitz, Statistics of correlations functions in the random Heisenberg chain, SciPost Phys. 7, 64 (2019).
  * 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).
  * Khemani  _et al._ [2017] V. Khemani, S. P. Lim, D. N. Sheng, and D. A. Huse, Critical properties of the many-body localization transition, Phys. Rev. X 7, 021013 (2017).
  * Serbyn  _et al._ [2017] M. Serbyn, Z. Papić, and D. A. Abanin, Thouless energy and multifractality across the many-body localization transition, Phys. Rev. B 96, 104201 (2017).
  * Kjäll [2018] J. A. Kjäll, Many-body localization and level repulsion, Phys. Rev. B 97, 035163 (2018).
  * 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).
  * Doggen and Mirlin [2019] E. V. H. Doggen and A. D. Mirlin, Many-body delocalization dynamics in long Aubry-André quasiperiodic chains, Phys. Rev. B 100, 104203 (2019).
  * Weiner  _et al._ [2019] F. Weiner, F. Evers, and S. Bera, Slow dynamics and strong finite-size effects in many-body localization with random and quasiperiodic potentials, Phys. Rev. B 100, 104204 (2019).
  * Sierant and Zakrzewski [2019] P. Sierant and J. Zakrzewski, Level statistics across the many-body localization transition, Phys. Rev. B 99, 104205 (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).
  * 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).
  * Schiulaz  _et al._ [2020] M. Schiulaz, E. J. Torres-Herrera, F. Pérez-Bernal, and L. F. Santos, Self-averaging in many-body quantum systems out of equilibrium: Chaotic systems, Phys. Rev. B 101, 174312 (2020).
  * Torres-Herrera  _et al._ [2020] E. J. Torres-Herrera, G. De Tomasi, M. Schiulaz, F. Pérez-Bernal, and L. F. Santos, Self-averaging in many-body quantum systems out of equilibrium: Approach to the localized phase, Phys. Rev. B 102, 094310 (2020).
  * Schiulaz  _et al._ [2019] M. Schiulaz, E. J. Torres-Herrera, and L. F. Santos, Thouless and relaxation time scales in many-body quantum systems, Phys. Rev. B 99, 174313 (2019).
  * Hopjan  _et al._ [2021] M. Hopjan, G. Orso, and F. Heidrich-Meisner, Detecting delocalization-localization transitions from full density distributions, Phys. Rev. B 104, 235112 (2021).
  * Š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).
  * Sierant  _et al._ [2020a] P. Sierant, D. Delande, and J. Zakrzewski, Thouless time analysis of anderson and many-body localization transitions, Phys. Rev. Lett. 124, 186601 (2020a).
  * 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).
  * 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._ [2020b] P. Sierant, M. Lewenstein, and J. Zakrzewski, Polynomially filtered exact diagonalization approach to many-body localization, Phys. Rev. Lett. 125, 156601 (2020b).
  * Š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).
  * Laflorencie  _et al._ [2020] N. Laflorencie, G. Lemarié, and N. Macé, Chain breaking and Kosterlitz-Thouless scaling at the many-body localization transition in the random-field Heisenberg spin chain, Phys. Rev. Research 2, 042033 (2020).
  * Sels and Polkovnikov [2023] D. Sels and A. Polkovnikov, Thermalization of dilute impurities in one-dimensional spin chains, Phys. Rev. X 13, 011041 (2023).
  * Morningstar  _et al._ [2022] A. Morningstar, L. Colmenarez, V. Khemani, D. J. Luitz, and D. A. Huse, Avalanches and many-body resonances in many-body localized systems, Phys. Rev. B 105, 174205 (2022).
  * Sels [2022] D. Sels, Bath-induced delocalization in interacting disordered spin chains, Phys. Rev. B 106, L020202 (2022).
  * Kiefer-Emmanouilidis  _et al._ [2020] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Evidence for unbounded growth of the number entropy in many-body localized phases, Phys. Rev. Lett. 124, 243601 (2020).
  * Kiefer-Emmanouilidis  _et al._ [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).
  * Luitz and Lev [2020] D. J. Luitz and Y. B. Lev, Absence of slow particle transport in the many-body localized phase, Phys. Rev. B 102, 100202 (2020).
  * Ghosh and Žnidarič [2022] R. Ghosh and M. Žnidarič, Resonance-induced growth of number entropy in strongly disordered systems, Phys. Rev. B 105, 144203 (2022).
  * Sierant and Zakrzewski [2022] P. Sierant and J. Zakrzewski, Challenges to observation of many-body localization, Phys. Rev. B 105, 224203 (2022).
  * Sierant  _et al._ [2023a] P. Sierant, T. Chanda, M. Lewenstein, and J. Zakrzewski, Slow dynamics of a mobile impurity interacting with an anderson insulator, Phys. Rev. B 107, 144201 (2023a).
  * Szołdra  _et al._ [2023] T. Szołdra, P. Sierant, M. Lewenstein, and J. Zakrzewski, Tracking locality in the time evolution of disordered systems, Phys. Rev. B 107, 054204 (2023).
  * Šuntajs and Vidmar [2022] J. Šuntajs and L. Vidmar, Ergodicity breaking transition in zero dimensions, Phys. Rev. Lett. 129, 060602 (2022).
  * De Roeck and Huveneers [2017] W. De Roeck and F. Huveneers, Stability and instability towards delocalization in many-body localization systems, Phys. Rev. B 95, 155129 (2017).
  * Luitz  _et al._ [2017] D. J. Luitz, F. Huveneers, and W. De Roeck, How a small quantum bath can thermalize long localized chains, Phys. Rev. Lett. 119, 150602 (2017).
  * Dumitrescu  _et al._ [2017] P. T. Dumitrescu, R. Vasseur, and A. C. Potter, Scaling theory of entanglement at the many-body localization transition, Phys. Rev. Lett. 119, 110604 (2017).
  * Goremykina  _et al._ [2019] A. Goremykina, R. Vasseur, and M. Serbyn, Analytically solvable renormalization group for the many-body localization transition, Phys. Rev. Lett. 122, 040601 (2019).
  * 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).
  * 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).
  * Haake [2010] F. Haake,  _Quantum Signatures of Chaos_ (Springer, Berlin, 2010).
  * Hopjan and Vidmar [2023] M. Hopjan and L. Vidmar, Scale invariant survival probability at eigenstate transitions, arXiv:2212.13888 (2023).
  * Delande and Orso [2014] D. Delande and G. Orso, Mobility edge for cold atoms in laser speckle potentials, Phys. Rev. Lett. 113, 060601 (2014).
  * Pasek  _et al._ [2017] M. Pasek, G. Orso, and D. Delande, Anderson localization of ultracold atoms: Where is the mobility edge?, Phys. Rev. Lett. 118, 170403 (2017).
  * Lüschen  _et al._ [2018] H. P. Lüschen, S. Scherg, T. Kohlert, M. Schreiber, P. Bordia, X. Li, S. Das Sarma, and I. Bloch, Single-particle mobility edge in a one-dimensional quasiperiodic optical lattice, Phys. Rev. Lett. 120, 160404 (2018).
  * Li and Das Sarma [2020] X. Li and S. Das Sarma, Mobility edge and intermediate phase in one-dimensional incommensurate lattice potentials, Phys. Rev. B 101, 064203 (2020).
  * Altshuler  _et al._ [1997] B. L. Altshuler, Y. Gefen, A. Kamenev, and L. S. Levitov, Quasiparticle lifetime in a finite system: A nonperturbative approach, Phys. Rev. Lett. 78, 2803 (1997).
  * Kjäll  _et al._ [2014] J. A. Kjäll, J. H. Bardarson, and F. Pollmann, Many-body localization in a disordered quantum Ising chain, Phys. Rev. Lett. 113, 107204 (2014).
  * Luitz  _et al._ [2015] D. J. Luitz, N. Laflorencie, and F. Alet, Many-body localization edge in the random-field Heisenberg chain, Phys. Rev. B 91, 081103 (2015).
  * Mondragon-Shem  _et al._ [2015] I. Mondragon-Shem, A. Pal, T. L. Hughes, and C. R. Laumann, Many-body mobility edge due to symmetry-constrained dynamics and strong interactions, Phys. Rev. B 92, 064203 (2015).
  * Villalonga  _et al._ [2018] B. Villalonga, X. Yu, D. J. Luitz, and B. K. Clark, Exploring one-particle orbitals in large many-body localized systems, Phys. Rev. B 97, 104406 (2018).
  * Enss  _et al._ [2017] T. Enss, F. Andraschko, and J. Sirker, Many-body localization in infinite chains, Phys. Rev. B 95, 045121 (2017).
  * Chanda  _et al._ [2020] T. Chanda, P. Sierant, and J. Zakrzewski, Many-body localization transition in large quantum spin chains: The mobility edge, Phys. Rev. Research 2, 032045 (2020).
  * De Roeck  _et al._ [2016] W. De Roeck, F. Huveneers, M. Müller, and M. Schiulaz, Absence of many-body mobility edges, Phys. Rev. B 93, 014203 (2016).
  * French and Wong [1970] J. French and S. Wong, Validity of random matrix theories for many-particle systems, Physics Letters B 33, 449 (1970).
  * Bohigas and Flores [1971] O. Bohigas and J. Flores, Spacing and individual eigenvalue distributions of two-body random hamiltonians, Physics Letters B 35, 383 (1971).
  * Hartmann  _et al._ [2004] M. Hartmann, G. mahler, and O. Hess, Gaussian quantum fluctuations in interacting many particle systems, Letters in Mathematical Physics 68, 103 (2004).
  * Hartmann  _et al._ [2005] M. Hartmann, G. Mahler, and O. Hess, Spectral densities and partition functions of modular quantum systems as derived from a central limit theorem, Journal of Statistical Physics 119, 1139 (2005).
  * Atas  _et al._ [2013] Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, Distribution of the ratio of consecutive level spacings in random matrix ensembles, Phys. Rev. Lett. 110, 084101 (2013).
  * Giraud  _et al._ [2022] O. Giraud, N. Macé, E. Vernier, and F. Alet, Probing symmetries of quantum many-body systems through gap ratio statistics, Phys. Rev. X 12, 011006 (2022).
  * [95] See Supplementary Material at [URL will be inserted by publisher] for details on calculation of spectral form factor, Thouless time, scaling of critical points and discussion of statistical uncertainties. It contains references [107], [108], [109], [110], [111], [112].
  * Page [1993] D. N. Page, Average entropy of a subsystem, Phys. Rev. Lett. 71, 1291 (1993).
  * Bianchi  _et al._ [2022] E. Bianchi, L. Hackl, M. Kieburg, M. Rigol, and L. Vidmar, Volume-Law Entanglement Entropy of Typical Pure Quantum States, PRX Quantum 3, 030201 (2022).
  * Das Sarma  _et al._ [1988] S. Das Sarma, S. He, and X. C. Xie, Mobility edge in a model one-dimensional potential, Phys. Rev. Lett. 61, 2144 (1988).
  * Biddle and Das Sarma [2010] J. Biddle and S. Das Sarma, Predicted mobility edges in one-dimensional incommensurate optical lattices: An exactly solvable model of anderson localization, Phys. Rev. Lett. 104, 070601 (2010).
  * Biddle  _et al._ [2011] J. Biddle, D. J. Priour, B. Wang, and S. Das Sarma, Localization in one-dimensional lattices with non-nearest-neighbor hopping: Generalized anderson and aubry-andré models, Phys. Rev. B 83, 075105 (2011).
  * Li  _et al._ [2017] X. Li, X. Li, and S. Das Sarma, Mobility edges in one-dimensional bichromatic incommensurate potentials, Phys. Rev. B 96, 085119 (2017).
  * Kondov  _et al._ [2011] S. S. Kondov, W. R. McGehee, J. J. Zirbel, and B. DeMarco, Three-dimensional anderson localization of ultracold matter, Science 334, 66 (2011).
  * Jendrzejewski  _et al._ [2012] F. Jendrzejewski, A. Bernard, K. Müller, P. Cheinet, V. Josse, M. Piraud, L. Pezzé, L. Sanchez-Palencia, A. Aspect, and P. Bouyer, Three-dimensional localization of ultracold atoms in an optical disordered potential, Nature Physics 8, 398 (2012).
  * Semeghini  _et al._ [2015] G. Semeghini, M. Landini, P. Castilho, S. Roy, G. Spagnolli, A. Trenkwalder, M. Fattori, M. Inguscio, and G. Modugno, Measurement of the mobility edge for 3d anderson localization, Nature Physics 11, 554 (2015).
  * Šuntajs  _et al._ [2024] J. Šuntajs, M. Hopjan, W. De Roeck, and L. Vidmar, Similarity between a many-body quantum avalanche model and the ultrametric random matrix model, Phys. Rev. Res. 6, 023030 (2024).
  * Colbois  _et al._ [2024] J. Colbois, F. Alet, and N. Laflorencie, Interaction-driven instabilities in the random-field xxz chain (2024), arXiv:2403.09608 [cond-mat.dis-nn] .
  * Gómez  _et al._ [2002] J. M. G. Gómez, R. A. Molina, A. Relaño, and J. Retamosa, Misleading signatures of quantum chaos, Phys. Rev. E 66, 036209 (2002).
  * Mehta [1990] M. L. Mehta,  _Random Matrices_ (Elsevier, Amsterdam, 1990).
  * Cotler  _et al._ [2017] J. Cotler, N. Hunter-Jones, J. Liu, and B. Yoshida, Chaos, complexity, and random matrices, Journal of High Energy Physics 2017, 48 (2017).
  * Sierant  _et al._ [2023b] P. Sierant, M. Lewenstein, A. Scardicchio, and J. Zakrzewski, Stability of many-body localization in floquet systems, Phys. Rev. B 107, 115132 (2023b).
  * Cowan [1998] G. Cowan,  _Statistical data analysis_ (Oxford University Press, USA, 1998).
  * Cardy [1996] J. Cardy,  _Scaling and renormalization in statistical physics_ , Vol. 5 (Cambridge University Press, 1996).

##  I Supplementary material for Many-body mobility edge in Quantum Sun models

##  II Intermediate statistics

Both average gap $\langle r\rangle$ ratio and entanglement entropy $S$ of
eigenstates in the QSM reveal, that in the transition point between ergodic
and localized regimes almost no system size dependence can be seen. On the
other hand, more detailed analysis of the spectral statistics shows that
probability distributions of both level spacing $s$ and gap ratio $r$ approach
respectively Poissonian prediction in the localized phase, and GOE prediction
in an ergodic one. At a finite system size those distributions are in fact
intermediate between both model cases of GOE and Poissonian statistics, thus
it constitutes a valid question what will be the fate of those probability
distributions in the critical point once the thermodynamic limit is
approached.

To investigate this question, we aim to calculate numerical distributions of
the spacing ratio $P(s)$, although doing so requires first getting rid of the
influence of the DoS present in the spectrum. After obtaining an ordered set
of eigenvalues $\\{E_{n}\\}$, this can be achieved through so-called unfolding
of the spectrum procedure and an ordered set of unfolded eigenvalues
$\\{\widetilde{E}_{n}\\}$ is obtained. Unfolding is done by calculating the
cumulative density of states function ${G(E)=\sum_{n}\Theta(E-E_{n})}$, where
$\Theta$ is the Heaviside step function, and then separating it into
continuous and fluctuating part ${G(E)=\overline{G}(E)+\delta G(E)}$, where
$\overline{G}(E)$ is taken as a polynomial of small degree $n_{p}$ fitted to
$G(E)$. Degree $n_{p}$ is chosen as the smallest one that gives average
spacing $\approx 1$ after the unfolding, which corresponds in all considered
cases to $n_{p}=3,4$. Unfolded eigenvalues are defined as
$\widetilde{E}_{n}=\overline{G}(E_{n})$. In principle, the continuous part
$\overline{G}$ might exhibit nonphysical regions, where density of states
would be decreasing. To get rid of this numerical artifact $\widetilde{E}_{n}$
is included in calculation of SFF only if, after adding it, the set
$\\{\widetilde{E}_{n}\\}$ is still ordered in nondecreasing manner.

After unfolding, spacing $s$ is calculated for $n_{\rm
ev}=\min\\{500,\dim\mathcal{H}/10\\}$ energies from the middle of the spectrum
to give a numerical distribution $P(s)$, shown for the critical value of
$\alpha$ in Fig. S.1. Numerical distributions were fitted with $\nu-h$
distribution [43] (originally called $\beta-h$, here renamed to avoid
confusion with parameter $\beta$ of the main text), interpolating between GOE
$(\nu=1,h\to\infty)$ and Poisson $(\nu=0)$ statistics, recently shown to well
describe the EBPT.

Figure S.1: Intermediate level statistics close to the transition point.
Colored lines show numerical distributions. (a) Original quantum sun model at
$\alpha=0.75$. (b) U(1) symmetric model at $\alpha=0.76$. Black lines show
GOE, Poisson and $\nu-h$ distribution (see text), the latter of which was
fitted to the numerical data, giving $\nu$ and $h$ shown in the corresponding
legends.

For both studied variants of the QSM, the level spacing distribution did not
show any signs of system size dependence at the critical point, and excellent
agreement with the theoretical $\nu-h$ distribution was observed. Moreover,
$0<\nu<1$ in both models proves that critical statistics is neither of chaotic
nor integrable system. Investigation of the tails of numerical distributions
reveals that they perfectly follow the $\nu-h$ fit.

##  III Spectral form factor

To further confirm presence of the localization transition, we use a different
indicator of ergodicity breaking, the spectral form factor (SFF) [75, 49, 50,
67], defined as:

$K(\tau)=\frac{1}{Z}\left\langle\left|\sum\limits_{n=1}^{\dim\mathcal{H}}\rho(\widetilde{E}_{n})e^{-i2\pi\widetilde{E}_{n}\tau}\right|^{2}\right\rangle,$
(S.1)

where $\widetilde{E}_{n}$ are eigenvalues after the unfolding procedure [107],
$\rho(\widetilde{E})$ is a Gaussian weight diminishing the importance of
spectral edges, $Z$ is a normalization constant, such that ${K(\tau\gg
1)\approx 1}$ and the average is performed over disorder realizations. The SFF
is known analytically for GOE [108] as
${K_{\mathrm{GOE}}(\tau)=2\tau-\tau\ln(1+2\tau)}$ for ${\tau\leq 1}$ and
${K_{\mathrm{GOE}}(\tau)=2-\tau\ln[(2\tau+1)/(2\tau-1)]}$ for ${\tau>1}$. For
many-body systems enjoying time reversal symmetry there exists a time
$\tau_{\mathrm{Th}}$, called a Thouless time, beyond which the SFF follows the
GOE prediction. It is worth stressing, that due to the unfolding procedure,
$\tau$ present in the definition (S.1) is in fact dimensionless and related to
the real physical time via the Heisenberg time, $t_{H}$ (the inverse of the
mean level spacing) as $\tau=t/t_{H}$.

In calculation of SFF, procedure of [67] is followed. After unfolding the
spectrum, SFF (S.1) is calculated using a spectral filter diminishing the
influence of spectral edges
$\rho(\widetilde{E})=\exp[{-(\widetilde{E}-\langle\widetilde{E}\rangle)^{2}/(2\eta^{2}\sigma_{\widetilde{E}}^{2})}]$,
with $\langle\widetilde{E}\rangle$ and $\sigma_{\widetilde{E}}^{2}$ being mean
and variance in a given realization respectively. Parameter $\eta$ controls
effective fraction of included eigenvalues, $\eta=0.5$ is taken. To carry out
the sum in (S.1), a nonuniform fast Fourier transform algorithm is used.
Normalization $K(\tau\gg 1)\approx 1$ in (S.1) is assured by putting
$Z=\left\langle\sum_{n=1}^{\dim\mathcal{H}}|\rho(\widetilde{E}_{n})|^{2}\right\rangle$.

Figure S.2: SFF in the U(1)-preserving QSM (5) for ${\alpha=0.765}$. Small
sizes data (a) are from full exact diagonalization, larger $L$ data (b)
obtained using POLFED algorithm [54] finding $2500$ eigenvalues in the middle
of the spectrum. Insets show zooms on corresponding main panels. The
horizontal dotted line corresponds to the value in the constant regime
${K(\tau)\approx 0.51}$. Circles indicate the obtained values of Thouless
times.

Following [67], $K(\tau)$ was calculated for $5000$ times
$\tau_{i}\in[1/(2\pi\dim\mathcal{H}),5]$, equidistant in logarithmic scale. As
shown in the literature [109], the SFF does not have self-averaging property,
so it has fluctuations that do not decrease with increasing number of
realizations. Because of that, an additional moving time average is performed
to obtain $\overline{K}(\tau_{i})=\sum_{j=i-w}^{i+w}K(\tau_{j})$ for $w=25$,
which reduces the list of $\tau_{i}$ to $4950$ points. The difference from GOE
is given by:

$\Delta
K(\tau)=\left|\ln\left(\frac{\overline{K}(\tau)}{K_{\mathrm{GOE}}(\tau)}\right)\right|.$
(S.2)

The Thouless time $\tau_{\mathrm{Th}}$ is defined to be the biggest time for
which $\Delta K(\tau)\geq\epsilon$, where, $\epsilon=0.06$ was chosen based on
the criterion that $\Delta K(\tau_{\mathrm{Th}})$ must be larger than the
noise.

The analysis of the SFF shape (Fig. S.2) for a value close to the critical
region ${\alpha=0.765}$ reveals that the SFF does not exhibit dependence on
system size. Here, apart from exact diagonalization, we use the POLFED
algorithm [54] to reach larger system sizes up to $L=20$. In the latter case,
we extract $2500$ eigenvalues in the middle of the spectrum, with ${1000,250}$
realizations for ${L=15-19,20}$ respectively. This does not allow determining
$K(\tau)$ for small $\tau$, as the regime before the initial drop (e.g.
$\tau<10^{-3}$ on right panel of Fig. S.2) carries no physical information.
However, this is not an interesting range of $\tau$ values. Right after the
initial drop $K(\tau)$ is constant and $L$ independent over a range of $\tau$
values.

The behavior of SFF observed closely resembles the one reported for an
original quantum sun model [67] showing, without any doubt, that the
conservation of global $S_{z}$ projection in the model (5) does not affect the
character of the ergodicity breaking transition.

In case of larger systems $L=16-20$ instead of the whole spectrum only
$n_{e}=2500$ eigenvalues from the middle of the spectrum are included in
(S.1), which in general for $n_{e}\ll\dim\mathcal{H}$ has the effect of
shifting the initial drop of SFF to bigger time $\tau_{n_{e}}$ [50], but for
$\tau\gg\tau_{n_{e}}$, the SFF is correctly reproduced. As seen in Fig. S.2,
$\tau_{n_{e}=2500}\approx 10^{-3}$, and since the region interesting for
extraction of Thouless time is $\tau>10^{-1}$, so the obtained value of
$\tau_{\mathrm{Th}}$ is not affected by limited number of eigenvalues.

Figure S.3: Values of $t_{H}/t_{\mathrm{Th}}=1/\tau_{\mathrm{Th}}$ for
different couplings $\alpha$ and number of outside spins $L$ for the
U(1)-preserving QSM (5).

The numerically obtained SFF easily allows for extracting the Thouless times
for all analyzed values of $\alpha\in\\{0.6,0.61,\ldots,0.9\\}$. The results
for different system sizes are presented in Fig. S.3. The data cross almost at
a single, size independent point ${\alpha\approx 0.76}$, in an excellent
agreement with the gap ratio analysis.

##  IV Calculation of statistical uncertainties

In our consideration of the MBL phase transition in the main text, we
calculate the gap ratios $r$ and average them over a small energy window
around the rescaled energy $\varepsilon$, obtaining a single realization
average $r_{S}$. It is known that single realization average is not self-
averaging quantity [46, 45], meaning that the variance $\langle(r_{S}-\langle
r_{S}\rangle)^{2}\rangle$, where $\langle\cdot\rangle$ denotes disorder
average, is not necessarily decreasing when system size is increased. Instead,
disorder average $\langle r_{S}\rangle$ is considered, denoted $\langle
r\rangle$ for short. Assuming that $r_{S}$ for different realizations are
uncorrelated, the standard deviation of $\langle r\rangle$ reads [110]:

$\sigma_{r}=\frac{\langle(r_{S}-\langle
r\rangle)^{2}\rangle^{1/2}}{N_{\mathrm{dis}}^{1/2}},$ (S.3)

where $N_{\mathrm{dis}}$ is the number of disorder realizations used in
calculation of $\langle r\rangle$.

Figure S.4: Function $\chi^{2}_{L=11}$ in the vicinity of the crossing point
of $\langle r\rangle$ vs $\alpha$ curves at $\varepsilon=0.5$ for the system
sizes $L=10$ and $L=12$ in the original QSM (4). Horizontal line segment shows
$\chi^{2}=1$ and positions of roots of this equation $\alpha_{\pm}$ are shown
with dotted lines.

Having uncertainties of data points on $\langle r\rangle$ vs $\alpha$ curves,
we can estimate statistical errors for $\alpha_{c}^{\varepsilon}(L)$, which
requires determining the crossing point of two polynomials $f_{L-1}(\alpha)$
and $f_{L+1}(\alpha)$, which approximate $\langle r\rangle$ vs $\alpha$ curves
close to the critical point for $L-1$ and $L+1$, respectively. Using standard
error propagation, we calculate errors of $f_{L-1},f_{L+1}$ as a function of
$\alpha$, obtaining $\sigma_{L-1}(\alpha)$ and $\sigma_{L+1}(\alpha)$
correspondingly. Based on those functions the $\chi^{2}_{L}$ function is
defined:

$\chi^{2}_{L}(\alpha)=\frac{[f_{L-1}(\alpha)-f_{L+1}(\alpha)]^{2}}{\sigma_{L-1}^{2}(\alpha)+\sigma_{L+1}^{2}(\alpha)}.$
(S.4)

Exactly at the intersection point $\chi^{2}_{L}=0$, and it increases
approximately quadratically with the distance from minimum (see Fig. S.4),
thus the equation $\chi^{2}_{L}(\alpha)=1$ has two solutions
$\alpha_{-}<\alpha_{+}$. The quantity $(\alpha_{+}-\alpha_{-})/2$ has
statistical interpretation [111] as the standard deviation of the position of
the crossing point $\alpha_{c}^{\varepsilon}(L)$. All other instances of
uncertainty calculation follow simply standard uncertainty propagation.

##  V Extraction and scaling of critical point

Considering a plot of mean gap ratios $\langle r\rangle$ around the middle of
the spectrum in the QSM (see Fig. 2 in the main text) it seems that all curves
intersect at a single point, thus suggesting total independence of system size
at the critical point between the ergodic and MBL phase. Instead of the middle
of the spectrum, gap ratio can be calculated around rescaled energy
$\varepsilon\in[0,1]$ (introduced in the main text), results are summarized in
Fig. S.5.

Figure S.5: Average gap ratio $\langle r\rangle$ as a function of position in
the spectrum $\varepsilon$. Left column (a,c,e): the original QSM (4). Right
column (b,d,f): the U(1)-preserving QSM (5). Upper panels (a,b) show localized
regime, panels (c,d) show data in the vicinity of the critical point while
panels (e,f) have the ergodic regime at the center and the localized regions
at the edges. Black dotted lines indicate limiting values for ergodic (GOE)
and localized (Poisson) systems.

The right column, panels (b,d,f), in Fig. S.5 shows data for the U(1)
symmetric QSM (5). For $\alpha=0.7<\alpha_{c}$, all states are localized, as
shown in Fig S.5(b). In that case, with increasing system size, the mean gap
ratio decreases towards the Poisson limit regardless of the energy. In Fig.
S.5(d), the results at the critical point ${\alpha_{c}=0.76}$ are shown.
Indeed, we observe that the mean gap ratio exhibits almost no size dependence
for ${\varepsilon=0.5}$, i.e., at the center of the spectrum, while the rest
of the spectrum shows tendency towards localization.

The situation is, however, manifestly different in Fig. S.5(f), for
${\alpha=0.8>\alpha_{c}}$. Here, at the center of the spectrum, the larger the
system size, the bigger the $\langle r\rangle$. This clearly indicates that
the system becomes ergodic with increasing $L$. On the other hand, for either
low ${\varepsilon<0.3}$ or large ${\varepsilon>0.7}$ energies the trend is
reversed, the bigger the system, the closer level statistics is to the many-
body localized, Poisson limit. This is a clear manifestation of the many-body
mobility edge in the system. As shown in the left column, panels (a,c,e), of
Fig. S.5 the situation is similar for the original QSM, (4), with all states
being localized for small $\alpha$ [cf. Fig. S.5(a)], no size dependence for
$\alpha=\alpha_{c}$ at $\varepsilon=0.5$ [cf. Fig. S.5(c)], and a clear
indication of the mobility edge for $\alpha>\alpha_{c}$ [cf. Fig. S.5(e)].

Figure S.6: Extraction of intersection points for the U(1)-preserving QSM (5)
close to the edge of the spectrum (a) and the middle (b). The dash-dotted line
is interpolation of data points with polynomial of degree $3$ and the shaded
region corresponds to the $1\sigma$ confidence interval. The black horizontal
line segment corresponds to the position and uncertainty of the crossing
point.

An investigation of $\langle r\rangle$ vs $\alpha$ curves for different values
of $\varepsilon$ yields, qualitatively, the same conclusions as for the center
of the spectrum. To see whether those preliminary conclusions are correct, a
more systematic approach needs to be applied. To that end, we extract
positions of intersections of $\langle r\rangle$ vs $\alpha$ curves for
different $L$, namely the intersection between $L-1$ and $L+1$ curves for a
given $\varepsilon$, and we define position of this crossing to be system size
and energy dependent critical point $\alpha_{c}^{\varepsilon}(L)$. Position of
the intersection is obtained by fitting a polynomial of degree $3$ to the
corresponding $\langle r\rangle$ curves in the interval
$\alpha\in[\alpha_{*}^{\varepsilon}-0.05,\alpha_{*}^{\varepsilon}+0.05]$,
where $\alpha_{*}^{\varepsilon}$ is the estimate of the intersection point for
a given $\varepsilon$. Then $\alpha_{c}^{\varepsilon}(L)$ is extracted as the
intersection point of fitted polynomials, see Fig. S.6.

Figure S.7: Boundary between the phases in the QSM. (a) The original QSM (4).
(b) The U(1)-preserving QSM (5). Transition points extracted from the data are
shown with solid lines. Dotted lines show the fitted dependence on
$\varepsilon$, according to Eq. (S.5).

Figure S.8: Extrapolation of parameters $A_{0}$ and $A_{1}$ fitted to the
mobility edge shape according to Eq. (S.5) shown with a dotted line. The gray
band represents $1\sigma$ the confidence interval. The original QSM (4) shown
on (a,c), while the U(1)-preserving QSM (5) is on (b,d).

From the extracted values of $\alpha_{c}^{\varepsilon}(L)$ it is clear that
the critical point systematically decreases with increasing system size,
although this change is small. To obtain prediction for the critical point at
the thermodynamic limit, we considered the following two parameter formula:

$\alpha_{c}(\varepsilon)=A_{0}e^{\frac{|\varepsilon-\varepsilon_{0}|^{q}}{A_{1}}}\;,$
(S.5)

and we fitted it to the numerically obtained values of
$\alpha_{c}^{\varepsilon}(L)$ for both models. The exponent $q$ was chosen
based on the form of the density of states, which was perfectly Gaussian $q=2$
for the original model, while for the U(1) symmetric version, $q=2.25$ gave
the best description of the density of states. Perfect agreement of the data
with those functional forms is shown in Fig. S.7.

Parameters of the fit $A_{0},A_{1}$ are then extrapolated to the thermodynamic
limit, as shown in Fig. S.8. The scaling used was of the form
$A_{1,2}(L)=a_{1,2}+b_{1,2}/L$, as the parameters changed almost perfectly
linearly in $1/L$, where $a_{1,2},b_{1,2}$ are constants being fitted. This
scaling is also well motivated analytically, as we found that both parameters
$A_{0},A_{1}$ have the leading contribution of order $\mathcal{O}(L^{0})$,
thus the biggest subleading contribution should be $\mathcal{O}(L^{-1})$. We
also considered a scaling linear in the parameter $1/(N+L)$, which gave the
same results as $1/L$, thus we are not discussing it in detail. As the shape
of the mobility edge in the thermodynamic limit we took the formula (S.5) with
parameters $A_{0}(\infty),A_{1}(\infty)$, while its confidence intervals were
estimated using standard error propagation.

##  VI Entanglement entropy of eigenstates

In the main text, it was shown, that using entanglement entropy of eigenstates
as an indicator of ergodicity breaking gives results in agreement with gap
ratio analysis. At $\alpha>\alpha_{c}$ entropy $S$ followed volume law and the
rescaled entropy $S_{r}$ (defined in the main text) was increased towards
value $1$, as expected from the ergodic system. On the other hand,
$\alpha<\alpha_{c}$ showed area law and $S_{r}$ was decreasing, clearly
indicating localization.

Figure S.9: Entanglement entropy of eigenstates in QSM. (a,b) Average entropy,
$S$, far from the center of the spectrum ($\varepsilon=0.15$) in the original
QSM (4) on (a), and in the U(1) variant of QSM (5) on (b). (c) Average
entropy, $S$, close to the center of the spectrum ($\varepsilon=0.5$) in the
original QSM (4). (d) Dependence of the rescaled entanglement entropies
$S_{r}$ (defined in the main text) for the original QSM (4) at
$\varepsilon=0.5$. Inset shows the data for $\varepsilon=0.15$.

The decrease of entropy with increasing system size, being the consequence of
$0-$dimensional nature of the system, was also understood, as the effect of
going further from the GOE sun, which is the source of entropy growth in the
system.

Here, we supplement those results of the main text with the equivalent of Fig.
4 for the original QSM, see Fig. S.9(c,d). Moreover, we present area-volume
law transition for eigenstates far from the middle of the spectrum in Fig.
S.9(a,b) for both studied models. The numerical procedure to obtain those
plots is the same as described for entanglement entropies in the main text.
The results of Fig. S.9 show the same behavior as the spectral analysis, with
clear shift of the critical point $\alpha_{c}$ between area and volume laws
for different values of $\varepsilon$, giving further confirmation of the
mobility edge. This proves that results obtained are robust, independently of
the signature of the ergodicity breaking being studied.

##  VII Critical exponents and data collapse

When the phase transition of the system is approached, the characteristic
length $\xi$ diverges [112] and indicators of the transition $X(L,\alpha)$
follow the scaling:

$X(L,\alpha)=f[\mathrm{sgn}(\alpha-\tilde{\alpha})L^{1/\nu}/\xi],$ (S.6)

with a critical exponent $\nu$ and transition point $\tilde{\alpha}$ and
$f$-continuous function. For the QSM, the characteristic length is [67]:

$\xi=\frac{1}{2\left|\ln\frac{\alpha}{\tilde{\alpha}}\right|}.$ (S.7)

Figure S.10: Data collapse for (a) the original QSM for system sizes $6-14$
and (b) U(1) QSM for system sizes $7-16$. Critical parameters obtained from a
cost function method are shown. Black dotted lines indicate limiting values
for ergodic (GOE) and localized (Poisson) systems. (c) The critical exponent
of the transition at different rescaled energies $\varepsilon$ for both
models.

We applied the cost function method [55] to find parameters
$\nu,\tilde{\alpha}$ that provide optimal collapse (S.6) of $X=\langle
r\rangle$. To this mean, we sorted $\langle r\rangle$ with respect to
nondecreasing values of $\mathrm{sgn}(\alpha-\tilde{\alpha})L^{1/\nu}/\xi$,
obtaining an ordered list $\\{X_{i}\\}_{i=1}^{n_{\max}}$. The cost function is
then given as:

$\mathcal{C}_{X}=\frac{\sum_{i=1}^{n_{\max}-1}|X_{i+1}-X_{i}|}{\max{X_{i}}-\min{X_{i}}}-1.$
(S.8)

Optimal collapse was obtained when $\alpha$ values were restricted to
$\alpha\in[\alpha_{*}-0.05,\alpha_{*}+0.05]$, where $\alpha_{*}$ is the
estimate for the transition point obtained from prior analysis, results are
presented in Fig. S.10. Critical point $\tilde{\alpha}$ is in a very good
agreement with previously obtained $\alpha_{c}$.

Analogous analysis was also performed for different rescaled energies
$\varepsilon$, which is summarized in Fig. S.10c. We have found that critical
exponent $\nu$ does not vary across the spectrum for both models. This
confirms an intuitive expectation that the universality class of the
ergodicity breaking transition in the two models does not change with the
energy. It is important to note, that the value of $\nu$ depends significantly
on specific choice of $\alpha$ and system sizes $L$ included in the cost
function, the error associated to this effect is estimated as 5% of $\nu$.
Moreover, to calculate the statistical error of $\nu$ we repeated its
calculation for $1000$ bootstrap resamples, which were created by adding to
each data point $\langle r\rangle$ a random number drawn from a Gaussian
distribution with $0$ mean and standard deviation equal to the uncertainty of
this point. Uncertainties shown on Fig. S.10c are a statistical combination of
both types of errors. This data is consistent with the claim, that the
ergodicity breaking transition belongs to the same universality class in the
two models. This is further demonstrated by the fact that characteristic
length (S.7) derived for the original QSM leads to perfect collapse of the
data also for the U(1) model.