# A Quantum Breakdown Model: from Many-body Localization to Chaos with Scars

Biao Lian  Department of Physics, Princeton University, Princeton, New Jersey
08544, USA

###### Abstract

We propose a quantum model of fermions simulating the electrical breakdown
process of dielectrics. The model consists of $M$ sites with $N$ fermion modes
per site, and has a conserved charge $Q$. It has an on-site chemical potential
$\mu$ with disorder $W$, and an interaction of strength $J$ restricting each
fermion to excite two more fermions when moving forward by one site. We show
the $N=3$ model with disorder $W=0$ show a Hilbert space fragmentation in all
charge $Q$ sectors and is exactly solvable except for very few Krylov
subspaces. The analytical solution shows that the $N=3$ model exhibits many-
body localization (MBL) as $M\rightarrow\infty$, which is stable against $W>0$
as our exact diagonalization (ED) shows. At $N>3$, our ED suggests a MBL to
quantum chaos crossover as $M/N$ decreases across $1$. At $W=0$, an exactly
solvable many-body scar flat band exists in many charge $Q$ sectors, leading
to measure nonzero number of quantum scar eigenstates in the thermodynamic
limit. We further calculate the time evolution of a fermion added to the
particle vacuum, which shows the model is in a breakdown (dielectric) phase
when $\mu/J<1/2$ ($\mu/J>1/2$). The breakdown is local when $M/N\gg 1$, and is
global when $M/N\ll 1$.

##  I Introduction

The breakdown of dielectrics in an electric field (the Townsend avalanche) [1]
is one of the most violent far-from-equilibrium systems. Generically, it
involves a spatially asymmetric process which creates more and more particles
towards a fixed spatial direction, as driven by certain forces (electric
field, enthalpy gradient, etc). As a result, a microscopic perturbation can be
amplified into a macroscopic signal in a short time, producing a particle
avalanche. Therefore, the breakdown process is often utilized in designing
quantum measurement apparatuses, such as the Geiger Counter which works near
the breakdown electric field of inert gases [2, 3]. Similar processes also
occur in chemical reactions and generic chain reactions. However, the study of
the breakdown phenomenon has been largely restricted to phenomenological and
semiclassical models, since the quantum problem inevitably involves an
enormous Hilbert space and strong many-body interaction.

The studies in the past decades have significantly advanced our knowledge on
non-equilibrium quantum dynamics of many-body systems. The thermalization is
revealed to be closely related to many-body quantum chaos by the eigenstate
thermalization hypothesis [4, 5, 6, 7], and solvable quantum chaotic models,
particularly the Sachdev-Ye-Kitaev (SYK) models [8, 9, 10, 11], manifest the
quantum butterfly effect of perturbation amplifications. On the other hand,
quantum systems lacking thermalization have been investigated extensively,
which are known as the many-body localization (MBL) phase [12, 13, 14, 15,
16]. The MBL has also been numerically studied in interacting fermions in an
electric field (interacting Wannier-Stark) [17, 18] and similar models [19,
20, 21, 22]. However, the electrical breakdown has not yet been explored.
Models showing Hilbert space fragmentations have also been studied as systems
lacking thermalization [19, 20, 21, 22, 23, 24, 25]. In addition, many-body
quantum scars have been studied as non-thermal eigenstates in quantum chaotic
systems [26, 27, 28, 29, 30, 31, 32, 33]. These developments have provided new
tools for us to study the quantum physics of the electrical breakdown.

In this paper, we introduce a _quantum breakdown model_ with a _spatially
asymmetric_ interaction, to give a simplified quantum modeling of the
breakdown process. The model consists of a 1D chain of $M$ sites, and each
site consists of $N$ fermion degrees of freedom. Except for an on-site
chemical potential, the model only contains an asymmetric nearest-neighbor
interaction that annihilates (creates) a fermion on the $m$-th site and
creates (annihilates) three fermions on the $(m+1)$-th site ($1\leq m<M$). The
model is to provide a simplified model for the dielectric gas breakdown in an
electric field (the Townsend avalanche [1]), where each site represents a
layer of the gas perpendicular to the electric field, and the layer separation
is roughly the mean free path. When an atom (or molecule) in a certain layer
is ionized into an ion and electron, the electron will be accelerated by the
electric field and excite more ionizations (two more in the model here) in the
next layer through collisions, as shown in Fig. 1(a). Therefore, although the
microscopic interaction in the gas is generically inversion symmetric, the
effective interaction in the electric field is inversion asymmetric.

We note that our model here is an oversimplified model for breakdown in the
following sense: first, the electric field (or other generalized force) is not
explicitly present, instead its effect is integrated into the asymmetric
interaction. Secondly, we have ignored the spatial structure within each layer
of gas, assuming an electron in one layer can excite any electrons in the next
layer. Despite these simplifications, we can approximately identify the
asymmetric interaction strength $J$ as the voltage difference between adjacent
layers, and view the chemical potential $\mu$ in each layer as the ionization
energy of an atom. We also include a disorder strength $W$ giving fluctuations
in $\mu$.

The quantum breakdown model exhibits a surprisingly rich physics. The model
has a conserved charge $Q$, which significantly reduces its numerical
complexity. We first show that $N=3$ model with disorder strength $W=0$ show a
Hilbert space fragmentation into exponentially large number of Krylov
subspaces in all charge $Q$ sectors, and is exactly solvable except for very
few Krylov subspaces of certain charge $Q$ sectors. The exact solution shows
the $N=3$ model is in a 1D MBL phase in the $M\rightarrow\infty$ limit (the 1D
limit), and we further verified the robustness of the MBL under nonzero
disorder $W>0$ via exact diagonalization (ED). For $N>3$, we perform ED
calculations, and from the level spacing statistics (LSS) and eigenstate
entanglement entropies, we find the model shows MBL features when $M/N>1$,
while is quantum chaotic when $M/N<1$. This lead us to the conjecture that the
model is in a 1D MBL phase with localization length of order $N$ in the limit
$M/N\rightarrow\infty$, and is quantum chaotic in the limit $M/N\rightarrow 0$
(the 3D limit), with a crossover around $M/N\approx 1$. Moreover, when $W=0$,
in many charge $Q$ sectors, there exists a set of exactly solvable degenerate
eigenstates forming a many-body scar flat band, leading to quantum scar states
with a nonzero measure in the thermodynamic limit.

We further investigate the time evolution of a fermion added to the first
layer of the model, to examine the breakdown phenomena. Irrespective of $M/N$,
our calculations show the system is in a breakdown phase when $\mu/J<1/2$
($\mu\geq 0$), and in a dielectric phase (no breakdown) when $\mu/J>1/2$. The
ratio $M/N$ affects the spatial ranges of the breakdown. When $M/N>1$, the
model can only have a local breakdown with the particles ionized within order
$N$ number of sites after the initial perturbation. When $M/N<1$, the model
will exhibit a global breakdown with particles ionized in the entire system.

The rest of this paper is organized as follows. The quantum breakdown model
and its conserved charge is defined in Sec. II. In Secs. III and IV, we first
give the exact and almost exact solution of the $N=3$ model at disorder $W=0$,
and then compare with the ED results at both $W=0$ and $W>0$. Sec. V shows the
ED calculations for $N>3$ models suggesting the MBL to chaos crossover around
$M/N\approx 1$, and gives the theoretical understanding of the many-body scar
flat bands. Sec. VI presents the time-evolution of a fermion added to the
particle vacuum of the model, which reveals the breakdown transition at
$\mu/J=1/2$. Lastly, remarks and future directions are discussed in Sec. VII.

##  II The quantum breakdown model

###  II.1 The Hamiltonian

Fig. 1(a) gives a simplified illustration the breakdown (Townsend avalanche)
process of a dielectric gas in a large enough electric field (to the left),
where we assume the gas consists of neutral atoms (the black dots). When an
atom is ionized into an electron and an ion, the electron will be accelerated
in the direction opposite to the electric field, and then hit and ionize more
atoms, which will eventually lead to a particle Townsend avalanche as
illustrated by the blue lines. In principle, the ions can induce a similar
particle avalanche, but it would happen in a much longer time scale because of
their much heavier mass, and may be terminated by the recombinations of ions
with electrons from the electrode on the left. Moreover, if the dielectric is
a solid instead of a gas, the ions will stay static. Therefore, we will ignore
the ions and focus on the electrons to simplify the quantum model of the
breakdown.

By ignoring ions, we can assume each ionized electron in layer $m$ can hit
layer $m+1$ and create more ionized electrons via interaction. The cost of
ignoring the ions is that one can only write down interactions with electron
number changing by an even number, since electrons are fermions. Therefore,
the simplest breakdown interaction is an inversion asymmetric interaction
annihilating one electron in layer $m$ and creating three electrons in layer
$m+1$ (and its Hermitian conjugate). The distance between two adjacent layers
is thus approximately $2\ell_{\text{mfp}}$ for 2 collisions to happen, where
$\ell_{\text{mfp}}$ is the mean free path of the electrons.

Figure 1: (a) Simplified illustration of the electrical breakdown, where an
electron in a layer hits the next layer and creates two more electrons,
leading to a particle avalanche. (b) The quantum breakdown model in a 1D
lattice with $M$ sites (open boundary), and $N$ fermion modes per site. The
only interaction annihilates a fermion in site $m$ and creates three fermions
in site $m+1$ (and the inverse process).

To further simplify the model, we ignore the spatial structure within each
layer, and treat each layer as a site of a 1D lattice with an open boundary
condition, as shown in Fig. 1(b). By assuming $M$ sites (layers) in total and
$N$ complex fermion (electron) modes in each site, we define the quantum
breakdown model Hamiltonian in the 1D lattice as containing three parts:

$H=H_{I}+H_{\mu}+H_{\text{dis}}\ .$ (1)

The first part is the asymmetric breakdown interaction taking the form

$H_{I}=\sum_{m=1}^{M-1}\sum_{i<j<k}^{N}\sum_{l=1}^{N}\left(J_{m,l}^{ijk}c^{\dagger}_{m+1,i}c^{\dagger}_{m+1,j}c^{\dagger}_{m+1,k}c_{m,l}+h.c.\right),$
(2)

where $c_{m,i}$ and $c_{m,i}^{\dagger}$ are the annihilation and creation
operators of the $i$-th fermion mode on the $m$-th site ($1\leq i\leq N$,
$1\leq m\leq M$), and $J_{m,l}^{ijk}$ are complex parameters which are totally
antisymmetric in indices $i,j,k$. We shall consider the cases with both
uniform and random interaction $J_{m,l}^{ijk}$. For $H_{I}$ not to vanish, we
require the number of fermion modes per site $N\geq 3$.

The second part of the Hamiltonian is an on-site chemical potential uniform
within each site:

$H_{\mu}=\sum_{m=1}^{M}\mu_{m}\hat{n}_{m}\
,\qquad\hat{n}_{m}=\sum_{i=1}^{N}c^{\dagger}_{m,i}c_{m,i}\ ,$ (3)

where $\hat{n}_{m}$ is the fermion number operator on site $m$, and to be
general, we assume the chemical potential $\mu_{m}$ can depend on site $m$.

The last part is a disorder potential within each site:

$H_{\text{dis}}=\sum_{m=1}^{M}\sum_{i}^{N}\nu_{m,i}c^{\dagger}_{m,i}c_{m,i}\
,$ (4)

where $\nu_{m,i}$ are Gaussian random potentials with mean value and variance
given by

$\langle\nu_{m,i}\rangle=0\ ,\qquad\langle\nu_{m,i}^{2}\rangle=W^{2}\ .$ (5)

The constant $W\geq 0$ characterizes the disorder potential strength. We note
that one could add a non-diagonal random hopping term
$\nu_{m,ij}c^{\dagger}_{m,i}c_{m,j}$ within a site, but one could always
rewrite the model in the diagonal fermion basis of the hopping matrix
$\nu_{m,ij}$ on each site $m$, after which the disorder potential will take
the diagonal form of Eq. (4) without loss of generality.

For later convenience, we define the interaction strength between sites $m$
and $m+1$ as

$J_{m}=\left(\frac{1}{N}\sum_{i<j<k}^{N}\sum_{l=1}^{N}\left|J_{m,l}^{ijk}\right|^{2}\right)^{1/2}\
.$ (6)

In particular, when the system is on average uniform, which is the case we
will focus on dominantly, we can assume both the chemical potential and the
interaction strength are site $m$ independent:

$\mu_{m}=\mu\ ,\qquad J_{m}=J\ .$ (7)

In the case of fully random interactions $J_{m,l}^{ijk}$, we assume the
interactions are generated from a Gaussian distribution with statistical mean
and variance

$\langle J_{m,l}^{ijk}\rangle=0\
,\quad\langle\left|J_{m,l}^{ijk}\right|^{2}\rangle=\frac{3!J^{2}}{N(N-1)(N-2)}\
,$ (8)

where $J$ is defined in consistency with Eq. (7). We note the similarity of
such an interaction to that in the SYK model [8, 9, 10, 11].

###  II.2 Conserved charge

With the open boundary condition as shown in Fig. 1(b), the quantum breakdown
model in Eq. (1) has a conserved global U(1) charge

$Q=\sum_{m=1}^{M}3^{M-m}\hat{n}_{m}\ .$ (9)

In other words, each fermion in layer $m$ carries an effective charge
$q_{m}=3^{M-m}$. The breakdown process can thus be understood as charged
particles keeping splitting into partons carrying a $1/3$ factor smaller
charge as the layer index $m$ increases by one. This allows us to study a
fixed charge $Q$ sector each time, which has a much smaller Hilbert space
dimension than that of the entire model.

We note that one could complicate the Hamiltonian of Eq. (1) by adding other
terms, such as the hopping between different sites. This will, however, lead
to the loss of the conserved charge $Q$, making the analytical and numerical
studies much more difficult. In this paper, we restrict ourselves to the model
in Eq. (1) with the conserved charge $Q$.

###  II.3 Chiral symmetry

We can define a unitary chiral transformation $C$ as

$C(c_{m,l},c_{m,l}^{\dagger})C^{-1}=(-1)^{m}(c_{m,l},c_{m,l}^{\dagger})\ ,$
(10)

and it is straightforward to see that the model Hamiltonian
$H=H_{I}+H_{\mu}+H_{\text{dis}}$ transforms as

$C(H_{I}+H_{\mu}+H_{\text{dis}})C^{-1}=-H_{I}+H_{\mu}+H_{\text{dis}}\ .$ (11)

Moreover, $[C,Q]=0$. Therefore, when the disorder potential $W=0$, namely
$H_{\text{dis}}=0$, the Hamiltonian $H$ at chemical potentials $\mu$ and
$-\mu$ are exactly opposite to each other via a chiral transformation, and
thus have exactly opposite energy spectra in the same charge $Q$ sector. We
call this property the _chiral symmetry_ between $\mu$ and $-\mu$. When $W>0$,
the chiral symmetry still approximately true due to randomness of
$H_{\text{dis}}$. Therefore, we will only focus on $\mu\geq 0$ hereafter.

###  II.4 Physical meaning of the parameters

In a realistic breakdown system, both the number of fermion modes $N$ per
layer and the number of layers $M$ are large. We can view the ratio $M/N$ as
proportional to the ratio of the length (in the direction of electric field)
and cross sectional area of the space containing the dielectric gas. The ratio
$M/N$ is thus tunable by the geometry of the container. In the 1D limit where
the container is a quasi-1D long tube with a small cross sectional area, we
have $M/N\gg 1$. In contrast, in the 3D limit where the container has similar
sizes in all the three directions, we have $M/N\ll 1$.

Physically, the interaction strength $J$ in Eq. (7) is approximately the
average energy a fermion can gain by turning into three fermions in the next
site. Therefore, one can identify $J$ with the electric potential energy
difference $2eE_{\text{ele}}\ell_{\text{mfp}}$ of an ionized electron between
neighboring layers, where $e$ is the electron charge, and $E_{\text{ele}}$ is
the electric field. The chemical potential $\mu\geq 0$ can be understood as
the ionization energy for creating an electron, with a disorder strength $W$
due to fluctuations, which will disfavor breakdown. One thus expects the
quantum dynamics of the system to undergo a dielectric breakdown transition as
the ratio $\mu/J$ decreases towards zero, which is studied in Sec. VI.

##  III The $N=3$ model with $W=0$

We first study the model with $N=3$, the minimal number of fermion modes per
site for the interaction in Eq. (2) not to vanish. In this section, we fix the
disorder potential strength $W=0$ (defined in Eq. (5)), namely,
$H_{\text{dis}}=0$. We will show that this $N=3$ model with $W=0$ shows a
Hilbert space fragmentation into an exponentially large number of Krylov
subspaces in each charge $Q$ sector, and is exactly solvable except for very
few Krylov subspaces in certain charge $Q$ sectors, regardless of the
parameters $\mu_{m}$ and $J_{m,l}^{ijk}$. The resulting many-body eigenstates
indicate that the model is in a 1D MBL phase as $M\rightarrow\infty$, and the
energy spectrum exhibits an interesting fractal-like structure. The $N=3$
model with disorder strength $W>0$ will be studied in Sec. IV, where we find
the MBL phase remains robust.

###  III.1 Model simplification

With disorder potential $W=0$, we are allowed to simplify the model
interactions by U(3) unitary rotations of the fermion basis $c_{m,l}$ in each
site. We define a new fermion basis $f_{m,l}$ by

$c_{m,l}=\sum_{l^{\prime}=1}^{3}U_{ll^{\prime}}^{(m)}f_{m,l^{\prime}}\ ,$ (12)

where the unitary matrices $U^{(m)}$ satisfy
$\sum_{l}J_{m,l}^{123}U_{ll^{\prime}}^{(m)}=\sqrt{3}J_{m}\delta_{1,l^{\prime}}\det\left(U^{(m+1)}\right)$,
with the interaction strength $J_{m}$ defined in Eq. (6), and
$\delta_{l,l^{\prime}}$ is the Kronecker delta function. Note that the product
operator in each site
$c_{m,1}c_{m,2}c_{m,3}=\det\left(U^{(m)}\right)f_{m,1}f_{m,2}f_{m,3}$ is
invariant up to a phase factor $\det\left(U^{(m)}\right)$. This transforms the
total Hamiltonian in Eq. (1) into the $f_{m,l}$ basis:

$\begin{split}&H=H_{I}+H_{\mu}\ ,\\\
&H_{I}=\sum_{m=1}^{M-1}\Big{[}\sqrt{3}J_{m}\Big{(}\prod_{j=1}^{3}f^{\dagger}_{m+1,j}\Big{)}f_{m,1}+h.c.\Big{]},\\\
&H_{\mu}=\sum_{m=1}^{M}\mu_{m}\hat{n}_{m}\ ,\end{split}$ (13)

where $\hat{n}_{m}=\sum_{l=1}^{3}f_{m,l}^{\dagger}f_{m,l}$ is equivalent to
the definition in Eq. (3). Note that only the first fermion mode $f_{m,1}$ in
each site can excite three fermions in the next site.

###  III.2 Eigenstates of the $Q=3^{M-1}$ charge sector

The Hamiltonian is block diagonal in each charge $Q$ sector, the Hilbert space
of which is spanned by the basis of all Fock states with integer fermion
number $n_{m}=\langle\hat{n}_{m}\rangle$ on each site $m$, satisfying
$Q=\sum_{m}3^{M-m}n_{m}$.

As a prototypical example, we investigate the $Q=3^{M-1}$ charge sector in
this subsection. This is the charge sector of adding one fermion to the first
site of the particle vacuum state $|0\rangle$ (vacuum state defined by
$\hat{n}_{m}|0\rangle=0$ for all $m$). The fermion numbers of the Fock state
basis are thus limited to the following string configurations:

$n_{m}=2\vartheta_{m_{s}-m+1}-2\delta_{m,1}+\delta_{m,m_{s}}\ ,$ (14)

where $m_{s}$ is an integer satifying $1\leq m_{s}\leq M$, and we have defined
the integer variable Heaviside function $\vartheta_{m}$ as

$\vartheta_{m}=\begin{cases}&1\ ,\quad(m>0)\\\ &0\ ,\quad(m\leq 0)\end{cases}\
,\qquad(m\in\mathbb{Z}).$ (15)

More explicitly, the fermion number configurations in Eq. 14 takes the form
$\\{n_{m}\\}=\\{0,2,\cdots,2,3,0,\cdots,0\\}$ when $m_{s}>1$, where the
$m_{s}$-th site has $3$ fermions and is the last site with a nonzero fermion
number (i.e., the site the fermions have spread to). When $m_{s}=1$, one has
$\\{n_{m}\\}=\\{1,0,\cdots,0\\}$, and we call this configuration the
_reference configuration_ of charge $Q=3^{M-1}$, which is the configuration
with the _smallest_ total number of fermions possible in this charge $Q$
sector.

The Fock states with $n_{m}$ given in Eq. (14) with all possible $m_{s}$
constitute the Hilbert space basis of the charge $Q=3^{M-1}$ sector, which has
a Hilbert space dimension $h_{Q}=\frac{3^{M-1}+5}{2}$. In the below, we will
show that this Hilbert space further fragments into an exponentially large
number of disjoint Krylov subspaces, each of which has its Hamiltonian
equivalent to a single-particle tight-binding model.

####  III.2.1 A Krylov subspace without blocking

A _Krylov subspace_ is defined as a closed Hilbert space spanned by states
$H^{n}|\psi\rangle$ of all nonnegative integers $n\geq 0$ from some root state
$|\psi\rangle$, where $H$ is the Hamiltonian. We first examine a special
Krylov subspace in the $Q=3^{M-1}$ charge sector, which has a basis of $M$
Fock states defined by

$|m_{s},1,1\rangle=f_{m_{s},1}^{\dagger}\prod_{m=2}^{m_{s}}f^{\dagger}_{m,2}f^{\dagger}_{m,3}|0\rangle\
,$ (16)

where $1\leq m_{s}\leq M$. Note that the first fermion mode $f_{m,1}$ on each
site $m$ is empty except for the $m_{s}$-th site. The 3-quantum number
notation of state $|m_{s},1,1\rangle$ will become clear in next subsection.
These $M$ Fock states form a closed Hilbert Krylov subspace, in which we
choose $|1,1,1\rangle$ to be the root state, and the Hamiltonian in Eq. (13)
acts in this Krylov subspace as

$\begin{split}H&|m_{s},1,1\rangle=\vartheta_{m_{s}-1}\sqrt{3}J_{m_{s}-1}|m_{s}-1,1,1\rangle\\\
&+\vartheta_{M-m_{s}}\sqrt{3}J_{m_{s}}|m_{s}+1,1,1\rangle+V_{m_{s}}|m_{s},1,1\rangle.\end{split}$
(17)

Here $V_{m_{s}}=\mu_{m_{s}}+2\sum_{m=2}^{m_{s}}\mu_{m}$, the coefficients
$J_{m}$ are given in Eq. 6, and the notation $\vartheta_{m}$ is defined in Eq.
(15). Such a Hamiltonian is thus equivalent to a single-particle tight-binding
Hamiltonian on sites $1\leq m_{s}\leq M$ with an open boundary condition,
which has an on-site potential $V_{m_{s}}$, and a nearest neighbor hopping
$\sqrt{3}J_{m_{s}}$. The many-body state $|m_{s},1,1\rangle$ is thus
equivalent to a single free effective particle on site $m_{s}$. We say this
subspace has no blocking (also defined in the next subsection), since all the
sites are connected by hopping.

It is straightforward to diagonalize such a single-particle Hamiltonian and
obtain the eigenstates in this subspace. Generically, if the parameters
$V_{m_{s}}$ and/or $J_{m_{s}}$ are disordered with respect to $m_{s}$, the
tight-binding Hamiltonian in Eq. 17 will exhibit the _Anderson localization_
[34], with all the eigenstates spatially localized around different sites
$m_{s}$. Since the basis in Eq. 16 are many-body fermion states, these
eigenstates should be viewed as many-body localized (MBL) in the original
fermion basis.

Figure 2: Illustration of the eigenstates of the $N=3$ model $W=0$. (a) The
Krylov subspace without blocking in the charge $Q=3^{M-1}$ sector, which is
conditionally extended. (b) Krylov subspaces with blocking site $m_{b}$ in the
charge $Q=3^{M-1}$ sector, which is localized on sites $m\leq m_{b}$ and
conditionally extended on sites $m>m_{b}$. (c) In a Krylov subspace of a
generic charge sector $Q$, the lattice is divided into alternating localized
and conditionally extended intervals, leading to MBL eigenstates.

In the case when Eq. 7 is satisfied, namely, when the chemical potential
$\mu_{m}=\mu$ and interaction strength $J_{m}=J$ are uniform, one has
$V_{m}=(2m-1)\mu$ that increases linearly with $m$. If $\mu\neq 0$, the
Hamiltonian of Eq. 17 resembles that of the non-interacting electron Bloch
oscillation in electric field, and thus the eigenstates will be _Wannier-Stark
localized_ [35] within roughly a localization length of order $J/|\mu|$.

In the special case when

$\mu_{m}=\frac{\mu_{1}}{3^{m-1}}\ ,\qquad J_{m}=J\ ,$ (18)

one has $V_{m_{s}}=\mu_{1}$ and $J_{m_{s}}=J$ independent of $m_{s}$. As a
result, the eigenstates will be extended in the entire 1D chain, given by:

$|\psi_{k,1,1}\rangle=\sum_{m_{s}=1}^{M}\sin\frac{\pi
km_{s}}{M+1}|m_{s},1,1\rangle\ ,$ (19)

where $k$ is an integer ($1\leq k\leq M$). The according eigenenergy is

$E_{k,1,1}=\mu_{1}+2\sqrt{3}J\cos\frac{\pi k}{M+1}\ .$ (20)

More generally, the eigenstates will be extended if $V_{m}$ and $J_{m}$ are
periodic in $m$.

Therefore, we say the eigenstates in this Krylov subspace of Eq. (16) are
_conditionally extended_ (on the hollow sites in Fig. 2(a)), meaning that they
are extended if and only if $V_{m}$ and $J_{m}$ are periodic in $m$. The
states in this Krylov subspace behave as a single free effective particle in
the entire lattice.

####  III.2.2 Other Krylov subspaces: blocking enforced localization

Following a similar idea, we can determine all the other Krylov subspaces in
the charge $Q=3^{M-1}$ sector and derive the eigenstates exactly.

The Fock basis in Eq. 16 in the Krylov subspace above has the first fermion
mode $f_{m,1}$ empty in all sites but site $m_{s}$. We now consider a set of
more generic Fock states defined as follows:

$\begin{split}&|m_{s},m_{b},\\{s_{m}^{j}\\}\rangle\\\
=&f_{m_{s},1}^{\dagger}\prod_{m=m_{b}+1}^{m_{s}}f^{\dagger}_{m,2}f^{\dagger}_{m,3}\prod_{m=2}^{m_{b}}f^{\dagger}_{m,s^{1}_{m}}f^{\dagger}_{m,s^{2}_{m}}|0\rangle,\end{split}$
(21)

where $2\leq m_{b}<m_{s}\leq M$, and $\\{s_{m}^{j}\\}$ are the indices of
occupied fermion modes on sites $2\leq m\leq m_{b}$ ($j=1,2$), with the
specific requirement on site $m_{b}$:

$s_{m_{b}}^{1}\equiv 1\ ,\qquad s_{m_{b}}^{2}=2\text{ or }3\ .$ (22)

For all the other sites $2\leq m<m_{b}$, the coefficients $s_{m}^{1}$ and
$s_{m}^{2}$ can take any two distinct values among $\\{1,2,3\\}$. Therefore,
$m_{b}$ and $m_{s}$ are the last two sites with the first fermion mode
$f_{m,1}$ occupied.

We call the site $m_{b}$ in Eq. 21 a _blocking site_ , for the reasons below.
Since the first fermion mode $f_{m_{b},1}$ on site $m_{b}$ is occupied, and
the $(m_{b}+1)$-th site is not empty, it is easy to see that the Hamiltonian
in Eq. 13 cannot change the occupation configuration on sites $1\leq m\leq
m_{b}$ by any number of actions on the Fock state in Eq. 21. Namely, the
coefficients $\\{s_{m}^{j}\\}$ are invariant. Therefore, the dynamics of the
Fock state in Eq. 21 is blocked at site $m_{b}$ and only happens on sites
$m>m_{b}$.

As a result, one can show that all the Fock states
$|m_{s},m_{b},\\{s_{m}^{j}\\}\rangle$ in Eq. 21 with the same blocking site
$m_{b}$ and occupation configuration $\\{s^{j}_{m}\\}$ form a closed Krylov
subspace, and we define the root state of this Krylov subspace to be the Fock
state with $m_{s}=m_{b}+1$. So we can simply label the Krylov subspaces by
quantum numbers $(m_{b},\\{s_{m}^{j}\\})$. The Hamiltonian in Eq. 13 acting on
Fock states in such a Krylov subspace in a similar way to Eq. 17, but only on
sites $m_{s}>m_{b}$:

$\begin{split}H&|m_{s},m_{b},\\{s_{m}^{j}\\}\rangle=V_{m_{s}}|m_{s},m_{b},\\{s_{m}^{j}\\}\rangle\\\
&\quad+\vartheta_{m_{s}-m_{b}-1}\sqrt{3}J_{m_{s}-1}|m_{s}-1,m_{b},\\{s_{m}^{j}\\}\rangle\\\
&\quad+\vartheta_{M-m_{s}}\sqrt{3}J_{m_{s}}|m_{s}+1,m_{b},\\{s_{m}^{j}\\}\rangle\
,\end{split}$ (23)

where as before, $V_{m_{s}}=\mu_{m_{s}}+2\sum_{1<m\leq m_{s}}\mu_{m}$, $J_{m}$
is defined in Eq. 6, and $\vartheta_{m}$ is given in Eq. 15. Therefore, in the
Krylov subspace of Fock states in Eq. 21 with fixed indices $m_{b}$ and
$\\{s_{m}^{j}\\}$, the Hamiltonian behaves as a single-particle tight-binding
model on sites $m_{b}+1\leq m_{s}\leq M$ with an open boundary condition, and
the states can be viewed as a single free effective particle in this interval.
In particular, the energy spectrum of Eq. 23 for a fixed blocking site $m_{b}$
is obviously independent of the indices $\\{s_{m}^{j}\\}$. Therefore, each
energy level with a blocking site $m_{b}$ in the charge $Q=3^{M-1}$ sector has
a degeneracy $2\cdot 3^{m_{b}-2}$, which is the number of possible index
configurations $\\{s_{m}^{j}\\}$, or the number of Krylov subspaces with the
same blocking site $m_{b}$.

The eigenstates of Eq. 23 in a Krylov subspace with blocking site $m_{b}$ are
thus conditionally extended on sites $m_{b}<m\leq M$, as illustrated in Fig.
2(b). If $V_{m_{s}}$ and $J_{m_{s}}$ are not periodic in $m_{s}$, one expects
the eigenstates to be Anderson localized (or Wannier-Stark localized if
$V_{m_{s}}$ is linear in $m_{s}$), or MBL in the fermion basis perspective. If
instead $V_{m_{s}}$ and $J_{m_{s}}$ are periodic, for instance if Eq. 18 is
satisfied, yielding $V_{m_{s}}=\mu_{1}$ and $J_{m_{s}}=J$, the eigenstates
will be extended within $m_{b}<m\leq M$, with wavefunctions

$|\psi_{k,m_{b},\\{s^{j}_{m}\\}}\rangle=\sum_{m_{s}=m_{b}+1}^{M}\sin\frac{\pi
k(m_{s}-m_{b})}{M-m_{b}+1}|m_{s},m_{b},\\{s_{m}^{j}\\}\rangle,$ (24)

where $k$ is an integer ranging within $1\leq k\leq M-m_{b}$. The energy of
this eigenstate is given by

$E_{k,m_{b},\\{s^{j}_{m}\\}}=\mu_{1}+2\sqrt{3}J\cos\frac{\pi k}{M-m_{b}+1}\ ,$
(25)

which is degenerate with respect to all the $2\cdot 3^{m_{b}-2}$ occupation
configurations $\\{s^{j}_{m}\\}$, as we have explained.

In addition to the above eigenstates in Krylov spaces with blocking site
$m_{b}\geq 2$, there are two Fock states decoupled from all the other states:

$|1,1,s\rangle=f^{\dagger}_{1,s}|0\rangle\ ,\quad(s=2,3)\ ,$ (26)

which are readily the eigenstates of the Hamiltonian in Eq. 13. Their energy
eigenvalues are simply

$E_{1,1,s}=\mu_{1}\ ,\qquad(s=2,3).$ (27)

We can view these two states as two one-dimensional Krylov subspaces with a
blocking site $m_{b}=1$. In this way, we can summarize the Fock states in Eqs.
16, 21 and 26 in the same notation $|m_{s},m_{b},\\{s_{m}^{j}\\}\rangle$, with
$1\leq m_{b}\leq m_{s}\leq M$, and $(m_{b},\\{s_{m}^{j}\\})$ label the Krylov
subspace. If $m_{b}\geq 2$, only those states with $m_{b}<m_{s}$ are valid.

One can verify that the total number of eigenstates we have found adds up to
the Hilbert space dimension $h_{Q}=\frac{3^{M-1}+5}{2}$ of the charge sector
$Q=3^{M-1}$.

####  III.2.3 Numerical spectrum and entanglement entropy

To verify the above exact solutions in the $Q=3^{M-1}$ sector, we perform
exact diagonalization (ED) calculations of the energy spectrum of the $N=3$
quantum breakdown model. For generality, our ED is performed in a full charge
$Q$ sector in the original fermion basis $c_{m,i}$, assuming not knowing any
Krylov subspaces.

We first show an example with $(N,M)=(3,10)$ and no randomness: we assume all
the interactions in Eq. 2 are uniformly $J^{123}_{m,l}=10$, the chemical
potential $\mu_{m}=\mu=0$, and the disorder potential $W=0$. The many-body
energy spectrum in the $Q=3^{M-1}=3^{9}$ sector is shown in Fig. 3(a), where
$\alpha$ is an integer labeling the energy levels $E(\alpha)$ sorted from low
to high ($1\leq\alpha\leq h_{Q}$). The energy levels agree well with Eqs. 20,
25 and 27. Not surprisingly, most energy levels have an extensive degeneracy,
making the spectrum fractal-like. This is consistent with the fact that the
energy levels with a blocking site $m_{b}$ have a degeneracy $2\cdot
3^{m_{b}-2}$ (the number of degenerate Krylov subspaces with the same $m_{b}$)
plus other accidental degeneracies, and the fractal-like spectrum is because
the degeneracies are proportional to powers of $3$. Moreover, we note that the
parameters of Fig. 3(a) satisfy Eq. 18 with $\mu_{1}=0$ and $J=10$, so the
eigenstates are extended on sites $m>m_{b}$.

Another example we show in the $(N,M)=(3,10)$ case has random complex
interactions $J^{ijk}_{m,l}$ satisfying Eq. 8, where the interaction strength
$J=10$, while the chemical potential $\mu_{m}=\mu=1$ is uniform, and $W=0$.
Fig. 3(c) shows the energy spectrum in the $Q=3^{M-1}=3^{9}$ sector. Clearly,
the fractal-like level degeneracies remain for random interactions, agreeing
well with the degeneracy $2\cdot 3^{m_{b}-2}$ for blocking site $m_{b}$. The
randomness of interaction does remove all accidental degeneracies.

Figure 3: ED calculation in a charge $Q$ sector of the model with
$(N,M)=(3,10)$ and disorder potential $W=0$. (a), (c), (e) shows the energy
levels $E(\alpha)$ sorted from low to high (labeled by $\alpha$), and the
corresponding LSS in the insets. (b), (d), (f) show the eigenstate
entanglement entropy in subregion $A$ (the first $\lfloor M/2\rfloor$ sites)
versus energy $E(\alpha)$ for (a), (c), (e), respectively (the darker the
higher density of dots). The parameters are labeled in the left panels. In
(a)-(b), the interaction $J^{123}_{m,l}=10$ is uniform, and $Q=3^{M-1}$. In
(c)-(f), the interaction $J^{123}_{m,l}$ is random with strength $J=10$, while
$Q=3^{M-1}$ in (c)-(d) and $Q=21924$ in (e)-(f).

The level spacing statistics (LSS) in a charge sector is known as an indicator
of integrability or quantum chaos, which is defined as the statistical
distribution $p(\delta_{E})$ of the nearest neighboring level spacing
$\delta_{E}(\alpha)=E(\alpha+1)-E(\alpha)$ ($1\leq\alpha\leq h_{Q}-1$).
Integrable systems show Poisson distribution LSS [36]

$p(\delta_{E})\propto e^{-\delta_{E}/\lambda_{0}}\ ,$ (28)

with some constant $\lambda_{0}$, while fully quantum chaotic systems show
Wigner-Dyson (WD) distribution LLS [37, 38, 39]

$p(\delta_{E})\propto\delta_{E}^{n}e^{-\delta_{E}^{2}/\lambda_{0}^{2}}\ ,$
(29)

where $n=1,2,4$ for models in the orthogonal, unitary and symplectic symmetry
classes, respectively.

The insets of Fig. 3(a) and (c) show the LSS in the charge $Q=3^{9}$ sector.
The extensive level degeneracies yield a delta function peak in the LSS at
$\delta_{E}=0$, while at $\delta_{E}>0$, the LSS shows a Poisson distribution
(although $\delta_{E}$ is large compared to $1/h_{Q}$). This is consistent
with the fact that all the eigenstates in the $Q=3^{M-1}$ sector are exactly
solvable, namely, integrable.

To tell how localized the eigenstates are, we calculate the entanglement
entropy

$S_{A}(\alpha)=-\text{tr}[\rho_{A}(\alpha)\ln\rho_{A}(\alpha)]$ (30)

of the reduced density matrix $\rho_{A}(\alpha)$ of each eigenstate
$|\alpha\rangle$ in a subregion $A$. The subregion $A$ is always defined as
sites $1\leq m\leq\lfloor M/2\rfloor$ hereafter, where $\lfloor x\rfloor$ is
the floor function giving the integer part of $x$. MBL and extended
eigenstates should exhibit area law and volume law entanglement entropies,
respectively. Fig. 3(b) and (d) show the entanglement entropy $S_{A}(\alpha)$
of the $Q=3^{M-1}$ sector eigenstates in Fig. 3(a) and (c), respectively. Each
dot represents an eigenstate, and the darker the color of the dots is, the
denser the dots are (to reflect the overlapping dots). In both cases, there
are a large number of eigenstates with low entanglement entropy irrespective
of the eigen-energy $E(\alpha)$, which can be identified as area-law MBL
eigenstates. The rest states, however, have higher entanglement entropy peaked
in the middle of the energy spectrum, showing a volume law. The presence of
volume law states has two reasons: first, as shown in Fig. 2, many eigenstates
have conditionally extended regions, which may show a relatively large
entanglement. In particular, compared to most other charge sectors, the
$Q=3^{M-1}$ has a relatively larger conditionally extended region (see Sec.
III.3 below), thus is more delocalized. In Fig. 3(b) where the extended
condition Eq. 18 is satisfied, the eigenstates are like an extended free
effective particle in the conditionally extended region. In Fig. 3(d) where
the chemical potential $\mu_{m}=\mu=1$, one expects Wannier-Stark localization
with localization length $\sim J/|\mu|$, thus entanglement entropies are
generically lower than those in Fig. 3(b). Second, although the eigenstates in
each Krylov subspace we give in Sec. III.2.2 are strictly localized on sites
$m<m_{b}$, one can linearly superpose the extensively degenerate eigenstates
from different Krylov subspaces into a new eigenbasis with high entanglement
entropies. This is generically the case in ED which ends up with a random
degenerate eigenbasis. Intriguingly, adding a small disorder potential $W$
(Eq. 5) is sufficient to slightly break the degeneracy and fix the eigenbasis
to a low-entanglement basis, as shown in Fig. 4(b), which has $W=10^{-6}$ and
the rest parameters the same as those in Fig. 3(b).

###  III.3 The generic charge sectors

We now study the Krylov subspaces and eigenstates in an arbitrary charge $Q$
sector of the $N=3$ model. For later purpose, we rewrite the charge $Q$ as

$Q=\sum_{m=1}^{M}n^{0}_{m}3^{M-m}\ ,\quad(n^{0}_{m}\in\\{0,1,2\\})\ ,$ (31)

with the restriction $0\leq n^{0}_{m}\leq 2$. The set of integers
$\\{n^{0}_{m}\\}$ has one-to-one unique correspondence with charge $Q$, which
is nothing but the digits of $Q$ in ternary number. For instance, the
$Q=3^{M-1}$ sector we studied in Sec. III.2 corresponds to
$\\{n^{0}_{m}\\}=\\{1,0,0,\cdots\\}$. We call $\\{n^{0}_{m}\\}$ the _reference
configuration_ of charge $Q$.

Similar to Eq. 14, the allowed fermion number configurations $\\{n_{m}\\}$ of
the on-site Fock basis are uniquely determined by charge $Q$. In particular,
the reference configuration $\\{n^{0}_{m}\\}$ of charge $Q$ in Eq. 31 is the
configuration with the smallest total number of fermions possible for charge
$Q$.

We assume there are $M_{Q}$ sites on which $n^{0}_{m}>0$ (i.e., $M_{Q}$ filled
sites in the reference configuration), and assume $m_{f,p}$ is the $p$-th such
filled site sorted from small to large ($1\leq p\leq M_{Q}$). We now show that
depending on the occupation numbers $n^{0}_{m_{f,p}}$ on these $M_{Q}$ sites,
the eigenstates in the charge $Q$ sector fractionalizes into an exponentially
large number of Krylov subspaces which are either exactly solvable or almost
solvable.

####  III.3.1 Fully exactly solvable charge sectors

If charge $Q$ in Eq. 31 has $n^{0}_{m_{f,p}}=2$ for any $2\leq p\leq M_{Q}$,
the eigenstates of all the Krylov subspaces in the charge $Q$ sector can be
exactly solved.

To see this, we first note that the $M_{Q}$ sites at positions $m_{f,p}$
separate the system into $M_{Q}+1$ decoupled subsystems. The first subsystem
is the interval of sites $1\leq m\leq m_{f,1}-1$ (which vanishes if
$m_{f,1}=1$), and stays trivially in the zero fermion state. The $(p+1)$-th
($1\leq p\leq M_{Q}$) subsystem is the interval of sites
$m\in[m_{f,p},m_{f,p+1}-1]$ (here we define $m_{f,M_{Q}+1}=M+1$). It is easy
to verify that the allowed Fock basis fermion number configurations within the
$(p+1)$-th subsystem ($1\leq p\leq M_{Q}$) are similar to Eq. 14, given by

$n_{m}=2\vartheta_{m_{s,p}-m+1}-(3-n^{0}_{m_{f,p}})\delta_{m,m_{f,p}}+\delta_{m,m_{s,p}}$
(32)

for sites $m\in[m_{f,p},m_{f,p+1}-1]$, where the integer
$m_{s,p}\in[m_{f,p},m_{f,p+1}-1]$ is a free parameter, and $\vartheta_{m}$ is
defined in Eq. 15. Namely, the allowed configurations within the $(p+1)$-th
subsystem typically take the form
$\\{n_{m}\\}=\\{n^{0}_{m_{f,p}}-1,2,\cdots,2,3,0,\cdots,0\\}$, where the site
with $3$ fermions is the $m_{s,p}$-th site. Note that the allowed charge
configuration in the $(p+1)$-th subsystem is solely determined by
$n^{0}_{m_{f,p}}$. This is because $n^{0}_{m_{f,p}}=2$ for any $p\geq 1$
guarantees that the $m_{f,p}$-th site is never empty, so the fermions in the
$p$-th subsystem is blocked and cannot spread into with the $(p+1)$-th
subsystem.

Because of the decoupling of the $M_{Q}+1$ subsystems, the many-body
eigenstates of the entire system decomposes into the direct products of
$M_{Q}+1$ subsystem eigenstates. Except for the first subsystem which is
trivially in the particle vacuum state, each subsystem has Krylov subspaces
and exactly solvable eigenstates by the same method as we used for the charge
$Q=3^{M-1}$ sector (see Sec. III.2). A Krylov subspace of the $(p+1)$-th
subsystem is labeled by a blocking site $m_{b,p}$, defined as the second last
site in the interval $[m_{f,p},m_{f,p+1}-1]$ with the first fermion mode
$f_{m_{b,p},1}$ occupied (similar to Eq. 21), and occupation indices
$\\{s_{m}^{j}\\}$ in the interval $[m_{f,p},m_{b,p}]$. Its root state is
defined as the Fock state with the smallest possible total fermion number in
this subspace. In this Krylov subspace, the Hamiltonian maps to a single-
particle tight-binding Hamiltonian similar to Eqs. 17 and 23 with $m_{s}$ and
$m_{b}$ replaced by $m_{s,p}$ and $m_{b,p}$ (see App. A), thus is exactly
solvable. The Krylov subspaces of the entire system are simply the direct
product of the Krylov subspaces of the $M_{Q}+1$ subsystems.

Fig. 2(c) illustrates the localization of a generic many-body eigenstate in a
Krylov subspace of such exactly solvable charge $Q$ sectors. Within the
$(p+1)$-th subsystem on sites $[m_{f,p},m_{f,p+1}-1]$, there is a blocking
site $m_{b,p}$. The subsystem is then frozen in the localized configuration
$\\{s_{m}^{j}\\}$ within the interval of sites $[m_{f,p},m_{b,p}]$, and maps
to a conditionally extended single free effective particle within sites
$[m_{b,p}+1,m_{f,p+1}-1]$, similar to the charge $Q=3^{M-1}$ sector (Fig.
2(b)). This separates the entire system into alternating localized and
conditionally extended intervals (Fig. 2(c)), which maps to $M_{Q}$ free
effective particles (one in each conditionally extended interval) decoupled
with each other. Such many-body eigenstates are therefore necessarily MBL in
the original fermion basis, with a localization length not exceeding the size
of each subsystem. In the $Q=3^{M-1}$ sector, there is only one subsystem and
some eigenstates can still be extended; while in most charge $Q$ sectors, the
number of subsystems $M_{Q}$ is usually a large number of the same order as
$M$, so the localization length is at most $M/M_{Q}$ which is small.

####  III.3.2 The other charge sectors: exactly solvable in most Krylov
subspaces

For a generic charge $Q$ sector where $n_{m}^{0}$ in Eq. 31 are not limited to
values $0$ and $2$, the above picture of $M_{Q}+1$ decoupled subsystems
separated by sites $m_{f,p}$ is no longer exact for all Krylov subspaces, but
is still a good approximation.

A Krylov subspace is still labeled by the blocking sites $\\{m_{b,p}\\}$ and
occupation indices $\\{s_{m}^{j}\\}$ in the $M_{Q}$ intervals
$[m_{f,p},m_{b,p}]$. We still define the $(p+1)$-th subsystem as the interval
$[m_{f,p},m_{f,p+1}-1]$. If one has $n_{m_{f,p}}^{0}=1$ on a site $m_{f,p}$
for some $p\geq 2$, there will be a weak interaction between the $p$-th
subsystem and the $(p+1)$-th subsystem on their boundary site $m_{f,p}$ in
certain Krylov subspaces. This is because $n_{m_{f,p}}^{0}=1$ allows site
$m_{f,p}$ to be either occupied or empty (by the actions of interaction
$H_{I}$ in Eq. 13). When site $m_{f,p}$ is empty, it can either gain $3$
fermions from the $p$-th subsystem on its left, or gain $1$ fermion from the
$(p+1)$-th subsystem on its right, which induces an interaction between the
two subsystems on site $m_{f,p}$. Note that the fermions in each of the two
subsystems cannot invade into each other beyond the boundary site $m_{f,p}$,
which restricts the interaction to be on the site $m_{f,p}$ only. However,
such an interaction is present only in the Krylov subspaces with the blocking
site $m_{b,p}=m_{f,p}$. This is because if $m_{b,p}>m_{f,p}$, the occupation
configuration in the interval $[m_{f,p}+1,m_{b,p}]$ will be fixed, and the
site $m_{f,p}$ can only either be empty or accept 3 fermions from the $p$-th
subsystem on its left, which eliminates the interaction. Therefore, there are
only very few Krylov subspaces having this inter-subsystem boundary
interaction, which must have $m_{b,p}=m_{f,p}$ for some $p\geq 2$. A more
detailed expression of such an inter-subsystem interaction can be found in
App. A.

If one view the Fock states in a Krylov subspace within each interval
$[m_{b,p},m_{f,p+1}]$ as a conditionally extended single effective particle on
site $m_{s,p}\in[m_{b,p},m_{f,p+1}]$ (the site with $3$ fermions) similar to
that in Eqs. 17 and 23, the above picture can be illustrated as follows. If
$n_{m_{f,p}}^{0}=1$ and $m_{b,p}=m_{f,p}$ for some $p\geq 2$, the $(p-1)$-th
effective particle and the $p$-th effective particle will have a hard core
repulsion on site $m_{f,p}$. Otherwise, the $(p-1)$-th effective particle and
the $p$-th effective particle are decoupled. Namely, among the $M_{Q}$
effective particles (one in each conditionally extended interval in Fig.
2(c)), two neighboring effective particles can have a boundary hard core
interaction with each other only when their conditionally extended intervals
touch each other and the boundary site $m_{f,p}$ between them has
$n_{m_{f,p}}^{0}=1$. Such an interaction obstructs us to analytically obtain
the eigenstates in these Krylov subspaces.

Since this interaction in any generic charge $Q$ sector only weakly exists on
the boundary between some neighboring subsystems in very few Krylov subspaces,
we expect the eigenstates to be well approximated by direct products of
eigenstates of $M_{Q}+1$ subsystems, which would still show MBL. As a
numerical test, in Fig. 3(e)-(f) we show the energy spectrum and subregion
entanglement entropy of the charge $Q=21924$ sector eigenstates of the
$(N,M)=(3,10)$ model calculated by ED, where the model parameters are the same
as Fig. 3(c)-(d), namely, uniform $\mu_{m}=\mu=1$, random interaction strength
$J=10$ (the same random $J^{ijk}_{m,l}$), and $W=0$. This charge $Q$ in the
form of Eq. 31 corresponds to
$\\{n^{0}_{1},n^{0}_{3},n^{0}_{7}\\}=\\{1,1,2\\}$ and all the other
$n^{0}_{m}=0$. Fig. 3(e) shows that there are still extensive degenerate
energy levels. The LSS shows a delta function peak at $\delta_{E}=0$ and a
Poisson distribution at $\delta_{E}>0$, in consistency with MBL. In Fig. 3(f),
we see an extensive number of low entanglement entropy eigenstates within the
entire energy range, indicating their MBL. The rest states have higher
entanglement entropy due to the energy level degeneracies, for the same reason
as described in the end of Sec. III.2.3. But the eigenstate entanglement
entropies of the $Q=21924$ sector are generically lower than those in the
$Q=3^{9}$ sector (in Fig. 3(d)), since the eigenstates are more localized with
multiple nearly decoupled subsystems.

As our Krylov subspace exact or almost solutions and ED calculations show, the
$N=3$ quantum breakdown model with $W=0$ shows a 1D MBL in generic charge $Q$
sectors when $M\rightarrow\infty$ (illustrated in Fig. 2). Note that such an
MBL can exist without any disorder and randomness: we can set $\mu_{m}=\mu$,
$M=0$, and $J^{123}_{m,l}=J$ all being constant and still have MBL.

##  IV The $N=3$ model with $W>0$

Figure 4: ED calculation in a charge $Q$ sector of the model with
$(N,M)=(3,10)$ and disorder potential $W>0$. The parameters except for $W$ are
exactly the same as Fig. 3 and labeled in the left panels. (a)-(b) has
$W=10^{-6}$, while (c)-(f) has $W=1$. (a), (c), (e) shows the energy levels
$E(\alpha)$, and their LSS in the insets. (b), (d), (f) show the eigenstate
entanglement entropy in subregion $A$ (the first $\lfloor M/2\rfloor$ sites)
for (a), (c), (e), respectively (the darker the higher density of dots).

The MBL nature of the eigenstates of the $N=3$ quantum breakdown model can be
strengthened when the disorder potential $H_{\text{dis}}$ in Eq. 4 is turned
on, namely, when the disorder strength $W>0$ (defined in Eq. 5).

In the presence of a nonzero $H_{\text{dis}}$ term, the basis transformation
in Eq. 12 can no longer simplify the model (there is no Krylov subspace
structure), and there is no obvious analytical method for solving the
eigenstates. We therefore calculate the eigenstates and spectrum of a charge
$Q$ sector numerically by ED.

Fig. 4 shows three examples of the ED energy spectrum and subregion
entanglement entropy of the model with $(N,M)=(3,10)$, where the model
parameters except for $W$ are exactly the same as those in the corresponding
panels of Fig. 3. Generically, the disorder potential $W$ further localizes
the eigenstates and strengthen the MBL, as can be seen from the lower
entanglement entropy of eigenstates in Fig. 4(b),(d),(e) as compared to Fig.
3(b),(d),(e) (where $W=0$). We discuss this result in more details below.

First, a tiny disorder potential can readily pin the originally degenerate
many-body eigenstates into a more local basis with considerably lower
entanglement entropy. Fig. 4(a)-(b) shows the numerical calculation in the
$Q=3^{9}$ sector by adding a disorder potential strength $W=10^{-6}$ to the
$(N,M)=(3,10)$ case with parameters in Fig. 3(a)-(b), after which the energy
spectrum is almost unchanged (see Figs. 3(a) and 4(a)). The inset of Fig. 4(a)
shows the LSS of in the range of $\delta_{E}\sim 10^{-8}$ (the energy scale of
level degeneracy breaking), which is still a Poisson distribution. In
contrast, the entanglement entropies of the degenerate levels in Fig. 4(b) are
much lower than those in Fig. 3(b), indicating the eigenstates are more
localized. This is because the disorder potential $W$ connects the Krylov
subspaces with different blocking sites $m_{b}$ and different occupation
configurations $\\{s^{j}_{m}\\}$, which induces a localization and degeneracy
breaking among all the originally degenerate eigenstates in different Krylov
subspaces, similar to the Anderson localization.

The energy level degeneracy breaking and MBL become stronger as the disorder
strength $W$ increases. Fig. 4(c)-(f) show the spectra in the $Q=3^{9}$ and
$Q=21924$ charge sectors at $(N,M)=(3,10)$ with $W=1$, where the other
parameters are exactly the same as those in Fig. 3(c)-(f). Compared with Fig.
3(c)-(f), the originally degenerate levels at distinct energies are broadened
and merged together by the disorder potential $W$, showing a clear Poisson
LSS. As shown in Fig. 4(d) and (f), extensive numbers of eigenstates show
extremely low entanglement entropy regardless of energy, indicating MBL. The
$Q=21924$ sector (Fig. 4(f), as a representative of generic charge sector) is
clearly much more localized than the $Q=3^{9}$ sector (Fig. 4(d)), due to the
breaking of the system into smaller decoupled or weakly coupled subsystems as
explained in Sec. III.3.2. In general, we find the MBL is robust for any value
$W>0$. The App. B Fig. 14(a)-(b) shows an example of the $N=3$ model with
$W=10$ comparable to the interaction strength $J=10$, where the MBL feature is
still clear. In the $W\gg J$ limit, the MBL is even stronger, as the disorder
potential $W$ also favors the localization of fermions.

##  V The $N>3$ models

The quantum breakdown model in Eq. 1 with $N>3$ becomes increasingly
complicated, and we do not find a method to solve it analytically. Therefore,
we employ ED to numerically compute the full energy spectrum and eigenstates
of the $N>3$ model. Within the system sizes (up to $NM\approx 25\sim 28$) and
charge $Q$ sectors calculable, we find the model exhibits a crossover from MBL
to quantum chaotic behaviors as $M/N$ decreases to $M/N\approx 1$ (see
subsection V.1 below). This suggests that the model has an MBL localization
length of order $N$ sites, and behave like a chaotic quantum dot within a
localization length. Therefore, we conjecture the model is in a 1D MBL phase
with localization length $M_{L}\sim N$ in the $M/N\rightarrow\infty$ limit
(the 1D limit), and is many-body chaotic in the $M/N\rightarrow 0$ limit (the
3D limit).

Moreover, for generic $(N,M)$ at disorder potential $W=0$, we can identify a
set of exactly solvable degenerate many-body eigenstates in many charge $Q$
sectors, which form a flat band of many-body quantum scar states. This is
discussed in subsection V.2.

###  V.1 MBL to chaos crossover as $M/N$ decreases

Similar to the $N=3$ case, we first discuss the $N>3$ quantum breakdown model
at disorder potential $W=0$, and then examine the robustness of the model
features at disorder potential $W>0$. In ED calculations, we set a uniform
chemical potential $\mu_{m}=\mu$, and random interactions with a strength $J$,
as we assumed in Eqs. 7 and 8. Throughout this section, we set $\mu=1$ and
$J=10$.

####  V.1.1 ED result

Fig. 5 shows the many-body energy spectrum and entanglement entropy of a
representative charge $Q$ sector in the $(N,M)=(4,7)$, $(5,5)$ and $(7,4)$
quantum breakdown model with $W=0$, respectively. The most obvious feature we
observe is the presence of an interval of flat dispersion in the middle of the
energy spectrum in Fig. 5(a),(c),(e), which indicates a highly degenerate
energy level in each of the charge $Q$ sectors. Such a degenerate level when
$W=0$ is observed in many charge $Q$ sectors at generic $(N,M)$ in our ED. We
call this degenerate level a _many-body scar flat band_ , the origin of which
is explained later in subsection V.2. For $N=3$, this is the level with the
highest degeneracy in Fig. 3. In addition, in the $(N,M)=(4,7)$ spectrum in
Fig. 5(a), other than the many-body scar flat band (the highest degeneracy
level in the middle), there are enormous other degenerate levels showing a
fractal-like structure (similar to Fig. 3), which cannot be understood by the
explanation later in subsection V.2. This suggests the existence of Hilbert
space fragmentation into degenerate Krylov subspaces at least also at $N=4$,
which we leave for future study.

Figure 5: ED calculation in a charge $Q$ sector of the models with $N=4,5$ and
$7$ and the largest $M$ calculable, respectively, and disorder potential
$W=0$. The parameters are labeled in the left panels. (a), (c), (e) shows the
energy levels $E(\alpha)$, and insets show their LSS. (b), (d), (f) show the
eigenstate entanglement entropy in subregion $A$ (the first $\lfloor
M/2\rfloor$ sites) for (a), (c), (e), respectively (the darker the higher
density of dots).

The insets of Fig. 5(a),(c),(e) show the LSS of the energy spectrum of the
corresponding charge $Q$ sector. First, all the LSS show a delta function peak
at $\delta_{E}=0$, which is due to the many-body scar flat band and other
degenerate levels if any. At $\delta_{E}>0$, the LSS is Poisson in Fig. 5(a)
when $(N,M)=(4,7)$, almost Wigner-Dyson when $(N,M)=(5,5)$, and clearly
Wigner-Dyson in Fig. 5(e) when $(N,M)=(7,4)$. The eigenstate entanglement
entropy of Fig. 5(a),(c),(e) in subregion A (the first $\lfloor M/2\rfloor$
sites) are given in Fig. 5(b),(d),(f), respectively. In Fig. 5(b) where
$(N,M)=(4,7)$, most eigenstates have low (area law) entanglement entropies
strongly fluctuating with energies, and the pattern resembles that of
$(N,M)=(3,10)$ in Fig. 3, showing MBL features. In contrast in Fig. 5(f) where
$(N,M)=(7,4)$, most eigenstates have high entanglement entropies, outlining a
volume law dome expected for quantum chaotic systems. In Fig. 5(d) where
$(N,M)=(5,5)$, the entanglement entropies show features in between the above
two cases. These results suggest there is a crossover from MBL behavior to
quantum chaos behavior as the ratio $M/N$ decreases to around $M/N\approx 1$.

This crossover behavior remains robust against the disorder potential $W>0$.
Fig. 6 shows calculations with disorder potential $W>0$ (equals $0.5$ or $1$,
see parameters in the panels), where all the other parameters are exactly the
same as Fig. 5. The degenerate scar level is lifted by $W>0$. When
$(N,M)=(4,7)$, most eigenstates still show low (area law) entanglement
entropies as shown in Fig. 6(b), which is similar to the $N=3$ case in Fig.
4(d) and suggests MBL. The LSS in Fig. 6(a) inset becomes semi-Poisson, which
we conjecture will tend to Poisson if $M/N$ could be larger, similar to the
$(N,M)=(3,10)$ case in Fig. 4(c) inset. This MBL-like behavior is robust for
any $W>0$. An example of the $(N,M)=(4,7)$ model with $W=10$ (comparable to
$J=10$) is shown in App. B Fig. 14(c)-(d), which show clear MBL features of
Poisson LSS and area law eigenstate entanglement entropies. The $(N,M)=(5,5)$
and $(7,4)$ cases with $W>0$ (Fig. 6(c)-(f)) are not too different from those
with $W=0$ (Fig. 5(c)-(f)), except that the eigenstate entanglement entropies
vary more smoothly with respect to the energy, showing clearer features of
volume law behaviors. Examples of more charge $Q$ sectors with $W=0$ or $W>0$
for the same $(N,M)$ can be found in the App. B Fig. 14. The LSS for smaller
$(N,M)$ are shown in App. B Fig. 16. All of these results suggest the
existence of an MBL to chaos behavior crossover as $M/N$ decreases to around
$M/N\approx 1$.

Figure 6: ED calculation in a charge $Q$ sector of the models with $N=4,5$ and
$7$ and the largest $M$ calculable, respectively, and disorder potential
$W>0$. The parameters except for $W$ are the same as Fig. 5, which are labeled
in the left panels. (a), (c), (e) shows the energy levels $E(\alpha)$, and
insets show their LSS. (b), (d), (f) show the eigenstate entanglement entropy
in subregion $A$ (the first $\lfloor M/2\rfloor$ sites) for (a), (c), (e),
respectively (the darker the higher density of dots).

Moreover, in the $(N,M)=(7,4)$ quantum chaotic case in Fig. 5(f) and Fig.
6(f), there are sharp dips of entanglement entropy in the middle of the energy
spectrum, which indicate these states are low-entanglement quantum scar
states. The dips are above and below the energy of the many-body scar flat
band in Fig. 5(e), where the density of states are low. Such quantum scar
states also occur in the $(N,M)=(5,5)$ case in Fig. 5(d) and Fig. 6(d). At
$W=0$, the many-body scar flat band states in Fig. 5(e) can also be understood
as low-entanglement quantum scar states, since adding a small disorder $W>0$
will fix the eigenbasis of the flat band into a low entanglement entropy basis
(see App. B Fig. 15). In Fig. 5(f) with $W=0$, the many-body scar flat band
states do not show low entanglement entropy, again because the ED calculation
arbitrarily picks a high entanglement entropy eigenbasis within the degenerate
Hilbert space. Remarkably, the fraction of such flat band quantum scar
eigenstates is not measure zero compared to the total number of eigenstates
(see subsection V.2). This is in stark contrast to previous quantum scar
models, where the number of scar states is usually measure zero in the
thermodynamic limit.

####  V.1.2 Argument for the crossover from MBL to chaos

To define the MBL and quantum chaos phases, certain thermodynamic limit need
be taken. As we explained in Sec. II.4, we define $M/N\rightarrow\infty$ as
the 1D thermodynamic limit (where $N$ can be a small number), and
$M/N\rightarrow 0$ as the 3D thermodynamic limit (where both $M$ and $N$ are
large). Our small size calculations then suggest the model is MBL in the 1D
thermodynamic limit, and quantum chaotic in the 3D thermodynamic limit, and
the crossover possibly happens at $M/N\approx 1$. We emphasize that this
conclusion may be subject to finite size effects. Nonetheless, we provide here
a theoretical argument for the crossover around $M/N\approx 1$.

The key is to identify the localization length $M_{L}$ (number of sites) of
the model. Recall that a fermion on the $m$-th site can move to the $(m+1)$-th
site and turn into three fermions. If two of these three fermions stay in the
$(m+1)$-th site, and the other fermion further move and turn into three
fermions on the $(m+2)$-th site (as is the example in Eq. 14), it appears this
process can go on to any distance away without localization. However, we argue
this is not true. This is because in this picture, only one fermion in the
last occupied site is mobile (subject to randomness induced localization),
while the mobilities of all the other fermions rely on this single fermion far
way, which would not be a bulk property. To define the localization length
$M_{L}$ as a bulk property, we need to allow at least one fermion per site to
be mobile independently within a localization length $M_{L}$. Since a fermion
needs to produce at least two more fermions per site to move forward, it
requires a Hilbert space of $2\cdot(M_{L}/2)=M_{L}$ fermion modes if it is to
move by $M_{L}/2$ sites, the mid-value of a localization interval $[0,M_{L}]$.
Within an interval of $M_{L}$ sites, to allow one mobile fermion per site, we
would need a Hilbert space of $M_{L}\cdot M_{L}=M_{L}^{2}$ number of fermion
modes; while the total number of fermion modes available in the interval is
$NM_{L}$. Equalizing these two numbers yields a MBL localization length

$M_{L}\sim N\ .$ (33)

Therefore, we expect the model to show MBL features when $M>M_{L}$, or
$M/N>1$. In contrast, when $M<M_{L}$, or $M/N<1$, the model should behave as a
random extended system, which is usually quantum chaotic.

###  V.2 Many-body scar flat band at $W=0$

We have seen in Fig. 5 (see also App. B Fig. 14) that at $W=0$, many charge
$Q$ sectors of the quantum breakdown model have an extensively degenerate
level in the middle of the spectrum, which we name as the _many-body scar flat
band_. Moreover, we find all the eigenstates $|\alpha\rangle$ in this many-
body scar flat band has integer number of fermions on each site:

$\langle\alpha|\hat{n}_{m}|\alpha\rangle=n_{m}\in\mathbb{Z}\ ,\quad(n_{m}\geq
0)$ (34)

and their eigen-energy is given by

$E(\alpha)=\sum_{m=1}^{M}\mu_{m}n_{m}\ .$ (35)

This indicates that the interaction $H_{I}$ in Eq. 2 has no contribution to
the eigen-energy of these states at all, so they are the zero modes of $H_{I}$
satisfying $H_{I}|\alpha\rangle=0$. We note that Eqs. 34 and 35 are
generically not true for energy levels not in the many-body scar flat band
(for instance, not true for any degenerate levels other than the many-body
scar flat band at energy $E=13$ in the $(N,M)=(4,7)$ spectrum in Fig. 5(a)).

Figure 7: (a)-(c) The classical degeneracy $d_{2}^{0}(n_{1},n_{2})$ (defined
in Eq. 37) for several charge $Q$ sectors $(N,M)=(15,2)$ as a function of
$n_{1}$ (recall that $n_{2}=Q-3n_{1}$). Each $Q$ sector has only one
$(n_{1},n_{2})$ configuration with $d_{2}^{0}(n_{1},n_{2})>0$. (d)-(f) The
maximal $d_{2}^{0}$ of each charge $Q$ sector, calculated for
$(N,M)=(15,2),(28,2)$ and $(29,2)$, respectively. The $y$ axis is
$\log_{10}(d_{2}^{0}+1)$ (such that zero corresponds to $d_{2}^{0}=0$).

The origin of such a many-body scar flat band here resembles that of the
single-particle flat bands of line-graph tight-binding models [40, 41] (e.g.,
the Lieb lattice [42] and the kagome lattice [43]), which is due to the
mismatch of dimensions of the Hilbert subspaces connected by the Hamiltonian
$H_{I}$, as we will show below.

####  V.2.1 $M=2$ sites

The simplest example is the quantum breakdown model with $M=2$ sites. Consider
the Hilbert space of Fock states with $n_{m}$ ($n_{m}\in\mathbb{Z}$) fermions
on site $m$ ($m=1,2$), which has a Hilbert space dimension

$h_{n_{1},n_{2}}=C_{N}^{n_{1}}C_{N}^{n_{2}}\ ,$ (36)

where $C_{N}^{n}=\frac{N!}{(N-n)!n!}$ is the binomial coefficient. Note that
$C_{N}^{n}=0$ if $n<0$ or $n>N$. The conserved charge $Q$ is given by
$Q=3n_{1}+n_{2}$. The interaction $H_{I}$ in Eq. 2 can hop from the Hilbert
space of fermion numbers $\\{n_{1},n_{2}\\}$ into the Hilbert space of fermion
numbers $\\{n_{1}\pm 1,n_{2}\mp 3\\}$. The $H_{I}$ hopping matrix between
these two Hilbert spaces thus has $h_{n_{1}+1,n_{2}-3}+h_{n_{1}-1,n_{2}+3}$
rows and $h_{n_{1},n_{2}}$ columns. If the Hilbert space dimensions for some
integers $\\{n_{1},n_{2}\\}$ satisfy

$d_{2}^{0}(n_{1},n_{2})=h_{n_{1},n_{2}}-h_{n_{1}+1,n_{2}-3}-h_{n_{1}-1,n_{2}+3}>0\
,$ (37)

the above $H_{I}$ hopping matrix will have a rank no larger than
$h_{n_{1}+1,n_{2}-3}+h_{n_{1}-1,n_{2}+3}$, and thus there will necessarily be

$d_{2}(n_{1},n_{2})\geq d_{2}^{0}(n_{1},n_{2})$ (38)

number of linearly independent states in the $\\{n_{1},n_{2}\\}$ Hilbert space
that can be annihilated by $H_{I}$. Namely, the interaction Hamiltonian
$H_{I}$ will have $d_{2}$ zero modes with fermion numbers $\\{n_{1},n_{2}\\}$,
which will be eigenstates with energy given by Eq. 35 when the disorder
potential $W=0$. These zero modes of $H_{I}$ therefore form a many-body scar
flat band. We call $d_{2}^{0}(n_{1},n_{2})$ in Eq. 37 the _classical
degeneracy_ (from simple dimension counting), and the actual degeneracy
$d_{2}(n_{1},n_{2})$ the _quantum degeneracy_.

Fig. 7(a)-(c) shows the classical degeneracy $d_{2}^{0}(n_{1},n_{2})$ defined
in Eq. 37 as a function of $n_{1}$ in several charge $Q$ sectors for
$(N,M)=(15,2)$ (recall $n_{2}=Q-3n_{1}$). Each of them has one configuration
$\\{n_{1},n_{2}\\}$ with $d_{2}^{0}(n_{1},n_{2})>0$, and thus has a many-body
scar flat band. Generically, for $N\leq 28$, most charge $Q$ sectors have one
(and usually only one) fermion number configuration $\\{n_{1},n_{2}\\}$ with
$d_{2}^{0}(n_{1},n_{2})>0$, as shown in Fig. 7(d)-(e). When $N\geq 29$, a
charge $Q$ sector can have $d_{2}^{0}(n_{1},n_{2})>0$ only if $Q\lesssim
N^{2/3}$ or $4N-Q\lesssim N^{2/3}$, and the example of $N=29$ is shown in Fig.
7(f). We give a theoretical understanding and estimation of such behaviors
with respect to $N$ in App. C.

Within the ED calculable system sizes with $M=2$, we find $d_{2}=d_{2}^{0}$ in
all charge sectors except for one case: the $Q=2N$ sector when $N$ is odd. In
this case, within $N\leq 7$, we find $d_{2}-d_{2}^{0}\leq 2$. Therefore, the
classical degeneracy $d_{2}^{0}$ gives a very accurate estimation of the
quantum degeneracy $d_{2}$.

####  V.2.2 Generic $M$ sites

The counting of many-body scar flat band degeneracy in Eq. 37 for $M=2$ sites
can be straightforwardly generalized to $M>2$ sites. Consider a fermion number
configuration $\\{n_{1},n_{2},\cdots,n_{M}\\}$ on the $M$ sites in a charge
sector $Q$. We assume the first $m$ sites ($1\leq m\leq M$) with the sub-
configuration $\\{n_{1},n_{2},\cdots,n_{m}\\}$ has $d_{M}^{(m)}$ zero modes of
$H_{I}$. When we add the $m$-th site to the first $m-1$ sites, the
$d_{M}^{(m)}$ zero modes of the first $m$ sites should span a subspace of the
direct product of the zero-mode Hilbert space of the first $m-1$ sites
(dimension $d_{M}^{(m-1)}$) and the $n_{m}$-fermion Hilbert space of the
$m$-th site (dimension $C_{N}^{n_{m}}$). This is because the terms in $H_{I}$
within the first $m-1$ sites have to annihilate these $d_{M}^{(m)}$ zero
modes. Then, these $d_{M}^{(m)}$ zero modes also need be annihilated by the
interaction $J_{m-1,l}^{ijk}$ between sites $m-1$ and $m$. From our discussion
in Sec. V.2.1, the matrix of the $J_{m-1,l}^{ijk}$ term within the two sites
$m-1$ and $m$ has a rank no larger than
$h_{n_{m-1}+1,n_{m}-3}+h_{n_{m-1}-1,n_{m}+3}$, where $h_{n_{1},n_{2}}$ is the
two-site Hilbert space dimension of fermion number $\\{n_{1},n_{2}\\}$ defined
in Eq. 36. Therefore, the number of modes that the $J_{m-1,l}^{ijk}$ term
cannot annihilate is at most its rank times the number of zero modes
$d_{M}^{(m-2)}$ of the first $m-2$ sites. This indicates the number of zero
modes $d_{M}^{(m)}$ of the first $m$ sites satisfy the recursion relation

$\begin{split}&d_{M}^{(m)}=d_{M}^{(m-1)}C_{N}^{n_{m}}\\\
&-d_{M}^{(m-2)}(h_{n_{m-1}+1,n_{m}-3}+h_{n_{m-1}-1,n_{m}+3})+r_{M}^{(m)},\end{split}$
(39)

where

$r_{M}^{(m)}\geq 0$ (40)

is a nonnegative quantum correction beyond the above rank counting, which may
arise from destructive quantum interferences. The two initial values of the
recursion relation are given by

$d_{M}^{(0)}=1\ ,\qquad d_{M}^{(1)}=C_{N}^{n_{1}}\ ,$ (41)

where we allowed $m$ to start from $0$ for convenience. The number of zero
modes, namely, quantum degeneracy $d_{M}(\\{n_{m}\\})$ of the many-body scar
flat band in the total $M$ sites is given by

$d_{M}=d_{M}^{(M)}$ (42)

from the recursion relation in Eq. 39.

Figure 8: Many-body scar flat band for $(N,M)=(4,3)$ with $W=0$ and random
interactions. (a) The quantum degeneracy $d_{M}$ from ED and the classical
degeneracy $d_{M}^{0}$ computed from Eq. 44 in all charge $Q$ sectors. (b) The
difference $d_{M}-d_{M}^{0}$ for all $Q$. (c)-(d) Fermion numbers $n_{1}$ and
$n_{2}$ in the many-body scar flat band in each charge $Q$ sector. Missing
points indicate the absence of a many-body scar flat band.

We have not found a general formula for the quantum correction $r_{M}^{(m)}$
in Eq. 40. From ED calculations, we find $r_{M}^{(m)}$ is zero in most charge
$Q$ sectors. In the rare cases where $r_{M}^{(m)}>0$, it is always much
smaller than $d_{M}^{(m)}$, so it only gives a small correction to
$d_{M}^{(m)}$. We therefore define a classical degeneracy
$d_{M}^{0}(\\{n_{m}\\})$ as

$d_{M}^{0}=d_{M}^{(M),0}\ ,$ (43)

where $d_{M}^{(M),0}$ is defined by the recursion relation ignoring the
quantum correction $r_{M}^{(m)}$, namely,

$\begin{split}&d_{M}^{(m),0}=d_{M}^{(m-1),0}C_{N}^{n_{m}}\\\ &\quad-
d_{M}^{(m-2),0}(h_{n_{m-1}+1,n_{m}-3}+h_{n_{m-1}-1,n_{m}+3})\ ,\end{split}$
(44)

with initial values $d_{M}^{(0),0}=1$ and $d_{M}^{(1),0}=C_{N}^{n_{1}}$. We
expect the classical degeneracy $d_{M}^{0}$ not to differ much from the actual
quantum degeneracy $d_{M}$.

Similar to the $M=2$ case in Sec. V.2.1, for small $N$ (at least for $N$
calculable in our ED) and $M>2$, most charge sectors $Q$ have one (or very
rarely two, see an example in App. D Fig. 17(e)) fermion number configuration
$\\{n_{m}\\}$ with $d_{M}>0$, leading to a many-body scar flat band. As an
example, Fig. 8(a) shows the many-body scar flat band quantum degeneracy
$d_{M}$ from ED and classical degeneracy $d_{M}^{0}$ from Eq. 44 in all charge
sectors $Q$ for $N=4,M=3$, and Fig. 8(b) shows their difference
$d_{M}-d_{M}^{0}$. We see that in most charge $Q$ sectors, $d_{M}=d_{M}^{0}$
precisely. There are only four charge sectors with $d_{M}-d_{M}^{0}>0$, but
this difference is much smaller than $d_{M}$. Fig. 8(c)-(d) shows the fermion
numbers on the first site and second site of the many-body scar flat band,
respectively (missing points meaning no many-body scar flat band in that
charge $Q$ sector).

##  VI The breakdown transition

To examine if the quantum breakdown model can describe the dielectric
breakdown phenomena, we numerically calculate the time-evolution of a fermion
added into the particle vacuum state of the system. We consider uniform
chemical potential $\mu_{m}=\mu$ and random interactions with strength $J>0$.
By the chiral symmetry (Sec. II.3), we can assume $\mu\geq 0$ without loss of
generality.

As explained in Sec. II.4, in the dielectric breakdown picture, $\mu$ can be
understood as the ionization energy for producing a fermion, while $J$ can be
viewed as the electric field induced potential energy difference between
neighboring layers (the energy a fermion can gain upon moving forward by one
site). When $J>2\mu$, the energy a fermion gains when moving a site forward
will be able to overcome the energy cost of producing two more fermions,
allowing the particle avalanche to happen. Therefore, we expect the breakdown
transition to happen at

$\frac{\mu}{J}=\frac{1}{2}\ ,\qquad(\mu\geq 0)\ .$ (45)

As we will show, this point indeed corresponds to a breakdown transition in
the quantum breakdown model.

###  VI.1 Energy spectrum dependence on $\mu/J$

We first investigate the dependence of many-body energy spectrum on the ratio
$\mu/J$. In Figs. 3 to 6, we have shown the calculations only at $\mu/J\leq
0.1$. Fig. 9 shows the energy spectrum with $W=0$ in the $Q=3^{M-1}$ sector at
larger $\mu/J$. We find that at small $\mu/J$, the energy spectrum does not
have large gaps (Fig. 9(a),(c)), while at large $\mu/J$, the energy spectrum
splits into many-body bands with gaps of order $2\mu$ (Fig. 9(b),(d)). This is
because as $\mu/J\rightarrow\infty$, or $J\rightarrow 0$, the energy spectrum
is simply given by $H_{\mu}$ in Eq. 3, which forms bands with an interband gap
$2\mu$. However, we observed no sharp changes but only an adiabatic crossover
in the energy spectrum as $\mu/J$ crosses $1/2$. Therefore, if a sharp
breakdown transition exists at $\mu/J=1/2$, the change in eigen-wavefunctions
must play an essential role.

Figure 9: The dependence of the energy spectrum on the chemical potential
$\mu_{m}=\mu$, for fixed random interactions with strength $J=10$ and disorder
potential $W=0$, and the charge sector is $Q=3^{M-1}$. (a)-(b) has
$(N,M)=(4,7)$, while (c)-(d) has $(N,M)=(7,4)$. At large $\mu$, the spectrum
splits into bands with gaps around $2\mu$. The LSS has no qualitative
dependence on $\mu$.

We also find that the LSS of the energy spectrum is insensitive to $\mu/J$. As
shown in the insets of Fig. 9, the LSS of the $(N,M)=(4,7)$ model at
$\mu/J=\frac{1}{2}$ and $2$ are both Poisson (Fig. 9(a)-(b)), while the LSS of
the $(N,M)=(7,4)$ model at $\mu/J=\frac{1}{2}$ and $2$ are both Wigner-Dyson
(Fig. 9(c)-(d)). This suggests that the MBL to chaos crossover at $M/N\approx
1$ is independent of $\mu/J$.

###  VI.2 Time evolution

To reveal the breakdown transition, we examine the dynamical evolution of an
initial state $|\psi(0)\rangle$ which adds one fermion to the first site of
the particle vacuum $|0\rangle$. Since the interactions are random, we can
assume the fermion is in the first fermion mode of site $m=1$ without loss of
generality, namely:

$|\psi(0)\rangle=c_{1,1}^{\dagger}|0\rangle\ ,$ (46)

where $c_{m,i}^{\dagger}$ is the fermion creation operator in Eq. 2. Such an
initial state resembles creating one ionized atom in an undisturbed dielectric
gas. We then employ ED to calculate the time-evolved state
$|\psi(t)\rangle=e^{-iHt}|\psi(0)\rangle$ at time $t$. Note that
$|\psi(t)\rangle$ is always in the charge $Q=3^{M-1}$ sector. Here we always
take random interactions with interaction strength $J=10$, and uniform
chemical potential $\mu_{m}=\mu$. Note that the wavefunctions
$|\psi(t)\rangle$ at chemical potentials $\mu$ and $-\mu$ are related by the
chiral symmetry $C$ (Eq. 10) followed by complex conjugation, giving the same
dynamics. We thus restrict ourselves to $\mu\geq 0$.

Figure 10: The long-time ($t\in[50,100]$) mean value (solid lines) and range
of fluctuation (yellow shaded regions) of the number of fermions $n_{1}$ on
site $m=1$ for the initial state in Eq. 46, as a function of chemical
potential $\mu_{m}=\mu$. The parameters are given in the panels. The vertical
dashed line shows the position where $\mu/J=1/2$.

We define the number of fermions on the $m$-th site at time $t$ as

$n_{m}(t)=\langle\psi(t)|\hat{n}_{m}|\psi(t)\rangle\ ,$ (47)

with the number operator $\hat{n}_{m}$ defined in Eq. 3. The breakdown can be
observed by measuring the number of fermions $n_{1}(t)$ on the first site. If
the breakdown happens, we expect $n_{1}(t)$ to approach zero as time
$t\rightarrow\infty$, since the fermion will turn into many fermions on the
$m>1$ sites. In contrast, if the breakdown does not happen, one expects a
nonzero $n_{1}(t)$ as $t\rightarrow\infty$, since the fermion cannot induce a
fermion avalanche and is trapped near the first site.

We calculate the time evolution within $t\in[0,100]$, which is sufficient for
examining the long time behavior of the evolution. Fig. 10 (a)-(h) show the
late time behavior of the model with different $(N,M)$ (with $W=0$ in the left
panels and $W=1$ in the right panels), where we vary $\mu\in[0,25]$ while
fixing $J=10$. The black solid line in each panel shows the long-time mean
value $\langle n_{1}\rangle$ of the fermion number $n_{1}(t)$ in the late time
window $t\in[50,100]$. The yellow shaded region shows the range of fluctuation
of $n_{1}(t)$ in the time window $t\in[50,100]$, which has boundaries at
$\langle n_{1}\rangle\pm\sqrt{\langle\delta n_{1}^{2}\rangle}$, where
$\sqrt{\langle\delta n_{1}^{2}\rangle}$ is the root mean square of the
fluctuation $\delta n_{1}(t)=n_{1}(t)-\langle n_{1}\rangle$. The model
parameters are given in each panel.

Figure 11: The number of fermions $n_{m}(t)/3^{m-1}$ on site $m$ as a function
of time $t$, for $(N,M)=(3,10)$ with the initial state in Eq. 46. $n_{1}(t)$
is highlighted as the thick solid lines. The parameters are given at the top
of each panel (the same $H_{I}$ and $H_{\text{dis}}$ as Fig. 10(a),(b)).

The most prominent feature is, regardless of the values of $N$ and $M$, the
long-time mean value $\langle n_{1}\rangle$ shows a transition around $\mu=5$,
or $\mu/J=1/2$, which is exactly the transition point we expected in Eq. 45.
On the $\mu/J>1/2$ side, the mean value $\langle n_{1}\rangle$ curve is
(mostly) smoothly approaching $1$ as $\mu$ grows, with the fluctuation
$\sqrt{\langle\delta n_{1}^{2}\rangle}$ becoming smaller and smaller. This is
the expected dielectric behavior before breakdown, since the fermion is frozen
on the first site when $\mu\rightarrow\infty$. The behavior of $\mu/J<1/2$
side depends on $N/M$ (for total fermion modes $NM$ around $25\sim 30$). In
Fig. 10 (a)-(d) where $N/M<1$, the value $\langle n_{1}\rangle$ shows
oscillations with respect to $\mu$. A special case is $(N,M)=(3,10)$ with
disorder potential $W=1$ in Fig. 10(b), where $\langle n_{1}\rangle$ shows a
sharp transition from intense oscillation to no oscillation as $\mu/J$ crosses
$1/2$. In contrast, in Fig. 10(e)-(h) with $N/M\geq 1$, $\langle n_{1}\rangle$
is nearly zero, consistent with our expectation for the breakdown phase. This
is the most obvious in the $(N,M)=(7,4)$ case with disorder potential $W=1$ in
Fig. 10(h), where $\langle n_{1}\rangle$ show a clear transition from almost
zero to nonzero at $\mu/J=1/2$. Hereafter, we call $\mu/J>1/2$ the _dielectric
regime_ , and $\mu/J<1/2$ the _breakdown regime_.

Figure 12: The number of fermions $n_{m}(t)/3^{m-1}$ on site $m$ as a function
of time $t$, for $(N,M)=(7,4)$ with the initial state in Eq. 46. $n_{1}(t)$ is
highlighted as the thick solid lines. The parameters are given on top of each
panel (the same $H_{I}$ and $H_{\text{dis}}$ as Fig. 10(g)).

We now investigate the time evolutions more closely by examining $n_{m}(t)$ as
a function of time $t$. Since our calculation is in the charge $Q=3^{M-1}$
sector, we have

$\sum_{m=1}^{M}\frac{n_{m}(t)}{3^{m-1}}=\frac{Q}{3^{M-1}}=1\ .$ (48)

The $N=3$ case is relatively special compared to the $N>3$ cases, due to its
almost analytically solvable nature. Fig. 11 shows $n_{m}(t)/3^{m-1}$ in the
$(N,M)=(3,10)$ case, where $n_{1}(t)$ is highlighted as the thick solid lines.
When $W=0$, the difference between $\mu/J<1/2$ and $\mu/J>1/2$ is not
significant, both having $\langle n_{1}\rangle$ far from zero, as Fig. 10(a)
implies. In both examples of $\mu/J=0.4$ and $\mu/J=0.6$ in Fig. 11(a)-(b),
$n_{m}(t)/3^{m-1}$ show persistent quasi-periodic oscillations with $t$. This
can be understood by noting that the initial state in Eq. 46 is in either the
Krylov subspace without blocking in Eq. 16 or the one-dimensional Krylov
subspaces in Eq. 26 in the $Q=3^{M-1}$ sector of the $N=3$ model with $W=0$,
which is equivalent to a tight-binding model with a single effective particle
(Eq. 17) and thus show no sharp changes as $\mu/J$ varies. When a finite
disorder potential $W=1$ is added, the $\mu/J>1/2$ side remains qualitatively
unchanged. The $\mu/J<1/2$ side, however, is significantly altered, which
develops a strong irregular dependence on chemical potential $\mu$ (see Fig.
10(b)). Fig. 11(c) and (d) show two examples with $W=1$ at $\mu/J=0.4$ and
$\mu/J=0.435$, respectively. In Fig. 11(c), the oscillation in $n_{1}(t)$ is
suppressed by the disorder potential $W$ (compared to Fig. 11(a)). In Fig.
11(d) (note the range of $t$ plotted is different from Fig. 11(a)-(c)), not
only the short period oscillation in $n_{1}(t)$ is not suppressed, but also a
persistent long period oscillation arises. The underlying mechanism of such
behaviors remains unclear, but these oscillations clearly indicate the absence
of thermal equilibrium in the $N=3$ model.

The $N>3$ cases do not show significant differences between disorder strength
$W=0$ and $W>0$, except that $\langle n_{1}\rangle$ in the breakdown regime
$\mu/J<1/2$ is generically smaller when $W>0$, as Fig. 10 shows. We therefore
only discuss the $W=0$ case here. In Fig. 12, we plot the time evolution of
$n_{m}(t)/3^{m-1}$ for $(N,M)=(7,4)$ with $W=0$, where $\mu/J$ is $0.4$, $0.5$
and $0.6$ in Fig. 12(a)-(c), respectively. It is clear that $n_{1}(t)$ decays
to zero in the breakdown regime $\mu/J<1/2$, and tends to certain oscillations
around a nonzero value in the dielectric regime $\mu/J>1/2$. We note that even
if the $(N,M)=(7,4)$ model is quantum chaotic as indicated by its Wigner-Dyson
LSS (see Fig. 9), the dielectric regime $\mu/J>1/2$ shows persistent
oscillations in fermion densities and does not thermalize.

Based on the above observations, we conjecture that, when both $N$ and $M$ are
large enough, irrespective of the ratio $M/N$, the breakdown transition
happens robustly at $\mu/J=1/2$, which is signaled by whether $n_{1}(t)$ tends
to zero at large $t$. This is because the breakdown picture only requires
sufficiently large number of fermion modes in the subsequent sites. However,
we expect the breakdown phase to have a difference between $M/N>1$ and
$M/N<1$. When $M/N>1$, our earlier analysis in Sec. V.1 implies that the
system shows 1D MBL with a localization length around $N$, so we expect the
fermion to spread only by a length around $N$ sites in the breakdown regime.
Since we expect this MBL to behave as a chaotic quantum dot within a
localization length $N$, we expect a local thermal equilibrium within the
length $N$ is reached. We call such a phase a _local breakdown_. When $M/N<1$,
the model is globally quantum chaotic within a charge $Q$ sector, so we expect
the spreading of the fermion to approach a global thermal equilibrium in the
breakdown regime, and we call this phase a _global breakdown_.

Figure 13: The conjectured phase diagram of the quantum breakdown model at
large $N$ and $M$, as a function of $\mu/J$ and $M/N$.

Our conjectures above are summarized into a phase diagram of the quantum
breakdown model at large $N$ and $M$, as shown in Fig. 13. A transition
between the dielectric phase and the breakdown phase happens at $\mu/J=1/2$
irrespective of $M/N$. When $M/N\gg 1$ (the 1D limit), the system is in a MBL
phase, and only local breakdown within order $N$ sites near the initial
perturbation can happen when $\mu/J<1/2$. When $M/N\ll 1$ (the 3D limit), the
system is quantum chaotic (with quantum scar states in certain charge $Q$
sectors), and a global breakdown in the entire space is expected when
$\mu/J<1/2$.

##  VII Discussion

We have shown evidences that our quantum breakdown model exhibit a crossover
from MBL to quantum chaos as $M/N$ decreases across $1$, and undergoes a
transition from dielectric to breakdown as $\mu/J$ decreases across $1/2$,
which are summarized in the phase diagram Fig. 13. At $N=3$ and with disorder
potential $W=0$, the model exhibit a Hilbert space fragmentation into
exponentially many Krylov subspaces, and is exactly solvable except for very
few Krylov subspaces in certain charge $Q$ sectors. The solution indicates
features of MBL. Remarkably, the MBL features can occur in the absence of any
translational symmetry breaking randomness, namely, when $\mu_{m}=\mu$,
$J^{ijk}_{m,l}=J$ are all uniformly constant and $W=0$. The MBL is robust
against a nonzero disorder potential $W>0$. We note that our model show rather
distinct features (the local and global breakdown, etc) compared to the
interacting Wannier-Stark type models [17, 18, 19, 20, 21, 22], although both
are describing particles in an electric field. This is because our underlying
settings are different: in our setup, the electrons do not feel the electric
field until they are ionized from the atoms. For generic $(N,M)$ with $W=0$,
we also show that a many-body scar flat band of quantum scar states exists in
many charge $Q$ sectors, the origin of which can be understood as the many-
body generalization of the line-graph tight-binding models such as the Lieb
lattice and kagome lattice models [40, 41, 42, 43]. It will be intriguing to
generalize this idea to construct more interacting models with many-body flat
bands or quantum scar states.

The $N=3$ quantum breakdown model, as an almost exactly solvable MBL model at
$W=0$ as we demonstrated, may allow one to further explore analytically the
MBL quantum dynamics, conserved quantities and quantum entanglement [44, 45,
46, 47, 48], in contrast to many MBL models such as the disordered magnetic
field XXZ model [12, 13, 14, 15, 16] which heavily relies on many-body
numerical calculations. Perturbation analysis on the $N=3$ model may give a
deeper understanding to the quantum avalanche transition from MBL to thermal
states studied recently [49, 50, 51, 52, 53, 54]. For the $N>3$ models, it may
also be possible to find exact or almost exact solutions for certain non-
random interactions, which may provide more insights to the phase diagram we
conjectured in Fig. 13. For $N>3$ with random interactions, Hilbert space
fragmentation into degenerate Krylov subspaces may still exist, as indicated
by the enormous fractal-like degenerate levels other than the many-body scar
flat band in the $N=4$ case in Fig. 5(a).

In the large $N$ limit, with random interactions satisfying Eq. 8, it is
possible to study our quantum breakdown model using the large $N$ expansion
techniques of the SYK models [8, 9, 10, 11]. It would be particularly
interesting to investigate the MBL to chaos crossover and breakdown transition
we conjectured in Fig. 13 in this SYK limit. This is, however, beyond the
scope of this paper, and will be studied in a separate future paper.
Nevertheless, we make a few general observations here. This SYK limit of our
model is obviously different from the previous SYK dot lattice models with
spatially symmetric interactions (upon random average) [55, 56, 57, 58, 59].
In our model, one expect the quantum chaos and scrambling to develop towards a
fixed direction, in analogy to the chiral SYK models [60, 61] and generic
chaotic chiral models [62, 63], although our model here is not chiral as a
lattice model. The conserved charge $Q$ in our model is also very different
from the on-site U(1) charge in complex SYK models [64]. Moreover, one expects
a richer phase diagram in the SYK limit of our model. A $1/N$ expansion may be
required for investigating the MBL to chaos crossover with respect to $M/N$.
We do note that MBL was studied before in SYK models with spatially symmetric
interactions [58].

Our quantum breakdown model further allows a variety of directions of
generalizations. First, the number of fermion modes $N$ on each site can be
generalized to a site-dependent number $N_{m}$, which allows one to more
freely change the geometry of the “dielectric gas container”. In this case,
the formula for the many-body scar flat band degeneracy in Eq. 39 simply
changes through the substitutions $C_{N}^{n_{m}}\rightarrow C_{N_{m}}^{n_{m}}$
and $h_{n_{m-1},n_{m}}\rightarrow C_{N_{m-1}}^{n_{m-1}}C_{N_{m}}^{n_{m}}$.
Secondly, one could assume certain spatial structure within each layer of
dielectric gas in Fig. 1(a), such that the interaction $J_{m,l}^{ijk}$ only
acts among nearby fermions, generalizing the model into higher spatial
dimensions. Moreover, the model can be generalized to a wide class of models
of fermions, bosons or spins with different types of asymmetric interactions,
which may give rise to a rich variety of non-equilibrium phenomena. It will
also be useful to investigate the effect of adding bilinear inter-site fermion
hopping terms to our model, allowing a fermion to move without creating
additional fermions, which may destabilize the MBL or the breakdown phases.
However, such hopping terms break the conserved charge $Q$, which complicates
the model and requires theoretical or numerical methods beyond this work. We
leave the studies of these questions to the future.

###### Acknowledgements.

The author thank David Huse, Yingfei Gu, Zhenbin Yang and Elliott Lieb for
insightful discussions. This work is supported by the Alfred P. Sloan
Foundation, the National Science Foundation through Princeton University’s
Materials Research Science and Engineering Center DMR-2011750, and the
National Science Foundation under award DMR-2141966. Additional support is
provided by the Gordon and Betty Moore Foundation through Grant GBMF8685
towards the Princeton theory program.

##  Appendix A The interaction between subsystems in arbitrary charge $Q$
sectors of the $N=3$ model

This appendix is to supplement our discussion in Sec. III.3.2, to write down
more explicitly the interaction between neighboring subsystems in certain
Krylov subspaces of an arbitrary charge $Q$ sector of the $N=3$ model with
$W=0$ given in Eq. 13.

Following the notations defined in and below Eq. 31, each charge $Q$ uniquely
corresponds to a reference configuration $\\{n_{m}^{0}\\}$ with $0\leq
n_{m}^{0}\leq 2$, and $m_{f,p}$ is the $p$-th site with $n_{m}^{0}>0$ ($1\leq
p\leq M_{Q}$). Generalizing Eq. 21, we define the following Fock state basis
$|\\{m_{s,p}\\},\\{m_{b,p}\\},\\{s_{m}^{j}\\}\rangle$ in a generic charge $Q$
sector:

$\begin{split}&|\\{m_{s,p}\\},\\{m_{b,p}\\},\\{s_{m}^{j}\\}\rangle\\\
&=\prod_{p=1}^{M_{Q}}\Big{(}f_{m_{s,p},s^{1}_{m_{s,p}}}^{\dagger}\prod_{m=m_{b,p}+1}^{m_{s,p}}f^{\dagger}_{m,2}f^{\dagger}_{m,3}\\\
&\times\prod_{m=m_{f,p}+1}^{m_{b,p}}f^{\dagger}_{m,s^{1}_{m}}f^{\dagger}_{m,s^{2}_{m}}\prod_{j=2}^{n_{m_{f,p}}}f^{\dagger}_{m_{f,p},s^{j}_{m_{f,p}}}\Big{)}|0\rangle,\end{split}$
(49)

where we restrict the blocking sites $m_{b,p}$ and spreading sites $m_{s,p}$
to satisfy

$m_{f,p}\leq m_{b,p}\leq m_{s,p}\leq m_{f,p+1}\ ,$ (50)

and we always fix

$s_{m_{b,p}}^{1}=1\ ,$ (51)

which is necessary for defining $m_{b}^{p}$ as a blocking site. In addition,
for the state we defined in Eq. 49 to be nonzero, we require
$s_{m_{s,p}}^{1}=1$ whenever $m_{s,p}>m_{f,p}$, and
$s^{2}_{m^{\text{nz}}_{p}}=2$ or $3$. All the other indices $s_{m}^{j}$ can
generically take values $1,2$ or $3$ (distinct for the same $m$).

In this representation, each Krylov subspace is labeled by
$\\{m_{b,p}\\},\\{s_{m}^{j}\\}$ for a given charge $Q$ sector (i.e., given
$\\{m_{f,p}\\}$). The action of the Hamiltonian in each Krylov subspace can be
generically written down as

$\begin{split}&H|\\{m_{s,p}\\},\\{m_{b,p}\\},\\{s_{m}^{j}\\}\rangle\\\
&=(H_{\text{hop}}+H_{\text{b-int}})|\\{m_{s,p}\\},\\{m_{b,p}\\},\\{s_{m}^{j}\\}\rangle\
,\end{split}$ (52)

where $H_{\text{hop}}$ is a part that can be viewed as hoppings of a non-
interacting tight-binding model:

$\begin{split}&H_{\text{hop}}|\\{m_{s,p}\\},\\{m_{b,p}\\},\\{s_{m}^{j}\\}\rangle=\sum_{p^{\prime}=1}^{M_{Q}}\Big{[}\delta_{s^{1}_{m_{s,p^{\prime}}},1}\Big{(}\vartheta_{m_{s,p^{\prime}}-m_{b,p^{\prime}}-1+\delta_{m_{b,p^{\prime}},m_{f,p^{\prime}}}}J_{m_{s,p^{\prime}}-1}|\\{m_{s,p}-\delta_{p,p^{\prime}}\\},\\{m_{b,p}\\},\\{s_{m}^{j}\\}\rangle\\\
&\qquad\qquad\qquad+\vartheta_{m_{f,p^{\prime}+1}-m_{s,p^{\prime}}-1}J_{m_{s,p^{\prime}}}|\\{m_{s,p}+\delta_{pp^{\prime}}\\},\\{m_{b,p}\\},\\{s_{m}^{j}\\}\rangle\Big{)}+V_{m_{s,p^{\prime}}}|\\{m_{s,p}\\},\\{m_{b,p}\\},\\{s_{m}^{j}\\}\rangle\Big{]}\
,\end{split}$ (53)

where $\vartheta_{m}$ is the integer variable Heaviside function defined in
Eq. 15, $J_{m}$ is defined in Eq. 6, and $V_{m}=\sum_{m=1}^{M}\mu_{m}n_{m}$
with $n_{m}$ being the site $m$ number of fermions of state
$|\\{m_{s,p}\\},\\{m_{b,p}\\},\\{s_{m}^{j}\\}\rangle$. Eq. 53 is nothing but
the generalization of Eq. 23, namely, an effective tight-binding Hamiltonian
for some effective particles on sites $m_{s,p}$ within subregions $m_{f,p}\leq
m_{s,p}<m_{f,p+1}$. The term $H_{\text{b-int}}$ is a boundary interaction term

$\begin{split}&H_{\text{b-int}}|\\{m_{s,p}\\},\\{m_{b,p}\\},\\{s_{m}^{j}\\}\rangle=\sum_{p^{\prime}=2}^{M_{Q}}\Big{(}\delta_{m_{s,p^{\prime}-1},m_{f,p^{\prime}}}\delta_{n_{m_{f,p^{\prime}}},3}J_{m_{f,p^{\prime}}-1}|\\{m_{s,p}-\delta_{p,p^{\prime}-1}\\},\\{m_{b,p}\\},\\{s_{m}^{j}\\}\rangle\\\
&+\delta_{m_{s,p^{\prime}-1},m_{f,p^{\prime}}-1}\delta_{n_{m_{f,p^{\prime}}},0}J_{m_{f,p^{\prime}}-1}|\\{m_{s,p}+\delta_{p,p^{\prime}-1}\\},\\{m_{b,p}\\},\\{s_{m}^{j}\\}\rangle\Big{)}\
,\end{split}$ (54)

where $n_{m_{f,p^{\prime}}}$ in the delta functions
$\delta_{n_{m_{f,p^{\prime}}},3}$ and $\delta_{n_{m_{f,p^{\prime}}},0}$ is the
number of fermions on site $m_{f,p^{\prime}}$ of the state
$|\\{m_{s,p}\\},\\{m_{b,p}\\},\\{s_{m}^{j}\\}\rangle$. This dependence on the
number of fermions $n_{m_{f,p^{\prime}}}$ makes Eq. 54 an interaction term for
the effective particles located at sites $m_{s,p}$, although it is confined to
the boundary sites $m_{f,p^{\prime}}$. Note that if
$m_{b,p^{\prime}}>m_{f,p^{\prime}}$, the fermion number $n_{m_{f,p^{\prime}}}$
cannot be changed by fermions on the right of site $m_{f,p^{\prime}}$, and one
will always have
$\delta_{n_{m_{f,p^{\prime}}},3}=\delta_{m_{s,p^{\prime}-1},m_{f,p^{\prime}}}$,
and
$\delta_{n_{m_{f,p^{\prime}}},0}=\delta_{m_{s,p^{\prime}-1},m_{f,p^{\prime}}-1}$,
in which case Eq. 54 becomes a free term for the effective particles on sites
$m_{s,p}$ as an additional part of Eq. 53. This interaction between the
$(p^{\prime}-1)$-th and $p^{\prime}$-th subsystem is only present when
$n_{f,p^{\prime}}^{0}=1$ and the Krylov subspace has
$m_{b,p^{\prime}}=m_{f,p^{\prime}}$, making it a very weak interaction. This
interaction is a hardcore repulsion on site $m_{f,p^{\prime}}$ between the two
effective particles on sites $m_{s,p^{\prime}-1}$ and $m_{s,p^{\prime}}$, as
is obvious from the fact that $m_{s,p^{\prime}-1}=m_{s,p^{\prime}}$ is
forbidden.

Figure 14: ED calculation of more charge $Q$ sectors of the models with
$N=3,4,5$ and $7$ and the largest $M$ calculable. The parameters are labeled
in the panels, where random interaction $J=10$, the chemical potential is
uniformly $\mu_{m}=\mu=1$. The first and third columns are the energy levels
$E(\alpha)$ and their LSS in the insets. The second and fourth columns show
the eigenstate entanglement entropy in subregion $A$ (the first $\lfloor
M/2\rfloor$ sites) for the energy spectra on their left (the darker the higher
density of dots). The disorder potential is $W=10$ in (a)-(d), to show that
the MBL behaviors in the $M/N>1$ case remain robust when $W$ is comparable to
$J$. We set $W=0$ in panels (e), (f), (i), (j), (m), (n), and $W=1$ or $0.5$
in panels (g), (h), (k), (l), (o), (p).

##  Appendix B Additional ED results for $N>3$ models

This appendix shows more examples of ED calculations in different charge $Q$
sectors of the $N>3$ models.

Fig. 14 show additional charge $Q$ sectors of the models with $N=3,4,5$ and
$7$ and the largest $M$ calculable, as a supplement to Figs. 3, 4, 5 and 6 in
the main text (the same system sizes). The detailed descriptions and
parameters are given in the caption. The results show the same trend as the
main text figures show, that the system tends to MBL if $M/N>1$ (with Poisson
or semi-Poisson LSS and many eigenstates with area law entanglement entropy,
subject to finite system sizes), and behaves quantum chaotically when $M/N<1$
(with Wigner-Dyson LSS and most eigenstates with volume law entanglement
entropy).

Figure 15: Changes of the eigenstate entanglement entropy of the charge
$Q=3^{3}$ (panels (a)-(b)) and $Q=30$ (panels (c)-(d)) sectors of the
$(N,M)=(7,4)$ model upon adding a small disorder potential $W=0.01$ (the right
two panels). The parameters are given in the panels. The many-body scar flat
band is at energy $E=7$ in the charge $Q=3^{3}$ sector, and is at energy $E=8$
in the $Q=30$ sector.  Figure 16: The energy spectrum $E(\alpha)$ (from low to
high) and LSS (the insets) of various small system sizes $(N,M)$ (see values
in the panels, with $NM\approx 16\sim 27$) with $W=0$, where the charge sector
is always $Q=3^{M-1}$, the interaction is random with interaction strength
$J=10$, and the chemical potential is uniformly $\mu_{m}=\mu=1$. The LSS shows
a clear crossover from Poisson to Wigner-Dyson as $M/N$ decreases across $1$.

The many-body scar flat band in a charge $Q$ sector can be immediately fixed
into a localized eigenbasis when a small disorder potential $W$ is added,
which breaks the degeneracy. Fig. 15 shows two examples of the eigenstate
entanglement entropies in the charge $Q=3^{3}$ and $Q=30$ sectors in the
$(N,M)=(7,4)$ model. When a disorder potential $W=0.01$ is added to the
system, while the entanglement entropies of most eigenstates are almost
unchanged, the entanglement entropy of the many-body scar flat band in the
middle immediately drops significantly, indicating the real space localization
of the eigenbasis.

To test the MBL to chaos behavior crossover around $M/N\approx 1$, we further
show in Fig. 16 the energy spectra and LSS of smaller system sizes $(N,M)$ at
$W=0$, where we always choose the $Q=3^{M-1}$ sector as an example. The
interaction is random with a strength $J$, and chemical potential is uniformly
$\mu_{m}=\mu=1$. One can clearly see a crossover from Poisson to Wigner-Dyson
LSS as $M/N$ decreases across $1$, in consistency with a MBL to chaos behavior
crossover.

##  Appendix C More on the classical degeneracy of many-body scar flat band
for $M=2$ sites

Here for the model with $M=2$ sites, we give a theoretical quantitative
understanding of the dependence of the classical degeneracy $d_{2}^{0}$ of
many-body scar flat band on $N$ as we showed in Fig. 7(d)-(f). In particular,
a transition is observed between $N=28$ and $N=29$ in Fig. 7: for $N\leq 28$,
most charge $Q$ sectors have a many-body scar flat band, while for $N\geq 29$,
only charge sectors with $Q$ close to $0$ and $4N$ have a many-body scar flat
band. We explain these two regimes in two subsections below, respectively.

By Eq. 37, classical degeneracy $d_{2}^{0}(n_{1},n_{2})$ is given by

$d_{2}^{0}(n_{1},n_{2})=h_{n_{1},n_{2}}-h_{n_{1}+1,n_{2}-3}-h_{n_{1}-1,n_{2}+3}\
,$ (55)

where $h_{n_{1},n_{2}}=C_{N}^{n_{1}}C_{N}^{n_{2}}$. We first note that for
$M=2$, the number of fermions $n_{2}$ on site $m=2$ can be represented by
charge $Q$ and the number of fermions $n_{1}$ on site $m=1$ as

$n_{2}=Q-3n_{1}\ .$ (56)

Therefore, we can express $h_{n_{1},n_{2}}$ as a function of $Q$ and $n_{1}$
as

$\begin{split}&h_{n_{1},n_{2}}=h_{n_{1},Q-3n_{1}}=C_{N}^{n_{1}}C_{N}^{Q-3n_{1}}\\\
&=\frac{N!^{2}}{n_{1}!(N-n_{1})!(Q-3n_{1})!(N-Q+3n_{1})!}\ .\end{split}$ (57)

We now examine the condition for a positive classical degeneracy
$d_{2}^{0}>0$.

###  C.1 Small $N$

As an approximation, as long as $N$ is not too small, we can regard $\ln
h_{n_{1},Q-3n_{1}}$ as a continuous function of $n_{1}$. The reason to take a
log is because $h_{n_{1},Q-3n_{1}}$ can change by several orders of magnitude
as a function of $n_{1}$. For a fixed charge $Q$, we can approximate the
classical degeneracy $d_{2}^{0}(n_{1},n_{2})=d_{2}^{0}(n_{1},Q-3n_{1})$ in Eq.
55 by Taylor expansion as

$\begin{split}&d_{2}^{0}(n_{1},Q-3n_{1})\approx h_{n_{1},Q-3n_{1}}\\\
&\times\left[-\frac{d^{2}\ln h_{n_{1},Q-3n_{1}}}{dn_{1}^{2}}-\left(\frac{d\ln
h_{n_{1},Q-3n_{1}}}{dn_{1}}\right)^{2}-1\right].\end{split}$ (58)

By a straightforward calculation, we find

$\begin{split}&\frac{d\ln h_{n_{1},Q-3n_{1}}}{dn_{1}}\approx\ln(N-n_{1})-\ln
n_{1}\\\ &\qquad\qquad+3\ln(Q-3n_{1})-3\ln(N-Q+3n_{1})\ ,\end{split}$ (59)

and

$\begin{split}&\frac{d^{2}\ln
h_{n_{1},Q-3n_{1}}}{dn_{1}^{2}}\approx-\frac{1}{N-n_{1}}-\frac{1}{n_{1}}\\\
&\qquad\qquad\qquad\qquad-\frac{9}{Q-3n_{1}}-\frac{9}{N-Q+3n_{1}}\
.\end{split}$ (60)

Figure 17: (a)-(d) The comparison between the classical degeneracy $d_{M}^{0}$
(blue dots) from Eq. 44 and quantum degeneracy $d_{M}$ (red circles) in ED,
for $(N,M)=(4,4)$ and $(N,M)=(5,3)$, respectively. The parameters are
$\mu_{m}=\mu=1$, random interaction strength $J=10$, and $W=0$. (e)-(f) The
energy spectrum of two example charge sectors in the $(N,M)=(5,3)$ model,
which has two and one many-body scar flat bands (red), respectively.

In order to obtain a positive classical degeneracy $d_{2}^{0}$ in Eq. 58, we
should have $\frac{d^{2}\ln h_{n_{1},Q-3n_{1}}}{dn_{1}^{2}}<0$, and at the
same time $\frac{d\ln h_{n_{1},Q-3n_{1}}}{dn_{1}}$ as close to zero as
possible. Numerically, it can be shown that $\frac{d\ln
h_{n_{1},Q-3n_{1}}}{dn_{1}}=0$ when

$n_{1}\approx\frac{Q}{4}\ ,$ (61)

at which one has

$\frac{d^{2}\ln
h_{n_{1},Q-3n_{1}}}{dn_{1}^{2}}\approx-40\left(\frac{1}{Q}+\frac{1}{4N-Q}\right)\
.$ (62)

Note that Eq. 61 implies $n_{2}=Q-3n_{1}\approx Q/4$. Therefore, we estimate
the largest classical degeneracy $d_{2}^{0}$ is reached around the fermion
number configuration $n_{1}\approx n_{2}\approx Q/4$, estimated by the formula

$\begin{split}&d_{2}^{0}(Q)\approx
h_{n_{1},Q-3n_{1}}\left[40\left(\frac{1}{Q}+\frac{1}{4N-Q}\right)-1\right]\
.\end{split}$ (63)

It is easy to see that $d_{2}^{0}$ is minimal when $Q=2N$ (note that $0\leq
Q\leq 4N$). Therefore, if $d_{2}^{0}>0$ in the $Q=2N$ sector, we would expect
almost all (not all because this estimation is very crude) charge $Q$ sectors
to have a many-body scar flat band. Since

$\begin{split}&d_{2}^{0}(Q=2N)\approx
h_{n_{1},Q-3n_{1}}\left(\frac{40}{N}-1\right)\ ,\end{split}$ (64)

which is positive only if $N<40$, we conclude that the criterion for having
many-body flat band in almost all charge $Q$ sectors is $N\lesssim 40$. In the
actual numerical calculation in Fig. 7(d)-(f), this criterion is $N\leq 28$,
which is reasonably not far from our crude estimation.

###  C.2 Large $N$

As we have shown above, for large $N$ ($N\geq 29$ from numerical calculation
in Fig. 7(d)-(f)), the many-body scar flat band will not occur in most charge
$Q$ sectors. However, here we show that there is always a many-body scar flat
band when $Q$ is close enough to $0$ or $4N$. Since the classical degeneracy
in Eq. 55 is symmetric between $Q$ and $4N-Q$, we only consider small $Q$
here.

The proof is easy by noting that the configuration $(n_{1},n_{2})=(0,Q)$
always has a positive classical degeneracy $d_{2}^{0}$ at sufficiently small
$Q$. By Eq. 55, one has

$\begin{split}&d_{2}^{0}(0,Q)=h_{0,Q}-h_{1,Q-3}\\\
&=C_{N}^{Q}\left[1-\frac{NQ(Q-1)(Q-2)}{(N-Q+1)(N-Q+2)(N-Q+3)}\right]\\\
&\approx C_{N}^{Q}\left(1-\frac{Q^{3}}{N^{2}}\right)\ .\end{split}$ (65)

Therefore, at large $N$, when $Q<N^{2/3}$, there is always a many-body scar
flat band with fermion number configuration $(n_{1},n_{2})=(0,Q)$. Similarly,
when $Q>4N-N^{2/3}$, there is always a many-body scar flat band with fermion
number configuration $(n_{1},n_{2})=(N,Q-3N)$.

##  Appendix D Additional results for many-body scar flat band with $M>2$
sites

In Fig. 17, we provide additional ED results as a supplementary to the main
text Fig. 8. Fig. 17(a)-(d) show the comparison between the quantum degeneracy
$d_{M}$ and classical degeneracy $d_{M}^{0}$ (calculated by Eq. 44) of many-
body scar flat band in all the charge $Q$ sectors, where the system sizes are
$(N,M)=(4,4)$ and $(N,M)=(5,3)$, respectively. Again, the difference between
$d_{M}$ and $d_{M}^{0}$ is zero in most charge $Q$ sectors, and small compared
to $d_{M}$ in all the charge $Q$ sectors.

In addition, we show a very rare charge $Q$ sector example in Fig. 17(e) which
has two many-body scar flat bands, in contrast to most charge $Q$ sectors
which has only one many-body scar flat band (e.g., Fig. 17(f)).

## References

  * Townsend [1911] J. S. Townsend, The Theory of Ionisation of Gases by Collision, Nature 85, 400 (1911).
  * Rutherford and Geiger [1908] E. Rutherford and H. Geiger, An electrical method of counting the number of $\alpha$-particles from radio-active substances, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 81, 141 (1908).
  * Geiger and Müller [1928] H. Geiger and W. Müller, Elektronenzählrohr zur Messung schwächster Aktivitäten, Naturwissenschaften 16, 617 (1928).
  * Jensen and Shankar [1985] R. V. Jensen and R. Shankar, Statistical behavior in deterministic quantum systems with few degrees of freedom, Phys. Rev. Lett. 54, 1879 (1985).
  * 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).
  * 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 .
  * Sachdev and Ye [1993] S. Sachdev and J. Ye, Gapless spin fluid ground state in a random, quantum Heisenberg magnet, Phys. Rev. Lett. 70, 3339 (1993), arXiv:cond-mat/9212030 [cond-mat] .
  * Polchinski and Rosenhaus [2016] J. Polchinski and V. Rosenhaus, The Spectrum in the Sachdev-Ye-Kitaev Model, JHEP 04, 001, arXiv:1601.06768 [hep-th] .
  * Maldacena and Stanford [2016] J. Maldacena and D. Stanford, Remarks on the Sachdev-Ye-Kitaev model, Phys. Rev. D94, 106002 (2016), arXiv:1604.07818 [hep-th] .
  * Kitaev and Suh [2018] A. Kitaev and S. J. Suh, The soft mode in the Sachdev-Ye-Kitaev model and its gravity dual, JHEP 05, 183, arXiv:1711.08467 [hep-th] .
  * 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–1205 (2006).
  * Gornyi  _et al._ [2005] I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Interacting electrons in disordered wires: Anderson localization and low-$t$ transport, Phys. Rev. Lett. 95, 206603 (2005).
  * Oganesyan and Huse [2007] V. Oganesyan and D. A. Huse, Localization of interacting fermions at high temperature, Phys. Rev. B 75, 155111 (2007).
  * Žnidarič  _et al._ [2008] M. Žnidarič, T. c. v. Prosen, and P. Prelovšek, Many-body localization in the heisenberg $xxz$ magnet in a random field, Phys. Rev. B 77, 064426 (2008).
  * Pal and Huse [2010] A. Pal and D. A. Huse, Many-body localization phase transition, Phys. Rev. B 82, 174411 (2010).
  * van Nieuwenburg  _et al._ [2019] E. van Nieuwenburg, Y. Baum, and G. Refael, From bloch oscillations to many-body localization in clean interacting systems, Proceedings of the National Academy of Sciences 116, 9269 (2019).
  * Schulz  _et al._ [2019] M. Schulz, C. A. Hooley, R. Moessner, and F. Pollmann, Stark many-body localization, Phys. Rev. Lett. 122, 040606 (2019).
  * Moudgalya  _et al._ [2021] S. Moudgalya, A. Prem, R. Nandkishore, N. Regnault, and B. A. Bernevig, Thermalization and its absence within krylov subspaces of a constrained hamiltonian, in  _Memorial Volume for Shoucheng Zhang_ (WORLD SCIENTIFIC, 2021) pp. 147–209.
  * Herviou  _et al._ [2021] L. Herviou, J. H. Bardarson, and N. Regnault, Many-body localization in a fragmented hilbert space, Phys. Rev. B 103, 134207 (2021).
  * Sala  _et al._ [2020] P. Sala, T. Rakovszky, R. Verresen, M. Knap, and F. Pollmann, Ergodicity breaking arising from hilbert space fragmentation in dipole-conserving hamiltonians, Phys. Rev. X 10, 011047 (2020).
  * Khemani  _et al._ [2020] V. Khemani, M. Hermele, and R. Nandkishore, Localization from hilbert space shattering: From theory to physical realizations, Phys. Rev. B 101, 174204 (2020).
  * Žnidarič [2013] M. Žnidarič, Coexistence of diffusive and ballistic transport in a simple spin ladder, Phys. Rev. Lett. 110, 070602 (2013).
  * Yang  _et al._ [2020] Z.-C. Yang, F. Liu, A. V. Gorshkov, and T. Iadecola, Hilbert-space fragmentation from strict confinement, Phys. Rev. Lett. 124, 207602 (2020).
  * Moudgalya  _et al._ [2020] S. Moudgalya, B. A. Bernevig, and N. Regnault, Quantum many-body scars in a landau level on a thin torus, Phys. Rev. B 102, 195150 (2020).
  * Bernien  _et al._ [2017] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Probing many-body dynamics on a 51-atom quantum simulator, Nature 551, 579 (2017).
  * Moudgalya  _et al._ [2018] S. Moudgalya, N. Regnault, and B. A. Bernevig, Entanglement of exact excited states of affleck-kennedy-lieb-tasaki models: Exact results, many-body scars, and violation of the strong eigenstate thermalization hypothesis, Phys. Rev. B 98, 235156 (2018).
  * Schecter and Iadecola [2018] M. Schecter and T. Iadecola, Many-body spectral reflection symmetry and protected infinite-temperature degeneracy, Phys. Rev. B 98, 035139 (2018).
  * Turner  _et al._ [2018a] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Weak ergodicity breaking from quantum many-body scars, Nature Physics 14, 745 (2018a).
  * Turner  _et al._ [2018b] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Quantum scarred eigenstates in a rydberg atom chain: Entanglement, breakdown of thermalization, and stability to perturbations, Phys. Rev. B 98, 155134 (2018b).
  * Choi  _et al._ [2019] S. Choi, C. J. Turner, H. Pichler, W. W. Ho, A. A. Michailidis, Z. Papić, M. Serbyn, M. D. Lukin, and D. A. Abanin, Emergent su(2) dynamics and perfect quantum many-body scars, Phys. Rev. Lett. 122, 220603 (2019).
  * Ho  _et al._ [2019] W. W. Ho, S. Choi, H. Pichler, and M. D. Lukin, Periodic orbits, entanglement, and quantum many-body scars in constrained models: Matrix product state approach, Phys. Rev. Lett. 122, 040603 (2019).
  * Bull  _et al._ [2019] K. Bull, I. Martin, and Z. Papić, Systematic construction of scarred many-body dynamics in 1d lattice models, Phys. Rev. Lett. 123, 030601 (2019).
  * Anderson [1958] P. W. Anderson, Absence of diffusion in certain random lattices, Phys. Rev. 109, 1492 (1958).
  * Wannier [1962] G. H. Wannier, Dynamics of band electrons in electric and magnetic fields, Rev. Mod. Phys. 34, 645 (1962).
  * Berry  _et al._ [1977] M. V. Berry, M. Tabor, and J. M. Ziman, Level clustering in the regular spectrum, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 356, 375 (1977), https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1977.0140 .
  * Bohigas  _et al._ [1984] O. Bohigas, M. J. Giannoni, and C. Schmit, Characterization of chaotic quantum spectra and universality of level fluctuation laws, Phys. Rev. Lett. 52, 1 (1984).
  * Wigner [1967] E. P. Wigner, Random matrices in physics, SIAM Review 9, 1 (1967).
  * Dyson [1970] F. J. Dyson, Correlations between eigenvalues of a random matrix, Communications in Mathematical Physics 19, 235 (1970).
  * Mielke [1991] A. Mielke, Ferromagnetic ground states for the hubbard model on line graphs, Journal of Physics A: Mathematical and General 24, L73 (1991).
  * Călugăru  _et al._ [2021] D. Călugăru, A. Chew, L. Elcoro, Y. Xu, N. Regnault, Z.-D. Song, and B. A. Bernevig, General construction and topological classification of crystalline flat bands, Nature Physics 18, 185 (2021).
  * Lieb [1989] E. H. Lieb, Two theorems on the hubbard model, Phys. Rev. Lett. 62, 1201 (1989).
  * Bergman  _et al._ [2008] D. L. Bergman, C. Wu, and L. Balents, Band touching from real-space topology in frustrated hopping models, Phys. Rev. B 78, 125104 (2008).
  * Serbyn  _et al._ [2013] 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 (2013).
  * 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. Müller, and A. Scardicchio, Integrals of motion in the many-body localized phase, Nuclear Physics B 891, 420–465 (2015).
  * Lian [2022] B. Lian, Conserved quantities from entanglement hamiltonian, Phys. Rev. B 105, 035106 (2022).
  * De Roeck and Huveneers [2017] W. De Roeck and F. m. c. Huveneers, Stability and instability towards delocalization in many-body localization systems, Phys. Rev. B 95, 155129 (2017).
  * Thiery  _et al._ [2018] T. Thiery, F. m. c. Huveneers, M. Müller, and W. De Roeck, Many-body delocalization as a quantum avalanche, Phys. Rev. Lett. 121, 140601 (2018).
  * Luitz  _et al._ [2017] D. J. Luitz, F. m. c. Huveneers, and W. De Roeck, How a small quantum bath can thermalize long localized chains, Phys. Rev. Lett. 119, 150602 (2017).
  * Goihl  _et al._ [2019] M. Goihl, J. Eisert, and C. Krumnow, Exploration of the stability of many-body localized systems in the presence of a small bath, Phys. Rev. B 99, 195145 (2019).
  * Crowley and Chandran [2020] P. J. D. Crowley and A. Chandran, Avalanche induced coexisting localized and thermal regions in disordered chains, Phys. Rev. Research 2, 033262 (2020).
  * 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).
  * Berkooz  _et al._ [2017] M. Berkooz, P. Narayan, M. Rozali, and J. Simón, Higher dimensional generalizations of the SYK model, Journal of High Energy Physics 2017, 10.1007/jhep01(2017)138 (2017).
  * Gu  _et al._ [2017] Y. Gu, X.-L. Qi, and D. Stanford, Local criticality, diffusion and chaos in generalized sachdev-ye-kitaev models, Journal of High Energy Physics 2017, 10.1007/jhep05(2017)125 (2017).
  * Davison  _et al._ [2017] R. A. Davison, W. Fu, A. Georges, Y. Gu, K. Jensen, and S. Sachdev, Thermoelectric transport in disordered metals without quasiparticles: The sachdev-ye-kitaev models and holography, Phys. Rev. B 95, 155131 (2017).
  * Jian and Yao [2017] S.-K. Jian and H. Yao, Solvable sachdev-ye-kitaev models in higher dimensions: From diffusion to many-body localization, Phys. Rev. Lett. 119, 206602 (2017).
  * Chen  _et al._ [2017] Y. Chen, H. Zhai, and P. Zhang, Tunable quantum chaos in the sachdev-ye-kitaev model coupled to a thermal bath, Journal of High Energy Physics 2017, 10.1007/jhep07(2017)150 (2017).
  * Lian  _et al._ [2019] B. Lian, S. L. Sondhi, and Z. Yang, The chiral SYK model, Journal of High Energy Physics 2019, 10.1007/jhep09(2019)067 (2019).
  * Hu and Lian [2021] Y. Hu and B. Lian, The Chiral Sachdev-Ye Model: Integrability and Chaos of Anyons in 1+1d, arXiv e-prints , arXiv:2109.13263 (2021), arXiv:2109.13263 [cond-mat.str-el] .
  * Hu and Lian [2022] Y. Hu and B. Lian, Integrability of the v=4/3 fractional quantum hall edge states, Quantum 6, 736 (2022).
  * Lian [2022] B. Lian, Interacting 1D Chiral Fermions with Pairing: Transition from Integrable to Chaotic, arXiv e-prints , arXiv:2202.06169 (2022), arXiv:2202.06169 [cond-mat.str-el] .
  * Gu  _et al._ [2020] Y. Gu, A. Kitaev, S. Sachdev, and G. Tarnopolsky, Notes on the complex sachdev-ye-kitaev model, Journal of High Energy Physics 2020, 10.1007/jhep02(2020)157 (2020).