Systematic construction of scarred many-body dynamics in 1D lattice models

TL;DR

This paper develops a systematic method to construct 1D lattice models exhibiting many-body scarring, leading to robust wave-function revivals and slow thermalization, expanding understanding of non-integrable quantum dynamics.

Contribution

It introduces a systematic construction of non-integrable 1D models with many-body scarring, generalizing phenomena observed in Rydberg atom systems.

Findings

Models exhibit robust wave-function revivals.

Scarring dynamics involve free precession and an interacting bottleneck.

Models include kinetically-constrained quantum clock models.

Abstract

We introduce a family of non-integrable 1D lattice models that feature robust periodic revivals under a global quench from certain initial product states, thus generalizing the phenomenon of many-body scarring recently observed in Rydberg atom quantum simulators. Our construction is based on a systematic embedding of the single-site unitary dynamics into a kinetically-constrained many-body system. We numerically demonstrate that this construction yields new families of models with robust wave-function revivals, and it includes kinetically-constrained quantum clock models as a special case. We show that scarring dynamics in these models can be decomposed into a period of nearly free clock precession and an interacting bottleneck, shedding light on their anomalously slow thermalization when quenched from special initial states.

Figures

Equations

$$H=\sum _j{P}_{j-1}^0{X}_j{P}_{j+1}^0,P_j^0\equiv |0_j⟩⟨0_j|,$$

$$h=i\sum _{n=0}^{N_c-1}\frac{2\pi i}{T}{k}_n|{\psi }_n⟩⟨{\psi }_n|.$$

$$H=\sum _{j=0}^{N-1}{P}_{j-1}{h}_j{P}_{j+1},P_j\equiv |{\varphi }_j⟩⟨{\varphi }_j|,$$

$$\mathcal{U}=e^{-iC}=\sum _{j=0}^{N_c-1}|j+1⟩⟨j|.$$

$$H_{\mathrm{clock}}=\sum _j{P}_{j-1}^0{C}_j{P}_{j+1}^0,$$

$${\mathcal{C}}_i=\left(\begin{array}{ccccc}1& 0& 0& ...& 0\\ 0& 0& ...& \pm 1& 0\\ 0& ...& \pm 1& 0& 0\\ ...& & & & \\ 0& \pm 1& 0& ...& 0\end{array}\right),$$

$$I:$$

$$T:$$

$$[H,I]$$

$$[H,\mathcal{C}T]$$

$$\{H,\mathcal{C}\}$$

$$\{H,\mathcal{C}I\}$$

$$\mathcal{U}=e^{-iC}=\sum _{j=0}^{N_c-1}|j+1⟩⟨j|,$$

$$H_{clock}$$

$$=\sum _n{P}_{n-1}^{{}^{\prime }}{X}_n{P}_{n+1}^{{}^{\prime }},P^{{}^{\prime }}=\left(\begin{array}{ccc}1/6& 1/3& -1/6\\ 1/3& 2/3& -1/3\\ -1/6& -1/3& 1/6\end{array}\right).$$

$$H=f\sum _{j=1}^N{\tau }_j^{†}{e}^{-i\varphi }+J\sum _{j=1}^{N-1}{\sigma }_j^{†}{\sigma }_{j+1}{e}^{-i\theta }+h.c,$$

$${\sigma }_j=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& \omega & 0& 0\\ 0& 0& {\omega }^2& 0\\ 0& 0& 0& {\omega }^3\end{array}\right),{\tau }_j=\left(\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\\ 1& 0& 0& 0\end{array}\right),$$

$$H^{site}=\left(\begin{array}{cccc}0& -i& 0& i\\ i& 0& -i& 0\\ 0& i& 0& -i\\ -i& 0& i& 0\end{array}\right).$$

$$H_{chiral}=\sum _j{P}_{j-1}{H}_j^{site}{P}_{j+1},$$

$$C_n^{(A)}$$

$$C_n^{(B)}$$

$$E^4-a{E}^2+b=0,$$

$$\text{(A)}a$$

$$b$$

$$\text{(B)}a$$

$$b$$

$$E_{1,\dots ,4}=\mp \frac{1}{\sqrt{2}}\sqrt{a\pm \sqrt{a^2-4b}}.$$

$$9a^2=100b,\Delta E_2=\Delta E_1,$$

$$4a^2=25b,\Delta E_2=2\Delta E_1,$$

$$25^2=676b,\Delta E_2=\frac{1}{2}\Delta E_1,$$

$$a^2=4b,\Delta E_1=\Delta E_3=0,$$

$$b=0,\Delta E_2=0$$

$$H_d^{PPXP}=\sum _n{P}_{n-1}{X}_n{P}_{n+1}(P_{n-d}+P_{n+d})$$

$$H_d^{PXXXP}=\sum _n{P}_{n-1}{X}_n{X}_{n+1},...,X_{n+d-1}{P}_{n+d},$$

Full text

Introduction.—The understanding of ergodicity and thermalization in isolated quantum systems is an open problem in many-body physics, with important implications for a variety of experimental systems Kinoshita et al. (2006); Schreiber et al. (2015); Smith et al. (2016); Kucsko et al. (2016); Kaufman et al. (2016). On the one hand, this problem has inspired important developments such as Eigenstate Thermalization Hypothesis (ETH) Deutsch (1991); Srednicki (1994); Rigol et al. (2008), which establishes a link between ergodicity and the properties of the system’s eigenstates. On the other hand, strong violation of ergodicity can result in rich new physics, such as in integrable systems Sutherland (2004), Anderson insulators Anderson (1958), and many-body localized phases Basko et al. (2006); Serbyn et al. (2013); Huse et al. (2014). In these cases, the emergence of many conservation laws prevents the system, initialized in a random state, from fully exploring all allowed configurations in the Hilbert space, causing a strong ergodicity breaking.

A recent experiment on an interacting quantum simulator Bernien et al. (2017) has reported a surprising observation of quantum dynamics that is suggestive of weak ergodicity breaking. Utilizing large 1D chains of Rydberg atoms Schauß et al. (2012); Labuhn et al. (2016); Bernien et al. (2017), the experiment showed that quenching the system from a Néel initial state lead to persistent revivals of local observables, while other initial states exhibited fast equilibration without any revivals. The stark sensitivity of the system’s dynamics to the initial states, which were all effectively drawn from an “infinite temperature” ensemble, appeared at odds with “strong” ETH D’Alessio et al. (2016); Gogolin and Eisert (2016); Mori et al. (2018).

In Ref. Turner et al., 2018a; Ho et al., 2019 the non-ergodic dynamics of a Rydberg atom chain was interpreted as a many-body generalization of the classic phenomenon of quantum scar Heller (1984). For a quantum particle in a stadium billiard, scars represent an anomalous concentration of the particle’s trajectory around (unstable) periodic orbits in the corresponding classical system, which has an impact on optical and transport properties Sridhar (1991); Marcus et al. (1992); Wilkinson et al. (1996). By contrast, in the strongly interacting Rydberg atom chain initialized in the Néel state, quantum dynamics remains concentrated around a small subset of states in the many-body Hilbert space, thus it is effectively “semiclassical” Ho et al. (2019). While recent works Khemani et al. (2019); Choi et al. (2018) have shown that revivals can be significantly enhanced by certain perturbations to the system, a general understanding of the conditions that allow scars to occur in a many-body quantum system is still lacking.

The observation of periodic dynamics was linked to the existence of atypical eigenstates at evenly spaced energies throughout the spectrum of the many-body system Turner et al. (2018a, b); Iadecola et al. (2019). Similar phenomenology was previously found in the non-integrable AKLT model Moudgalya et al. (2018a, b), where highly-excited eigenstates with low entanglement have been analytically constructed. A few of such exact eigenstates are now also available for the Rydberg atom chain model Lin and Motrunich (2019). In a related development, it was proposed that atypical eigenstates of one Hamiltonian can be “embedded” into the spectrum of another, ETH-violating, Hamiltonian Shiraishi and Mori (2017). However, although the collection of models that feature atypical eigenstates is rapidly expanding Kormos et al. (2016); James et al. (2019); Robinson et al. (2019); Iadecola and Znidaric (2018); Surace et al. (2019), including recent examples of topological phases Ok et al. (2019) and fractons Pai et al. (2019), their relation to periodic dynamics remains largely unclear.

In this paper we systematically construct interacting lattice models that exhibit periodic quantum revivals when quenched from a particular product state. The basic building block has a Hilbert space containing $N_{c}$ states (“colors”) and a time-independent Hamiltonian that yields periodic unitary dynamics, $\mathcal{U}(t+T)=\mathcal{U}(t)$. The interacting models are defined by coupling these building blocks under a kinetic constraint. Intriguingly, the dynamics in these models can be decomposed into periods of nearly free precession, in which the local degrees of freedom coherently cycle through the available states on a single site, followed by an interacting segment of dynamical evolution, reminiscent of a kicked quantum top Haake (2001). In all cases, the existence of atypical scarred eigenstates underpins the revivals. We show that our construction includes known models, such as chiral clock models Fendley (2014), which are shown to support scars, and also gives a way of enhancing the revivals in spin-$s$ generalisations of the Rydberg chain Ho et al. (2019). Moreover, in selected cases with small values of $N_{c}$, we numerically explore general deformations of the models, verifying that our construction yields optimal models, i.e., with the highest amplitude of the wave function revivals.

PXP model.—We start by briefly reviewing the model of a 1D Rydberg atom chain Sun and Robicheaux (2008); Olmos et al. (2009, 2012); Lesanovsky and Katsura (2012). The system can be modelled as coupled two level systems (with states $|0\rangle$, $|1\rangle$) described by an effective “PXP” Hamiltonian

[TABLE]

where $X_{j}=|0_{j}\rangle\langle 1_{j}|+|1_{j}\rangle\langle 0_{j}|$ denotes the Pauli matrix. The model in Eq. (1) describes a kinetically constrained paramagnet Ostmann et al. (2018): each atom can flip only if both its neighbors are in $|0\rangle$ state.

The Hamiltonian in Eq. (1) is interacting and non-integrable Turner et al. (2018a), yet it exhibits unconventional thermalization. For example, the model has atypical (ETH-violating) eigenstates with low entanglement at high energy densities Turner et al. (2018b). Moreover, when the system is quenched from the Néel initial state, $|\psi_{0}\rangle=|0101\ldots\rangle$, local observables such as domain wall density Bernien et al. (2017) and even the many-body wave function fidelity, $F(t)=|\langle\psi_{0}|\psi(t)\rangle|^{2}$, all revive with the same frequency Turner et al. (2018a); Schecter and Iadecola (2018); Surace et al. (2019). At the same time, quenches from other initial states, such as $|0000\ldots\rangle$, do not lead to observable revivals Bernien et al. (2017). The revival frequency from the Néel state is set by the energy separation between atypical eigenstates, as the same eigenstates also maximize the overlap with the Néel state Turner et al. (2018a). Thus, the quench dynamics from the Néel state is largely restricted to a subset of many-body eigenstates, and can be viewed as precession of a large spin, which traces a periodic orbit that can be accurately captured by time-dependent variational principle (TDVP) on a manifold spanned by matrix product states with low bond dimension Ho et al. (2019).

Construction of scarred models.— Consider now a system with a local basis $|0\rangle$, $|1\rangle$, …, $|N_{c}-1\rangle$, and an arbitrary time independent Hamiltonian $h$ whose unitary dynamics is periodic, such that $\mathcal{U}^{T}\equiv\exp(-ihT)=\mathds{I}$ for arbitrary $T$ (not necessarily integer). The eigenvalues of $\mathcal{U}$ are then $\lambda_{n}=\exp(i2\pi k_{n}/T)$, with the corresponding eigenvectors $|\psi_{n}\rangle$, where $k_{n}$ are arbitrary integers. We can then obtain candidate Hamiltonians $h$ by choosing particular $\{\lambda_{n}\}$ which guarantee a periodic $\mathcal{U}$ and taking its logarithm:

[TABLE]

To define a many-body lattice Hamiltonian, we take a tensor product of $h$ and impose the kinetic constraint that $h$ only acts on sites whose neighbors are in some unlocking state $|\phi\rangle$:

[TABLE]

where $N$ is the number of lattice sites. The only other condition we place on $h$ is that the many-body system possesses a particle-hole symmetry $\rho$, which anticommutes with $H$, $\{H,\rho\}=0$, leading to the symmetry $E\leftrightarrow-E$ of the energy spectrum. This is motivated by the fact that PXP model possesses such a symmetry, and its revivals are improved by the addition of perturbations which preserve this symmetry Khemani et al. (2019); Choi et al. (2018). Precise form of $\rho$ is unimportant here and can be found in SOM . We thus focus on cases where $\{k_{n}\}$ are symmetric around zero, resulting in $h$ being off diagonal and compatible with $\rho$. A particularly illustrative example of this construction is when $\mathcal{U}$ is interpreted as the shift operator of a quantum clock Fendley (2014); Vernier et al. (2018), as we explain next.

Scars in clock models.—The scarred clock models are defined by choosing $T=N_{c}$, which gives

[TABLE]

In this case, $\lambda_{n}=\exp(2\pi ik_{n}/N_{c})$ and $|\psi_{n}\rangle=\sum_{j=0}^{N_{c}-1}(1/\lambda_{n}^{j})|j\rangle$. For odd $N_{c}$, $k_{n}$ takes the values $-\frac{N_{c}-1}{2},\ldots,0,\ldots,\frac{N_{c}-1}{2}$. For $N_{c}$-even, we need to double the period, $T=2N_{c}$, in order to make $h$ off-diagonal in the $|j\rangle$ basis. This allows to choose $k=-\frac{N_{c}-1}{2},\ldots,-\frac{1}{2},\frac{1}{2},\ldots,\frac{N_{c}-1}{2}$, and Eq. (S8) continues to be valid for $N_{c}$-even after performing a gauge transformation, $|j\rangle\to e^{i\pi j/N_{c}}|j\rangle$.

The inspiration behind Eq. (S8) is that local dynamics is a cyclic rotation around the basis of $N_{c}$ “clock” states $|j\rangle$, Fig. 1(a). The many-clock “PCP” Hamiltonian from Eq. (3) then becomes

[TABLE]

with $C$ represents $h$ in Eq. (2). Without loss of generality, we can choose the projector onto any of the clock basis states, e.g., $P^{0}=|0\rangle\langle 0|$. Thus, each site precesses around the clock of basis states if both its neighbors are in $|0\rangle$ state, otherwise the site remains frozen, Fig 1(a).

We have studied the PCP model in Eq. (5) using exact diagonalization with periodic boundary conditions. For any $N_{c}\leq 12$ accessible to us numerically, we find long-lived oscillatory dynamics when the system is quenched from any Néel-like state, $|0101...\rangle$, $|0202...\rangle$, etc. Fig. 1(b) summarizes the result for $N_{c}=4$. The dynamics proceeds in two steps. First, each unfrozen clock nearly freely cycles through its states, $|1\rangle\to|2\rangle\to\ldots|N_{c}-1\rangle$. After this coherent process is complete, the many-clock state shifts, $|N_{c}-1,0,N_{c}-1,0\ldots\rangle\to|0101\ldots\rangle$. In this second step, interactions kick in and some fidelity is lost to thermalization. This cycle of free precession followed by a short interacting segment is reminiscent of a kicked quantum top. We now see that the PXP model (which is equivalent to $N_{c}=2$ clock) is special in that it lacks free-precession dynamics. On the other hand, similar to the PXP case, in scarred clock models coherence also remains protected to a large degree during the interacting part of the process, allowing the wave function to keep returning to the initial state at later times.

In order to visualize the dynamics, in Fig. 1(b) we plot the fidelity $|\langle\phi|\exp(-itH)|1010\ldots\rangle|^{2}$ w.r.t. several product states $|\phi\rangle$ corresponding to either the initial state, the internal shift of each clock, or to the overall translation of the initial state. The duration of individual clock ticks (e.g., $|1010\ldots\rangle\to|2020\ldots\rangle$) matches that of the unconstrained clock model. Following the convention that $C$ is rescaled such that nearest neighbor hoppings have magnitude one, the frequency of the putative free precession is found to be $\approx 0.902$ (in units $\hbar=1$) while the frequency of the single site precession (in the absence of a constraint) is $\approx 0.900$. We note that time evolution of local observables is consistent with the presented picture of the underlying dynamics SOM .

In Fig. 1(c) we show the overlap of all eigenstates with the Néel state $|0101\ldots\rangle$, while Fig. 1(d) shows the bipartite entanglement entropy $S=-\mathrm{tr}\rho_{A}\ln\rho_{A}$, where $\rho_{A}$ is the reduced density matrix of one half of the chain. The scar states are easily identifiable as a band of special eigenstates (circled in red) that extend throughout the spectrum. Total number of special states is $(N_{c}-1)N+1$. Similar to the PXP model, the special eigenstates are distinguished by their high overlap with the Néel state, or alternatively as ones with atypically low entanglement. We note that some of the eigenstates with smallest entanglement belong to a different band of scarred states associated with a “defected $\mathds{Z}_{4}$” state $|20002030103000\rangle$ (blue circles in Fig. 1(d)). Apart from these special states, we also observe tower structures in the spectrum, which reflect the clustering of neighboring eigenstates around the energies of the scarred eigenstates. Deep in the bulk of the spectrum, the density of states [indicated by color scheme in Fig. 1(c)] appears to be uniform, as expected from the ETH. Indeed, at $N=14$ we find a mean level spacing ratio Oganesyan and Huse (2007) of $\langle r\rangle=0.5218$, consistent with Wigner-Dyson statistics and the ETH. We have confirmed that the frequency of the revival to the initial state matches the energy separation between special eigenstates in Fig. 1(c).

Relation to spin-$s$ and chiral clock models.—In Ref. Ho et al. (2019) the TDVP approach was generalized to spin-$s$ PXP models with the kinetic constraint $P^{0}$. Periodic revivals were numerically demonstrated for $s=1,2$. Both spin-$s$ PXP model and $N_{c}=2s+1$ colored PCP clock models are obtained from our construction in Eq. (3) by taking $k=-s,...,s$. Thus by performing a basis rotation, the clock Hamiltonian can be expresed in the spin basis, $H_{\rm clock}=\sum_{j}P_{j-1}^{\prime}X_{j}P_{j+1}^{\prime}$, where $P^{\prime}$ is a deformation of the projector $P^{0}$ in Eq. (5) SOM . We have numerically found that the number of scarred states remains the same for PXP models expressed in terms of either the spin $P^{0}$ or $P^{\prime}$; however, the amplitude of the revivals is always higher when using $P^{\prime}$ instead of spin $P^{0}$ for $N_{c}$-odd SOM . Thus, our construction shows how to improve the revivals in the standard PXP models. In addition, mapping to the clock representation allows to clearly delineate nearly-free precession from the interacting part of the dynamics, which is not transparent in the spin representation.

Furthermore, our construction includes models for which $C$ is not related to spin matrices via a change of basis. One family of models for even $N_{c}$ is obtained by choosing $k=-\frac{N_{c}}{2},\ldots,-1,1,\ldots,\frac{N_{c}}{2}$, with the projector $P^{0}$ as above. For $N_{c}=4$, the resulting $C$ is equivalent to the 4-color Chiral Clock Model (CCM) at the fixed point in the disordered phase Baxter (1989); Fendley (2014); SOM . This model exhibits two types of oscillatory behaviour. Quenches from $|0202...\rangle$ result in decaying fidelity revivals. Quenches from $|1010...\rangle$, $|3030...\rangle$ essentially freeze out the sublattice containing [math]’s, such that the system oscillates like a nearly free paramagnet SOM .

General phase diagram of scarred models.—We now perform an extensive search for scarred models with the fixed kinetic constraint $P^{0}$. By varying elements of $C$, we scan all models of the form Eq. (5). We map out the phase diagram of these models based on the quality of scars, i.e., the first revival maximum of the fidelity from the Néel-like states. We restrict the matrix $C$ to be purely imaginary and off diagonal, as this leads to a symmetric spetrum SOM , thus preserving the desired particle-hole symmetry.

Consider the $N_{c}=4$ case. Allowed distortions involve varying 5 matrix elements in $C$, so we take slices where only two parameters are simultaneously varied. We consider two cases, (a) vary the next-nearest-neighbor hoppings $C_{02}=C_{13}=\alpha i$, while also varying $C_{03}=-\beta i$, or (b) switch off next-nearest-neighbor hoppings, while varying $C_{12}=-\alpha i$ and $C_{03}=-\beta i$. The corresponding phase diagrams are shown in Fig. 2. On these diagrams there are several limiting cases at special values of $(\beta,\alpha)$. For variation (a), we have: (i) $(1,1/\sqrt{2})$ is $N_{c}=4$ clock; (ii) $(-1,0)$ is $N_{c}=4$ CCM model. For variation (b): (iii) $(0,2/\sqrt{3})$ is spin-$\frac{3}{2}$ PXP; (iv) $(1,-1)$ is also $N_{c}=4$ CCM; (v) at $(0,0)$, we have $C=i\sum_{j=0,2}|j\rangle\langle j+1|-\mathrm{h.c.}$, which (with $P^{0}$) can be viewed as the sum of a spin-$\frac{1}{2}$ PXP and a free $s=\frac{1}{2}$ paramagnet. Points marked F correspond to decoupled free paramagnets.

The maximum fidelity at first revival for $N_{c}$-even is generally comparable between clock and spin-$s$ PXP models. For example, for $N_{c}=4$ in Fig. 2, $F_{max}\approx 0.761$ (clock) and $F_{max}\approx 0.783$ for spin-$\frac{3}{2}$ PXP. For $N_{c}=6$ and $N=8$, we obtain $F_{max}\approx 0.813$ (spin) and $F_{max}\approx 0.802$ (clock), while for $N_{c}=8$, $N=8$ we find $F_{max}\approx 0.793$ (spin) and $F_{max}\approx 0.806$ (clock). On the other hand, for $N_{c}$-odd, we find a considerable improvement in the fidelity of a clock compared to the spin-$s$ PXP model. For example, for $N_{c}=3$, the maximum fidelity of the clock model is $F_{max}\approx 0.724$ versus $F_{max}\approx 0.653$ for spin-$1$; for $N_{c}=5$, $N=10$, the improvement is even bigger, $F_{max}\approx 0.563$ vs. $F_{max}\approx 0.766$ (clock) SOM . Thus, our construction for odd $N_{c}$ gives a way to improve the revivals over corresponding $s=(N_{c}-1)/2$ spin-$s$ PXP models.

Since the phase diagram in Fig. 2 is evidently quite rich, we look for a simple guiding principle that predicts the most robust scarring models. The commensurability of the eigenvalue spectrum of $C$ provides such a criterion, see lines and dots in Fig. 2. White lines mark the models for which $C$ has equidistant energy levels, i.e., its eigenvalues are $E_{n}=k\epsilon$, $k\in\mathds{Z}$. Our $N_{c}=4$ clock model lies on one of these lines, as shown in Fig. 2(a). We can consider further commensurability conditions where the energy spacings of $C$ are in simple ratios, e.g., 1:2 (pink lines). Finally, red points mark the areas where $C$ contains one pair of degenerate eigenvalues (with trivially commensurate spacing in our case). One of these points hosts the $N_{c}=4$ CCM at its fixed point in the disordered phase. Another one, along the diagonal in Fig. 2(b), hosts a combination of the free paramagnet and spin-$\frac{1}{2}$ PXP model. We note, however, that our simple criterion based on the non-interacting spectrum of $C$ only serves as a rough indicator of scarring models. The precise parameter values where such models are realized are determined by the non-trivial interplay between this condition and the kinetic constraint, i.e., $P^{0}$.

Conclusion.—We have presented a systematic construction of non-integrable PCP models exhibiting many-body revivals and quantum scars. The construction is based on embedding local unitary precession, $U^{T}=e^{-iCT}=\mathbb{I}$, into an interacting quantum system. The obtained models are expressed in terms of kinetic constraints which are present in quantum simulators in the Rydberg blockade regime Bernien et al. (2017); Surace et al. (2019); Keesling et al. (2019). Kinetic constraints of this kind also emerge naturally in lattice gauge theories, which have recently been realized in periodically driven optical lattices Schweizer et al. (2019). The strongest reviving models are predicted by considering the commensurability of $C$’s eigenvalues. For odd $N_{c}$ and equidistant eigenvalues for $C$, the obtained models revive better than the corresponding spin $s=(N_{c}-1)/2$ PXP model. Rotating $C\rightarrow X$, $P\rightarrow P^{{}^{\prime}}$, our construction thus provides a prescription for improving PXP revivals. If we do not restrict to equidistant eigenvalues of $C$, our construction yields further families of scarred models not related to PXP by rotation. Further, clock models provide a simple physical picture of the underlying dynamics – a period of nearly free precession followed by an interacting bottleneck. This “effective drive” is reminiscent of kicked systems, where mixed phase space dynamics (both recurrent and thermalizing behaviour), can emerge due to the presence of a continuous spectrum in the Floquet operator Milek and Seba (1990). Taking the same constraint $U^{T}=\mathbb{I}$, one can also engineer time-translation symmetry breaking in driven systems Schaefer et al. (2019); Else et al. (2017). These observations suggest a deeper connection between oscillatory scarred models and time crystals, complementing recent description of scarred PXP states as $\pi$ magnon condensates which possess long range order in both space and time Iadecola et al. (2019).

SI Revivals of local observables

In the main text we analyzed the dynamics in the clock models from the perspective of fidelity revivals, i.e., the many-body wave function returning to a particular product state after time $t$. By inspecting different product states, we have concluded that the dynamics proceeds in two steps: (i) precession of every clock between the available basis states; (ii) transition to a product state that is related to the initial state by an overall translation of the system by one lattice site, e.g., $|1010\ldots\rangle\to|0101\ldots\rangle$. While the process (i) is nearly free (and hence coherent), remarkably even the process (ii) results in a relatively small loss of fidelity.

Revival of a many-body wave function, of course, is formidable to measure experimentally, and here we demonstrate that the same picture of the dynamics can be inferred from studying revivals in local observables, $Q$. For example, measuring a color on a given site, $Q\equiv P^{n}=|n\rangle\langle n|$ will show the same pattern of dynamics, see Fig. S1. Consider the two sublattices, corresponding to even and odd sites. When the system is quenched from the Néel state $|0101...\rangle$, even sites initially cycle through the basis set $1\rightarrow 2\rightarrow 3...$ while odd sites are locked at [math] by the constraint. These sublattices exchange and the process repeats after a full precession $1\rightarrow 0$ has been completed.

SII Discreet Symmetries in Clock Models

In the main text we noted our construction of the clock models yielded systems possessing a desired (anti-unitary) particle-hole symmetry. This was motivated by the fact that the original PXP model possesses a particle-hole symmetry $\rho=\prod_{i}Z_{i}$, which anticommutes with H, $\{H,\rho\}=0$. Our clock model actually possesses several discrete unitary and anti-unitary symmetries, originating from charge conjugation, parity and time reversal. We define single-site charge conjugation as the following $N_{c}\times N_{c}$ operator Whitsitt et al. (2018):

[TABLE]

where $+$ corresponds to odd $N_{c}$ while $-$ corresponds to even $N_{c}$. Further, we define parity and time reversal operators in the conventional way:

[TABLE]

where $\hat{K}$ corresponds to complex conjugation. It then follows that the PCP models possess the following discreet symmetries:

[TABLE]

where $\mathcal{C}=\prod_{i}\mathcal{C}_{i}$.

SIII Relation between clock and spin-$s$ PXP models

In the main text we pointed out that spin-$s$ PXP models are included in our construction as a special case. From the definition of $\mathcal{U}$ governing the free precession of our clock models, which we reproduce here,

[TABLE]

the eigenvalues are found to be $\lambda_{n}=\exp(2\pi ik_{n}/N_{c})$, with $k_{n}\in\{-(N_{c}-1)/2,\ldots,0,\ldots,(N_{c}-1)/2\}$ for odd $N_{c}$, while $k_{n}\in\{-(N_{c}-1)/2,\ldots,-1/2,1/2,\ldots,(N_{c}-1)/2\}$ for even $N_{c}$. This then implies the eigenvalues of $C\propto i\ln(\mathcal{U})$ are proportional to $k_{n}$, which we immediately recognize as angular momenta quantum numbers for spin $2s+1=N_{c}$. This tells us that operators $C$ and $X$ are equivalent by a change of basis.

The distinction between the clock and spin models thus lies entirely in the kinetic constraint, i.e., the choice of the projector $P=|\phi\rangle\langle\phi\rangle$. For example, at $N_{c}=3$, we find:

[TABLE]

For $N_{c}=3$, it is straightforward to visualize the entire phase diagram with the fixed projector $P^{0}$. There are two couplings to vary, as by a rescaling of $C$ we are free to fix the value of one of the three couplings. The resulting phase diagram for $N_{c}=3$ is shown in Fig. S2 as a function of $C_{02}$ and $C_{12}$ matrix elements.

All the regions in the diagram that display strong scarring can be interpreted as deformations of one of our previously discussed models. Specifically, in Fig. S2 we identify the representative points as: (i) and ii) are spin-$\frac{1}{2}$ PXP, (iii) is spin$-1$ PXP, with $P=|S_{z}=-1\rangle\langle S_{z}=-1|$, (iv) spin-$1$ PXP, with $P=|S_{z}=0\rangle\langle S_{z}=0|$, (v) spin-$1$ PXP, with $P=|S_{z}=1\rangle\langle S_{z}=1|$, (vi) $N_{c}=3$-color clock model, (vii) free paramagnet. From this we conclude that for the on-site Hilbert space dimension $N_{c}=3$, the only scarred models are obtained by various ways one can constrain a free paramagnet. As mentioned in the main text, the maximum fidelity of the clock model, $0.724$, is greater than the spin-$1$ model, $0.653$. This improvement in the revival fidelity is seen in all $N_{c}$-odd models. For example, at $N_{c}=5$, $N=10$, spin: $0.563$, clock: $0.766$. Thus, by constraining the free paramagnet with $P^{{}^{\prime}}$ as opposed to $P^{0}$ leads to better revivals for $N_{c}$ odd. Finally, although still expressable as spin models, the clock basis provides a much simpler interpretation of the dynamics, clearly showing a period of free precession followed by an interacting segment.

In Fig. S3 we compare the 5-color clock with the PXP spin-2 model and the 6-color clock with the PXP spin $5/2$ model. The different choice of the projector ($P^{0}$ or $P^{{}^{\prime}}$) leads to a rather different structure of the scar states near the edge of the spectrum. While the corresponding clock and spin models both contain the same number of scar states, we see the band of scars remains better separated from the bulk near the edge of the spectrum in the clock models. Fidelity decays similarly in both models, but the spin model lacks the period of free precession observed in clock models in the main text, instead just bouncing back and forth between Néel and anti-Néel states, $|0404...\rangle$ and $|4040...\rangle$, as in spin-$1/2$ PXP. Clock models were found to “shed” less fidelity for odd $N_{c}$, where the difference in the first revival peak maximum of $|\langle\psi(t)|\psi(0)\rangle|^{2}$ was found to be 0.226 for $N_{c}=5$. For even $N_{c}$, however, the maximum fidelity is comparable between clock and spin models.

SIV Relation to chiral clock models

By analogy with the quantum Ising model in a transverse field, Fendley has studied generalizations of models that involve clock operators Fendley (2014), originally introduced by Baxter Baxter (1989). Instead of focusing on critical properties of such models, we mainly focus on their fixed point, where only the shift operator is present in the Hamiltonian. The first non-trivial case that is distinct from our clock operators occurs for $N_{c}=4$ Chiral Clock Model (CCM) Fendley (2014), with the Hamiltonian defined as:

[TABLE]

where the operators $\sigma$ and $\tau$ generalize the spin-1/2 Pauli matrices, $\sigma_{z}$ and $\sigma_{x}$, respectively. They are given by

[TABLE]

where (for $N_{c}=4$) we have $\omega=i$. Taking $J=0$, $\phi=-\pi/2$, the Hamiltonian (S11) becomes non-interacting, with the on-site Hamiltonian reducing to:

[TABLE]

We construct an interacting many-clock model by introducing the same Rydberg-like constraint:

[TABLE]

with $P_{j}=|0_{j}\rangle\langle 0_{j}|$. This Hamiltonian can be obtained as the deformation of our clock model introduced in the main text.

In Fig. S4 we plot the overlap of all eigenstates of the Hamiltonian in Eq. (S14) with two product states states: $|1010\ldots\rangle$ (top) and $|2020\ldots\rangle$ (bottom). System size is $N=14$ with periodic boundary conditions. These plots show that the model is not completely ergodic and supports a band of scarred states. However, the scarred bands look rather different in the two cases, with a stronger clustering of eigenstates into towers for the case of $|1010\ldots\rangle$ state. We note that the structure of the scarred band in the bottom panel of Fig. S4 is quite reminiscent of the PXP spin-1/2 model Turner et al. (2018b). In fact, as we show next, the dynamics in this case reduces to the same type of oscillation that was found in the PXP spin-1/2 model.

Fig. S5 shows the single site precession of $H^{site}$ in the absence of a constraint. We see that the dominant dynamics is $0\rightarrow 2$ flips and $1\rightarrow 3$ flips. Thus, when the constrained model, dressed with $P^{0}$, is quenched from $|0101...\rangle$, the dominant dynamics is just to oscillate between $|0101...\rangle\rightarrow|0303...\rangle$, as $0\rightarrow 2$ is blocked due to the constraint. Neglecting leakage, this is just a decoupled free paramagnet precessing at every other site. However, when quenched from $|0202...\rangle$, a transition through the polarized state $|0000...\rangle$ is permitted and we see transitions like $|0202...\rangle\rightarrow|2020...\rangle$, reminiscent of spin-$1/2$ PXP.

We confirm this by studying the fidelity time series, shown in Fig. S6, for two initial states, $|\psi_{0}\rangle=|1010\ldots\rangle$ (top) and $|\psi_{0}\rangle=|2020\ldots\rangle$. As in the main text, to get a clearer insight into the dynamics, we plot the generalized fidelity, $|\langle\phi|\exp(-itH)|\psi_{0}\rangle|^{2}$, where for $|\phi\rangle$ we do not necessarily pick the initial state, but any $|n0n0\ldots\rangle$, $n=1,2,3$, and their translated copies, $|0n0n\ldots\rangle$. We notice that the predominant dynamics of this Hamiltonian is to drive oscillations between the $|1\rangle$, $|3\rangle$ states and the $|0\rangle$, $|2\rangle$ states, with a small leakage into the remaining states. When quenched from $|0101...\rangle$, as $0\rightarrow 2$ is blocked due to the Rydberg constraint, the resulting dynamics is primarily just a free precession of the decoupled spin-$1/2$ chain, with a small leakage. This is observed in Fig. S6, where we see transitions between $|0101...\rangle\rightarrow|0303...\rangle$ without passing through $|0000....\rangle$. In contrast, $0\rightarrow 2$ flips are permitted when quenching from the state $|2020\ldots\rangle$. This allows the state to evolve through the polarized state $|0000...\rangle$ and exchange sublattices, undergoing oscillations between $|0202...\rangle\rightarrow|2020...\rangle$. This is precisely the same dynamics of spin-$1/2$ PXP, neglecting the leakage. As we see in Fig. S6 (bottom), we obtain non-zero fidelity only with the states $|2020\ldots\rangle$ and $|0202\ldots\rangle$, suggesting that the 4-color chiral clock model contains an embedded PXP spin-1/2.

Finally, we note that CCMs, as defined in Ref. Fendley (2014), appear to have scars only in $N_{c}\leq 4$ cases. For odd $N_{c}\geq 5$, we did not find revivals from Néel-type initial states. However, we note that for $N_{c}$-even it is possible to construct a special case of our clock models by choosing $k_{n}\in\{-N_{c}/2,...,-1,1,...N_{c}/2\}$. These models do appear to support scars, and for $N_{c}=4$ this choice of $k$’s reduces to CCM in Eq. (S11). Nevertheless, for $N_{c}>4$ longer range couplings are introduced and the models are no longer equivalent to Ref. Fendley (2014). The scarring in these models can similarly be explained by studying the single-site free precession, where dominant cyclic transitions between basis states emerge, with cycles shorter than $N_{c}$. Ignoring leakage away from these paths, they are effectively the same as the corresponding clock model with a matching cycle length.

SV Analytic prediction of parameters of 4-color scarred models

In the main text we presented a phase diagram of $N_{c}=4$ models of the form $H=\sum_{n}P^{0}_{n-1}C_{n}P^{0}_{n+1}$, which were shown to contain strong revivals. We noted that optimal models in these diagrams could be roughly predicted by determining the values of parameters that yield a single-site Hamiltonian $C$ with evenly spaced eigenvalues, such that the single-site dynamics is periodic. Further models with weaker revivals are obtained for $C$ whose eigenvalue spacings differed by a ratio of 2. In the following we derive analytical expressions for the lines in Fig. 2 in the main text.

Varying the matrix elements of $C$, we consider all Hermitian, purely imaginary, non-diagonal $C$, in analogy with the clock and spin Hamiltonians discussed in the main text. We are free to rescale $C$ such that one coupling has magnitude $1$, so for $N_{c}=4$, $5$ couplings can be varied. We consider two slices of the phase diagram: (A) vary the next-nearest neighbor hoppings, $C_{02}=C_{13}=\alpha i$, while also varying $C_{03}=-\beta i$; (B) switch off next-nearest-neighbor hoppings while varying $C_{12}=-\alpha i$ and $C_{03}=-\beta i$. Explicitly:

[TABLE]

The characteristic eigenvalue equation for both these matrices reduces to the form

[TABLE]

with $a$ and $b$ being:

[TABLE]

It follows that the eigenvalues $E$ take the form

[TABLE]

These eigenvalues are symmetric about zero, such that the spacings $\Delta E_{1}=\Delta E_{3}$, with $\Delta E_{n}=E_{n}-E_{n-1}$. The analytic prediction for robust scarred models, based on the single-site analysis, then reads:

• •

Equidistant eigenvalues:

[TABLE]

• •

Commensurate eigenvalues with ratio 2:

[TABLE]

• •

Degenerate eigenvalues:

[TABLE]

These equations are plotted in terms of parameters $\alpha,\beta$ in Fig S7. These lines give the band of scarred models in our $N_{c}=4$ phase diagram in Fig. 2 in the main text.

SVI Open Boundary Conditions

In the PXP spin-1/2 model Turner et al. (2018a, b), open boundary conditions have little effect on the scarring physics. Special scarred eigenstates still exist in the open PXP chain, and oscillatory quench dynamics still follow from certain initial states (Néel), while rapid thermalization is observed from generic quenches. Fig S8 summarizes the dynamics of the $N_{c}=3$ colored clock, for the largest accesible system size with open boundary conditions ($N=14,dim(\mathcal{H})=21845$). Numerically, we observe the same scarring phenomenology as with periodic boundary conditions, with the same pattern of transitions $|0101...\rangle\rightarrow|0202...\rangle\rightarrow|1010...\rangle$. The only difference, as in PXP spin-1/2, is that revivals decay slightly faster in an open chain.

SVII Finite Size Scaling of Numerical Results

The system sizes accesible to us are quite limited, due to relying on exact diagonalization techniques to obtain all eigenstates of systems with fairly large local Hilbert space dimensions. While our models could be studied using MPS/DMRG techniques White (1992); Vidal (2007), we note such techniques can only be used to simulate dynamics at short times. By contrast, important information about the atypical eigenstates of the system are currently only obtainable via exact diagonalization methods.

Fig S9 presents a finite size scaling analysis of our clock models, illustrating that our main results are well-converged in the accessible system sizes. Bottom panel of Fig S9 shows $t_{0}(N)$, the time at which the first maximum of the fidelity $F(t)=|\langle 1010...|e^{-iHt}|1010...\rangle|^{2}$ occurs, while $f_{0}$ is the amplitude of this revival, $f_{0}=F(t_{0})$. More precisely, top panel of Fig S9 shows the fidelity density, $\mathcal{F}(N)=\ln(f_{0}(N))/N$, which has a well-defined thermodynamic limit. The largest system sizes accesible to us were $N=20$ for $N_{c}=3$ and $N=16$ for $N_{c}=4$. Resolving both translational and parity symmetry, the Hilbert space dimension of the relevant sectors are of the order of $2\times 10^{4}$ for both $N_{c}=3$ and $N_{c}=4$.

SVIII Approximating scarred eigenstates

In the PXP spin-1/2 model Turner et al. (2018a, b) an approximation scheme called “the forward scattering approximation” (FSA) was introduced in order to approximate the $N+1$ scarred eigenstates in a chain of size $N$. This approximation was shown to accurately capture the eigenstates that have highest overlaps with the Néel state and whose energies are equally spaced. The FSA is formally based on the splitting of the Hamiltonian into two parts, $H=H^{+}+H^{-}$, where $H^{+}$ increases the Hamming distance from the Néel state by one (and $H^{-}$, conversely, decreases this distance by one). The FSA is defined by projecting the PXP Hamiltonian into the $(N+1)$-dimensional basis spanned by vectors $(H^{+})^{n}|1010\ldots\rangle$, $n=0,1,\ldots,N$. By construction, these vectors are orthogonal to each other. Diagonalizing the PXP Hamiltonian in this basis produces approximate eigenenergies and eigenvectors, which were found to match the exact ones with high accuracy. Moreover, the projection of the Hamiltonian into the FSA basis can be performed analytically, allowing one to reach much larger sizes than in exact diagonalization. Finally, the FSA scheme also provided an intuitive picture for the quench dynamics in terms of the wave function spreading over the Hilbert space graph Turner et al. (2018a).

Unfortunately, direct generalization of the FSA to $N_{c}>2$ clock models is not trivial. In retrospective, in the PXP case, the splitting of the Hamiltonian into $H^{+}$ and $H^{-}$ is straightforward because we can rely on the conventional definition of the Hamming distance. In the higher color models, this is no longer sufficient; moreover, it may be necessary to split the Hamiltonian into more than two parts. Nevertheless, in Fig. S10 we demonstrate that a more general approximation scheme yields accurate results for the $N_{c}=3$ clock model. Instead of splitting $H$, we follow the standard Lanczos procedure and build the Krylov subspace using full $H$: $\mathcal{K}=\{H^{n}|1010\ldots\rangle,n=0,1,2,\ldots\}$. Compared to the FSA, the Krylov subspace may contain an arbitrary number of vectors and we need to explicitly impose a cutoff. Furthermore, the Krylov vectors are no longer necessarily orthogonal, and Gram-Schmidt orthogonalization needs to be explicitly performed. We note that this approximation scheme is identical to the standard Lanczos algorithm for approximating the ground state of a system. In that case, however, the size of the Krylov subspace can be as large as as required to reach a target precision, while in our case we stop the algorithm after a given number of steps. Moreover, we obtain approximations to all scarred states in a single run of the algorithm.

Fig. S10(top) shows the result of the Krylov approximation for the 3-color clock model. We see that the approximation accurately predicts the energies of scarred eigenstates, as well as the eigenstate overlap with $|1010\ldots\rangle$ product state. Fig. S10(bottom) also shows the overlap of the predicted scarred wave functions with those obtained by exact diagonalization. The approximation is very accurate for the ground state, which is expected given that we are effectively using the Lanczos algorithm. Neglecting the zero energy state (which is not uniquely defined due to large degeneracy at $E=0$), the average overlap with the exact states for $N_{c}=3,N=12$ is $0.7410$. This is considerably smaller compared to PXP $1/2$, where the same procedure results in an average accuracy of $0.9173$. Moreover, the explicit orthogonalization that needs to be performed numerically restricts the Krylov approximation scheme to much smaller system size compared to the FSA.

SIX Improving Revivals

We noted in the main text that revivals are slightly more robust in PCP models over PXP models for odd $N_{c}$. Revivals can be substantially improved further, for any $N_{c}$, by the addition of certain perturbations, which we discuss here.

For PXP spin-1/2, it has previously been shown Khemani et al. (2019); Choi et al. (2018) that adding the family of perturbations

[TABLE]

can substantially improve revivals. Following the above discussion on the FSA approximation, the Hamming operators $H^{+},H^{-}$ approximately obey an su(2) algebra, in the sense that $[H^{+},H^{-}]=2iH^{Z}+Q$, where $Q$ is the deviation from a perfect su(2) algebra. One can show the perturbations $H_{d}^{PPXP}$ can reduce this deviation, such that oscillations from the Néel state can be improved and interpreted as precession of a large ($s=N/2$) spin Choi et al. (2018).

Further, we have numerically found a different family of perturbations

[TABLE]

with $d=3,5,7,\ldots$, which also improves revivals. Taking the same form of perturbations, but letting $X\rightarrow C$, revivals in clock models can also be substantially improved. We find the optimal coefficients for the 1st order perturbations using golden section search, taking as a cost function the first revival fidelity deficit, $1-F=1-|\langle\psi(0)|\psi(t)\rangle|^{2}$. Optimal values for $N_{c}=2,N=22$ were $0.10856$ (PPXP), $0.12296$ (PXXXP), while for $N_{c}=3,N=16$ optimal values were $0.09443$ (PPCP), $0.09112$ (PCCCP).

We note here that PPCP perturbations appear to improve the first fidelity revival to some threshold, and revivals no longer decay after this. Intuitively, this suggest all components of the wave function not residing in some subspace is rapidly shed, followed by perfect oscillations in this subspace. Via a rotation, we know $H=PCP+\lambda(PPCP+PCPP)=P^{{}^{\prime}}XP^{{}^{\prime}}+\lambda(P^{{}^{\prime}}P^{{}^{\prime}}XP^{{}^{\prime}}+P^{{}^{\prime}}XP^{{}^{\prime}}P^{{}^{\prime}})$. Thus we also expect an emergent $SU(2)$ subspace, whose algebra is specified by $H^{+}$,$H^{-}$ now given in terms of $P^{{}^{\prime}}$. However, in the rotated basis, the initial state $|0101...\rangle_{clock}$ is now a highly entangled state. Thus it appears any component not in the emergent $SU(2)$ subspace will dephase away and we are left with perfect revivals in this subspace.

In contrast, PCCCP perturbations can improve fidelity revivals to a greater degree for short times, but fidelity can still be observed to decay at longer times. This suggests the improvements comes from an orthogonal mechanism, not an emergent SU(2) symmetry. We leave this question for future work.

SIX.1 Proximity to Integrability

It has been suggested Khemani et al. (2019) that $H_{d}^{PPXP}$ perturbations, Eq. (S27), drive the spin-1/2 PXP Hamiltonian to an integrable point. The argument for this was the numerically observed minimum in the mean level statistics parameter $\langle r\rangle$ Oganesyan and Huse (2007), tending towards the Poisson value $\approx 0.39$. This can be observed in Fig S12, where we scan $\langle r\rangle$, resolved in the $k=0$, inversion symmetric sector, while varying the strength of perturbation. Vertical lines mark the optimal perturbation strength when optimizing for fidelity revivals. For PXP, both perturbations, PPXP and PXXXP, result in a dip of $\langle r\rangle\approx 0.39$, although the optimal perturbation strength with respect to fidelity moves the system further from this minimum Choi et al. (2018).

In contrast, no similar features is observed for the $N_{c}=3$ colored clock (Fig S13), where $\langle r\rangle\approx 0.53$ for all values scanned, including no perturbation and the optimal perturbation. While we cannot rule out a proximate integrable point for the clock models after only scanning along two axes of the vector space of all possible perturbations, it seems that optimal revivals are not correlated with integrability.