Digital quantum simulation of Floquet symmetry-protected topological phases

Quantum many-body systems away from equilibrium host a rich variety of exotic phenomena that are forbidden by equilibrium thermodynamics. A prominent example is that of discrete time crystals1–8, in which time-translational symmetry is spontaneously broken in periodically driven systems. Pioneering experiments have observed signatures of time crystalline phases with trapped ions9,10, solid-state spin systems11–15, ultracold atoms16,17 and superconducting qubits18–20. Here we report the observation of a distinct type of non-equilibrium state of matter, Floquet symmetry-protected topological phases, which are implemented through digital quantum simulation with an array of programmable superconducting qubits. We observe robust long-lived temporal correlations and subharmonic temporal response for the edge spins over up to 40 driving cycles using a circuit of depth exceeding 240 and acting on 26 qubits. We demonstrate that the subharmonic response is independent of the initial state, and experimentally map out a phase boundary between the Floquet symmetry-protected topological and thermal phases. Our results establish a versatile digital simulation approach to exploring exotic non-equilibrium phases of matter with current noisy intermediate-scale quantum processors21.

In general, SPT states are characterized by non-trivial edge states that are confined near the boundaries and protected by certain symmetries [21,22].In a clean system without disorder, these edge states in general only occurs for the ground states of the system with a bulk energy gap.On raising temperature, they are typically destroyed by mobile thermal excitations.However, adding strong disorder can make the system many-body localized (MBL) [23,24], allowing for a sharply defined topological phase and stable edge states even at infinite temperature [25,26].Strikingly, the topological phase and corresponding edge states can even survive external drives, as long as the driving frequency is large enough so that the localization persists [27].
The interplay between symmetry, topology, localization, and periodic driving gives rise to various peculiar phases of matter that exist only out of equilibrium [27].Understanding and categorizing these unconventional phases poses a notorious scientific challenge.On the theoretical side, topological classifications of periodically driven (Floquet) systems with [28][29][30] and without [31] interactions have already been ob-tained through a range of mathematical techniques (such as group cohomology), revealing a number of mysterious phases with no equilibrium counterparts [27].Yet, we still lack powerful analytical tools or numerical algorithms to thoroughly address these phases and their transitions to other ones.On the experimental side, observations of signatures for time crystals [1][2][3][4][5][6][7][8], which is a paradigmatic example of exotic phases beyond equilibrium [32], have been reported in a wide range of systems [9][10][11][12][13][14][15][16][17][18][19][20].However, none of these experiments encompasses topology as a key ingredient and the realization of a topological time crystal [33,34], which demands delicate concurrence of topology and localization as well as periodic driving, still remains a notable experimental challenge so far.
In this paper, we report the observation of a SPT time crystal with an array of programmable superconducting qubits.With high controllability and long coherence time, we successfully implement the dynamics of the system under a prototypical Floquet Hamiltonian, which we theoretically predict to exhibit a SPT time crystalline phase.We measure the temporal correlations and the local spin magnetizations, and demonstrate that both of them show a subharmonic response at the boundaries but not in the bulk of the chain.The subharmonic response, which manifests the spontaneous breaking of discrete time-translational sysmmetry carried by the Floquet Hamiltonian, maintains for an extended parameter region and is robust to various experimental imperfections, independent of the initial states.In addition, the SPT time crystal is further explored from the perspectives of entanglement dynamics, entanglement spectrum, and the dynamics of stabilizers that underlies its topological nature.Through measuring the variance of the subharmonic peak height in the Fourier spectrum, we also experimentally map out the phase boundary be- Symmetry protected topological time crystal and schematics of the experimental setup.a, Fourteen qubits used in our experiment are coupled to their neighbors with capacitive couplers.b, A chain of spins are periodically driven with the stroboscopic Floquet Hamiltonian H(t), giving rise to a SPT time crystal characterized by spontaneous time translational symmetry breaking at the boundaries.c, The quasienergy spectrum ( ) of the Floquet unitary, which is the time evolution operator over one period.For the SPT time crystal, every eigenstate is twofold degenerate and has a cousin separated by quasienergy π. d, A schematic illustration of the experimental circuits used to implement the time dynamics governed by the Floquet Hamiltonian H(t).We randomly sample the Hamiltonians and prepare the initial states onto random product states or random SPT states.After running a sequence of quantum gates, we measure the local magnetization for each qubit or stabilizers at discrete time points.e, Illustration of the quantum processor, with the 14 qubits highlighted in green.
tween the time crystalline and thermal phases.The observed SPT time crystal in our experiment differs drastically from all other conventional ones with trivial topology, which opens a door for harnessing this exotic phase of matter for practical quantum information processing [35].

Model Hamiltonian
We consider an one-dimensional (1D) spin- 1  2 chain governed by the following Floquet Hamiltonian (see Fig. 1b): where σx,z k is the Pauli matrix acting on the k-th spin; J k , V k , and h k are random parameters drawn independently from uniform distributions over and [h−∆ h , h+∆ h ], respectively.For simplicity, we fix T = 2T 1 = 2 throughout this paper.We note that H(t) has a Z 2 ×Z 2 symmetry.For a suitable parameter region, it has been shown that H 2 can be in a MBL phase, where topological edge states can survive as coherent degree of freedom at arbitrarily high energies [26].The localization and edge states carry over to the case of periodic driving with the Hamiltonian H(t), giving rise to a SPT time crystalline phase (see Supplementary Materials I. B).
The Floquet unitary that fully characterizes the SPT time crystal reads U F = U 2 U 1 , where U 1 = e −iH1 and U 2 = e −iH2 are the unitary operators generated by the Hamiltonian H 1 and H 2 , respectively.To understand why time translational symmetry only breaks at the boundary but not in the bulk, we consider the idealized cluster limit (V k = h k = 0) and set δ = 0. We suppose that the system is initially prepared to a random product state in the computational basis, and we use the dynamics of local magnetization as a diagnosis.In this simple scenario, the non-trivial structure of the cluster states (eigenstates of U 2 ) gives rise to edge modes that behave as free spins.At each driving period, the unitary operator U 1 basically flips all spins.As a result, the edge spins are reversed after one period and return to their initial configuration after two, leading to the doubling response periodicity for the local magnetization at the boundaries.For spins in the bulk, however, the unitary operator U 2 plays a role and evolves the random product state to a state with vanishing magnetization, resulting in no double periodicity or symmetry breaking.A closer look at the quasienergy specturm of U F reveals that its every eigenstate is two-fold degenerate with a long-range spatial order and has a cousin eigenstate separated by quasienergy π (see Fig. 1c).This is essential for the robustness of the subharmonic response of the edge spins against local perturbations.

Experimental setup
Our experiment is performed on a flip-chip superconducting quantum processor designed to encapsulate a square array of 6 × 6 transmon qubits with adjustable nearest-neighbor couplings (Fig. 1e), on which a chain of up to L = 14 qubits, denoted as Q 1 through Q L , that alternate with L − 1 couplers, denoted as C 1 through C L−1 , are selected to observe the SPT time crystal (Fig. 1a).All L qubits can be individually tuned in frequency with flux biases, excited by microwaves, and measured using on-chip readout resonators; all couplers are also of transmon type with the characteristic transition frequencies higher than those of the qubits, which can be controlled with flux biases to tune the effective nearest-neighbor couplings.During an experimental sequence (Fig. 1d), we first initialize each qubit, Q j , in |0 at its idle frequency ω j , following which we alternate the single-qubit gates at ω j with the two-qubit controlled-π (CZ) gates realized by biasing Q j and its neighboring qubit to the pairwise frequencies listed in ω A (B) j , ω A (B) j+1 for a fixed interaction time (see Supplementary Materials III.C).Meanwhile, each coupler is dynamically switched between two frequencies [36][37][38][39][40][41]: one is to turn off the effective coupling where the neighboring two qubits can be initialized and operated with single-qubit gates; the other one is to turn on the nearest-neighbor coupling to around 11 MHz for a CZ gate.After n number of layers of the alternating single-and two-qubit gates, we finally tune all qubits to their respective ω m j for simultaneous quantum-state measurement.Qubit energy relaxation times measured around ω j are in the range of 11 to 37 µs.More characteristic qubit parameters, including the above mentioned frequencies, anharmonicities, and readout fidelities, can be found in Supplementary Materials Tab.S1.
We explore a quantum digital simulation scheme to implement the dynamics of the system under the driven Hamiltonian H(t).More specifically, we decompose the evolution op-erators to the experimentally feasible single-qubit gates [X(θ), Y(θ), and Z(θ)], and the two-qubit gates [CR z (±π)], where X(θ), Y(θ), and Z(θ) are rotations around x-, y-, and z-axis by the angle θ, respectively, and CR z (±π) are the z-axis rotations of the target qubit by ±π conditioned on the state of the control qubit (see Fig. 1d and Supplementary Materials III.A for the ansatz that generating the gate sequences).X(θ) and Y(θ) are realized by applying 50 ns-long microwave pulses with full width half maximum of 25 ns, whose quadrature correction terms are optimized to minimize state leakages to higher levels [42].Simultaneous randomized benchmarkings indicate that the single-qubit gates used in this experiment have reasonably high fidelities, averaging above 0.99 (see Supplementary Materials Tab.S1).Z(θ) is realized using the virtual-Z gate, which encodes the information θ in the rotation axes of all subsequent gates [43], and is combined with CZ to assemble CR z (±π).Here, we adopt the strategy reported elsewhere [40,44] to realize the CZ gate, i.e., diabaticly tune the coupler frequency while keeping |11 and |02 (or |20 ) for the subspace of the two neighboring qubits in near resonance.The 40 ns-long CZ gate for a pair of neighboring qubits can be individually optimized to be around 0.99 in fidelity as calibrated by interleaved randomized benchmarking; when simultaneously running the CZ gates for multiple pairs of neighboring qubits as required in the experimental sequence, the averaged CZ gate fidelities can be around 0.985 as obtained by simultaneous randomized benchmarkings (see Supplementary Materials Tab.S1).

Symmetry breaking at boundaries
The characteristic signature of a SPT time crystal is the breaking of the discrete time-translational symmetry at the boundaries of the chain but not in the bulk.This can be manifested by the persistent oscillation with period 2T of local magnetizations at the boundaries.In Fig. 2, we plot the time evolution of disorder-averaged local magnetizations σ z j (t) for different phase regions.From Fig. 2a, it is evident that in the SPT time crystal region, the disorder-averaged magnetizations at the two ends of the chain, namely σ z 1 (t) and σ z L (t), oscillate with a 2T periodicity, up to over 20 driving cycles.In stark contrast, the magnetizations at the bulk of the chain (σ z j (t) with 2 ≤ j ≤ L − 1) decay quickly to zero and do not show periodic-doubled oscillations.This unconventional behavior is independent of disorder average.Even for a single random disorder instance, the magnetizations show similar dynamical features, as shown in Fig. 2b.The distinction between the dynamics of boundary and bulk magnetizations can also be clearly seen by examining σ z j (t) in the frequency domain.As shown in Fig. 2d, the edge spins lock to the subharmonic frequency of the drive period ω/ω 0 = 1/2 whereas the bulk spins show no such peak.We stress that the subharmonic response for the edge spins obtained in our experiment is notably robust to various perturbations and experimental imperfections (see Supplementary Materials I. B for a more in-depth discussion).For comparison, we also experimentally measure the dynam- While the bulk magnetization decays quickly to zero, the edge spins oscillate with a stable subharmonic response up to 20 cycles.b, The evolution dynamics of local magnetizations for different random instances.Here, each layer corresponds to a specific random instance.c, Magnetization dynamics deep in the thermal region (J = ∆J = 1, V = h = ∆V = ∆ h = 0, and δ = 0.8).d, Fourier transform of experimentally measured σ z j (t) in the time crystal region.The edge spins lock to the subharmonic frequency, which is in sharp contrast to that for bulk spins.e, Fourier spectra of σ z j (t) in the thermal region.No robust subharmonic frequency peak appears for both edge and bulk spins in this case.f, Time-dependence of the autocorrelator Aj = Ψ0|σ z j (t)σ z j (0)|Ψ0 up to 40 cycles, obtained from averaging over 20 random instances deep in the SPT time crystalline phase, with the initial states prepared at random product states |{0, 1} ⊗L .The black solid lines show the results of "echo" circuits for two boundary qubits.ics of the magnetizations in the thermal region.Our results are shown in Fig. 2c and Fig. 2e, where we see that the magnetizations for both the edge and bulk spins decay quickly to zero and no subharmonic response appears at all.
The spontaneous breaking of the discrete time translational symmetry at the boundaries can also be detected by the autocorrelators defined as A j = Ψ 0 |σ z j (t)σ z j (0)|Ψ 0 .Our experimental result for autocorrelators up to 40 driving cycles is plotted in Fig. 2f, showing the breaking of time translational symmetry at the boundaries but not in the bulk again.We mention that in the SPT time crystal region the local magnetizations for the edge spins exhibit a gradually decaying envelope, which is attributed to either the external circuit errors or slow internal thermalization.To distinguish these two mechanisms, we carry out an additional experiment about the "echo" circuit U echo ≡ (U † F ) t U t F , whose deviation from the identity operation measures the effect of circuit errors [18].The square root at the output of U echo (black solid lines shown in Fig. 2f) fits well with the decaying envelop of U F , indicating that the decay of the envelope is due to circuit errors rather than thermalization, which assures that the system is indeed in the MBL region.

Localization-protected topological states
In the above discussion, the initial states are random product states.To establish SPT time crystal, additional experiments on other initial states and other local observables are advisable.In this section, we show that the stabilizers in the bulk do not break the discrete time translational symmetry, but at the boundaries do.To understand this, we consider the idealized cluster and spin-flip limit, i.e., V k = h k = 0 and δ = 0.In this limit, H 2 reduces to a summation of stabilizers: To break the degeneracy of H 2 , we consider adding two boundary terms 2 ) and J L S L (S L ≡ σx L σz L−1 ), which are commuting with all bulk stabilizers, into the Hamiltonian H 2 and form a new cluster Hamiltonian We note that the eigenstates of H 2 are also eigenstates of H 2 with degeneracy split by the boundary terms.We choose the initial states to be random eigenstates of H 2 and evolve the system with the Floquet Hamiltonian H(t).We measure the time-dependence of local stabilizers.In Fig. 3a, we show a sketch of the quantum circuit used in our experiment to prepare the desired random eigenstates of H 2 .To manifest the topological nature of these eigenstates, we study their entanglement spectra [45], which are widely used as a crucial diagnostic for universal topological properties of quantum phases [45][46][47][48].In our experiment, we prepare random eigenstates of H 2 with both open and periodic boundary conditions and obtain the reduced density matrices of half of the system ρ half through quantum state tomography.Figure 3b displays the entanglement spectra (− ln v where v stores the eigenvalues of ρ half ) of two experimentally prepared eigenstates for the open and periodic boundary condition, re-spectively.From this figure, a clear two-fold degeneracy for the low entanglement energy Schmidt states is obtained for the open boundary condition, which corresponds to the spin half degrees of freedom released at the bipartitioning boundary.Whereas, for the periodic boundary condition, the spectrum is four-fold degenerate, corresponding to two effectively decoupled spins at the two boundaries of the bipartition.The degeneracy of the entanglement spectra and its dependence on boundary conditions marks a characteristic feature of the SPT state prepared in our experiment.We remark that the degeneracy disappears at high entanglement energy.This is due to the finite size effects and experimental imperfections.
In Fig. 3c, we plot the time-dependence of local stabilizers.In the SPT time crystal region, we observe that the stabilizers at the boundaries oscillate with a 2T periodicity, indicating again a spontaneous breaking of the discrete time translational symmetry at the boundaries.In the bulk, the stabilizers oscillate with a T periodicity and are synchronized with the driving frequency, thus no symmetry breaking occurs.This is in sharp contrast to the dynamics of bulk magnetizations, which decays rapidly to zero and exhibits no oscillation.In fact, in the SPT time crystal region the system is MBL and there exist a set of local integrals of motion, which are the "dressed" version of the stabilizers with exponentially small tails [26].The persistent oscillations of the bulk stabilizers observed in our experiment originates from these local integrals of motion and is a reflection of the fact that the system is indeed in a MBL phase.

Phase transition
We now turn to the phase transition between the SPT time crystalline phase and the trivial thermal phase.For simplicity and concreteness, we fix other parameters and vary the drive perturbation δ and strength of the interactions V .For the chosen parameter regions, theoretically the system would exhibit a SPT time crystalline phase with small δ and V .Whereas, with increasing δ and V the strong interaction diminishes localization and eventually thermalizes the system.At some critical values of δ and V , a transition between these two phases occurs.In Fig. 4a, we plot the δ-V phase diagram obtained from numerical simulations, where the phases boundary, although not very sharp due to the finite-size effect, can be located and visualized approximately.
To experimentally examine this phase transition, we further fix the interaction strength V = 0. We probe the transition point by measuring the variance of the subharmonic spectral peak height, i.e., Fourier spectrum of σ z 1 at ω = ω 0 /2 for the boundary spin.Figure 4b shows the subharmonic peak height as a function of the drive perturbation δ.With small δ, the system is in the SPT time crystalline phase and the peak height remains at a notable value about 0.5.As we increase δ to a large value, the crystal 'melts' and peak height vanishes.This confirms the theoretical analysis above.The largest variance of the peak height corresponds to the phase transition point.The inset of Fig. 4b shows the measured standard de-

Conclusion and outlook
In summary, we have experimentally observed signatures of a SPT time crystal with a programmable superconducting quantum processor.In contrast to previously reported conventional time crystals, for our observed SPT time crystal the discrete time translational symmetry only breaks at the boundries but not in the bulk.We measured the persistent oscillations of edge spins with a subharmonic frequency and experimentally demonstrated that the SPT time crystal phase is robust against to perturbations in the drive and imperfections in the experiment.In addition, we also demonstrated that subharmonic response of the boundary observables is independent of the initial states.
The controllability and scalability of the superconducting platform demonstrated in our experiment open up several new avenues for future fundamental studies and potential applications.In particular, it would be interesting and important to explore other exotic non-equilibrium phases beyond classical simulability, with the superconducting or other platforms.In practice, the observed SPT time crystal may have applications in some quantum information processing tasks, such as quantum metrology or implementing a robust quantum memory with topological protection.
Supplementary Materials: Observation of a symmetry-protected topological time crystal with superconducting qubits I. THEORETICAL UNDERSTANDING A. Introduction to the conventional Floquet time crystals In order to obtain a better intuitive understanding of the symmetry-protected topological (SPT) time crystal, we first introduce the basic concepts of the traditional Floquet time crystals and give a prototypical model as a concrete example.
Spontaneous symmetry-breaking is an important concept in modern physics.This occurs when the steady states of a physical system do not respect the symmetries of the Hamiltonian governing this system.An important example manifests spontaneous symmetry-breaking is the ordinary crystals, which break the continuous spatial translational symmetry.More precisely, in a crystal, the system is not invariant under continuous translation operators, which is respected by the system's Hamiltonian.Analogously, systems that spontaneously break the time translational symmetry are named as time crystals [1,49].Although there is a no-go theory for the continuous time crystals Bruno [50], Floquet time crystals manifests themselves in many physical systems.There are two equivalent definitions of Floquet time crystals in Ref [2], which characterize this concept from the perspective of the expectation value of an operator and from the perspective of the eigenstates of the Floquet evolution unitary respectively.The first definition states that time translational symmetry breaking occurs if for every state |ψ(t) at arbitrary time t, there exists an operator O satisfying that ψ(t , where |ψ(t + T ) = U F (T ) |ψ(t) with U F (T ) being the Floquet evolution unitary of one period T .This definition implies how to observe time crystals experimentally and is used in our paper.The second definition states that the time translational symmetry breaking occurs if all eigenstates of the Floquet evolution unitary are long-range correlated.This concept is used in our theoretical analysis.
To be more concrete, we introduce the following prototypical time-dependent Hamiltonian previously studied in Ref [2,3] as an example of Floquet time crystals: where J k and h z k are uniformly chosen from For simplicity, we choose T = 2T 1 = 2.The Floquet evolution operator in one period can be written as U F = exp(−iH 2 ) exp(−iπ/2 k σx k ).We consider the eigenstates of H 2 , which are the product states polarized in z direction |Θ = |{s k } with s k = ±1, and these states can be prepared in experiments.
As evolution unitary in the first time interval has the effect of flipping all spins, the state |Θ is related to another state |−Θ = |{−s k } , which is also an eigenstate of H 2 .As for the time evolution in the second time interval, operating unitary e −iH2 on the two states, we will have are the expectation value of the first and the second part of H 2 for the state |Θ .Therefore, under Floquet operator U F , we have relations: From the above equations, we see the Floquet evolution unitary U F mixes the two states |Θ and |−Θ .In the subspace formed by |±Θ , U F has the matrix form: Diagonalize the matrix, we can get eigenvalues ± exp(−iE + (Θ)), and eigenstates |Θ ±exp Obviously, the eigenstates of U F are paired cat states with long-range correlations.Thus, in the localized region, this model satisfies the second definition of a Floquet time crystal in Ref. [2] and discrete time translational symmetry breaking could occur in such a system.Besides, as the Floquet operator has eigenvalue ± exp(−iE + (Θ)), if we diagonalize the effective Hamiltonian of the Floquet operator, we can get two eigenvalues with quasi-energy difference π.This is just the π-spin glass phase introduced in [4].

B. Our model: the SPT time crystal
Different from the model mentioned above, our SPT time crystal only breaks the discrete time translational symmetry at the boundaries.Precisely speaking, our model manifests the subharmonic response at frequency 2π/2T only at the edges, but not in the bulk of the system.Here T is the period of the Floquet driving.Now, let us give more theoretical analysis of this model.

Localized and SPT quantum states
Our SPT time crystal has two distinct governing Hamiltonians in different time intervals as showing in the main text.In the first time interval, this governing Hamiltonian is the sum of one-body Pauli operators on different sites.And in the second time interval, the governing Hamiltonian includes interaction among neighboring sites, which introduces the subtle many-body properties in this system.Let us begin with the static Hamiltonian H 2 [26]: where the parameters are chosen as explained in the main text.This Hamiltonian has a Z 2 × Z 2 symmetry, corresponding to σz,y In the extreme case V k = h k = 0, the eigenstates of this Hamiltonian are the mutual eigenstates of all stabilizers, which are named as cluster states.They are SPT states with Z 2 × Z 2 symmetry.The SPT phase manifests itself in the open boundary case: there is one effective free spin at each end of the chain.The topological property of those states is encapsulated by the string-order parameter: which takes random values O st (l, j) = ±1 for different eigenstates and different disorder realizations.Thus, we can define a non-local analogue of the Edwards-Anderson glassorder parameter to characterize the SPT time crystal phase: where [[]] denotes an average over sites, states and random realizations.Furthermore, the entanglement spectra of the SPT states are degenerate.This degeneracy can serve as another manifestation of the topological property.Besides, in this case, all energy levels are exactly four-fold degenerate.The corresponding degenerate eigenstates can be divided into four groups: Here we are working in the σz basis for simplicity, and the two arrows represent the effective boundary spins and the ... denotes the bulk spins.They are related by σx the two body terms and the one-body terms make the eigenstates of this Hamiltonian departing from the cluster states.However, if we keep this Hamiltonian deep in the topological phase (the region we are interested in): we can also interpret this model from a manybody localized (MBL) perspective.Different from the Hamiltonian with strictly localized stabilizers being integrals of motion, in this region, the system obtains another set of mutually commuting quasi-local integrals of motion.Each of them has a finite interaction length, and thus, is not strictly local.Similarly, in open boundary condition, there exists a quasi-local effective free spin at each edge, which contains the bulk components decaying exponentially with the distance.In this case, the string-order parameter and the degeneracy of the entanglement spectra can also manifest the topological property of the corresponding eigenstates.Moreover, while the energy spectra stop to be exactly four-fold degenerate for a finite system For the edge spins, it is clear that σ z 1 and σ z 100 displays persistent oscillations with a 2T periodicity, manifesting the breaking of discrete time-translational symmetry.In stark contrast, in the bulk region σ z k decays rapidly to zero and no symmetry breaking is observed.This is the defining feature of the SPT time crystals: time-translational symmetry only breaks at the boundary, not in the bulk.b.Fourier spectra of σ k .We find that σ1(ω) has a peak at ω/ω0 = 1/2, where ω0 = 2π/T is the driving frequency.We mention that this peak is robust and rigidly locked at ω0/2, a manifestation of the robustness of the SPT time crystals.For the bulk spins, there is no such peak, consistent with no symmetry breaking in the bulk.c.Logarithmic entanglement entropy growth.In the SPT time crystal region, the system is many-body localized.We thus expect a logarithmic entanglement growth, which is showing in this figure.We note that the entanglement has a initial quick rise till time Jt ∼ 1, corresponds to expansion of wave packets to a size of the order of the localization length.size, it is still nearly four-fold degenerate and the corresponding eigenstates can be divided into four groups as mentioned above.

The emergence of the SPT time crystal
With those understandings of the static Hamiltonian H 2 , now let us consider the Floquet case, wherein we periodically drives the above SPT Hamiltonian as discussed in the paper: Let us begin with the perfect case, where λ = π/2, δ = 0 and V k , h k = 0.So the energy spectrum of H 2 is perfectly four-fold degenerate.The eigenstates can be divided into four groups, i.e.
There exists a local effective free spin at each boundary.Since the effect of U 1 = e −iπ/2 k σx k is perfectly flipping the spins at all sites, we can obtain the following relations of the Floquet operator U F = exp (−iH 2 ) exp (−iπ/2 k σx k ): From this, we see that U F mixes the states |A k with |D k , and mixes the states |B k with |C k .Writing in the subspace of |A k and |D k , U F has the matrix form: In the subspace of |B k and |C k , U F has the matrix form: Therefore, in the subspace formed by Thus, the floquet effective Hamiltonian has eigen-energies: As the energy spectrum is four-fold degenerate: Therefore, the original four-fold degeneracy breaks into two-fold degeneracy after the drive.This two-fold degeneracy is the reminiscence of the original topological order.As for the Floquet eigenstates, by diagonizing U F in the two subspaces, we will get that these Floquet eigenstates are cat-like linear combination of topological eigenstates: |A k ± |D k , and |B k ± |C k .These states are long-range correlated, and thus, the mutual information between the two boundaries is 2 log 2. The long-range correlated property of these eigenstates indicates that discrete time translational symmetry breaking can occur in this system, implying a Floquet time crystal.
When we turn on the two-body terms and the one-body terms, but still constrain the Hamiltonian H 2 deep in the topological region: J k h k , V k , it has four nearly degenerate eigenstates related by the symmetry operations for a finite system size.The effective free spin in each boundary becomes quasi-local.Under the Floquet driving, the nearly four-fold degeneracy breaks into nearly two-fold degeneracy: Following the similar discussion, we can find that the eigenstates of the Floquet unitary are also cat-like states, and thus, time translational symmetry breaking can occur in this case.This discussion makes it clear that our model satisfies the definition of Floquet time crystals in this case, as long as the governing Hamiltonian in the second time interval is also in the deep topological region.The stability of the SPT time crystal will be discussed more precisely in the later subsection.

Dynamical properties of the SPT time crystal
Next, we will consider the evolution of this system, and explicitly manifest the behavior of the SPT time crystal.
Let us start from a product state: |ψ 0 = |↓ ... ↑ , here ... denotes the product states of bulk spins.Because the boundary spins fall into the group of {|B k }, we can expand the initial state as |ψ 0 = k b k |B k .Under the time evolution of U F for one driving period: where |B k = |↓ ... ↑ , |C k = |↑ ... ↓ .So if we measure the edge spins at the starting point, we have Similarly, after two Floquet periods, the state becomes Thus, we see the edge spins exhibit breaking of the time-translational symmetry.
As for the bulk spins, assume that one bulk spin σz j has expectation value for the initial product state: ψ 0 | σz j |ψ 0 = 1.Writing |ψ 0 in the |B k basis, we have: Since the spins of |C k are opposite to the spins of |B k at all sites, we immediately have that : However, the expectation value of σz j in the state |ψ 1 should be expressed as: Comparing the two above equations, we see that because of the extra phase factor exp −iE C k + iE C k before each component, the σz j will not have definite value after the Floquet time evolution and will decay to zero quickly after random averaging.Thus, the bulk spins do not exhibit breaking of the time-translational symmetry.
The above derivations tell that for our model, with the time evolution, the edge spins exhibit discrete time-transnational symmetry breaking, while the bulk spins relax very fast.Thus, the time translational symmetry breaking only occurs at the boundaries as showing in Fig. S1a.We stress the importance of topology here.It protects the edge spins, ensuring the robustness of the edge spins against local perturbations that respect the underlying symmetry.
Besides, for a localized system, the entanglement entropy will grow logarithmically [51].In the deep SPT time crystal region, the system represented by the static Hamiltonian H 2 has a series of quasi-local integrals of motion.So, spins far away from each other can build significant entanglement only after exponentially long evolving time.Thus, under the Floquet time evolution, the entanglement entropy of this system will exhibit a logarithmic growth, which is showing in Fig. S1b, and will finally saturate to a value proportional to the system size.
Furthermore, in the deep topological region, when the system size is finite, the Floquet time evolution can lead the states of the system to the thermal states after long enough evolving time.Actually, the quasi-local effective free spins at the boundaries have exponentially decaying tails, containing the bulk components.When the system size is finite, the tails of the quasi-local effective free spins have an exponentially small overlap, which will relax the two effective free spins.Due to the fact that the overlap is exponentially small with the system size, the lifetime of the time translational symmetry breaking states will exponentially diverge with system size.This phenomenon is demonstrated in our numerical simulation Fig. S2a.The omitted parameters are chosen as in Fig. S1.We find there is an initial quick decay of the σ z 1 , followed by a plateau that extends up to a time diverging exponentially with the system size.The inset shows the exponential scaling of τ * with the system size, where τ * is the time when the spin at edge decays to 1/2.b.The phase diagram of the SPT time crystal with respect to the imperfect drive δ and the average strength of the interactions V .Here, we adapt the string order parameter Osg (averaged over 100 random realizations) as the indicator.It shows that when the imperfections are not very large, the string order parameter is approximately equal to one, indicating the topological property.The other parameters are chosen as J = ∆J = 1, h = ∆h = 0.01.

The stability of the SPT time crystal
The above considerations rely on the fact that during each period, we perfectly flip spins at all sites.To show that this is indeed a time crystal, we should make sure this phase is stable even for imperfect drive.We follow the similar arguments introduced in Ref. [2].
In above subsection, we already showed that in the perfect drive case, the eigenstates of the Floquet evolution operator are cat-like states ψ AD are all long-range correlated, with different eigen-energies.Then for any experimentally prepared short-range states (such as product states), they can only be formed by taking superposition of those long-range Floquet eigenstates with different energies.Thus, after one period of Floquet evolution, local observables at the edge will not be invariant and has to break the discrete time translational symmetry.Now we add local perturbations into the system, such as imperfect drives and other local interactions.Since the system is in the MBL phase, a local perturbation will only significantly affect those near sites.Thus, we expect that the long-range correlations in the eigenstates of the Floquet unitary will not disappear.Actually, there exists one quasi-local unitary operator U , which relates the perturbed Floquet eigenstates with the un-perturbed Floquet eigenstates.Since U is a quasi-local operator, it cannot destroy the long-range correlations of the un-perturbed Floquet eigenstates.Therefore, time translational symmetry breaking can also occur in the locally perturbed system.To explicitly show that the SPT time crystal is indeed a physical phase, we use the string order parameter O sg as the indicator to plot the phase diagram with respect to the imperfect drive (δ) and the average strength of the two-body interactions (V ), which is showing in Fig. S2b.

II. DETAILS OF TEBD METHOD
We numerically simulate the time evolution process of the SPT time crystals using the time-evolving block decimation (TEBD) methods.This method was proposed for the time evolution of the matrix product states (MPS) [52,53] and is variant of the density matrix renormalization group (DMRG) algorithm [54,55].At the heart of the TEBD method lies the Trotter-Suzuki decomposition of the time evolution operator U (∆t) of a short-range interaction system in a small time interval ∆t.Usually, we can represent the operator U (∆t) in the matrix product operator (MPO) form with small Trotter error, and then repeatedly apply it on the MPS representing current quantum state |ψ(t) to implement the time evolution.
Our SPT time crystal has two distinct Hamiltonian operators in different time intervals as showing above.For the first time interval, the corresponding Hamiltonian is the sum of one-body operators on different sites.So the evolution operator is a direct product of one-body evolution operators which can be represented as MPO directly.To obtain the corresponding expectation values of some local observables at different times, we also decompose the time evolution operator of an entire time interval T 1 into several small time intervals ∆t.We show the implementation of U 1 (∆t) in Fig. S3a.
For the second time interval, the Hamiltonian H 2 consists of multiple short-range interaction terms: . Thus, we can approximate the time evolution operator using Trotter-Suzuki decomposition U 2 (t) ≈ [U 2 (∆t)] t/∆t = e −i∆tH2 t/∆t with ∆t t.To efficiently construct the MPO representation of U 2 (∆t), we group together terms in H 2 which are commuting with each other.For the three-body operators, all of them are stabilizer operators and they are commuting with each other.For the two-body terms, we can group terms without overlap and construct corresponding MPO operators.For one-body terms, all of them are operators on different sites and thus, they are commuting with each other.For simplicity, we denote imation as following Thus, the evolution operator of H 2 of time interval t has ap-proximation t/∆t = e −i∆tD e −i∆tC e −i∆tB e −i∆tA t/∆t + O(∆t).
Furthermore, considering the efficiency, the implementation of three-body terms can be accomplished layer by layer, wherein each layer only contains three-body operators with no overlap with each other, so that they can be applied to the MPS parallelly.We emphasize that since the Trotter error is of order ∆t, the time interval ∆t should be small enough to avoid large TEBD error.The implementation of U 2 (∆t) is showed in Fig. S3b.The circuit ansatz for the time evolution unitary in the second time interval, where the system is in deep topological phase.c.The circuit ansatz for the time evolution unitary in the second time interval, where the system contains no two-body operators and no one-body operators: Algorithm 2: Optimization for a quantum circuit Output: Optimal parameters of the given quantum circuit.Input: A quantum circuit C, evolution unitary U 2 (∆t) and learning rate α.Randomly initialize θ; To observe the SPT time crystal on a real superconducting quantum computer, we need to decompose the time evolution unitary into a quantum circuit consisting of a series of experimentally implementable quantum gates.Due to the direct product form of the evolution unitary in the first time interval U 1 (t) = e −itH1 , this unitary can be represented as a quantum circuit using a layer of rotation gates along x axis.Thus, it can be constructed and implemented relatively easily.As for the second time interval, the interaction among different sites leads the time evolution unitary far away from a direct product form, making things a little different.
With the progress of the research on the variational quan-tum circuits, we are able to adapt this method to construct the quantum circuit of the second time evolution unitary in our time crystal model.Variational quantum circuit is a powerful tool and is intensively investigated in recent years.Algorithms based on variational quantum circuits are of great potential in the noisy intermediate scale quantum (NISQ) There are many algorithms based on variational parameterized quantum circuits, such as variational quantum eigen-solver [56], quantum neural network [57], etc.They are proposed to solve practical problems.The major distinction between usually used quantum circuits and variational quantum circuits is that the gates contained in a variational quantum circuit are not fixed.They can be modified by switching their parameters using different parameter-updating algorithms.With the updating of those parameters, the unitary implemented by the variational circuit is varying.We end this updating procedure until a satisfactory result is obtained.
Here, our target is to find a variational quantum circuit with some fine-tuned parameters nearly equivalent to the evolution unitary in the second time interval U 2 (t) = e −itH2 .We accomplish this target within two steps: finding an implementable variational quantum circuit ansatz which can be used to represent the target unitary, and updating those parameters contained in this circuit ansatz according to the evolution unitary of a particular random realization of the hamiltonian.
We use the neuroevolution method [58] to find a suitable variational circuit architecture.The elementary gates used in our experiments are single qubit rotation gates (X, Y, Z) and control-rotation gate along z axis.And each of them contains a variational parameter.They can form various quantum circuit layers.i.e. quantum circuits with depth equal to one.We construct the direct graph of those layers consisting of those elementary gates with respect to the rules in [58].Then, a quantum circuit can be represented as a path in this graph.The later problem can be solved with the following procedures: 1) Randomly generate several variational quantum circuits with fixed length based on the directed graph; 2) Update parameters contained in those quantum circuits using gradient based algorithm to minimize the loss function L(θ) = 1 − Tr U 2 (∆t) † U circuit (θ) /d, where U 2 (∆t) is the evolution unitary in the second time interval, U circuit (θ) is the unitary represented by the current quantum circuit with variational parameters θ, and d is the dimensionality of the corresponding Hilbert space; 3) Chose quantum circuits with small loss values and extend them based on the direct graph to generate new circuits; 4) Iterate processes 2) and 3) until the loss is lower than a fixed threshold.The circuit ansatz giving the smallest loss value is regarded as the optimal ansatz representing the evolution unitary and is adapted in our experiments.We show the pseudo-code of this algorithm in Algorithm 1.
The quantum circuit ansatzs used in our experiments are showing in Fig. S4.We notice that the quantum circuit for the evolution unitary in the second time interval has a sandwich form U 2 (∆t) ≈ W D(θ)W † , where D(θ) is a layer of single qubit rotation gates with θ being the evolution-time dependent parameters.So for one driving period, the depth of the corre-sponding quantum circuit is fixed, i.e. for n driving periods, the depth of the corresponding quantum circuit is 6n, which is independent with the number of time slices in a period.
With this circuit ansatz in hand, we can then use it to construct the experiment circuits.For a particular disorder realization of the second hamiltonian deep in the topological phase (J k V k , h k ), we begin with this ansatz containing randomly generated variational parameters θ.From this, the gradients of the loss function L(θ) with respect to those variational parameters are computed, and are used to update the current parameters θ (n+1) = θ (n) − α∇ θ (n) L, where α is the given learning rate (we usually chose 0.001 ≤ α ≤ 0.01).
In our calculation, we iterate this optimization procedure until the operator fidelity [59] Tr U 2 (∆t) † U circuit (θ) /d ≥ 0.999 (L(θ) ≤ 0.001).Then, we take the quantum circuit with the final parameters as the approximation of the evolution unitary U 2 (∆t) in our experiments.We show the pseudo-code of this algorithm in Algorithm 2.
We emphasize that this optimization procedure is suitable for small systems.Nevertheless, with the exponential growth of the dimensionality of the Hilbert space, the optimization for large systems is unpractical.It is helpful that the quantum circuit ansatz found by the neuroevolution method can exactly represent the evolution unitary U 2 (t) = e −itH2 , when H 2 consists no two-body operators and no one-body operators (as showing in Fig. S4c).This indicates that we can theoretically construct the corresponding quantum circuits for arbitrary many qubits.In our simulations and experiments, for systems of L ≤ 8, the two-body terms and one-body terms are considered and the parameters in the corresponding quantum circuits are obtained using gradient-based optimization method.For 14-qubit systems, we only consider the stabilizer terms and exactly construct those corresponding quantum circuits.

B. Device overview and measurement setup
To illustrate the idea of SPT time crystal, we select a chain of up to L = 14 qubits in a superconducting quantum processor, which is a flip-chip device hosting an array of 6 × 6 qubits distributed in a square lattice.To realize high-fidelity controlled-Z (CZ) gates, we adopt the tunable coupler architecture [36] to mediated nearest neighbor qubit-qubit interactions, i.e., individual couplers are inserted between neighboring qubits with the qubit-coupler coupling strengths designed to be around 130 MHz for qubits at 6.5 GHz.All qubits (couplers) are of transmon type, with anharmonicities around 250 (350) MHz and maximum resonance frequencies around 7 (10.5)GHz.Each qubit has its own control line, which takes microwave (XY) inputs for rotating the qubit state around xor y-axis and flux bias (Z) pulses for tuning the qubit frequency and rotating the qubit state around z-axis; each coupler is frequency tunable via its own flux bias (Z) line, which guarantees that the effective coupling strength between two neighboring qubits at 6.5 GHz can be dynamically turned on, up to −25 MHz, or off, ≤ 0.25 MHz.Each qubit capacitively couples to its own readout resonator, designed in the frequency range from 4.1 to 4.4 GHz, for qubit state measurement.9 readout resonators share one readout transmission line (TL) running across the processor chip, and 4 readout TLs can cover all 36 qubits in the processor.
processor was fabricated using the flip-chip recipe: all qubits and couplers are located on sapphire substrate (top chip); most of the control/readout lines and readout resonators are located on silicon substrate (bottom chip).These two chips have lithographically defined base wirings, junction loops, and airbridges made of aluminum, and are galvanically connected via indium bumps with titanium under-bump metallization, as described elsewhere [60].The indium bumps were formed by the lift-off method with 9 µm-thick indium deposited on both chips, after which these two chips were aligned and bonded together at room temperature to complete the flip-chip device.The indium bumps in our processor are not only for ground connectivity, but also for passing through control signals from the bottom chip to the top chip where the qubits are located.
The processor was loaded into a multilayer printed circuit board (PCB) enclosure, which was then mounted inside a dilution refrigerator (DR) with the base temperature down to 15 mK. Figure S5 shows the schematics of the control/readout electronics and wiring setup.In this setup, the XY microwave signals and fast Z pulses synthesized by digital-to-analog converters (DACs) are first joined together at room temperature, which are attenuated and filtered at multiple cold stages of DR, and later combined with the slow Z (DC) pulses via homemade bias-tees at the mixing chamber stage of DR before being transmitted into the qubit control lines.The multiplexed readout signals are also heavily attenuated and filtered before going into the readout TLs of the processor to retrieve the qubit state information.To boost the signal-to-noise ratio (SNR), output signals from TLs are sequentially amplified by Josephson parametric amplifier (JPA), high electron mobility transistor (HEMT) amplifier, and room temperature (RT) amplifiers before being demodulated by analog-to-digital converters (ADCs) with 10-bit vertical resolution and 1.0 GS/s sampling rate.All arbitrary microwave signals are generated by mixing the DAC outputs with continuous microwaves using IQ mixers.DACs used to synthesize XY microwave signals and fast Z pulses in this experiment has 14-bit vertical resolution and 300 MHz output bandwidth.Slow Z (DC) pulses are generated by commercial 16-bit DACs with maximum outputs of ±2.5 V.

C. Single-and Two-qubit gates
Single-qubit gates used in this experiment include X(θ), Y(θ), and Z(θ), which rotate the qubit state by an arbitrary angle θ around x-, y-, and z-axis, respectively.We realize X(θ) and Y(θ) by controlling the amplitude and phase of XY microwave pulses, and implement Z(θ) via the virtual Z gate [61].Single-qubit gate errors are characterized by simulatneous randomized benchmarkings, yielding an average gate fidelity above 0.99 (see Tab. S1).
The basic structure to implement the CZ gate consists of two flux tunable qubits and one flux tunable coupler, which are, respectively, noted as Q 1 , Q 2 , and C here for the clarity of description.The effective coupling strength is composed of a direct coupling strength between two qubits and a part mediated by the coupler, which can be continuously adjusted by controlling the flux or frequency of the coupler.The Hamiltonian of this three body system is written as The effective coupling strength between qubits is In Fig. S6a we plot the dynamic range of g (bottom panel) processed using the two-qubit swap dynamics after initializing Q 1 -Q 2 in |10 , which shows that the effective coupling strength is tunable in the range from −25 MHz to ≤ 0.25 MHz.Experimentally we can apply single-qubit gates while locating the coupler at around 10.5 GHz to turn off g.
To realize the CZ gate, we apply a flux bias (fast Z) pulse to steer the coupler's frequency in a trajectory of 10.5→7.3→10.5 GHz, and meanwhile we turn on the fast Z pulses to bring Q 1 and Q 2 from their idle frequencies to the pair of values in ω A(B) j , ω

A(B) j+1
depending on which pair we use, where |11 and |02 in the two-qubit subspace are in near resonance (see Tab. S1).After a finite period for this diabatic interaction, an unitary two-qubit gate equivalent to CZ up to trivial single-qubit phase factors can be obtained as A sine decorated square pulse with the amplitude A = z 0 × 1 − r + r sin π t tgate is used for the coupler in order to minimize state leakages.Experimentally we fix r = 0.3 and only fine-tune the parameter z 0 .All pulses are digitally smoothed by convolving them via a Gaussian window with σ = 2 ns before applying our pulse calibration routines [65].The CZ gate pulse duration is 30 ns, and there are additional 5 ns padding times before and after the 30-ns gate in compensation for the finite small tails of the smoothed pulse.
Individual CZ gates are calibrated following the procedure below: TABLE S1.Device parameters calibrated during the experiment.ω 0 j is the maximum frequency of Qj at zero flux bias.ωj is the idle frequency where we initialize Qj in |0 and subsequently apply single-qubit gates.ηj is Qj's anharmonicity, which is approximately a constant within the frequency range relevant to this experiment.ω A (B) j , ω A (B) j+1 list the pairwise frequency values for neighboring two qubits where |11 and |02 in the two-qubit subspace are in near resonance for a CZ gate; the CZ gates for qubit pairs in the same group A (B) are implemented simultaneously when executing the multilayer quantum circuit to simulate SPT time crystal.ω m j is the readout frequency of Qj where we apply readout pulses to excite Qj's readout resonator for quantum state measurement.ω r j is the resonant frequency of Qj's readout resonator.T1,j and T * 2,j are the energy relaxation time and Ramsey dephasing time of Qj, respectively.F0,j and F1,j are the readout fidelity values for Qj prepared in |0 and |1 , respectively, which are used to correct raw probabilities to eliminate readout errors as done previously [62].esq list the single-qubit gate errors obtained by simultaneous randomized benchmarkings.e A(B) CZ lists the CZ gate errors obtained by both individual and simultaneous randomized benchmarkings for qubit pairs in group A (B).We note that the qubit parameters may slowly drift over time [63,64] Unlike thermal phases without disorder or Anderson localized phases without interaction, where entanglement grows ballistically [66][67][68] or saturates to an area law at long time, the entanglement entropy of a MBL system grows logarithmically and saturates to an volume law in the long time limit [51].In our experiment, we also extract the entanglement dynamics, through a full quantum state tomography of the reduced density state of half the system.In Fig. S7a, we plot the reduced density matrix ρ half for a single random instance of the Hamiltonian at the end of one driving period.With the tomography data of the reduced density states at different times, we extract the desired information about the entanglement growth for a SPT time crystal.Our results are plotted in Fig. S7c.From this figure, it is clear that in the thermal region, the entanglement grows quickly and saturate to a maximal volume law (∼ L 2 ln 2) after about three driving period.In contrast, in the SPT time crystalline region, the entanglement grows much slower.Due to the decoherence and other imperfections in our experiment, we are not able to observe the logarithmic growth of the entanglement.This not only demands a significant improvement of the gate fidelity and a substantial increase of the coherence time, but also a more efficient and scalable approach to measure entanglement for a many-body system.
FIG.1.Symmetry protected topological time crystal and schematics of the experimental setup.a, Fourteen qubits used in our experiment are coupled to their neighbors with capacitive couplers.b, A chain of spins are periodically driven with the stroboscopic Floquet Hamiltonian H(t), giving rise to a SPT time crystal characterized by spontaneous time translational symmetry breaking at the boundaries.c, The quasienergy spectrum ( ) of the Floquet unitary, which is the time evolution operator over one period.For the SPT time crystal, every eigenstate is twofold degenerate and has a cousin separated by quasienergy π. d, A schematic illustration of the experimental circuits used to implement the time dynamics governed by the Floquet Hamiltonian H(t).We randomly sample the Hamiltonians and prepare the initial states onto random product states or random SPT states.After running a sequence of quantum gates, we measure the local magnetization for each qubit or stabilizers at discrete time points.e, Illustration of the quantum processor, with the 14 qubits highlighted in green.

FIG. 2 .
FIG.2.Observation of a SPT time crystal.a, Time evolution of disorder-averaged local magnetizations deep in the SPT time crystal region (J = ∆J = 1, V = h = ∆V = ∆ h = 0, and δ = 0.01).The initial states are |0 ⊗L and data shown are averaged over 20 random disorder instances.While the bulk magnetization decays quickly to zero, the edge spins oscillate with a stable subharmonic response up to 20 cycles.b, The evolution dynamics of local magnetizations for different random instances.Here, each layer corresponds to a specific random instance.c, Magnetization dynamics deep in the thermal region (J = ∆J = 1, V = h = ∆V = ∆ h = 0, and δ = 0.8).d, Fourier transform of experimentally measured σ z j (t) in the time crystal region.The edge spins lock to the subharmonic frequency, which is in sharp contrast to that for bulk spins.e, Fourier spectra of σ z j (t) in the thermal region.No robust subharmonic frequency peak appears for both edge and bulk spins in this case.f, Time-dependence of the autocorrelator Aj = Ψ0|σ z j (t)σ z j (0)|Ψ0 up to 40 cycles, obtained from averaging over 20 random instances deep in the SPT time crystalline phase, with the initial states prepared at random product states |{0, 1} ⊗L .The black solid lines show the results of "echo" circuits for two boundary qubits.

8 FIG. 3 .
FIG.3.Dynamics of stabilizers with random initial SPT states.a, Schematic of the experimental circuit for preparing random SPT states.To prepare the system into the ground state of the cluster Hamiltonian H 2 , we apply a Hardmard gate (H) on each qubit and then run CZ gates in parallel for all neighboring qubit pairs in two steps.Then we act Z(π) on some random sites to create excitations and transfer the ground state to a highly excited eigenstate of H 2 .This procedure enables the preparation of random SPT states at high energy.We then evolve these states with the Floquet Hamiltonian H(t) to study the dynamics of stabilizers.b, Entanglement spectra of random SPT states prepared in our experiment, with both open and periodical boundary conditions.The dashed lines indicate the lowest degenerate entanglement energies obtained by simulation in error-free case.The two-and four-fold degeneracy at low entanglement energy is a characterizing feature for the topological nature of these states.c, The time dependence of stabilizers in the SPT time crystal region, averaged over 20 random circuit instances.

FIG. 4 .
FIG.4.Phase diagram and detection of phase transition.a, The numerical δ-V phase diagram obtained by examining the central subharmonic peak height for the edge spins in the Fourier spectrum, averaged over 1000 disorder instances.The dashed line corresponds to the maximal height variances for varying V with each fixed δ point, which gives a rough estimation of the phase boundary.Here, the parameters are chosen as L = 8, J = ∆J = 1, and h = ∆ h = ∆V = 0.05.b, Experimental result of the subharmonic peak height as a function of the drive perturbation δ with fixed h = V = ∆ h = ∆V = 0 and J = ∆J = 1, averaged over 20 disorder instances uniformly sampling from the interval [J −∆J , J +∆J ] for 14 qubits, with the shadow outlining the standard deviation.Inset, the standard deviation of the central peak height as a function of δ, which signals a phase transition from the SPT time crystalline phase to a thermal phase.
FIG. S1.The evolution of the SPT time crystal with system size L = 100, computed using the time-evolving block decimation methods.Results showed here are averaged over 1000 random realizations, with parameters J = ∆J = 1, h = ∆h = V = ∆V = δ = 0.05.a.Time evolution of disorder-averaged local observables.For the edge spins, it is clear that σ z 1 and σ z 100 displays persistent oscillations with a 2T periodicity, manifesting the breaking of discrete time-translational symmetry.In stark contrast, in the bulk region σ z k decays rapidly to zero and no symmetry breaking is observed.This is the defining feature of the SPT time crystals: time-translational symmetry only breaks at the boundary, not in the bulk.b.Fourier spectra of σ k .We find that σ1(ω) has a peak at ω/ω0 = 1/2, where ω0 = 2π/T is the driving frequency.We mention that this peak is robust and rigidly locked at ω0/2, a manifestation of the robustness of the SPT time crystals.For the bulk spins, there is no such peak, consistent with no symmetry breaking in the bulk.c.Logarithmic entanglement entropy growth.In the SPT time crystal region, the system is many-body localized.We thus expect a logarithmic entanglement growth, which is showing in this figure.We note that the entanglement has a initial quick rise till time Jt ∼ 1, corresponds to expansion of wave packets to a size of the order of the localization length.
FIG. S2.The decay of spins at boundaries and the phase diagram of the SPT crystal.a.The decay of the random realization averaged magnetization for the spin at the edge.Here, the number of the disorder realization ranges from 3 × 10 4 (L = 6) to 10 3 (L = 14).The omitted parameters are chosen as in Fig.S1.We find there is an initial quick decay of the σ z 1 , followed by a plateau that extends up to a time diverging exponentially with the system size.The inset shows the exponential scaling of τ * with the system size, where τ * is the time when the spin at edge decays to 1/2.b.The phase diagram of the SPT time crystal with respect to the imperfect drive δ and the average strength of the interactions V .Here, we adapt the string order parameter Osg (averaged over 100 random realizations) as the indicator.It shows that when the imperfections are not very large, the string order parameter is approximately equal to one, indicating the topological property.The other parameters are chosen as J = ∆J = 1, h = ∆h = 0.01.

±=
|A k ± |D k , and ψ BC ± = |B k ± |C k .We say that an effective short-range correlated topological state is that ψ| σ1 σN |ψ − ψ| σ1 |ψ ψ| σN |ψ → 0. Obviously, |A k , |B k , |C k , |D k are short-range correlated topological states, but the Floquet eigenstates ψ AD ± and ψ BC ± FIG. S3.Pictorial illustration of the implementation of the time evolution unitary, where a connected wire between different blocks means contraction of indices.a. Implementation of U1(∆t).The gray blocks represent the current state in MPS form and the green blocks represent the time evolution unitary consisting of one-body operators (α = −(π/2 − δ)).b.Implementation of U2(∆t).The gray blocks represent the current state in MPS form; the orange blocks represent the time evolution unitary consisting of three-body operators arranged in three groups; the blue blocks represent the time evolution unitary consisting of two-body operators arranged in two groups and the green blocks represent the time evolution unitary consisting of one-body operators (α = h k ).They are applied to the the current state layer by layer.
FIG. S4.Quantum circuit ansatz used in experiments.a.The circuit ansatz for the time evolution unitary in the first time interval.b.The circuit ansatz for the time evolution unitary in the second time interval, where the system is in deep topological phase.c.The circuit ansatz for the time evolution unitary in the second time interval, where the system contains no two-body operators and no one-body operators: U2 = exp (i k J k σz FIG. S6.Two-qubit CZ gate.a, two-qubit swap dynamics showing the dynamical range of the effective coupling strength as tuned by the coupler.b, |1 -state population landscape for Q1 resulting from the |11 and |02 interaction after initializing Q1-Q2 in |11 .The white star marks the vicinity of the gate parameters used for CZ and lines indicate how we sweep Z parameters to approach this vicinity.
FIG. S7.Entanglement dynamics.a, Tomography of the reduced density matrix for the second half of a six-qubit chain after one driving period of one instance.The real and imaginary parts are shown in the left and right panel, respectively.b, Entanglement growth in both the thermal and time crystal regions, averaging from 9 instances.