# Floquet Engineering of Hilbert Space Fragmentation in Stark Lattices

Li Zhang  Institute of Quantum Precision Measurement, State Key Laboratory of
Radio Frequency Heterogeneous Integration, Shenzhen University, Shenzhen
518060, China  College of Physics and Optoelectronic Engineering, Shenzhen
University, Shenzhen 518060, China  Yongguan Ke  Laboratory of Quantum
Engineering and Quantum Metrology, School of Physics and Astronomy, Sun Yat-
Sen University (Zhuhai Campus), Zhuhai 519082, China  Institute of Quantum
Precision Measurement, State Key Laboratory of Radio Frequency Heterogeneous
Integration, Shenzhen University, Shenzhen 518060, China  College of Physics
and Optoelectronic Engineering, Shenzhen University, Shenzhen 518060, China
Ling Lin  Institute of Quantum Precision Measurement, State Key Laboratory of
Radio Frequency Heterogeneous Integration, Shenzhen University, Shenzhen
518060, China  College of Physics and Optoelectronic Engineering, Shenzhen
University, Shenzhen 518060, China  Chaohong Lee  Email<EMAIL_ADDRESS>Institute of Quantum Precision Measurement, State Key Laboratory of Radio
Frequency Heterogeneous Integration, Shenzhen University, Shenzhen 518060,
China  College of Physics and Optoelectronic Engineering, Shenzhen University,
Shenzhen 518060, China  Quantum Science Center of Guangdong-Hongkong-Macao
Greater Bay Area (Guangdong), Shenzhen 518045, China

###### Abstract

The concept of Hilbert space fragmentation (HSF) has recently been put forward
as a routine to break quantum ergodicity. While HSF widely exists in dynamical
constraint models, it is still challenging to tune HSF. Here, we propose a
scheme to tune HSF in a one-dimensional tilted lattice of interacting spinless
fermions with periodically driven tunneling. The dynamics is governed by
effective Hamiltonians with kinetic constraints, which appear as density-
dependent tunneling in the weak-tunneling perturbation expansion. The kinetic
constraint can be tuned via changing the driving frequency, and three
different kinds of strong HSF can be engineered. In general, the system is
strongly constrained and exhibits a strong HSF. Two partial resonance
frequencies are _analytically_ given by a time-dependent perturbation theory
for Floquet systems, at which some kinetic constraints are released and the
system exhibits another two different strong HSF. We demonstrate the
perturbation analysis with exact numerical simulation of the entanglement
entropy, the density correlation functions and the saturated local density
profiles. Our result provides a promising way to control HSF through Floquet
engineering.

##  I Introduction

Recently, there has been much interests in exploring whether and how an
isolated quantum many-body system can reach thermal equilibrium under unitary
dynamics. While ergodic systems can reach thermal equilibrium rapidly via
eigenstate thermalization hypothesis Deutsch (1991); Srednicki (1994, 1999);
Rigol and Dunjko (2008); Abanin  _et al._ (2016), the failure of
thermalization due to ergodicity breaking has been predicted in several
systems Nandkishore and Huse (2015); Abanin  _et al._ (2019); Moudgalya  _et
al._ (2022). Two well-known examples for ergodicity breakdown are quantum
integrable systems and many-body localized systems, both of which possess an
extensive number of conserved quantities Kinoshita  _et al._ (2006); Rigol
_et al._ (2007); Polkovnikov  _et al._ (2011); Serbyn  _et al._ (2013); Huse
_et al._ (2014); Nandkishore and Huse (2015); Abanin  _et al._ (2019).
Moreover, exotic violation of ergodicity has been found in systems subjected
to strong tilting potentials Schulz  _et al._ (2019); van Nieuwenburg  _et
al._ (2019); Taylor  _et al._ (2020); Yao and Zakrzewski (2020); Chanda  _et
al._ (2020); Zhang  _et al._ (2021); Yao  _et al._ (2021); Wang  _et al._
(2021); Guo  _et al._ (2021); Morong  _et al._ (2021), which support quantum
many-body scars Turner  _et al._ (2018a, b); Ho  _et al._ (2019); Serbyn  _et
al._ (2021); Moudgalya  _et al._ (2022); Chandran  _et al._ (2023) and
fractured Hilbert space Pai  _et al._ (2019); Sala  _et al._ (2020); Khemani
_et al._ (2020); Moudgalya  _et al._ (2021, 2022).

The Hilbert space fragmentation (HSF) describes the phenomenon that the
Hilbert space of a system is exponentially split into many dynamically
disjoint invariant subspaces (referred to as Krylov subspaces), which can not
be captured by the conventional symmetry Pai  _et al._ (2019); Sala  _et al._
(2020); Khemani  _et al._ (2020); Moudgalya  _et al._ (2021, 2022). Thus,
large parts of the Hilbert space are inaccessible to certain initial states
and the ergodicity breaks down regarding the full system. The HSF arises due
to dynamical constraints, such as fractonlike constraints of conservation of
U(1) charge and its associated dipole Pai  _et al._ (2019); Sala  _et al._
(2020); Khemani  _et al._ (2020); Moudgalya  _et al._ (2021) and local
tunneling constraints from strong interaction or strong tilting potentials De
Tomasi  _et al._ (2019); Yang  _et al._ (2020); Frey  _et al._ (2022); Ghosh
_et al._ (2023); Scherg  _et al._ (2021); Kohlert  _et al._ (2023). Depending
on whether the largest Krylov subspace can span almost the entire symmetry
space or not in the thermodynamic limit, the HSF is categorized as weak and
strong respectively Sala  _et al._ (2020); Moudgalya  _et al._ (2022).
Although the HSF can be different when changing the filling number Morningstar
_et al._ (2020); Pozderac  _et al._ (2023), to our best knowledge, there is
still no way to engineer the HSF by tuning the physical parameters.

Floquet engineering, the coherent control of a quantum system via periodic
driving, is a powerful tool for engineering synthetic Hamiltonians with novel
properties Eckardt (2017). It has been used to realize topological Bloch
bands, dynamical localization and synthetic gauge field, etc Eckardt (2017).
Specifically, density-dependent tunneling has been engineered by suitable
driving schemes Eckardt (2017); Greschner  _et al._ (2014); Meinert  _et al._
(2016); Zhao  _et al._ (2019). Recently, it have been proposed to generate
quantum many-body scars and HSF via tuning local tunneling constraints to some
specific parameters Zhao  _et al._ (2020); Ghosh  _et al._ (2023). However,
the way to adjust the HSF via Floquet engineering is still lacking.

In this article, we provide a way to tune the HSF in a periodically driven
Stark chain of interacting spinless fermions. The fermions with nearest
neighboring interactions are subjected to a tilted field, and they tunnel
between nearest neighboring sites with periodically varying tunneling
strength. When the tunneling strength is much smaller than the tilting
strength, the interaction strength and the driving frequency, the system at
stroboscopic times can be viewed as transitions between multiple energy levels
by absorbing or emitting “photons” from the driving field. The energy levels
are formed by the eigenergies of the tilting and interacting Hamiltonian. The
transitions between different energy levels happen through the action of the
hopping Hamiltonian, with frequencies and amplitudes given by the Fourier
expansion of the tunneling strength. Thus, there appears a resonant transition
when the Fourier frequencies match the energy difference before and after the
hopping, which depends on the particle numbers near the hopping sites. The
non-resonant transitions can be neglected in the presence of resonant ones.
Inspired by this picture, we can tune the relation between the driving
frequency and the interaction and tilting strength to design different
Hamiltonians with different kinetic constraints, in which some particular
hoppings are non-resonant thus forbidden. As a concrete example, we set the
interaction and tilting strength being equal, at which an intrinsic resonance
is always present, and tune the driving frequency to obtain other resonant
transitions. In the absence of the driving, the hopping is constrained to
conserve the total energy of interaction and tilting and leads to strong HSF.
In the presence of the periodic driving, we engineer two extra resonant
hoppings at two particular resonant frequencies respectively, which are
analytically given by a perturbation method for Floquet system. The released
kinetic constraint changes the structure of the Hilbert space and leads to two
different kinds of HSF. At non-resonant frequencies, the driving does not
assist the hopping, and the HSF is the same as that in the absence of the
periodic driving. We demonstrate the perturbation analysis through exact
numerical calculation of the dynamics of entanglement entropy, the density
autocorrelation function and the saturated density profiles. Our result
provides a route to controlled switching between different HSF through tuning
the driving frequency.

The rest of the paper is organized as follows. In Sec. II, we introduce the
model for our physical system and present a time-dependent perturbation theory
for Floquet system to derive the effective Hamiltonians. In Sec. III, we
present how to engineer the HSF by tuning the driving frequency. At last, in
Sec. IV, we give a brief summary of our results.

##  II Model and Method

###  II.1 Model

We consider an ensemble of interacting spinless fermions in a one-dimensional
Stark lattice under periodic driving, see Fig. 1. It is described by the
Hamiltonian $\hat{H}(t)=\hat{H}_{0}+\hat{H}_{\mathrm{d}}(t)$ with the time-
independent part

$\hat{H}_{0}=\hat{H}_{J}+U\sum_{j=0}^{L-2}\hat{n}_{j}\hat{n}_{j+1}-g\sum_{j=0}^{L-1}j\hat{n}_{j},$
(1)

where the hopping term
$\hat{H}_{J}=J\sum_{j=0}^{L-2}(\hat{c}^{\dagger}_{j}\hat{c}_{j+1}+\hat{c}^{\dagger}_{j+1}\hat{c}_{j})$,
and the time periodic part

$\hat{H}_{\mathrm{d}}(t)=u(t)\hat{H}_{J}.$ (2)

Here, $\hat{c}^{\dagger}_{j}(\hat{c}_{j})$ creates (annihilates) a fermion at
site $j$, $\hat{n}_{j}=\hat{c}^{\dagger}_{j}\hat{c}_{j}$ is the particle
number operator, and $u(t)=u\cos(\omega t)$ is periodic in time with frequency
$\omega$ and amplitude $u$. The parameters $J$, $U$ and $g$ are the nearest-
neighbor tunneling strength, the nearest-neighbor interaction strength and the
tilting field strength, respectively. $L$ is the total number of lattice
sites. In Appendix A, we provide a discussion of the realization of our model
in Rydberg atom platform. In the following study, we consider open boundary
condition and focus on the half-filling sector with particle number $N=L/2$.
The energy unit is set as $J=1$.

Figure 1:  Scheme diagram of interacting spinless fermions in a periodically
driven tilted lattice. The tilting and nearest neighboring interaction
strength are $g$ and $U$ respectively. The fermions tunnel between nearest
neighboring sites with periodic varying strength $J[1+u(t)]$, where
$u(t)=u\cos\omega t$.

###  II.2 Time dependent perturbation theory in Floquet systems

Generally, the dynamics of a periodically driven system $\hat{H}(t)$ after a
period is governed by its Floquet Hamiltonian

$\hat{H}_{F}=\frac{i}{T}\ln\hat{F}.$ (3)

Here the Floquet operator $\hat{F}$ is the unitary evolution operator over a
period:

$\hat{F}=\mathcal{T}\mathrm{exp}\left[-i\int_{0}^{T}\hat{H}(t)dt\right],$ (4)

with $\mathcal{T}$ denoting time ordering. In general cases, $\hat{H}_{F}$ is
untractable, as $\hat{H}(t)$ at different times does not commute with each
other. Here, we introduce a time-dependent perturbation theory (TDPT) in
Floquet systems to obtain the Floquet Hamiltonian for small driving amplitudes
Soori and Sen (2010); Sen  _et al._ (2021).

Consider a Hamiltonian $\hat{H}(t)=\hat{H}_{s}+\lambda\hat{H}_{V}(t)$, where
$\hat{H}_{s}$ is time independent and solvable with eigenequation
$\hat{H}_{s}|n\rangle=E_{n}|n\rangle$, $\hat{H}_{V}(t)$ is time periodic with
frequency $\omega$ and period $T=2\pi/\omega$, and $\lambda\ll 1$. We now find
solution to the time-dependent Schrödinger equation (TDSE)
$i\partial_{t}\left|{\Psi(t)}\right\rangle=\hat{H}(t)\left|{\Psi(t)}\right\rangle$
with nonzero $\lambda$, where $\left|{\Psi(t)}\right\rangle$ should satisfy

$|\Psi(T)\rangle=\hat{F}|\Psi(0)\rangle,$ (5)

according to the definition of the Floquet operator $\hat{F}$. When
$\lambda=0$, $|n(t)\rangle\equiv e^{-iE_{n}t}|n\rangle$ satisfies the TDSE
$i\partial_{t}|n(t)\rangle=\hat{H}_{s}|n(t)\rangle$. When $0<\lambda\ll 1$, we
expand $|\Psi(t)\rangle$ in the basis of $|\\{n(t)\rangle\\}$

$|\Psi(t)\rangle=\sum_{n}c_{n}(t)e^{-iE_{n}t}|n\rangle.$ (6)

Substituting it to the TDSE yields

$c_{n}(t)=c_{n}(0)-i\lambda\sum_{l}\int_{0}^{t}c_{l}(t^{\prime})e^{-i\Delta_{l,n}t^{\prime}}V_{n,l}(t^{\prime})dt^{\prime},$
(7)

where $\Delta_{l,n}=E_{l}-E_{n}$ is the energy difference between the states
$|l\rangle$ and $|n\rangle$, and $V_{n,l}(t)=\langle
n|\hat{H}_{V}(t)|l\rangle$ is the matrix elements of $\hat{H}_{V}(t)$ in the
basis $\\{|n\rangle\\}$. When $\lambda\ll 1$, we expand $c_{n}(t)$ in powers
of $\lambda$ as

$c_{n}(t)=c_{n}^{(0)}+\lambda
c_{n}^{(1)}(t)+\lambda^{2}c_{n}^{(2)}(t)+\mathcal{O}(\lambda^{3}).$ (8)

Through substituting it into Eq.(7) and comparing the coefficients before the
$k$th order of $\lambda$ of both sides of the equation, we obtain the series
of coefficients

$\displaystyle
c_{n}^{(1)}(t)=c_{n}^{(1)}(0)-i\sum_{l}\int_{0}^{t}c_{l}^{(0)}e^{-i\Delta_{l,n}t_{1}}V_{n,l}(t_{1})dt_{1},$
$\displaystyle
c_{n}^{(2)}(t)=c_{n}^{(2)}(0)-i\sum_{l}\int_{0}^{t}c_{l}^{(1)}(t_{1})e^{-i\Delta_{l,n}t_{1}}V_{n,l}(t_{1})dt_{1}$
$\displaystyle\qquad\quad=c_{n}^{(2)}(0)-i\sum_{l}\int_{0}^{t}c_{l}^{(1)}(0)e^{-i\Delta_{l,n}t_{1}}V_{n,l}(t_{1})dt_{1}-\sum_{l,m_{1}}\int_{0}^{t}\left[\int_{0}^{t_{1}}c_{l}^{(0)}e^{-i\Delta_{l,m_{1}}t_{2}}V_{m_{1},l}(t_{2})dt_{2}\right]e^{-i\Delta_{m_{1},n}t_{1}}V_{n,m_{1}}(t_{1})dt_{1},$
$\displaystyle\cdots.$ (9)

Thus,

$\displaystyle c_{n}(t)$ $\displaystyle=$ $\displaystyle c_{n}^{(0)}+\lambda
c_{n}^{(1)}(0)+\lambda^{2}c_{n}^{(2)}(0)-i\lambda\sum_{l}\int_{0}^{t}\left[c_{l}^{(0)}+\lambda
c_{l}^{(1)}(0)\right]e^{-i\Delta_{l,n}t_{1}}V_{n,l}(t_{1})dt_{1}$ (10)
$\displaystyle\qquad-\lambda^{2}\sum_{l,m_{1}}\int_{0}^{t}\left[\int_{0}^{t_{1}}c_{l}^{(0)}e^{-i\Delta_{l,m_{1}}t_{2}}V_{m_{1},l}(t_{2})dt_{2}\right]e^{-i\Delta_{m_{1},n}t_{1}}V_{n,m_{1}}(t_{1})dt_{1}+\mathcal{O}(\lambda^{3})$
$\displaystyle=$ $\displaystyle
c_{n}(0)-i\lambda\sum_{l}c_{l}(0)\int_{0}^{t}e^{-i\Delta_{l,n}t_{1}}V_{n,l}(t_{1})dt_{1}$
$\displaystyle\qquad\quad-\lambda^{2}\sum_{l,m_{1}}c_{l}(0)\int_{0}^{t}\left[\int_{0}^{t_{1}}e^{-i\Delta_{l,m_{1}}t_{2}}V_{m_{1},l}(t_{2})dt_{2}\right]e^{-i\Delta_{m_{1},n}t_{1}}V_{n,m_{1}}(t_{1})dt_{1}+\mathcal{O}(\lambda^{3}).$

Labeling the column vector $[c_{1}(t);c_{2}(t);\cdots;c_{\mathcal{D}}(t)]$ as
$|c(t)\rangle$ with $\mathcal{D}$ being the total Hilbert space dimension,
Eq.(10) at $t=T$ can be written in the matrix form

$|c(T)\rangle=[\hat{I}+\lambda\hat{F}_{1}+\lambda^{2}\hat{F}_{2}+\mathcal{O}(\lambda^{3})]|c(0)\rangle,$
(11)

where $\hat{I}$ is the identity matrix and

$\displaystyle\hat{F}_{1}$ $\displaystyle=$
$\displaystyle-i\sum_{n,l}\int_{0}^{T}e^{-i\Delta_{l,n}t_{1}}V_{n,l}(t_{1})dt_{1}|n\rangle\langle
l|,$ $\displaystyle\hat{F}_{2}$ $\displaystyle=$
$\displaystyle-\sum_{n,l}\sum_{m_{1}}\int_{0}^{T}e^{-i\Delta_{m_{1},n}t_{1}}V_{n,m_{1}}(t_{1})\left[\int_{0}^{t_{1}}e^{-i\Delta_{l,m_{1}}t_{2}}V_{m_{1},l}(t_{2})dt_{2}\right]dt_{1}|n\rangle\langle
l|.$ (12)

Writing Eq.(5) in the matrix form

$e^{-i\hat{E}T}|c(T)\rangle=\hat{F}|c(0)\rangle,$ (13)

where $e^{-i\hat{E}T}=\mathrm{diag}([e^{-iE_{1}T}\ e^{-iE_{2}T}\ \cdots\
e^{-iE_{\mathcal{D}}T}])$, one can directly write down the Floquet operator

$\hat{F}=e^{-i\hat{E}T}\left[\hat{I}+\lambda\hat{F}_{1}+\lambda^{2}\hat{F}_{2}+\mathcal{O}(\lambda^{3})\right].$
(14)

The $k$th order Floquet operator $\lambda^{k}\hat{F}_{k}$ can be viewed as the
transition between levels $l$ and $n$ intermediated by $k$ virtual transitions
assisted by the driving. For each transition between levels $m_{i}$ and
$m_{i+1}$, expanding the periodic matrix element in the Fourier series as
$V_{m_{i+1},m_{i}}(t)=\sum_{q_{i}}f_{q_{i}}(m_{i},m_{i+1})e^{iq_{i}\omega t}$
($q_{i}$ being integer numbers), the driving provides multiple frequencies
$q_{i}\omega$ if $f_{q_{i}}$ is nonzero. If there exists a $\tilde{q}_{i}$
making $\tilde{q}_{i}\omega=\Delta_{m_{i},m_{i+1}}$, single-photon resonance
occurs between $m_{i}$ and $m_{i+1}$. For a $p$ order transition between
levels $m_{i}$ and $m_{i+p}$, if all intermediate single-photon processes are
non-resonant but there exists a set of
$\\{\tilde{q}_{i},\tilde{q}_{i+1},\cdots,\tilde{q}_{i+p-1}\\}$ making
$\Delta_{m_{i},m_{i+p}}=(\tilde{q}_{i}+\tilde{q}_{i+1}+\cdots+\tilde{q}_{i+p-1})\omega$,
$p$-photon resonance between $m_{i}$ and $m_{i+p}$ occurs. If $\hat{F}_{k}$
contains $q_{r}(0\leq q_{r}\leq k)$ resonance processes, no matter single-
photon or multi-photon resonance, the integrals over time lead to a factor
$T^{q_{r}}$. Thus, in a tight manner, if $\lambda T\ll 1$, one can expand the
Floquet Hamiltonian in the Taylor series as

$\displaystyle\hat{H}_{F}=\frac{i}{T}\ln\left\\{e^{-i\hat{E}T}\left[\hat{I}+\lambda\hat{F}_{1}+\lambda^{2}\hat{F}_{2}+\mathcal{O}(\lambda^{3})\right]\right\\}$
$\displaystyle=\hat{H}_{\omega}+\frac{i}{T}\left[\lambda\hat{F}_{1}+\lambda^{2}\hat{F}_{2}-\frac{\lambda^{2}}{2}\hat{F}_{1}^{2}+\mathcal{O}(\lambda^{3}T^{3})\right],$
(15)

where $\hat{H}_{\omega}=\frac{i}{T}\ln
e^{-i\hat{E}T}=\mathrm{diag}\\{\mathrm{mod}([E_{1}\ E_{2}\ \cdots
E_{\mathcal{D}}],\omega)\\}$ is the effective on-site potential by folding the
unperturbed energies into one Floquet Brilloun zone with width $\omega$.

For our model, we treat $\lambda\hat{H}_{V}(t)=[1+u(t)]\hat{H}_{J}$ as
perturbation to
$\hat{H}_{s}=U\sum_{j=0}^{L-2}\hat{n}_{j}\hat{n}_{j+1}-g\sum_{j=0}^{L-1}j\hat{n}_{j}$,
assuming $J,u\ll g,U,\omega/2\pi$, and $g$ and $U$ are comparable. The
eigenstates of $\hat{H}_{s}$ are the Fock states
$|\vec{n}\rangle=|n_{0}n_{1}\cdots n_{L-2}n_{L-1}\rangle$ with eigenergy
$E_{\vec{n}}=U\sum_{j=0}^{L-2}n_{j}n_{j+1}-g\sum_{j=0}^{L-1}jn_{j}$, where
$n_{j}=0,1$ denotes the particle number on site $j$. We consider perturbation
to first order. The nonzero tunneling matrix element between $|\vec{n}\rangle$
and $|\vec{l}\rangle$ can be expanded in the Fourier series as
$V_{\vec{n},\vec{l}}(t)=J(1+0.5ue^{i\omega t}+0.5ue^{-i\omega t})$. Thus,
resonance occurs if $\Delta_{\vec{l},\vec{n}}=0$ and $\pm\omega$. In our
model, there are three different pairs of conjugate tunneling processes: a
particle tunnels to the right/left site with the total interaction strength
unchanged, increased/decreased and decreased/increased by $U$, which cause
energy difference $|\Delta|=g,|g-U|$ and $|g+U|$ respectively. Thus, we can
tune the relation between the interaction, tilting strength and the driving
frequency to make particular tunneling process resonant.

To reduce flexibility, we fix $U=g$, which ensures the intrinsic resonant
tunneling without driving. The first order Floquet Hamiltonian in this case
reads

$\displaystyle\hat{H}_{F}^{(1)}$ $\displaystyle=$ $\displaystyle\frac{i}{T}\ln
e^{-i\hat{E}T}+J\sum_{j}\hat{\mathcal{P}}^{(0)}_{j-1,j+2}\hat{c}^{\dagger}_{j+1}\hat{c}_{j}+h.c.$
(16) $\displaystyle+$
$\displaystyle\left\\{\frac{iJ}{T}\left(e^{-igT}-1\right)\left[\frac{1}{g}+\frac{u\left(1-\delta_{g,\omega}\right)}{2(g-\omega)}+\frac{u}{2\left(g+\omega\right)}\right]+\frac{uJ}{2}\delta_{g,\omega}\right\\}\sum_{j}\hat{P}^{(g)}_{j-1,j+2}\hat{c}^{\dagger}_{j+1}\hat{c}_{j}+h.c.$
$\displaystyle+$
$\displaystyle\left\\{\frac{iJ}{T}\left(e^{-i2gT}-1\right)\left[\frac{1}{2g}+\frac{u\left(1-\delta_{2g,\omega}\right)}{2(2g-\omega)}+\frac{u}{2\left(2g+\omega\right)}\right]+\frac{uJ}{2}\delta_{2g,\omega}\right\\}\sum_{j}\hat{P}^{(2g)}_{j-1,j+2}\hat{c}^{\dagger}_{j+1}\hat{c}_{j}+h.c.,$

where

$\displaystyle\hat{\mathcal{P}}_{j-1,j+2}^{(0)}$ $\displaystyle=$
$\displaystyle\hat{n}_{j+2}(1-\hat{n}_{j-1}),$
$\displaystyle\hat{\mathcal{P}}_{j-1,j+2}^{(g)}$ $\displaystyle=$
$\displaystyle 1-(\hat{n}_{j-1}-\hat{n}_{j+2})^{2},$
$\displaystyle\hat{\mathcal{P}}_{j-1,j+2}^{(2g)}$ $\displaystyle=$
$\displaystyle\hat{n}_{j-1}\left(1-\hat{n}_{j+2}\right)$ (17)

accounts for the tunneling processes with energy difference being 0, $g$ and
$2g$ respectively. Apart from the intrinsic resonance proportional to
$\hat{\mathcal{P}}^{(0)}_{j-1,j+2}$, there can be two partial resonances
between the driving and the system at frequencies $\omega_{1}=g$ and
$\omega_{2}=2g$ respectively, where the tunneling processes proportional to
$\hat{\mathcal{P}}_{j-1,j+2}^{(g)}$ and $\hat{\mathcal{P}}_{j-1,j+2}^{(2g)}$
is resonant correspondingly. At other frequencies, the system does not
resonate with the driving. Particularly, at $\omega=g/q$ (with $q$ being
integers larger than 1), the driving does not assist the tunneling and the two
tunneling processes vanish. In the following study, we show that the HSF can
be tuned by changing the driving frequency to and away from the resonance
frequencies.

##  III Floquet Engineering of Hilbert space fragmentation

In this section, we explore the tuning of HSF by controlling the driving
frequency within the half-filling sector. We explore the splitting of the
Hilbert space based on $\hat{H}_{F}^{(1)}$, and supplement the results by the
exact dynamical evolution governed by the evolution operator
$\hat{U}(t)=\exp[-i\hat{H}_{0}t]$ in the absence of driving, and the Floquet
operator $\hat{F}=\mathcal{T}\int_{0}^{T}\mathrm{exp}[-i\hat{H}(t)t]dt$ in the
presence of driving. The system size is set as $L=16$ unless otherwise
specified.

Some of the key diagnostics used to identify the HSF include the saturated
values of bipartite entanglement entropy (EE), infinite-temperature
autocorrelation functions and local observables.

The bipartite von Neumann EE between a subsystem $A$ and the rest of the
system $B$ is defined as

$S=-\mathrm{Tr}\left(\hat{\rho}_{A}\ln\hat{\rho}_{A}\right),$ (18)

where $\hat{\rho}_{A}=\mathrm{Tr}_{B}(|\psi\rangle\langle\psi|)$ is the
reduced density matrix of the subsystem $A$. When evolving from a low-
entangled initial state $|\psi_{0}\rangle$ in a finite system, the EE will
saturate to a value $S_{\infty}$. If the Krylov subspace to which the initial
state belongs is ergodic, $S_{\infty}$ is consistent with the Page value of
that subspace Page (1993); Yang  _et al._ (2020); Ghosh  _et al._ (2023). For
a system without HSF, the Krylov subspace is trivially the symmetry sector.
The Page value $S_{p}$ of the half-filling sector can be calculated as

$S_{p}=-\sum_{n_{A}=0}^{N}\frac{d_{A}(n_{A})d_{B}(N-n_{A})}{\mathcal{D}}\ln\frac{d_{B}(N-n_{A})}{\mathcal{D}}-\frac{1}{2},$
(19)

where $\mathcal{D}$ is the dimension of the half-filling sector and
$d_{A}(n_{A})[d_{B}(N-n_{A})]$ is the dimension of the subsystem $A[B]$ with
particle number $n_{A}[N-n_{A}]$ Ghosh  _et al._ (2023). For a system with
HSF, the Page value $S_{p}[\mathcal{K}]$ in a Krylov subspace $\mathcal{K}$
can be calculated via averaging the von Neumann EE of the random canonical
states in $\mathcal{K}$

$|\phi\rangle=\sum_{n=1}^{\mathcal{D}_{\mathcal{K}}}\frac{z_{n}}{\sqrt{\mathcal{N}}}|n\rangle,$
(20)

where $\\{|n\rangle\\}$ is the set of basis in $\mathcal{K}$,
$\mathcal{D}_{\mathcal{K}}$ is the dimension of $\mathcal{K}$, $z_{n}$ are
normally distributed real random numbers with mean 0 and variance 1 and
$\sqrt{\mathcal{N}}$ is the normalized coefficient Ghosh  _et al._ (2023). In
the following study, we consider the half-chain von Neumann EE with $A$ being
the left half chain of the system, and calculate $S_{p}[\mathcal{K}]$ in a
Krylov subspace $\mathcal{K}$ by averaging the EE over
$\mathcal{D}_{\mathcal{K}}$ random canonical states.

The infinite-temperature autocorrelation function of an operator $\hat{O}$ is
defined as
$C(t)=\langle\psi_{\mathrm{inf}}|\hat{O}(t)\hat{O}(0)|\psi_{\mathrm{inf}}\rangle$,
where $|\psi_{\mathrm{inf}}\rangle$ is a random infinite-temperature state
distributing evenly on the whole Hilbert space. It captures the spreading of
operator $\hat{O}$ and the symmetry and transport properties of the system
Moudgalya  _et al._ (2022). In the following study, we consider the infinite-
temperature density autocorrelation function at site $j$:

$C_{j}(t)=\langle\psi_{\mathrm{inf}}|[\hat{n}_{j}(t)-1/2][\hat{n}_{j}(0)-1/2]|\psi_{\mathrm{inf}}\rangle.$
(21)

At long times, $C_{j}(t)$ is larger than a lower bound Sala  _et al._ (2020);
Moudgalya  _et al._ (2022); Ghosh  _et al._ (2023)

$C_{j}^{(f)}=\frac{1}{\mathcal{D}}\sum_{\mathcal{K}}\frac{\left[\mathrm{Tr}(\hat{P}_{\mathcal{K}}\hat{n}_{j}\hat{P}_{\mathcal{K}}-\frac{1}{2})\right]^{2}}{\mathcal{D}_{\mathcal{K}}},$
(22)

where $\hat{P}_{\mathcal{K}}$ is the projector to the Krylov subspace
$\mathcal{K}$ and the summation sums over all Krylov subspaces. In the absence
of HSF, $\mathcal{K}$ equals the symmetry sector and $C_{j}^{(f)}=0$ for the
half-filling case. In the presence of HSF, $C_{j}^{(f)}$ provides a finite
lower bound for $C_{j}(t)$ Sala  _et al._ (2020); Moudgalya  _et al._ (2022);
Ghosh  _et al._ (2023).

When the Hilbert space of a system is fragmented, one can redefine ergodicity
within the Krylov subspaces Moudgalya  _et al._ (2021); Herviou  _et al._
(2021). For an ergodic system without HSF, when evolving from an out-of-
equilibrium state $|\psi_{0}\rangle$, the local observables should relax to
the equilibrium values predicted by the ensemble of the symmetry sector. In
contrast, for a system with HSF, they should relax to the equilibrium values
predicted by the ensemble of the Krylov subspace which the initial state
belongs to, if the Krylov subspace is ergodic. In the following, we study the
evolution of the local densities $n_{j}(t)=\langle\hat{n}_{j}\rangle$. For our
system, $n_{j}$ will relax to 0.5 if there is no HSF. Otherwise, it will relax
to the value

$n_{j}[{\mathcal{K}}]=\frac{1}{\mathcal{D}_{\mathcal{K}}}\mathrm{Tr}(\hat{P}_{\mathcal{K}}\hat{n}_{j}\hat{P}_{\mathcal{K}}),$
(23)

where $\mathcal{K}$ is the Krylov subspace which the initial state belongs to.

###  III.1 Strong Hilbert space fragmentation without driving

When $u=0$, the Hamiltonian of the system reduces to $\hat{H}_{0}$ and the
effective Hamiltonian up to first order (using Schrieffer-Wolff transformation
Bravyi  _et al._ (2011)) in the large tilting limit is given by

$\hat{H}_{u=0}^{\mathrm{eff}}=J\sum_{j}\hat{\mathcal{P}}_{j-1,j+2}^{(0)}(\hat{c}^{\dagger}_{j}\hat{c}_{j+1}+\hat{c}^{\dagger}_{j+1}\hat{c}_{j}),$
(24)

where only the intrinsic resonance tunneling exists. This effective
Hamiltonian conserves the sum of the dipole moment and the number of pairs of
occupied adjacent sites
$\hat{E}=-\sum_{j}j\hat{n}_{j}+\sum_{j}\hat{n}_{j}\hat{n}_{j+1}$, see Fig.
2(a).

Figure 2:  (a) Schematic of the allowed tunneling processes in the Hamiltonian
$\hat{H}^{\mathrm{eff}}_{u=0}$, which conserve $\hat{E}$. (b) Schematic of the
HSF in the half-filling sector under the action of
$\hat{H}^{\mathrm{eff}}_{u=0}$. There are $N^{2}$ symmetry spaces
$\mathcal{S}^{(0)}_{q}$ with good quantum number $E_{q}$, denoted by the empty
blue squares. Apart from
$\mathcal{S}^{(0)}_{0},\mathcal{S}^{(0)}_{1},\mathcal{S}^{(0)}_{N^{2}-2}$ and
$\mathcal{S}^{(0)}_{N^{2}}$, the other symmetry spaces are further split into
Krylov subspaces $\mathcal{K}^{(0)}_{i}$, denoted by the filled orange
squares. The two largest Krylov subspaces $\mathcal{K}^{(0)}_{\mathrm{cdw1}}$
and $\mathcal{K}^{(0)}_{\mathrm{cdw2}}$ lie in the symmetry spaces
$\mathcal{S}^{(0)}_{q_{\mathrm{cdw1}}}$ and
$\mathcal{S}^{(0)}_{q_{\mathrm{cdw2}}}$ respectively. Inset in (b): the ratio
between the dimension of the largest Krylov subspaces and that of the half-
filling sector, as a function of the system size. The blue line is the fitting
function.

The action of Hamiltonian(24) fragments the Hilbert space beyond the simple
conservation of $\hat{E}$ to Krylov subspaces $\mathcal{K}^{(0)}_{i}$, with
the superscript denoting the splitting way under Hamiltonian(24) and the
subscript denoting the different Krylov subspaces, see Fig. 2(b). There are
totally $N^{2}$ symmetry spaces $\mathcal{S}^{(0)}_{q}$ in the half-filling
sector with good quantum numbers $E_{q}=E_{0}+q$, where
$E_{0}=\frac{-3N(N-1)}{2}-1$ and $q=0,1,2,\cdots N^{2}-2$ and $N^{2}$. Apart
from the first and last two symmetry spaces (with dimensions 2, 3, 1, 1 for
all system sizes), the Hilbert space in all the other symmetry spaces further
fractures into disconnected Krylov subspaces. There are two largest Krylov
subspaces with the same dimension, labelled as
$\mathcal{K}_{\mathrm{cdw1}}^{(0)}$ and $\mathcal{K}_{\mathrm{cdw2}}^{(0)}$,
which respectively contain the charge density wave states
$|\mathrm{CDW1}\rangle=|0101\cdots 01\rangle$ and
$|\mathrm{CDW2}\rangle=|1010\cdots 10\rangle$, and lie in the symmetry spaces
$\mathcal{S}_{q_{\mathrm{cdw1}}}^{(0)}$ with $E_{q_{\mathrm{cdw1}}}=-N^{2}$
and $\mathcal{S}_{q_{\mathrm{cdw1}}}^{(0)}$ with
$E_{q_{\mathrm{cdw2}}}=-N(N-1)$. In the inset of Fig. 2(b), we show
$\max[\mathcal{D}_{\mathcal{K}^{(0)}}]/\mathcal{D}$, the dimension ratio
between the largest Krylov subspace and the half-filling sector, as a function
of the system size. This ratio decays exponentially with $L$ and tends to 0 in
the thermodynamic limit, which implies strong HSF regarding the half-filling
sector. We also calculate the dimension ratio between the largest Krylov
subspace $\mathcal{K}_{\mathrm{cdw1}}^{(0)}$
[$\mathcal{K}_{\mathrm{cdw2}}^{(0)}$] and the corresponding symmetry space
$\mathcal{S}_{q_{\mathrm{cdw1}}}^{(0)}$
[$\mathcal{S}_{q_{\mathrm{cdw2}}}^{(0)}$], which also decays exponentially
with the system size and tends to 0 in the thermodynamic limit and implies
strong HSF regarding $\mathcal{S}_{q_{\mathrm{cdw1}}}^{(0)}$
[$\mathcal{S}_{q_{\mathrm{cdw2}}}^{(0)}$].

To demonstrate the HSF of $\hat{H}_{0}$ in the large tilting limit, we study
the growth of EE starting from Fock states by the exact time evolution.
According to the random states(20) in different spaces, we can calculate the
Page value of EE, $S_{p}[\mathcal{S}^{(0)}_{q}]$, in the symmetry space
$\mathcal{S}^{(0)}_{q}$, and $S_{p}[\mathcal{K}^{(0)}_{i}]$, in the Krylov
subspace $\mathcal{K}^{(0)}_{i}$. We calculate the growth of EE starting from
all Fock states in $\mathcal{K}_{\mathrm{cdw1}}^{(0)}$, and show the average
over all the initial states at different values of $g$ in Fig. 3(a). For small
values of $g$, $S(t)$ saturates to $S_{p}$ but not
$S_{p}[{\mathcal{S}_{q_{\mathrm{cdw1}}}^{(0)}}]$, which indicates that the
system is not fragmented and thermalizes within the half-filling sector. For a
larger $g$, $S(t)$ saturates to $S_{p}[\mathcal{K}_{\mathrm{cdw1}}^{(0)}]$,
which demonstrates the splitting of Hilbert space under the action of
$\hat{H}_{u=0}^{\mathrm{eff}}$. The results in
$\mathcal{K}_{\mathrm{cdw2}}^{(0)}$ are qualitatively the same, which are not
shown here.

The HSF at large tilting strength is further evidenced by the saturated
density profile. In Fig. 3(c), we plot the saturated density profile
$\bar{n}_{j}$ at different values of $g$. It is obtained by averaging
$n_{j}(t)$ over the times $t\in[700,800]\frac{2\pi}{g}$, during which
$n_{j}(t)$ already saturates, and over all initial Fock states in
$\mathcal{K}_{\mathrm{cdw1}}^{(0)}$. For comparison, we also plot the density
profile $n_{j}[\mathcal{K}^{(0)}_{\mathrm{cdw1}}]$, predicted by the infinite-
temperature ensemble of the Krylov subspace
$\mathcal{K}^{(0)}_{\mathrm{cdw1}}$. For small values of $g$, $\bar{n}_{j}$
concentrate on 0.5 for all lattice sites, which is consistent with the result
of thermalization within the half-filling sector. For a larger $g$, the
density profile is consistent with
$n_{j}[{\mathcal{K}_{\mathrm{cdw1}}^{(0)}}]$, again implying that the Hilbert
space splits according to $\hat{H}_{u=0}^{\mathrm{eff}}$. The results in
$\mathcal{K}_{\mathrm{cdw2}}^{(0)}$ are qualitatively the same, which are not
shown here.

Figure 3:  (a) The growth of the normalized von Neumann EE $S(t)/S_{p}$ for
different values of $g$, averaged over all initial states (34 in total) in the
largest Krylov subspace $\mathcal{K}_{\mathrm{cdw1}}^{(0)}$. The black dotted
line and blue dashed lines denote
$S_{p}[\mathcal{S}_{q_{\mathrm{cdw1}}}^{(0)}]/S_{p}$ and
$S_{p}[\mathcal{K}_{\mathrm{cdw1}}^{(0)}]/S_{p}$, respectively. (b) At $g=50$,
the fidelity dynamics from the frozen state
$|\psi_{\mathrm{frozen}}\rangle=|1111111010000000\rangle$ (blue line), which
lies in $\mathcal{S}^{(0)}_{N^{2}-2}$, and the domain state
$|\psi_{\mathrm{d}}\rangle=|0000000011111111\rangle$ (red line), which lies in
$\mathcal{S}^{(0)}_{0}$. Red dotted line shows the amplitude of state transfer
between $|\psi_{\mathrm{d}}\rangle$ and
$|\bar{\psi}_{\mathrm{d}}\rangle=|0000000101111111\rangle$. (c) The saturated
density profile $\bar{n}_{j}$ for different values of $g$, averaged over time
$t\in[700,800]\frac{2\pi}{g}$ and all initial Fock states in
$\mathcal{K}_{\mathrm{cdw1}}^{(0)}$. The blue dashed line denotes
$n_{j}[\mathcal{K}_{\mathrm{cdw1}}^{(0)}]$. (d) The evolution of $C_{L}(t)$
for different values of $g$, starting from a random infinite-temperature
state. The blue dashed line denotes the lower bound $C_{L}^{(0)}=0.25$,
predicted by the HSF due to Hamiltonian(24). In (a), (c) and (d), the data
with $g=15$ and $50$ almost collapse with each other. The other parameters are
$u=0$ and $L=16$ in all plots. The energies are scaled in units of $J$.

Apart from the large Krylov subspaces, there are numerous frozen states and
small Krylov subspaces. The frozen states correspond to the Fock states which
are zero-energy eigenstates of $\hat{H}_{u=0}^{\mathrm{eff}}$. Any Fock states
are frozen under the action of $\hat{H}_{u=0}^{\mathrm{eff}}$, if they lack
both the configurations $``\cdots 0011\cdots"$ and $``\cdots 0101\cdots"$,
where the sequences $011$ and $101$ are away from the left particles by at
least one site. In Fig. 3(b), we plot the evolution of the fidelity
$|\langle\psi_{0}|\hat{U}(t)|\psi_{0}\rangle|^{2}$ from a frozen state
$|\psi_{0}\rangle=|\psi_{\mathrm{frozen}}\rangle$ at $g=50$. The fidelity
keeps near 1 for all the time considered, evidencing that the state stays to
the initial state. Furthermore, we study the dynamics in a Krylov subspace
with two elements, which contains a state $|\psi_{\mathrm{d}}\rangle$ and its
partner
$|\bar{\psi}_{\mathrm{d}}\rangle=\hat{H}_{u=0}^{\mathrm{eff}}|\psi_{\mathrm{d}}\rangle$.
In Fig. 3(b), we plot the fidelity dynamics from
$|\psi_{0}\rangle=|\psi_{\mathrm{d}}\rangle$ at $g=50$. The fidelity
oscillates periodically, and the state transfers almost perfectly between
$|\psi_{\mathrm{d}}\rangle$ and $|\bar{\psi}_{\mathrm{d}}\rangle$, see the red
solid and dotted lines in Fig. 3(b). The oscillation period is
$\tilde{T}_{2}=\pi/J$, which is consistent with the prediction by projecting
the effective Hamiltonian $\hat{H}_{u=0}^{\mathrm{eff}}$ into the Krylov
subspace.

At last, we study the autocorrelation function $C_{j=L}(t)$ starting from a
random infinite-temperature state. The results make little difference for
different infinite-temperature states. The lower bound predicted by the
present HSF is $C_{L}^{(0)}=0.25$. In Fig. 3(d), we plot the evolution of
$C_{L}(t)$ for different values of $g$. For small values of $g$, the
perturbation theory does not work, and $C_{L}(t)$ decays to zero quickly,
which is the result of thermalization in the half-filling sector. As $g$
increases, $C_{L}(t)$ saturates to finite values, and tends to the bound
$C_{L}^{(0)}$. All these results show that the half-filling sector is split
strongly according to $\hat{H}^{\mathrm{eff}}_{u=0}$ in the large tilting
limit and without driving.

###  III.2 Strong Hilbert space fragmentation at the partial resonance
frequency $\omega_{1}=g$

When the driving frequency $\omega=\omega_{1}$, the tunneling processes
projected by $\hat{P}^{(0)}_{j-1,j+2}$ and $\hat{P}^{(g)}_{j-1,j+2}$ are
allowed, while the one projected by $\hat{P}^{(2g)}_{j-1,j+2}$ vanishes, see
schematic in Fig. 4(a). The first order effective Floquet Hamiltonian in this
case reads

$\displaystyle\hat{H}_{\omega_{1}}^{\mathrm{eff}}$ $\displaystyle=$
$\displaystyle
J\sum_{j}\hat{\mathcal{P}}_{j-1,j+2}^{(0)}(\hat{c}^{\dagger}_{j}\hat{c}_{j+1}+h.c.)$
(25) $\displaystyle+$
$\displaystyle\frac{uJ}{2}\sum_{j}\hat{\mathcal{P}}_{j-1,j+2}^{(g)}(\hat{c}^{\dagger}_{j}\hat{c}_{j+1}+h.c.).$

The additional tunneling term relaxes some kinetic constraints, thus changes
the structure of the Hilbert space. The Hamiltonian(25) does not show obvious
symmetry, and splits the half-filling sector into $N$ Krylov subspaces
$\mathcal{K}_{i}^{(1)}$ $(i=1,\cdots,N)$, with the superscript denoting the
splitting way under $\hat{H}_{\omega_{1}}^{\mathrm{eff}}$. Each Krylov
subspace $\mathcal{K}_{i}^{(1)}$ contains the state
$|\psi_{i}\rangle=|(0\cdots 0)^{i-1}(1\cdots 1)^{N}(0\cdots 0)^{N-i+1}\rangle$
with $(0\cdots 0)^{q}$ and $(1\cdots 1)^{q}$ denoting $q$ contiguous empty
sites and occupied sites respectively. The dimension of
$\mathcal{K}_{i}^{(1)}$ increases with $i$ with $\mathcal{K}_{1}^{(1)}$ being
the minimal subspace with dimension 1 for all system sizes, i.e., $|(1\cdots
1)^{N}(0\cdots 0)^{N}\rangle$ is a frozen state. In Fig. 4(b), we plot the
dimension ratio $\max[\mathcal{D}_{\mathcal{K}^{(1)}}]/\mathcal{D}$ between
the largest Krylov subspace and the half-filling sector as a function of the
system size. We numerically find that
$\max[\mathcal{D}_{\mathcal{K}^{(1)}}]/\mathcal{D}=\frac{8(L+1)}{(L+2)(L+4)}$,
with the analytical derivation deserving further study. This ratio means that
the largest Krylov subspace occupies a vanishing fraction of the half-filling
sector in the thermodynamic limit and the system is strongly fragmented.

Figure 4:  (a) Schematic of the allowed tunneling processes in the Hamiltonian
$\hat{H}^{\mathrm{eff}}_{\omega_{1}}$. (b) Schematic of the HSF in the half-
filling sector under the action of $\hat{H}^{\mathrm{eff}}_{\omega_{1}}$. It
is fragmented into $N$ Krylov subspaces $\mathcal{K}^{(1)}_{i}$ with the
dimension increasing with $i$, denoted by the filled orange squares. Inset in
(b): the ratio between the dimension of the largest Krylov subspace and that
of the half-filling sector, as a function of the system size. The blue circles
denote the numerical data and the blue line is the function
$y=\frac{8(L+1)}{(L+2)(L+4)}$.

To verify the strong HSF, we calculate the EE dynamics, the saturated density
profile and the density autocorrelation function dynamics by the exact time
evolution at $\omega=g$, $u=1$ and system size $L=16$. In Fig. 5(a), we plot
$S(kT)/S_{p}$ versus the driving cycles $k$ at different values of $g$ for the
largest Krylov subspace $\mathcal{K}^{(1)}_{8}$. The EE are averaged over 50
initial random Fock states in $\mathcal{K}^{(1)}_{8}$. For small values of
$g$, $S(kT)$ saturates to $S_{p}$ rapidly. As $g$ increases, $S(kT)$ starts to
saturate to smaller values, and reaches $S_{p}[\mathcal{K}^{(1)}_{8}]$
finally. At a large value of $g$, the EE for Fock states in the other Krylov
subspaces also saturates to the Page value of the corresponding Krylov
subspaces. In Fig. 5(b), we plot $S(kT)/S_{p}$ at $g=50$ in Krylov subspaces
$\mathcal{K}^{(1)}_{7}$, $\mathcal{K}^{(1)}_{5}$, $\mathcal{K}^{(1)}_{3}$ and
$\mathcal{K}^{(1)}_{1}$. For $\mathcal{K}^{(1)}_{7,5,3}$, the EE are averaged
over 50 initial random Fock states, and they saturate to values
$S_{p}[\mathcal{K}^{(1)}_{7,5,3}]$ respectively. For $\mathcal{K}^{(1)}_{1}$,
the EE keeps near 0, which is consistent with the frozen state prediction.
Fig. 5(c) shows the saturated density profile $\bar{n}_{j}$ at different
values of $g$, for the same initial states as that in Fig. 5(a). It is
obtained by averaging $n_{j}(kT)$ over the driving cycles $k\in[2900,3000]$,
during which $n_{j}(kT)$ already saturates, and over all the 50 initial
states. For small values of $g$, $\bar{n}_{j}$ shows a uniform profile around
0.5. For large values of $g$, $\bar{n}_{j}$ is consistent with
$n_{j}[\mathcal{K}^{(1)}_{8}]$. At last, Fig. 5(d) shows the evolution of
$C_{L}(kT)$ for different values of $g$, starting from the same infinite-
temperature state as that in Fig. 3(d). The results make little difference for
different infinite-temperature states. For small $g$, $C_{L}(kT)$ decays to
zero quickly. For a large $g$, $C_{L}(kT)$ saturates to $C_{L}^{(1)}$
predicted by the present HSF. All these features demonstrate that the half-
filling sector is split strongly according to
$\hat{H}^{\mathrm{eff}}_{\omega_{1}}$ at $\omega=\omega_{1}$ in the large
tilting limit.

Figure 5:  (a) $S(kT)/S_{p}$ versus the driving cycles $k$ for different
values of $g$, averaged over 50 randomly chosen Fock states in the largest
Krylov subspace $\mathcal{K}_{8}^{(1)}$. The blue dashed line denotes
$S_{p}[{\mathcal{K}_{8}^{(1)}}]/S_{p}$. (b) $S(kT)/S_{p}$ versus the driving
cycles $k$ at $g=50$, averaged over 50 randomly chosen Fock states in
$\mathcal{K}_{7}^{(1)}$ (blue line), $\mathcal{K}_{5}^{(5)}$ (red line) and
$\mathcal{K}_{3}^{(1)}$ (yellow line), and starting from the Fock state in
$\mathcal{K}_{1}^{(1)}$ (purple line). The dashed lines with different colors
denote $S_{p}[\mathcal{K}_{i}^{(1)}]/S_{p}$ ($i=7,5,3$ and 1) of the
corresponding Krylov subspaces $\mathcal{K}_{i}^{(1)}$. For
$\mathcal{K}_{1}^{(1)}$, the Fock state is frozen and $S(kT)$ stays near 0.
(c) The saturated density profile $\bar{n}_{j}$ for different values of $g$,
averaged over the driving cycles $k\in[2900,3000]$ and 50 randomly chosen Fock
states in $\mathcal{K}_{8}^{(1)}$. The blue dashed line denotes
$n_{j}[\mathcal{K}_{8}^{(1)}]$. (d) $C_{L}(kT)$ versus the driving cycles $k$
for different values of $g$, starting from the same infinite-temperature state
in Fig. 3(d). The blue dashed line denotes the lower bound $C_{L}^{(1)}$
predicted by the HSF due to $\hat{H}_{\omega_{1}}^{\mathrm{eff}}$. In (a), (c)
and (d), the numerical data with $g=30$ and $50$ almost collapse with each
other. The other parameters are $u=1,\omega=g$ and $L=16$ in all plots. The
energies are scaled in units of $J$.

###  III.3 Strong Hilbert space fragmentation at the partial resonance
frequencies $\omega_{2}=2g$

When the driving frequency $\omega=\omega_{2}$, the tunneling processes
projected by $\hat{P}^{(0)}$ and $\hat{P}^{(2g)}$ are resonant. The tunneling
process projected by $\hat{P}^{(g)}$ is off-resonant, but yet does not vanish.
The first order effective Floquet Hamiltonian in this case reads

$\displaystyle\hat{H}_{\omega_{2}}^{{}^{\prime}\mathrm{eff}}$ $\displaystyle=$
$\displaystyle\frac{g}{2}\left[1-(-1)^{\sum_{j}\left(-j\hat{n}_{j}+\hat{n}_{j}\hat{n}_{j+1}\right)}\right]$
(26) $\displaystyle+$ $\displaystyle
J\sum_{j}\hat{\mathcal{P}}_{j-1,j+2}^{(0)}(\hat{c}^{\dagger}_{j}\hat{c}_{j+1}+\hat{c}^{\dagger}_{j+1}\hat{c}_{j})$
$\displaystyle-$
$\displaystyle\frac{2iJ(3-u)}{3\pi}\sum_{j}\hat{\mathcal{P}}_{j-1,j+2}^{(g)}(\hat{c}^{\dagger}_{j+1}\hat{c}_{j}-\hat{c}^{\dagger}_{j}\hat{c}_{j+1})$
$\displaystyle+$
$\displaystyle\frac{uJ}{2}\sum_{j}\hat{\mathcal{P}}_{j-1,j+2}^{(2g)}(\hat{c}^{\dagger}_{j}\hat{c}_{j+1}+\hat{c}^{\dagger}_{j+1}\hat{c}_{j}).$

At first sight, all the kinetic constraints are released now. However, due to
the effective on-site potential, the tunneling process projected by
$\hat{\mathcal{P}}^{(g)}$ is biased by energy $g$, which is much larger than
the tunneling amplitude $2J(3-u)/3\pi$. Thus, to first order in $J/g$, we can
neglect this tunneling process and reduce the effective Hamiltonian to

$\displaystyle\hat{H}_{\omega_{2}}^{\mathrm{eff}}$ $\displaystyle=$
$\displaystyle
J\sum_{j}\hat{\mathcal{P}}_{j-1,j+2}^{(0)}(\hat{c}^{\dagger}_{j}\hat{c}_{j+1}+\hat{c}^{\dagger}_{j+1}\hat{c}_{j})$
(27) $\displaystyle+$
$\displaystyle\frac{uJ}{2}\sum_{j}\hat{\mathcal{P}}_{j-1,j+2}^{(2g)}(\hat{c}^{\dagger}_{j}\hat{c}_{j+1}+\hat{c}^{\dagger}_{j+1}\hat{c}_{j}),$

where the effective on-site potential has been dropped, since the two
tunneling processes conserve the parity $\hat{\mathcal{P}}_{E}=(-1)^{\hat{E}}$
of $\hat{E}=\sum_{j}\left(-j\hat{n}_{j}+\hat{n}_{j}\hat{n}_{j+1}\right)$, see
Fig. 6(a).

The action of $\hat{H}^{\mathrm{eff}}_{\omega_{2}}$ splits the half-filling
sector into Krylov subspaces $\mathcal{K}^{(2)}_{i}$ beyond the conservation
of $\hat{\mathcal{P}}_{E}$, where the superscript of $\mathcal{K}^{(2)}_{i}$
denotes the splitting way under $\hat{H}_{\omega_{2}}^{\mathrm{eff}}$ and the
subscript denotes the different Krylov subspaces, see Fig. 6(b). It should be
noted that the two allowed processes are related by a spatial-reflection
transformation, i.e.
$\hat{R}[\hat{\mathcal{P}}^{(0)}_{j-1,j+2}(\hat{c}^{\dagger}_{j}\hat{c}_{j+1}+h.c.)]\hat{R}^{-1}=\hat{\mathcal{P}}^{(2g)}_{L-j-1,L-j+2}(\hat{c}^{\dagger}_{L-j+1}\hat{c}_{L-j}+h.c.)$,
where the spatial-reflection operator $\hat{R}$ is defined by
$\hat{R}|n_{0}n_{1}\cdots n_{L-2}n_{L-1}\rangle=|n_{L-1}n_{L-2}\cdots
n_{1}n_{0}\rangle$. This leads to an odd-even effect of the particle number on
the splitting of the Hilbert space (see Appendix B for details). When the
particle number $N$ is odd, the Krylov subspaces in the even- and odd-parity
symmetry space $\mathcal{S}^{(2)}_{e}$ and $\mathcal{S}^{(2)}_{o}$ are in one-
to-one correspondence through the spatial-reflection transformation. Thus,
there are two largest Krylov subspace with the same dimension, labelled as
$\mathcal{K}^{(2)}_{e1}$ and $\mathcal{K}^{(2)}_{o1}$, which lie in
$\mathcal{S}^{(2)}_{e}$ and $\mathcal{S}^{(2)}_{o}$ respectively. When $N$ is
even, there exists reflection-invariant Krylov subspace and the largest Krylov
subspace lies in $\mathcal{S}^{(2)}_{e}$. So, the scalings of the dimension
ratio between the largest Krylov subspace and the half-filling sector is
different for quadruple and non-quadruple lattice sites, and we need two
functions to fit them. In the inset of Fig. 6(b), we plot this ratio as a
function of the system size. It shows exponential decay with the system size
and vanishes in the thermodynamic limit, which implies strong HSF regarding
the half-filling sector. We also calculate the dimension ratio between the
largest Krylov subspace $\mathcal{K}_{e1}^{(2)}$ [$\mathcal{K}_{o1}^{(2)}$]
within the even (odd)-parity symmetry space $\mathcal{S}_{e}^{(2)}$
[$\mathcal{S}_{o}^{(2)}$] and that of $\mathcal{S}_{e}^{(2)}$
[$\mathcal{S}_{o}^{(2)}$]. They also decays exponentially with the system size
and tend to 0 in the thermodynamic limit, and implies strong HSF regarding
$\mathcal{S}_{e}^{(2)}$ and $\mathcal{S}_{o}^{(2)}$.

Figure 6:  (a) Schematic of the allowed tunneling processes in the Hamiltonian
$\hat{H}^{\mathrm{eff}}_{\omega_{2}}$, which conserve the parity of $\hat{E}$.
(b) Schematic of the HSF in the half-filling sector under the action of
$\hat{H}^{\mathrm{eff}}_{\omega_{2}}$. The even and odd parity symmetry spaces
$\mathcal{S}^{(2)}_{e}$ and $\mathcal{S}^{(2)}_{o}$ are denoted by the empty
blue squares. $\mathcal{S}^{(2)}_{e}$ and $\mathcal{S}^{(2)}_{o}$ are further
fragmented into Krylov subspaces $\mathcal{K}^{(2)}_{i}$, denoted by the
filled orange squares. Inset in (b): the ratio between the dimension of the
largest Krylov subspace and that of the half-filling sector, as a function of
the system size. The blue circles denote the numerical data. The blue dashed
and dotted lines denote the fitting functions for quadruple and non-quadruple
lattice sites, respectively.  Figure 7:  (a) $S(kT)/S_{p}$ versus the driving
cycles $k$ for different values of $g$, averaged over 50 randomly chosen Fock
states in $\mathcal{K}^{(2)}_{o1}$. The blue dashed line denotes
$S_{p}[\mathcal{K}^{(2)}_{o1}]/S_{p}$. (b) $S(kT)/S_{p}$ versus the driving
cycles $k$ at $g=50$, averaged over 50 randomly chosen Fock states in
$\mathcal{K}^{(2)}_{e1}$ (blue line), $\mathcal{K}^{(2)}_{e2}$ (red line) and
$\mathcal{K}^{(2)}_{o2}$ (yellow line). The dashed lines with different colors
denote the normalized Page value of EE of the corresponding Krylov subspaces.
(c) The saturated density profile $\bar{n}_{j}$ for different values of $g$,
averaged over the driving cycles $k\in[2900,3000]$ and 50 randomly chosen Fock
states in $\mathcal{K}^{(2)}_{o1}$. The blue dashed line denotes
$n_{j}[\mathcal{K}^{(2)}_{o1}]$. (d) $C_{L}(kT)$ versus $k$ for different
values of $g$, starting from the same infinite-temperature state in Fig. 3(d).
The blue dashed line denotes the lower bound $C_{L}^{(2)}$ predicted by the
HSF due to $\hat{H}^{\mathrm{eff}}_{\omega_{2}}$. In (a), (c) and (d), the
numerical data with $g=15,30$ and $50$ almost collapse with each other. The
other parameters are $u=1,\omega=2g$ and $L=16$ in all plots. The energies are
scaled in units of $J$.

To verify this strong HSF, we calculate the EE dynamics, the saturated density
profile and the density autocorrelation function dynamics by the exact time
evolution at $\omega=2g$ and $u=1$. In Fig. 7(a), we plot $S(kT)/S_{p}$ versus
the driving cycles $k$ at different values of $g$ for the largest Krylov
subspace $\mathcal{K}^{(2)}_{o1}$ within the odd-parity space. The EE are
averaged over 50 initial random Fock states in $\mathcal{K}^{(2)}_{o1}$. For
small values of $g$, $S(kT)$ saturates to $S_{p}$ rapidly. As $g$ increases,
$S(kT)$ starts to saturate to smaller values, and reaches
$S_{p}[\mathcal{K}^{(2)}_{o1}]$ finnaly. At a large value of $g$, the EE for
Fock states in other Krylov subspaces also saturates to the Page value of the
corresponding Krylov subspaces. In Fig. 7(b), we plot $S(kT)/S_{p}$ at $g=50$
within the largest and second largest Krylov subspace $\mathcal{K}^{(2)}_{e1}$
and $\mathcal{K}^{(2)}_{e2}$ within the even-parity space, and the second
largest Krylov subspace $\mathcal{K}^{(2)}_{o2}$ within the odd-parity space.
For all the Krylov subspaces, the EE are calculated by averaging 50 initial
random Fock states, and they saturate to the Page values of the corresponding
Krylov subspaces. Fig. 7(c) shows the saturated density profile $\bar{n}_{j}$
at different values of $g$, for the same initial states as that in Fig. 7(a).
It is obtained by averaging $n_{j}(kT)$ over the driving cycles
$k\in[2900,3000]$, during which $n_{j}(kT)$ already saturates, and over all
the initial states. For small values of $g$, $\bar{n}_{j}$ shows a uniform
profile around 0.5. For large values of $g$, $\bar{n}_{j}$ is consistent with
$n_{j}[\mathcal{K}^{(2)}_{o1}]$. At last, Fig. 7(d) shows the evolution of
$C_{L}(kT)$ for different values of $g$, starting from the same infinite-
temperature state as that in Fig. 3(d). The results make little difference for
different infinite-temperature states. For small $g$, $C_{L}(kT)$ decays to
zero quickly. For a large $g$, $C_{L}(kT)$ saturates to $C_{L}^{(2)}$
predicted by the present strong HSF. All these features demonstrate the half-
filling sector is split strongly under $\hat{H}^{\mathrm{eff}}_{\omega_{2}}$
at $\omega=\omega_{2}$ in the large tilting limit.

###  III.4 Strong Hilbert space fragmentation at non-resonant frequencies

When the driving frequency $\omega$ is away from the partial resonant values
$\omega_{1}$ and $\omega_{2}$, the tunneling strength of the processes
projected by $\hat{\mathcal{P}}^{(g)}$ and $\hat{\mathcal{P}}^{(2g)}$ is much
smaller than the effective potential bias, and the two tunneling processes can
be neglected. One then expects the same HSF as the case without driving.
Specifically, when $\omega=g/q$ with $q$ being integers larger than 1, the
tunneling strength of the two processes equals 0. We demonstrate the strong
HSF at the non-resonant frequencies by studying the EE dynamics, the fidelity
dynamics, the saturated density profile and the density autocorrelation
function dynamics through the exact time evolution, see Fig. 8 for
$\omega=g/2$ and $u=1$. It should be noted that for the two-dimensional Krylov
subspace, the period of fidelity oscillation in units of the driving period
$T$ is predicted as $\tilde{T}_{2}/T=\omega/2J$. All the observables show
qualitatively similar behaviour as that in Fig. 3, thus demonstrating that in
the large tilting limit and at non-resonant frequencies, the Hilbert space is
split in the same way as that in Sec. III.1.

Figure 8:  (a) The growth of the averaged and normalized von Neumann EE
$S(kT)/S_{p}$ for different values of $g$. (b) The fidelity dynamics at
$g=50$, from the frozen state $|\psi_{\mathrm{frozen}}\rangle$ and the domain
state $|\psi_{\mathrm{d}}\rangle$. Red dotted line shows the amplitude of
state transfer between $|\psi_{\mathrm{d}}\rangle$ and
$|\bar{\psi}_{\mathrm{d}}\rangle$. (c) The saturated density profile
$\bar{n}_{j}$ for different values of $g$, averaged over the driving cylcles
$k\in[900,1000]$ and all the initial states. (d) The evolution of $C_{L}(kT)$
for different values of $g$. The initial states in (a)-(d) are the same as
that in Fig. 3(a)-(d) respectively. In (a), (c) and (d), the data with
$g=15,30$ and $50$ almost collapse with each other. The other parameters are
$u=1$, $\omega=g/2$ and $L=16$ in all plots. The energies are scaled in units
of $J$.

###  III.5 Tuning the Hilbert space fragmentation by controlling the driving
frequency at large tilting strength

Now, we have established that the Hilbert space will fragment in three
different ways when tuning the driving frequency in the large tilting limit.
To make it clearer, we study the evolution of $S(kT)$ from one Fock state
$|\psi_{0}\rangle$ and of $C_{L}(kT)$ from one infinite-temperature state.

In Fig. 9(a), we plot the growth of the normalized EE from the initial state
$|\psi_{0}\rangle=|1010011001101001\rangle$, for $\omega=g/2,g/1.5,2g$ and
$g$, $g=50$ and $u=1$. The initial state is chosen in some rand way, but
requiring that the Page values of EE of the three Krylov subspaces
$\mathcal{K}^{(0)}_{\psi_{0}}$, $\mathcal{K}^{(1)}_{\psi_{0}}$ and
$\mathcal{K}^{(2)}_{\psi_{0}}$ to which $|\psi_{0}\rangle$ belongs make
sufficient difference for better illustration. For $\omega=g$ and $2g$, the EE
saturates to $\mathcal{S}_{p}[\mathcal{K}^{(1)}_{\psi_{0}}]$ and
$\mathcal{S}_{p}[\mathcal{K}^{(2)}_{\psi_{0}}]$ respectively. For the other
two frequencies, the EE saturates to
$\mathcal{S}_{p}[\mathcal{K}^{(0)}_{\psi_{0}}]$. In Fig. 9(b), we plot the
normalized saturated EE $\bar{S}/S_{p}$ versus $g/\omega$ for
$|\psi_{0}\rangle$ at $g=50$ and $u=1$. $\bar{S}$ is calculated by averaging
$S(kT)$ over the driving cycles $k\in[1900,2000]$. We see two clear peaks at
$\omega=\omega_{1}$ and $\omega_{2}$, where the saturated EE is consistent
with $\mathcal{S}_{p}[\mathcal{K}^{(1)}_{\psi_{0}}]$ and
$\mathcal{S}_{p}[\mathcal{K}^{(2)}_{\psi_{0}}]$ respectively. At other
frequencies, the saturated EE is consistent with
$\mathcal{S}_{p}[\mathcal{K}^{(0)}_{\psi_{0}}]$.

The saturated EE depends on the initial state strongly. The saturated
autocorrelation function barely depends on the specific form of the initial
infinite-temperature state, thus serves as a better criterion to distinguish
the HSF. In Fig. 9(c), we plot the saturated density autocorrelation function
$\bar{C}_{L}$ versus $g/\omega$ from the infinite-temperature state in Fig.
3(d), at $g=50$ and $u=1$. $\bar{C}_{L}$ is obtained by averaging $C_{L}(kT)$
over the driving cycles $k\in[5900,6000]$. It shows two dips at
$\omega=\omega_{1}$ and $\omega_{2}$, which is consistent with the predicted
value $C_{L}^{(1)}$ and $C_{L}^{(2)}$ respectively. For other values of
$\omega$, $\bar{C}_{L}$ is consistent with $C_{L}^{(0)}$.

Finally, in Fig. 9(d), we plot a phase diagram on the parameter plane
$(g,g/\omega)$ at $u=1$ and $L=12$, which shows the thermal and three kinds of
HSF features of $\bar{C}_{L}$. For $\omega=\omega_{1}$ and $\omega_{2}$,
$\bar{C}_{L}$ respectively crossovers from 0 to $C_{L}^{(1)}$ and
$C_{L}^{(2)}$ as $g$ increases, while for other values of $\omega$, it
crossovers to $C_{L}^{(0)}$. All these features demonstrate that one can tune
the HSF by controlling the driving frequency when the tilting strength is
sufficiently strong.

Figure 9:  (a) $S(kT)/S_{p}$ versus the driving cycles $k$ for different
driving frequency $\omega=g/2$ (blue line), $g/1.5$ (red line), $g$ (yellow
line) and $2g$ (purple line). The initial state is
$|\psi_{0}\rangle=|1010011001101001\rangle$. (b) $\bar{S}/S_{p}$ versus
$g/\omega$, with the same initial state in (a). In (a) and (b), the blue
dotted, dashed and dot-dashed lines denote the three normalized Page values
$S_{p}[\mathcal{K}_{\psi_{0}}^{(0)}]/S_{p}$,
$S_{p}[\mathcal{K}_{\psi_{0}}^{(1)}]/S_{p}$ and
$S_{p}[\mathcal{K}_{\psi_{0}}^{(2)}]/S_{p}$ respectively. (c) $\bar{C}_{L}$
versus $g/\omega$, starting from the same infinite-temperature state in Fig.
3(d), 5(d) and 7(d). In (a)-(c), the other parameters are $u=1$, $g=50$ and
$L=16$. (d) $\bar{C}_{L}$ as a function of $(g,g/\omega)$ at $u=1$ and system
size $L=12$. In all plots, the energies are scaled in units of $J$.

##  IV Summary

In conclusion, we have proposed a novel scheme to tune the HSF within the
half-filling sector, through Floquet engineering of constrained hoppings in
strongly tilted lattices. In the absence of driving, the system is kinetically
constrained and stays in a strong fragmented phase. In the presence of
driving, we engineer two effective Hamiltonians which release some kinetic
constraints in different ways, by utilizing the partial resonance between the
driving and the system. We give the two partial resonance frequencies
analytically by using the time-dependent perturbation theory in Floquet
system. The releasing of the kinetic constraints changes the structure of the
Hilbert space and leads to another two kinds of strong Hilbert space
fragmentation. At the non-resonant frequencies, the system remains in the
strong fragmented phase without driving. We demonstrate the three kinds of
Hilbert space fragmentation by studying the EE dynamics, the saturated density
profile and the density autocorrelation functions. Our results provide an
efficient way to tune the Hilbert space fragmentation by controlling the
driving frequency of particle tunneling in tilted and interacting lattice
systems. The possibility to control the hopping channels through tuning the
resonance between the driving and the system provides a new opportunity to
study the relation between kinetic constraints and other ergodicity breaking
phenomenon, such as the quantum many body scars Serbyn  _et al._ (2021); Guo
_et al._ (2023).

###### Acknowledgements.

This study is supported by the National Natural Science Foundation of China
(Grants No. 12025509, No. 12305048) and the National Key Research and
Development Program of China (Grant No. 2022YFA1404104). Y.K. is supported by
the National Natural Science Foundation of China (Grant No. 12275365) and the
Natural Science Foundation of Guangdong (Grant No. 2023A1515012099).

##  Appendix A Proposal for the cosine driving protocol

We discuss a scheme to realize our model
$\hat{H}(t)=\hat{H}_{J}(t)+\hat{H}_{U}+\hat{H}_{g}$ in Rydberg atom platform,
where $\hat{H}_{J}(t)=J[1+u\cos(\omega
t)]\sum_{j}(\hat{c}^{\dagger}_{j}\hat{c}_{j+1}+h.c.)$,
$\hat{H}_{U}=U\sum_{j}\hat{n}_{j}\hat{n}_{j+1}$ and
$\hat{H}_{g}=-g\sum_{j}j\hat{n}_{j}$.

We consider a one-dimensional chain of $L$ Rydberg atoms, with lattice
constant $a$ and each trapped in optical tweezers Barredo  _et al._ (2016);
Endres  _et al._ (2016); Browaeys  _et al._ (2016); Barredo  _et al._ (2018);
Browaeys and Lahaye (2020). Two Rydberg states are chosen to simulate the
empty state and the fermions occupied state with $|0\rangle=|nS\rangle$ and
$|1\rangle=|n^{\prime}P\rangle$, where $|nS\rangle$ and $|n^{\prime}P\rangle$
are the two Rydberg states with principal quantum numbers $n$ and $n^{\prime}$
($n\simeq n^{\prime}$) and angular momentum $S$ and $P$ respectively (we note
that $n$ labels the principal quantum number in this appendix and the fermions
number in the main text). The Rydberg atom chain is subjected to a gradient
magnetic field which causes an effective Zeeman splitting $M_{j}$ between the
two Rydberg states with $M_{j}-M_{j+1}=\tilde{g}$, as described the
Hamiltonian
$\hat{H}_{\tilde{g}}=-\tilde{g}\sum_{j}j(|n^{\prime}P\rangle\langle
n^{\prime}P|)_{j}$.

Figure A1:  The two Raman-assisted three-order transition processes coupling
$|nS\rangle_{j}|n^{\prime}P\rangle_{j+1}$ to
$|n^{\prime}P\rangle_{j}|nS\rangle_{j+1}$. The left energy level scheme
denotes the three-order transition process
$|nS\rangle_{j}|n^{\prime}P\rangle_{j+1}\to|n^{\prime}P\rangle_{j}|nS\rangle_{j+1}\to|n^{\prime}P\rangle_{j}|\alpha\rangle_{j+1}\to|n^{\prime}P\rangle_{j}|nS\rangle_{j+1}$,
in which the first process (i) is the bare exchange interaction and the second
(ii) and third processes (iii) induces a two-photon detuning
$\omega_{21}=\omega_{\mathrm{R}_{2}}-\omega_{\mathrm{R}_{1}}$. The right
energy level scheme denotes the three-order transition process
$|nS\rangle_{j}|n^{\prime}P\rangle_{j+1}\to|\alpha\rangle_{j}|n^{\prime}P\rangle_{j+1}\to|nS\rangle_{j}|n^{\prime}P\rangle_{j+1}\to|n^{\prime}P\rangle_{j}|nS\rangle_{j+1}$,
in which the first and second processes (I) and (II) induces the two-photon
detuning $\omega_{21}$ and the third process (III) is the exchange
interaction.

Utilizing the gradient magnetic field and a laser-assisted dipole-dipole (DD)
interaction between Rydberg atoms Yang  _et al._ (2022), we can realize the
term $\hat{H}_{g}+\hat{H}_{J}(t)$. For two atoms in different Rydberg states,
there is DD interaction between them, which scales as $1/R^{3}$ with $R$ being
the distance between the two atoms Browaeys  _et al._ (2016); Browaeys and
Lahaye (2020). For sufficiently large lattice constant $a$, one can only
consider the nearest-neighboring DD interaction. The direct DD interaction
between site $j$ and $j+1$ exchanges the Rydberg states in the two sites and
reads $V_{d}[(|nS\rangle\langle n^{\prime}P|)_{j}(|n^{\prime}P\rangle\langle
nS|)_{j+1}+h.c.]$ with $V_{d}=C_{3}/a^{3}$ and $C_{3}\sim n^{4}$ being the DD
interaction coefficient Browaeys  _et al._ (2016). It simulates the tunneling
term in our model. To make the tunneling strength vary in time periodically,
we propose a Raman coupling scheme, which contains two Raman lasers coupling
the state $|nS\rangle$ to a low-lying state $|\alpha\rangle$, as described by

$V_{R}(t)=\sum_{j}\\{[\Omega_{1}\cos(\omega_{\mathrm{R}_{1}}t)+\Omega_{2}\cos(\omega_{\mathrm{R}_{2}}t)](|nS\rangle\langle\alpha|)_{j}+h.c.\\}.$
(28)

Here $\Omega_{1,2}$ are the Rabi frequency and
$\omega_{\mathrm{R}_{1},\mathrm{R}_{2}}$ are the frequency of the two Raman
lasers. We suppose that the tilting strength $\tilde{g}$ is much larger than
$V_{d}$ and the bare exchange process is suppressed due to this energy offset.
In the same spirit of Ref. Yang  _et al._ (2022), we utilize the two-photon
Raman process to compensate part of the energy offset with
$\tilde{g}-(\omega_{\mathrm{R}_{2}}-\omega_{\mathrm{R}_{1}})=g$ and recover
the exchange coupling.

There are two intermediate states
$|n^{\prime}P\rangle_{j}|\alpha\rangle_{j+1}$ and
$|\alpha\rangle_{j}|n^{\prime}P\rangle_{j+1}$ that help the coupling between
$|nS\rangle_{j}|n^{\prime}P\rangle_{j+1}$ and
$|n^{\prime}P\rangle_{j}|nS\rangle_{j+1}$, with the energy compensation
happening at site $j+1$ and $j$ respectively, see Fig. A1. Labeling
$|a\rangle_{j}=|nS\rangle_{j}|n^{\prime}P\rangle_{j+1}$,
$|b\rangle_{j}=|\alpha\rangle_{j}|n^{\prime}P\rangle_{j+1}$,
$|c\rangle_{j}=|n^{\prime}P\rangle_{j}|\alpha\rangle_{j+1}$ and
$|d\rangle_{j}=|n^{\prime}P\rangle_{j}|nS\rangle_{j+1}$, the laser-assisted DD
interaction in the four-level system under rotating wave approximation reads

$\hat{H}^{\mathrm{LDD}}(j,t)=\begin{pmatrix}E_{a}&\frac{\Omega_{1}}{2}e^{-i\omega_{\mathrm{R}_{1}}t}+\frac{\Omega_{2}}{2}e^{-i\omega_{\mathrm{R}_{2}}t}&0&V_{d}\\\
\frac{\Omega_{1}^{*}}{2}e^{i\omega_{\mathrm{R}_{1}}t}+\frac{\Omega_{2}^{*}}{2}e^{i\omega_{\mathrm{R}_{2}}t}&E_{b}&0&0\\\
0&0&E_{c}&\frac{\Omega_{1}^{*}}{2}e^{i\omega_{\mathrm{R}1}t}+\frac{\Omega_{2}^{*}}{2}e^{i\omega_{\mathrm{R}2}t}\\\
V_{d}&0&\frac{\Omega_{1}}{2}e^{-i\omega_{\mathrm{R}_{1}}t}+\frac{\Omega_{2}}{2}e^{-i\omega_{\mathrm{R}_{2}}t}&E_{d}\end{pmatrix},$
(29)

where $E_{d}-E_{a}=\tilde{g}$, $E_{a}-E_{b}=E_{d}-E_{c}=\Delta$ and $\Delta$
is the energy difference between $|nS\rangle$ and $|\alpha\rangle$.
Transforming to the interaction picture, the Hamiltonian reads

$\hat{H}^{\mathrm{LDD}}_{\mathrm{I}}(j,t)=\begin{pmatrix}0&\frac{\Omega_{1}}{2}e^{i\delta_{1}t}+\frac{\Omega_{2}}{2}e^{i\delta_{2}t}&0&V_{d}e^{-i\tilde{g}t}\\\
\frac{\Omega_{1}^{*}}{2}e^{-i\delta_{1}t}+\frac{\Omega_{2}^{*}}{2}e^{-i\delta_{2}t}&0&0&0\\\
0&0&0&\frac{\Omega_{1}^{*}}{2}e^{-i\delta_{1}t}+\frac{\Omega_{2}^{*}}{2}e^{-i\delta_{2}t}\\\
V_{d}e^{i\tilde{g}t}&0&\frac{\Omega_{1}}{2}e^{i\delta_{1}t}+\frac{\Omega_{2}}{2}e^{i\delta_{2}t}&0\end{pmatrix},\\\
$

with $\delta_{1}=\Delta-\omega_{\mathrm{R}_{1}}$ and
$\delta_{2}=\Delta-\omega_{\mathrm{R}_{2}}$ being the single-photon detunings.
In the interaction picture, the evolution operator can be expanded in the
Dyson series as

$\displaystyle\hat{U}_{\mathrm{I}}(j,t)$ $\displaystyle=$ $\displaystyle
1-i\int_{0}^{t}\hat{H}^{\mathrm{LDD}}_{\mathrm{I}}(j,t_{1})dt_{1}+(-i)^{2}\int_{0}^{t}\hat{H}^{\mathrm{LDD}}_{\mathrm{I}}(j,t_{1})\int_{0}^{t_{1}}\hat{H}^{\mathrm{LDD}}_{\mathrm{I}}(j,t_{2})dt_{2}dt_{1}$
(30) $\displaystyle+$
$\displaystyle(-i)^{3}\int_{0}^{t}\hat{H}^{\mathrm{LDD}}_{\mathrm{I}}(j,t_{1})\int_{0}^{t_{1}}\hat{H}^{\mathrm{LDD}}_{\mathrm{I}}(j,t_{2})\int_{0}^{t_{2}}\hat{H}^{\mathrm{LDD}}_{\mathrm{I}}(j,t_{3})dt_{3}dt_{2}dt_{1}+\cdots.$

Suppose
$\tilde{g},|\delta_{1}|,|\delta_{2}|,|\tilde{g}\pm\delta_{1}|,|\tilde{g}\pm\delta_{2}|,|\delta_{1}-\delta_{2}|\gg
V_{d},|\Omega_{1}|,|\Omega_{2}|$ and
$\tilde{g}\sim\omega_{\mathrm{R2}}-\omega_{R1}$. The first- and second-order
processes are non-resonant. In the third-order processes, the couplings
between $\\{|a\rangle_{j},|d\rangle_{j}\\}$ and
$\\{|b\rangle_{j},|c\rangle_{j}\\}$ are non-resonant, and there are near-
resonant couplings between $|a\rangle_{j}(|b\rangle_{j})$ and
$|d\rangle_{j}(|c\rangle_{j})$. Thus, the evolution operator in the
interaction picture is block diagonalized. According to the near-resonant
processes in Fig. A1, we can calculate the element matrix

$U_{\mathrm{I}}^{(ad)}=\ _{j}\langle
d|\hat{U}_{\mathrm{I}}(j,t)|a\rangle_{j}\\\
=J_{\mathrm{eff}}\frac{e^{-i(\omega_{21}-\tilde{g})t}-1}{\omega_{21}-\tilde{g}},$

where $\omega_{21}=\omega_{\mathrm{R}_{2}}-\omega_{\mathrm{R}_{1}}$, and

$J_{\mathrm{eff}}=\frac{V_{d}\Omega_{1}^{*}\Omega_{2}}{4}\left[\frac{1}{\delta_{1}\omega_{21}}+\frac{1}{\tilde{g}(\tilde{g}-\delta_{1})}\right].$
(31)

The evolution operator in the bases $\\{|a\rangle_{j},|d\rangle_{j}\\}$ up to
third-order Dyson expansion then reads

$\hat{U}_{\mathrm{I}}^{(ad)}(j,t)\approx
1+\left[U_{\mathrm{I}}^{(ad)}(|d\rangle\langle a|)_{j}+h.c.\right].$ (32)

Transforming to the Schrödinger picture, we can obtain the Hamiltonian

$\displaystyle\hat{H}^{(ad)}$
$\displaystyle(j,t)=\hat{H}_{0}(j)+ie^{-i\hat{H}_{0}(j)t}\frac{\partial\hat{U}^{(ad)}_{\mathrm{I}}}{\partial
t}{\hat{U}}^{(ad)\dagger}_{\mathrm{I}}e^{i\hat{H}_{0}(j)t}$ $\displaystyle=$
$\displaystyle\hat{H}_{0}(j)+J_{\mathrm{eff}}e^{-i(\omega_{21}-\tilde{g}+E_{d}-E_{a})t}(|d\rangle\langle
a|)_{j}+h.c.,$

where $\hat{H}_{0}(j)=E_{a}(|a\rangle\langle a|)_{j}+E_{d}(|d\rangle\langle
d|)_{j}$. At last, we perform a rotating transformation and obtain the
effective Hamiltonian

$\displaystyle\hat{H}^{(ad)}_{\mathrm{eff}}(j)$ $\displaystyle=$
$\displaystyle\hat{S}(t)\hat{H}^{(ad)}(j,t)\hat{S}^{\dagger}(t)-i\hat{S}(t)\frac{\partial\hat{S}^{\dagger}(t)}{\partial
t}$ (34) $\displaystyle=$
$\displaystyle-\frac{\tilde{g}-\omega_{21}}{2}(|a\rangle\langle
a|)_{j}+\frac{\tilde{g}-\omega_{21}}{2}(|d\rangle\langle d|)_{j}$
$\displaystyle+$ $\displaystyle J_{\mathrm{eff}}(|d\rangle\langle
a|)_{j}+h.c.,$

where

$\hat{S}(t)=\begin{pmatrix}e^{i(E_{a}-\frac{\omega_{21}-\tilde{g}}{2})t}&0\\\
0&e^{i(E_{d}+\frac{\omega_{21}-\tilde{g}}{2})t}.\end{pmatrix}$ (35)

If we set $\Omega_{1}=\Omega_{2}=\Omega$, and tune the Rabi frequency as
$|\Omega|^{2}=|\Omega_{0}|^{2}|[1+u\cos(\omega t)]$, with the cosinoidal
driving been realized experimentally Omran  _et al._ (2019), the sum of
$\hat{H}_{\mathrm{eff}}^{(ad)}(j)$ over $j$ simulates
$\hat{H}_{g}+\hat{H}_{J}(t)$ with
$J=V_{d}|\Omega_{0}|^{2}[\frac{1}{4\delta_{1}\omega_{21}}+\frac{1}{4\tilde{g}(\tilde{g}-\delta_{1})}]$
and $g=\tilde{g}-\omega_{21}$.

The van der Waals interaction between Rydberg atoms can be used to simulate
$\hat{H}_{U}$. The atoms in the same Rydberg states can interact through the
van der Waals interaction, which scales as $1/R^{6}$ Browaeys  _et al._
(2016); Browaeys and Lahaye (2020). If the lattice constant is sufficiently
large, one can only consider the nearest-neighboring interaction. The sum of
the van der Walls interactions of the two kinds of Rydberg states read

$\displaystyle\hat{H}_{\mathrm{vdW}}$ $\displaystyle=$
$\displaystyle\frac{C_{6}^{(2)}+C_{6}^{(1)}}{a^{6}}\sum_{j=0}^{L-2}\left(|n^{\prime}P\rangle\langle
n^{\prime}P|\right)_{j}\left(|n^{\prime}P\rangle\langle
n^{\prime}P|\right)_{j+1}$ (36) $\displaystyle+$
$\displaystyle\frac{C_{6}^{(1)}}{a^{6}}\left[\left(|n^{\prime}P\rangle\langle
n^{\prime}P|\right)_{0}+\left(|n^{\prime}P\rangle\langle
n^{\prime}P|\right)_{L-1}\right],$

where $C_{6}^{(1)}\sim n^{11}$ and $C_{6}^{(2)}\sim n^{\prime 11}$ are the van
der Waals coefficients of $|nS\rangle$ and $|n^{\prime}P\rangle$ Browaeys  _et
al._ (2016); Browaeys and Lahaye (2020), and a constant energy
$-LC_{6}^{(1)}/a^{6}$ has been dropped. The defect on the boundary sites can
be compensated by light shift. Thus, the van der Walls interaction simulates
$\hat{H}_{U}$ with $U=[C_{6}^{(2)}+C_{6}^{(1)}]/a^{6}$.

##  Appendix B The odd-even effect of particle number on the largest Krylov
subspace when $\omega=\omega_{2}$

The spatial-reflection operator reflects the particles on a Fock state
$|n_{0}n_{1}\cdots n_{L-2}n_{L-1}\rangle$, such that
$\hat{R}|\vec{n}\rangle=|n_{L-1}n_{L-2}\cdots
n_{1}n_{0}\rangle\equiv|\underline{\vec{n}}\rangle$. It is easy to see that

$\displaystyle\hat{\mathcal{P}}^{(2g)}_{L-j-1,L-j+2}\hat{c}^{\dagger}_{L-j+1}\hat{c}_{L-j}$
$\displaystyle=$
$\displaystyle\hat{R}\hat{\mathcal{P}}^{(0)}_{j-1,j+2}\hat{c}^{\dagger}_{j}\hat{c}_{j+1}\hat{R}^{-1},$
$\displaystyle\hat{\mathcal{P}}^{(2g)}_{L-j-1,L-j+2}\hat{c}^{\dagger}_{L-j}\hat{c}_{L-j+1}$
$\displaystyle=$
$\displaystyle\hat{R}\hat{\mathcal{P}}^{(0)}_{j-1,j+2}\hat{c}^{\dagger}_{j+1}\hat{c}_{j}\hat{R}^{-1}.$

Thus, corresponding to a nonzero
$\langle\vec{l}|\hat{H}^{\mathrm{eff}}_{\omega_{2}}|\vec{n}\rangle$, the
matrix element
$\langle\vec{\underline{l}}|\hat{H}^{\mathrm{eff}}_{\omega_{2}}|\vec{\underline{n}}\rangle$
is also nonzero. It is obvious that if a Krylov subspace contains a pair of
mutually reflected states $\\{|\vec{n}\rangle,|\underline{\vec{n}}\rangle\\}$
(including the case when $|\vec{n}\rangle$ is a reflection invariant state
with $|\underline{\vec{n}}\rangle=|\vec{n}\rangle$), then all the Fock states
in this Krylov subspace can be paired by the spatial reflection
transformation. Otherwise, the Krylov subspace does not contain any pair of
mutually reflected Fock states. Under the spatial reflection, the Krylov
subspace of the former kind is invariant, in the sense that the set of bases
in the Krylov subspace does not change. The Krylov subspace of the the latter
kind will be transformed to its reflection partner.

Depending on the parity of the total particle number $N$, there will exist
reflection invariant Krylov subspace or not. The expectation values of
$\hat{E}$ in a Fock state $|\vec{n}\rangle$ and its reflected state read
$E_{\vec{n}}=\langle\vec{n}|\hat{E}|\vec{n}\rangle=-\sum_{j}jn_{j}+\sum_{j}n_{j}n_{j+1}$
and $E_{\vec{\underline{n}}}=E_{\vec{n}}-2\sum_{j}jn_{L+1-j}+N(L+1)$,
respectively. If $N$ is odd, the parity of $\hat{E}$ in the two states are
opposite to each other, thus $|\vec{n}\rangle$ and
$|\vec{\underline{n}}\rangle$ belong to different parity symmetry spaces. In
this way, the Krylov subspaces in the even- and odd-parity symmetry spaces are
in one-to-one correspondence through the spatial-reflection transformation,
and there is no reflection invariant Krylov subspace. If $N$ is even, the
parity of $\hat{E}$ in a pair of mutually reflected Fock states are the same,
and there exists reflection invariant Krylov subspaces. One such Krylov
subspace is the one that contains the charge density wave state
$|\mathrm{CDW}1\rangle=|0101\cdots 01\rangle$ and its reflected partner
$|\mathrm{CDW}2\rangle=|1010\cdots 10\rangle$. Viewing four sites as a cell,
both the sequences $0101$ in $|\mathrm{CDW}1\rangle$ and $1010$ in
$|\mathrm{CDW}2\rangle$ can be changed to $0110$ under the action of
$\hat{H}^{\mathrm{eff}}_{\omega_{2}}$. Thus, under $N/2$ times of action of
$\hat{H}^{\mathrm{eff}}_{\omega_{2}}$, both $|\mathrm{CDW}1\rangle$ and
$|\mathrm{CDW2}\rangle$ can be couple to the same state $|01100110\cdots
0110\rangle$, which is reflection invariant. So, $|\mathrm{CDW}1\rangle$ and
$|\mathrm{CDW2}\rangle$ are in the same Krylov subspace and this subspace is
reflection invariant.

Numerically, we find that the largest Krylov subspace at all system sizes
contains the charge density wave states. Thus, when $N$ is odd, there are two
largest Krylov subspace in the half-filling sector, which lie in the even- and
odd-parity symmetry space of $\hat{E}$ respectively. When $N$ is even, there
is only one largest Krylov subspace, which lies in the even-parity symmetry
space and is reflection invariant [note that $E_{\mathrm{CDW}1}=-N^{2}$ and
$E_{\mathrm{CDW}2}=-N(N-1)$]. This odd-even effect leads to different fitting
functions for quadruple and non-quadruple lattice sites.

## References

  * Deutsch (1991) J. M. Deutsch, Quantum statistical mechanics in a closed system, Phys. Rev. A 43, 2046 (1991)
  * Srednicki (1994) M. Srednicki, Chaos and quantum thermalization, Phys. Rev. E 50, 888 (1994)
  * Srednicki (1999) M. Srednicki, The approach to thermal equilibrium in quantized chaotic systems, J. Phys. A 32, 1163 (1999)
  * Rigol and Dunjko (2008) M. Rigol, M. Dunjko and V. Olshanii, Thermalization and its mechanism for generic isolated quantum systems, Nature 452, 854 (2008)
  * Abanin  _et al._ (2016) D. A. Abanin, W. De Roeck, and F. Huveneers, Theory of many-body localization in periodically driven systems, Ann. Phys. 372, 1 (2016)
  * Nandkishore and Huse (2015) R. Nandkishore and D. A. Huse, Many-Body Localization and Thermalization in Quantum Statistical Mechanics, Annu. Rev. Condens. Matter. Phys. 6, 15 (2015)
  * Abanin  _et al._ (2019) D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Colloquium: Many-body localization, thermalization, and entanglement, Rev. Mod. Phys. 91, 021001 (2019)
  * Moudgalya  _et al._ (2022) S. Moudgalya, B. A. Bernevig, and N. Regnault, Quantum many-body scars and Hilbert space fragmentation: a review of exact results, Rep. Prog. Phys. 85, 086501 (2022)
  * Kinoshita  _et al._ (2006) T. Kinoshita, T. Wenger, and D. S. Weiss, A quantum Newton’s cradle, Nature 440, 900 (2006)
  * Rigol  _et al._ (2007) M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, Relaxation in a Completely Integrable Many-Body Quantum System: An Ab Initio Study of the Dynamics of the Highly Excited States of 1D Lattice Hard-Core Bosons, Phys. Rev. Lett. 98, 050405 (2007)
  * Polkovnikov  _et al._ (2011) A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Colloquium: Nonequilibrium dynamics of closed interacting quantum systems, Rev. Mod. Phys. 83, 863 (2011)
  * 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)
  * Schulz  _et al._ (2019) M. Schulz, C. A. Hooley, R. Moessner, and F. Pollmann, Stark Many-Body Localization, Phys. Rev. Lett. 122, 040606 (2019)
  * van Nieuwenburg  _et al._ (2019) E. van Nieuwenburg, Y. Baum, and G. Refael, From Bloch oscillations to many-body localization in clean interacting systems, Proc. Natl. Acad. Sci. U.S.A. 116, 9269 (2019)
  * Taylor  _et al._ (2020) S. R. Taylor, M. Schulz, F. Pollmann, and R. Moessner, Experimental probes of Stark many-body localization, Phys. Rev. B 102, 054206 (2020)
  * Yao and Zakrzewski (2020) R. Yao and J. Zakrzewski, Many-body localization of bosons in an optical lattice: Dynamics in disorder-free potentials, Phys. Rev. B 102, 104203 (2020)
  * Chanda  _et al._ (2020) T. Chanda, R. Yao, and J. Zakrzewski, Coexistence of localized and extended phases: Many-body localization in a harmonic trap, Phys. Rev. Research 2, 032039(R) (2020)
  * Zhang  _et al._ (2021) L. Zhang, Y. Ke, W. Liu, and C. Lee, Mobility edge of Stark many-body localization, Phys. Rev. A 103, 023323 (2021)
  * Yao  _et al._ (2021) R. Yao, T. Chanda, and J. Zakrzewski, Many-body localization in tilted and harmonic potentials, Phys. Rev. B 104, 014201 (2021)
  * Wang  _et al._ (2021) Y.-Y. Wang, Z.-H. Sun, and H. Fan, Stark many-body localization transitions in superconducting circuits, Phys. Rev. B 104, 205122 (2021)
  * Guo  _et al._ (2021) Q. Guo, C. Cheng, H. Li, S. Xu, P. Zhang, Z. Wang, C. Song, W. Liu, W. Ren, H. Dong, R. Mondaini, and H. Wang, Stark Many-Body Localization on a Superconducting Quantum Processor, Phys. Rev. Lett. 127, 240502 (2021)
  * Morong  _et al._ (2021) W. Morong, F. Liu, P. Becher, K. S. Collins, L. Feng, A. Kyprianidis, G. Pagano, T. You, A. V. Gorshkov, and C. Monroe, Observation of Stark many-body localization without disorder, Nature 599, 393 (2021)
  * 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, Nat. Phys. 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)
  * 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)
  * Serbyn  _et al._ (2021) M. Serbyn, D. A. Abanin, and Z. Papić, Quantum many-body scars and weak breaking of ergodicity, Nat. Phys. 17, 675 (2021)
  * Chandran  _et al._ (2023) A. Chandran, T. Iadecola, V. Khemani, and R. Moessner, Quantum Many-Body Scars: A Quasiparticle Perspective, Annu. Rev. Condens. Matter. Phys. 14, 443 (2023)
  * Pai  _et al._ (2019) S. Pai, M. Pretko, and Rahul M. Nandkishore, Localization in Fractonic Random Circuits, Phys. Rev. X 9, 021003 (2019)
  * 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)
  * 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,  _Memorial Volume for Shoucheng Zhang_ (2021) Chap. Chapter 7, pp. 147
  * De Tomasi  _et al._ (2019) G. De Tomasi, D. Hetterich, P. Sala, and F. Pollmann, Dynamics of strongly interacting systems: From Fock-space fragmentation to many-body localization, Phys. Rev. B 100, 214313 (2019)
  * Yang  _et al._ (2020) Zhi-Cheng Yang, Fangli Liu, Alexey V. Gorshkov, and Thomas Iadecola, Hilbert-space fragmentation from strict confinement, Phys. Rev. Lett. 124, 207602 (2020)
  * Frey  _et al._ (2022) P. Frey, L. Hackl, and S. Rachel, Hilbert space fragmentation and interaction-induced localization in the extended Fermi-Hubbard model, Phys. Rev. B 106, L220301 (2022)
  * Ghosh  _et al._ (2023) Somsubhra Ghosh, Indranil Paul, and K. Sengupta, Prethermal fragmentation in a periodically driven fermionic chain, Phys. Rev. Lett. 130, 120401 (2023)
  * Scherg  _et al._ (2021) S. Scherg, T. Kohlert, P. Sala, F. Pollmann, B. Hebbe Madhusudhana, I. Bloch, and M. Aidelsburger, Observing non-ergodicity due to kinetic constraints in tilted Fermi-Hubbard chains, Nat. Commun. 12, 4490 (2021)
  * Kohlert  _et al._ (2023) T. Kohlert, S. Scherg, P. Sala, F. Pollmann, B. Hebbe Madhusudhana, I. Bloch, and M. Aidelsburger, Exploring the Regime of Fragmentation in Strongly Tilted Fermi-Hubbard Chains, Phys. Rev. Lett. 130, 010201 (2023)
  * Morningstar  _et al._ (2020) A. Morningstar, V. Khemani, and D. A. Huse, Kinetically constrained freezing transition in a dipole-conserving system, Phys. Rev. B 101, 214205 (2020)
  * Pozderac  _et al._ (2023) C. Pozderac, S. Speck, X. Feng, D. A. Huse, and B. Skinner, Exact solution for the filling-induced thermalization transition in a one-dimensional fracton system, Phys. Rev. B 107, 045137 (2023)
  * Eckardt (2017) A. Eckardt, Colloquium: Atomic quantum gases in periodically driven optical lattices, Rev. Mod. Phys. 89, 011004 (2017)
  * Greschner  _et al._ (2014) S. Greschner, L. Santos, and D. Poletti, Exploring Unconventional Hubbard Models with Doubly Modulated Lattice Gases, Phys. Rev. Lett. 113, 183002 (2014)
  * Meinert  _et al._ (2016) F. Meinert, M. J. Mark, K. Lauber, A. J. Daley, and H.-C. Nägerl, Floquet Engineering of Correlated Tunneling in the Bose-Hubbard Model with Ultracold Atoms, Phys. Rev. Lett. 116, 205301 (2016)
  * Zhao  _et al._ (2019) H. Zhao, J. Knolle, and F. Mintert, Engineered nearest-neighbor interactions with doubly modulated optical lattices, Phys. Rev. A 100, 053610 (2019)
  * Zhao  _et al._ (2020) H. Zhao, J. Vovrosh, F. Mintert, and J. Knolle, Quantum many-body scars in optical lattices, Phys. Rev. Lett. 124, 160604 (2020)
  * Soori and Sen (2010) A. Soori and D. Sen, Nonadiabatic charge pumping by oscillating potentials in one dimension: Results for infinite system and finite ring, Phys. Rev. B 82, 115432 (2010)
  * Sen  _et al._ (2021) A. Sen, D. Sen, and K. Sengupta, Analytic approaches to periodically driven closed quantum systems: methods and applications, J. Phys.: Condens.Matter 33, 443003 (2021)
  * Page (1993) Don N. Page, Average entropy of a subsystem, Phys. Rev. Lett. 71, 1291 (1993)
  * 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)
  * Bravyi  _et al._ (2011) S. Bravyi, D. P. DiVincenzo, and D. Loss, Schrieffer-Wolff transformation for quantum many-body systems, Ann. Phys. 326, 2793 (2011)
  * Guo  _et al._ (2023) Z. Guo, B. Liu, Y. Gao, A. Yang, J. Wang, J. Ma, and L. Ying, Origin of Hilbert-space quantum scars in unconstrained models, Phys. Rev. B 108, 075124 (2023)
  * Barredo  _et al._ (2016) D. Barredo, S. de Léséleuc, V. Lienhard, T. Lahaye, and A. Browaeys, An atom-by-atom assembler of defect-free arbitrary two-dimensional atomic arrays, Science 354, 1021 (2016)
  * Endres  _et al._ (2016) M. Endres, H. Bernien, A. Keesling, H. Levine, E. R. Anschuetz, A. Krajenbrink, C. Senko, V. Vuletic, M. Greiner, and M. D. Lukin, Atom-by-atom assembly of defect-free one-dimensional cold atom arrays, Science 354, 1024 (2016)
  * Browaeys  _et al._ (2016) A. Browaeys, D. Barredo, and T. Lahaye, Experimental investigations of dipole-dipole interactions between a few rydberg atoms, J. Phys. B: At. Mol. Opt. Phys. 49, 152001 (2016)
  * Barredo  _et al._ (2018) D. Barredo, V. Lienhard, S. de Léséleuc, T. Lahaye, and A. Browaeys, Synthetic three-dimensional atomic structures assembled atom by atom, Nature (London) 561, 79 (2018)
  * Browaeys and Lahaye (2020) A. Browaeys and T. Lahaye, Many-body physics with individually controlled rydberg atoms, Nat. Phys. 16, 132 (2020)
  * Yang  _et al._ (2022) T.-H. Yang, B.-Z. Wang, X.-C. Zhou, and X.-J. Liu, Quantum Hall states for Rydberg arrays with laser-assisted dipole-dipole interactions, Phys. Rev. A 106, L021101 (2022)
  * Omran  _et al._ (2019) A. Omran, H. Levine, A. Keesling, G. Semeghini, T. T. Wang, S. Ebadi, H. Bernien, A. S. Zibrov, H. Pichler, S. Choi, J. Cui, M. Rossignolo, P. Rembold, S. Montangero, T. Calarco, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Generation and manipulation of Schrödinger cat states in Rydberg atom arrays, Science 365, 570 (2019)