Observation of a Floquet symmetry-protected topological phase with superconducting qubits

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 crystals [1–8], where time translational symmetry is spontaneously broken in periodically driven systems. Pio-neering experiments have observed signatures of time crystalline phases with trapped ions [9, 10], spins in nitrogen-vacancy centers [11–13], ultracold atoms [14, 15], solid spin ensembles [16, 17], and superconducting qubits [18–20]. Here, we report the observation of a distinct type of intrinsically non-equilibrium state of matter, a Floquet symmetry-protected topological phase, which is implemented through digital quantum simulation with an array of programmable superconducting qubits. Unlike the discrete time crystals reported in previous experiments, where spontaneous breaking of the discrete time translational symmetry occurs for local observables throughout the whole system, the Floquet symmetry-protected topological phase observed in our experiment breaks the time translational symmetry only at the boundaries and has trivial dynamics in the bulk. More concretely, we observe robust long-lived temporal correlations and sub-harmonic temporal response for the edge spins over up to 40 driving cycles using a circuit whose depth exceeds 240. We demonstrate that the sub-harmonic response is independent of whether the initial states are random product states or symmetry-protected topological states, and experimentally map out the phase boundary between the time crystalline and thermal phases. Our work paves the way to exploring novel non-equilibrium phases of matter emerging from the interplay between topology and localization as well as periodic driving, with current noisy intermediate-scale quantum processors [21] . Here we provide more details on the theory of the Floquet symmetry-protected topological phase (Sec. I), on our numerical simulations (Sec. II), and on our experimental setup (Sec. III). We also provide additional experimental data (Sec. IV).

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 crystals [1][2][3][4][5][6][7][8], where time translational symmetry is spontaneously broken in periodically driven systems. Pioneering experiments have observed signatures of time crystalline phases with trapped ions [9,10], spins in nitrogen-vacancy centers [11][12][13], ultracold atoms [14,15], solid spin ensembles [16,17], and superconducting qubits [18][19][20]. Here, we report the observation of a distinct type of intrinsically non-equilibrium state of matter, a Floquet symmetry-protected topological phase, which is implemented through digital quantum simulation with an array of programmable superconducting qubits. Unlike the discrete time crystals reported in previous experiments, where spontaneous breaking of the discrete time translational symmetry occurs for local observables throughout the whole system, the Floquet symmetry-protected topological phase observed in our experiment breaks the time translational symmetry only at the boundaries and has trivial dynamics in the bulk. More concretely, we observe robust long-lived temporal correlations and sub-harmonic temporal response for the edge spins over up to 40 driving cycles using a circuit whose depth exceeds 240. We demonstrate that the sub-harmonic response is independent of whether the initial states are random product states or symmetry-protected topological states, and experimentally map out the phase boundary between the time crystalline and thermal phases. Our work paves the way to exploring novel non-equilibrium phases of matter emerging from the interplay between topology and localization as well as periodic driving, with current noisy intermediate-scale quantum processors [21].
Symmetry-protected topological (SPT) phases are characterized by non-trivial edge states that are confined near the boundaries of the system and protected by global symmetries [22][23][24][25][26]. In a clean system without disorder, these edge states typically only occur for the ground states of systems with a bulk energy gap. At finite temperature, they are typically destroyed by mobile thermal excitations. However, adding strong disorder can make the system many-body localized (MBL) [27][28][29][30][31], allowing for a sharply defined topological phase and stable edge states even at infinite temperature [32][33][34][35][36]. Strikingly, the topological phase and corresponding edge states can even survive external periodic driving, as long as the driving frequency is large enough so that the localization persists [37,38].
The interplay between symmetry, topology, localization, and periodic driving gives rise to various peculiar phases of matter that exist only out of equilibrium [38]. Understanding and categorizing these unconventional phases poses a notorious scientific challenge. On the theoretical side, topological classifications of periodically driven (Floquet) systems with [4, [39][40][41] and without [42] interactions have already been obtained through a range of mathematical techniques (such as group cohomology), revealing a number of "Floquet SPT" (FSPT) phases with no equilibrium counterparts [38]. 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, signatures of discrete time crystals (DTCs) [1][2][3][4][5][6][7][8], which are paradigmatic examples of exotic phases beyond equilibrium [43], 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 encompass topology as a key ingredient. A recent experiment simulating an FSPT phase on a trapped-ion quantum computer found that the phase was short-lived due to the presence of coherent errors in the device [44]. Realizing a long-lived FSPT phase, which demands a delicate concurrence of topology, localization, and periodic driving, thus still remains a notable experimental challenge.
Here, we report the observation of an intrinsically nonequilibrium FSPT phase with a programmable array of superconducting qubits ( Fig. 1)   , giving rise to an FSPT phase characterized by time translational symmetry breaking at the boundaries. c, The quasienergy spectrum ( ) of the Floquet unitary UF , which is the time evolution operator over one period. For the FSPT phase, every eigenstate is two-fold degenerate and has a cousin separated by quasienergy π. Here, |A ± |D and |B ± |C denote the eigenstates of UF . d, A schematic illustration of the experimental circuits used to implement the time dynamics governed by the time-periodic Hamiltonian H(t). We randomly sample the Hamiltonians and prepare the initial states as random product states or random static SPT states. After running a sequence of quantum gates, we measure the local magnetization or stabilizer operators at discrete time points. e, Illustration of the quantum processor, with the 14 qubits highlighted in green.
long coherence time. We successfully implement the dynamics of the system under a prototypical time-periodic Hamiltonian, which we theoretically predict to exhibit an FSPT phase. This requires us to digitally simulate a three-body-interacting Hamiltonian exhibiting an SPT phase via a large-depth quantum circuit obtained using a novel neuroevolution algorithm [45]. We then measure local spin magnetizations and their temporal correlations and demonstrate that both quantities show a subharmonic response at the boundaries but not in the bulk of the chain. This situation differs drastically from the case of DTCs, which exhibit subharmonic response everywhere in the bulk. This contrast stems from a fundamental distinction between DTC and FSPT phases: the former exhibit conventional long-range order in the bulk intertwined with the spontaneous breaking of discrete time-translational symmetry [43,46,47] while the latter exhibit SPT order that can only be revealed through boundary effects or nonlocal "string operators" in the bulk [39,41,48]. The observed boundary subhar-monic response persists over an extended range of parameters and is robust to various experimental imperfections, independent of the initial states. We further explore the FSPT phase experimentally from the perspectives of entanglement dynamics, entanglement spectrum, and the dynamics of stabilizer operators that underlies its topological nature. By measuring the variance of the subharmonic peak height in the Fourier spectrum, we also experimentally map out the phase boundary between the time crystalline and thermal phases. This work opens the door to harnessing this exotic phase of matter for practical quantum information processing.

Model Hamiltonian and its implementation
We consider a one-dimensional spin-1 2 chain governed by the following time-periodic 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 uni- 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 regime, it has been shown that H 2 can be in an MBL phase, where topological edge states can survive as coherent degrees of freedom at arbitrarily high energies [34]. The localization and edge states carry over to the case of periodic driving with the Hamiltonian H(t), giving rise to an FSPT phase. In this FSPT phase, the time translational symmetry only breaks at the boundary but not in the bulk. The Floquet unitary that fully characterizes the FSPT phase reads U F = U 2 U 1 , where U 1 = e −iH1 and U 2 = e −iH2 are the unitary operators generated by the Hamiltonians H 1 and H 2 , respectively. The quasienergy specturm of U F reveals that every eigenstate is two-fold degenerate and has a cousin eigenstate separated by quasienergy π (see Fig. 1c). The degenerate eigenstates also exhibit longrange mutual information between the boundary spins; this is essential for the robustness of the subharmonic response of the edge spins against local perturbations, including finite δ and V k , that respect the protecting Z 2 × Z 2 symmetry. We refer to the Methods and Supplementary Information for a more in-depth discussion about the model.
To implement H(t) with superconducting qubits, the threebody term in H 2 , which is crucial for the SPT phase at high energy, poses an apparent challenge since no three-body interaction appears naturally in the superconducting system. We thus employ the idea of digital quantum simulation [49] to implement H(t) with quantum circuits (Fig. 1d). For V k = h k = 0, we find optimal circuits in an analytical fashion that can implement H(t) with arbitrary J k and δ, while for nonvanishing V k and h k we utilize a neuroevolution algorithm [45] to design suitable quantum circuits (see Methods). With the obtained quantum circuits, we perform our experiment on a flip-chip superconducting quantum processor (Fig. 1e) with a chain of L = 14 transmon qubits denoted as Q 1 through Q L (Fig. 1a). See the Method and Supplementary Information for the details of the experimental setup.

Symmetry breaking at boundaries
The characteristic signature of an FSPT phase 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 phases. From Fig. 2a, it is evident that in the FSPT phase, the disorder-averaged magnetizations at the two ends of the chain, namely σ z 1 (t) and σ z L (t), oscillate with a 2T periodicity, for up to 20 driving cycles. In stark contrast, the local magnetizations in the bulk of the chain (σ z j (t) with 2 ≤ j ≤ L − 1) decay quickly to zero and do not show period-doubled oscillations. This unconventional behavior is independent of disorder averaging. Even for a single random disorder instance, the magnetizations exhibit 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 Information I. B for a more in-depth discussion). For comparison, we also experimentally measure the dynamics of the magnetizations in the thermal phase. 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 breaking of the discrete time translational symmetry at the boundaries can also be detected by the disorder-averaged autocorrelators defined as A j = Ψ 0 |σ z j (t)σ z j (0)|Ψ 0 . Our experimental measurements of autocorrelators for up to 40 driving cycles are plotted in Fig. 2f, again showing the breaking of time translational symmetry at the boundaries but not in the bulk. We mention that, in the FSPT phase, the local magnetizations for the edge spins exhibit a gradually decaying envelope, which could be attributed to either external circuit errors or slow internal thermalization. To distinguish these two mechanisms, we carry out an additional experiment on the "echo" circuit U echo ≡ (U † F ) t U t F , whose deviation from the identity operator measures the effect of circuit errors [18]. The square root of the output of U echo (black solid lines shown in Fig. 2f) fits well with the decaying envelope of the results obtained by evolution under U F . This indicates that the decay of the envelope is due to circuit errors rather than thermalization, which corroborates that the system is indeed in the localized phase.

Localization-protected topological states
In the above discussion, the initial states are random product states. To establish the FSPT phase, additional experiments on other initial states and other local observables are necessary. In this section, we show that the stabilizers in the bulk do not break the discrete time translational symmetry, but at the boundaries they do. To understand this, we consider the idealized cluster-state and spin-flip limit, i.e., V k = h k = 0 and δ = 0. In this limit, H 2 reduces to a summation of stabilizers: x kσ z k+1 . To break the degeneracy of H 2 , we consider adding two bound- The initial states are |0 ⊗L and data shown is averaged over 20 random disorder instances. While the bulk magnetization decays quickly to zero, the edge spins oscillate with a stable subharmonic response for 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 phase (J = ∆J = 1, V = h = ∆V = ∆ h = 0, and δ = 0.8). d, Fourier transform of experimentally measured σ z j (t) in the FSPT phase. The edge spins lock to the subharmonic frequency, which is in sharp contrast to the bulk spins. e, Fourier spectra of σ z j (t) in the thermal phase. No robust subharmonic frequency peak appears for either edge spins or bulk spins in this case. f, Time-dependence of the autocorrelator Aj = Ψ0|σ z j (t)σ z j (0)|Ψ0 for up to 40 cycles, obtained from averaging over 20 random instances deep in the FSPT phase, with the initial states prepared as random product states in the computational basis. The black solid lines show the results of "echo" circuits for two boundary qubits. ary terms J 1 S 1 (S 1 ≡σ x 1σ z 2 ) and J L S L (S L ≡σ x Lσ z L−1 ), which are commuting with all bulk stabilizers, to the Hamiltonian H 2 , forming a new cluster-state Hamiltonian H 2 = H 2 + J 1 S 1 + J L S L . 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 time-periodic Hamiltonian H(t) to 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 [50], which are widely used as a crucial diagnostic for universal topological properties of quantum phases [50][51][52][53]. 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} are the eigenvalues of ρ half ) of two experimentally prepared eigenstates for open and periodic boundary conditions, respectively. From this figure, a clear two-fold degeneracy for the low-lying Schmidt states is obtained for the open boundary conditions. This degeneracy corresponds to an effectively decoupled spin-half degree of freedom at the boundary of the bipartition. For periodic boundary conditions, the spectrum is four-fold degenerate, corresponding to two effectively decoupled spins at the two boundaries of the bipartition. The degeneracy of the entanglement spectrum and its dependence on boundary conditions marks a characteristic feature of the SPT state prepared in our experiment. We remark that the degeneracy disappears above the entanglement gap. This is due to finite size effects and experimental imperfections.
In Fig. 3c, we plot the time-dependence of local stabilizers in the FSPT phase. We observe that the stabilizers at the boundaries oscillate with a 2T periodicity, indicating again the breaking of discrete time translational symmetry at the bound- aries. In the bulk, the stabilizers oscillate with a T periodicity and are synchronized with the driving frequency, showing that no symmetry breaking occurs. This is in sharp contrast to the dynamics of bulk magnetizations, which decay rapidly to zero and exhibit no oscillation. In fact, in the FSPT phase, the system is many-body-localized and there exist a set of local integrals of motion, which are the "dressed" versions of the stabilizers with exponentially small tails [34]. The persistent oscillations of the bulk stabilizers observed in our experiment originate from these local integrals of motion and are a reflection of the fact that the system is indeed in an MBL phase.

Phase transition
We now turn to the phase transition between the FSPT phase and the trivial thermal phase. For simplicity and concreteness, we fix other parameters and vary the drive perturbation δ and the interaction strength V . Theoretically, the system is expected to exhibit an FSPT phase for small δ and V . 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 phase boundary, although not very sharp due to finite-size effects, 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., the amplitude of the 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 δ. At small δ, the system is in the FSPT phase, and the peak height remains at a value around 0.5. As we increase δ to a large value, the system transitions out of the topological phase and the peak height vanishes. This is consistent with 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 deviation as a function of δ, indicating a phase transition point around δ ≈ 0.2.

Conclusion and outlook
In summary, we have experimentally observed signatures of an intrinsically non-equilibrium Floquet SPT phase with a programmable superconducting quantum processor. In contrast to previously reported conventional time crystals, for our observed FSPT phase, the discrete time translational symmetry only breaks at the boundaries but not in the bulk. We measured the persistent oscillations of edge spins with a subharmonic frequency and experimentally demonstrated that the FSPT phase is robust to symmetry-respecting perturbations in the drive and imperfections in the experiment. In addition, we also demonstrated that the subharmonic response of boundary observables is independent of the initial state.
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 [54], with superconducting or other near-term quantum platforms. In practice, the observed FSPT phase may have applications in some quantum information processing tasks, such as quantum metrology or implementing a robust quantum memory with topological protection.

Characterization of the model Hamiltonian
To understand why time translational symmetry breaks at the boundary but not in the bulk, we consider the idealized "cluster-model" limit (V k = h k = 0) and set δ = 0. We suppose that the system is initially prepared in a random product state in the computational basis, and we use the dynamics of local magnetization as a diagnostic. In this simple scenario, the topologically 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 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 period-doubled dynamics of 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 period doubling. When V k = 0, the Hamiltonian (1) can be mapped to free Majorana fermions (see Supplementary Information I. B and, e.g., Refs. [55,56]). Further setting δ = h k = 0, we find that Eq. (1) maps onto two decoupled copies of the fixedpoint model of a Z 2 FSPT phase considered in Ref. [39]. The robustness of the subharmonic responses of the topologically protected edge spins to perturbations respecting the Z 2 × Z 2 symmetry is discussed in depth in the Supplementary Information I. B.

Logarithmic entanglement growth
For a many-body localized system, the entanglement entropy will feature a logarithmic growth [57], which is in sharp contrast to the case of Anderson localization without interactions. For the model Hamiltonian H(t) studied in this paper, we also expect a logarithmic growth of the entanglement entropy inside the FSPT phase with V k = 0. We numerically simulate the entanglement dynamics of the system deep in the FSPT phase with the time-evolving block decimation algorithm up to a system size L = 100 (see Supplementary Information II.). Our results clearly verify the logarithmic entanglement growth, which again implies that the FSPT phase in indeed many-body localized with nonvanishing V k . In our experiment, we also study the entanglement dynamics for a small system size through quantum tomography (see Supplementary Information IV.). We find that in the thermal phase, the entanglement grows much faster than that for the FSPT phase. However, due to decoherence and other experimental imperfections, we are not able to observe the logarithmic entanglement growth.

Quantum circuits for implementing H(t)
Direct implementation of the Floquet Hamiltonian H(t) with superconducting qubits faces a notable difficulty: the natural interactions hosted by the superconducting qubits are only two-body, so the three-body terms in H 2 cannot emerge directly. Fortunately, programmable superconducting qubits are universal for quantum computation; thus we can explore the idea of digital quantum simulation to emulate the dynamics of H(t). However, due to inevitable experimental imperfections, the depth of the quantum circuits is limited. As a result, obtaining well-performing circuits with an optimal depth which can implement H(t) (or equivalently the Floquet unitary U F ) is of crucial importance for the success of our experiment.
To find the desired quantum circuits, we utilize a neuroevolution method introduced in Ref. [45] which outputs a nearoptimal architecture for a family of variational quantum circuits that can implement H(t) with different random disorder instances. For a given instance of J k , V k , and h k , we use the gradient decent method to tune the variational parameters of the ansatz circuits to minimize the distance between the unitary represented by the circuit and the unitary generated by H(t) within a samll time interval. In the idealized "clustermodel" limit (V k = h k = 0), we can find a simple exact oneto-one correspondence between J k and the variational parameters, independent of the system size and the values of J k and δ. Thus, we are able to construct an analytical quantum circuit that can implement H(t) precisely, and at the same time are experimentally friendly and practical. The details of how to obtain the desired quantum circuits are given in Sec. III of the Supplementary Information.

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 FSPT phase (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 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 for a fixed interaction time (see Supplementary Information III. C). Meanwhile, each coupler is dynamically switched between two frequencies [58][59][60][61][62][63]: 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 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 Information 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 operators into the experimentally feasible single-qubit gates [X(θ), Y(θ), and Z(θ)] and two-qubit gates [CR z (±π)], where X(θ), Y(θ), and Z(θ) are rotations around x-, y-, and zaxis 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 Information III. A for the ansatz that generates 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 [64]. Simultaneous randomized benchmarkings indicate that the single-qubit gates used in this experiment have reasonably high fidelities, averaging above 0.99 (see Supplementary Information Tab. S1). Z(θ) is realized using the virtual-Z gate, which encodes the information θ in the rotation axes of all subsequent gates [65], and is combined with CZ to assemble CR z (±π). Here, we adopt the strategy reported elsewhere [62,66] to realize the CZ gate, i.e., we diabatically 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 benchmarking (see Supplementary Informations Tab. S1).
Data availability The data presented in the figures and that support the other findings of this study are available upon reasonable request from the corresponding authors.
Code availability The data analysis and numerical simulation codes are available from the corresponding authors on reasonable request.  Here we provide more details on the theory of the Floquet symmetry-protected topological phase (Sec. I), on our numerical simulations (Sec. II), and on our experimental setup (Sec. III). We also provide additional experimental data (Sec. IV).

A. Introduction to Floquet time crystals
In order to obtain a better intuitive understanding of the Floquet symmetry-protected topological (FSPT) phase, we first introduce the basic concepts behind Floquet time crystals and present a prototypical model as a concrete example.
Spontaneous symmetry breaking is an important concept in modern physics. It occurs when the steady state of a physical system does not respect the symmetries of the Hamiltonian governing this system. An important example that manifests spontaneous symmetry breaking is an ordinary crystal, which breaks the continuous spatial translational symmetry. More precisely, in a crystal, the state of the system, unlike its Hamiltonian, is not invariant under continuous translation operators. Analogously, systems that spontaneously break time-translation symmetry are named time crystals [S1, S2]. Although there is a no-go theorem for continuous time crystals at equilibrium [S3, S4], Floquet time crystals manifest themselves in many physical systems. There are two equivalent definitions of Floquet time crystals in Ref. [S5], 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-translation symmetry breaking occurs if, for every short-range-correlated state |ψ(t) at arbitrary time t, there exists an operator O satisfying ψ(t + T )| O |ψ(t + T ) = ψ(t)| O |ψ(t) , where |ψ(t + T ) = U F (T ) |ψ(t) , with U F (T ) the Floquet evolution unitary corresponding to one period T . This definition implies how to observe time crystals experimentally and is used in our paper. The second definition states that timetranslation 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 Refs. [S5, S6] as an example of a Floquet time crystal: where J k and h z k are uniformly chosen from the following intervals: J k ∈ [J/2, 3J/2] and h z k ∈ [0, h z ]. We set T = 2T 1 = 2. The Floquet evolution operator for one period can then be written as U F = exp(−iH 2 ) exp(−iπ/2 kσ x k ). We consider eigenstates of H 2 , which are product states in the computational z basis: |Θ = |{s k } with s k = ±1. Such states are easy to prepare experimentally. Since U F has the effect, up to a global phase, of flipping all spins, the state |Θ is related to another state |−Θ = |{−s k } , which is also an eigenstate of H 2 . Defining E + (Θ) and Therefore, in the subspace formed by |±Θ , U F has the matrix form Diagonalizing this matrix gives eigenvalues ± exp(−iE + (Θ)) and eigenstates |Θ ± exp [iE − (Θ)] |−Θ . The eigenstates of U F are thus paired cat states with longrange correlations. Thus, this model satisfies the second definition of a Floquet time crystal in Ref. [S5], so discrete time-translation symmetry breaking occurs in this system. (Note that, in order for these correlations to be stable to perturbations, disorder in the couplings J k and h z k that is sufficiently strong to render U F many-body localized is required.) Futhermore, as the Floquet operator has eigenvalues ± exp(−iE + (Θ)), if we diagonalize the effective Hamiltonian of the Floquet operator, we will get two eigenvalues with quasi-energy difference π. This model therefore corresponds to the π-spin-glass phase introduced in Ref. [S7].

B. Our model: the Floquet SPT phase
Unlike Floquet time crystals introduced above, the Floquet SPT phase breaks discrete time-translation symmetry only at the boundaries. To be specific, our model of the Floquet SPT phase exhibits 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. We will now present additional theoretical analysis of our model.

Localized and SPT quantum states
Our Floquet SPT phase has two distinct governing Hamiltonians during different time intervals as shown in the main text. In the first time interval, this governing Hamiltonian H 1 is the sum of one-body Pauli operators on different sites. In the second time interval, the governing Hamiltonian H 2 includes interaction among neighboring sites, which introduces the subtle many-body properties in this system. Let us begin by studying the static Hamiltonian H 2 [S8], where the parameters are chosen as in the main text. This Hamiltonian has a Z 2 × Z 2 symmetry, corresponding tô Here we are working in theσ z , and the two arrows represent the effective boundary spins and the ... denotes the bulk spins. These states are related by σ , and σ x all ≡ kσ x k . When V k , h k = 0, the two-body terms and the one-body terms make the eigenstates of this Hamiltonian depart from cluster states. However, if we keep the Hamiltonian deep in the topological phase (the phase we are interested in), i.e. J k V k , h k , we can also interpret this model from a many-body localized (MBL) perspective. Unlike the V k = h k = 0 Hamiltonian with strictly localized stabilizers as the integrals of motion, in the MBL phase, the system posses a set of mutually commuting quasi-local integrals of motion. Similarly, for open boundary conditions, there exists a quasi-local effective free spin at each edge, which contains bulk components decaying exponentially with the distance from the edge. In this case, the string-order parameter and the degeneracy of the entanglement spectra can still manifest the topological nature of the eigenstates. Moreover, while the energy spectrum is no longer exactly four-fold degenerate in a finite system, it is still nearly four-fold degenerate, and the corresponding eigenstates can still be divided into the four groups introduced above.

The emergence of the FSPT phase
Having reviewed the properties of the static Hamiltonian H 2 , let us now consider the Floquet case, wherein we periodically drive the above SPT Hamiltonian as discussed in the a. Time evolution of disorder-averaged local observables. For the edge spins, it is clear that σ z 1 and σ z 100 (which are right on top of each other, so that k = 1 is not visible) display persistent oscillations with 2T periodicity, manifesting the breaking of discrete time-translation symmetry. In stark contrast, bulk spin σ z k (plotted for k = 2, 50, and 99) decays rapidly to zero and no symmetry breaking is observed. This is the defining feature of the FSPT phase: time-translation symmetry only breaks at the boundary, not in the bulk. b. Fourier spectra of σ k . We find that σ1(ω) and σ100(ω) (which are on top of each other, so that k = 1 is not visible) have a peak at ω/ω0 = 1/2, where ω0 = 2π/T is the driving frequency. This peak is robust and rigidly locked to ω0/2, a manifestation of the robustness of the FSPT phase. For bulk spins, there is no such peak, consistent with no symmetry breaking in the bulk. c. Logarithmic entanglement entropy growth. In the FSPT phase, the system is many-body localized. We thus expect logarithmic entanglement growth, which is shown in this figure. The entanglement has an initial quick rise until time t ∼ 1/J. This initial rise corresponds to the expansion of wave packets to a size on the order of the localization length. main text: where T = 2T 1 = 2.
There exists a local effective free spin at each boundary. Since the effect of U 1 = e −iπ/2 kσ x k is to perfectly flip the spins at all sites, we obtain the following properties 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 and |D k and mixes the states |B k and |C k . Within the subspace of |A k and |D k , U F has matrix form Within the subspace of |B k and |C k , U F has matrix form Therefore, in the subspace formed by |A k , |B k , |C k , and |D k , U F has eigenvalues ± exp[−i(E A k + E D k )/2] and ± exp[−i(E B k +E C k )/2]. Thus, the Floquet effective Hamiltonian has eigen-energies (E A k + E D k )/2, (E B k + E C k )/2, (E A k + E D k )/2 + π, (E B k + E C k )/2 + π mod 2π. As the energy spectrum is four-fold degenerate ( Therefore, the original four-fold degeneracy breaks into two-fold degeneracy in the presence of the drive. This two-fold degeneracy is a remnant of the original topological order. As for the Floquet eigenstates, they are cat-like linear combinations of topological eigenstates: |A k ± |D k and |B k ± |C k . The mutual information between the two boundary spins is 2 log 2, indicating that there are long-range correlations between the boundaries. When we turn on the two-body terms and the one-body terms in H 2 , but still keep the system deep in the topological phase (J k h k , V k ), H 2 has four nearly degenerate eigenstates related by the symmetry operations. The effective free spin at each boundary becomes quasilocal. Under Floquet driving, the near-four-fold degeneracy breaks into near-two-fold degeneracy: Similarly, the eigenstates of the Floquet unitary are still cat-like states, and thus time-translation symmetry breaking can occur in this case. The stability of the FSPT phase will be discussed in more detail in Sec. I B 4.

Dynamical properties of the FSPT phase
Next, we will consider the evolution of this system and explicitly demonstrate the behavior of the FSPT phase.
As for bulk spins, assume that one bulk spinσ z j has the following expectation value in 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 state |ψ 1 can be expressed as (S19) Comparing the last two 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, bulk spins do not exhibit breaking of the timetranslation symmetry.
The above derivations tell that, for our model, the edge spins exhibit discrete time-translation symmetry breaking, while bulk spins relax very fast. Thus, the time-translation symmetry breaking only occurs at the boundaries as showing S2. The decay of boundary-spin magnetization and the phase diagram of the system. a. The decay of the first-spin magnetization, averaged over random disorder realizations. Here, the number of disorder realizations ranges from 3 × 10 4 (L = 6) to 10 3 (L = 14). The omitted parameters are chosen as in Fig. S1. We see an initial quick decay of σ z 1 , followed by a plateau that extends up to a time diverging exponentially with system size. The inset shows the exponential scaling of τ * with system size, where τ * is the time when the edge spin decays to 1/2. b. The phase diagram of the system as a function of the parameter δ in the definition of H1 and the average strength V of the two-body interactions. 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 phase. The other parameters are chosen as J = ∆J = 1, h = ∆h = 0.01.
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.
Deep in the FSPT phase, the system represented by the static many-body-localized Hamiltonian H 2 has a complete set of quasi-local integrals of motion [S9]. Therefore, spins far away from each other can build significant entanglement only after exponentially long evolution time [S10]. Thus, under the Floquet time evolution, the entanglement entropy of our system exhibits logarithmic growth, as shown in Fig. S1c, and will eventually saturate to a value proportional to the system size. Furthermore, when system size is finite, even deep in the topological phase, the Floquet time evolution will eventually lead to the decay of the spin signal at the boundaries. Indeed, the quasi-local effective free spins at the boundaries have tails that decay exponentially into the bulk. When system size is finite, these tails have an exponentially small overlap, which leads to the relaxation of the two effective free spins, with the lifetime diverging exponentially with system size. We demonstrate this phenomenon numerically in Fig. S2a.

The stability of the FSPT phase
The above considerations rely on the fact that, during each period, we perfectly flip spins at all sites. To show that the FSPT phase is indeed stable, we should make sure its defin-ing properties hold even for an imperfect drive. We follow arguments similar to those introduced in Ref. [S5].
We showed above that, in the perfect-drive case (δ = 0), the eigenstates of the Floquet evolution operator are cat-like states ψ AD ± = |A k ± |D k and ψ BC ± = |B k ± |C k . We say that an effective short-range correlated topological state satisfies ψ|σ 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 ± are all long-range correlated, with different quasienergies. Then, any experimentally prepared short-range correlated state (such as a product state) can only be formed by taking a superposition of those long-range-correlated Floquet eigenstates with different quasienergies. Thus, after one period of Floquet evolution, local observables at the edge will not be invariant, signaling a breaking of discrete time-translation symmetry. Now we add local Z 2 × Z 2 -symmetric perturbations into the system, such as an imperfect drive (δ = 0), two-body interactions (V k = 0), and single-body terms (h k = 0). As long as the system is in an MBL phase, a local perturbation will significantly affect only nearby sites. Thus, we expect that the long-range correlations in the eigenstates of the Floquet unitary will not disappear. In fact, there exists a quasi-local Z 2 × Z 2 -symmetric unitary operator U , which constructs the perturbed Floquet eigenstates from the unperturbed Floquet eigenstates. Since U is quasi-local and symmetric, it cannot destroy the long-range boundary correlations of the unperturbed Floquet eigenstates. (Note, however, that perturbations that break the protecting symmetry but maintain MBL can destroy the FSPT phase, as discussed in Ref. [S11].) Therefore, time-translation-symmetry breaking can also occur in the locally perturbed system. To explicitly show that the FSPT phase is indeed a phase, we use the string order parameter O sg as the indicator to plot in Fig. S2b the phase diagram with respect to the drive imperfection δ and the average strength V of two-body interactions.

Mapping to free fermions when V k = 0
In this section, we review the mapping of the time-periodic Hamiltonian H(t) defined in Eq. (1) in the main text [equivalently Eq. (S7)] to free fermions when the two-body interactions V k are set to zero. This is achieved by a Jordan-Wigner transformation whereby a spin operator on site k is represented in terms of two Majorana operators,α k andβ k . The Majorana operators are defined via the nonlocal mappinĝ Under this transformation, we havê The mapping thus results in redefined Hamiltonians and Note that H 2 can be rewritten as which corresponds to two decoupled Kitaev chains [S12], one on the odd and one on the even sublattice. The Z 2 × Z 2 symmetry of H(t) then manifests itself as the separate conservation of the two fermion parity operators with eigenvalues ±1. When T 1 = 1 and δ = h k = 0, the time-dependent Hamiltonian H(t) maps onto two copies of the fixed-point model for the nontrivial class D FSPT phase studied in Ref. [S13]. To see this, note that, when δ = 0 and T 1 = 1, we have (up to an unimportant overall phase factor) exp (−iT 1 H 1 ) = P even P odd .
Thus, we obtain the Floquet operator (setting T = 2T 1 = 2) If we additionally set h k = 0, this Floquet operator is just a product of two decoupled copies of the class D model considered in Ref. [S13]. The model studied in this work is thus expected to remain in this universality class for any small perturbations that respect the Z 2 × Z 2 symmetry of Eq. (S26), including finite δ, h k , and V k .

II. DETAILS OF THE TEBD METHOD
We numerically simulate the time evolution process of the FSPT phase using the time-evolving block decimation (TEBD) method. This method was proposed for the time evolution of matrix product states (MPS) [S14, S15] and is a variant of the density matrix renormalization group (DMRG) algorithm [S16, S17]. At the heart of the TEBD method lies the a b Trotter-Suzuki decomposition of the time-evolution operator U (∆t) of a short-range interacting system over 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 the current state |ψ(t) of the system to implement the time evolution.
Our FSPT phase has two distinct Hamiltonian operators in different time intervals as shown above. For the first time in-terval, 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 an MPO directly. To obtain the corresponding expectation values of 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 with ∆t t. To efficiently construct the MPO representation of U 2 (∆t), we group together terms in H 2 that commute with each other. The three-body operators are all stabilizer operators and commute with each other. For the two-body terms, they also commute with each other. For one-body terms, all of them are act on different sites and thus commute with each other. For simplicity, we denote Thus, the time-evolution operator for H 2 over time interval t is approximated by Furthermore, to make the numerics more efficient, the implementation of three-body terms and the two-body terms can be accomplished layer by layer, wherein each layer only contains operators with nonoverlapping support, so that they can be applied to the MPS in parallel. 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.

III. EXPERIMENTAL DETAILS
A. Quantum circuit ansatz Algorithm 1: Neuroevolution Method Output: Quantum circuit ansatz approximating target unitary. Input: Elementary gate set S, evolution unitary U 2 (∆t) and threshold β. G = Direct Graph(S); C = Random Generation of Quantum Circuit(G); L = Optimization(C, U 2 (∆t)); while min{L} > β do C = Quantum Circuit Extension(C, G); L = Optimization(C, U 2 (∆t)); end return argmin C {L}; 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 θ; U circuit (θ) = Unitary(C, θ); To observe the FSPT phase on a digital 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 directproduct structure of the evolution unitary U 1 (t) = e −itH1 in the first time interval, this unitary can be represented as a quantum circuit using a layer of rotation gates along the x axis. Thus, it can be constructed and implemented relatively easily. As for the second time interval, the interaction among different sites takes the time-evolution unitary far away from a direct product form, making things a little different.
With the progress of research on variational quantum circuits, we are able to adapt this method to construct the quan- tum circuit of the second-time-interval unitary. Variational quantum circuits are a powerful tool that has been intensively investigated in recent years. Algorithms based on variational quantum circuits hold great potential in the noisy intermediate-scale quantum (NISQ) era. There are many algorithms based on variational parameterized quantum circuits, such as the variational quantum eigen-solver [S18], the quantum neural network [S19], etc.. The major distinction between standard quantum circuits and variational quantum circuits is that the gates composing a variational quantum circuit are not fixed. They can be modified by tuning their parameters using different parameter-updating algorithms. As these parameters are updated, the unitary implemented by the variational circuit is also updated. We terminate the updating procedure when a satisfactory result is obtained. Our target is to find a variational quantum circuit, with some fine-tuned parameters, that approximates to high precision the evolution unitary U 2 (t) = e −itH2 in the second time interval. We accomplish this target in two steps: find an implementable variational quantum circuit ansatz that can be used to represent the target unitary and keep updating the parameters contained in this circuit ansatz to find a good approximation of the desired unitary.
We use the neuroevolution method [S20] to find a suitable variational circuit architecture. The elementary gates used in our experiments are variable-angle single-qubit rotation gates, X(θ), Y (θ), Z(θ), and a variable-angle control-rotation gate along the z axis. Each of these gates contains a variational parameter, the rotation angle. These gates can form various quantum circuit layers. i.e. quantum circuits with depth equal to one. Using the method of Ref. [S20], we construct out of these layers a directed graph, so that a quantum circuit can be represented as a path in this graph. To find the desired circuit, we follow the following procedure: 1) Randomly generate several variational quantum circuits of fixed depth based on the directed graph; 2) Update parameters contained in those quantum circuits using a 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 over the second time interval, U circuit (θ) is the unitary represented by the current quantum circuit with variational parameters θ, and d is the dimension of the corresponding Hilbert space; 3) Chose quantum circuits with small values of the loss function and extend them based on the directed graph to generate new circuits; 4) Iterate processes 2) and 3) until the loss function is below a desired threshold. The circuit ansatz giving the smallest value of the loss function 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 ansatzes used in our experiments are shown in Fig. S4. We notice that the quantum circuit for the evolution unitary over 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-timedependent parameters. So, for one driving period, the depth of the corresponding quantum circuit is fixed, i.e. for n driving periods, the depth of the corresponding quantum circuit is 6n.
With this circuit ansatz in hand, we can then use it to construct the circuits for our experiments. For a particular disorder realization of H 2 deep in the topological phase (J k V k , h k ), we begin with this ansatz containing randomly generated variational parameters θ. The the gradient of the loss function L(θ) with respect to those variational parameters is computed and is used to update the current parameters θ (n+1) = θ (n) − γ∇ θ (n) L, where γ is a given learning rate (we usually chose 0.001 ≤ γ ≤ 0.01). In our calculation, we iterate this optimization procedure until the operator fidelity [S21] satisfies Tr U 2 (∆t) † U circuit (θ) /d ≥ 0.999 (L(θ) ≤ 0.001). We then 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. On the other hand, because of the exponential growth of the dimension of the Hilbert space, the optimization for large systems is impractical. It is helpful that the quantum circuit ansatz found using the neuroevolution method can exactly represent the evolution unitary U 2 (t) = e −itH2 when H 2 has no two-body operators and no one-body operators (as shown in Fig. S4c). This indicates that we can analytically construct the corresponding quantum circuits for arbitrarily many qubits when V k = h k = 0, regardless of what values J k and δ have. In fact, in this case, we can find an exact simple one-to-one mapping between J k and the variational rotation angles in Fig. S4c. 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 correspond-ing quantum circuits are obtained using the above-described gradient-based optimization method. For 14-qubit systems, we only consider the stabilizer terms in H 2 and exactly construct the corresponding quantum circuits.

B. Device overview and measurement setup
To illustrate the idea of the FSPT phase, 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 [S22] 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 the xor y-axis and flux-bias (Z) pulses for tuning the qubit frequency and rotating the qubit state around the 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.
The processor was fabricated using the flip-chip recipe: all qubits and couplers are located on the sapphire substrate (top chip); most of the control/readout lines and readout resonators are located on the 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 [S23]. 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, then attenuated and filtered at multiple cold stages of DR, and later combined with the slow Z (DC) pulses via home- made 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 a Josephson parametric amplifier (JPA), a 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. An arbitrary microwave signal can be generated by mixing the DAC outputs with continuous microwaves using IQ mixer. DACs used to synthesize XY microwave signals and fast Z pulses in this experiment have 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 [S24]. Single-qubit gate errors are characterized by simultaneous randomized benchmarking, yielding an average gate fidelity above 0.99 (see Tab. S1).
The basic structure to implement the CZ gate consists of 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.
is a list of two frequencies for two neighboring qubits in group A (group B), chosen such that |11 and |02 in the two-qubit subspace have nearly the same energy 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 the FSPT phase. ω 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; these fidelity values are used to correct raw probabilities to eliminate readout errors as done previously [S25]. esq lists the single-qubit gate errors obtained by simultaneous randomized benchmarking. e A(B) CZ list the CZ gate errors obtained by both individual and simultaneous randomized benchmarking for qubit pairs in group A (B). We note that the qubit parameters may slowly drift over time [S26, S27].
ω 0 j /2π (GHz) 7.021 6.970 7.000 6.864 6.840 7.028 6.819 6.879 6.770 6.854 6.818 6.962 6.925 6.970 ωj/2π (GHz) 6.450 6.730 6.890 6.651 6.565 6.750 6.676 6.600 6.520 6.620 6.721 6.893 6.838 6.960 two flux tunable qubits and one flux tunable coupler, which are, respectively, denoted as Q 1 , Q 2 , and C here for 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 where a † i and a i are raising and lowering operators, and g ij is the coupling strength between each pair in {Q 1 , Q 2 , C}. The effective coupling strength between qubits is In Fig. S6a, we plot the dynamic range ofg (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 tuning the frequency of the coupler to around 10.5 GHz to turn offg. To realize the CZ gate, we apply a flux bias (fast Z) pulse to steer the coupler's frequency along the following trajectory: 10.5→7.3→10.5 GHz. Meanwhile, we turn on the fast Z pulses to bring a pair of qubits from their idle frequencies to the pair of values ω (see Tab. S1), chosen such that |11 and |02 in the two-qubit subspace have nearly the 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 leakage. 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 [S28]. 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: 1. Optimize coupler Z bias amplitude for minimum state leakage: We initialize Q 1 -Q 2 in |11 and fix their fre-quency detuning at ω 1 − ω 2 ≈ −2π × 250 MHz, following which we apply the sine-decorated square pulse with a total length of 40 ns to the coupler. We search for the optimized pulse amplitude z 0 which maximizes the |1 -state population for Q 1 , i.e., minimum state leakage. In Fig. S6b, we plot the whole landscape of state leakage as functions of the Z bias amplitudes of both the coupler and Q 1 , where the black solid line indicates how we sweep the coupler Z pulse amplitude.
2. Optimize phase factors: We fix the coupler Z pulse and sweep Q 1 's Z pulse amplitude using different initial states to calculate the three phase factors in Eq. (S32), aiming at the condition φ 3 − φ 2 − φ 1 = π. The black dashed line in Fig. S6b shows the routine of how we sweep the qubit Z pulse amplitude. We apply virtual Z gates to remove the trivial single-qubit phases.
3. Fine-tune gate parameters according to randomized benchmarking: We choose the randomized benchmarking sequence fidelity as a goal function to optimize relevant gate parameters, including the Z pulse amplitudes of both qubits and the coupler, and the single-qubit phases. We use the Nelder-Mead method to speed up the parameter optimization process.

IV. DYNAMICS OF ENTANGLEMENT
Unlike thermal phases without disorder or Anderson localized phases without interaction, where entanglement grows ballistically [S29-S31] or saturates to an area law at long times, respectively, the entanglement entropy of an MBL system grows logarithmically and saturates to a volume law in the long-time limit [S10]. In our experiment, we also extract the entanglement dynamics, through a full quantum state tomography of the reduced density matrix describing one half of 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. Using the tomographically extracted ρ half at different times, we extract the desired information about entanglement growth for the FSPT phase. Our results are plotted in Fig. S7b. From this figure, it is clear that, in the thermal phase, entanglement grows quickly and saturates to a maximal volume law (∼ L 2 ln 2). In contrast, in the FSPT phase, entanglement grows much slower. Due to decoherence and other imperfections in our experiment, we are not able to observe the logarithmic growth of entanglement. This not only demands a significant improvement of gate fidelities and a substantial increase of the coherence time, but also a more efficient and scalable approach to measure entanglement for a many-body system.