Dynamical localization and slow thermalization in a class of disorder-free periodically driven one-dimensional interacting systems

We study if the interplay between dynamical localization and interactions in periodically driven quantum systems can give rise to anomalous thermalization behavior. Specifically, we consider one-dimensional models with interacting spinless fermions with nearest-neighbor hopping and density-density interactions, and a periodically driven on-site potential with spatial periodicity $m=2$ and $m=4$. At a dynamical localization point, these models evade thermalization either due to the presence of an extensive number of conserved quantities (for weak interactions) or due to the kinetic constraints caused by drive-induced resonances (for strong interactions). Our models therefore illustrate interesting mechanisms for generating constrained dynamics in Floquet systems which are difficult to realize in an undriven system.

Periodic driving of quantum systems can give rise to a host of interesting phenomena which have no equilibrium counterparts, such as the generation of drive-induced topological phases [49][50][51][52][53][54], Floquet time crystals [55][56][57], dynamical localization [58][59][60][61][62][63][64], dynamical freezing [65][66][67][68][69][70], tuning between ergodic and non-ergodic behaviors [71][72][73], and dynamical transitions [74][75][76][77][78][79].The outof-equilibrium dynamics of a wide class of closes quantum systems is believed to be governed by the eigenstate thermalization hypothesis (ETH) [80][81][82][83][84].According to ETH, all the eigenstates near the middle of the energy spectrum of a closed, non-integrable and disorder-free quantum system are thermal; the thermal nature of such states guarantees the ergodicity of the system.However, some instances are known where ETH is violated, for example, in integrable quantum systems and in manybody localized phases in one dimension in the presence of disorder and interactions [85,86].In recent years, it has been found that ETH can be broken in some quantum systems which are not integrable and have no disorder.The breaking of ETH may be weak or strong.In the case of weak ergodicity breaking, systems evade ergodicity due to the presence of quantum many-body scars [87][88][89][90][91][92][93][94][95][96][97][98].Quantum many-body scars are states which lie near the middle of the spectrum and have anomalously low entanglement entropy between two halves of the system.The number of scar states is typically much smaller than the full Hilbert space dimension.Moreover, the scar states form a subspace which is almost decoupled from the thermal subspace.Hence they are protected from thermalization for a long period of time and show persistent long-time coherent oscillations in their dynamics; this has been observed recently in Rydberg atoms simulators.Furthermore, the interplay between quantum manybody scar states and periodic driving can generate rich dynamical phase diagrams which have been studied in a number of papers [73,[99][100][101][102][103][104].One possible mechanism for strong ergodicity breaking is Hilbert space fragmentation [105][106][107][108] which occurs due to the presence of certain kinetic constraints in the dynamics.These kinetic constraints lead to the fragmentation of the Hilbert space into many disconnected sectors which can give rise to non-ergodic behavior in such systems.The Hilbert space fragmentation in Floquet systems has been examined recently [109].
It has been shown in a series of theoretical works that quantum many-body scar states can appear in systems hosting flat bands supported by compact localization [94][95][96].Motivated by this idea, we will pose a similar question in the context of Floquet systems undergoing dynamical localization (DL) [74][75][76][77][78][79].DL can be achieved in Floquet systems by tuning some of the system parameters which makes the effective hopping amplitudes zero or very small.DL can thus be a powerful tool for generating flat bands.The effects of interactions then become predominant which makes such systems highly promising platforms for investigating correlated out-of-equilibrium phases of matter.DL is special to periodically driven systems since this phenomenon has no equilibrium analog.Furthermore, the interplay between DL and quasiperiodic driving can prevent thermalization in a quantum system; this has been studied recently [110].Having understood the rich possibilities that DL can offer in Floquet systems, the natural question that motivates our work is as follows.Can the interplay between DL and interactions in a one-dimensional disorder-free periodically driven closed quantum system induce anomalous thermalization behavior?Model of Hilbert space fragmentation but different from the period-2 case Table 1: Schematic of main results obtained for the period-2 and period-4 models.µ and ω denote the driving amplitude and frequency respectively.
Another motivation for our work is as follows.As mentioned above, it has been known for many years that periodic driving of non-interacting systems can be used to generate quantum systems with a wide variety of band structures.It is therefore natural to ask if periodic driving of interacting systems can produce new kinds of quantum many-body systems whose parameters can be readily tuned.
The plan of this paper is as follows.In Sec. 2, we introduce our general model which consists of a one-dimensional system of spinless fermions with nearest-neighbor hopping, an on-site potential which is periodic in space [111] and is also driven periodically in time, and a density-density interaction between nearest-neighbor sites.In Sec. 3, we study in detail a model in which the potential has a periodicity of 2 sites.We first use first-order Floquet perturbation theory to derive an effective Floquet Hamiltonian for a non-interacting system when the driving amplitude and frequency are much larger than the hopping amplitude.We find that the system shows DL for certain values of the system parameters.We look at a two-point correlation function as a function of the stroboscopic time t = nT , where T = 2π/ω is the time period.We find that the correlator can decay as a power-law, where the power depends on the structure of the quasienergy dispersion around zero momentum.An interesting dynamical phase transition is found to occur when the dispersion changes, and a crossover between different powers occurs.Next, we look at the effects of interactions on DL.Exactly at a DL point, we find that there is an emergent integrability with a large number of conserved quantities; these quantities are just the fermion occupation numbers (0 or 1) at the different sites.An exact numerical calculation of the time evolution of the system confirms this result from the Floquet Hamiltonian, namely, we find that the two-point correlation function and Loschmidt echo (which is the overlap between an initial state and its time-evolved state) oscillate in time and show almost perfect revivals with a frequency which depends on the interaction strength.The integrability disappears as we go away from a DL point, and the Loschmidt echo then decays rapidly with time.In Sec.3.4, we consider the combined effects of resonances [112] and DL for the period-2 model.The first-order Floquet Hamiltonian then has a remarkable structure in which we have both the large number of conserved quantities as well as a density-dependent hopping in which the hopping between sites j + 1 and j + 2 depends on the occupation numbers on their neighboring sites j and j + 3. We thus obtain dynamical constraints on the hopping.This leads to the appearance of an exponentially large number of zero quasienergy states (an expression for this number is derived in Appendix A using a transfer matrix method), a highly fragmented Hilbert space, and several states with low entanglement entropy which lie near the middle of the quasienergy spectrum.
In Sec. 4, we study a model in which the potential has a periodicity of 4 sites.The potential has an amplitude which is driven periodically in time and a phase φ.The system has a mirror symmetry for two values of the phase, 0 and 7π/4.At a DL point, a period-4 model with φ = 0 behaves similarly as the period-2 model.But a period-4 model with φ = 7π/4 exhibits a different and remarkably rich set of behaviors.First, in the absence of interactions, the Floquet Hamiltonian has the form of the Su-Schrieffer-Heeger (SSH) model in which the nearest-neighbor hoppings have a staggered structure; this leads to the appearance of modes near the ends of an open system.At a DL point, we obtain an extreme limit of the SSH model in which hoppings on alternate bonds vanish.This leads to a large number of conserved quantities which is the total fermion occupation number on two sites between which the hopping is non-zero.Labeling the unit cell with two such sites as j, the conserved occupation number n j can take the values 0, 1 and 2. It is convenient to map the two possible states of a unit cell with n j = 1 to the states of a spin-1/2 object.We then discover that when interactions are introduced, a set of consecutive unit cells all of which have n j = 1 is described by the transverse field Ising model, in which neighboring spin-1/2's have σ x j σ x j+1 interactions and there is a transverse magnetic field term σ z j .In addition, the two boundary sites of this model have a longitudinal magnetic field term σ x j .The exact spectrum of this model can be found by mapping the spin-1/2 model to a model of fermions using the Jordan-Wigner transformation.Once again we examine the spectrum of the entanglement entropy versus the quasienergy and the time evolution of the Loschmidt echo.We find a clear fragmentation of the Hilbert space in terms of the quasienergy spectrum, and the Loschmidt echo shows oscillations for a long period of time.Both of these are consequences of the conserved quantities.We then study the effects of a staggered on-site potential at a DL point; we find that the fragments of the Hilbert space further break up into secondary fragments.Finally, we study what happens in the period-4 model when both resonances and DL are simultaneously present.The Floquet Hamiltonian again consists of a density-dependent nearest-neighbor hopping.We summarize our main results in Sec. 5.

Hamiltonian of period-m model
In this paper we will discuss a class of periodically driven one-dimensional models of interacting spinless fermions.The general form of the Hamiltonian is where J is a uniform time-independent nearest-neighbor hopping, µ(t) is the strength of an on-site potential which varies in space with period m (we will call this a period-m model), φ is a phase whose significance will be discussed later, V is the strength of a nearest-neighbor density-density interaction, n j = c † j c j , and N denotes the number of sites (we will use periodic boundary conditions unless otherwise specified).We will take µ(t) to be a periodic function of time with a time period T and the form and f (t+T ) = f (t).Since H(t) = H(t+T ), we will use the Floquet formalism to examine this model.The model can be analytically studied by performing a Floquet-Magnus expansion [1][2][3][4][5][6][7][8][9][10][11][12][13] which is valid in the large ω limit, where ω = 2π/T is the driving frequency.However, we will analytically study the model by finding an effective Floquet Hamiltonian H F using Floquet perturbation theory (FPT) [11,73,[113][114][115]; this approach is valid when both ω and µ are much larger than all the other parameters of the model.We will examine in detail two classes of models, period-2 and period-4, both analytically using the Hamiltonian H F and numerically by computing the Floquet operator U which evolves the system through one time period.We will set ħ h = 1 in this paper.

Period-2 model
In this section we will consider the period-2 model, whose form can be obtained by putting m = 2 and setting φ = 0 in Eq. ( 1), This Hamiltonian can be written in the following form in the language of a unit cell with two sites, where N /2 denotes the number of unit cells (we assume N is even), and each unit cell consists of two sites labeled a and b with a † j and b † j denoting the creation operator for a particle at odd-and even-numbered sites, respectively.We will assume periodic boundary conditions.
We will first study the non-interacting case with V = 0.The Hamiltonian in momentum space is then given by where k takes N /2 equally spaced values lying in the range of (−π/2, π/2].For each value of k, we define the Floquet operator where T denotes time ordering.The Floquet operator can be re-written in the form where H F k is the time-independent effective Floquet Hamiltonian.Assuming µ >> J, we can write We will now calculate H (1) F k using FPT to first order in the hopping J.We see from Eq. ( 2) that the two instantaneous eigenvalues of H 0 given by E ± k (t) = ±µ(t) satisfy the condition Hence we need to use degenerate FPT [11,73,113,115].
The eigenfunctions corresponding to E ± k are given by To construct the first-order Floquet Hamiltonian, we start with the Schrödinger equation where we assume that |ψ(t)〉 has the form and |n〉 = |±〉 k in our case.Using Eq. ( 11) and keeping terms up to first order in H 1 , we find that This can be re-written as follows where I denotes the identity matrix and H (1) F k refers to the first-order effective Hamiltonian.We then find that where To first order in H 1 , therefore, H F is given by Before discussing further, we note a symmetry of the Floquet operator, which follows from the expression in Eq. ( 7) and the fact that the driving protocol satisfies µ(T −t) = −µ(t).
Eq. ( 18) implies [77] that Hence H F can only have terms with odd powers of J.This implies that the next order term after the first order will be third order since there cannot be a term of second-order in J. Hence, the first-order effective Hamiltonian will be a very good approximation to the exact Hamiltonian in the limit µ J.
We also note that the Floquet quasienergy E F k must be an even function of k if we hold J, µ fixed.To prove this, we do the unitary transformation Thus we obtain With this Hamiltonian, it is clear that Hence the Floquet quasienergy must be an even function of k.

Dynamical localization for a single-particle system
It is evident from the form of the Hamiltonian obtained by first-order FPT that the system will (approximately) exhibit DL [58][59][60][61][62][63][64] when where n is a non-zero integer.We note that this condition for DL becomes more and more exact as µ/J → ∞ and the higher order corrections to the first-order effective Hamiltonian become negligible.In this limit, H F k vanishes for all values of k which produces a flat band with zero quasienergy.We can see in Fig. 1 (a) and (b), that for a system with J = 1, µ = 20, and ω = 20 and ω = 10, respectively, the quasienergy band is almost flat with a bandwidth ∆ ∼ 0.02.In both cases, we have taken µ >> J, and the higher order corrections to the first-order effective Hamiltonian are therefore very small.However, if we decrease the value of µ holding J fixed, the first-order Floquet Hamiltonian becomes less and lass accurate, and the DL begins to fail as can be seen in Fig. 1 (c), for a system with J = 1, and µ = ω = 10.In this case, the bandwidth, ∆ ∼ 0.08, is significantly larger than in the first two cases (∆ ∼ 0.02).

Dynamical phase transition
Motivated by our previous work [78], we now study the relaxation behavior of some correlators for the non-interacting model (V = 0).We will examine the correlation function a † j b j (where j denotes a particular unit cell) at stroboscopic times t = nT , where |Ψ 0 〉 is an initial state.For simplicity, we will take |Ψ 0 〉 to be a product state in momentum space with the form Since this state is translation invariant, the correlator C n will not depend on the unit cell index j.We will now investigate the relaxation behavior of C n for large values of n, particularly to see if there is any crossover behavior.We observe numerically that generally in the µ >> J limit (with a few exceptions), the correlation function exhibits a n −1/2 decay with oscillations and there is no crossover behavior.To explain this, we first note that the Floquet quasienergy obtained from the first-order Floquet Hamiltonian has the form where The stationary point of E k within the first Brillouin zone lies at k = 0. Next, we can show that the time-dependent part of the correlation function can be expressed as where For N → ∞ and assuming f (k) to be real, the above equation takes the integral form For large values of n, this integral gets a dominant contribution from the regions close to the stationary points of E k [75,77,78].Therefore, expanding Eq. ( 28) around the stationary point at k = 0, we find that where ζ = 4J T (sin A/A).Assuming f (k) = 0 for a generic initial state, we find that the correlation function for large values of n will decay as a power n −1/2 , and there will be oscillations due to the term cos(ζn).This is the usual decay behavior unless there are competing terms which come from higher-order corrections.
We encounter such higher-order terms slightly away from the DL points (we recall that a DL point is where A is an integer multiple of π, i.e., µ/ω is an integer).Since an analytical calculation of the third-order effective Hamiltonian is a tedious task, we analyze this regime numerically.We first find the Floquet quasienergy from the numerically exact calculation and then do a fitting analysis of it.The results we obtain from such an analysis are as follows.We consider the parameter values J = 1, µ = 10, and ω = 10.6.Taking into account the structure of the stationary point obtained from the first-order effective Hamiltonian, we fit the numerically computed Floquet quasienergy around k = 0 as a function of k up to sixth order in k.We find the following functional form All terms with odd powers of k are found to be zero which is expected by the symmetry discussed in Eq. (22).Further, the coefficients of k 2 and k 6 are also found to be zero at these particular parameter values.Hence Eq. (30) shows that E k goes as k 4 near k = 0.An analysis similar to the one following Eq.( 29) then shows that Assuming f (k) = 0, we see that for large values of n, the correlation function will decay as n −1/4 with oscillations due to the cos(2nT p 0 ) term.This is what we see in the Fig. 2 (c): for a system with the parameter values mentioned above, the correlation function function decays as an oscillatory term times n −1/4 .If we plot |δC n | (rather than δC n ) versus n, the period of oscillations ∆n will be given by the condition, 2∆np 0 T = π, which implies that ∆n = ω/(4p 0 ).Putting ω = 10.6, and p 0 = 0.07865, we find that ∆n ∼ 34), which agrees very well with the oscillation period seen in Fig. 2. Interestingly, we observe a crossover from n −1/4 to n −1/2 as we move slightly away from ω = 10.6.To see this crossover, we consider J = 1, µ = 10, and ω = 10.7 as an example.We again follow the same fitting procedure, and find the following functional form The terms with odd powers of k vanish as before.In contrast to the previous case, however, there is now a competition between the k 2 and k 4 terms, which can be seen in Eq. ( 32).Close to k = 0, we see that E k has a leading contribution coming from the k 4 term, followed by a subleading correction due to the k 2 term very close to k = 0. Expanding the integrand in Eq. ( 28) around k = 0 then gives Defining a scaled variable k = kn 1/4 , and assuming f (k) = 0, we find Since |p 2 | << |p 4 |, it is clear from Eq. ( 33) that there will be a n −1/4 scaling (along with oscillations) when |ε|n 1/2 << 1, where ε = p 2 /p 4 .However, when |ε|n 1/2 >> 1, the n −1/4 scaling breaks down and we then encounter a different scaling law, namely, n −1/2 (along with oscillations due to the cos(2nT p 0 ) term).We can extract the crossover scale n c from this analysis; a crossover between the n −1/4 and n −1/2 power-laws occurs when |ε|n 1/2 c ∼ 1, which implies that n c ∼ 1/|ε| 2 .

Effects of interactions on dynamical localization
In this section, we will look at the effects of density-density interactions [61,63], H I = j V n j n j+1 on DL.For this case, we take the system to be at half-filling with periodic or antiperiodic boundary conditions (c N +1 = c 1 or − c 1 ) depending on whether the particle number is even or odd, respectively, to avoid any degeneracy in the spectrum.To get an analytical understanding of this system, we will first compute the effective Hamiltonian up to first order in V .Since H I is diagonal in the position basis and commutes with the unperturbed Hamiltonian, H 0 (t), the effective Hamiltonian to first order in V simply reads as Hence the full effective Hamiltonian to first order in J and V is as follows We note that the Floquet evolution operator U(T ) satisfies the same condition as mentioned in Eq. ( 19) [77], and therefore H F possesses the symmetry Hence the first-order effective Hamiltonian will be a very good approximation to the exact Hamiltonian for µ >> J, V , since the higher order corrections will be negligible compared to the first-order term.
We see from Eq. ( 36) that H (1) F, J = 0 at the DL points where A is an integer multiple of π, and then H (1) F I .Consequently, the spectrum of the Floquet quasienergies becomes easy to compute at any filling due to the diagonal form of the effective Hamiltonian in the position basis.
We now consider a system at half-filling with antiperiodic boundary conditions for J = 1, µ = 20, ω = 20, V = 1, and L = 16, and calculate the spectrum of the Floquet quasienergies and half-chain entanglement entropy.The dimension of the Hilbert space for N = L/2 and L = 16 is 16 C 8 = 12870.We first discuss the Floquet quasienergies and the degeneracies for some of the Floquet eigenstates which can be obtained analytically.We observe that there are eight equally spaced quasienergies lying between 0 and 8V with a energy spacing of V for L = 16, as shown in Fig. 3 (a).Next, there are exactly two Floquet eigenstates, |±〉 = 1/ 2 |1 〉 ± |2 〉 with E F = 0, where |1 〉 and |2 〉 are equal to |1010101010101010〉 and |0101010101010101〉 respectively in the number basis.Furthermore, there are 16 Floquet eigenstates with E F = 8V .We note that the DL induces an emergent integrability which leads to the appearance of many flat bands and several low-entanglement states (with S L/2 << S pa g e , where S pa g e = (L/2) ln 2 − 1/2 [116,117]) near the middle of the spectrum as can be seen in Figs. 3 (a) and (b) respectively.However, this emergent integrability starts to break down as we move away from a DL point, as we see in Fig. 3 In Figs. 3 (c) and (d), we observe that the flatness of the bands begins to disappear as we move away from the DL limit.
In Fig. 4, we show some dynamical properties of the system, namely, the two-point correlation function and the Loschmidt echo at a DL point and away from a DL point.The parameter values chosen for this figure are the same as in Fig. 3.We will take the initial state to be the ground state of the undriven Hamiltonian.In Fig. 4 (a), we see almost perfect revivals of the correlation function in time, which is expected due to the emergent integrability at a DL point.Further, the oscillation period of the revivals can be calculated by using the condition e iV δt = 1.Hence the time period is given by δt = 2π/V , where V = 1 in our case.Away from the DL point, the correlation function decays rapidly in time due to a breakdown of the integrable structure.In Fig. 4 (c), we show the Loschmidt echo, |〈ψ(t)|ψ(0)〉|, for the same parameter values as in Figs. 4 (a) and 4 (b) and with the same initial states.For the first case with ω = 20, we again see perfect revivals with a period of 2π in Fig. 4 (c), upper panel.For the second case with ω = 22, the amplitude of Loschmidt echo decays rapidly as shown in Fig. 4 (c), lower panel.These results indicate that the system evades thermalization for a long time at a DL point but thermalizes quickly as we go away from a DL point [58][59][60][61][62][63][64].

Effects of resonances
In this section we will examine the effects of resonances [5,112] on DL for the period-2 model.To obtain an analytical insight for this case, we will consider parameter values with µ = ω = V >> J and derive an effective Floquet Hamiltonian using first-order FPT.We first consider a four-site system so that we can easily identify various non-trivial processes, and we will then generalize it to larger system sizes.We find that a four-site system only offers four distinct non-trivial processes for a particular choice of a periodic potential pattern due to the constraints imposed by DL.We further note that there are two possible potential patterns available for such a system consisting of four sites.Therefore, we need to consider a total of eight distinct non-trivial processes while formulating the first-order effective Floquet Hamiltonian.These eight processes and the corresponding time-dependent effective Hamiltonians are listed below.
Pattern of periodic potential Process Effective time-dependent Hamiltonian Table 2: Allowed processes and their corresponding effective time-dependent Hamiltonians for a four-site system with all possible patterns of periodic on-site potential in the case of resonance and dynamical localization for a period-2 model.
As an example, we consider the first process listed above and calculate the first-order Floquet Hamiltonian for this case.We first note that the Hamiltonian can be recast as where I, σ x and σ z denote the identity and two of the Pauli matrices, and H 0 and H 1 are the unperturbed Hamiltonian and perturbation, respectively.Assuming µ = ω and V >> J, the instantaneous eigenvalues of H 0 are given by E ± k = ±(µ(t) + V /2).The eigenfunctions corresponding to E ± k are given by |+〉 = 0 1 and |−〉 = 1 0 .These two eigenvalues satisfy the condition given in Eq. ( 9), and we therefore use degenerate FPT.This gives

〈+|H
(1) Putting µ = V = ω, we find that I(µ, V, ω) = −4i/(3π).Thus, the effective Hamiltonian for this particular process is where we have set J = 1.Following similar procedures, we can compute the effective Hamiltonians for all the other processes.These are given below.
Pattern of periodic potential Process First-order Floquet Hamiltonian + -+ -1100 ↔ 1010 H (1) Table 3: First-order effective FPT Hamiltonians for the allowed processes with all possible patterns of periodic potential in the case of resonance and dynamical localization for a period-2 model.
Taking all these processes into account, the complete effective Hamiltonian for the case where a resonance and DL occur simultaneously is given by We can perform the unitary transformation c 2 j → c 2 j and c 2 j+1 → ic 2 j+1 to obtain a simpler form of the effective Hamiltonian This form implies that hoppings between two nearest-neighbor sites are forbidden whenever their neighboring sites are both completely empty or completely occupied.Hence, these forbidden processes act as kinetic constraints in the dynamics [105][106][107][108][109], and these constraints can, in principle, lead the system towards an anomalous thermalization behavior.We also find that the Hamiltonian in Eq. ( 42) has several zero-energy states.The number of such states can be found using a transfer matrix method as shown in Appendix A. We discover that the number grows exponentially with system size as 1.466 L .In Figs. 5 (a) and (b), we show the variation of the half-chain entanglement entropy S L/2 as a function of E F , obtained from exact numerical calculations and from first-order FPT for J = 1, µ = ω = 20, and V = 20.Both these cases point towards many low-entanglement states near the middle of the spectrum; these arise due to the kinetic constraints simultaneously imposed by DL and the resonance condition.Before proceeding further, we note that the effective Hamiltonian described in Eq. ( 42) supports many fragmented Hilbert space sectors [105][106][107][108][109], which can be shown as follows.One of the simplest fragments consists of the four states, |0000111111110000〉, |0001011111110000〉, |0000111111101000〉, and |0001011111101000〉, and their translated partners (we have written all the states in the occupation number basis).The action of the effective Hamiltonian on these four states is schematically shown in Fig. 6.Taking into account the action of the effective Hamiltonian on these four states, we find an The eigenvalues of this Hamiltonian are given by E 1,2 = ±0.85,E 3 = 0, and E 4 = 0. Hence these four eigenvalues offer two distinct difference in energies, i.e., ∆E = 0.85 and ∆E = 1.7, which will be important later in the discussion of dynamics.In Figs. 7 (a) and (b), the variation of the Loschmidt echos with time, t = nT , is shown as found from the exact numerical calculation and the first-order FPT, respectively, for the parameter values, J = 1, µ = ω = V = 20, and L = 16, taking the initial state to be |0000111111110000〉.We can show analytically that the Loschmidt echo for this particular choice of initial state takes the form |a + b cos(∆E t)|, which implies that it oscillates with a period ∆t = 2π/∆E.Putting ∆E = 0.85, the period of oscillation in the revival pattern turns out to be ∆t ≈ 7.4, which almost perfectly captures the numerically obtained value.In Fig. 7 (c), we show the overlaps of the same initial state with the Floquet eigenstates (obtained from the exact numerical calculation) as a function of E F , where the color bar indicates the variation of S L/2 of the Floquet eigenstates.Interestingly, we observe that the overlap is highest for the Floquet eigenstates with E ± 0.85 and 0, which almost identically agrees with the analytically predicted values.
4 Period-4 model We will now discuss a model with the second type of periodic potential, namely, the period-4 model.The general form of the Hamiltonian for this class is given by where J denotes the nearest-neighbor hopping, µ is the amplitude of the periodic potential with m = 4, V defines the nearest-neighbor density-density interaction, and φ refers to a generalized phase.This  model possesses a mirror symmetry [118] for certain special values of φ.Since we are interested in mirror-symmetric configurations for our analysis, it turns out that we can only have two possible realizations, namely, φ = 0 and φ = 7π/4.We will denote these as type-1 and type-2 cases respectively.

Type-1 mirror-symmetric case
Taking µ(t) to be proportional to µ in Eq. ( 44), we find that the on-site potential pattern for φ = 0 is given by (µ, 0, −µ, 0) on four consecutive sites numbered (4n, 4n + 1, 4n + 2, 4n + 3); this is shown in Fig. 8 (a).Assuming µ >> J, and using the results presented in Appendix B, we find that the first-order FPT Hamiltonian is given by Since the interaction term commutes with the unperturbed Hamiltonian, H 0 , that part of the Hamiltonian will just be given by to first order in V .
The symmetry property discussed in Eq. ( 37) again holds for this model, and therefore, similar to the period-2 case, this system will also exhibit DL when M = 0, i.e., when µ = 2nω where n = 1, 2, 3, • • • .Thus, this mirror-symmetric configuration [118] of the period-4 model with φ = 0 and the period-2 model are identical to each other exactly at a DL point.In Figs. 9 (a) and 9 (b), we show the Floquet quasienergy spectrum E F and the variation of the entanglement entropy S L/2 with E F for J = 1, µ = 20, ω = 20, and V = 0.5.The spectrum looks almost identical to the spectrum of the period-2 model at DL with many low-entanglement states near the middle of the spectrum.

Type-2 mirror-symmetric case
The period-4 model with φ = 7π/4 is another mirror-symmetric configuration with many interesting properties, which we will now discuss in detail.Taking µ(t) to be proportional to µ 2 in Eq. ( 44), the on-site potential pattern for φ = 7π/4 is given by (µ, µ, −µ, −µ) on four consecutive sites numbered (4n, 4n + 1, 4n + 2, 4n + 3).Assuming µ >> J, we obtain the following first-order FPT Hamiltonian, Remarkably, we see that the non-interacting part of the Hamiltonian exactly describes the SSH model [119], with nearest-neighbor hoppings which have alternating strengths given by J and J|M 1 |.(It is clear that |M 1 | is always smaller than 1.The phase of M 1 can be removed by doing a unitary transformation of the form c j → c j e iα j with appropriately chosen α j 's).We thus see that the periodicity of the model has effectively reduced from 4 to 2.
The expression for M 1 implies that this model will exhibit DL for µ = nω, where n = 1, 2, • • • .Exactly at these points, the effective first-order Hamiltonian is given by The non-interacting part of this Hamiltonian is an extreme limit of the SSH model with alternating nearest-neighbor hoppings γ 1 = 1, and γ 2 = 0. Our model therefore inherits the property of the SSH model that a system with open boundary conditions has topologically protected zero-energy edge modes provided that the hopping strength on the leftmost or rightmost bond is weaker than the strength of the bond next to it.For φ = 7π/4, we find that the leftmost bond (between sites numbered 0 and 1) has a hopping strength which is larger than the strength of the next bond, and therefore, the system has no edge modes.However, the system has edge modes for φ = π/4, i.e., when the bonds are shifted by one unit cell and the stronger and weaker bonds get interchanged (see the schematic pictures in Fig. 10).In Figs.Since the interaction part again commutes with unperturbed Hamiltonian, H 0 , the effective Hamiltonian to first order in V again reads as H (1) F 2 = V n j n j+1 .The interplay between interaction and DL in this case gives rise to various intriguing phenomena.We will study this in the next few sections using an effective spin model based on the first-order Floquet Hamiltonian.

Effective spin model based on first-order effective Hamiltonian
To derive the effective spin model, it is convenient to recast the effective Hamiltonian in terms of unit cells where a j and b j denote the annihilation operators on the even and odd numbered sites of the j-th unit cell.(Henceforth we will refer to the j-th unit cell as the j-th site for convenience).The schematic of the periodic-4 model at a DL point is shown in Fig. 12.This above form suggests that the particle number n j at the j-th site (unit cell) commutes with H (1) F has L/2 conserved quantities.(Note that these will be only approximately conserved quantities.The exact effective Hamiltonian will have higher order terms which do not commute with these quantities).For a system consisting of two sites with n max j = 2, we can have nine possible effective Hamiltonians which are shown in the table below.
Table 4: Allowed configurations and the corresponding effective Hamiltonians for two unit cells at a DL point with µ >> V for a period-4 model.
In the table above, n 1 and n 2 are the occupation numbers of the first and second unit cells, respectively.We can see that eight out of the nine possibilities shown above can be mapped to a noninteracting problem.The only instance when the effect of the interaction is non-trivial is when both the unit cells are singly occupied.For this case, an effective spin degrees of freedom can be defined as where |10〉 defines a unit cell with the left and the right sites being occupied and empty, respectively, and |01〉 means the other way around.With this definition, the correlated two-site problem takes the following form where σ z and σ x are Pauli matrices.This two-site problem can now easily be generalized to larger system sizes.To do so, we first consider the case where all the sites are singly occupied.The effective spin Hamiltonian for this case is given by which is essentially the transverse field Ising model with the interaction term −V /4 and the transverse field J.The other four cases for a system with L/2 − 2 unit cells being singly occupied and with the two boundary unit cells being either empty or doubly occupied have effective spin Hamiltonians as follows.
In Table 5, n L and n R denote the occupation numbers of the leftmost and rightmost unit cells labeled j = 1 and L/2, respectively.Therefore, we see from Table 5 that the effective spin Hamiltonian has the form of the transverse field Ising model with additional longitudinal magnetic field terms of strength ±V /4 at the boundary sites, depending on the adjoining sites having n = 0 or 2. Before Table 5: The four possible effective spin Hamiltonians emerging from a system with L/2−2 unit cells being singly occupied and two boundary unit cells, each of them either completely occupied or completely empty, at a DL point with µ >> V for a period-4 model.
proceeding further, we perform the transformation σ x j → σ z j , σ z j → −σ x j , and σ y j remains unchanged.The Hamiltonian then takes the form It may appear that the longitudinal field terms at the boundary sites 2 and L/2 − 1 would make it difficult to find the energy spectrum analytically for this model.To overcome this problem, we add two more sites, labeled 1 and L/2, with Pauli operators σ x 1 and σ x L/2 , which couple to σ x 2 and σ x L/2−1 respectively [120].The Hamiltonian then becomes where we have ignored some constants.Note that σ x 1 and σ x L/2 commute with the Hamiltonian.Hence, there are four decoupled sectors of states corresponding to σ x 1 = ±1 and σ x L/2 = ±1.These four sectors precisely cover the four possible combinations of ± signs in Eq. ( 52).
The Hamiltonian in Eq. ( 53) can now be solved analytically by writing it in terms of Majorana fermion operators using the Jordan-Wigner transformation, where α j , β j are Majorana operators.In terms of these operators, the Hamiltonian takes the form Since this Hamiltonian is quadratic in terms of Majorana operators, it describes a non-interacting system and its spectrum can be found exactly [120].To examine the effects of DL on the thermalization of the system, we consider the variation of the half-chain entanglement entropy S L/2 with the quasienergy E F , which gives a static measure of ergodicity.As shown in Figs. 13 (a) and 13 (b), we consider the system at a DL point with J = 1, µ = 20, and ω = 20, and take V to be 0.5 and 2, respectively.In the first case, we observe many finger-like structures [121] in the entanglement spectrum, which are due to the presence of an extensive numbers of approximate conserved quantities arising due to the DL.Furthermore, the DL offers many frozen states with extremely low-entanglement values, i.e., S L/2 = ln 2 0.693 or 2 ln 2 1.386 near the middle of the spectrum.Some of the frozen states can be found easily from the effective Hamiltonian, such as |22220000〉, |20220200〉, |22202000〉 and their translated partners.Nevertheless, as shown in Fig. 13 (b), these finger-like structures are absent for V = 2 due to the disappearance of these approximate conserved quantities with increasing interaction strength.The low-entanglement states near the middle of the spectrum are still present, which again indicates that this system would thermalize very slowly.We can, therefore, conclude that this slow thermalization occurs due to two possible mechanisms: (i) the existence of extensive numbers of conserved quantities arising due to the DL, which grows exponentially with the system size as 3 L/2 (which grows less rapidly than the Hilbert space dimension which goes as 2 L ) [122].We emphasize again that these quantities are conserved to a good approximation only for µ >> J, V .(ii) the presence of many frozen state configurations, which do not participate in the dynamics at a DL point.In Fig. 13 (c), we consider a system away from a DL point with µ = 10, ω = 20 and V = 0.5, and we see that the low-entanglement states have disappeared, signaling that the system should thermalize quickly.To see a dynamical signature of slow thermalization, we study the time evolution of the Loschmidt echo precisely at a DL point with two non-trivial initial states, i.e., |1〉=|21012101〉 and |2〉=|21102110〉.As shown in Fig. 14 (a), the Loschmidt echo exhibits long-time revivals, which indicates that the system is evading a thermalized behavior.To gain an analytical understanding, we first derive an effective Hamiltonian for the state |210〉, H e f f = V J J 0 .This Hamiltonian has the energy eigenvalues . Taking this into account, we see that the initial state |1〉 would have the highest overlap with 4 × 2 4 = 64 such Floquet eigenstates with , and . From the numerically obtained data for a system with µ = ω = 20, and V = 0.5, we find that two Floquet eigenstates with E F = 4V = 2 and 4V + 2 J 2 + V 2 4 = 4.06 have the highest overlaps with this particular initial state, which agrees quite well with the analytically predicted values.Within the approximation of these two highest overlapping states, the Loschmidt echo will have a time-dependence of the form be |1+e which oscillates with a period given by ∆t = 2π/(E 1 − E 2 ).For the parameter values given above, we find ∆t ∼ 3, which agrees with what we see in Fig. 14 (a).In a similar manner, we can find the period of oscillation for the other three DL points with ω = 5, 6.67 and 10.For Fig. 14 (b), we choose the initial state |2〉 which has the highest overlaps with 64 such Floquet eigenstates.As shown in Fig. 14 (b), we see long-time oscillating behaviors in the Loschmidt echo at the four DL points, which again indicates that the system will not thermalize for a long time.In Figs. 14 (c) and (d), we plot the overlaps of these two initial states, |1〉 and |2〉, with all the Floquet eigenstates at a DL point with µ = ω = 20, where the color bar shows the variation of the entanglement entropy of the Floquet eigenstates.In both cases, we observe that these two initial states have overlaps with multiple Floquet eigenstates, some of them having extremely low entanglement entropy which possibly causes the long-time persistent oscillations in the Loschmidt echo.
To further confirm the above findings, we investigate the dynamics of S L/2 and the unequal-time two-point density-density correlation function with the initial state taken to be |21012101〉.As shown in Figs. 15 (a), we see that the entanglement entropy slowly increases before reaching a saturation value for µ = 2ω = 20 and V = 0.5, as expected.Further, the saturation value (∼ 3.5) is much less than S pa g e ∼ 5.1 signaling a deviation from a volume law behavior.In Fig. 15 (b), we examine the two-point correlation function with the same initial state and see a behavior similar to S L/2 .Thus, the dynamics also confirms the behaviors predicted from the static signatures.In Figs. 15 (c) and 15 (d), we repeat the same analysis for the system away from a DL point, namely, for µ = 20, ω = 11, and V = 0.5.As opposed to Fig. 15 (a), we see in Fig. 15 (c) that the entanglement entropy saturates quite quickly.In Fig. 15 (d), the two-point correlation function exhibits a similar behavior with the saturation value of 〈n j 〉 2 = 1/4 at half-filling, which suggests that thermalization has occurred.
Finally, we compare the results obtained from the exact numerics and the first-order FPT calculation for µ >> J, V .In Figs.16 (a) and 16 (b), we see that the quasienergies obtained from exact numerics and from first-order FPT for J = 1, µ = ω = 20 (lying at a DL point), and V = 0.5 agree very well with each other.However, the results for the entanglement entropy do not agree that well.We believe that this disagreement is due to the corrections in FPT which are higher than first-order.The correction terms have a relatively small effect on the Floquet quasienergies but have a noticeable effect on the Floquet eigenfunctions.Consequently, the entanglement entropy deviates significantly from the exact numerical values.Nevertheless, we note that the qualitative behavior of the entanglement entropy is the same in the two cases.In Figs.16 (c) and 16 (d), we compare the same plots as in Since the first-order term dominates over higher order corrections away from a DL point, we see that the Floquet quasienergy and entanglement entropy agree with each other almost identically.

Effects of staggered on-site potential
It is interesting to incorporate the effects of a staggered on-site potential with amplitude w in the period-4 model with φ = 7π/4 since such a potential commutes with the unperturbed Hamiltonian H 0 .Hence, in the presence of such a potential, the first-order FPT Hamiltonian becomes H (1) F,stagg = w L j=1 (−1) j n j .This term can also be incorporated within the effective spin model.Assuming that all the unit cells are singly occupied (so that there is no boundary field), we obtain (In general, the σ z term only appears for unit cells with single occupation, and has no effect on unit cells with n j = 0 or 2).The other sectors where only some of the unit cells have single occupation get modified in a similar way.In Fig. 17, we show the variation of the entanglement entropy S L/2 versus E F for a system with J = 1, µ = ω = 20, V = 0.5, and (a) w = 1 and (b) w = 3.In both cases the spectrum contains multiple finger-like primary structures with further secondary fragments [123], and the secondary fragments become more prominent with increasing value of the staggered potential.In Fig. 18, the dynamics of the Loschmidt echo is shown for a system with w = 3, and ω = 5, 6.67, 10 and 20, respectively.For all four DL points, the Loschmidt echo for the initial state |21012101〉 exhibits an oscillatory behavior for a long period of time, which implies that the system thermalizes very slowly.

Effects of resonances
We have so far discussed cases with V << µ where the system can evade thermalization at a DL point.In this section, we will address the effects of resonances at two DL points, i.e., µ = ω = V and µ = 2ω = V .For both cases, we note that µ, V >> J, which is the opposite limit to the previously examined cases, and we want to investigate whether this limit can give rise to non-ergodic behavior similar to the period-2 model.To do so, we will first derive an effective Hamiltonian based on the first-order FPT.Similar to the period-2 case, we will first identify the non-trivial processes for a system with four sites.Due to the periodic pattern of the on-site potential, we need to consider a total of sixteen independent processes to construct the time-dependent effective Hamiltonian.These are shown in Table 6.
As an example, we will derive the effective time-independent first-order FPT Hamiltonian for the first process shown in Table 6.In this case, the effective time-dependent Hamiltonian can be written Table 6: Allowed processes and their corresponding effective time-dependent Hamiltonians for a four-site system with all possible patterns of the periodic potential in a period-4 model in which both resonances and dynamical localization are present. as where H 0 and H 1 are the unperturbed Hamiltonian and the perturbation, respectively, and we will assume that µ, V >> J.The instantaneous eigenvalues of H 0 are E ± = ±(µ(t) + V /2).The eigenfunctions corresponding to E ± k are given by |+〉 = 1 0 and |−〉 = 0 1 .These two eigenvalues again satisfy the condition given in Eq. ( 9).Therefore, following the usual steps of degenerate FPT, we obtain

〈+|H
(1) Substituting µ = V = ω, I(µ, V, ω) turns out to be 4i/(3π).Thus, the effective Hamiltonian for a system consisting of four sites with this specific choice of periodic potential pattern becomes (setting J = 1) Following the same procedure, we can compute the effective Hamiltonian for all the other processes.
Pattern of periodic potential Process First-order Floquet Hamiltonian + + --1100 ↔ 1010 H (1) Table 7: First-order effective FPT Hamiltonians corresponding to the allowed correlated processes for a four-site system with all possible patterns of periodic potential in the case of resonance and dynamical localization for a period-4 model.
These results enable us to deduce the complete first-order effective Hamiltonian for µ = ω = V , namely, The form of this Hamiltonian suggests that certain nearest-neighbor hoppings are forbidden, as elaborated below, when the DL and resonance condition are simultaneously satisfied; in principle this can lead to an anomalous thermalization behavior.In Fig. 19, we show the variation of S L/2 with E F obtained from exact numerical calculations and from the first-order FPT analysis, for J = 1, and µ = ω = V = 20.As anticipated, the middle of the entanglement spectrum consists of many lowentanglement states possibly due a fragmented nature of the Hilbert space as described below.For example, we can see that this effective Hamiltonian supports a simple fragment which consists of only one state, |1100110011001100〉, and its translated partners.To construct this single fragment, we note the following constraints following from Eq. ( 60).The hopping on the bonds (2 j + 1, 2 j + 2) is only possible if the neighboring sites have (n 2 j , n 2 j+3 ) = (0, 1) or (1, 0).However, the hopping on the bonds (2 j+2, 2 j+3) requires the neighboring sites to have (n 2 j+1 , n 2 j+4 ) = (0, 0) or (1, 1).These two constraints enables us to show that the above state forms a fragment on its own, and it does not mix with other states due to the action of the Hamiltonian.In Fig. 20 (a) and (b), we study the dynamics of the Loschmidt echo for the initial state |1100110011001100〉 using the exact Floquet dynamics and the first-order effective Hamiltonian, respectively.In both cases, we find that the Loschmidt echo stays very close to 1, with some small oscillations in (a).In Fig. 20 (c), we see that these initial states have highest overlap with two mid-spectrum Floquet eigenstates with S L/2 = ln 2. Since S L/2 for this initial state is much smaller than S pa g e ( 5.1 for L = 16), this state is likely to retain its initial memory for a long period of time.Now we will consider the case µ = 2ω = V .Although both µ = ω and µ = 2ω give rise to DL for a non-interacting system, the presence of V makes a significant difference in the effective Hamiltonian description.Going through the same procedure as before, we obtain the effective Hamiltonian Note that there is no hopping on the bonds (n 2 j+1 , n 2 j+2 ).This implies that the occupation number n 2 j + n 2 j+1 in the j-th unit cell commutes with the effective Hamiltonian for all values of j.Hence there are L/2 approximately conserved quantities, which can protect some of the mid-spectrum states from thermalization for a long time.In Fig. 21, the entanglement entropy spectrum obtained by (a) exact numerics and (b) first-order FPT are shown for J = 1, µ = 2ω = V = 20, and L = 16.The fragmentation in the spectrum points towards the existence of conserved charges following from the first-order effective Hamiltonian.

Discussion
The central results presented in our paper are as follows.We have unraveled an intricate dynamical behavior of a class of disorder-free one-dimensional interacting spinless fermionic models with a  The color intensity indicates the density of states, revealing that the majority of states do not attain the thermal value of te entanglement entropy which is given by the upper envelope of the plots.
periodically driven on-site potential which is also periodic in space.In the absence of interactions, this class of models exhibits DL for a particular set of parameter values, giving rise to one or more flat bands.We have focused in detail on two models, corresponding to potentials with period-2 and period-4 on the lattice.For the period-2 model, we describe a dynamical phase transition which can be observed in the relaxation behavior of correlators in the absence of any interactions.Our investigation shows that a crossover behavior between different power laws of the decay of correlations generally occurs away from the DL points.We find that in the period-2 model, the flat bands which arise due to DL are stable in the presence of a comparatively weak interaction strength due to an emergent integrability.Further, the spectrum of half-chain entanglement entropy as a function of the Floquet quasienergy for a weakly interacting system reveals that there are many low-entanglement states near the middle of the quasienergy spectrum, implying that the system may evade thermalization for a long time.The persistent oscillations in the correlation functions and in the Loschmidt echo which survive for a long period of time support the above statement about thermalization.However, these oscillations decay rapidly in time when we move away from these fine-tuned parameter values.Remarkably, our model also exhibits Hilbert space fragmentation due to the presence of kinetic constraints when the DL and resonance condition are simultaneously satisfied and the interaction is strong.In the case of period-4, the behavior appears to be much more intriguing, even in the regime of weak interaction.The period-4 model possesses two mirror-symmetric configurations, corresponding to two values of the phase, φ = 0 and φ = 7π/4.Our study reveals that the φ = 0 case is identical to the period-2 model at the DL, although the conditions for DL for the two cases are slightly different.The φ = 7π/4 model is much more rich compared to the earlier models.In the strong driving limit, the first-order Floquet perturbation theory suggests that the φ = π/4 model at the DL points reduces to the SSH Hamiltonian with perfect dimerization, with the hopping amplitudes alternating as γ 1 = 0 and γ 2 = 1 respectively.Hence the system supports robust zero-energy edge modes, which are topologically protected.The entanglement spectrum in this regime demonstrates a finger-like structure with many low-entanglement states near the middle of the quasienergy spectrum.We find that there are some initial states which either show long-time persistent oscillations or do not participate in the dynamics at all.We put forward two possible mechanisms for this ergodicity breaking.(i) There exists an extensive numbers of conserved quantities at the DL points giving rise to sectors which are decoupled from the each other.The number of sectors grows exponentially with the system size as 3 L/2 , which is slower than the growth of the Hilbert space dimension 2 L .We would like to emphasize that these quantities are only approximately conserved quantities; the conservation becomes more and more exact as the driving amplitude is increased.(ii) Another possibility is that there are many configurations of frozen states which do not evolve with time.However, these states quickly thermalize as we move away from a DL point.
In the presence of interactions, the period-4 model with φ = 7π/4 is found to have a Floquet Hamiltonian which describes the transverse field Ising model with longitudinal fields at the boundaries.We find it surprising and remarkable that periodic driving of a period-4 model can be tuned to generate well-known systems like the SSH model and the transverse field Ising model which have been extensively studied for many years.
We also discuss the effects of a staggered on-site potential on the DL.In this case, we find that the finger-like structure in the entanglement spectrum further breaks up into secondary fragments, and the Loschmidt echo produces long-time coherent oscillations at the same fine-tuned parameter values.Next, we examine the stability of this non-ergodic behavior whenever the DL and resonance condition are simultaneously satisfied.To study this regime, we choose two different sets of parameter values, µ = ω = V >> J, and µ = 2ω = V >> J.In both cases, the non-interacting part of the effective Hamiltonian supports DL.For µ = ω = V >> J, the entanglement spectrum again demonstrates many low-entanglement states and slow thermalization of the system.In this regime, the effective Floquet Hamiltonian found using first-order perturbation theory shows that DL and resonance together put strict restrictions on the allowed hopping processes, and these restrictions protect some of the midspectrum states from thermalization.These kinetic constraints on the dynamics generate dynamically disconnected sectors, a phenomenon called Hilbert space fragmentation.For µ = 2ω = V >> J, we see a fractured entanglement spectrum with many segments and with many low-entanglement states near the middle of the quasienergy spectrum.The first-order effective Hamiltonian for this case shows that there are an extensive numbers of conserved quantities.Furthermore, as in the previous case, some processes are again strictly forbidden due to the combination of DL and resonance.These two mechanisms can, in principle, lead the system towards non-ergodic behavior.
We would like to emphasize that the results obtained from the first-order Floquet Hamiltonian (such as the appearance of a large number of conserved quantities) agree well with the results from an exact numerical calculation of the Floquet operator only when (i) the driving amplitude and frequency are much larger than all the other parameters of the system, and (ii) the time scale of observation of correlation functions and Loschmidt echos is not very large.The two sets of results are expected to deviate from each other at very long times because the effects of higher-order terms in the FPT then become important.
In summary, we have presented a number of models in this paper which can be tailored by Floquet engineering to exhibit rich topological and dynamical phase diagrams which have no counterparts in a time-independent (undriven) model.

Appendix A
In this appendix we will use a transfer matrix method to determine the number of zero-energy modes for the Hamiltonian given in Eq. (42).We see from that Hamiltonian that there cannot be any hopping between sites j + 1 and j + 2 if either (i) the occupation numbers at sites (n j , n j+3 ) are either (0, 0) or (1, 1), or (ii) the occupation numbers at sites (n j+1 , n j+2 ) are either (0, 0) or (1, 1).Hence, any configuration which satisfies any of the above conditions for all values of j must necessarily be a zero-energy state.This leads us to define an 8 × 8 transfer matrix T 1 whose rows correspond to the eight possible occupation numbers (n j , n j+1 , n j+2 ), i.e., (111), ( 110), ( 101), ( 100), ( 011), (010), ( 001) and (000), and the columns correspond in a similar way to the eight possible occupation numbers (n j+1 , n j+2 , n j+3 ).The conditions given above imply that T 1 must have the form 1 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 1 We now have to find the eigenvalues of this matrix.We first note that the matrix has two 4 × 4 blocks which are not coupled to each other: the blocks consist of the rows and columns numbered (1235) and (4678), and the two blocks have identical eigenvalues.We can therefore look at either of the two blocks and find the four eigenvalues; the eigenvalues of T 1 will then be given by these four eigenvalues, each repeated twice.The block corresponding to (1235) takes the form We find that one of the eigenvalues of T 1 is zero.The other three eigenvalues must therefore be solutions of a cubic equation which turns out to be The solutions of this equation are found to be where k can take the values 0, 1, 2. We then find that the three eigenvalues are given by 1.466 and −0.233± 0.793i approximately.The eigenvalues of the original transfer matrix T 1 are therefore given by 1.466, −0.233 ± 0.793i and 0, each repeated twice.Hence, for a system with a large number of sites L, the number of zero-energy states grows exponentially as 1.466 L (compared to the total number of states which grows as 2 L ).
It is interesting to compare the number of zero-energy states in this model with the number of such states in a different model which also has a kinetic constraint on the nearest-neighbor hoppings [105].The Hamiltonian of that model is given by It has been shown recently that this model can appear as the effective Hamiltonian of a periodically driven system for some special values of the driving parameters [109].In this model, we see that there cannot be any hopping between sites j + 1 and j + 2 if either (i) the occupation numbers at sites (n j , n j+3 ) are either (0, 1) or (1, 0), or (ii) the occupation numbers at sites (n j+1 , n j+2 ) are either (0, 0) or (1, 1).

Appendix B
In this appendix, we will derive the first-order FPT for a two-site model with nearest-neighbor hopping J, and on-site potentials µ 1 and µ 2 , which can be generalized easily to a system with larger system sizes.We will further assume that the on-site potential is periodically driven in time with a square pulse protocol.Assuming µ 1 , µ 2 >> J, we can recast the Hamiltonian as and The eigenfunctions corresponding to E 1,2 = µ 1 f (t) and µ 2 f (t) are given by |+〉 = 1 0 and respectively.The instantaneous eigenvalues E 1,2 satisfy the degeneracy condition in >> J m = 2 µ = nω (n = 1, 2, • • • ) Many-body flat bands, slow thermalization due to an emergent integrability Model of Hilbert space fragmentation m = 4 (φ = 0) µ = 2nω (n = 1, 2, • • • ) Same as period-2 model Same as the period-2 model m = 4 (φ = 7π/4) µ = nω (n = 1, 2, • • • ) Many lowentanglement states near the middle of the spectrum due to the presence of an extensive number of conserved quantities

Figure 1 :
Figure 1: Dispersion of quasienergies E Fk at DL points: Plots showing E F k as a function of k obtained from the exact numerical calculation for J = 1, µ = 20, and (a) ω = 20 and (b) ω = 10.For both cases the bandwidth ∆ ∼ 0.02, which implies an almost perfectly flat band at zero energy.(c) the same plot for J = 1, µ = 10, and ω = 10.In this case, the bandwidth ∆ = 0.08 is significantly larger than in the first two cases.

Figure 2 :
Figure 2: Quasienergies and crossover behaviors of correlation functions: Plots showing E F k as a function of k obtained from the exact numerical calculation for J = 1, µ = 10, and (a) ω = 10.6 and (b) ω = 10.7.In plot (a), E k ∼ k 4 around k = 0 as can be seen in Eq. (31).In plot (b), E k ∼ k 2 + εk 4 , where |ε| << 1, as seen in Eq. (33).Log-log plots of the absolute value of the n-dependent part of the correlation function δC n as a function of the time nT , for J = 1, µ = 10, and (c) ω = 10.6 and (d) ω = 10.7.(c) The correlation function decays as n −1/4 along with oscillations.(d) The plot shows a crossover between an oscillatory term times n −1/4 and an oscillatory term times n −1/2 .

Figure 3 :Figure 4 :
Figure 3: Quasienergy and entanglement entropy spectrum of the period-2 model at a DL point and away from a DL point: (a-b) Plots of the quasienergy spectrum E F and half-chain entanglement entropy S L/2 as a function of E F exactly at a dynamical localization point with J = 1, µ = ω = 20, and V = 1, where the system exhibits many-body flat bands with many low-entanglement states near the middle of the spectrum.(c-d) Plots showing the same quantities as in plots (a-b) but away from a DL point, with J = 1, µ = 20, ω = 22 and V = 1.For this case, the many-body flat bands start to disappear as we tune the system away from a DL point.

Figure 5 :
Figure 5: Entanglement entropy spectrum of the period-2 model for the case of DL and resonance: Plots showing the entanglement entropy S L/2 as a function of the Floquet quasienergy E F obtained from (a) the exact numerical calculation and (b) the first-order effective FPT Hamiltonian, for J = 1, µ = ω = V = 20, and L = 16.In both cases, we see many low-entanglement states near the middle of the spectrum.In both plots, the color intensity indicates the density the states, implying that the majority of Floquet eigenstates show almost thermal entanglement.Moreover, the quasienergies obtained from the firstorder FPT agree quite well with the exact numerically computed values.However, there are a large number of Floquet eigenstates for which the entanglement estimated from FPT is much smaller than the exact numerically obtained values.

Figure 6 :
Figure 6: Schematic of a Hilbert space fragment for the period-2 case: Figure showing a particular Hilbert space sector consisting of four states, |0000111111110000〉, |0001011111110000〉, |0000111111101000〉, and |0001011111101000〉.The black and white dots indicate occupied and empty sites respectively.

Figure 7 :
Figure 7: Dynamics of the Loschmidt echo and overlaps with Floquet eigenstates in a resonant case: (a-b): Plots of the Loschmidt echo versus time as obtained from (a) the exact numerical calculation and (b) the first-order effective FPT Hamiltonian, for the same parameter values as in Fig. 5, taking the initial state to be |0000111111110000〉.Both plots show show long-time oscillations in time indicating an anomalous thermalization behavior.(c): Overlaps of the same initial state with the Floquet eigenstates as a function of E F computed from the exact numerical calculation for the same parameter values, with a color bar indicating the variation of S L/2 .The quasienergies of the Floquet eigenstates having the highest overlaps with the initial state agree with the analytically predicted values.

Figure 8 :
Figure 8: Schematic of the mirror-symmetric periodic potential pattern for the period-4 model: Schematic diagrams showing the potential patterns for the two mirror-symmetric configurations of the period-4 model corresponding to φ = 0 and 7π/4.

Figure 9 :
Figure 9: Quasienergy and entanglement entropy spectrum of the period-4 model with φ = 0: Plots showing the spectrum of (a) E F and (b) S L/2 as a function of E F at a DL point with µ = 2ω = 20 and V = 0.5.Both the spectra look identical to those of the period-2 case at a DL point.

Figure 10 :
Figure 10: Schematic of the topologically protected edge modes for the period-4 model at a DL point: Schematic picture showing topologically protected zero-energy edge modes for the φ = π/4 at a DL point and no edge modes for φ = 7π/4.
11 (a) and 11 (b), we see two zero-energy edge modes and no edge modes for φ = π/4 and φ = 7π/4, respectively, for J = 1, µ = ω = 20 and L = 2000 with open boundary conditions.As shown in Figs.11(c) and 11 (d), the two modes are localized at the two edges of the system, which can be seen from a plot of the probability |ψ( j)| 2 versus the site number j.Note that the parameter values J = 1 and µ = ω = 20 imply that the system is at a DL point.Another interesting point to observe is that a static model with a nearest-neighbor hopping J and a period-4 time-independent potential with strength µ does not have any such zero-energy end modes.

Figure 11 :
Figure 11: Quasienergies and wave function probabilities of the edge modes: (a-b): Plots of E F in increasing order versus state number for a system with open boundary conditions with J = 1, µ = ω = 20.(a) A system with φ = π/4 supports two zero-energy edge modes, while (b) a system with φ = 7π/4 does not.(c-d): Plots showing the probability versus site number for the edge modes showing that the modes are localized at the ends of the system.

Figure 12 :
Figure 12: Schematic of period-4 model with interaction at a DL point: The schematic picture of the period-4 model at a DL point, with a nearest-neighbor hopping which alternates between J and zero, and an interaction V .

Figure 13 :
Figure 13: Entanglement entropy spectrum for the period-4 case at a DL point and away from a DL point: Plots showing the half-chain entanglement entropy S L/2 versus the quasienergy E F .For all three cases, we take J = 1 and ω = 20.For the first two cases, we consider a DL point with µ = 20, and (a) V = 0.5 and (b) V = 2, respectively.(a): The entanglement entropy consists of many finger-like structures with multiple lowentanglement states near the middle of the spectrum.(b):No finger-like structure is present; however, the system still exhibits multiple low-entanglement states near the middle of the spectrum.(c): We consider a point away from DL, with µ = 10, and V = 0.5.As opposed to the behavior at a DL point, we see the system exhibits low-entanglement states (except at the end points of the quasienergy spectrum where the entanglement is always low).In all the plots, the color intensity indicates the density of states.In plot (c) we see that the majority of the Floquet eigenstates show thermal entanglement.

Figure 14 :
Figure 14: Dynamics of the Loschmidt echo starting from an initial state and the overlap with Floquet eigenstates at a DL point for the period-4 model: (a-b): Variation of Loschmidt echo with time for two different initial states, |210120101〉, and |21102110〉, respectively, for four DL points with µ = 20, V = 0.5, and ω = 5, 6.67, 10 and 20 (satisfying ω = µ/n).In all these cases, the initial states show perfect revivals due to the presence of approximate conserved charges.(c-d): The overlaps of these two initial states with the Floquet eigenstates as a function of E F , with the color bar indicating the variation of S L/2 .The initial states have significant amounts of overlap with a large number of Floquet eigenstates, and some of these eigenstates have extremely low entanglement.

Figure 15 :Figure 16 :
Figure 15: Dynamics of entanglement entropy and correlation function at a DL point and away from a DL point for the period-4 model: (a-b): Plots showing the dynamics of S L/2 and the correlation function 〈n j (nT )n 0 (0)〉 for the initial state |21012101〉 at a DL point with J = 1, µ = 20, ω = 10, and V = 0.5.(c-d): Same plots away from a DL point with J = 1, µ = 20, ω = 11, and V = 0.5.(a): At a DL point, we see that the entanglement entropy increases extremely slowly before reaching a saturation value which is less than the thermal value.(b): The correlation function shows a behavior similar to S L/2 with a longtime revival pattern.(c): Away from a DL point, S L/2 reaches a saturation value soon after an initial growth in time.(d): The correlation function demonstrates a similar behavior, suggesting thermalizing behavior away from a DL point.Note that the time scales in (c-d) are much shorter than in (a-b)

Figure 17 :
Figure 17: Entanglement entropy spectrum with a staggered on-site potential at a DL point: Plots of S L/2 versus E F at a DL point with J = 1, µ = ω = 20, V = 0.5, and (a) w = 1 and (b) w = 3.In both cases we see many finger-like structures with multiple secondary fragments, with the effect becoming clearer with increasing strength of the staggered potential w.The color intensity indicates the density of states, revealing that the majority of states demonstrate athermal behavior with S L/2 being much smaller than the thermal value.

Figure 18 :
Figure 18: Dynamics of Loschmidt echo in the presence of a staggered potential: Plots showing the dynamics of Loschmidt echo for the initial state |21012101〉 at four DL points with J = 1, V = 0.5, w = 3, µ = 20, and ω = 5, 6.67, 10 and 20.In all four cases, the Loschmidt echos show an oscillatory behavior with an extremely slow decay, suggesting that the system retains the information of the initial state for a long period of time.

Figure 19 :Figure 20 :
Figure 19: Entanglement entropy spectrum for a resonant case at a DL point of the period-4 model: Plots showing the entanglement entropy S L/2 versus the quasienergy E F obtained from (a) an exact Floquet calculation and (b) the first-order FPT Hamiltonian, for a DL point exhibiting a resonance with J = 1 and µ = ω = V = 20.The color intensity indicates the density of states, suggesting that most of the states attain the thermal value.However, there are also many low-entanglement states present near the middle of the spectrum.The quasienergy spectrum obtained from FPT (b) agrees well with the exact numerically computed spectrum (a).However, the entanglement entropy obtained from the first-order FPT is much less than the exact numerically obtained values for many of the states.

Figure 21 :
Figure 21: Entanglement spectrum at another DL point and at resonance for the period-4 model: Plots showing the entanglement entropy S L/2 obtained from (a) exact numerical calculations and (b) first-order Floquet Hamiltonian, at a DL point with J = 1, and µ = 2ω = V = 20.Both calculations show that the quasienergy spectrum consists of multiple fragments with several low-entanglement states lying near the middle of the spectrum.The color intensity indicates the density of states, revealing that the majority of states do not attain the thermal value of te entanglement entropy which is given by the upper envelope of the plots.