Prolonging a discrete time crystal by quantum-classical feedback

Nonequilibrium phases of quantum matter featuring time crystalline eigenstate order have been realized recently on noisy intermediate-scale quantum (NISQ) devices. While ideal quantum time crystals exhibit collective subharmonic oscillations and spatiotemporal long-range order persisting for infinite times, the decoherence time of current NISQ devices sets a natural limit to the survival of these phases, restricting their observation to a shallow quantum circuit. Here we propose a time-periodic scheme that leverages quantum-classical feedback protocols in subregions of the system to enhance a time crystal signal significantly exceeding the decoherence time of the device. As a case of study, we demonstrate the survival of the many-body localized discrete time crystal phase in the one-dimensional periodically kicked Ising model, accounting for decoherence of the system with an environment. Based on classical simulation of quantum circuit realizations we find that this approach is suitable for implementation on existing quantum hardware and presents a prospective path to simulate complex quantum many-body dynamics that transcend the low depth limit of current digital quantum computers.


I. INTRODUCTION
The exploration of non-equilibrium phases of matter has become a central focus in current many-body physics research, attracting significant attention and driving numerous research endeavors [1,2].The fast development in the field has been motivated by the stunning degree of control achieved in experimental settings involving trapped ions [3,4], ultra-cold atoms [5,6], light-induced superconductors [7] and ultrafast topology switching [8], to name a few.More recently, the growing interest on making use of available noisy-intermediate scale quantum (NISQ) devices [9] has opened the path for intensive research in the field of digital quantum simulation applied to complex non-equilibrium phenomena.
Understanding the capabilities of current digital quantum simulators for realizing non-equilibrium phases of matter is a problem on demand [10][11][12][13][14][15].It is widely accepted that NISQ devices suffer from a major limitation to perform quantum computations due to their intrinsic decoherence times.In digital quantum simulation, where the study of many-body quantum dynamics is based on generic Trotterization schemes, this translates into a pronounced loss of fidelity as the number of gates composing a quantum circuit is increased [16,17].These findings naturally lead to questions about the potential of NISQ machines beyond the constraints of low circuit depth, particularly in terms of simulating quantum many-body systems and detecting emergent phases of matter that are out of equilibrium.
Recent experiments performed on quantum computers suggest that error mitigation techniques can successfully push the applicability of current NISQ devices to perform quantum simulations [18], while the long-term goal for achieving full quantum error correction [19] is still under investigation.In parallel, research has been conducted in the field of monitored quantum dynamics [14] to understand how the application of local measurement operations affects the evolution of a quantum system, particularly concerning its entanglement dynamics [20][21][22][23][24][25][26][27][28][29][30][31][32].A less explored approach in this domain is to employ the outcomes obtained from measurements to construct feedback protocols that steer the evolution in a desired direction.Such quantum-classical feedback schemes are beginning to be considered in superconducting qubit platforms as a way to improve current qubit fidelities [33].Along these lines, we pose the following question: can we use in-circuit quantum-classical feedback protocols to simulate quantum many-body dynamics beyond the intrinsic decoherence time of the device?
In this work, we propose a stabilization protocol directly implementable in quantum circuits, that combines measurements on the system with a classical error correction code.The protocol is based on the periodic application of projective measurements in randomly located domain walls across subregions of the system, with the incircuit classical computation of correlations taking place after a measurement event.More explicitly, we employ a correction scheme that identifies single qubit flips based on the construction of the correlation matrix of measurement outcomes and the comparison with its value at initialization.The outcome obtained from the classical processing is then employed to determine the target qubit undergoing a correction operation (see Fig. 1 (a) and Appendix section A 5 for details on the specific correction protocol employed).
We focus on the particular case of the MBL-DTC phase recently observed in Ref. [94], and whose model is schematically represented in Fig. 1 (b).We test and verify our scheme employing classical simulation of quantum circuits realizations, comparing the unitary protocol  (i), the open quantum system (ii) and the stabilization protocol in the open system (iii), all of which are represented in Fig. 1 (c).

II. MODEL AND PROTOCOL
As a typical example for an MBL-DTC we consider the one dimensional periodically kicked Ising model.The model consists of a one dimensional chain of L qubits with open boundary conditions and nearest-neighbour interactions, subjected to an external periodic drive of period T .The Floquet unitary operator has the form [34,60,94,95]: where σ x j , σ z j are the x, z Pauli matrices at site j of the chain, respectively.The parameter T represents the Floquet period of the external drive, which we set to T = 1.The parameter g is the pulse parameter, with g = 1 representing a perfect π-pulse on the Bloch sphere.The coupling parameters J j and the magnetic fields h j are randomly sampled from uniform distributions [94] with J j ∈ [−1.5π, −0.5π] and h j ∈ [−π, π].With this choice of the parameters, the model in Eq. ( 1) is exactly equivalent to the one in Ref. [94], and has been implemented recently on NISQ hardware [94,95].The model has three different phases (paramagnetic, thermal, MBL-DTC) depending on the values of g [34]; here we focus on the MBL-DTC phase of the model, which has been reported to exist for values of g ⪆ 0.82 [34,94].
Protocol (i) corresponds to unitary evolution of the state in a closed quantum system.For protocols (ii) and (iii), we assume that the system is open and coupled to a dissipative Markovian environment modelled as a quantum channel Φ depending on a single noise parameter p (see Appendix section A 3 for details on the noise model employed).The coherent errors that might take place on a real setup can generally be mapped to incoherent error models through the Kraus formalism, in the spirit of Ref. [103].We thus assume here an incoherent noise model, noting that the application of more sophisticated and realistic noise models are possible by employing such approaches.All time evolution protocols are carried out employing tensor networks and by representing the state of the system as a Matrix Product Density Operator (MPDO), see Appendix sections A 1, B 4 and B 5 for details.Due to the existence of an infinite number of local conservation laws [88,104], closed systems featuring MBL effects generically obey an unbounded logarithmic growth of the entanglement entropy [105,106] in the thermodynamic limit.For a one-dimensional system, such area-law behavior of entanglement translates into an efficient representation of these systems employing Matrix Product State (MPS) methods [107].

A. Stabilization for different initial states
Characterizing a true MBL-DTC phase requires independent results regardless of the choice of any initial product state [34,37].Here we study the time evolution of the model in Eq. (1) under protocols (i), (ii) and (iii), for a set of different initially prepared product states: the Néel state |10⟩ ⊗L/2 , the zeros state |0⟩ ⊗L , the wall state |1⟩ ⊗L/2 |0⟩ ⊗L/2 , and a random bit-string |random⟩.
Our first observation is that results are independent of the initial preparation of the state.In Fig. 2 (a), we have represented the disorder averaged magnetization for the four product states considered.The application of protocol (iii) leads to an overall restoring of the MBL-DTC phase in all cases, which is otherwise rapidly lost under protocol (ii) due to repeated action of the quantum channel [34] at each step n.This is our main result: The correction scheme enhances the subharmonic response beyond the intrinsic decoherence time of the device.We note that with the proposed correction scheme in protocol (iii), qubits located at the edges of any finite chain have a lower probability of correction than qubits deep in the bulk, an effect that is appreciated in Fig. 2 (a).This effect is due to the open boundary conditions employed in the correction protocol (iii), for which an edge qubit 0.6 0.7 0.8 0.9 1.0 Glassy spatial order enhancement.The spin glass order parameter χ SG as a function of the pulse parameter g, obtained by averaging over random initial product states and disorder realizations for L = 24 qubits, and p = 0.025 for protocols (ii) and (iii).We generate a total of N dis = 1280 disorder realizations for all protocols.Following Ref. [94], we employ a total of 60 time steps.χ SG is obtained by averaging the values between steps 48 and 60.The horizontal dashed line corresponds to χ SG = 1.The straight line joining data points in protocol (i) has been included as a guide to the eye.
has a probability 1/(L − M ) of being part of a domain wall.On the other hand, bulk qubits have a probability of M/(L−M ) of belonging to a randomly located domain wall, and thus of being corrected.This effect is expected to be absent for bigger system sizes, i.e. when addressing the thermodynamic limit.
In order to quantify the effect of the correction scheme onto the subharmonic response, we show in Fig. 2 (b) the spatially and disorder averaged autocorrelation ⟨⟨C(n)⟩⟩ (see Appendix section A 6). Employing protocol (iii) greatly enhances the sub-harmonic response of the system when compared to the pure dissipative case of protocol (ii).The qubit-qubit correlations between the midchain qubit and the right-most qubit of the chain have been represented in Fig. 2 (c).The results for protocol (iii) show a saturation to a finite value for increasing number of steps, a behavior in accordance with the unitary protocol (i), independent of the initial state.
To further investigate the robustness of the results with respect to different initial (product states) preparations, we proceed to calculate the Edwards-Anderson spin-glass order parameter χ SG .This extensive quantity identifies the amount of spatial spin-glass order developing in the system.As the initial condition, we prepare a set of random product states in the computational basis, keeping the same values for L, p and M in all cases.
The spin glass order parameter is represented in Fig. 3.In accordance with previous results [94], protocol (i) shows a clear transition from the thermal to the MBL-DTC phase of model Eq. ( 1), when averaged over dif-ferent initial (short-ranged) states and different disorder realizations.Under protocol (ii), glassy spatial order is absent for the considered noise model due to decoherence, showing a value χ SG ∼ 1.The application of protocol (iii) shows a restoring of the χ SG order parameter to values indicating a high degree of spatial correlation between qubits, i.e. the appearance of glassy spatial patterns as in the unitary case of protocol (i).Tuning the value of g towards the thermal phase becomes computationally expensive due to the growth of entanglement; thus, we limit ourselves to results with g ≥ 0.9.

B. Scaling and Relevance of the correction scheme
A key question to address concerns the scaling of protocol (iii) with increasing system size L and noise level p, as well as the verification of the relevant role played by the correction scheme in the observed behavior.Hence, we prepare the initial state to |10⟩ ⊗L/2 , studying variations of the L and p parameters, as well as a different number of qubits undergoing a correction operation.
In Figs. 4 (a),(b) ((c),(d)), we represent the autocorrelation (the qubit-qubit correlation) for two different noise levels.With increasing system size, we observe that for moderate noise levels p = 0.025, employing the correction protocol leads to robust sub-harmonic response and stabilization of correlations compared to protocol (ii).The enhancement of correlations in the p = 0.025 case is a consequence of the specific design of the correction scheme employed in protocol (iii), which relies on the classical processing of the correlation matrix between qubits forming the domain wall.As the correlation matrix increases its size, the identification of a defective qubit based on the total number of sign flips becomes more accurate (see Appendix section A 5 for details on the correction protocol).The finite value of correlations between qubits indicates a restoring of spatial order.For strong noise p = 0.1, we observe that protocol (iii) is incapable of overcoming noise, and the system evolves very close to the dissipative case in protocol (ii).
In Figs. 4 (e),(f), we observe a pronounced growth of the mid-chain operator entanglement entropy (see Appendix section A 6 for the definition of this quantity) for the two noise levels of p under protocol (ii).The application of protocol (iii) leads to an overall slow growth in the dynamics of entanglement, with the case p = 0.025 showing closer behavior to the unitary case.We note that the employed correction scheme does not lead to a complete saturation of entanglement in any case.In Figs. 4 (g),(h) we include statistics on the single-qubit probability of correction under protocol (iii).Increasing the value of the noise level leads to a higher probability of correction, with all curves saturating to finite values in the long-time limit.The relation between saturation values for the correction probability and system size for both noise levels is included in the inset in Fig. 4 (g).For moderate values of p, the probability eventually sta-bilizes to a finite value dependent on the system size L, with bigger L showing higher correction probabilities, as expected due to the increased possibility of detecting a bit-flip event after measuring bigger domain walls.
The results indicate that the correction scheme in protocol (iii) is able to correct to a certain degree erroneous outcomes of the qubits experiencing a measurement operation M. In the context of monitored quantum dynamics in random circuits with repeated local measurements, a field where quantum-classical feedback remains uncharted [14], we now seek to address whether the measurement operations in M alone can restore the subharmonic response without the necessity of a correction scheme.To further verify that the restoring of the DTC behavior takes place due to the applied correction scheme, we have represented in Fig. 5 the autocorrelation and the qubit-qubit correlations for protocol (iii) employing a different number of corrected qubits in C. We observe that the absence of a correction protocol (zero corrected qubits) is not enough to restore the sub-harmonic response of the system, i.e. local measurement operations M alone do not on average preserve the MBL-DTC phase.Moreover, we observe that the repeated application of measurements alone lead to a slightly faster decay than in the pure dissipative case of protocol (ii), compare Fig. 4 (a).In contrast, the composed operation C • M with at least one corrected qubit in C is able to partially restore the subharmonic response of the system.We conclude that the correction scheme employed here outperforms a pure Zeno effect in the presence of dissipation.

C. Long time limit
Here we address the question of the survival of the MBL-DTC under protocol (iii) in the large circuit depth limit, under the tuning of several parameters, namely the system size, the noise level and the pulse parameter.
In Figs. 6 (a),(b) we have represented the absolute value for both the autocorrelation and the qubit-qubit correlation function for different system sizes and domain wall sizes.We observe that under protocol (iii), there is a pronounced enhancement of the subharmonic response of the system, with a much slower decay in ⟨⟨C(n)⟩⟩ with increasing value of the domain wall M , confirming that the proposed correction protocol is able to accurately identify qubit flips that destroy the subharmonic pattern following a local measurement operation for a noise level p = 0.01.
In Figs. 6 (c),(d) we report results for the long-time behaviour under variations of the parameters p and g, fixing the system and domain wall sizes.For a noise level of p = 0.01, protocol (iii) is able to restore spatiotemporal order for times way beyond the decoherence scale in the model, with increasing values of p posing a strong limitation in the application of the correction protocol.The tuning of g towards values approaching the thermal region plays a minor role compared to the  effect of noise.
The results in Figs. 6 indicate that even in the presence of small noise, the protocol approaches the classical limit at the point g = 1, where perfect subharmonic response is expected, but no quantum entanglement develops in the system.This occurs as a consequence of the application of local projective measurements in protocol (iii), leading to an overall decrease of entanglement of the state evolution, while at the same time, increasing the purity.

IV. DISCUSSION
In this work, we have proposed a scheme to overcome dissipation effects in the MBL-DTC phase recently observed in NISQ devices [94,95].Our approach is based on the application of a periodic quantum-classical feedback protocol, where single qubit correction based on classical processing of measurement outcomes is able to enhance spatio-temporal order, eventually restoring the discrete time crystal signal well beyond the decoherence time of the device.We demonstrate the capabilities of the protocol by systematically investigating the dependence on initial state, system size and noise level, while addressing the question of how do feedback schemes based on local projective measurements affect the monitored dynamics of open quantum systems [14].
The increase in the subharmonic response with increasing system size, Fig. 6, suggests a complete stabilization of the MBL-DTC in the thermodynamic limit.This can not be proven rigorously with the system sizes we can simulate classically.Demonstrating the stabilization effect on a real NISQ device therefore presents a rare opportunity to obtain a deeper understanding about quantum many-body dynamics from a digital quantum computer.
There are several directions and extensions that can be explored using similar stabilizing schemes.From the theory perspective of DTC, we motivate further studies that employ more sophisticated noise models, differ-ent measurement schemes and alternative classical (or quantum) error correction protocols, all in the context of how spatio-temporal order gets affected by such local controlled operations.We believe that a particularly interesting question to address concerns the extension of the method to multi-qubit correction schemes employing more than a single domain wall in the measurement process.Providing further insight from a fundamental point of view in order to better understand the observed collective behaviour is also highly desirable, particularly regarding the general conditions necessary for the observation of discrete time crystals in more general open quantum systems [38,96].In this case achieving full control over the phase is essential for developing efficient and long-lived quantum memory devices.
Due to its conceptual design the approach is suitable for immediate implementation on current quantum hardware.The correction protocol discussed in this work needs sufficient idle times in the qubits in order to per-form the readout and classical processing following a projective measurement.We note that such functionality is at the heart of quantum error correction schemes, constituting a very active area of research regarding improved coherence qubit times [33].In superconducting qubit platforms [108] as well as trapped ion platforms [109], employing feedback schemes is now possible, making our protocol feasible in state-of-the-art devices.
A potential application concerns the stabilization of phases of quantum matter in digital quantum devices, where the use of current NISQ devices finds strong limitations due to decoherence processes.An interesting direction concerns the investigation of related in-circuit feedback protocols in the simulation of diverse many-body quantum systems, in particular regarding the study of recently realized measurement-induced quantum phases in trapped ion quantum computers [110].We highlight the potential application of these methods in simulating quantum systems defined in higher dimensional spaces, with qudit-based processors [111] standing out as candidates where such stabilization protocols could lead to the observation of novel out of equilibrium phases of quantum matter.

DATA AVAILABILITY
The data supporting the findings of this study are available from Zenodo at https://doi.org/10.5281/zenodo.8318707and also upon request from the corresponding author.

ACKNOWLEDGMENTS
This project was made possible by the DLR Quantum Computing Initiative and the Federal Ministry for Economic Affairs and Climate Action of Germany; qci.dlr.de/projects/ALQU.

Appendix A 1. Matrix Product Density Operators (MPDO)
We represent the state of the system ρ in terms of a finite number of variational parameters, in the form of a Matrix Product Density Operator [112,113] (MPDO): where the matrices M σ1σ 1 ′ Dj Dj+1 have dimensions D j = χ j χ ′ j and D j+1 = χ j+1 χ ′ j+1 , with r representing an auxiliary space index.To carry out the time evolution of the state, we employ the Time Evolving Block Decimation (TEBD) method [113][114][115][116][117] extended to open quantum systems, in the spirit of Ref. [118], where local updates of the tensor ρ are performed employing a series of Singular Value Decompositions (SVD).The dimension χ plays the role of the bond dimension employed in the unitary evolution of Matrix Product States (MPS) [113].The error upon SVD decompositions after application of unitary entangling gates is controlled by a TEBD tolerance parameter ε.We increase dynamically the value of the bond dimension χ whenever the tolerance ε is exceeded after every SVD decomposition.
The dimension K, which is identified with the Kraus dimension of the channel, generally grows as a result of subsequent SVD decompositions [118,119] performed in the auxiliary indices.Upon application of a non-unitary gate, only the lowest K singular values are retained.Thus, there are two main sources of error, one associated with the growth of entanglement in the system, and one associated with the statistical noise error due to the channel employed.For the data displayed in the main text, we fix ε = 10 −4 , K = 10.We have checked that choosing smaller (larger) values of ε (K) do not affect the results in a significant way (see Appendix section B 4).As a technical detail, projective measurement operations are realized as non-unitary local quantum gates, which are known to destroy the Schmidt decomposition in the standard TEBD method [114].This can be overcome by the application of identity operations sweeping across the chain, in analogy with the approach followed in the imaginary time evolution TEBD.

State initialization and time evolution protocols
Any initially prepared state is represented by a product state with factorized density matrix ρ 0 in the computational basis: To study disorder ensembles with different couplings in Eq. ( 1), we generate a total of N dis disorder realizations, each realization corresponding to a different circuit of the form represented in Fig. 1 (a) and having its own initial state ρ 0 .We study each disorder realization independently, with results being averaged over all disordered samples in the end.Note that for a single realization of the coupling parameters, any measurements carried out in real quantum hardware would need to be repeated several times for those same values of the couplings, in order to reduce shot noise.We focus on the time evolution of the state ρ at stroboscopic times n under three different time evolution protocols schematically represented in Fig. (1) (c) in the main text.
In protocol (i), the system is closed, and the unitary evolution is carried out by direct action of the Floquet unitary U F .In protocol (ii) the system is open, and we model decoherence effects in the system by applying a dissipation layer Φ immediately after the action of the unitary layer U F , with each qubit experiencing a quantum channel.The full quantum channel acts over the state as a mapping Φ(ρ) chosen to be a completely positive and trace-preserving map (CPTP).Under protocol (iii), the dissipative layer is followed by the application of a non-unitary layer formed by a measurement operation M on a subset of the qubits, which is chosen to be a domain wall of size M that is randomly located at any time step n.The second operation C consists of a (classical) correction protocol where the outcome of the measured qubits is processed in a classical way (see the correction operation section below).
The three different time evolution protocols and their corresponding operations over the state of the system at a given time step n are: Protocol (i) Unitary evolution of the state: Protocol (ii) Unitary step followed by quantum channel Φ: Protocol (iii) Unitary step, followed by quantum channel Φ, measurement M and a correction operation C: Note that the order of operations is taken from right to left, with each composed operation taking place at any step n.

Dissipative layer
We focus on arguably one of the simplest noise channel to represent the dissipative layer Φ, namely the bit-flip channel.We consider that each of the qubits experiences a bit-flip channel immediately after application of the Floquet unitary (see Fig .(1)(a) in the main text).This noise model assumes only local decoherence effects in the qubits.For the DTC considered in this work, the assumption of local bit-flip effects is a natural choice; however, extensions of the method to more realistic noise models is straightforward.
The quantum channel is described by a Completely Positive Trace Preserving (CPTP) map Φ, whose action over the state ρ is represented by a set of Q Kraus operators K i : where I is the identity matrix.The single-qubit bit-flip channel depends on a single noise parameter p, and it is defined by two Kraus operators: This quantum channel corresponds to a bit-flip occurring with probability p on each qubit.Since all qubits experience a local bit-flip channel in our case, this translates into a set of Q = 2 L Kraus operators by virtue of Eq. (A3).For single-qubit gates, current technology estimates [34] of the noise parameter are p ∼ 10 −3 .Further generalizations for the noise models to be employed are possible.We have included the case of a fully depolarizing channel model in Appendix section B 3.

Measurement operation
A set of M < L qubits forming the domain wall are measured at a given time step n under the action of the M operation.Calling x (n) 0 the coordinate of the leftmost qubit forming the domain wall at time step n, the domain wall indices at that time step is given by the set of integers: Note that the domain wall S (n) is dynamic and changes coordinates at any time step n in a random fashion.We distinguish between two different projective measurement protocols, each represented by projectors P0 and P1 in the computational basis: To simulate a measurement operation at time step n, we proceed as follows: For a selected qubit with index i ∈ S (n) , we calculate the expectation value of the zcomponent of the Pauli matrix at that particular qubit location , and generate a random number r i ∈ [0, 1] from a uniform distribution.If r i > m i , we apply the measurement protocol with P0 , simulating a collapse of the qubit over the |0⟩ state; otherwise, we employ the protocol with P1 , simulating a collapse over the |1⟩ state.For each disorder realization we then follow the trajectory of measurement outcomes that is determined by the generated random numbers.Note that these operations are non-unitary.The resulting measurement outcomes are recorded into a classical register with M values, where each value being either +1 or −1 for any index in S (n) .The measurement values are recorded into an M -entry vector: (A7) Note that the vector depends implicitly on the coordinates of the domain wall.For that particular domain wall indices set, there is an associated vector at initialization, ⃗ v(n = 0), which is known from the initial preparation of the product state.The above measurement procedure is repeated at any time step n and with the associated domain wall S (n) .

Correction operation
Here we present the correction operation C. Due to the probabilistic nature of a measurement, the value obtained for a given qubit might be erroneous (i.e. it breaks the subharmonic response pattern characterizing the DTC).To overcome this, the correction protocol C is applied immediately after the measurement operation M, with C being solely based on the initial state configuration ρ 0 and the obtained measurement outcomes at the time step n; we refer again to Fig.
(1) (a) for a schematic representation of this operation.
The initial correlation matrix between qubits at positions given by S (n) is known at time n = 0, and it is defined as the outer product of the associated vector at initialization: The associated correlation matrix of the recorded outcomes for the domain wall S (n) at time n is: Both matrices, which are the result of a classical computation, are compared by their element-wise product: where int represents the integer part, J ij is the matrix of ones of dimension M × M , and * indicates elementwise multiplication.Note that δ ij (n) is a real, symmetric matrix with integer entries.In the case where noise levels are moderate and we are deep in the DTC phase with g ∼ 1, it is expected that repeated measurements yield entries for ⃗ v(n) that follow the subharmonic pattern characterizing the MBL-DTC phase.Thus, we mask those elements that do not experience a sign change in the components δ ij (n).This way, the non-zero entries of the matrix correspond to qubits experiencing a sign flip respect to C 0 (i, j).We define the position of the qubit to be corrected at that time step n as i(n), given by: Unless otherwise stated, we correct one qubit out of the total M measured qubits forming the domain wall.If more than one qubit is corrected, the indices follow the criteria employed in Eq. (A11), from bigger to smaller values of i(n).

Observables
We identify the relevant observables that characterize the MBL-DTC phase following Refs.[34,37].We represent a spatially and disordered averaged operator by a double angle bracket.For any operator that can be decomposed as a sum of local terms O = L j=1 O j , the averaging over L lattice sites and N dis disorder realizations at stroboscopic times n is given by: where tr denotes the trace and ρ (m) n represents the state of the system for the m-th disorder realization of the circuit at time step n.
A quantity of experimental relevance [39,41,42,44,94,95,120] to describe the DTC phase is the spatially and disorder averaged autocorrelation, defined by: The disorder averaged, equal-time qubit-qubit correlations between qubits at positions L/2 and L of the chain is given by: A quantity of relevance to determine the existence of glassy spatial patterns is the Edwards-Anderson spinglass order parameter [34,94]: This parameter is averaged over the total number of disorder realizations N dis .The mid-chain entanglement entropy measures the amount of entanglement for a bipartition of the system into two equal size subregions A, B. Averaged over disorder realizations, it is given at stroboscopic times n by: For pure states, this quantity can be computed easily by representing the state as a Matrix Product State (MPS) in canonical form [113,114] (inner sums indicating contractions are implied): The entanglement entropy of the bipartition in this case is associated with the singular values encoded in Λ L/2 as: For mixed states represented by a MPDO, calculating the entanglement of the bipartition poses a challenging computational task [121,122].Since local updates of the MPDO only occur in the A σj ,r χj χj+1 in Eq. (A1), we associate, in analogy with the unitary case, all updates of singular values in the middle bond of the chain as a measure of quantum correlations in the state, i.e. as a measure to characterize entanglement dynamics in the MPDO.We therefore extend the definition of entanglement to mixed states, with the mid-chain entanglement entropy in Eq. (A18) corresponding to the limiting case under unitary evolution for pure states.We call this generalization the operator entanglement entropy.We stress that this definition does not necessarily represent the entanglement of the bipartition of the chain in the MPDO case [121,122], but rather serves as a proxy to the development of quantum correlations in the state even in the presence of non-unitary operations.Following the discussion on the initial state independence of results in the main text, here we report additional data regarding the mid-chain operator entanglement entropy, the probability of correction and the averaged qubit purities, for the four different product states considered in Sec.III A.
The average amount of entanglement for a single qubit with index j respect to the rest of a quantum system is represented by the purity of the reduced density matrix ρ n,j at stroboscopic times n, defined by the partial trace over the system degrees of freedom excluding qubit j: where Ω represents the system excluding the qubit at the lattice index j.The spatially and disorder averaged qubit purities are defined by: In Fig. 7 (a) we represent the mid-chain operator entanglement entropy.We observe nearly identical behaviour of the entanglement growth for the set of selected initial states.Due to the non-unitary operations in protocol (iii), the operator entanglement entropy grows at a much slower rate compared to protocol (ii); we argue that such slow growth could be logarithmic.The (unbounded) logarithmic growth of entanglement is a characteristic feature of quenched product states in one dimensional MBL phases in the thermodynamic limit.In Fig. 7 (b), we have represented the probability of performing a single qubit correction in protocol (iii) at different time steps.The results suggest that this probability is independent of the initially prepared state.

Long time behaviour of operator entanglement entropy and probability of correction
In Fig. 8. we have represented the mid-chain operator entanglement entropy, along with the probability of performing a correction operation in protocol (iii), for different tuning of the parameters.In Fig. 8 (a), (b), increasing system size leads to a decrease in the growth rate of the operator entanglement entropy, while at the same time the probability of correction saturates.The behavior of the operator entanglement entropy is in accordance with the behavior observed in the autocorrelator in the main text.In panels (c), (d), the major role played by the noise level p respect to the pulse parameter g is manifest.
Colors in (d) apply those in the legend in (c).

Fully depolarizing noise model
We have in addition tested the protocol under a fully depolarizing noise model, in which each qubit experiences a depolarizing channel given by Kraus operators: We note that under this definition of the quantum channel, the value of the noise parameter p is not equivalent to the one employed for the bit-flip channel.
The results are represented in Fig. 9.We observe that the main results remain robust under depolarizing noise channel for the values of noise parameter p reported, with spatial and temporal correlations showing a clear amplification for different system sizes under protocol (iii).

Numerical convergence
Here we provide evidence that variation of the tolerance threshold parameter ε and the Kraus auxiliary dimension K in the time evolution of the MPDO does not lead to different results (see Appendix A 1 for a definition of these parameters).
To test this, we focus on a single circuit realization of the system for L = 24 qubits, n s = 60 and g = 0.97.The same realization of the circuit is employed in all simulations for variations of both ε and K.For protocol (iii) we must, in addition of using the same circuit realization, employ the same measurements and corrections in all the cases; i.e. for any pair (ε, K), both the locations of the domain wall at any step n and the simulated measurement outcomes are the same.We exclude protocol (i) because this protocol does not depend on the Kraus dimension K due to unitary evolution of the state.
In Fig. 10, results for the autocorrelator and the qubitqubit correlation show nearly identical values in all cases, with the exception of K = 8, which is slightly deviated from the reference case ε = 10 −5 , K = 12.Note that results reported in the main text correspond to ε = 10 −4 and K = 10.The jumps observed under protocol (iii) correspond to the successive application of projective measurements.

Exact evolution benchmark
In order to benchmark the MPDO approach and the application of non-unitary gates as measurements, here we show comparison of the MPDO evolution of the state against the exact evolution of the density matrix operator (ED) for small chain sizes.We fix L = 8 as our system size, and create a single disorder realization for the couplings.We also create all measurement events at any time step before executing both algorithms, so that the evolution of the state is exactly equivalent under both methods.
The method is benchmarked employing the noise model advertised in the main text (see Appendix section A 3).This noise model applies a mapping over the density matrix ρ with a total of 2 L possible outcomes for the state.In the exact simulation, constructing the associated Kraus operators acting on C 2 L can be done employing the computational basis states.For a given state in the computational basis, any site is associated with a Kraus operator K 0 = √ 1 − pI (K 1 = √ pσ x ) for the zero state (one state), where I, σ x are the identity and the Pauli sigma x matrices acting on the site Hilbert space of local dimension C 2 .Thus, a single outcome ρ {s} out of the 2 L possibilities is given by: The comparison between both methods is shown in Figs.11,

FIG. 1 .
FIG.1.Quantum-classical feedback scheme for the stabilization of an MBL-DTC on a quantum computer.(a) Quantum circuit realization for a chain of L = 8 qubits.At initialization, an arbitrary product state with factorized density matrix ρ0 is prepared in the computational basis.The evolution of the state is studied under three different protocols (i), (ii), (iii), with each protocol step containing different layers.The layers are applied in a periodic fashion during n steps (circuit depth).We distinguish three different layers: a unitary layer UF , a dissipative layer in the form of a quantum channel Φ, and a non-unitary layer composed of a measurement M and a correction C operation, highlighted in grey background.For the measurement operation M, a domain wall of M qubits at an arbitrary position along the chain is measured in the z-basis; M = 3 in the example.Measurement outcomes are recorded and processed classically to compute the correlation matrix between all measured qubits Cn(i, j), which is compared with its value at initialization C0(i, j) for the same subset of qubits forming the domain wall.The qubit experiencing the most sign flips is identified and corrected.(b) A discrete time crystal (DTC) responds sub-harmonically in the presence of a external drive.The DTC can be regarded as an infinitely extended spatio-temporal lattice with perioddoubled oscillations along the time dimension, with arbitrary couplings Jj and local fields hj in space (different colors represent different coupling and field strengths).(c) The three different time evolution protocols: (i) A closed quantum system; (ii) An open quantum system; and (iii) An open quantum system subjected to a measurement and a correction operation M and C, respectively.In the panels, we show the time evolution for the z-projection magnetization of the qubits ⟨σ z j ⟩, for a single realization of the circuit with L = 12, M = 3 in (a) under protocols (i), (ii) and (iii).

FIG. 2 .
FIG.2.Time evolution simulation for different initial product states.The system size is fixed to L = 24 sites, g = 0.97, p = 0.025 (protocols (ii) and (iii)), employing a total of 80 steps on each protocol, and results are averaged over N dis = 1280 disorder realizations in all cases, with (b) and (c) only showing data for the initial state |10⟩ ⊗L/2 .The domain wall size in protocol (iii) is fixed to M = 6 qubits, where at most one of them experiences a correction (see Fig.1(a) and Appendix sections A 4 and A 5 for details).(a) The disorder averaged magnetization for the four different initial states (columns), for the three different time evolution protocols (rows), for different qubit positions and time steps n.(b) Disorder and spatially averaged autocorrelation.Legends correspond to those in (c).(c) The qubit-qubit correlations along the z-projection, between the qubit in the middle of the chain and that located at the right-most edge.

LL 1 FIG. 4 .FIG. 5 .
FIG.4.System size scaling and noise dependence.Comparison of the different time evolution protocols starting from the initial state |10⟩ ⊗L/2 , for g = 0.97 and p = 0.025, 0.1, and different system sizes L. All results were obtained averaging over N dis = 5120 disorder realizations.Straight lines between data points in all figures are included as a guide to the eye.Panels (a),(b): The spatially and disordered averaged autocorrelation at stroboscopic times.Panels (c),(d): The qubit-qubit correlations between the mid-chain qubit and the qubit located at the right-most edge of the chain.Panels (e),(f): The midchain operator entanglement entropy (see Appendix section A 6).Under protocol (iii), the measurement protocol operations lead to an overall slow growth of the operator entanglement entropy.Panels (g), (h): Probability of performing a single qubit correction under protocol (iii), for different system sizes and noise levels.The inset in (g) shows the values marked by the dashed lines as a function of 1/L, for the two noise levels considered.

1 FIG. 6 .
FIG.6.Long time limit under different tuning of parameters.Panels (a),(b): Absolute value of (a) the autocorrelation and (b) the qubit-qubit correlation function, for different L, M , g = 0.97, p = 0.01 and a total of 300 time steps.Data corresponding to protocol (iii) was averaged over N dis = 1280 disorder realizations; for protocol (ii), results are averaged over N dis = 256 for 80 time steps and L = 24.For protocol (iii), data markers have been included every ten time steps.The inset in (a) shows the subharmonic oscillations obtained for ⟨⟨C(n)⟩⟩ in the cases M = 3 (orange) and M = 6 (purple) under protocol (iii).Legend labels in (b) apply to (a).Panels (c),(d): Absolute value of (c) the autocorrelation and (d) the qubit correlation function under protocol (iii), for L = 24, M = 6, averaging over N dis = 1280 disorder realizations, for different values of g and p.Data markers have been included every ten steps.The inset in (d) shows the subharmonic response obtained in the autocorrelation in (c), for the cases p = 0.05, g = 0.99 (purple) and p = 0.01, g = 0.97 (orange).Legend labels in (c) apply to (d).

Appendix B 1 .
Operator entanglement growth and probability of correction for different initial states

FIG. 7 .
FIG.7.: Operator entanglement and probability of correction dynamics: (a) The disorder averaged mid-chain operator entanglement entropy S L/2 (n).Protocol (i) shows a very slow growth of entanglement compared to the pure dissipative case of protocol (ii).Application of non-unitary operations under protocol (iii) leads to a much slower growth rate for the mid-chain operator entanglement entropy.The inset represents the averaged qubit purities.(b) The probability of applying the correction scheme described in protocol (iii).We observe saturation to a value close to 0.5 independent of the initial state.All parameters correspond to those of Fig.(2) in the main text.

97 (M = 3 , 75 FIG. 8 .
FIG. 8. Long time dynamics of operator entanglement entropy and probability of correction: (a) Mid-chain operator entanglement entropy S L/2 and (b) the probability of performing a correction under protocols (ii) (L = 24, M = 6, N dis = 256) and (iii) (N dis = 1280 and different system sizes), for fixed p = 0.01, g = 0.97.(c) Mid-chain operator entanglement entropy S L/2 and (d) the probability of performing a correction under protocol (iii), for different values of p and g, fixing L = 24, M = 6.

FIG. 11 .
FIG. 11.Comparison of the MPDO method against the exact evolution of the state for a chain of size L = 8, with g = 0.97, p = 0.025 and M = 3 for the size of the domain wall.The figure represents the local magnetization ⟨σ z i ⟩ at different sites i of the chain.

FIG. 12 .
FIG.12.Comparison of the MPDO method against the exact evolution of the state for a chain of size L = 8, with g = 0.97, p = 0.025 and M = 3 for the size of the domain wall.The figure represents the qubit-qubit correlations ⟨σ z L/2 σ z L/2+i ⟩ at different distances i from the middle of the chain.