Collapse and revival of quantum many-body scars via Floquet engineering

TL;DR

This paper demonstrates how periodic driving in a Rydberg atom chain can induce transitions between ergodic and non-ergodic quantum regimes, revealing a method to control quantum many-body scars and long-term coherence.

Contribution

It introduces Floquet engineering as a tool to tune quantum many-body scars and ergodic behavior in a driven Rydberg system, showing reentrant transitions.

Findings

Drive frequency controls ergodic and non-ergodic regimes.

Floquet spectrum contains scars in non-ergodic regime.

Long-time coherence is linked to Floquet scars.

Abstract

The presence of quantum scars, athermal eigenstates of a many-body Hamiltonian with finite energy density, leads to absence of ergodicity and long-time coherent dynamics in closed quantum systems starting from simple initial states. Such non-ergodic coherent dynamics, where the system does not explore its entire phase space, has been experimentally observed in a chain of ultracold Rydberg atoms. We show, via study of a periodically driven Rydberg chain, that the drive frequency acts as a tuning parameter for several reentrant transitions between ergodic and non-ergodic regimes. The former regime shows rapid thermalization of correlation functions and absence of scars in the spectrum of the system's Floquet Hamiltonian. The latter regime, in contrast, has scars in its Floquet spectrum which control the long-time coherent dynamics of correlation functions. Our results open a new possibility of drive frequency-induced tuning between ergodic and non-ergodic dynamics in experimentally realizable disorder-free quantum many-body systems.

Figures

Equations

$$H_{\mathrm{RYD}}$$

$$H_{\mathrm{RYD}}^{\prime }$$

$$H_{\mathrm{spin}}$$

$$U(T,0)$$

$$O_{ij}$$

$$H_0$$

$$H_1$$

$$H_F=-w\frac{\sin \gamma }{\gamma }\sum _j[\cos \gamma {\tilde{\sigma }}_j^x+\sin \gamma {\tilde{\sigma }}_j^y].$$

$$f_{1[2]}$$

$$f_1=w^2\frac{{\sin }^2(\gamma )}{{\gamma }^2},\text{mbox}{f}_2=0.$$

$$H$$

$$H_d=-w\sum _{\ell }(d_{\ell }+d_{\ell }^{†})+\lambda \sum _{\ell }{n}_{\ell },$$

$$H_{\mathrm{RYD}}$$

$$H_{\mathrm{RYD}}^{\prime }$$

$${\sigma }_{\ell }^z$$

$$H_{\mathrm{spin}}$$

$$U=e^{-i{H}_{+}T/(2\hslash )}{e}^{-i{H}_{-}T/(2\hslash )}=e^{X_{+}}{e}^{X_{-}}=U_{+}{U}_{-},$$

$$X_{1[2]}$$

$$\ln [e^{X_{+}}{e}^{X_{-}}]$$

$$H_F^0$$

$$H_F^2$$

$$C_{1a}$$

$$C_{1b}$$

$$H_{F3}^{(1)}$$

$$H_{F3}^{(2)}$$

$$H_{F3}^{(3)}$$

$$H_{F3}^{(4)}$$

$$H_2$$

$$i\hslash \frac{\partial |n(t)⟩}{\partial t}=H(t)|n(t)⟩$$

$$|n(T)⟩=e^{-i{\theta }_n}|n(0)⟩.$$

$$e^{-i{\theta }_n}=e^{-(i/\hslash ){\int }_0^{T}dt{E}_n(t)}.$$

$$|n(t)⟩=\sum _m{c}_m(t)e^{-(i/\hslash ){\int }_0^{t}d{t}^{\prime }{E}_m(t^{\prime })}|m⟩,$$

Full text

The eigenstate thermalization hypothesis (ETH) is one of the central paradigms for understanding out-of-equilibrium dynamics of closed non-integrable quantum systems rev1a ; rev1b ; rev1c ; rev1d ; rev2 ; deutsch1 ; srednicki1 ; rigol1 . It posits that all bulk eigenstates of a generic quantum many-body Hamiltonian are thermal; their presence ensures ergodicity and leads to eventual thermalization for out-of-equilibrium dynamics of a generic many-body state rev2 . This hypothesis is strongly violated in certain cases, the most famous example being one-dimensional (1D) disordered electrons in their many-body localized phase mblref1 ; mblref2 . More recently another example of a weaker failure of ETH, due to the presence of quantum many-body scar states, has been studied extensively in disorder-free systems scarrefqm1 ; scarref1 ; scarref2a ; scarref2b ; scarref2c ; scarref2d ; scarref2e ; scarref3a ; scarref3b ; scarref3c ; scarref3d ; scarref3e . Scars are eigenstates with finite energy density but anomalously low entanglement scarref1 ; scarref2b ; scarref2c ; scarref2d ; scarref3b which form an almost closed subspace in the system’s Hilbert space under the action of its Hamiltonian. Their presence leads to persistent coherent oscillatory dynamics of correlation functions starting from initial states that have a high overlap with scars. This provides an observable consequence of their presence as verified in recent experiments on quench dynamics of a chain of ultracold Rydberg atoms scarref1 .

Here we study the fate of such ergodicity violation in a periodically driven Rydberg chain. It is well known that the stroboscopic dynamics of a periodically driven quantum system is controlled by its Floquet Hamiltonian $H_{F}$ rev3 which is related to its unitary evolution operator $U$ through $U(T,0)=\exp[-iH_{F}T/\hbar]$, where $T=2\pi/\omega_{D}$ is the time period of the drive and $\omega_{D}$ is the drive frequency. For generic disorder-free systems, such driving is expected to cause thermalization to a featureless “infinite temperature” state steady1b ; steady1c ; steady2 ; steady3 . In what follows, we will study the possibility of the existence of scars in the eigenstates of $H_{F}$ as a function of $\omega_{D}$ and relate their influence on the dynamics of correlation functions. Our initial state will be an experimentally realized ${\mathbb{Z}}_{2}$ symmetry broken many-body state which has one Rydberg excitation in alternate lattice sites scarref1 ; rydramp1 ; rydramp2 ; rydramp3 .

The central results of this study are as follows. First, for large $\omega_{D}$ and starting from a initial $\mathbb{Z}_{2}$ state, we show the presence of long-time persistent oscillations of the density-density correlator of Rydberg atoms. Such oscillations have characteristic frequencies which are different from $\omega_{D}$ indicating a lack of synchronization (a hallmark of thermalization in periodically driven systems). We relate this oscillation frequency to the quasienergy separation between the scar states of the Floquet Hamiltonian indicating the central role of these states in the dynamics. Second, at ultra-low drive frequencies, we find that there are no persistent oscillations, and the behavior of the correlator agrees with that expected from ETH. We show numerically that in this regime, there are no scars in the eigenspectrum of $H_{F}$ and the dynamics is controlled by a set of thermal states. Finally, we find several drive-frequency-induced transitions between thermal and coherent regimes at intermediate frequencies. These transitions, that have no analogs in the non-driven systems studied earlier scarref1 ; scarref2a ; scarref2b ; scarref2c ; scarref2d ; scarref2e ; scarref3a ; scarref3b ; scarref3c ; scarref3d ; scarref3e , provide a route to controlled switching between ergodic and non-ergodic dynamics of the Rydberg atoms. We chart out the critical drive frequencies at which these transitions occur, provide an analytic understanding of their occurrence, and suggest experiments which can test our theory.

Model: The low-energy properties of an ultracold Rydberg atom chain can be described by an effective two-state Hamiltonian on each site given by scarref1 ; rydramp1 ; rydramp2 ; rydramp3

[TABLE]

The two states correspond to the ground ($|g_{i}\rangle$) and Rydberg excited states ($|e_{i}\rangle$) of the atoms on site $i$. The dipole blockade in these systems ensures that there is at most one Rydberg excitation per site: $n_{r}^{i}\leq 1$, where $n^{r}_{i}=(1+\sigma^{z}_{i})/2$ is the Rydberg excitation number operator, and $\sigma_{i}^{x}=|g_{i}\rangle\langle e_{i}|+|e_{i}\rangle\langle g_{i}|$ denotes a Pauli matrix on site $i$ which couples the ground and excited states. In Eq. (1), $\Delta$ is the detuning parameter which can be used to excite an atom to its Rydberg state, $V_{ij}$ denotes an interaction between two Rydberg excitations, and $\Omega$ is the coupling strength between ground and excited states. In experiments, it is possible to reach a regime where $V_{i,i+1}\gg\Omega,\Delta\gg V_{i,i+2}$ so that the Hamiltonian of the model becomes scarref1 ; dipoleexp1 ; dipoleexp2 ; subir1 ; subir2

[TABLE]

and this is to be supplemented by the constraint that $n^{r}_{i}n^{r}_{i+1}=0$ for all sites $i$.

This constrained model can be easily mapped into an Ising-like spin model in the presence of both longitudinal and transverse fields of strength $\Delta$ and $\Omega$ respectively. Within the constrained Hilbert space of the system, one can represent $H_{\rm RYD}$ as scarref2c ; scarref2d

[TABLE]

where $P_{i}=(1-\sigma_{i}^{z})/2$ is a local projection operator, ${\tilde{\sigma}}^{\alpha}_{i}=P_{i-1}\sigma_{i}^{\alpha}P_{i+1}$ and $\alpha=x,y,z$, $\Omega\equiv-w$, and $\lambda\equiv\Delta$. Our analysis will be based on this model. Eq. (3) also provides a low-energy description for the 1D tilted Bose-Hubbard model as detailed in the Supplementary Information. For $\lambda=0$, $H_{\rm spin}$ reduces to the “PXP model” studied in Refs. scarref2c ; scarref2d which is known to host quantum scars among its eigenstates.

Analysis: We analyze the periodic dynamics of $H_{\rm spin}$ for a square pulse protocol: $\lambda(t)=-(+)\lambda$ for $t\leq(>)T/2$. The unitary evolution operator at the end of a drive cycle can be written as $U(T,0)=e^{-iH_{\rm spin}[\lambda]T/2}e^{-iH_{\rm spin}[-\lambda]T/2}$. The evolution operator can then be expressed as

[TABLE]

where $\epsilon_{p}^{+(-)}$ and $|p^{+(-)}\rangle$ are eigenstates and eigenfunctions of $H_{\rm spin}[+(-)\lambda]$ and $c_{pq}^{-+}=\langle p^{-}|q^{+}\rangle$. The spin correlation function $O_{ij}=\langle(1+\sigma^{z}_{i})(1+\sigma^{z}_{i+j})\rangle/4$ of the spins between any two sites $i$ and $i+j$ can then be obtained, after $n$ drive cycles, as

[TABLE]

where $c_{\psi_{0}p}^{a}=\langle\psi_{0}|p^{a}\rangle$ for $a=\pm$, and $|\psi_{0}\rangle$ is the initial state. Unless explicitly stated otherwise, we will choose $|\psi_{0}\rangle=|\mathbb{Z}_{2}\rangle=|\cdots\downarrow\uparrow\downarrow\uparrow\cdots\rangle$ to be a $\mathbb{Z}_{2}$ symmetry broken state, with $\langle\psi_{0}|\sigma^{z}_{j}|\psi_{0}\rangle=(-1)^{j+1}$. We note that $O_{ij}$ provides direct information of the density-density correlation function of the Rydberg atoms after $n$ cycles of the drive. Our numerical analysis will involve computation of $\epsilon_{p}^{\pm}$ and $|p^{\pm}\rangle$ using exact diagonalization for finite chains of size $L\leq 26$ and subsequent evaluation of $O_{ij}$ using Eq. (5).

To obtain an analytical understanding of the nature of the dynamics, we derive the Floquet Hamiltonian corresponding to $U=\exp[-iH_{F}T/\hbar]$ using a Magnus expansion which is expected to yield an accurate description of the dynamics for high drive frequencies rev3 . Further details are provided in Supplementary Information. To $\mathcal{O}(1/\omega_{D}^{3})$, this calculation leads to $H_{F}^{\rm Magnus}=H_{0}+H_{1}$, where

[TABLE]

Here we have defined dimensionless quantities $\gamma=\lambda T/4\hbar$ and $\delta=wT/(4\hbar)$, $C_{1}=1-2\gamma^{2}/3$, and $C_{2}=\gamma[1-(\gamma^{2}-4\delta^{2})/3]$. We find that the terms in $H_{0}$ (which we denote as PXP terms) constitute a renormalized PXP model (up to a global spin rotation); consequently, for $\hbar\omega_{D}\gg\lambda,\delta$, where the effect of $H_{1}$ can be ignored, we expect $H_{F}^{\rm Magnus}$ to host scar states similar to those in the PXP model. However at moderate $\omega_{D}$, the terms in $H_{1}$ (which we denote as non-PXP terms), are expected to become important. The competition between these two classes of terms can be tuned using the drive frequency and will be discussed in detail below.

In what follows, we will be interested in large drive amplitudes for which $\lambda\gg w$ ($\gamma\gg\delta$). In this regime, as detailed in Supplementary Information, the Floquet Hamiltonian can be perturbatively calculated to $\mathcal{O}(w)$ for an arbitrary $\omega_{D}$ and gives only PXP terms with all non-PXP terms (and further PXP terms) generated at $\mathcal{O}(w^{2}/\lambda)$ and beyond. To $\mathcal{O}(w)$, the Floquet Hamiltonian equals

[TABLE]

Eq. (7) will be used to understand the transitions between ergodic and non-ergodic regimes.

Results: To demonstrate the presence of ergodic to non-ergodic transitions as a function of the drive frequency $\omega_{D}$, we first compute the dynamics of the correlators starting from $|\mathbb{Z}_{2}\rangle$. For this, we perform exact diagonalization and compute $O_{22}$ from Eq. (5) as a function of the stroboscopic time $n$ (number of drive cycles) for several $\omega_{D}$. In addition we also compute the half-chain entanglement entropy $S_{L/2}$ for the eigenstates of the Floquet spectrum by obtaining these via numerical diagonalization of $U$ (Eq. (4)). This is followed by a computation of the reduced density matrix for these eigenstates for a half chain from which $S_{L/2}$ can be obtained using a standard procedure (see Supplementary Information). Quantum scars have $S\sim\ln L$ and are thus expected to have much lower entanglement compared to thermal states for which $S\sim L$. Thus $S_{L/2}$ provides a reliable way to distinguish between thermal states and scars for a finite-size many-body system.

The results of these calculations are shown in Fig. 1. Panel a [b] of Fig. 1 shows the behavior of $O_{22}$ as a function of $n$ for $\hbar\omega_{D}\gg[\ll]\lambda,w$. We find that for $\hbar\omega_{D}\gg\lambda$, the dynamics exhibits long-time coherent oscillations as expected from the quench dynamics of the PXP model at $\lambda=0$ studied earlier scarref2c ; scarref2d . The frequency of these oscillations differs from $\omega_{D}$ indicating a clear lack of synchronization. This behavior is expected from Eq. (6) where the non-PXP terms appear in ${\mathcal{O}}(1/\omega_{D}^{3})$ and hence are small. The presence of scars in the Floquet Hamiltonian in this regime may be confirmed from Fig. 1 (c)which shows $S_{L/2}$ for Floquet eigenstates (denoted by $|\Phi_{n}\rangle$ henceforth) as a function of the Floquet quasienergies $E_{F}$. The scar states are seen as clear outliers in this plot. The eigenstates $|\Phi_{n}\rangle$ with large overlaps with $\mathbb{Z}_{2}$ ($|\langle\mathbb{Z}_{2}|\Phi_{n}\rangle|^{2}\geq 0.01$) are circled in red; from this we find that the scars have maximal overlap with $|\mathbb{Z}_{2}\rangle$ and thus control the dynamics leading to violation of ETH scarref2d .

In contrast, for $\hbar\omega_{D}/w\ll 1$, all the states including those controlling the dynamics are thermal (Fig. 1 (d)). Consequently, there are no persistent oscillations for $O_{22}$ (Fig. 1 (b)) and one finds thermalization consistent with ETH. We also note that the oscillatory behavior seen in Fig. 1 is a property of the initial $|\mathbb{Z}_{2}\rangle$ state; a similar study of the dynamics for any drive frequency starting from the Rydberg vacuum state $|0\rangle=|\downarrow\downarrow\downarrow...\rangle$ always provides fast thermalization consistent with ETH (red curves in Fig. 1 (a) and (b)).

The quasienergy separation $w_{R}$ of the scars for $\hbar\omega_{D}\geq\lambda$ as a function of $\omega_{D}$ is shown in Fig. 2. We note that $w_{R}$ starts to decrease when $\gamma\to 1$; such a behavior follows from the decrease of the norm of the PXP terms in $H_{F}$ at $\mathcal{O}(w)$ with increasing $\gamma$ (Eq. (7)) which gives $w_{R}=(\sin(\gamma)w_{\infty})/\gamma$, where $w_{\infty}$ is the scar quasienergy separation for the undriven PXP model. Fig. 2 shows the exact match of the numerical result with that obtained analytically. The lowering of $w_{R}$ implies a sharp decrease in the oscillation frequency $\omega_{\rm res}$ of $O_{22}$ as a function of $\omega_{D}$ in this regime. To check this, we extract $\omega_{\rm res}$ as a function of $\omega_{D}$ from the Fourier transform of $O_{22}$ (upper inset of Fig. 2) which matches the corresponding values of $w_{R}/\hbar$ almost perfectly (lower inset of Fig. 2) and shows a clear decrease with $\omega_{D}$. This provides a drive-induced control over the quasienergy separation of the scars and hence on the oscillation frequency which has no analog in earlier quench studies. (Interestingly, the bottom inset of Fig. 2 shows a large number of states with zero quasienergy. They arise due to a symmetry of the Floquet operator as discussed in Supplementary Information).

Next, we analyze the regime $\hbar\omega_{D}<\lambda$ where we encounter the reentrant transitions between coherent and thermal regimes. Here, we follow Ref. scarref2c ; scarref2d and use the state $|\Psi_{0}\rangle=(|\mathbb{Z}_{2}\rangle+|{\bar{\mathbb{Z}}}_{2}\rangle)/\sqrt{2}$ as our initial state, where $|{\bar{\mathbb{Z}}}_{2}\rangle$ denotes the spin-flipped version of $|\mathbb{Z}_{2}\rangle$. This allows us access to larger chain length $L\leq 26$ since $|\psi_{0}\rangle$ has weight only in the sector with zero total momentum and spatial inversion (parity) symmetry.

The result of evolution of $O_{22}$ in this subspace is shown in Fig. 4 near the first reentrant transition. Fig. 4 (a) shows non-ergodic persistent oscillatory dynamics at $\omega_{D}=8.25$. As we reduce $\omega_{D}$, these oscillations dampen (Fig. 3 (b)); such a behavior can be interpreted as a precursor to ergodic dynamics and thermalization. Upon further reduction of $\omega_{D}$, ergodic dynamics consistent with ETH sets in and the fastest thermalization is seen around $\omega_{D}=7.75$ (Fig. 3 (c)). Finally, at lower $\omega_{D}$, the persistent oscillations return (Fig. 3 (d)).

Such transitions between non-ergodic and ergodic dynamics as a function of drive frequency can be tied to the presence or absence of scars in the spectrum of $H_{F}$. This is shown in Fig. 4. Figs. 4 (a)($\omega_{D}=8.25$) and (d) ($\omega_{D}=7.5$) clearly indicate the presence of scars having high overlap with $(|\mathbb{Z}_{2}\rangle+|{\bar{\mathbb{Z}}}_{2}\rangle)/\sqrt{2}$. This is consistent with the presence of non-ergodic dynamics characterized by persistent long-time oscillations (Figs. 3 (a, d)). These scar states start to merge with the thermal band around $\omega_{D}=8.0$ (Fig. 4 (b)) indicating precursor to the thermal behavior (Fig. 3 (b)). Fig. 4 (c) at $\omega_{D}=7.75$ shows complete absence of scars resulting in ergodic dynamics of $O_{22}$ and fast thermalization predicted by ETH (Fig. 3 (c)).

To obtain a qualitative understanding of these transitions, we compute the norm of the PXP-like terms in the Floquet Hamiltonian numerically using $L=14$. To this end, we write the matrix representation of $H_{F}$ in the basis states $|\phi_{n}\rangle$ of $\sigma^{z}$ and identify the matrix elements that have $\langle\phi_{n}|\sum_{\ell}\tilde{\sigma_{\ell}}^{x/y}|\phi_{m}\rangle\neq 0$. Let us denote this set as ${\mathcal{N}}_{0}$ which has $N_{0}=2LF_{L-1}$ elements (where $F_{L}$ is a Fibonacci number defined by $F_{n}+F_{n+1}=F_{n+2}$ with $F_{1}=F_{2}=1$). We then define $f_{1}$ and $f_{2}$ as

[TABLE]

Clearly, $f_{2}$ represents the contribution from the non-PXP type of terms in $H_{F}$. We note that in general $f_{1}$ will also have contributions from non-PXP terms since such terms may have non-zero matrix elements for some states included in ${\mathcal{N}}_{0}$. However, at large $\lambda/w$, the contributions from these terms are expected to be small by at least ${\mathcal{O}}(w^{2}/\lambda^{2})$. In fact from Eq. (7), to leading order, we find that

[TABLE]

To numerically verify that this is indeed the case, we plot $f_{1}$ as a function of $\hbar\omega_{D}/\lambda$ for several $\lambda/w$ (Fig. 5 (a)). These curves coincide indicating that $f_{1}$ is almost independent of $\delta$. Thus in this regime $f_{1}$ receives negligible contributions from the non-PXP terms in $H_{F}$ which necessarily depend on $\delta$ (Fig. 5 (a) also shows that $\lambda\sim 10$ lies outside this regime).

Fig. 5 (a) and Eq. (9) also show that $f_{1}$ is an oscillatory function of $\omega_{D}$. From the inset of Fig. 5 (a), we find that $f_{1}\simeq f_{2}$ near the minima of $f_{1}$ at $\hbar\omega_{D}=\lambda/(2n_{0})$ (Eq. (9)) where $n_{0}$ is a positive integer; in other regions, $f_{1}\gg f_{2}$. The ergodic dynamics of $O_{22}$ always occur in a finite frequency interval around $\omega_{D}^{c}$ (where $f_{1}=f_{2}$). This observation sheds light on the cause of the transitions. We note that the PXP Hamiltonian supports scars; consequently, for $\omega_{D}$ where $H_{F}$ is PXP-like ($f_{1}\gg f_{2}$), one expects the presence of scars among its eigenstates. These scars lead to coherent non-ergodic dynamics. In contrast, near the minima of $f_{1}$, where $f_{1}\simeq f_{2}$, $H_{F}$ receives significant contributions from the non-PXP terms. In the presence of such terms which can be long-ranged at low $\omega_{D}$, $H_{F}$ does not support scars. The bulk of its eigenstates around $\omega_{D}=\omega_{D}^{c}$ are thermal; consequently the dynamics of $O_{22}$ displays ergodic behavior consistent with the prediction of ETH. In fact, the level spacing statistics of the eigenvalues of $H_{F}$ follow a Gaussian orthogonal ensemble for $\omega_{D}\approx\omega_{D}^{c}$ as expected for an ergodic system (see Supplementary Information). The behavior of the fidelity $|\langle\Psi_{n}|\Psi_{0}\rangle|^{2}$, where $|\Psi_{n}\rangle$ is the state after $n$ drive cycles, across the transition is also discussed in the Supplementary Information.

The schematic phase diagram for these reentrant ergodic to non-ergodic transitions is summarized in Fig. 5 (b). The ergodic regions, where $O_{22}$ displays either complete ergodic behavior or a precursor to thermalization, are located in a small frequency window (shown schematically in grey) around $\omega_{D}^{c}$ (indicated as red dots joined by blue dashed lines in Fig. 5 (b)). The yellow lines indicate the positions of the transitions as obtained from Eq. (7) (i.e., $\lambda/\hbar\omega_{D}$ is an even integer), while the white regions denote the ranges of $\omega_{D}$ where $O_{22}$ shows non-ergodic oscillatory dynamics due to the presence of scars in $H_{F}$. The thermal regions become denser with decreasing $\omega_{D}$ and ultimately merge into a continuum at sufficiently small $\omega_{D}$ where non-ergodic coherent dynamics of $O_{22}$ ceases to exist.

Discussion: To summarize, we have studied the kinematically constrained PXP model, a paradigmatic model for many-body eigenstates called quantum scars that violate ETH, in the presence of a pulsed transverse magnetic field that varies periodically in time. In the limit of large drive amplitude of the field, the instantaneous Hamiltonian does not have any scars but the corresponding Floquet Hamiltonian that controls the stroboscopic dynamics of local quantities can still host them depending on the drive frequency. We find (a) the presence of several non-ergodic (characterized by a coherent oscillatory behavior of density-density correlators and scars in the Floquet Hamiltonian) and ergodic (characterized by a thermal non-oscillatory behavior of density-density correlators and absence of scars) regimes as a function of the driving frequency, and (b) the possibility of tuning the quasienergy spacing of the scars in the non-ergodic regime as a function of the drive frequency to control the frequency of oscillations of the correlators. Both these features are entirely absent in the undriven PXP model and can be tested by standard experiments using finite-size Ryberg chains and starting from an initial state that has one Rydberg excitation at each alternate lattice site. In particular, these reentrant transitions from non-ergodic to ergodic regimes by tuning the drive frequency are possibly the first example of this kind in a system without any spatial disorder.

The mechanism for these transitions is also rather transparent in the large drive amplitude regime as discussed here. Using a Floquet perturbation theory, the Floquet Hamiltonian can be grouped into PXP and non-PXP type terms, with the non-PXP terms being suppressed by at least the inverse of the drive amplitude. The leading PXP terms can be resummed to all orders in the drive frequency which shows that these can diminish and become comparable to the non-PXP terms in the neighborhood of special drive frequencies leading to the emergence of the thermal regime. Lastly, on a theoretical front, it would be interesting to explore the influence of the small non-PXP terms in the non-ergodic regime to see whether they lead to eventual thermalization on a much longer time scale.

Acknowledgments: The work of A.S. is partly supported through the Max Planck Partner Group program between the Indian Association for the Cultivation of Science (Kolkata) and the Max Planck Institute for the Physics of Complex Systems (Dresden). D.S. thanks DST, India for Project No. SR/S2/JCB-44/2010 for financial support.

I Supplementary Material

II Model Hamiltonian

The Hamiltonian of the Bose-Hubbard model with a tilt is given by subir1

[TABLE]

where $b_{i}$ ($b^{\dagger}_{i}$) denotes the boson annihilation (creation) operator on site $i$ of a 1D chain, $n^{b}_{i}=b_{i}^{\dagger}b_{i}$ is the boson number operator, $w^{\prime}$ is the hopping amplitude of the bosons, $E_{0}$ denotes the magnitude of the tilt, $\mu$ is the chemical potential, $w$ is the hopping amplitude, and $U$ is the on-site interaction between the bosons. The tilt can be generated either by shifting the center of the parabolic trap confining the bosons or by applying a linearly varying Zeeman field which couples to the spin of the bosons. The latter variation can be made time dependent by using a magnetic field which varies periodically in time. It is well known that the low-energy physics of these bosons deep inside the Mott phase, whose occupation number is denoted by $n_{0}$, (here we focus on the case $n_{0}=1$) and where $U,E_{0}\gg w^{\prime},|U-E_{0}|$, is given by

[TABLE]

where $d_{\ell}=b_{i}^{\dagger}b_{j}/\sqrt{n_{0}(n_{0}+1)}$ denotes a dipole annihilation operator on link $\ell$ between sites neighboring $i$ and $j$ on a 1D lattice, $w=\sqrt{2}w^{\prime}$ for $n_{0}=1$, $n_{\ell}=d_{\ell}^{\dagger}d_{\ell}$ is the dipole number operator on link $\ell$, $w$ is the amplitude of spontaneous dipole creation or destruction, and $\lambda$ is the chemical potential for the dipoles. This dipole model is to be supplemented by two constraints which make it non-integrable: $n_{\ell}\leq 1$ and $n_{\ell}n_{\ell+1}=0$ for all links. The phase diagram of this model has been studied theoretically in Ref. subir1, and has also been experimentally verified dipoleexp1 . It is well-known that $H_{d}$ support a quantum phase transition at $\lambda_{c}=-1.31w$ separating a $\mathbb{Z}_{2}$ symmetry broken ground state ($|\mathbb{Z}_{2}\rangle$) for $\lambda<\lambda_{c}$ and a featureless dipole vacuum ($|0\rangle$) for $\lambda>\lambda_{c}$. The non-equilibrium dynamics of this model has also been studied for quench, ramp and periodic protocols dipoledyn .

The dipole model described in Eq. (S2) also serves as an effective model for describing the low energy physics of the Rydberg atoms. To see this we first consider the Hamiltonian of such atoms given by scarref1 ; scarref1p5 ; scarref2 ; scarref3 ; rydramp

[TABLE]

Here $n^{r}_{i}\leq 1$ denotes the number of Rydberg excitations on a given site, $\Delta$ is the detuning parameter which can be used to excite an atom to its Rydberg state, $V_{ij}$ denotes the interaction between two Rydberg excitations, $\sigma_{i}^{x}=|g_{i}\rangle\langle r_{i}|$ denotes the coupling between the ground ($|g\rangle$) and Rydberg excited ($|r\rangle$) states, and $\Omega$ is the corresponding coupling strength. In experiments scarref1 , it is possible to reach a regime where $V_{i,i+1}\gg\Omega,\Delta\gg V_{i,i+2}$; in this case, the Hamiltonian the model becomes equivalent to that of

[TABLE]

supplemented by the constraint that $n^{r}_{i}n^{r}_{i+1}=0$ for all sites $i$. Clearly, this model is equivalent to Eq. (S2) with the identification $n^{r}_{i}\to n_{\ell}$, $\Delta\to\lambda$ and $\Omega\to-w$.

Furthermore it is also easy to see that the dipole Hamiltonian (Eq. (S2)) is identical to the PXP model studied in Ref. scarref1, for $\lambda=0$. The simplest way to see this involves mapping of the dipole operators to Ising spins via the transformation

[TABLE]

where $\sigma^{\alpha}$ denote the Pauli matrices for $\alpha=x,y,z$. Moreover, the constraint of not having dipoles on adjacent links can be implemented via a local projection operator $P_{\ell}=(1-\sigma_{\ell}^{z})/2$ scarref1p5 ; scarref2 . Using these, one finds the spin Hamiltonian

[TABLE]

where $\tilde{\sigma}_{\ell}^{\alpha}=P_{\ell-1}\sigma_{\ell}^{\alpha}P_{\ell+1}$ for $\alpha=x,y,z$, and we have ignored an unimportant constant term while writing down the expression for $H_{\rm spin}$. The physics of $H_{\rm spin}$ within the constrained dipole Hilbert space is identical to that of $H_{d}$ and $H_{\rm RYD}$. At $\lambda=0$, $H_{\rm spin}$ reduced to the PXP model studied in Ref. scarref1, . Note that both the constraints of the dipole model are incorporated in $H_{\rm spin}$ via the local projection operators $P_{\ell}$.

Eq. (S6) has been used for the analysis in the main text.

III Magnus expansion

We consider the Hamiltonian $H_{\rm spin}$ given by Eq. (S6) in the presence of a periodic drive characterized by a square pulse protocol with time period $T=2\pi/\omega_{D}$, where $\omega_{D}$ is the drive frequency: $\lambda(T)=-(+)\lambda$ for $t\leq(>)T/2$. In what follows, we will chart out the details of the computation of the Floquet Hamiltonian of such a driven system using a high-frequency Magnus expansion.

For this protocol, the unitary matrix governing the evolution by a time period is given by

[TABLE]

where $X_{\pm}=(-iT/2\hbar)H_{\pm}$. For future reference, we also define $X_{1,2}$ given by

[TABLE]

such that $X_{\pm}=X_{1}\pm X_{2}$. The Floquet Hamiltonian can be obtained from $U$ as $H_{F}=(i\hbar/T)\ln U$. Using the Baker-Campbell-Hausdorff formula, one can express

[TABLE]

From Eq. (S9) we can find terms of different order in the Floquet Hamiltonian. The computation of these terms up to $O(1/\omega_{D}^{2})$ is straightforward and yields

[TABLE]

where $\gamma=\lambda T/(4\hbar)$. Note that these terms lead to a renormalized PXP model; it amounts to a change of magnitude of the coefficient $w$ of the standard PXP Hamiltonian and also a rotation in spin-space which depends on $\omega_{D}$. Note that the second order term in the expansion has an opposite sign compared to the zeroth order term. As we will see later, this is a general feature of the model; any term $\sim\gamma^{n}$ in the renormalized PXP model always comes with a opposite sign compared to a term $\sim\gamma^{n-2}$.

The first non-trivial longer-ranged terms in $H_{F}$ arises in $\mathcal{O}(1/\omega_{D}^{3})$. Its derivation involves some subtle issues. To see this, let us consider the commutator $C_{1}=[X_{+},[X_{+},X_{-}]]=C_{1a}+C_{1b}$. It is easy to see after a straightforward calculation that

[TABLE]

We note that within the constrained Hilbert space, any term with $\sigma_{\ell-1}^{\alpha}\sigma_{\ell}^{\beta}P_{\ell}$ identically vanishes for $\alpha,\beta=x,y$. Furthermore, the projection operators $P_{\ell}$ on any link satisfy $(1-P_{\ell\pm 1})\sigma_{\ell}^{\alpha}=0$ for $\alpha=x,y$. Using these results, we can simplify $C_{1b}$ to obtain

[TABLE]

Using Eq. (S12) and evaluating the necessary commutators, we finally get $H_{F}^{3}=H_{F3}^{(1)}+H_{F3}^{(2)}+H_{F3}^{(3)}+H_{F3}^{(4)}$, where

[TABLE]

Here $\delta=wT/(4\hbar)$, and we note that $\delta/\gamma=w/\lambda\ll 1$ in the limit of large $\lambda/w$. Thus Eqs. (S10) and (S13) yield $H_{0}$ in the main text while Eqs. (S14 - S16) yield $H_{1}$. This completes our derivation of the Magnus expansion to ${\mathcal{O}}(1/\omega_{D}^{3})$.

Before ending this section we note that if we concentrate on the large $\lambda/w$ limit, it is possible to compute higher-order corrections to the coefficients in $H_{F}^{0}$ and $H_{F}^{1}$. This can be seen by noting that at each order the contribution to such terms comes from $[X_{2},[X_{2},[X_{2},\cdots[X_{2},X_{1}]]]\cdots]$, i.e., the $n$-th order contribution involves commutator of $n-1$ terms involving $\sigma_{j}^{z}$ and one $\sigma_{j}^{x}$. These commutators provide the leading contribution in the large $\lambda/w$ limit. This structure allows us to compute leading higher-order terms in the Magnus expansion which contribute to the coefficients of the PXP term. A straightforward but cumbersome computation yields

[TABLE]

It will be shown in Sec. IV that the coefficients in $H_{2}$ can be resummed to yield a closed form valid for arbitrary $\omega_{D}$: $H_{F}=-(w\sin\gamma/\gamma)\sum_{j}(\cos\gamma~{}\tilde{\sigma}^{x}_{j}+\sin\gamma~{}\tilde{\sigma}^{y}_{j})$.

IV Floquet perturbation theory

We will now present a perturbation theory for a periodically driven system soori . We consider a Hamiltonian $H(t)=H_{0}(t)+V$, where $H_{0}(t)$ varies in time with a period $T=2\pi/\omega$, and $V$ is a small time-independent perturbation. We will assume that $H_{0}(t)$ commutes with itself at different times, and will work in the basis of eigenstates of $H_{0}(t)$ which are time-independent and will be denoted as $|n\rangle$, so that $H_{0}(t)|n\rangle=E_{n}(t)|n\rangle$, and $\langle m|n\rangle=\delta_{mn}$. We will also assume that $V$ is completely off-diagonal in this basis, namely, $\langle n|V|n\rangle=0$ for all $n$. We will now find solutions of the Schrödinger equation

[TABLE]

which satisfy

[TABLE]

For $V=0$, we have $|n(t)\rangle=e^{-(i/\hbar)\int_{0}^{t}dt^{\prime}E_{n}(t^{\prime})}|n\rangle$, so that the eigenvalue of the Floquet operator $U$ is given by

[TABLE]

We will now develop a perturbation theory to first order in $V$. We assume that the $n$-th eigenstate can be written as

[TABLE]

where $c_{n}(t)\simeq 1$ for all $t$, while $c_{m}(t)$ is of order $V$ for all $m\neq n$ and all $t$. Eq. (S18) implies that

[TABLE]

where the dot over $c_{m}$ denotes $d/dt$. Taking the inner product of Eq. (S22) with $\langle n|$ and using $\langle n|V|n\rangle=0$, we find that $\dot{c}_{n}=0$. We can therefore choose $c_{n}(t)~{}=~{}1$ for all $t$. We thus have

[TABLE]

Hence Eq. (S19) implies that the Floquet eigenvalue is still given by Eq. (S20) up to first order in $V$.

Next, taking the inner product of Eq. (S22) with $\langle m|$, where $m\neq n$, and integrating from $t=0$ to $T$, we get

[TABLE]

Since we know that Eq. (S23) satisfies

[TABLE]

we must have

[TABLE]

for all $m\neq n$. Eqs. (S24-S26) imply that we must choose

[TABLE]

We see that $c_{m}(t)$ is indeed of order $V$ provided that the denominator on the right hand side of Eq. (S27) does not vanish; we will call this case non-degenerate.

The above analysis breaks down if

[TABLE]

for a pair of states $|m\rangle$ and $|n\rangle$. We then have to develop a degenerate perturbation theory. Suppose that there are $p$ states $|m\rangle$ ($m=1,2,\cdots,p$) which have energy eigenvalues $E_{m}(t)$ satisfying Eq. (S28) for every pair of states. Ignoring all the other states of the system for the moment, we will assume that a solution of the Schrödinger equation is given by

[TABLE]

where we now allow all the $c_{m}(t)$’s to be order 1. Then we again obtain an equation like Eq. (S22) except that the sum over $m$ only goes over $p$ states. To first order in $V$, we can replace $c_{m}(t)$ by the time-independent constants $c_{m}(0)$ on the right hand side of Eq. (S22). Upon integrating from $t=0$ to $T$, this gives

[TABLE]

This can be written as a matrix equation

[TABLE]

where $c(t)$ denotes the column $(c_{1}(t),c_{2}(t),\cdots,c_{p}(t))^{T}$ (where the superscript $T$ denotes transpose), $I$ is the $p$-dimensional identity matrix, and $M$ is a $p$-dimensional Hermitian matrix with matrix elements given by

[TABLE]

Let the eigenvalues of $M$ be given by $\mu_{n}$ ($n=1,2,\cdots,p$). To first order in $V$, $I-iM$ is a unitary matrix and will therefore have eigenvalues of the form $e^{-i\mu_{n}}$; the corresponding eigenstates satisfy $c(T)=e^{-i\mu_{n}}~{}c(0)$. Next, we want the wave function in Eq. (S29) to satisfy $|\psi(T)\rangle=e^{-i\theta_{n}}|\psi(0)\rangle$. This implies that the Floquet eigenvalues are related to the eigenvalues of $M$ as

[TABLE]

Given a Floquet operator $U$, we can define a Floquet Hamiltonian $H_{F}$ as $U=e^{-iH_{F}T/\hbar}$. Comparing this with Eqs. (S31) and (S32), we see that the matrix elements of $H_{F}$ are given by

[TABLE]

We will now apply the above analysis to our model, where

[TABLE]

with $\lambda(t)=-\lambda$ for $0<t<T/2$ and $+\lambda$ for $T/2<t<T$. We will do Floquet perturbation theory assuming that $w\ll\lambda$. To do this, we consider states in the $\sigma_{l}^{z}$ basis. According to the Hamiltonian $H_{0}(t)$ in Eq. (S35), all such states $|n\rangle$ have an instantaneous energy eigenvalue $E_{n}(t)=(\lambda(t)/2)\sum_{l}\sigma^{z}_{l}$, which implies that the $\int_{0}^{T}dtE_{n}(t)=0$. Thus the unperturbed Floquet eigenvalue $e^{-i\theta_{n}}$ is equal to 1 for all states; we therefore have to do degenerate perturbation theory.

If $|m\rangle$ and $|n\rangle$ are two states which are connected by the perturbation $V$ in Eq. (S35), they differ by the value of $\sigma_{l}^{z}$ at only one site and therefore $E_{m}(t)-E_{n}(t)=\lambda(t)$, assuming that $|m\rangle$ and $|n\rangle$ have spin-up and spin-down respectively at that site. We find that the integral in Eq. (S30) is given by

[TABLE]

We therefore see that if

[TABLE]

where $q$ is an integer, then the expression in Eq. (S36) vanishes. This means that even in degenerate perturbation theory, there is no change in the Floquet eigenvalues and they remain equal to 1.

We will now use Eqs. (S34) and (S36). If $|m\rangle$ and $|n\rangle$ are two states which are connected by the perturbation $V$, we have $\langle m|V|n\rangle=-w$. We then obtain

[TABLE]

where $\gamma=\pi\lambda/(2\hbar\omega)$. The Floquet Hamiltonian is therefore given by

[TABLE]

This vanishes if Eq. (S37) is satisfied; we will then have to go to higher order in perturbation theory.

We note that Eq. (S39) can also be obtained by a straightforward expansion of the evolution operators $U_{\pm}$. To see this, we first note that for any two different sites $j$ and $j^{\prime}$ we have

[TABLE]

Thus as long as we are interested in terms of ${\mathcal{O}}(w)$, we can write

[TABLE]

One can then carry out a straightforward expansion of $U_{\pm}$. A few lines of algebra, required to gather terms of ${\mathcal{O}}(w/\lambda)$, yield

[TABLE]

Using Eq. (S42), one can compute $U=U_{+}U_{-}$. A re-exponentiation of terms of ${\mathcal{O}}(w/\lambda)$ then yields Eq. (S39) in a straightforward manner.

Next, we note that the magnitude of the right hand side of Eq. (S38) is independent of $|m\rangle$ and $|n\rangle$. Hence

[TABLE]

where the sum runs over all all pairs of states $m,n\in{\cal N}_{0}$ which are connected by $V$.

Given a system of size $L$ and periodic boundary condition ($PBC$), and the constraint that two up-spins cannot be on neighboring sites in any state, we can find the number of states and the value of $\sum_{mn}1$ in Eq. (S43). To this end, we first define the Fibonacci numbers which satisfy $F_{n}+F_{n+1}=F_{n+2}$, with $F_{1}=F_{2}=1$; as $n$ increases, $F_{n}$ quickly approaches $\tau^{n}/\sqrt{5}$, where $\tau=(\sqrt{5}+1)/2$ is the golden ratio. Keeping the up-spin constraint in mind, we define the transfer matrix

[TABLE]

where the indices $(ij)$ of $A_{ij}$ can take values 1 (spin-up) and 2 (spin-down). The number of states in an $L$-site system is then given by ${\rm Tr}(A^{L})=F_{L-1}+F_{L+1}$. To calculate $\sum_{mn}1$, we note that a spin at, say, site 2, can flip between up and down only if the spins at sites 1 and 3 are both down. The contribution of all such states to $\sum_{mn}1$ is 2 times the number of all possible states for sites 4 to $L$ with open boundary condition ($OBC$) and the up-spin constraint; the factor of 2 is because the spin at site 2 can be up or down, giving two states. Since the number of bonds for an $OBC$ system with sites 4 to $L$ is $L-4$, the number of states for such a system is given by the sum of all the matrix elements of $A^{L-4}$. This is equal to $F_{L-1}$. In the above argument, we assumed that there the spin which can flip between up and down is at site 2. However, the site 2 could have been anywhere else in the $L$-site system. We therefore see that $N_{0}\equiv\sum_{mn}1=2LF_{L-1}$. This leads us to define the normalized quantity

[TABLE]

V Symmetry of Floquet operator and zero modes

We will now discuss an exact symmetry of the Floquet operator for our driving protocal. We will then see that this symmetry implies that there will be a large number of states with zero quasienergy.

We define an operator

[TABLE]

which is unitary and satisfies $Q^{-1}=Q$. The eigenvalues of $Q$ are $\pm 1$, and the corresponding eigenstates have an even (odd) number of down spins. Next, we see that $Q$ anticommutes with the first term and commutes with the second term in Eq. (S6). As a result, the Floquet operator defined in Eq. (S7) satisfies

[TABLE]

This means that if $|\psi_{n}\rangle$ is an eigenstate of $U$ with eigenvalue $e^{-i\theta_{n}}$, $Q|\psi_{n}\rangle$ will be an eigenstate of $U$ with eigenvalue $e^{i\theta_{n}}$. Hence, all the quasienergies must come in $\pm$ pairs, except for quasienergies 0 and $\pi$ which correspond to Floquet eigenvalues equal to $\pm 1$. We also see that eigenstates of $U$ with eigenvalues $\pm 1$ can be simultaneously chosen to be eigenstates of $Q$; hence they will be superpositions of states all of which have an even or an odd number of down spins.

Given the Floquet Hamiltonian $H_{F}$ defined through $U=e^{-iH_{F}T/\hbar}$, Eq. (S47) implies that

[TABLE]

Note that this is an exact symmetry, independent of the Magnus expansion or Floquet perturbation theory. We will now use the arguments given in Ref. scarref1p5, to argue that there will be a large number of eigenstates of $H_{F}$ with zero eigenvalue; we will call these zero modes. Eq. (S48) implies that $H_{F}$ can be thought of as defining a tight-binding model of a particle moving on a bipartite lattice, where the two sublattices correspond to eigenvalues of $Q$ being equal to $\pm 1$. For such a tight-binding model, it is known that a lower bound on the number of zero modes is given by the difference of the number of states with $Q$ equal to $\pm 1$.

Next, we use the parity symmetry of our system corresponding to a reflection about the middle of the bond connecting sites labeled $L/2$ and $L/2+1$ (we will assume that $L$ is even). We define a parity transformation $P$ which does this reflection. Given an arbitrary state $|\psi\rangle$, the parity transformed state is $P|\psi\rangle$. The superpositions $|\psi\rangle\pm P|\psi\rangle$ then give two states with even (odd) parity respectively. However, states of the product form

[TABLE]

where each $\sigma_{n}$ can be spin-up or down, clearly satisfy $P|\psi\rangle=|\psi\rangle$. Such states therefore lie in the even parity sector; further, they have eigenvalue of $Q$ equal to $+1$ since each value of $\sigma_{i}$ appears twice in Eq. (S49). Let the number of such states be $N^{\prime}$ (we will calculate this number below). Now, let $N_{PQ}$ denote the number of orthonormal states with $P=\pm$ (denoting even/odd) and $Q=\pm$ (denoting $\pm 1$). We then see that $N_{++}=N_{-+}+N^{\prime}$ and $N_{+-}=N_{--}$. A lower bound on the number of zero modes in the odd parity sector is $|N_{-+}-N_{--}|$, and in the even parity sector is $|N_{++}-N_{+-}|=|N_{-+}-N_{--}+N^{\prime}|$. Since

[TABLE]

regardless of the values of $N_{-+}$ and $N_{--}$, we see that $N^{\prime}$ gives a lower bound on the number of zero modes.

To calculate the number $N^{\prime}$, we note that the string $(\sigma_{1},\sigma_{2},\cdots,\sigma_{L/2})$ for the first $L/2$ sites in Eq. (S49) must begin and end with a down spin to ensure that the neighboring sites $(L/2,L/2+1)$ and $(L,1)$ do not both have spin-up. The number of of such strings is given by the $22$ matrix element of $A^{L/2-1}$; this gives $N^{\prime}=F_{L/2}$. Thus the number of zero modes increases exponentially with $L$, as $\tau^{L/2}/\sqrt{5}$.

VI Half-chain entanglement

Here we detail out the procedure for computation of the half-chain entanglement $S_{L/2}$ used in the main text. The procedure could be applied to equilibrium or Floquet eigenstates or to an arbitrary driven state of the given model in Eq. S6.

We first note that the full density matrix ($DM$) is given by $\rho_{AB}=|{\psi}\rangle\langle{\psi}|$ where $AB$ is the whole system of size $L$ with $PBC$. We divide this system into two parts $A$ and $B$ of size $L/2$ each with $OBC$ as schematically shown in Fig. S1.

Next, we calculate the reduced density matrix of any of the subsystems by integrating out the other subsystem. This leads to

[TABLE]

This procedure involves taking a partial trace over the environment degrees of freedom. For example, the $(ij)$-th element, $\rho_{A}(i,j)$, of the reduced density matrix is given by

[TABLE]

where $i$ and $j$ represent product states in $A$, and $k$ represent product states in $B$ with OBC. However, since the full system ($AB$) had PBC, the Hilbert space dimension (HSD) of $AB$ does not match the HSD of ${\mathcal{H}}_{A}\otimes{\mathcal{H}}_{B}$ with OBC. To see this, we note that $HSD_{L}^{PBC}=F_{L-1}+F_{L+1}$ and $HSD_{L}^{OBC}=F_{L+2}$, where $F_{L}$ denotes the $L$-th Fibonacci number. Since $F_{L+2}=F_{L}+F_{L+1}>F_{L-1}+F_{L+1}$ for any $L$, one has $F_{L_{1}+2}F_{L_{2}+2}>F_{L_{1}+L_{2}-1}+F_{L_{1}+L_{2}+1}$ for any $L_{1},\,L_{2}$. Thus, while taking the summation in Eq. (S51), one has to exclude the states for which at any one of the junctions marked in red in Fig. S1, the end points of both $A$ and $B$ are occupied by a dipole.

Since the Hamiltonian and hence the full density matrix is in the configuration (product) basis, any matrix element of the reduced density matrix can be expressed as a sum of some of the elements of full DM, and we can rewrite Eq. (S52) as

[TABLE]

We stress here that one needs to be careful about this procedure if the density matrix is expressed using some other basis. For example in our case $\rho_{AB}=\sum_{m,n}\rho(m,n)|{m}\rangle\langle{n}|$ where $|{m}\rangle$,$|{n}\rangle$ are states in the zero total momentum ($K=0$) and even parity ($P=+1$) sector. In this case

[TABLE]

This indicates that one needs to search for states $|{m^{i;k}}\rangle$, $|{n^{j;k}}\rangle$ having non-zero overlap with the states $|{i;k}\rangle$ and $|{j;k}\rangle$ given by

[TABLE]

In this case, $(i,j)$-th element of $\rho_{A}$ is given by

[TABLE]

and the diagonalization of $\rho_{A}$ gives the Von-Neumann entropy since $S_{A}=-\sum_{i=1}^{HSD_{A}^{OBC}}p_{i}\ln(p_{i})$, where $p_{i}$ are the eigenvalues of $\rho_{A}$. This procedure is used to compute $S_{L/2}$ in the main text.

VII Fidelity and Level statistics

In this section, we note that the transitions from the ergodic to non-ergodic behaviors will also manifest themselves in the fidelity and eigenvalue statistics.

The fidelity of the driven Rydberg chain computed after $n$ cycles of the drive is defined as

[TABLE]

where $|\Psi(n)\rangle$ denotes the state of the system after $n$ cycles of the drive, and $|\psi_{0}\rangle$ is the initial state. In the regime where the dynamics is dominated by scars, the presence of long-term coherent oscillations indicates that one could expect periodic revival of ${\mathcal{F}}(n)$ to values near unity; in contrast, the thermal region, ${\mathcal{F}}(n)$ is expected to decay rapidly to zero and never revive. These behaviors are numerically confirmed in Fig. S2 near the ergodic to non-ergodic transition around $\omega_{D}=7.75$. The top left ($\omega_{D}=8.25$) and the bottom right ($\omega_{D}=7.5$) panels display periodic persistent revivals of ${\mathcal{F}}(n)$ before and after the first transition. Such revivals are completely absent in the bottom left panel ($\omega_{D}=7.75$) in the thermal region where ${\mathcal{F}}(n)$ decays to zero without any subsequent revival. The top right panel ($\omega_{D}=8$) shows an intermediate behavior displaying a few (smaller) revivals and subsequent decay. This indicates a crossover from a coherent to a thermal regime.

Next we discuss the characteristics of eigenvalue statistics across the transition. To this end, we first arrange the Floquet eigenvalues ($E_{F}^{n}$) (excluding the zero modes) in ascending order in the range $[-\omega_{D},\omega_{D})$, and then calculate the gaps $\delta_{n}=E_{F}^{n}-E_{F}^{n-1}$ between adjacent eigenvalues. This allows us to compute the distribution of the ratio of successive gaps in the energy spectrum huse

[TABLE]

The distribution of $r_{n}$ for an ergodic (thermal) system obeys the Gaussian orthogonal ensemble (GOE) and can be computed using random matrix theory Atas to be

[TABLE]

with an average value $\langle r\rangle_{GOE}\approx 0.535$. In contrast, for a fully non-ergodic (localized) system, the distribution of $r_{n}$ is Poissonian $P_{POI}(r)=2/(1+r)^{2}$ with $\langle r\rangle_{POI}\approx 0.386$. In Fig. S3, we plot $\langle r\rangle$ vs $\omega_{D}$ computed using the eigenvalues of $H_{F}$. The plot indicates that $\langle r\rangle$ reaches its GOE form precisely at the transition points; it remains lower than this value for other $\omega_{D}$. It is to be noted that $\langle r\rangle$ never reaches its Poisson value which indicates non-integrability of the system for all $\omega_{D}$. This confirms the signature of the transition in the level statistics.