# Discrete time crystals enforced by Floquet-Bloch scars

Biao Huang<EMAIL_ADDRESS>Kavli Institute for Theoretical
Sciences, University of Chinese Academy of Sciences, Beijing 100190, China
Tsz-Him Leung  Department of Physics, University of California, Berkeley, USA
Dan Stamper-Kurn  Department of Physics, University of California, Berkeley,
USA  Materials Sciences Division, Lawrence Berkeley National Laboratory, USA
W. Vincent Liu<EMAIL_ADDRESS>Department of Physics and Astronomy, University
of Pittsburgh, Pittsburgh, PA 15260, USA  Shenzhen Institute for Quantum
Science and Engineering, Southern University of Science and Technology,
Shenzhen 518055, China

###### Abstract

We analytically identify a new class of quantum scars protected by
spatiotemporal translation symmetries, dubbed Floquet-Bloch scars. They
distinguish from previous (quasi-)static scars by a rigid spectral pairing
only possible in Floquet systems, where strong interaction and drivings
equalize the quasienergy corrections to all scars and maintain their spectral
spacings against generic bilinear perturbations. Scars then enforce the
spatial localization and rigid discrete time crystal (DTC) oscillations as
verified numerically in a trimerized kagome lattice model relevant to recent
cold atom experiments. Our analytical solutions offer a potential scheme to
understand the mechanisms for more generic translation-invariant DTCs.

Introduction — Systems far from equilibrium have become a fertile ground
cultivating unexpected phenomena recently. Among them, discrete time crystals
(DTC) [1, 2, 3, 4, 5, 6, 7] constitute an intriguing example. As foundational
concepts of ground state and temperature fall apart in the absence of thermal
equilibrium, Landau’s theory of symmetry breaking [8] is replaced by new
principles like spectral pairing and eigenstate orders [1, 2] in handling time
translation symmetries. That results in the DTC phenomena where Hamiltonians
$\hat{H}(t+T)=\hat{H}(t)$ give rise to observables $O(t+NT)=O(t)$
($1<N\in\mathbb{Z}$) oscillating like a temporal charge/spin density wave.
Crucially, the periodicity $NT$ demands no fine-tuning and withstands generic
perturbations.

The concept of DTCs has been considered in several physical realizations [9,
10, 11, 12, 13, 14, 15]. While the strongly disordered cases are relatively
well understood [1, 2, 3, 16], the possibility of DTCs in translation-
invariant ordered systems is less clear. Empirical evidence for DTCs is
accumulating in both physical and numerical experiments [9, 10, 17, 18, 19,
20, 21, 22, 23]. However, analytical explanations based on many-body
localization (MBL) [24] or prethermalization [25] do not seem to apply to
these cases. Recently, it was indicated that quasi-conservation laws [26, 27],
which can be enhanced by single-particle terms, may help protect phenomena
pertinent to DTCs. Meanwhile, the initial state dependence of clean DTCs [28,
26] has been reexamined in terms of scar physics [29] in recent numerics [30,
31]. Altogether, continued investigation on DTCs in non-disordered systems,
with the objectives of uncovering the underlying mechanism that supports the
DTC and the specific role of many-body (vs. single-particle) effects, is
warranted.

In this Letter, we gain insights on these two research objectives by studying
a small cluster of soft-core bosons on driven trimerized kagome lattices,
relevant to recent experiments [32] and feasible for numerical verifications.
We find analytically that it is a class of quantum scars protected by
spatiotemporal translation invariance, dubbed Floquet-Bloch scars (FBS), that
gives rise to DTC behaviors for sublattice density oscillations. FBS’s
identified here neither exploit a static scar (i.e. “PXP” model [29, 33, 34])
nor end up with engineered static Hamiltonians. Instead, these FBS’s exhibit a
unique DTC feature. Specifically, each scar quasienergy may be shifted
considerably under perturbation. However, the interplay of strong interactions
and drivings equalizes the scar level shifts, which is proved to all orders in
our perturbative treatment. Then, the quasienergy difference $\omega_{0}$
between FBSs remains invariant and enforces the persisting
$2\pi/\omega_{0}$-periodic DTC. Rigid scar level spacing here resembles the
“spectral pairing rigidity” for all Floquet eigenstates in MBL DTCs [1, 16].
Also, such a mechanism allows for rather generic perturbations compared with
preexisting scar models typically relying on microscopic details to achieve
configuration separations [35, 36, 37, 38, 39, 34, 40, 41, 29, 33, 42, 43,
44]. Thus, our analytical solutions not only offer a more definitive
understanding of clean DTC mechanisms, but also point out a new way of
constructing scars showing peculiar spectral orders characteristic of Floquet
systems.

Model and phenomena — We consider bosons evolving under a Hamiltonian
$\hat{H}(t+T)=\hat{H}(t)$ that is toggled between two settings repetitively
within each period $T$:

$\displaystyle\hat{H}_{1}T/2\hbar$
$\displaystyle=\phi_{1}\sum_{\boldsymbol{r},\mu\neq\nu}if_{\mu\nu}\left[\hat{\psi}_{\boldsymbol{r}\mu}^{\dagger}\hat{\psi}_{\boldsymbol{r}\nu}+\lambda\hat{\psi}_{\boldsymbol{r}+\boldsymbol{e}_{\mu},\mu}^{\dagger}\hat{\psi}_{\boldsymbol{r}+\boldsymbol{e}_{\nu},\nu}\right],\,t\in[0,T/2)$
$\displaystyle\hat{H}_{2}T/2\hbar$
$\displaystyle=\sum_{\boldsymbol{r},\mu}\left[\phi_{2}\hat{n}_{\boldsymbol{r}\mu}(\hat{n}_{\boldsymbol{r}\mu}-1)+\theta_{\mu}\hat{n}_{\boldsymbol{r}\mu}\right],\,t\in[T/2,T).$
(1)

Here, $\hat{H}_{1}$ describes the hopping of non-interacting bosons in a
trimerized kagome lattice with complex hopping amplitudes, as shown in Fig.
1(a), while $\hat{H}_{2}$ describes the combination of on-site single-particle
and interaction energy shifts. Dimensionless parameters
$(\phi_{1},\lambda,\phi_{2},\theta_{\mu})$ characterize the Floquet operator
$\hat{U}_{F}=P_{t}e^{-(i/\hbar)\int_{0}^{T}dt\hat{H}(t)}=e^{-i\hat{H}_{2}T/2\hbar}e^{-i\hat{H}_{1}T/2\hbar}$.
$\hat{\psi}_{\boldsymbol{r}\mu}$ and
$\hat{n}_{\boldsymbol{r}\mu}=\hat{\psi}_{\boldsymbol{r}\mu}^{\dagger}\hat{\psi}_{\boldsymbol{r}\mu}$
are annihilation and particle number operators respectively, for $L^{2}$ unit
cells $\boldsymbol{r}=m_{1}\boldsymbol{e}_{1}+m_{2}\boldsymbol{e}_{2}$
($m_{1,2}=0,1,\dots,L-1$) and three sublattices $\mu,\nu=0,1,2$. Here
$\boldsymbol{e}_{1,2}$ are Bravais vectors for kagome lattices and
$\boldsymbol{e}_{0}\equiv\boldsymbol{0}$. $if_{\mu\nu}=(1+2e^{2\pi
i(\mu-\nu)/3})/\sqrt{3}=\pm i$ specifies the $+i$ hopping directions in Fig. 1
(a). Note that $\sum_{\mu}\theta_{\mu}=0$ can always be achieved by
subtracting $(N_{b}/3)\sum_{\mu}\theta_{\mu}$ from $\hat{H}_{2}T/2\hbar$,
where total bosons
$N_{b}=\sum_{\boldsymbol{r}\mu}\hat{n}_{\boldsymbol{r}\mu}$.

(a) Lattice and phenomena

(c) Eigenstate correlations

(b) Dynamics ($m\in\mathbb{Z}$)

Figure 1:  (a) Trimerized kagome lattice and the schematic illustration of DTC
dynamics. Triangles with strong/weak bonds are denoted by black/gray colors,
with $+i$ hopping directions indicated by arrows. (b) Particle dynamics
$n_{\mu}(NT)=\sum_{\boldsymbol{r}}\langle\psi_{\text{ini}}|(\hat{U}_{F}^{\dagger})^{N}n_{\boldsymbol{r}\mu}\hat{U}_{F}^{N}|\psi_{\text{ini}}\rangle$,
with the initial state $|\psi_{\text{ini}}\rangle$ at $t/T=1$ that all
particles locate at a single site $\boldsymbol{r}=\boldsymbol{0},\mu=0$. To
facilitate reading, data is grouped into 3 sets at $t\mod 3T=0,1,2$
respectively. For comparison, the non-interacting $\phi_{2}=0,L=3$ case is
plotted at all $t$ in the upper panel as translucent dots. (c) Temporal
correlation functions indicating infinite-time response frequencies ($L=3$).
Unless denoted otherwise, $N_{b}=5$, $\phi_{1}=2\pi/3\sqrt{3}$, $\lambda=0.1$,
$\phi_{2}=1.1$, $\theta_{1,2,3}=(0.1,0.2,-0.3)$.

DTC dynamics obtained by exact diagonalization is briefly shown in Fig. 1.
When $\lambda\rightarrow 0$, $\hat{H}_{1}$ enters the strongly trimerized
regime composed of disconnected triangles, where $\pi/2$-fluxes equalize the
spacing between single-particle flat bands $\omega_{n}=0,\pm\sqrt{3}\phi_{1}$
($\hat{U}_{F}|\omega_{n}\rangle=e^{i\omega_{n}}|\omega_{n}\rangle$). Then,
$\phi_{1}=2\pi/3\sqrt{3}$ leads to $3T$ ballistic oscillations for particles
$\hat{U}_{F}^{\dagger}\hat{\psi}^{\dagger}_{\boldsymbol{r},\mu=0,1,2}\hat{U}_{F}=\hat{\psi}_{\boldsymbol{r},\mu=1,2,0}$
breaking the Hamiltonian time translation symmetry of $T$, as in Fig. 1 (a).
Frequencies given by single-particle physics are, of course, unstable against
perturbations. It is then the hallmark for DTC where strong interactions
$\phi_{2}$ stabilize the $3T$ periodicity without fine-tuning, see Fig. 1 (b).
Late time dynamics can be further confirmed by the temporal correlation
functions $C(\omega)=\sum_{N=-\infty}^{\infty}\frac{e^{i\omega
N}}{2\pi}\sum_{n}\langle\omega_{n}|\hat{P}(N)\hat{P}(0)|\omega_{n}\rangle=\sum_{mn}\delta(\omega-\omega_{mn})A(\omega_{mn})$
for the sublattice density bias, i.e.
$\hat{P}(N)=\left(\hat{U}_{F}^{\dagger}\right)^{N}N_{b}^{-1}\sum_{\boldsymbol{r}}(\hat{n}_{\boldsymbol{r}0}-\hat{n}_{\boldsymbol{r}1})\hat{U}_{F}^{N}$.
Note that the summation $N$ is over infinite time without truncation. The
spectral weight
$A(\omega_{mn})=|\langle\omega_{m}|\hat{P}|\omega_{n}\rangle|^{2}$,
$\omega_{mn}=\omega_{m}-\omega_{n}$ in Fig. 1 (c) showing strong peaks at
frequencies $\omega_{0}\rightarrow\pm 2\pi/3+O(1/D)$ verifies long-time
oscillation periods $2\pi T/|\omega_{0}|\rightarrow 3T$. The small deviation
$O(1/D)$ suppressed by Hilbert space dimension $D$ gives an envelop modulation
in Fig. 1 (b) as noticed previously for both MBL [2, 16] and clean [17] DTCs.

The above phenomena may be viewed from several angles. Particularly, in the
case of complete trimerization, $\lambda=0$, the two-dimensional lattice
breaks up into isolated trimers. DTCs observed in this case is then explained
simply as that of a microscopic three-site chiral system similar to Ref. [45].
If we were to regard intertrimer coupling as simply opening up each one-trimer
DTC to an external bath composed of other trimers, we might expect the overall
DTC dynamics to be destroyed over short time at $\lambda\neq 0$ [46]. Yet,
such expectations contradict results in Fig. 1 (b) (c). Below, we offer an
explanation that DTCs in the coupled-trimer regime is stabilized by a special
class of scar Floquet eigenstates each spanning over the entire two-
dimensional lattice.

(a) Level statistics

(b) EE cut

(c) $\phi_{2}=0.05$

(d) $\phi_{2}=1.1$

(e) Lowest $S_{\text{ent}}$ for $\phi_{2}=1.1$

Figure 2:  (a) $\langle r\rangle$ shows generic ergodicity. (b) Subsystem for
computing EE. (c) – (d) Eigenstate EE in (c) proximate-integrable and (d) DTC
regimes, where low-entropy scars in (d) are highlighted by larger dots. (e)
Lowest EE approaching size-insensitive values at small $\lambda$ and volume
law $S_{\text{ent}}\sim(N_{\text{sub}}/3L^{2})\ln(D_{L})$ at large $\lambda$.
Here $D_{L}$ is the total Hilbert space dimension and $N_{\text{sub}}$ the
subsystem site number enclosed in (b). Inset shows $\Delta
S_{\text{ent}}=S_{\text{ent}}^{(L=3)}-S_{\text{ent}}^{(L=2)}$ near the
crossing $\lambda_{0}\approx 0.135$. Unless specified otherwise, in all plots
parameters are the same as in Fig. 1. Blue (or red) colors denote $L=2$ (or
$L=3$) respectively.

Identifying scars — Quantum scars are rare non-ergodic eigenstates within an
eigenstructure that is otherwise thermalizing [47]. Numerical calculations
confirm the overall thermalizing, non-integrable nature of our model system.
Specifically, we point to two signatures of non-integrability: level-spacing
statistics and entanglement entropy.

Consider first the level spacing. Ordering quasi-energies as
$\omega_{n}<\omega_{n+1}$, following Ref. [48], we test for ergodicity by
calculating the level spacing ratios
$r_{n}=\min(\delta_{n},\delta_{n+1})/\max(\delta_{n},\delta_{n+1})$ for
consecutive gaps $\delta_{n}=\omega_{n+1}-\omega_{n}$. Clearly from Fig. 2
(a), except for a vanishingly small region in proximity to single particle
limit $\phi_{2}\rightarrow 0$, our model is generically far from the
integrable Poissonian case $\langle r\rangle\rightarrow 0.39$. We also note a
crossover between two ergodic Gaussian orthogonal/unitary ensembles (GOE/GUE)
purely by different drivings, an interesting feature previous seen in spin
models [49].

We next exploit the entanglement entropy (EE) to examine each Floquet
eigenstate $|\omega_{n}\rangle$. Reduced density matrices
$\rho_{A}=\text{Tr}_{B}(|\omega_{n}\rangle\langle\omega_{n}|)$ for subsystem
$A$ (region enclosed by highlighted paths in Fig. 2 (b)) can be formed by
tracing out the remaining part $B$ in real space. The EE
$S_{\text{ent}}=-\text{Tr}\left(\rho_{A}\ln\rho_{A}\right)$ then shows that in
both proximate-integrable (Fig. 2 (c)) and DTC (Fig. 2 (d)) regimes, majority
eigenstates do exhibit the typical arch shape for $S_{\text{ent}}$ whose
values increase with Hilbert space dimensions [50]. The narrow distribution of
EE for eigenstates of similar quasi-energy in the DTC regime confirms that
majority arch eigenstates are ergodic [50], in consistent with $\langle
r\rangle$ results previously.

However, in the DTC regime, additional non-ergodic states are observed. As
exhibited in Fig. 2(d), we identify precisely $3L^{2}$ low $S_{\text{ent}}$
scar states (each scar dot in the figure is $L^{2}$-fold degenerate). Each set
of scars separates from the others by quasienergy $|\Delta E|\rightarrow
2\pi/3$, corresponding to exactly the DTC frequency in Fig. 1 (c). The scaling
of lowest EE in Fig. 2 (e) shows a system size $L$ insensitive scar EE for
$\lambda\rightarrow 0$. With increasing $\lambda$, a possible transition is
observed around $\lambda_{0}\approx 0.135$ 111Due to limited sizes accessible
here, we would postpone a more comprehensive examination of criticality to
future work and only take $\lambda_{0}$ as a reference scar vanishing point.,
after which all eigenstates approach the volume law ergodic limit.

We have confirmed numerically that parameters in Fig. 2 (c) give rather short
DTC lifetime, unlike the lifetime shown in Fig. 1 (b) for parameters in Fig. 2
(d). It strongly indicates that the DTC behaviors here are intimately
associated with scars rather than (approximate) overall integrability.

Analytical results for FBS — To characterize these quantum scars further, we
work in the many-body momentum basis [52]
$|\boldsymbol{k},\\{n_{\boldsymbol{r},\mu}\\}\rangle=(1/L)\sum_{m_{1},m_{2}=0}^{L-1}e^{(2\pi
i/L)(k_{1}m_{1}+k_{2}m_{2})}\left|\\{n_{\boldsymbol{r}+m_{1}\boldsymbol{e}_{1}+m_{2}\boldsymbol{e}_{2},\mu}\\}\right\rangle$
constructed from Fock basis
$|\\{n_{\boldsymbol{r}\mu}\\}\rangle=\prod_{\boldsymbol{r}\mu}\left[\left(\hat{\psi}_{\boldsymbol{r}\mu}^{\dagger}\right)^{n_{\boldsymbol{r}\mu}}/\sqrt{n_{\boldsymbol{r}\mu}!}\right]|0\rangle$.
Here $\\{n_{\boldsymbol{r}\mu}\\}$ specifies occupation numbers at different
sites, and $\boldsymbol{k}\sim k_{1,2}=0,1,\ldots,L-1$. Then, translation-
invariant $\hat{U}_{F}$ are block-diagonalized
$\langle\boldsymbol{k},\\{n_{\boldsymbol{r}\mu}\\}|U_{F}|\boldsymbol{k}^{\prime},\\{n_{\boldsymbol{r}\mu}^{\prime}\\}\rangle\sim\delta_{\boldsymbol{k},\boldsymbol{k}^{\prime}}$.
Each $\boldsymbol{k}$ sector would be shown later to host $3$ scar states,
leading to the $3\times L^{2}$-fold scars in Fig. 2 (d).

It is helpful to write down the solution
$U_{F}|\boldsymbol{k},\ell,\\{n_{\boldsymbol{r}\mu}\\}\rangle=e^{iE(\ell,\\{n_{\boldsymbol{r}\mu}\\})}|\boldsymbol{k},\ell,\\{n_{\boldsymbol{r}}\\}\rangle$
to Eqs. (1) at the anchor point $\lambda=0$,

$\displaystyle\left|\boldsymbol{k},\ell,\\{n_{\boldsymbol{r}\mu}\\}\right\rangle$
$\displaystyle=\,\frac{1}{\sqrt{3}}\sum_{m=0,1,2}e^{-i(\frac{2\pi
m}{3}\ell-\alpha_{m})}|\boldsymbol{k},\\{n_{\boldsymbol{r},\mu+m\text{ mod
}3}\\}\rangle,$ (2) $\displaystyle
E\left(\ell,\\{n_{\boldsymbol{r}\mu}\\}\right)$
$\displaystyle=\frac{2\pi}{3}\ell+\phi_{2}\sum_{\boldsymbol{r}\mu}n_{\boldsymbol{r}\mu}(n_{\boldsymbol{r}\mu}-1),\qquad\ell=0,\pm
1,$ (3)

where $\alpha_{0}=0$,
$\alpha_{1}=\sum_{\boldsymbol{r}\mu}\theta_{\mu}n_{\boldsymbol{r}\mu}$, and
$\alpha_{2}=-\sum_{\boldsymbol{r}\mu}\theta_{\mu}n_{\boldsymbol{r},\mu+2\text{
mod }3}$. Each eigenstate populates 3 sublattices $\mu=0,1,2$ coherently, and
therefore an arbitrary Fock state $|\\{n_{\boldsymbol{r}\mu}\\}\rangle$,
usually taken as initial states, will simultaneously overlap with all three
branches $\ell=0,\pm 1$ separating from each other by quasi-energy $|\Delta
E|=2\pi/3$. Then, observables diagonal in the Fock basis, such as
$\hat{O}=\hat{n}_{\boldsymbol{r}\mu}$ or
$\hat{O}=\hat{P}=N_{b}^{-1}\sum_{\boldsymbol{r}}(\hat{n}_{\boldsymbol{r}0}-\hat{n}_{\boldsymbol{r}1})$,
will demonstrate an oscillation $\langle\hat{O}\rangle(t)\sim
c_{1}^{*}c_{2}\langle\boldsymbol{k}_{1},\ell_{1},\\{n_{\boldsymbol{r}\mu}\\}|\hat{O}|\boldsymbol{k}_{2},\ell_{2},\\{n_{\boldsymbol{r}\mu}\\}\rangle
e^{-i\Delta Et/T}+c.c.$ with periodicity $2\pi T/\Delta E=3T$.

Spectral pairing $\Delta E$ for majority eigenstates in Eqs. (2) (3) is, as
expected, unstable against perturbations. The crucial difference here from the
disordered case [1, 16, 2, 3] is the uniform interaction strength $\phi_{2}$
in Eq. (3), which results in an enormous Floquet emergent degeneracy.
Specifically, consider the combination
$Q_{a}=\\{(q^{(a)}_{j},N^{(a)}_{j})|j=1,2,\ldots,M\\}$ for, i.e. $q^{(a)}_{j}$
copies of sites each hosting $N^{(a)}_{j}\geq 0$ particles. In terms of the
Hubbard interaction
$\phi_{2}\sum_{\boldsymbol{r}\mu}n_{\boldsymbol{r}\mu}(n_{\boldsymbol{r}\mu}-1)=\phi_{2}\sum_{j}q^{(a)}_{j}N^{(a)}_{j}(N^{(a)}_{j}-1)$,
each $Q_{a}$ manifold contains degenerate levels of different
$\\{n_{\boldsymbol{r}\mu}\\}$ as
$\text{deg}(Q_{a})=(3L^{2})!/\prod_{j=1}^{M}q^{(a)}_{j}!$. The degeneracy,
though partially lifted by $2\pi\ell/3$ in Eq. (3), leads to the instability
that a small perturbation could generally trigger a reconstruction for
extensive numbers of eigenstates in Eq. (2) with different configurations
$\\{n_{\boldsymbol{r}\mu}\\}$, leading to the ergodicity as indicated by Fig.
2. Correspondingly, a Fock initial state would overlap with large numbers of
eigenstates with different quasienergies without rigid spectral pairings.

To identify FBS, we then seek for manifolds with low degeneracy. Except for a
homogeneous distribution $n_{\boldsymbol{r}\mu}=N_{b}/3L^{2}$ (deg=$1$)
without dynamical signatures, the lowest degenerate
$\\{n_{\boldsymbol{r}\mu}\\}$ deposit all $N_{b}$ bosons into a single site
$n_{\boldsymbol{r}\mu}=\delta_{\boldsymbol{r},\boldsymbol{r}_{0}}\delta_{\mu,\mu_{0}}N_{b}$.
There are apparently $3L^{2}$ such $\\{n_{\boldsymbol{r}\mu}\\}$ with $N_{b}$
bosons allocated into different sites $(\boldsymbol{r}_{0},\mu_{0})$. They
compose the FBS eigenstates

$\displaystyle|\boldsymbol{k},\ell,N_{b}\rangle=\frac{1}{\sqrt{3}}\sum_{m=0,1,2}e^{-i(\frac{2\pi
m}{3}\ell-\alpha_{m})}|\boldsymbol{k},\\{n_{\boldsymbol{r}\mu}=\delta_{\boldsymbol{r},\boldsymbol{0}}\delta_{\mu,m}N_{b}\\}\rangle,$
(4)

with quasienergy $E_{\text{scar}}(\ell)=2\pi\ell/3+\phi_{2}N_{b}(N_{b}-1)$.
The $3L^{2}$ FBSs equally partition into $L^{2}$ conserved many-body momentum
$\boldsymbol{k}$ sectors, each hosting 3 scars with $\ell=0,\pm 1$. Spatial
translation symmetry then forbids hybridizing eigenstates of different
$\boldsymbol{k}$, and temporal translation symmetry protects the conserved
quasienergy separating different $\ell$ by $|\Delta
E|=|E_{\text{scar}}(\ell+1)-E_{\text{scar}}(\ell)|=2\pi/3$. Therefore, FBS’s
experience no degenerate-level perturbations.

It still remains to consider non-degenerate perturbations. In particular, the
periodicity $2\pi$ of Floquet quasienergy constrains Hubbard-interaction gap
for different $Q_{a}$ to be of the order unity. Then, one may expect each scar
level to receive an energy correction $\sim\lambda^{2}$ (of Fermi golden rule
type), resulting in fast detuning within $t\sim T/\lambda^{2}\sim 100T$ for
$\lambda=0.1$. However, such estimations directly contradict Fig. 1 (b).

The resolution turns out to be that all three $\ell=0,\pm 1$ scars are shifted
identically, such that their quasienergy difference, dubbed spectral pairing
gap $|\Delta E|=2\pi/3$ 222In strongly disordered cases [1, 16, 2], spectral
pairing happens for all eigenstates due to many-body localization. Here we use
the terminology to describe similar behaviors for scars but due to different
reasons., is unchanged. In Supplemental Materials (SM) [54], we construct the
strong-drive perturbation theory. For conciseness, we illuminate the essential
physics below by elaborating results up to the second order in the
perturbation series, while higher orders cases are left to SM [54].

Arrange a Floquet operator in the form $U_{F}=U_{0}U^{\prime}$, where $U_{0}$
corresponds to Eq. (1) at $\lambda=0$, and perturbations are factored into
$U^{\prime}\equiv e^{i\lambda H^{\prime}}$. For our purposes, it is more than
enough to take $H^{\prime}$ as a generic hopping Hamiltonian $H^{\prime}=\sum
J_{\boldsymbol{r}\mu\neq\boldsymbol{r}^{\prime}\mu^{\prime}}\hat{\psi}_{\boldsymbol{r}\mu}^{\dagger}\hat{\psi}_{\boldsymbol{r}^{\prime}\mu^{\prime}}$.
(See SM [54] for factorization process). Scar quasienergy corrections
$e^{i\tilde{E}_{\text{scar}}(\ell)}=e^{i(E_{\text{scar}}(\ell)+\sum_{\alpha=1}^{\infty}\lambda^{\alpha}E_{\ell}^{(\alpha)})}$
up to the second order read
$E^{(1)}_{\ell}=\langle\boldsymbol{k},\ell,N_{b}|H^{\prime}|\boldsymbol{k},\ell,N_{b}\rangle$,
$E^{(2)}_{\ell}=-\frac{1}{2}\sum_{(\ell^{\prime},\\{n_{\boldsymbol{r}\mu}\\})}^{\prime}|\langle\boldsymbol{k},\ell^{\prime},\\{n_{\boldsymbol{r}\mu}\\}|H^{\prime}|\boldsymbol{k},\ell,N_{b}\rangle|^{2}\cot\frac{E_{\text{scar}}(\ell)-E(\ell^{\prime},\\{n_{\boldsymbol{r}\mu}\\})}{2}$,
where summation $\sum^{\prime}$ excludes the scar eigenstate in consideration.
Here, $E^{(1)}_{\ell}=0$ is trivially identical for all $\ell$. Importantly,
Eqs. (2)–(4) show that each term for $E^{(2)}_{\ell}$ depends only on the
difference $(\ell-\ell^{\prime})$. Due to $2\pi$ quasienergy periodicity,
quantum numbers $\ell$ in Eqs. (2) and (3) are only defined modulo $3$. That
allows for shifting dummy indices $\ell^{\prime}$ in the summation
$E^{(2)}_{\ell_{1}}\equiv\sum^{\prime}_{\ell^{\prime},\\{n_{\boldsymbol{r}\mu}\\}}\varepsilon(\ell_{1}-\ell^{\prime})=\sum^{\prime}_{\ell^{\prime},\\{n_{\boldsymbol{r}\mu}\\}}\varepsilon(\ell_{2}-(\ell^{\prime}-\ell_{1}+\ell_{2}))=\sum^{\prime}_{\ell^{\prime\prime},\\{n_{\boldsymbol{r}\mu}\\}}\varepsilon(\ell_{2}-\ell^{\prime\prime})=E^{(2)}_{\ell_{2}}$,
proving the equality of energy corrections for all scars. SM [54] also
numerically verifies spectral pairing rigidity for Eq. (1) and against more
generic bilinear perturbations.

Importantly, it is exactly the Floquet spectrum periodicity that allows for
shifting all three $\ell$’s in Eq. (3) by the same integer and end up with an
identical set of levels, which is crucial for the above proof. In SM [54], we
prove that the spectral pairing rigidity persists to all perturbation orders
for FBS’s. Therefore, $O(L^{2})$ initial states overlapping with multiple
FBS’s separating by a rigid $\omega_{0}=2\pi/3T$ will exhibit persisting
$2\pi/\omega_{0}=3T$ DTC oscillations.

Analytical identification of FBS’s and proof for their spectral pairing
rigidity are the main results of our work. They rely on three pivotal factors.
First, strong interactions validate the starting point from Eqs. (2) and (3)
for kicked Fock states. Second, strong Floquet drivings produce three
identical $\ell=0,\pm 1$ spectral plethora at $\lambda=0$, and the $2\pi$
quasienergy periodicity intrinsic of Floquet nature enables the rigid spectral
pairing for FBS against perturbations. Third, spatiotemporal translation
symmetry prevents FBS from mutual hybridization. Therefore, FBS’s describe
genuine strongly interacting Floquet matters in clean systems.

Numerical verification — Revisiting previous numerics can now be illuminating.
Spectral function peaks in Fig. 1 (c) derive from pairs of FBS’s in Eq. (4),
$A(\omega_{0})|_{\lambda\rightarrow
0}=|\langle\boldsymbol{k},\ell_{1},N_{b}|\hat{P}|\boldsymbol{k},\ell_{1}\pm
1,N_{b}\rangle|^{2}=1/3,|\omega_{0}|=2\pi/3$. The spectral pairing rigidity
then stabilizes $|\omega_{0}|$ against perturbation up to finite size effects,
resulting in DTC oscillations in Fig. 1 (b). Also, Eq. (4) prescribes an
$L$-independent EE for FBS at $\lambda\rightarrow 0$ (see SM [54] for
analytical calculation) as in Fig. 2 (e).

(a) IPR for eigenstates

(b) Component

(c) 3 dominant components in (b).

(d) IPR${}_{\text{max}}$ scaling

Figure 3:  Structure of Floquet eigenstates in $\boldsymbol{k}=\boldsymbol{0}$
sector. Other $\boldsymbol{k}$ sectors show essentially the same results. (a)
Most eigenstates involve extensive number of basis
$|\boldsymbol{k},\\{n_{\boldsymbol{r}\mu}\\}\rangle$ leading to vanishing IPR,
except for the 3 FBSs. (b) Expand for instance one FBS in the basis
$|\boldsymbol{k},\\{n_{\boldsymbol{r}\mu}\\}\rangle$, we see it is dominated
by 3 components depicted in (c), exactly as given by Eq. (4). (d) Scaling of
the maximal IPR, where
$\Delta\text{IPR}_{\text{max}}=\text{IPR}_{\text{max}}^{(L=3)}-\text{IPR}_{\text{max}}^{(L=2)}$
in the inset. Parameters are the same as in Fig. 1, and $L=3$ for (a) (b).

Finally, we offer an efficient way to benchmark FBS by exploiting their
peculiar $\boldsymbol{k}$ space localization. A natural measure is then the
momentum space inverse participation ratio
IPR$\,=\sum_{\\{n_{\boldsymbol{r}}\\}}|\langle\boldsymbol{k},\\{n_{\boldsymbol{r}\mu}\\}|\omega_{n}\rangle|^{4}$,
where scars would show exceptionally large IPR as in Fig. 3 (a). Due to the
absence of degenerate level hybridization, the original scar components in Eq.
(4) still dominate upon perturbation as in Fig. 3 (b) and (c). The scaling of
largest IPRs in Fig. 3 (d) reproduces the reference transition
$\lambda_{0}\approx 0.135$ as in Fig. 2 (e).

Experimental relevance — Small clusters studied above can be readily realized
using the latest technology of quantum gas microscopes [55, 56, 57], which
allows for manipulation and detection with single-site resolutions. We now
further discuss cases with finite filling fractions relevant to wider ranges
of experiments.

In principle, previous analytical results show that initial states populating
more than one unit cell will chiefly overlap with non-scar ergodic
eigenstates. Therefore, a finite filling fraction among all unit cells will
eventually lead to a thermalizing behavior without dynamical signatures.
However, there could exist a finite and predictable time window before decay
to observe the scar DTCs due to scar localization.

To show it, we first take a closer look at Fig. 1 (b). The initial state of
putting $N_{b}$ bosons on one site overlaps with all FBS’s (perturbed Eq. (4))
in different $(\boldsymbol{k},\ell)$ sectors; they interfere destructively
everywhere except for the unit cell $\boldsymbol{r}=\boldsymbol{0}$, resulting
in a real-space localization. As such, two scar DTCs localized in different
regions will take time to sense the presence of and affect each other
$\phi_{2}n_{\text{ol}}(n_{\text{ol}}-1)t_{0}\sim 1$ by interactions, giving
rise to the characteristic time scale $t_{0}$ to observe DTC’s before decays.
Here $n_{\text{ol}}$ is the density overlap for two scar DTCs hypothetically
left alone in a lattice. Then, one can predict that larger distance gives a
smaller density overlap $n_{\text{ol}}$, which prolongs the scar DTC lifetime.
Such expectations are verified numerically in SM [54] for two lattice settings
relevant to the Berkeley platform. It confirms the possibility of observing
DTC signatures with finite filling fractions over the experimentally
accessible time, and further point out theoretically the controlling parameter
for DTC lifetime therein: the distance of initially populated cells.

Conclusion — We show a distinct DTC phenomenon enforced by the analytically
discovered FBS’s. Its intrinsic Floquet and many-body nature stabilizes
spectral pairings against translation-invariant bilinear perturbations.
Moreover, the new scheme of checking Floquet emergent degeneracy and scar
spectral pairing indicates a possible procedure to unveil the long-sought
universal mechanism behind clean DTCs in arbitrary dimensions. It is also
tantalizing to incorporate more intricate crystalline spacegroup symmetries
aside translations into designing DTCs with unique structures and phenomena in
clean systems.

Acknowledgment — This work is supported by the National Natural Science
Foundation of China Grant No. 12174389 (BH), the NSF Grant No. PHY-1806362
(THL, DS), the MURI-ARO Grant No. W911NF17-1-0323 through UC Santa Barbara
(BH, THL, DS, WVL), AFOSR Grant No. FA9550-16-1-0006 (WVL), and the Shanghai
Municipal Science and Technology Major Project (Grant No. 2019SHZDZX01) (WVL).

## References

  * Khemani  _et al._ [2016] V. Khemani, A. Lazarides, R. Moessner, and S. Sondhi, Phase structure of driven quantum systems, Phys. Rev. Lett. 116, 250401 (2016).
  * Else  _et al._ [2016] D. V. Else, B. Bauer, and C. Nayak, Floquet time crystals., Phys. Rev. Lett. 117, 090402 (2016).
  * Yao  _et al._ [2017] N. Y. Yao, A. C. Potter, I.-D. Potirniche, and A. Vishwanath, Discrete time crystals: Rigidity, criticality, and realizations., Phys. Rev. Lett. 118, 030401 (2017).
  * Ho  _et al._ [2017] W. W. Ho, S. Choi, M. D. Lukin, and D. A. Abanin, Critical time crystals in dipolar systems, Phys. Rev. Lett. 119, 010602 (2017).
  * Sacha [2015] K. Sacha, Modeling spontaneous breaking of time-translation symmetry, Phys. Rev. A 91, 033617 (2015).
  * Zhang  _et al._ [2017] J. Zhang, P. W. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I.-D. Potirniche, A. C. Potter, A. Vishwanath, N. Y. Yao, and C. Monroe, Observation of a discrete time crystal., Nature 543, 217 (2017).
  * Choi  _et al._ [2017] S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, C. von Keyserlingk, N. Y. Yao, E. Demler, and M. D. Lukin, Observation of discrete time-crystalline order in a disordered dipolar many-body system, Nature 543, 221 (2017).
  * Landau and Lifshitz [1980] L. D. Landau and E. M. Lifshitz,  _Statistical Physics, third edition_ (Butterworth-Heinemann, 1980).
  * Rovny  _et al._ [2018] J. Rovny, R. L. Blum, and S. E. Barrett, Observation of discrete-time-crystal signatures in an ordered dipolar many-body system, Phys Rev Lett 120, 180603 (2018).
  * Pal  _et al._ [2018] S. Pal, N. Nishad, T. Mahesh, and G. Sreejith, Temporal order in periodically driven spins in star-shaped clusters, Phys Rev Lett 120, 180602 (2018).
  * Mi and et. al. [2022] X. Mi and et. al., Time-crystalline eigenstate order on a quantum processor, Nature 601, 531 (2022).
  * Randall  _et al._ [2021] J. Randall, C. E. Bradley, F. V. van der Gronden, A. Galicia, M. H. Abobeih, M. Markham, D. J. Twitchen, F. Machado, N. Y. Yao, and T. H. Taminiau, Many-body–localized discrete time crystal with a programmable spin-based quantum simulator, Science 374, 1474 (2021).
  * Kyprianidis  _et al._ [2021] A. Kyprianidis, F. Machado, W. Morong, P. Becker, K. S. Collins, D. V. Else, L. Feng, P. W. Hess, C. Nayak, G. Pagano, N. Y. Yao, and C. Monroe, Observation of a prethermal discrete time crystal, Science 372, 1192 (2021).
  * Estarellas  _et al._ [2020] M. P. Estarellas, T. Osada, V. M. Bastidas, B. Renoust, K. Sanaka, W. J. Munro, and K. Nemoto, Simulating complex quantum networks with time crystals, Sci. Adv. 6, eaay8892 (2020).
  * Frey and Rachel [2022] P. Frey and S. Rachel, Realization of a discrete time crystal on 57 qubits of a quantum computer, Sci. Adv. 8, 10.1126/sciadv.abm7652 (2022).
  * von Keyserlingk  _et al._ [2016] C. W. von Keyserlingk, V. Khemani, and S. L. Sondhi, Absolute stability and spatiotemporal long-range order in floquet systems, Phys. Rev. B 94, 085112 (2016).
  * Huang  _et al._ [2018] B. Huang, Y.-H. Wu, and W. V. Liu, Clean floquet time crystals: Models and realizations in cold atoms, Phys. Rev. Lett. 120, 110603 (2018).
  * Russomanno  _et al._ [2017] A. Russomanno, F. Iemini, M. Dalmonte, and R. Fazio, Floquet time crystal in the lipkin-meshkov-glick model, Phys Rev B 95, 214307 (2017).
  * Zeng and Sheng [2017] T.-S. Zeng and D. N. Sheng, Prethermal time crystals in a one-dimensional periodically driven floquet system, Phys. Rev. B 96, 094202 (2017).
  * Lyu  _et al._ [2020] C. Lyu, S. Choudhury, C. Lv, Y. Yan, and Q. Zhou, Eternal discrete time crystal beating the heisenberg limit, Phys. Rev. Research 2, 033070 (2020).
  * Yu  _et al._ [2019] W. C. Yu, J. Tangpanitanon, A. W. Glaetzle, D. Jaksch, and D. G. Angelakis, Discrete time crystal in globally driven interacting quantum systems without disorder, Phys Rev A 99, 033618 (2019).
  * Mizuta  _et al._ [2018] K. Mizuta, K. Takasan, M. Nakagawa, and N. Kawakami, Spatial-translation-induced discrete time crystals, Phys Rev Lett 121, 093001 (2018).
  * Barfknecht  _et al._ [2019] R. E. Barfknecht, S. E. Rasmussen, A. Foerster, and N. T. Zinner, Realizing time crystals in discrete quantum few-body systems, Phys Rev B 99, 144304 (2019).
  * 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).
  * Else  _et al._ [2017] D. V. Else, B. Bauer, and C. Nayak, Prethermal phases of matter protected by time-translation symmetry, Phys. Rev. X 7, 011026 (2017).
  * Luitz  _et al._ [2020] D. J. Luitz, R. Moessner, S. Sondhi, and V. Khemani, Prethermalization without temperature, Phys. Rev. X 10, 021046 (2020).
  * [27] W. W. Ho and W. D. Roeck, A rigorous theory of prethermalization without temperature, arXiv:2011.14583 .
  * Khemani  _et al._ [2019a] V. Khemani, R. Moessner, and S. L. Sondhi, A brief history of time crystals (2019a).
  * Turner  _et al._ [2018] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Weak ergodicity breaking from quantum many-body scars, Nat Phys 14, 745 (2018).
  * Pizzi  _et al._ [2020] A. Pizzi, D. Malz, G. D. Tomasi, J. Knolle, and A. Nunnenkamp, Time crystallinity and finite-size effects in clean floquet systems, Phys Rev B 102, 214207 (2020).
  * Yarloo  _et al._ [2020] H. Yarloo, A. E. Kopaei, and A. Langari, Homogeneous floquet time crystal from weak ergodicity breaking, Phys Rev B 102, 224309 (2020).
  * Barter  _et al._ [2020] T. H. Barter, T.-H. Leung, M. Okano, M. Block, N. Y. Yao, and D. M. Stamper-Kurn, Spatial coherence of a strongly interacting bose gas in the trimerized kagome lattice, Phys Rev A 101, 011601 (2020).
  * Khemani  _et al._ [2019b] V. Khemani, C. R. Laumann, and A. Chandran, Signatures of integrability in the dynamics of rydberg-blockaded chains, Phys Rev B 99, 161101 (2019b).
  * 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).
  * 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).
  * Bluvstein  _et al._ [2021] D. Bluvstein, A. Omran, H. Levine, A. Keesling, G. Semeghini, S. Ebadi, T. T. Wang, A. A. Michailidis, N. Maskara, W. W. Ho, S. Choi, M. Serbyn, M. Greiner, V. Vuletić, and M. D. Lukin, Controlling quantum many-body dynamics in driven rydberg atom arrays, Science 371, 1355 (2021).
  * Maskara  _et al._ [2021] N. Maskara, A. Michailidis, W. Ho, D. Bluvstein, S. Choi, M. Lukin, and M. Serbyn, Discrete time-crystalline order enabled by quantum many-body scars: Entanglement steering via periodic driving, Phys Rev Lett 127, 090602 (2021).
  * Sugiura  _et al._ [2021] S. Sugiura, T. Kuwahara, and K. Saito, Many-body scar state intrinsic to periodically driven system, Phys. Rev. Research 3, 012010 (2021).
  * Mizuta  _et al._ [2020] K. Mizuta, K. Takasan, and N. Kawakami, Exact floquet quantum many-body scars under rydberg blockade, Phys. Rev. Research 2, 33284 (2020).
  * 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).
  * Mukherjee  _et al._ [2020] B. Mukherjee, S. Nandy, A. Sen, D. Sen, and K. Sengupta, Collapse and revival of quantum many-body scars via floquet engineering, Phys Rev B 101, 245107 (2020).
  * Desaules  _et al._ [2021] J.-Y. Desaules, A. Hudomal, C. J. Turner, and Z. Papić, Proposal for realizing quantum scars in the tilted 1d fermi-hubbard model, Phys Rev Lett 126, 210601 (2021).
  * Scherg  _et al._ [2021] S. Scherg, T. Kohlert, P. Sala, F. Pollmann, B. H. Madhusudhana, I. Bloch, and M. Aidelsburger, Observing non-ergodicity due to kinetic constraints in tilted fermi-hubbard chains, Nat Commun 12, 4490 (2021).
  * [44] G.-X. Su, H. Sun, A. Hudomal, J.-Y. Desaules, Z.-Y. Zhou, B. Yang, J. C. Halimeh, Z.-S. Yuan, Z. Papić, and J.-W. Pan, Observation of unconventional many-body scarring in a quantum simulator,  arXiv:2201.00821 .
  * Pizzi  _et al._ [2019] A. Pizzi, J. Knolle, and A. Nunnenkamp, Period-$n$ discrete time crystals and quasicrystals with ultracold bosons, Phys Rev Lett 123, 150601 (2019).
  * Lazarides and Moessner [2017] A. Lazarides and R. Moessner, Fate of a discrete time crystal in an open system, Phys Rev B 95, 195135 (2017).
  * 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).
  * Atas  _et al._ [2013] Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, Distribution of the ratio of consecutive level spacings in random matrix ensembles, Phys Rev Lett 110, 084101 (2013).
  * Regnault and Nandkishore [2016] N. Regnault and R. Nandkishore, Floquet thermalization: Symmetries and random matrix ensembles, Phys Rev B 93, 104203 (2016).
  * 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, Adv Phys 65, 239 (2016).
  * Note [1] Due to limited sizes accessible here, we would postpone a more comprehensive examination of criticality to future work and only take $\lambda_{0}$ as a reference scar vanishing point.
  * Sandvik  _et al._ [2010] A. W. Sandvik, A. Avella, and F. Mancini, Computational studies of quantum spin systems, in  _AIP Conference Proceedings_ (AIP, 2010).
  * Note [2] In strongly disordered cases [1, 16, 2], spectral pairing happens for all eigenstates due to many-body localization. Here we use the terminology to describe similar behaviors for scars but due to different reasons.
  * [54] See Supplemental Materials for details of Floquet perturbation treatment, entanglement entropy calculations, and more experimental proposals. Additional Refs. [58, 59, 60, 61, 62, 63, 64] are included therein.
  * 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).
  * Tai  _et al._ [2017] M. E. Tai, A. Lukin, M. Rispoli, R. Schittko, T. Menke, D. Borgnia, P. M. Preiss, F. Grusdt, A. M. Kaufman, and M. Greiner, Microscopy of the interacting harper–hofstadter model in the two-body limit, Nature 546, 519 (2017).
  * Semeghini  _et al._ [2021] G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, A. Omran, S. Sachdev, A. Vishwanath, M. Greiner, V. Vuletić, and M. D. Lukin, Probing topological spin liquids on a programmable quantum simulator, Science 374, 1242 (2021).
  * Thomas  _et al._ [2017] C. K. Thomas, T. H. Barter, T.-H. Leung, M. Okano, G.-B. Jo, J. Guzman, I. Kimchi, A. Vishwanath, and D. M. Stamper-Kurn, Mean-field scaling of the superfluid to mott insulator transition in a 2d optical superlattice, Phys Rev Lett 119, 100402 (2017).
  * Leung  _et al._ [2020] T.-H. Leung, M. N. Schwarz, S.-W. Chang, C. D. Brown, G. Unnikrishnan, and D. Stamper-Kurn, Interaction-enhanced group velocity of bosons in the flat band of an optical kagome lattice, Phys Rev Lett 125, 133001 (2020).
  * [60] C. D. Brown, S.-W. Chang, M. N. Schwarz, T.-H. Leung, V. Kozii, A. Avdoshkin, J. E. Moore, and D. Stamper-Kurn, Direct geometric probe of singularities in band structure, arXiv:2109.03354 .
  * Jo  _et al._ [2012] G.-B. Jo, J. Guzman, C. K. Thomas, P. Hosur, A. Vishwanath, and D. M. Stamper-Kurn, Ultracold atoms in a tunable optical kagome lattice, Phys. Rev. Lett. 108, 045305 (2012).
  * Quelle  _et al._ [2017] A. Quelle, C. Weitenberg, K. Sengstock, and C. M. Smith, Driving protocol for a floquet topological phase without static counterpart, New J Phys 19, 113010 (2017).
  * Taie  _et al._ [2015] S. Taie, H. Ozawa, T. Ichinose, T. Nishio, S. Nakajima, and Y. Takahashi, Coherent driving and freezing of bosonic matter wave in an optical lieb lattice, Sci. Adv. 1, e1500854 (2015).
  * Polkovnikov [2010] A. Polkovnikov, Phase space representation of quantum dynamics, Ann Phys-new York 325, 1790 (2010).

Supplemental materials: Discrete time crystals enforced by Floquet-Bloch scars

###### Contents

  1. S-1 Floquet perturbation treatment
    1. S-1.1 Preliminary: Factoring out perturbations
    2. S-1.2 Formalism and explicit results up to the second order
    3. S-1.3 Scar spectral pairing rigidity up to the second order
    4. S-1.4 Numerical verification
    5. S-1.5 Higher orders
  2. S-2 Entanglement entropy at the anchor point $\lambda\rightarrow 0$
  3. S-3 More details for experimental proposals
    1. S-3.1 Scheme I: Trimerized kagome lattice (TKL)
    2. S-3.2 Scheme II: Shaken honeycomb lattice
    3. S-3.3 Initial state preparation
    4. S-3.4 Detection of sublattice particle number by band mapping
    5. S-3.5 Simulations of results

##  S-1 Floquet perturbation treatment

In this section, we will derive a strongly-driven Floquet perturbation theory
for analyzing the scar stability. A concrete series up to the second order is
firstly presented for intuition. Then, we would obtain the formal structure to
all higher-orders and show the spectral pairing rigidity for FBS.

###  S-1.1 Preliminary: Factoring out perturbations

The Floquet driving in Eq. (1) of the main text has the form

$\displaystyle
U_{F}=\hat{U}_{2}\hat{U}_{1}=e^{-i\hat{H}_{2}T/2\hbar}e^{-i(\hat{h}_{1}+\lambda\hat{h^{\prime}})},$
(S1)

where $\hat{h}_{1}+\lambda\hat{h^{\prime}}=\hat{H}_{1}T/2\hbar$. For later
analysis, it will be convenient to factor out the perturbation related to
$\lambda\hat{h^{\prime}}$ through the Baker-Campbell-Hausdorff (BCH) formula,

$\displaystyle\hat{U}_{F}$
$\displaystyle=\left(e^{-i\hat{H}_{2}T/2\hbar}e^{-i\hat{h}_{1}}\right)\left(e^{i\hat{h}_{1}}e^{-i(\hat{h}_{1}+\lambda\hat{h}^{\prime})}\right)\equiv\hat{U}_{0}\hat{U}^{\prime},$
(S2) $\displaystyle\hat{U}^{\prime}$
$\displaystyle=e^{i\hat{h}_{1}}e^{-i(\hat{h}_{1}+\lambda\hat{h}^{\prime})}=e^{(i\hat{h}_{1}-i(\hat{h}_{1}+\lambda\hat{h}^{\prime}))+\frac{1}{2}[i\hat{h}_{1},-i(\hat{h}_{1}+\lambda\hat{h}^{\prime})]+\frac{1}{12}[i\hat{h}_{1},[i\hat{h}_{1},-i(\hat{h}_{1}+\lambda\hat{h}^{\prime})]]+\frac{1}{12}[-i(\hat{h}_{1}+\lambda\hat{h}^{\prime}),[-i(\hat{h}_{1}+\lambda\hat{h}^{\prime}),i\hat{h}_{1}]]+\dots}$
$\displaystyle=e^{i\lambda(-\hat{h}^{\prime}+\frac{-i}{2}[\hat{h}_{1},\hat{h}^{\prime}]+\frac{1}{6}[\hat{h}_{1},[\hat{h}_{1},\hat{h}^{\prime}]]+\dots)+O(\lambda^{2})}.$
(S3)

In our case, both $\hat{h}_{1}$ and $\hat{h}^{\prime}$ consist of bilinear
terms, so their mutual commutators to any orders are also of bilinear forms.
Meanwhile, on the exponential part of Eq. (S3), all remaining terms involve
commutators between $i\hat{h}_{1}$ and
$-i(\hat{h}_{1}+\lambda\hat{h}^{\prime}d)$, and hence are at least of the
order $\lambda^{1}$. Therefore, we can group all terms into an effective
perturbing static Hamiltonian

$\displaystyle\hat{U}^{\prime}=e^{i\hat{h}_{1}}e^{-i(\hat{h}_{1}+\lambda\hat{h}^{\prime})}=e^{i\lambda\hat{H}^{\prime}},\qquad\hat{H}^{\prime}=\sum_{\boldsymbol{r}\mu,\boldsymbol{r}^{\prime}\mu^{\prime}}J_{\boldsymbol{r}\mu,\boldsymbol{r}^{\prime}\mu^{\prime}}\hat{\psi}^{\dagger}_{\boldsymbol{r}\mu}\hat{\psi}_{\boldsymbol{r}^{\prime}\mu^{\prime}}=\hat{H}^{\prime}_{d}+\hat{H}^{\prime}_{nd}$
(S4)
$\displaystyle\hat{H}_{d}=\sum_{\boldsymbol{r}\mu}J_{\boldsymbol{r}\mu,\boldsymbol{r}\mu}\hat{\psi}_{\boldsymbol{r}\mu}^{\dagger}\hat{\psi}_{\boldsymbol{r}\mu}=\sum_{\boldsymbol{r}\mu}J_{\boldsymbol{r}\mu,\boldsymbol{r}\mu}\hat{n}_{\boldsymbol{r}\mu},\qquad\hat{H}_{nd}=\sum_{\boldsymbol{r}\mu}J_{\boldsymbol{r}\mu\neq\boldsymbol{r}^{\prime}\mu^{\prime}}\hat{\psi}_{\boldsymbol{r}\mu}^{\dagger}\hat{\psi}_{\boldsymbol{r}\mu}$
(S5)

For our purposes, we would further factor out the onsite energy offsets in
$\hat{H}^{\prime}_{d}$, as these terms can be grouped into $\hat{U}_{0}$:

$\displaystyle\hat{U}_{0}e^{i\lambda\hat{H}^{\prime}_{d}}=\hat{U}_{2}e^{-i\hat{h}_{1}}e^{i\lambda\sum_{\boldsymbol{r}\mu}J_{\boldsymbol{r}\mu,\boldsymbol{r}\mu}\hat{\psi}_{\boldsymbol{r}\mu}^{\dagger}\hat{\psi}_{\boldsymbol{r}\mu}}e^{i\hat{h}_{1}}e^{-i\hat{h}_{1}}=\hat{U}_{2}e^{i\lambda\sum_{\boldsymbol{r}\mu}J_{\boldsymbol{r}\mu,\boldsymbol{r}\mu}\hat{\psi}_{\boldsymbol{r},\mu+1}^{\dagger}\hat{\psi}_{\boldsymbol{r}\mu+1}}e^{-i\hat{h}_{1}}=\hat{\tilde{U}}_{0}.$
(S6)

Here, we use the fact that under the unperturbed $\hat{h}_{1}$,
$e^{-i\hat{h}_{1}}\hat{\psi}_{\boldsymbol{r}\mu}e^{i\hat{h}_{1}}=\hat{\psi}_{\boldsymbol{r},\mu+1}$,
and that onsite terms $\sim n_{\boldsymbol{r}\mu}$ commute with the Hubbard
interaction and chemical potentials in $\hat{U}_{2}$. Then,
$\hat{\tilde{U}}_{0}$ is related to $\hat{U}_{0}$ by a renormalized chemical
potential $\theta_{\mu}\rightarrow\theta_{\mu}-\lambda
J_{\boldsymbol{r}\mu-1,\boldsymbol{r}\mu-1}$. To factor out onsite terms, we
can perform an iteration:

  1. 1.

The perturbing Hamiltonian $\hat{H}^{\prime}$ can be obtained from Eq. (S4),
and then one could obtain the diagonal terms $\hat{H}_{d}^{\prime}$.

  2. 2.

Factor out onsite terms up to $\lambda^{1}$ by
$\hat{U}^{\prime}=e^{i\lambda\hat{H}_{d}^{\prime}}\left(e^{-i\lambda\hat{H}_{d}^{\prime}}e^{i\lambda(\hat{H}_{d}^{\prime}+\hat{H}_{nd})}\right)=e^{i\lambda\hat{H}_{d}^{\prime}}e^{i\lambda\hat{H}_{nd}^{\prime}+i\lambda^{2}(\hat{H}_{d}^{(2)}+\hat{H}_{nd}^{(2)})}$.
Here we have used the BCH formula in the last step, where
$\lambda\hat{H}_{d}^{\prime}$ in the exponential part are canceled, and the
remaining terms other than $\lambda\hat{H}^{\prime}_{nd}$ would involve
commutators between $\lambda\hat{H}_{d}^{\prime}$ and
$\lambda\hat{H}_{nd}^{\prime}$ and therefore are at least of the order
$\lambda^{2}$.

  3. 3.

Iterating the process in step 2 for $n$ times will result in
$\hat{U}^{\prime}=e^{i\left(\lambda\hat{H}_{d}+\sum_{m=2}^{n-1}\lambda^{m}\hat{H}_{d}^{(m)}\right)}e^{i\left(\lambda\hat{H}_{nd}^{\prime}+\sum_{m=2}^{n-1}\lambda^{m}\hat{H}_{nd}^{(m)}+\lambda^{n}(\hat{H}_{d}^{(n)}+\hat{H}_{nd}^{(n)})\right)}$.
That means one can factor out the diagonal energy offset terms up to arbitrary
accuracy.

In summary, for a bilinear type of Hamiltonian, one can transform the
perturbation in Eq. (1) of the main text concerning $\hat{H}_{1}$ such that
the Floquet unitary reads $\hat{U}_{F}=\hat{U}_{0}\hat{U}^{\prime}$. Here,
$\hat{U}_{0}$ is the unperturbed Eq. (1) with $\lambda=0$ and onsite chemical
potentials renormalized. $\hat{U}^{\prime}=e^{i\lambda\hat{H}^{\prime}}$
contains all the perturbations with
$\hat{H}^{\prime}=\sum_{\boldsymbol{r}\mu\neq\boldsymbol{r}\mu}J_{\boldsymbol{r}\mu,\boldsymbol{r}\mu}\hat{\psi}_{\boldsymbol{r}\mu}^{\dagger}\hat{\psi}_{\boldsymbol{r}^{\prime}\mu^{\prime}}$.

###  S-1.2 Formalism and explicit results up to the second order

Consider a Floquet operator involving a strongly-driven but exactly solvable
part $U_{0}$, and a perturbation $U_{1}$,

$\displaystyle U_{F}=U_{0}U^{\prime},\qquad
U_{0}|\omega_{n}\rangle=e^{i\omega_{n}}|\omega_{n}\rangle,\quad
U^{\prime}=e^{i\lambda
H^{\prime}}=\sum_{\alpha=0}^{\infty}\frac{(i\lambda)^{\alpha}}{\alpha!}(H^{\prime})^{\alpha}$
(S7)

Note that $U_{0}$ could be far from an identity operator unlike conventional
high-frequency expansions around static limits. Now, to solve the eigenproblem
perturbatively,

$\displaystyle
U_{F}|\tilde{\omega}_{n}\rangle=e^{i\tilde{\omega}_{n}}|\tilde{\omega}_{n}\rangle,$
(S8)

we expand the quasienergy and Floquet eigenstates into series

$\displaystyle
e^{i\tilde{\omega}_{n}}=e^{i\omega_{n}}\exp\left(i\sum_{\alpha=1}^{\infty}\lambda^{\alpha}\omega_{n}^{(\alpha)}\right),\qquad\qquad|\tilde{\omega}_{n}\rangle$
$\displaystyle=\dots e^{i\lambda^{3}S_{3}}e^{i\lambda^{2}S_{2}}e^{i\lambda
S_{1}}|\omega_{n}\rangle,$ (S9)

such that $e^{i\lambda^{\alpha}S_{\alpha}}$ should diagonalize $U_{F}$ up to
the order of $\lambda^{\alpha}$. Note that unlike previous perturbation
treatments [2] expressing $(\alpha+1)$-th order results by $\alpha$-th order
ones, we would aim at obtaining series of any orders in terms of the zeroth-
order solutions $\omega_{n},|\omega_{n}\rangle$. This is necessary for
quantitative evaluations like spectral pairing rigidity, as only the zeroth-
order results are exactly solvable.

To $\lambda^{1}$ order,

$\displaystyle\langle\omega_{m}|e^{-i\lambda S_{1}}U_{0}e^{i\lambda
H^{\prime}}e^{i\lambda S_{1}}|\omega_{n}\rangle=$
$\displaystyle\,Q^{(1)}_{n}\left(\langle\omega_{m}|U_{0}|\omega_{n}\rangle+i\lambda\langle\omega_{m}|(U_{0}H^{\prime}+[U_{0},S_{1}])|\omega_{n}\rangle\right)+O(\lambda^{2})$
$\displaystyle=$
$\displaystyle\,Q^{(1)}_{n}\left(e^{i\omega_{m}}\delta_{mn}+i\lambda(e^{i\omega_{m}}[H^{\prime}]_{mn}+(e^{i\omega_{m}}-e^{i\omega_{n}})[S_{1}]_{mn})\right)+O(\lambda^{2}).$
$\displaystyle=$
$\displaystyle\,Q^{(1)}_{n}e^{i\omega_{m}}\left(\delta_{mn}+i\lambda([H^{\prime}]_{mn}+(1-e^{i\omega_{nm}})[S_{1}]_{mn})\right)+O(\lambda^{2}).$
(S10)

where $\omega_{nm}=\omega_{n}-\omega_{m}$ and matrix elements i.e.
$[S_{1}]_{mn}=\langle\omega_{m}|S_{1}|\omega_{n}\rangle$. The normalization
constant $Q^{(1)}_{n}$ makes sure that the diagonal term ($m=n$) on the right-
hand-side for perturbed eigenvalue, after keeping terms only up to
$\lambda^{1}$, is still unitary. To diagonalize $U_{F}$ up to $\lambda^{1}$,
choose

$\displaystyle[S_{1}]_{mn}=\frac{[H^{\prime}]_{mn}}{e^{i\omega_{nm}}-1},\quad(m\neq
n);\qquad\qquad[S_{1}]_{nn}=0.$ (S11)

Then, compare the diagonal term in Eq. (S10) with Eq. (S9), we have a similar
result as in static case

$\displaystyle\omega_{n}^{(1)}=[H^{\prime}]_{nn}.$ (S12)

A qualitative difference from static perturbation can be observed in the
$\lambda^{2}$ order,

$\displaystyle\langle\omega_{m}|e^{-i\lambda
S_{1}}e^{-i\lambda^{2}S_{2}}U_{0}e^{i\lambda
H^{\prime}}e^{i\lambda^{2}S_{2}}e^{i\lambda S_{1}}|\omega_{n}\rangle$
$\displaystyle=$
$\displaystyle\,Q_{n}^{(2)}\left(e^{i\omega_{n}}(1+i\lambda[H^{\prime}]_{nn})\delta_{mn}+\lambda^{2}\langle\omega_{m}|[-\frac{1}{2}(U_{0}S_{1}^{2}+S_{1}^{2}U_{0})+i[U_{0},S_{2}]-\frac{1}{2}U_{0}(H^{\prime})^{2}+[S_{1},U_{0}H^{\prime}]+S_{1}U_{0}S_{1}]|\omega_{n}\rangle\right)+O(\lambda^{3}).$
(S13)

To diagonalize it, we similarly require the $\lambda^{2}$ term for off-
diagonal elements $m\neq n$ to vanish,

$\displaystyle 0=$
$\displaystyle\,\lambda^{2}\left(-\frac{[S_{1}^{2}]_{mn}}{2}(e^{i\omega_{m}}+e^{i\omega_{n}})+i[S_{2}]_{mn}(e^{i\omega_{m}}-e^{i\omega_{n}})-\frac{[(H^{\prime})^{2}]_{mn}}{2}e^{i\omega_{m}}-[H^{\prime}S_{1}]_{mn}e^{i\omega_{m}}+\sum_{l}[S_{1}]_{ml}[H^{\prime}]_{ln}e^{i\omega_{l}}+\sum_{l}[S_{1}]_{ml}[S_{1}]_{ln}e^{i\omega_{l}}\right),$

A straightforward calculation simplifies it to

$\displaystyle i[S_{2}]_{mn}(e^{i\omega_{nm}}-1)$
$\displaystyle=\sum_{l}\frac{[H^{\prime}]_{ml}[H^{\prime}]_{ln}}{(e^{i\omega_{lm}}-1)(e^{i\omega_{nl}}-1)}\left(-\frac{1+e^{i\omega_{nm}}}{2}-\frac{(e^{i\omega_{lm}}-1)(e^{i\omega_{nl}}-1)}{2}-(e^{i\omega_{lm}}-1)+e^{i\omega_{lm}}(e^{i\omega_{nl}}-1)+e^{i\omega_{lm}}\right)$
$\displaystyle=\sum_{l}\frac{[H^{\prime}]_{ml}[H^{\prime}]_{ln}}{(e^{i\omega_{lm}}-1)(e^{i\omega_{nl}}-1)}\frac{e^{i\omega_{nl}}-e^{i\omega_{lm}}}{2}$
(S14)

Therefore, we can similarly choose

$\displaystyle[S_{2}]_{mn}$
$\displaystyle=\frac{-i}{e^{i\omega_{nm}}-1}\sum_{l}\frac{[H^{\prime}]_{ml}[H^{\prime}]_{ln}}{(e^{i\omega_{lm}}-1)(e^{i\omega_{nl}}-1)}\frac{e^{i\omega_{nl}}-e^{i\omega_{lm}}}{2},\quad(m\neq
n);\qquad\qquad[S_{2}]_{mm}=0,$ (S15)

Then, the diagonal term of Eq. (S13), where $\lambda^{2}$ terms corresponding
to the right-hand-side of Eq. (S14) with $m=n$, gives the second order
correction to quasienergy

$\displaystyle{\color[rgb]{0,0,1}\omega_{n}^{(2)}}$
$\displaystyle=-\sum_{l\neq
n}\left|\frac{[H^{\prime}]_{nl}}{(e^{i\omega_{ln}}-1)}\right|^{2}\sin\omega_{ln}={\color[rgb]{0,0,1}-\frac{1}{2}\sum_{l\neq
n}|[H^{\prime}]_{nl}|^{2}\cot\frac{\omega_{ln}}{2}}$ (S16)

We can benchmark the results by considering high frequency limits where all
$\omega_{nl}$’s are very small. Then, we can expand
$\cot(\omega_{nl}/2)\approx 2/\omega_{nl}$, so $\omega_{n}^{(2)}=\sum_{l\neq
n}|[H^{\prime}]_{nl}|^{2}/\omega_{nl}$, which recovers the familiar static
perturbation result.

It is worth noting that in the Floquet case, level differences $\omega_{ln}$
contribute a periodic $2\pi$ correction $\cot(\omega_{ln}/2)$ in contrast to
the $1/\omega_{ln}$ factor for static cases without periodicity. Such a
difference is crucial for the spectral pairing rigidity only possible in
Floquet systems.

In sum, up to the second order, the perturbed quasienergy reads

$\displaystyle{\color[rgb]{0,0,1}e^{\tilde{\omega}_{n}}\approx
e^{i\omega_{n}}\exp\left(i\lambda[H^{\prime}]_{nn}-i\frac{\lambda^{2}}{2}\sum_{l\neq
n}|[H^{\prime}]_{nl}|^{2}\cot\frac{\omega_{ln}}{2}\right),}$ (S17)

For completeness, we also write the dressed eigenstates

$\displaystyle|\tilde{\omega}_{n}\rangle$
$\displaystyle\approx(1+i\lambda^{2}S_{2})(1+i\lambda
S_{1}-\frac{\lambda^{2}}{2}S_{1}^{2})|\omega_{n}\rangle$
$\displaystyle=|\omega_{n}\rangle+i\lambda\sum_{m\neq
n}[S_{1}]_{mn}|\omega_{m}\rangle+\lambda^{2}\left(\sum_{m}\sum_{l\neq
m,n}-\frac{[S_{1}]_{ml}[S_{1}]_{ln}}{2}+i[S_{2}]_{mn}\right)|\omega_{m}\rangle.$
$\displaystyle=|\omega_{n}\rangle+i\lambda\sum_{m\neq
n}\frac{[H^{\prime}]_{mn}}{e^{i\omega_{mn}}-1}|\omega_{m}\rangle-\frac{\lambda^{2}}{2}\left(\sum_{m}\sum_{l\neq
m,n}[S_{1}]_{ml}[S_{1}]_{ln}+\sum_{m\neq n}\sum_{l\neq
m,n}\frac{e^{i\omega_{ml}}-e^{i\omega_{ln}}}{e^{i\omega_{mn}}-1}[S_{1}]_{ml}[S_{1}]_{ln}\right)|\omega_{m}\rangle$
(S18)

up to a normalization factor. Here $[S_{1}]_{mn}$ etc are given by Eq. (S11)

###  S-1.3 Scar spectral pairing rigidity up to the second order

Now, we apply the previously developed strong-drive Floquet perturbation
theory to analyze the scar spectral pairings. Recall Eq. (3) (4) in the main
text for the eigenstates at $\lambda=0$, which we will take as eigenstates for
$U_{0}$,

$\displaystyle\left|\boldsymbol{k},\ell,\\{n_{\boldsymbol{r}\mu}\\}\right\rangle$
$\displaystyle=\,\frac{1}{\sqrt{3}}\sum_{m=0,1,2}e^{-i(\frac{2\pi
m}{3}\ell-\alpha_{m})}|\boldsymbol{k},\\{n_{\boldsymbol{r},\mu+m\text{ mod
}3}\\}\rangle,$ (2) $\displaystyle
E\left(\ell,\\{n_{\boldsymbol{r}\mu}\\}\right)$
$\displaystyle=\frac{2\pi}{3}\ell+\phi_{2}\sum_{\boldsymbol{r}\mu}n_{\boldsymbol{r}\mu}(n_{\boldsymbol{r}\mu}-1),\qquad\ell=0,\pm
1.$ (3)

Each level labeled by
$|\boldsymbol{k},\ell,\\{n_{\boldsymbol{r}\mu}\\}\rangle$ corresponds to the
zeroth order $|\omega_{n}\rangle$ here, and quasienergy
$e^{iE(\ell,\\{n_{\boldsymbol{r}\mu}\\})}$ is denoted $e^{i\omega_{n}}$. Among
them, the special set of scars are defined in Eq. (4) of the main text,

$\displaystyle|\boldsymbol{k},\ell,N_{b}\rangle=\frac{1}{\sqrt{3}}\sum_{m=0,1,2}e^{-i(\frac{2\pi
m}{3}\ell-\alpha_{m})}|\boldsymbol{k},\\{n_{\boldsymbol{r}\mu}=\delta_{\boldsymbol{r},\boldsymbol{0}}\delta_{\mu,m}N_{b}\\}\rangle,\qquad
E_{\text{scar}}(\ell)=\frac{2\pi\ell}{3}+\phi_{2}N_{b}(N_{b}-1),\qquad\ell=0,\pm
1.$ (4)

In this subsection, we would aim at proving that all the three scars
$\ell=0,\pm$ in each $\boldsymbol{k}$ sector will receive the same quasienergy
corrections in Eq. (S17), with the perturbing Hamiltonian
$H^{\prime}=\sum_{\boldsymbol{r}\mu\neq\boldsymbol{r}^{\prime}\mu^{\prime}}J_{\boldsymbol{r}\mu,\boldsymbol{r}^{\prime}\mu^{\prime}}\hat{\psi}^{\dagger}_{\boldsymbol{r}\mu}\hat{\psi}_{\boldsymbol{r}^{\prime}\mu^{\prime}}$
being a generic bilinear hopping one respecting translation invariance.

First, as onsite chemical potential terms are all factored out into $U_{0}$
(see subsection S-1.1), the first order corrections
$E_{\ell}^{(1)}\propto\langle\boldsymbol{k},\ell,N_{b}|H^{\prime}|\boldsymbol{k},\ell,N_{b}\rangle=0$
vanish — in $|\boldsymbol{k},\ell,N_{b}\rangle$ (for $N_{b}>1$) all particles
are allocated onto a single site and no hopping happens. Next, for the second
order corrections, we write explicitly

$\displaystyle
E_{\ell}^{(2)}=-\frac{1}{2}\sum_{\\{n^{\prime}_{\boldsymbol{r}\mu}\\},\ell^{\prime}}|\langle\boldsymbol{k},\ell,N_{b}|H^{\prime}|\boldsymbol{k},\ell^{\prime},\\{n^{\prime}_{\boldsymbol{r}\mu}\\}\rangle|^{2}\cot\left(\frac{E(\ell^{\prime},\\{n^{\prime}_{\boldsymbol{r}\mu}\\})-E_{\text{scar}}(\ell)}{2}\right).$
(S19)

As the perturbing Hamiltonian only involves hopping terms, the non-vanishing
matrix elements have $\\{n_{\boldsymbol{r}\mu}^{\prime}\\}$ with $(N_{b}-1)$
particles on one site, and $1$ particle on another. Written in the main text
notation, $\\{n_{\boldsymbol{r}\mu}^{\prime}\\}\in
Q_{a}=\\{(q^{(a)}_{1},N^{(a)}_{1})=(1,N_{b}-1),(q^{(a)}_{2},N^{(a)}_{2})=(1,1),(q^{(a)}_{3},N^{(a)}_{3})=(3L^{2}-2,0)\\}$.
Thus,
$E(\ell^{\prime},\\{n^{\prime}_{\boldsymbol{r}\mu}\\})=2\pi\ell^{\prime}/3+\phi_{2}\sum_{j}q^{(a)}_{j}N^{(a)}_{j}(N^{(a)}_{j}-1)=2\pi\ell^{\prime}/3+\phi_{2}(N_{b}-1)(N_{b}-2)$,
and all

$\displaystyle
E(\ell^{\prime},\\{n_{\boldsymbol{r}\mu}\\})-E_{\text{scar}}(\ell)=-(2\pi(\ell-\ell^{\prime})/3+2\phi_{2}).$
(S20)

Further, within the manifold $\\{n_{\boldsymbol{r}\mu}^{\prime}\\}\in Q_{a}$,

$\displaystyle\langle\boldsymbol{k},\ell,N_{b}|H^{\prime}|\boldsymbol{k},\ell^{\prime},\\{n^{\prime}_{\boldsymbol{r}\mu}\\}\rangle$
$\displaystyle=\frac{1}{3}\sum_{m,m^{\prime}=0}^{2}\langle\boldsymbol{k},\\{n_{\boldsymbol{r}\mu}=\delta_{\boldsymbol{r},\boldsymbol{0}}\delta_{\mu,m}N_{b}\\}|H^{\prime}|\boldsymbol{k},\\{n^{\prime}_{\boldsymbol{r}\mu+m^{\prime}}\\}\rangle
e^{i\frac{2\pi}{3}(m\ell-m^{\prime}\ell^{\prime})}e^{-i(\alpha_{m}-\tilde{\alpha}_{m^{\prime}})},$
(S21)

the site with $N_{b}$ particles for $\\{n_{\boldsymbol{r}\mu}\\}$ must be the
same as the site with $(N_{b}-1)$ particles for
$\\{n_{\boldsymbol{r}\mu}^{\prime}\\}$, see Fig. S1 for example.

Figure S1: Exemplary non-vanishing matrix elements in the second order
perturbation for scar states $|\boldsymbol{k},\ell,N_{b}\rangle$. The non-
vanishing matrix elements would have $m^{\prime}=m+\Delta m$ mod $3$, where
$\Delta m=2$ for all $m,m^{\prime}$.

That means a non-vanishing term would have
$\langle\boldsymbol{k},\\{n_{\boldsymbol{r}\mu}=\delta_{\boldsymbol{r},\boldsymbol{0}}\delta_{\mu,m}N_{b}\\}|H^{\prime}|\boldsymbol{k},\\{n^{\prime}_{\boldsymbol{r}\mu+m^{\prime}}\\}\rangle\propto\delta_{m^{\prime},m+\Delta
m}$, with a fixed $\Delta m$ (mod 3) for all $m,m^{\prime}$. Then,

$\displaystyle|\langle\boldsymbol{k},\ell,N_{b}|H^{\prime}|\boldsymbol{k},\ell^{\prime},\\{n^{\prime}_{\boldsymbol{r}\mu}\\}\rangle|^{2}$
$\displaystyle=\left|\frac{1}{3}e^{-i\frac{2\pi\Delta
m}{3}\ell^{\prime}}\sum_{m=0}^{2}\langle\boldsymbol{k},\\{n_{\boldsymbol{r}\mu}=\delta_{\boldsymbol{r},\boldsymbol{0}}\delta_{\mu,m}N_{b}\\}|H^{\prime}|\boldsymbol{k},\\{n^{\prime}_{\boldsymbol{r}\mu+m+\Delta
m}\\}\rangle
e^{i\frac{2\pi}{3}m(\ell-\ell^{\prime})}e^{-i(\alpha_{m}-\tilde{\alpha}_{m^{\prime}})}\right|^{2}.$

In sum, the second order quasi-energy correction

$\displaystyle E_{\ell}^{(2)}$
$\displaystyle=\frac{1}{18}\sum_{\\{n_{\boldsymbol{r}\mu}^{\prime}\\},\ell^{\prime}}\left|\sum_{m=0}^{2}\langle\boldsymbol{k},\\{n_{\boldsymbol{r}\mu}=\delta_{\boldsymbol{r},\boldsymbol{0}}\delta_{\mu,m}N_{b}\\}|H^{\prime}|\boldsymbol{k},\\{n^{\prime}_{\boldsymbol{r}\mu+m+\Delta
m}\\}\rangle
e^{i\frac{2\pi}{3}m{\color[rgb]{1,0,0}(\ell-\ell^{\prime})}}e^{-i(\alpha_{m}-\tilde{\alpha}_{m^{\prime}})}\right|^{2}\cot\left(\frac{\pi{\color[rgb]{1,0,0}(\ell-\ell^{\prime})}}{3}+\phi_{2}\right)$
(S22)

indeed only depend on the difference $(\ell-\ell^{\prime})$, but not
individual $\ell$ and $\ell^{\prime}$. Thus, we can generically write
$E_{\ell_{1}}^{(2)}\equiv\sum_{\ell^{\prime}}\varepsilon(\ell_{1}-\ell^{\prime})=\sum_{\ell^{\prime}}\varepsilon(\ell_{2}-(\ell^{\prime}-\ell_{1}+\ell_{2}))=\sum_{\ell^{\prime\prime}}\varepsilon(\ell_{2}-\ell^{\prime\prime})=E_{\ell_{2}}^{(2)}$,
proving the spectral pairing rigidity up to the second order.

###  S-1.4 Numerical verification

Now we verify the previous analysis numerically. First, we test the model in
Eq. (1) of the main text. Quasienergy for the 3 scar states, obtained
according to maximal momentum space IPRs as in Fig. 3 of the main text, is
shown in Fig. S2 (a). At the anchor point $\lambda\rightarrow 0$, mutual
spacing for FBS’s approaches $2\pi/3$ giving the $3T$-periodic oscillations.
Then, under perturbations $\lambda$, each individual FBS indeed receive an
energy correction dominated by $\sim\lambda^{2}$ as expected. However, the key
feature is that all three FBS’s receive identical quasienergy corrections, as
shown by Fig. S2 (b), leading to a rigid $2\pi/3$ spectral pairing between
pairs of FBS’s schematically illustrated in Fig. S2 (c).

(a) Quasienergy for largest IPR eigenstates

(b) Deviations and spectral pairing rigidity

(c) Spectral pairing for FBS’s

Figure S2:  Test of spectral pairing for the main text model. (a) Based on the
data of Fig. (3) in the main text, we obtain the quasi-energy of the 3
eigenstates with largest IPR’s. Below the transition $\lambda\approx 0.135$
indicated by the green dashed line, they correspond to the FBS’s, while for
$\lambda>0.135$ all eigenstates have similarly vanishing IPR’s. We take
$\lambda=0.005$ as a reference point for $\lambda\rightarrow 0$ limits
(corresponding to gray horizontal lines), and measure quasienergy deviations
$\Delta E_{n}\equiv E_{n}(0.005)-E_{n}(\lambda)$. (b) The spectral deviations
of each FBS. We see that individual scar quasienergy $\Delta E_{n}$ indeed
exhibits a $\sim\lambda^{2}$ deviation of Fermi-Golden rule type. However,
different scars demonstrate almost identical deviations $\Delta
E_{n_{1}}-\Delta E_{n_{2}}\rightarrow 0$, such that their mutual spacing of
$2\pi/3$ in (a) is preserved, leading to a stable $2\pi T/(2\pi/3)=3T$
periodic DTC oscillations. (c) Schematic illustrate for scar spectral
pairings. All data is for the $\boldsymbol{k}=\boldsymbol{0}$ sector, and
other $\boldsymbol{k}$ sectors demonstrate essentially the same characters.
$L=3$, and parameters are the same as Fig. 1 in the main text.

Further, we test the spectral pairing rigidity against more generic
perturbations using the second order perturbation results in Eq. (S16). Here
we generalize $H^{\prime}$ to involve all possible hoppings up to nearest
neighbors, where each bond can possess different hopping matrix elements as
shown by Fig. S3 (a). These bonds possess random strengths $J_{n}\in[0,1]$ and
also carry random fluxes $\Theta_{n}\in[0,2\pi]$. Each set of
$\\{(J_{n},\Theta_{n})|n=1,\dots,6)\\}$ consists of a sample. From the results
in Fig. S3 (b) for different sets of samples, it is clear that each scar
quasienergy can receive notable second order corrections $\omega_{n}^{(2)}$ as
given by Eq. (S16). However, all 3 scars receive equal amount of corrections
$|\omega_{n+1}^{(2)}-\omega_{n}^{(2)}|\rightarrow 0$ (up to numerical errors),
such that their mutual spectral pairing remain rigidly $2\pi/3$ just like that
for $U_{0}$, reproducing again the scheme in Fig. S2 (c).

(a) Random bonds

(a) Second order quasienergy corrections (in units of $\lambda^{2}$)

Figure S3: Generalized model allowing for random bond parameters (preserving
translation symmetry) up to nearest neighbor hoppings. Floquet operator takes
the form in Eq. (S7), where $N_{b}=5,L=3$, $U_{0}$ is for the main text model
at the anchor point $\phi_{1}=2\pi/(3\sqrt{3}),\lambda=0$, and the
perturbation $U^{\prime}=e^{i\lambda H^{\prime}}$ . (a) Random hopping
parameters for $H^{\prime}$. (b) Results for the
$\boldsymbol{k}=\boldsymbol{0}$ sector. Each sample means one set of
parameters in (a). To compare with Fig. S2, note the rescaling for each point,
$\omega_{n}^{(2)}\sim\Delta E_{n}/\lambda^{2}$, where $\omega_{n}^{(2)}$ is
given by Eq. (S16).

###  S-1.5 Higher orders

Formally, as can be seen from Eq. (S17) and Eqs. (S7) – (S9), the perturbed
quasi-energy at the $\alpha$-th order takes the most generic form satisfying
Floquet quasienergy periodicity as

$\displaystyle\omega_{n}^{(\alpha)}=\sum_{\beta}\left.\sum_{l_{1}\dots
l_{\alpha-1}}\right.^{\prime}g_{\beta}([H^{\prime}]_{nl_{1}}[H^{\prime}]_{l_{1}l_{2}}\dots[H^{\prime}]_{l_{\alpha-2}l_{\alpha-1}}[H^{\prime}]_{l_{\alpha-1}n})h_{\beta}(e^{i\omega_{nl_{1}}},e^{i\omega_{l_{1}l_{2}}},\dots,e^{i\omega_{l_{\alpha-1}n}}),$
(S23)

where $\beta=1,2,\dots$ denotes a set of functions. The major difference
between the second and higher orders is that one would encounter, i.e. $l_{1}$
and $l_{2}$ denoting unperturbed eigenstates within the same degenerate
manifold. To be concrete, we give the explicit third order result

$\displaystyle\omega_{n}^{(3)}$ $\displaystyle=\frac{1}{3!}\sum_{l_{1}\neq
l_{2}\neq
n}[H^{\prime}]_{nl_{1}}[H^{\prime}]_{l_{1}l_{2}}[H^{\prime}]_{l_{2}n}\left[\frac{1}{2}\frac{\cos(\omega_{l_{1}n}/2)}{\sin(\omega_{l_{2}l_{1}}/2)\sin(\omega_{nl_{2}}/2)}+\frac{1}{2}\frac{\cos(\omega_{nl_{2}}/2)}{\sin(\omega_{l_{1}n}/2)\sin(\omega_{l_{2}l_{1}}/2)}-\frac{\cos(\omega_{l_{2}l_{1}}/2)}{\sin(\omega_{l_{1}n}/2)\sin(\omega_{nl_{2}}/2)}\right]$
(S24)

Here, two different levels denoted by $l_{1},l_{2}$ could be different
configurations in the same degenerate manifold
$\ell,Q_{a}=\\{(q^{(a)}_{1},N^{(a)}_{1})=(1,N_{b}-1),(q^{(a)}_{2},N^{(a)}_{2})=(1,1),(q^{(a)}_{3},N^{(a)}_{3})=(3L^{2}-2,0)\\}$,
so the two different levels $l_{1},l_{2}$ both have unperturbed quasienergy
$2\pi\ell/3+\sum_{j}q_{j}N_{j}(N_{j}-1)=2\pi\ell/3+(N_{b}-1)(N_{b}-2)$. As
such, we need to first perform a degenerate level perturbation for non-scar
eigenstates within each subspace $\ell,Q_{a}$.

Figure S4:  Non-vanishing submatrix elements within the
$Q_{a}=\\{(q^{(a)}_{j},N^{(a)}_{j})=(1,N_{b}-1),(1,1),(3L^{2}-2,0)\\}$
degenerate manifolds. Similar to the case in Fig. S1, the site with majority
particle $N_{b}-1$ forces the matrix element to only depend on the relative
$\ell-\ell=0$ — so there is no dependence on $\ell$ for degenerate sub-space
corrections.

Importantly, the submatrix spanned by
$\\{|\boldsymbol{k},\ell,\\{n_{\boldsymbol{r}\mu}\\}\rangle|\\{n_{\boldsymbol{r}\mu}\in
Q_{a}\\}\\}$ would contribute a degeneracy-lifting energy independent of
$\ell$. This is due to the same reason as before: the site with $N_{b}-1$
particles must be the same in order for the hopping matrix element to be non-
vanishing. One can choose the gauge that the unperturbed eigenstates within
$Q_{a}$ have the $(N_{b}-1)$-particle site being the same for each $m$ in
different $\\{n_{\boldsymbol{r}\mu}\\}$, as shown in Fig. S4. Then, the
degenerate space corrections in different $\ell$ sectors are the same. That
means for degenerate eigenstates within the manifold
$Q_{a}=\\{(q_{j}^{(a)},N_{j}^{(a)})|j=1,2,\dots,M\\}$, one only needs to
replace Eq. (2), (3) in the main text with

$\displaystyle|\boldsymbol{k},\ell,\\{\\{n_{\boldsymbol{r}\mu}\\}\in
Q_{a}\\},\gamma\rangle=\sum_{\\{n_{\boldsymbol{r}\mu}\\}\in
Q_{a}}A_{\\{n_{\boldsymbol{r}\mu}\\}}^{(\gamma)}|\boldsymbol{k},\ell,\\{n_{\boldsymbol{r}\mu}\\}\rangle=\frac{1}{\sqrt{3}}\sum_{m=0,1,2}e^{-i\frac{2\pi
m}{3}\ell}\left(\sum_{\\{n_{\boldsymbol{r}\mu}\\}\in
Q_{a}}A_{\\{n_{\boldsymbol{r}\mu}\\}}^{(\gamma)}e^{i\alpha_{m}}|\boldsymbol{k},\\{n_{\boldsymbol{r}\mu}\\}\rangle\right)$
($2^{\prime}$) $\displaystyle E\left(\ell,\\{\\{n_{\boldsymbol{r}\mu}\\}\in
Q_{a}\\},\gamma\right)=\frac{2\pi}{3}\ell+\phi_{2}\sum_{j=1}^{M}q_{j}N_{j}(N_{j}-1)+E_{\gamma}(\\{n_{\boldsymbol{r}\mu}\\}\in
Q_{a}).$ ($3^{\prime}$)

Due to the absence of $\ell$ in the submatrix within degenerate $Q_{a}$
manifold, the quasienergies in each $\ell$ sector are lifted by the
$\ell$-independent $E_{\gamma}$, and the coefficients
$A_{\\{n_{\boldsymbol{r}\mu}\\}}^{(\gamma)}$ are also independent of $\ell$.

Now, as all degeneracies are lifted, and the resulting levels still has the
structure of identical plethoras of $\ell=0,\pm 1$, we can use the previous
analysis to prove spectral pairing rigidity. Specifically, in Eq. (S23), when
we change the scar level

$\displaystyle n\sim(\ell,N_{b})\quad\rightarrow\quad n^{\prime}\sim(\ell\pm
1,N_{b})$ (S25)

the right-hand-side should remain the same. Specifically, we note that quantum
numbers $\ell=0,\pm 1$ is only defined modulo 3. That means summing over
$\sum_{\ell=-1,0,1}$ is equivalent to
$\sum_{\ell+1=0,1,2}=\sum_{\ell+1=0,1,-1}$ and
$\sum_{\ell-1=-2,-1,0}=\sum_{1,-1,0}$. Then, Eq. (S25) would be compensated by
a simultaneous change of the dummy index in Eq. (S23)

$\displaystyle
l_{1}\sim(\ell_{1},\\{n_{\boldsymbol{r}\mu}\\}_{1})\quad\rightarrow\quad
l_{1}^{\prime}\sim(\ell_{1}\pm 1,\\{n_{\boldsymbol{r}\mu}\\}_{1}),\qquad\qquad
l_{\alpha-1}\sim(\ell_{\alpha-1},\\{n_{\boldsymbol{r}\mu}\\}_{\alpha-1})\quad\rightarrow\quad
l_{\alpha-1}^{\prime}\sim(\ell_{\alpha-1}\pm
1,\\{n_{\boldsymbol{r}\mu}\\}_{\alpha-1}),$ (S26)

while all the remaining indices are unaffected.

In summary, we have proved that for perturbations $H^{\prime}$ of a generic
bilinear form conserving translation symmetries, all scar levels
$|\boldsymbol{k},\ell,N_{b}\rangle$ will receive the same amount of energy
correction in the perturbation series, leading to the spectral pairing
rigidity $|\Delta E|=2\pi/3$ for FBS’s.

##  S-2 Entanglement entropy at the anchor point $\lambda\rightarrow 0$

To compute the entanglement entropy, we first rewrite Eq. (4) in the main text
in the real space representation

$\displaystyle|\boldsymbol{k},\ell,N_{b}\rangle=\frac{1}{\sqrt{3}}\sum_{m=0,1,2}e^{-i(\frac{2\pi
m}{3}\ell-\alpha_{m})}\frac{1}{L}\sum_{m_{1},m_{2}=1}^{L}e^{-(2\pi
i/L)(k_{1}m_{1}+k_{2}m_{2})}\frac{(\hat{\psi}_{m_{1}\boldsymbol{e}_{1}+m_{2}\boldsymbol{e}_{2},m}^{\dagger})^{N_{b}}}{\sqrt{N_{b}!}}|0\rangle$
(S27)

The full density matrix for each FBS is defined as
$\rho=|\boldsymbol{k},\ell,N_{b}\rangle\langle\boldsymbol{k},\ell,N_{b}|$, and
in the real space representation

$\displaystyle\rho$
$\displaystyle=\frac{1}{3L^{2}}\sum_{m,m^{\prime}=0,1,2}e^{i(\frac{2\pi(m-m^{\prime})}{3}\ell-(\alpha_{m}-\alpha_{m^{\prime}}))}\sum_{m_{1},m_{2},m_{1}^{\prime},m_{2}^{\prime}=1}^{L}e^{\frac{2\pi
i}{L}(k_{1}(m_{1}-m_{1}^{\prime})+k_{2}(m_{2}-m_{2}^{\prime}))}\left(\frac{(\hat{\psi}_{m_{1}\boldsymbol{e}_{1}+m_{2}\boldsymbol{e}_{2},m}^{\dagger})^{N_{b}}}{\sqrt{N_{b}!}}|0\rangle\langle
0|\frac{(\hat{\psi}_{m_{1}^{\prime}\boldsymbol{e}_{1}+m_{2}^{\prime}\boldsymbol{e}_{2},m^{\prime}})^{N_{b}}}{\sqrt{N_{b}!}}\right)$
(S28)

Then, we obtain the reduced density matrix by enclosing $N_{s}$ subsystem unit
cells ($3N_{s}$ sites), dubbed region $A$, among totally $3L^{2}$ sites. The
remaining part is denoted as $B$. Due to the form of Eq. (S28), $\rho_{A}$
only involves tracing over configurations of $0$ or $N_{b}$ bosons in region
$B$,

$\displaystyle\rho_{A}=\text{Tr}_{B}(\rho)$ $\displaystyle=\langle
0_{B}|\rho|0_{B}\rangle+\sum_{m^{\prime\prime}_{1},m^{\prime\prime}_{2},m^{\prime\prime}\in
B}\langle
0_{B}|\frac{(\hat{\psi}_{m^{\prime\prime}_{1}\boldsymbol{e}_{1}+m^{\prime\prime}_{2}\boldsymbol{e}_{2},m^{\prime\prime}})^{N_{b}}}{\sqrt{N_{b}!}}\rho\frac{(\hat{\psi}_{m^{\prime\prime}_{1}\boldsymbol{e}_{1}+m^{\prime\prime}_{2}\boldsymbol{e}_{2},m^{\prime\prime}}^{\dagger})^{N_{b}}}{\sqrt{N_{b}!}}|0_{B}\rangle$
$\displaystyle=\frac{1}{3L^{2}}\sum_{\scriptsize\begin{array}[]{l}m_{1},m_{2},m,\\\
m_{1}^{\prime},m_{2}^{\prime},m^{\prime}\end{array}\in
A}e^{i(\frac{2\pi(m-m^{\prime})}{3}\ell-(\alpha_{m}-\alpha_{m^{\prime}}))}e^{\frac{2\pi
i}{L}(k_{1}(m_{1}-m_{1}^{\prime})+k_{2}(m_{2}-m_{2}^{\prime}))}\left(\frac{(\hat{\psi}_{m_{1}\boldsymbol{e}_{1}+m_{2}\boldsymbol{e}_{2},m}^{\dagger})^{N_{b}}}{\sqrt{N_{b}!}}|0_{A}\rangle\langle
0_{A}|\frac{(\hat{\psi}_{m_{1}^{\prime}\boldsymbol{e}_{1}+m_{2}^{\prime}\boldsymbol{e}_{2},m^{\prime}})^{N_{b}}}{\sqrt{N_{b}!}}\right)$
(S31) $\displaystyle\quad+\frac{3L^{2}-3N_{s}}{3L^{2}}|0_{A}\rangle\langle
0_{A}|.$ (S32)

Here, the first term would allow for arbitrary indices
$(m_{1},m_{2},m),(m_{1}^{\prime},m_{1}^{\prime},m^{\prime})\in A$ because they
all correspond to zero particles in region $B$. For the second term, there
must be $(m_{1},m_{2},m)=(m_{1}^{\prime},m_{2}^{\prime},m^{\prime})$ equaling
to the trace indices
$(m_{1}^{\prime\prime},m_{2}^{\prime\prime},m^{\prime\prime})\in B$, so all
the phase factors vanish, and there are $3L^{2}-3N_{s}$ sites in region $B$
giving rise to the prefactor. Denote

$\displaystyle F_{m_{1},m_{2},m}=e^{-i(\frac{2\pi
m}{3}\ell-\alpha_{m})}e^{-\frac{2\pi
i}{L}(k_{1}m_{1}+k_{2}m_{2})},\qquad|F\rangle=(F_{0,0,0},F_{0,0,1},F_{0,0,2},F_{0,1,0},\dots,F_{L_{x}^{(A)},L_{y}^{(A)},2})^{T},$
(S33)

the reduced density operator, written in the matrix form, has

$\displaystyle\rho_{A}$
$\displaystyle=\frac{1}{3L^{2}}\begin{pmatrix}\begin{pmatrix}|F\rangle\langle
F|\end{pmatrix}_{3N_{s}\times 3N_{s}}&0\\\ 0&3L^{2}-3N_{s}\end{pmatrix}$ (S34)

Apparently, there are only 2 nonzero eigenvalues. The first one is for the
$(3N_{s}\times 3N_{s})$ matrix whose eigenvector is
$(1/\sqrt{3N_{s}})|F\rangle$ (note $\langle F|F\rangle=3N_{s}$) corresponding
to the eigenvalue $N_{s}/L^{2}$. The second one is the $1\times 1$ part
obviously corresponding to eigenvalue $(L^{2}-N_{s})/L^{2}$. That yields the
entanglement entropy

$\displaystyle
S_{\text{ent}}=-\text{Tr}\left(\rho_{A}\ln\rho_{A}\right)=-\gamma\ln\gamma-(1-\gamma)\ln(1-\gamma),\qquad\gamma\equiv\frac{N_{s}}{L^{2}}.$
(S35)

Therefore, for the choices of region $A$ in the main text Fig. 2 (b), we have
$S_{\text{ent}}=\ln 2\approx 0.6931$ for $L=2(\gamma=1/2)$, and
$S_{\text{ent}}=(4/9)\ln(4/9)+(5/9)\ln(5/9)\approx 0.6870$ for
$L=3(\gamma=4/9)$. They lead to slight differences of $S_{\text{ent}}$ at
$\lambda\rightarrow 0$ due to different subsystem portions
$\gamma=N_{s}/L^{2}$. However, the reference scar vanishing point
$\lambda_{0}\approx 0.135$ is unlikely to be dominated by such differences, as
significant deviations of $S_{\text{ent}}$ already takes place there compared
with $S_{\text{ent}}$ at $\lambda\rightarrow 0$. This is also confirmed by the
IPR scaling in Fig. 3 (d) of main texts (irrelevant of subsystem size) that
also gives $\lambda_{0}\approx 0.135$.

##  S-3 More details for experimental proposals

This section gives a more detailed account for the experimental proposals.
They are most relevant to the kagome lattice platform at Berkeley [58, 32, 59,
60, 61], while similar schemes can be generalized into other lattices. The
laser system we consider consists of 6 beams, with 3 being red
($\lambda_{R}=1064$nm) and 3 green ($\lambda_{G}=532$nm). They are directed
along the same plane, where each set of monochromatic laser beams form
$120^{\circ}$ angles with respect to each other, as in Fig. S5 (a). The laser
system exhibits good tunability in forming different lattices in a unified
setting, including the (trimerized) kagome, honeycomb, stripe, and Su-
Schrieffer-Heeger types of lattices. We would discuss two exemplary choices of
experimental setup in the following.

###  S-3.1 Scheme I: Trimerized kagome lattice (TKL)

The TKL setting in Ref. [32] uses all 6 bichromatic lasers beams, where green
beams are polarized in-plane, while red ones are along $z$. The lattice
potential is given by

$\displaystyle V(x,y)$ $\displaystyle=V_{R}(x,y)+V_{G}(x,y),$ $\displaystyle
V_{R}(x,y)$ $\displaystyle=-V_{0}^{(R)}\left|\sum_{j=1}^{3}e^{\frac{2\pi
i}{\lambda_{R}}((x-x_{0})\cos\Phi_{j}+(y-y_{0})\sin\Phi_{j})}\right|^{2},$
$\displaystyle V_{G}(x,y)$
$\displaystyle=V_{0}^{(G)}\left(+\left|\sum_{j=1}^{3}\cos(\Phi_{j}+\pi/2)e^{\frac{2\pi
i}{\lambda_{G}}(x\cos\Phi_{j}+y\sin\Phi_{j})}\right|^{2}+\left|\sum_{j=1}^{3}\sin(\Phi_{j}+\pi/2)e^{\frac{2\pi
i}{\lambda_{G}}(x\cos\Phi_{j}+y\sin\Phi_{j})}\right|^{2}\right),$
$\displaystyle(\Phi_{1},\Phi_{2},\Phi_{3})=\left(-\frac{\pi}{2},\frac{\pi}{6},\frac{5\pi}{6}\right).$
(S36)

(a) Lasers

(b) $V_{0}^{(R)}=V_{0}^{(G)}$

(b) $V_{0}^{(R)}=0.8V_{0}^{(G)}$

Figure S5: Schemes for the driven trimerized kagome lattice, where
$x_{0}\approx 405\text{nm},y_{0}=0$. The extent of trimerization can be tuned
continuously by the relative strength between red and green laser beams.

TKLs correspond to our main text model, and the associated driving protocols
are

  * •

Shake the lattice circularly in order to endow a $\pi/2$ flux per triangle.
This shaking is kept on throughout the whole period.

  * •

Drive the lattice potential strength periodically, such that the Hamiltonian
switches between hopping-dominant terms for $\hat{H}_{1}$ and Hubbard
interaction (plus possible sublattice energy offset) dominant terms for
$\hat{H}_{2}$. Then, the Hamiltonians (in ideal situations) take the form as
in Eq. (1) of the main text.

  * •

For specific parameter control, one can fix the duration $t_{1}$ for the first
half of a period according to $\phi_{1}=Jt_{1}/\hbar\approx 2\pi/3\sqrt{3}$,
where $J$ is the hopping strength of strong bonds. Then,
$\lambda=J^{\prime}/J$, with $J^{\prime}$ the hopping strength for weak bonds.
Similarly, the “interaction strength” for the Floquet parameter
$\phi_{2}=Ut_{2}/\hbar$ can be controlled by the duration $t_{2}$ in the
second half of a period, where $U$ is the Hubbard interaction strength. One
driving period $T=t_{1}+t_{2}$ here.

###  S-3.2 Scheme II: Shaken honeycomb lattice

The honeycomb lattice [60] only requires monochromatic lasers, while for our
purposes the lattice should be dimerized as described later. Here, we propose
to use the green laser beams to generate a honeycomb lattice, while the red
beams would be used later to engineer the initial state. So the green beams
here should be polarized along $z$, giving

$\displaystyle
V(x,y)=V_{G}^{\prime}(x,y)=+V_{0}\left|\sum_{j=1}^{3}e^{\frac{2\pi
i}{\lambda_{G}}(x\cos\Phi_{j}+y\sin\Phi_{j})}\right|^{2},\qquad(\Phi_{1},\Phi_{2},\Phi_{3})=\left(-\frac{\pi}{2},\frac{\pi}{6},\frac{5\pi}{6}\right).$
(S37)

The lattice potentials are illustrated in Fig. S6. The corresponding driving
protocols here are

  1. 1.

One can smoothly control the extent of dimerization, shown in Fig. S6 (c) (d),
by shaking linearly the whole lattice at all time [62].

  2. 2.

Similar to the trimerized kagome lattice, we can add a driving in terms of
laser intensity to produce the relatively slow Floquet driving, where the
first and second step of a Floquet driving produce the Hamiltonians

Hopping to three neighbors with unequal strength:
$\displaystyle\frac{\hat{H}_{1}T}{2\hbar}=\phi_{1}\sum_{\boldsymbol{r}}\hat{\psi}^{\dagger}_{\boldsymbol{r}}(\hat{\psi}_{\boldsymbol{r}+\boldsymbol{d}_{1}}+\lambda(\hat{\psi}_{\boldsymbol{r}+\boldsymbol{d}_{2}}+\hat{\psi}_{\boldsymbol{r}+\boldsymbol{d}_{3}}))$
(S38) Onsite interactions and possible sublattice energy offsets:
$\displaystyle\frac{\hat{H}_{2}T}{2\hbar}=\phi_{2}\sum_{\boldsymbol{r}}n_{\boldsymbol{r}}(n_{\boldsymbol{r}}-1)+\theta_{\boldsymbol{r}}n_{\boldsymbol{r}}$
(S39)

(a) Lasers and shaking scheme

(b) Potential contours

(c) Honeycomb lattice with dimerized hopping strength

Three nearest neighbor bonds
$\boldsymbol{d}_{1}=\boldsymbol{e}_{x},\boldsymbol{d}_{2}=-\frac{1}{2}\boldsymbol{e}_{x}+\frac{\sqrt{3}}{2}\boldsymbol{e}_{y},\boldsymbol{d}_{3}=-\frac{1}{2}\boldsymbol{e}_{x}-\frac{\sqrt{3}}{2}\boldsymbol{e}_{y}$

$\gamma_{m}\rightarrow\gamma_{m}J_{0}\left(\frac{m\omega\boldsymbol{d}_{m}\cdot\boldsymbol{e}_{y}}{\hbar}\right)$

Strong bond $1$ not affected, weak bonds $\rightarrow 0$ when shaking
frequency $\frac{m\omega\sqrt{3}}{2\hbar}\rightarrow 2.405$

Figure S6:  The dimerized honeycomb lattice. The extent of dimerization is
achieved by tuning the shaking frequencies. $J_{0}(x)$ is the Bessel function
of the first kind.

###  S-3.3 Initial state preparation

Now, we show the scheme to realize initial states of depositing particles in
one sublattice. That can be achieved by taking advantage of the highly tunable
kagome optical lattice platform, for both the trimerized kagome and dimerized
honeycomb settings, as illustrated in Fig. S7.

(a) Trimerized kagome case

(b) Dimerized honeycomb case

Figure S7: Laser schemes to prepare for initial states. (a) For the trimerized
kagome lattice case, the scheme is realized by slightly modifying the setting
in Fig. S5 by tuning the relative phase of the red beam, such that
$x_{0}\approx 0.5nm$ and the minima of the red laser potential matches one
sublattice site. After loading the atoms to one sublattice as in (a), one
could quench $x_{0}\rightarrow 0.405nm$ and recover the scheme in Fig. S5. (b)
For the dimerized honeycomb case, in addition to the scheme in Fig. S6 using
purely green beams, one could further apply the red beams polarized along $z$
forming triangular lattices. Matching the potential minima to one of the
honeycomb site would produce the plotted potentials for preparing initial
states.

###  S-3.4 Detection of sublattice particle number by band mapping

The projection measurement of particle number in individual sublattice sites
of a non-primitive unit cell of an optical lattice can be performed in three
steps. First, the lattice depth is suddenly increased to quench tunneling
between lattice sites and project the subsystem in each site to approximately
a number Fock state. Then, the superlattice potential is adiabatically
deformed, say by changing the relative position of the two underlying
sublattices that add up to create the trimerized kagome lattice or the
dimerized honeycomb lattice, to energetically detune all subllatice sites. If
the detuning is sufficiently large and all the sublattice sites are decoupled,
then each energy band of the system is predominantly associated with one
sublattice site only. Finally, band mapping, a standard technique where the
lattice potential is addiabatically turned off to map quasimomentum to free-
particle momentum, is performed. Band population can thus be measured in time-
of-flight imaging.

Sublattice-site particle number measurement performed with a similar but
slightly different technique can be found in [63].

###  S-3.5 Simulations of results

To simulate concrete experimental situation with large lattices and filling
fractions, we resort to a semiclassical numerical method, the truncated Wigner
approximation (TWA) [64]. Roughly speaking, this method goes beyond a mean
field analysis by sampling over different initial states

$\displaystyle
W[\Phi_{\boldsymbol{r}\mu}^{(0)}]=\frac{1}{2\pi\sigma^{2}}e^{-|\Phi_{\boldsymbol{r}\mu}^{(0)}-\sqrt{n_{\boldsymbol{r}\mu}^{(0)}}|^{2}/2\sigma^{2}}$
(S40)

where quantum operators $\hat{\psi}_{\boldsymbol{r}\mu}$ are replaced by their
mean field values $\psi_{\boldsymbol{r}\mu}$, and their could deviation from
the initial state values $\\{n_{\boldsymbol{r}\mu}^{(0)}\\}$ in different
samplings. The fluctuation $\sigma=1/2$ means there is on average one half
excessive particles per site $\langle\delta
n_{\boldsymbol{r}\mu}\rangle=\int_{-\infty}^{\infty}|\Phi_{\boldsymbol{r}\mu}^{(0)}-\sqrt{n^{(0)}_{\boldsymbol{r}\mu}}|^{2}W[\Phi_{\boldsymbol{r}\mu}^{(0)}]d(Re\Phi_{\boldsymbol{r}\mu})d(Im\Phi_{\boldsymbol{r}\mu})=2\sigma^{2}=1/2$,
which will be canceled by the symmetrization process for transforming
operators into Weyl symbols, i.e.
$\hat{n}_{\boldsymbol{r}\mu}=\frac{1}{2}(\hat{\psi}_{\boldsymbol{r}\mu}^{\dagger}\hat{\psi}_{\boldsymbol{r}\mu}+\hat{\psi}_{\boldsymbol{r}\mu}\hat{\psi}_{\boldsymbol{r}\mu}^{\dagger})-1/2\rightarrow|\Phi_{\boldsymbol{r}\mu}|^{2}-1/2$
[64]. Meanwhile, the evolutions are still prescribed by classical differential
equations, which can be obtained by first using the Heisenberg’s equation of
motion
$i\partial_{t}\hat{\psi}_{\boldsymbol{r}\mu}=[\hat{\psi}_{\boldsymbol{r}\mu},H]$
and then replacing $\hat{\psi}_{\boldsymbol{r}\mu}$ with Weyl symbols (complex
numbers) $\Phi_{\boldsymbol{r},\mu}$. For instance, the main text model of
trimerized kagome lattice in Eq. (1) prescribes the equations of motion in
each period $T$ as

$\displaystyle t\in[0,T/2):\qquad$
$\displaystyle\qquad\partial_{t}\Phi_{\boldsymbol{r}\mu}=J\sum_{\nu\neq\mu}f_{\mu\nu}(\Phi_{\boldsymbol{r}\nu}+\lambda\Phi_{\boldsymbol{r}-\boldsymbol{e}_{\nu}+\boldsymbol{e}_{\mu},\nu}),$
(S41) $\displaystyle t\in[T/2,T):\qquad$ $\displaystyle\qquad
i\partial_{t}\Phi_{\boldsymbol{r}\mu}=2U|\Phi_{\boldsymbol{r}\mu}|^{2}\Phi_{\boldsymbol{r}\mu}$
(S42)

Here the parameter $J,U$ are related to those in Eq. (1) by
$JT/2\hbar=\phi_{1},UT/2\hbar=\phi_{2}$. This way, the semiclassical results
are expected to capture the quantum fluctuations during early time of
evolution, which is most relevant to experimental observations. As a side
remark, we notice that if only a pure mean field simulation is adopted (no
initial state sampling), one would observe a deceptive infinite time DTC
oscillation without decay at all for a rather wide range of initial states.
That contradicts exact diagonalization and analytical results, and the generic
thermalizing nature confirmed by level spacing statistics. Therefore, it is of
vital importance to incorporate fluctuations at least for the initial states
so as to simulate a realistic situation at early time.

(a) Numerical results

(experimentally accessible) (b) No separation

(c) Separation by 1 cell

(d) Separation by 2 cells

Figure S8: TWA simulation of scar enforced DTC dynamics in trimerized kagome
lattice. We simulate a lattice containing $12\times 12$ unit cells under
periodic boundary conditions. Each populated site denoted by colored dots in
(b) – (d) contains 5 bosons as the initial state. Evolutions from these states
are simulated in (a) for the corresponding colors. Two sublattices are denoted
as $\mu=0$ and $\mu=1$ following the convention in main text. $n_{0}$ is the
total particle number in sublattice $0$, and $N_{b}$ the total number of
bosons in all sites. Parameters for the model written in Eq. (S38) and (S39)
are $\phi_{1}=\pi/2,\phi_{2}=1.1,\theta_{\boldsymbol{r}}=0,\lambda=0.05$.
Simulations contain 3000 initial state Monte-Carlo samples for each site and
error bars denote the standard deviation when all data are group into 10 bins.

(a) Numerical results

(b) No separation

(experimentally accessible) (c) Separation by 1 cell

(d) Separation by 2 cells

Figure S9:  Similar simulations as in Fig. S8 for the dimerized honeycomb
lattice. System size is again $12\times 12$ unit cells, and each initially
populated site in (b) – (d) contains 5 bosons. $n_{0}$ means the total
particle number in sublattice $\mu=0$ and $N_{b}$ is the total boson number.
Parameters for the model written in Eq. (S38) and (S39) are
$\phi_{1}=\pi/2,\phi_{2}=1.1,\theta_{\boldsymbol{r}}=0,\lambda=0.05$.
Simulations contain 3000 initial state Monte-Carlo samples for each site and
error bars denote the standard deviation when all data are group into 10 bins.

The results for sublattice density dynamics is shown in Fig. S8 and S9, where
we compare the DTC decay rates starting from different initial states. It is
clear that with larger spatial separations for the initially populated sites
lead to prolonged oscillations due to longer time needed for FBS localized at
different unit cells to interact with each other.