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] , in which time-translational symmetry is spontaneously broken in periodically driven systems. Pioneering experiments have observed signatures of time crystalline phases with trapped ions 9,10 , solid-state spin systems [11][12][13][14][15] , ultracold atoms 16,17 and superconducting qubits [18][19][20] . Here we report the observation of a distinct type of non-equilibrium state of matter, Floquet symmetryprotected 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 symmetryprotected 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 processors 21 .

Symmetry-protected topological (SPT) phases are characterized by nontrivial 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 in general 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 well-known scientific challenge. On the theoretical side, topological classifications of periodically driven (Floquet) systems with 4,[39][40][41][42] and without 43 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 44 , 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 owing to the presence of coherent errors in the device 45 . 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 non-equilibrium FSPT phases with a programmable array of 26 superconducting qubits (Fig. 1) with high controllability and long coherence time. We successfully implement the dynamics of prototypical time-(quasi)periodic Hamiltonians with × 2 2

Z Z , 2 Z , or no microscopic symmetries, and observe subharmonic temporal responses for the edge spins. In particular, we focus on a one-dimensional (1D) time-periodic Hamiltonian with three-body interactions and × 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 44,47,48 whereas the latter exhibit SPT order that can only be revealed through boundary effects or non-local 'string operators' in the bulk 39,41,49 . The observed boundary subharmonic 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, the 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 experimentally map out the phase boundary between the FSPT and thermal phases.

We mainly consider a 1D spin-1 2 chain governed by the following time-periodic Hamiltonian (Fig. 1b):

where δ denotes the drive perturbation; σ ˆk x z , is the Pauli matrix acting on the kth 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′ = 2, which roughly corresponds to 0.3 μs for running the corresponding quantum circuit in our experiment. We note that H(t) has a Z Z × 2 2 symmetry. For a suitable parameter regime, it has been shown that H 2 can be in an MBL phase, in which 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 = e H 1 -i 1 and U = e H 2 -i 2 are the unitary operators generated by the Hamiltonians H 1 and H 2 , respectively. The quasi-energy spectrum of U F reveals that every eigenstate is two-fold degenerate and has a cousin eigenstate separated by the quasi-energy π (Fig. 1c). The degenerate eigenstates also exhibit long-range 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 Z Z × 2 2 symmetry (Methods and Supplementary Information I).

To implement H(t) with superconducting qubits, the three-body term in H 2 , which is crucial for the SPT phase at high energy, poses an apparent challenge because no three-body interaction appears naturally in the superconducting system. We thus use the idea of digital quantum simulation 50 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 δ, whereas for non-vanishing V k and h k we use a neuroevolution algorithm 46 to design suitable quantum circuits (Methods). With the obtained quantum circuits, we perform our experiment on a flip-chip superconducting quantum processor (Fig. 1e) with a chain of L = 26 transmon qubits denoted as Q 1 to Q L (Fig. 1a). See Methods and Supplementary Information for the details of the experimental setup, and for experimental results from another processor with a chain of 14 qubits. (Supplementary Information I.B). d, A schematic illustration of the experimental circuits used to implement the time dynamics governed by 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 26 qubits used in the experiment highlighted in green. Yellow circles are functional qubits, but not used owing to limited gate fidelity. The remaining lattice sites (denoted as red circles) are non-functional qubits.

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 the disorder-averaged local magnetizations σ t ⟨ ( )⟩ j z 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 σ t

, oscillate with a 2T periodicity, for up to 40 driving cycles. In stark contrast, the local magnetizations in the bulk of the chain ( σ t ⟨ ( )⟩ j z with j L 2 ≤ ≤ -1) decay quickly to zero and do not show period-doubled oscillations. This unconventional behaviour 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 σ t ⟨ ( )⟩ j z 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 (including non-zero δ) 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,e, 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 σ t σ = ⟨Ψ | ( ) (0)|Ψ ⟩ j j z j z 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 (that is, experimental imperfections such as decoherence, pulse distortions and cross-talk effects) or slow internal thermalization (namely, an intrinsic tendency towards thermalization in the model). To distinguish these two mechanisms, we carry out an additional experiment on the echo circuit U U U ≡ ( )

, the deviation of which 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.

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, that is, V k = h k = 0 and δ = 0. In this limit, H 2 reduces to a summation of stabilizers:

+1 . We choose the initial states to be random eigenstates of H s 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 s . To manifest the topological nature of these eigenstates, we study their entanglement spectra 51 , which are widely used as a crucial diagnostic for universal topological properties of quantum phases [51][52][53][54] . To show that H(t) preserves the topological nature of the SPT states, we prepare random eigenstates of H s with both open and periodic boundary conditions, evolve the system for one driving period with H(t) and then measure the reduced density matrix ρ half of half of the system through quantum-state tomography. Figure 3b displays the entanglement spectra (eigenvalues of ρ -ln( ) half ) 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 generated in our experiment. We note 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 boundaries. 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, as shown in Fig. 2a. In fact, in the FSPT phase, the system is MBL 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.

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 x 8

Sim. (noisy) Exp.

simulations, in which the phase boundary, although not very sharp because of finite-size effects (for a small system size L = 8, the coupling between the two edge modes is not negligible and thus will decrease the central subharmonic peak height and result in a blurred boundary), 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, that is, the amplitude of the Fourier spectrum of σ t ⟨ ( )⟩ 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.30, which is consistent with the numerically predicted value of 0.34. The small deviation between the numerical prediction and experimental result is mainly attributed to finite-size effects, experimental noise and the limited number of disorder instances implemented in the experiment.

The digital simulation approach used in our experiment is generally applicable for quantum simulations of various exotic phases of matter. The model Hamiltonian in equation ( 1) possesses a Z Z × 2 2 symmetry, which can also support robust edge modes in the static equilibrium setting. For driven non-equilibrium systems, however, the edge modes may be stabilized by emergent dynamical symmetries. To demonstrate this and illustrate the general applicability of our approach, we also digitally simulate two other models with our quantum device, namely a periodically driven Ising chain with Z 2 symmetry and a quasiperiodically driven model without any microscopic symmetry (Methods and Supplementary Information VI). Our results are summarized in Extended Data Figs. 1 and2, in which robust subharmonic edge oscillations are also observed.

In summary, we have experimentally observed signatures of nonequilibrium Floquet SPT phases with a programmable superconducting quantum processor. In contrast to previously reported conventional time crystals, for our observed FSPT phases, the discrete time-translational symmetry only breaks at the boundaries and not in the bulk. We measured the persistent oscillations of edge spins with a subharmonic frequency and experimentally demonstrated that the FSPT phases are 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 digital quantum simulation approach explored in our experiment is generally applicable to the simulation of a wide range of non-equilibrium systems hosting unconventional topological phases, including those with multi-body interactions.

Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-022-04854-3.

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

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 part and evolves the random product state to a state with vanishing magnetization, resulting in no period doubling. When V k = 0, the Hamiltonian in equation ( 1) can be mapped to free Majorana fermions (Supplementary Information I.B and, for example, refs. 55,56 ). Further setting δ = h k = 0, we find that equation (1) maps onto two decoupled copies of the fixed-point model of a 2 Z FSPT phase considered in ref. 39 . The robustness of the subharmonic responses of the topologically protected edge spins to perturbations respecting the Z Z × 2 2 symmetry is discussed in depth in Supplementary Information I.B.

For a MBL 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 Article, 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 (Supplementary Information II). Our results clearly verify the logarithmic entanglement growth, which again implies that the FSPT phase is indeed MBL with non-vanishing V k . In our experiment, we also study the entanglement dynamics for a small system size (L = 6) through quantum tomography (Supplementary Information V). We find that in the thermal phase the entanglement grows much faster than that in the FSPT phase. However, because of the small system size and experimental imperfections (such as decoherence, pulse distortions and cross-talk effects), we are not able to observe the logarithmic entanglement growth (Supplementary Information VI).

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, because of inevitable experimental imperfections, the depth of the quantum circuits is limited. As a result, obtaining well-performing circuits with an optimal depth that 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 use a neuroevolution method introduced in ref. 46 , which outputs a near-optimal 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 small time interval. In the idealized 'cluster-model' limit (V k = h k = 0), we can find a simple exact one-to-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 (see Supplementary Fig. 4c for an explicit illustration of the circuit for L = 6) that can implement H(t) precisely and, at the same time, in a way that is experimentally friendly and practical. The details of how to obtain the desired quantum circuits are given in Supplementary Information III.

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-neighbour couplings (Fig. 1e), on which a chain of up to L = 26 qubits, denoted as Q 1 to Q L , that alternate with L - 1 couplers, denoted as C 1 to 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-neighbour 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 neighbouring qubit to the pairwise frequencies of group A(B) listed in ω ω ( , )

A(B) for a fixed interaction time (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 neighbouring two qubits can be initialized and operated with single-qubit gates; the other one is to turn on the nearestneighbour 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 ω j m (here, the superscript 'm' stands for 'measurement') for simultaneous quantum-state measurement. Qubit energy relaxation times measured around ω j are in the range of 7-41 μs, averaging above 30 μs. More characteristic qubit parameters, including the above mentioned frequencies, anharmonicities and readout fidelities, can be found in Supplementary Table 1. The parameters for another processor with 14 qubits used are displayed in the Supplementary Table 2.

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 the x, y and z axes 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 (Fig. 1d and Supplementary Information III.A for the ansatz that generates the gate sequences). Here X(θ) and Y(θ) are realized by applying 50-ns-long microwave pulses with a full-width half-maximum of 25 ns, for which the 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 (Supplementary Table 1). Then 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, that is, we diabatically tune the coupler frequency while keeping 11⟩ and 02⟩ (or 20⟩) for the subspace of the two neighbouring qubits in near resonance. When simultaneously running the 40-ns-long CZ gates for multiple pairs of neighbouring qubits as required in the experimental sequence, the average CZ gate fidelities can be above 0.98, as obtained by simultaneous randomized benchmarking (Supplementary Table 1).

The digital simulation strategies of our experiments are capable of simulating a wide range of models hosting unconventional non-equilibrium topological phases. To illustrate this, we also implement two other dynamical SPT phases with our superconducting quantum processor: an FSPT phase in a periodically driven random Ising chain 4 and an emergent dynamical SPT (EDSPT) phase in a quasiperiodically driven chain 67 .

The first model has a 2 Z (Ising) symmetry. For the FSPT phase (ref. 4 and Supplementary Information VI.A), the evolution is realized by applying two unitaries in an alternating fashion (Extended Data Fig. 1a) to random initial states. For the parameters chosen in our experiments, the corresponding Floquet unitary U = e e symmetry (where Z describes discrete time-translation symmetry), despite the fact that the original static Hamiltonian only possesses a Z 2 symmetry. This enlarged dynamical symmetry protects the edge modes of this phase, one at quasi-energy 0 and the other at quasi-energy π. This leads to unusual dynamics of the edge spins. If one applies this evolution to a product state in the x basis, the edge spins will return to their initial states only at even periods. In our experiments, we measure the random disorder-averaged local magnetization σ k

x during the evolution (Extended Data Fig. 1b). Persistent subharmonic oscillations are observed for the edge spins, whereas the averaged magnetizaiton in the bulk is synchronized with the driving frequency and shows no breaking of the discrete time-translational symmetry.

The EDSPT model has no microscopic symmetry (see refs. 45,67 and Supplementary Information VI.B). The evolution of an initial state is realized by applying on it a sequence of evolution unitaries at Fibonacci times, U U t F = ( = )

, with F v being the vth element of the Fibonacci sequence. Although the underlying Hamiltonian of this model includes random fields breaking all microscopic symmetries, the evolution unitary possesses a locally dressed Z Z × 2 2 symmetry emergent from the quasiperiodic drive 45,67 . The emergent symmetry hosts two non-trivial edge modes, which can be manifested by the distinct dynamics of the edge spins. In particular, the edge spins would exhibit 3v-periodic oscillations when measured at Fibonacci times t v = F v , whereas the magnetization of the bulk spins will decay to zero rapidly. In our experiment, we prepare random initial states and use the circuits shown in Extended Data Fig. 1a,b to implement the quasiperiodic driving of the system. We measure the random disorder-averaged magnetizations σ ⟨ ⟩ j z and σ ⟨ ⟩ j x at Fibonacci times. Our experimental results are summarized in Extended Data Fig. 2c, in which persistent quasiperiodic oscillations for edge spins are indeed observed. The results shown in Extended Data Figs. 1 and 2 were obtained using 12 qubits on a third device with slightly improved performance. We note that an experimental implementation of the EDSPT model with ten trapped-ion qubits has recently also been reported 45 .