Homogeneous Floquet time crystal protected by gauge invariance

Angelo Russomanno, 2 Simone Notarnicola, Federica Maria Surace, 4 Rosario Fazio, 2 Marcello Dalmonte, 4 and Markus Heyl Abdus Salam ICTP, Strada Costiera 11, I-34151 Trieste, Italy NEST, Scuola Normale Superiore & Istituto Nanoscienze-CNR, I-56126 Pisa, Italy Dipartimento di Fisica e Astronomia “Galileo Galilei” & INFN, via Marzolo 8, I-35131, Padova, Italy SISSA, Via Bonomea 265, I-34136 Trieste, Italy Max-Planck-Institut für Physik Komplexer Systeme, Nöthnitzer Strasse 38, D-01187, Dresden, Germany

Introduction.-Isolated quantum matter can feature phases with long-range order in highly excited states that cannot be captured by thermodynamic ensembles [1,2]. This crucially relies on ergodicity breaking and a failure of the Eigenstate Thermalization Hypothesis (ETH) [3]. One robust mechanism for achieving such nonergodic behavior is to impose strong disorder giving rise to the many-body localized (MBL) phase [2,[4][5][6][7][8], which can host long-range ordered phases such as the MBL-spin glass [1,9] or discrete time crystals [10][11][12][13][14]. Recently, it has been realized that lattice gauge theories (LGTs) entail another robust mechanism for nonergodic dynamics in short-ranged systems protected by gauge invariance instead of disorder [6][7][8] due to the specific structure of their Hilbert spaces, which are built up of disconnected superselection sectors [7]. However, it has remained an open question to which extent they can also accommodate nonequilibrium phases with long-range order and therefore to which extent they can contribute to the open quest of realizing robust nonequilibrium ordered phases of homogeneous quantum matter.
In this work we focus on the time-crystal physics [10, and show that homogeneous LGTs can feature robust time-crystalline phases protected by gauge invariance instead of disorder. We identify a period-doubling discrete Floquet time crystal in a periodically kicked Z 2 LGT, see Fig. 1. This 'gauge time crystal' is characterized by both spatial and temporal long-range order. We solve the kicked Z 2 LGT exactly by a mapping onto a free fermionic theory using a Jordan-Wigner (J-W) transformation, which allows us to explore the phase diagram for large system sizes. We observe that the Floquet eigenstates appear in pairs with a quasienergy difference of π, so that our system shares many of the features of the π-spin glass in a periodically kicked Ising chain with quenched disorder [10]. Importantly, we find that this gauge time crystal represents a robust phase which does not require fine tuning and persists over a wide range of parameters. In particular, we also study the influ-ence of perturbations breaking the exact solvability and preserving the Z 2 gauge symmetry, where we find numerical evidence for stability by means of exact diagonalization. The mechanism behind this time-crystalline phase relies solely on gauge invariance and can therefore be directly extended to other LGTs with discrete gauge symmetries. Importantly, our observation of a robust timecrystalline phase in a homogeneous short-ranged system goes beyond recent approaches which lead to prethermal spatiotemporal order [35][36][37][38][39], consider fermionic systems with strong short-range interactions giving rise to emergent Floquet integrability [40], or require long-range interactions in driven unitary [24,29] and dissipative dynamics [41][42][43][44][45].
The model.-We consider a Z 2 Higgs-LGT in one spatial dimension. The theory describes the dynamics of Higgs fields -defined by Pauli-matrix operatorsξ α j at vertex j on the lattice -coupled to Z 2 gauge fields -defined by Z 2 parallel transportersτ x j,j+1 at the bond (j, j + 1) as illustrated in the box of Fig. 1. The system Hamiltonian reads [46,47] The Higgs-field operators can also be interpreted as hardcore bosonsb j withξ x j =b † j +b j . The first two terms denote mass and gauge interactions, while the third describes the coupling between the Higgs and gauge fields. We drive the Z 2 Higgs-LGT out of equilibrium by periodically kicking the strength of the Higgs-gauge coupling, leading to the following time-dependent Hamiltonian x jτ z j, j+1ξ x j+1 (2) This system exhibits a local symmetry:Ĥ(t) commutes with the operatorsĜ j = −τ  be understood as the complex exponentials of the local Gauss' operators). Thus, the Hilbert space of size 2 2L−1 is partitioned in N = 2 L superselection sectors, where all the states |Ψ { qα} in a given sector are identified by the same set of local static charges q j = ±1 viâ In the following we consider initial product states of the form |Ψ = |ϕ H ⊗ |ψ g where |ϕ H is a product state which satisfies H ϕ|ξ z j |ϕ H = 0 for all j = 1, . . . , N and |ψ g is the initial condition for the gauge degrees of freedom, that we will specify later in the text. Such initial conditions, which represent superpositions over many superselection sectors, can yield robust nonergodic behavior for LGTs and disorder-free localization [6][7][8]48]. Concretely, for our Z 2 LGT the dynamics in a given superselection sector specified by the charges {q α } is determined by an effective Hamiltonian (see SM Sec. I) . This integration is related to the duality between Ising models and Ising LGTs [49,50]. As a result the Hamiltonian becomes a kicked transverse-field Ising chain with binary bond disorder due to q j = ±1, which can be solved exactly via a J-W transformation for large systems. We will also study the influence of a perturbation of the formĤ K = 4K breaking the J-W solvability. After the integration it adds a transverse interaction term for the gauge fieldŝ We solve the dynamics of the LGT in a set of N real randomly chosen superselection sectors and finally perform an average when computing observables. In the shown data we include error bars resulting from the finiteness of N real . But let us emphasize again that the overall problem is homogeneous both in the initial condition and in the Hamiltonian. Initial conditions and observables.-In order to reveal both the temporal and spatial order we use two complementary setups.
On the one hand, we take initial conditions which explicitly break the Z 2 symmetry of the model yielding a nonzero magnetization m x for the gauge degrees of freedom which we then monitor in the subsequent evolution: where we have defined · · · t ≡ g ψ(t)| · · · |ψ(t) g and the overline marks the average over the N real pseudodisorder realizations [51]. In this way we obtain direct access to the time-crystalline period-doubling dynamics. In Fig. 1(b) we show results for m x (t) in the fully interacting case K = 0 obtained through exact diagonalization. We see the existence of period-doubling oscillations which are persistent for an infinite time in the thermodynamic limit. We show this fundamental property of persistence [11,28] in Fig. 1(c), where we see that the decay time t * of the period doubling oscillations exponentially scales to infinity with the system size. We determine t * as the time after which (−1) t/T m x (t) changes sign [30,52] averaged over disorder.
On the other hand we can choose initial conditions which are Z 2 -symmetric with a vanishing magnetization m x (t), which allows us to address the spatial long-range ordering in the system. For that purpose we study the correlation parameter with · · · t defined as above. Whenever S xx t > 0 while at the same time m x (t) = 0, the system exhibits long-range spatial order.
Exactly solvable case.-Let us first focus on the case with K = 0, where the model can be mapped onto a system of non-interacting fermions by means of a J-W transformation. In each superselection sector {q α } we initialize the dynamics with the same initial state |ψ g chosen as the ground state of the Hamiltonian . This state has a non-vanishing correlation parameter if h 0 < 1 and is symmetric under Z 2 which allows us to address the long-range spatial ordering in the system; for a study of the temporal order we perform a spectral analysis, as we are going to detail below. In the J-W framework it is well known how to numerically study the dynamics and how to evaluate the correlation parameter as a Pfaffian (see [53][54][55][56][57]). Here it is enough to say that the dynamics is induced by an effective 2(L − 1) × 2(L − 1) timeperiodic single-particle Hamiltonian. This is important to mention because we can compute the 2(L − 1) singleparticle Floquet states and the 2(L − 1) single-particle quasi-energies α (see for instance [58]). These quantities will play an important role in what follows.
We find that the correlation order parameter reaches an asymptotic value S xx asy after a transient (see the discussion below Eq. (7)). We plot the long-time value of S xx t as a function of kicking strength φ for different values of L in the main panel of Fig. 2. We observe three regimes whose separating phase boundaries we indicate by the colored zones. In the regimes i) and iii) S xx asy con- verges to a nonzero value as L → ∞, while in regime ii) S xx asy vanishes as the L is increased (see also the inset of Fig. 2). Both regions i) and iii) mark the existence of an eigenstate phase [1,2], where eigenstates exhibit long range spatial order (as in [9], for instance). This eigenstate order is protected by nonergodicity since in a thermalizing system such order is impossible due to the Mermin-Wagner theorem.
Although the behaviour of S xx asy is qualitatively similar in both i) and iii), these two regions mark different phases since i) in addition also supports temporal order. An example of this property for φ = 1.02π can be seen in Fig. 1(c) (curve with K = 0): The system is initialized in a state explicitly breaking the Z 2 symmetry and the decay time t * exponentially increases with the system size. This fact can be understood by an analysis of the Floquet spectrum [10]. The presence of a temporal time-crystal ordering corresponds to spectral pairing, where each Floquet state has a partner with quasienergy shifted by π. This situation is realized if there is a singleparticle quasienergy exactly at π with a marked gap separating it from the rest of the spectrum. In this way it does not hybridize with the bulk, and each many-body Floquet state has a π-shifted partner obtained by adding the quasiparticle with quasienergy π. We evaluate this gap as δ π = 1 [59] and plot it in Fig. 3. We see that it is non-vanishing in all the regime i). Moreover, as we show in the SM (Sec. II), in this regime 2L−2 averaged over the disorder is exactly equal to π. In SM-Sec. II we show also that the single-particle bulk Floquet states are always Anderson localized. This is very important, because without localization it is possible to have a gap in the Floquet spectrum at π and still observe no time crystal (see for instance [58]): In the absence of localization, local operators expand in time obeying the Lieb-Robinson bound and no time-periodic behaviour whatsoever is possible [25]. Of course, the transition to localization and the one to glassy order of the excited eigenstates are independent [9], and this is the reason why the transition from regime i) and ii) occurs at a value of φ different from the one where δ π vanishes.
In Fig. 2 we have initialized with a specific value of h 0 , but we have checked that the presented phenomenology doesn't depend on this choice.
General case.-At this point we break J-W solvability by considering the term of Eq. (4), with K = 0. We consider a value of φ for which we see this phenomenon at K = 0; then we take K = 0 and we study the properties of the asymptotic correlation parameter. An interval of K where this quantity does not scale with the size would mark the persistence of the time crystal. We now perform a conventional exact-diagonalization simulation of the system, up to size L = 13. To evaluate the asymptotic correlation parameter, we can resort to the diagonal ensemble and we get where |φ β g are the many-body Floquet states, N is the dimension of the Hilbert space and R β ≡ g ψ(0) |φ β g denotes the overlap with the initial state. We remark that we can use Eq. (7) even if the many-body Floquet quasienergies µ β appear in degenerate pairs, due to the Z 2 symmetry. The point is that the operatorsτ x j, j+1τ x i, i+1 commute with the same Z 2 symmetry and hence have no matrix elements between states with different parity (see SM-Sec. III for details).
We plot the dependence of S xx asy versus K for different L in Fig. 4. We take two different initial conditions, in the upper panel we take the state with all the spins pointing down along the x axis (|ψ g = |s x 1,2 = −1 . . . s x j,j+1 = −1 . . . s x L−1,L = −1 g ), in the lower panel we take the uniform superposition of all the eigenstates ofτ x j, j+1 ∀ j obeying the condition L−1 j=1 s x j,j+1 ≤ −(L − 1)f with f = 0.8. We see that for K 0.2 there is no decrease with L, marking the persistence of the time-crystal behaviour. This persistence can be seen also in Fig. 1(c) where the t * introduced above exponentially increases with L.
Time crystallinity in Abelian lattice gauge theories.
-We now investigate more generally if time crystallinity can appear in disorder-free Abelian LGTs in (1+1)-d. We consider the generic Hamiltonian coupling Higgs fields to Abelian gauge fields [60]: wheren j =φ † jφ j is the Higgs occupation on site j and E j,j+1 ,Û j,j+1 are respectively the electric field and the parallel transporter, andĤ(t) is defined analogously to the Z 2 case above. The electric-field interaction energy is local in these theories, differently from the Z 2 term involving at least two neighboring sites. For a Z N LGT (i.e. a theory where nowÛ j,j+1 andÊ j,j+1 are not Pauli matrices but the more general clock operators), we can use a similar approach as the one used in the Z 2 LGT. We consider an initial state where matter is in an equalweight superposition of all possible eigenvalues of the Higgs number operator, and the gauge fields are in a generic state. The evolution of such states can be mapped exactly into the one of Z N clock models under the effect of quasi-random local fields: since the latter class of models has been shown to display time-crystal behavior for small values of N and random disorder [30], it is natural to expect that the mechanism discussed above holds true also for N > 2. This mechanism does not work for continuous U (1) LGTs (see SM-Sec. IV for details), which, however, doesn't exclude other ones for the generation of time crystals in such theories.
Concluding discussion.-In this work we have demonstrated that homogeneous LGTs can realize timecrystalline phases, where the protecting nonergodicity is enforced by the local constraints imposed by gauge invariance. In more general terms, our results show that homogeneous LGTs can realize eigenstate order, which naturally leads to the question to which extent also other eigenstate phases can occur in homogeneous LGTs, e.g., analogues of the MBL-spin glass [1,9] or topological order at elevated energy densities [61].
Further, our results can be directly extended to Z N LGTs which opens up the possibility, in principle, of generating period N -tupling time-crystals. While our approach cannot be immediately applied to LGTs with continuous groups, it would be intriguing to see whether discrete non-Abelian symmetries can also support the formation of defect-free time crystals.

ACKNOWLEDGMENTS
We acknowledge fruitful discussions with R. Khasseh and M. Wauters. A. R. and R. F. thank the European Union for partial financial support through QUIC project (under Grant Agreement 641122). A. R. thanks the Max-Planck-Institut für Physik Komplexer Systeme for partial financial support and the warm hospitality received during the preparation of this work. S. N. and M. D. acknowledge partial support by the H2020 Project -QUANTUM FLAGSHIP -PASQUANS (2019-2022). We are indebted with G. E. Santoro for the subroutine performing the diagonalization of unitary operators. This work is partly supported by the ERC under grant number 758329 (AGEnTh) and by the QUANTERA project QTFLAG.

I. DERIVATION OF THE EFFECTIVE HAMILTONIANĤ {qα} (t)
The derivation of the effective HamiltonianĤ {qα} (t) fromĤ(t), defined respectively in Eqs. (3) and (2) of the main paper, needs two steps. In the first, we restrict ourselves to one of the superselection sectors defined by a set of static charges {q α }. To do so, for a generic state |Ψ , we consider its component |Ψ {qα} on the chosen sector, defined as where P {qα} = j P j (q j ) is the projector on the chosen superselection sector and P j (q j ) = (1 − q jτ x j−1, jτ x j, j+1ξ z j )/2 projects on the sector with static charge q j on site j. It follows that for each state |Ψ {qα} in the chosen sector we have In the second step, we exploit the above relation in order to cancel the matter field operatorsξ α j from the Hamilto-nianĤ(t). The derivation is now straightforward. Then, we redefine the operatorsτ α j, j+1 in order to cancel the matter field from the second part ofĤ(t) Note that the proper commutation relations are still satisfied, in particular τ z j, j+1 2 = 1 and [τ z j, j+1 , τ z k, k+1 ] = 0 ∀ k, j. By applying this substitution to the HamiltonianĤ(t) we have that where the prime will be henceforth omitted. The same substitution allows to cancel out the matter field in the termĤ K , which is introduced in Eq. (4) of the main paper as a function of the gauge field only. We thus obtain the effective Hamiltonian (Eq. The state for the Higgs field is a product state of spins, each one living on the equator of the Bloch sphere ( ξ z j = 0 for every j). For each j we can find the unit vector n j = cos θ jx + sin θ jŷ giving its position on the Bloch sphere, such that We now consider a generic sector {q α } and we see that for every site i we have where we used the fact that {( ξ i ·n i ),ξ z i } = 0. From the unitarity of ( ξ i ·n i ) we derive the relation between the norms which implies that the projections on sectors which only differ for one local charge have equal norms. Since the last relation holds for every set of {q α } and for every i, we find that all the probabilities p {qα} have to be equal, with p {qα} = 1/N .

II. SINGLE-PARTICLE FLOQUET SPECTRUM
Fig. S1 shows the average over N real pseudo-disorder realizations of 2(L−2) . We can see that it stays at π in an interval of φ larger than the one where δ π is nonvanishing (see Fig. 3 of the main paper). Fig. S2 shows the bulk-averaged single-particle Floquet inverse participation ratio defined as where we define · · · as the average over the N real realizations of pseudo-disorder, w α = u 1, α · · · u L, α | v 1, α · · · v L, α are the singleparticle Floquet states (see for instance [58] for more details) and α runs over the 2(L − 2) values corresponding to the bulk single-particle Floquet states (the ones with α = ±π). For all the considered values of φ we can clearly see that IPR bulk does not scale with the system size and the same occurs for the error bars (evaluated as the r. m. s. fluctuation over the pseudo-disorder realizations). This marks the fact that all the single-particle Floquet states are localized and therefore the model shows Anderson localization for all the considered values of φ.

III. CONVERGENCE TO THE FLOQUET DIAGONAL ENSEMBLE
In the text we have claimed that the correlation order parameter converges for long times towards the Floquet  diagonal ensemble value given by Eq. 6 of the main text. We have numerically verified this point; we show an example of this convergence in Fig. S3. We would like also to better discuss this point from the theoretical point of view. Let us expand S xx t in the Floquet basis, we find (S11) The indexes σ and σ can take values + and − and mark the Z 2 symmetry sector. The system is symmetric under the Z 2 symmetry, so the Floquet states are doubly degenerate µ + β = µ − β ∀β. Moreover, also the operators τ x j, j+1 τ x i, i+1 are symmetric under this symmetry and we have therefore g φ + β | τ x j, j+1 τ x i, i+1 |φ + β g = g φ − β | τ x j, j+1 τ x i, i+1 |φ − β g . We can therefore rewrite Eq. (S11) as (S12) The off-diagonal term vanishes in the long time after the disorder average, due to the destructive interference between the oscillating phase factors. Only the blockdiagonal term survives. Thanks to the degeneracy of the expectations with respect to the index σ, we can compute this term directly using the Floquet states given by the numerical diagonalization (which, for each β, are in general superpositions of σ = + and σ = −). Moreover, |R + β | 2 +|R − β | 2 takes the same value whichever basis in the degenerate Floquet subspace is considered. That's why in the main text we do not write the index σ.

IV. NO TIME CRYSTAL FOR CONTINUOUS GAUGE SYMMETRY
Here we briefly discuss the case of a continuous gauge symmetry. In order to show that the time-crystal question in this case is very delicate (and most probably a dis-crete time-crystal behaviour in this form is impossible) let us consider the lattice Schwinger model in the Wilson formulation as an example. As discussed in Ref. [7], this model displays an extremely slow dynamics, which is qualitatively less ergodic than conventional MBL -in particular, with entanglement entropy growing as ln(ln(t)). However, due to the absence of any periodicity in the gauge field Hilbert space, it is not possible to identify a clear time-dependent Hamiltonian whose dynamics could lead to a time crystal: for instance, applying the same recipe as above would lead to an infinite period. This case immediately illustrates that the absence of ergodic dynamics -typical of gauge theories due to superselection sectors -is by far not enough for engineering translationinvariant time crystals: Identifying the proper gauge symmetry is absolutely key and, in the present context, doable thanks to rather direct analogies with inhomogeneous clock models emerging from a specific class of initial states.