Realization of two-dimensional discrete time crystals with anisotropic Heisenberg coupling

Domain: physics
Published: 2026-01-28
Editorial outcome: minor_revision

Abstract
A discrete time crystal (DTC) is an out-of-equilibrium phase of matter that spontaneously breaks discrete time-translation symmetry. Previous studies have been limited to a set of models with Ising-like couplings-and mostly only in one dimension-thus precluding our understanding of the existence (or not) of DTCs in models with more realistic interactions. In this work, by combining the latest generation of IBM quantum processors with state-of-the-art tensor network methods, we demonstrate the existence of a DTC in a two-dimensional system governed by anisotropic Heisenberg interactions. We uncover a rich phase diagram encompassing spin-glass, ergodic, and time-crystalline phases, and identify the interplay of initialization, interaction anisotropy, and driving protocols in stabilizing the DTC phase. By extending the study of Floquet matter beyond simplified models, we lay the groundwork for exploring how driven systems bridge the gap between quantum coherence and emergent non-equilibrium thermodynamics.

Introduction
The mechanisms by which the laws of thermodynamics arise from the underlying unitary evolution of quantum mechanics is a topic of fundamental interest in the physical sciences. For example, it is now well understood that certain out-of-equilibrium quantum systems can indeed thermalize and hence act as their own heat bath such that local observables converge to thermal values described by a Gibbs ensemble1,2 or a generalized Gibbs ensemble3,4. These out-of-equilibrium phenomena are typically addressed by considering toy-models in the context of a “quantum quench”, whereby non-trivial dynamics are induced by a single rapid change of the Hamiltonian of the model5,6. This immediately raises the question of what might be different for the thermodynamic behavior of quantum systems outside of the quench paradigm, such as driven quantum systems7, in which energy is constantly exchanged with the environment - a situation which may more closely reflect naturally occurring thermodynamic processes. A striking example of the fundamentally different nature of quantum systems with and without a drive is provided by the emergence of discrete time crystals (DTC)8–11. Unlike equilibrium phases, DTCs exhibit macroscopic manifestations of coherent quantum dynamics, challenging the conventional narrative that thermodynamic behavior universally erases quantum signatures12.

Discrete time crystals emerge in certain periodically driven many-body-localized (MBL) spin systems13,14, where periodic driving combines with stable long-range ordering. When the drive (nearly) flips every spin, the ordering imposes a structure on the single-period unitary evolution operator, resulting in pairs of eigenstates separated by an exact π-phase rotation that remains stable against small perturbations8–10. Consequently, local observables exhibit oscillations with a period that is a multiple of the driving frequency, persisting indefinitely in perfectly isolated systems. This subharmonic response represents a spontaneous breaking of discrete time-translation symmetry, analogous to the breaking of continuous spatial symmetry in conventional solid-state crystals.

In this work, we explore a new setup for observing discrete time crystals. Previously, most studies focused on one-dimensional Ising-XXZ type couplings between spins due to their lower classical simulation complexity and straightforward implementation in experimental setups, which often rely on linear arrays15–26. Building on prior work27–30, we extend these experimental investigations to a two-dimensional system on a “decorated” (heavy) hexagonal lattice and introduce spin-flip coupling. Exploring such Heisenberg-type models is motivated by their applicability to several physical systems. Heisenberg-type interactions naturally arise in a wider range of physical systems, from single-molecule magnets and metallic chains to quantum dot-based architectures31–38. The ability to realize such models in different experimental settings underscore the broader relevance of our findings beyond the specific implementation on superconducting qubits.

Our quantum simulations demonstrate that, despite weaker localization and theoretical arguments challenging the existence of MBL in two dimensions39, disordered systems exhibit resilience to weak spin-flip coupling across various initial states, at least for short times. At the same time, such strong coupling drives a transition into an intermediate ergodic phase similar to that in one-dimensional systems40,41, see Fig. 1a. Unlike the Ising model, which is largely insensitive to the choice of initial state24,40 and governed by discrete symmetry, the addition of the Heisenberg terms to the Ising Hamiltonian, even as a perturbation, exhibits continuous U(1) symmetry. As a result, we find initial states which exhibit strikingly different behavior to the prototypical example of a Néel initial state or a randomly chosen product state. We show that for certain values of the spin-flip coupling, oscillations originating from a Néel state decay rapidly, whereas a fully polarized state displays a remarkably stable subharmonic response reminiscent of quantum many-body scars42,43, including cat-scar DTCs44,45. This asymmetry, reflected in the phase diagrams we present in this work, underscore the fundamentally different dynamical behaviors arising from Ising versus Heisenberg interactions.

Results
Similar to one-dimensional systems40, we anticipate probing transitions between three distinct phases: an ordered localized spin-glass phase, a discrete time crystal phase, and a featureless ergodic phase. To differentiate these phases, we define specific observables. A defining characteristic of the discrete time crystal phase is the appearance of persistent, period-doubled oscillations in spin dynamics, accompanied by long-range spatial order. This behavior is quantified through the evolution of the time correlation function of the same qubit.

Consider an initial product state in the computational basis, i.e., |ψ0⟩ = |s1⟩ ⊗ · · · ⊗ |sN⟩, where si ∈ {−1, 1} specifies the spin projection along the z-axis. Then the time correlation function takes the form

Δ(t) = 1/N ∑i=1^N ⟨ψ0|Zi(0)Zi(t)|ψ0⟩ = 1/N ∑i=1^N si ⟨ψt|Zi|ψt⟩, (5)

where N is the total number of qubits, the integer t denotes the number of Floquet cycles and the operator Zi(t) = UtF† Zi UtF represents the Pauli-Z operator in the Heisenberg picture. This measure effectively quantifies the individual spin memory of the initial state that persists over time.

The results of the time-correlation measurements in the time-crystalline regime (specifically at the point ϕ = 0.45π, ϵ = 0.05) are presented in Fig. 2a. Chosen to lie between the suspected transition to the ergodic regime and the Clifford (integrable) point, this point yields a robust DTC signal while keeping entanglement moderate and enabling reliable tensor-network simulations. The measured data exhibit distinct double-periodic oscillations, which gradually decay over time due to hardware noise, in contrast to the idealized classical simulations. The presence of double-periodic oscillations is further corroborated by the Fourier transform plot shown in the right panel of Fig. 2a, which reveals a pronounced peak at the frequency ω = π. To mitigate the effects of noise, we employ the renormalization technique described in Eq. (9) to recover the noiseless value of Δ(t). Specifically, parameters c(t) and c0(t) in Eq. (9) were extracted from data obtained on a 3 × 3 system and subsequently applied to the full 3 × 7 system. The recovery of the original peak is also apparent in the frequency components. Notably, we find that a single disorder realization suffices to capture the mean of single-qubit observables, owing to the large number of qubits employed. To illustrate this, we plot the error of the mean from Eq. (5) as a shaded band in Fig. 2a, b, which remains relatively small compared to the signal.

In addition to examining time-correlation functions, we also analyze spatial correlations. We define a same-time spin correlation order parameter

χ(t) = 1/M ∑⟨i,j⟩ ⟨ψt|Zi Zj|ψt⟩, (6)

and the sum is taken over all M nearest neighbors. The results are shown in Fig. 2b. Similar to the time correlations, the nearest-neighbor correlations decay over time due to hardware noise. This behavior is demonstrated by a comparison with the noiseless result obtained from classical simulations. To recover the noiseless values for the 3 × 7 system, we apply a renormalization technique analogous to that described in Eq. (9). This approach uses noise parameters extracted from a smaller 3 × 3 system (see Supplementary Information, Section III B for details). The results demonstrate that the system retains spin-spin ordering over time, which is a key feature for the emergence of the ordered regime. We also find that by expanding the set M to all site pairs (which is equivalent to determining the Edwards-Anderson spin glass order parameter24), the noisy hardware results are similar in behavior to the results of Fig. 2b, see Supplementary Information, Section IV C for more details.

In addition to low-order correlations, it is instructive to examine a more complex measure that characterizes the circuit output. One such metric is the Hamming distance, which quantifies the minimum number of bit flips required to transform the measured bit string {s0i}(t) = ⟨ψt|Zi|ψt⟩ into the initial configuration {si}. The distribution of Hamming distances Φ(d, t) can be employed to characterize various correlations, and the variance of the distribution corresponds to the quantum Fisher information (QFI). The full distribution Φ(d) is obtained by sampling the bit strings {s0i} with appropriate probabilities, a process that is carried out using quantum hardware, see Supplementary Information, Section III C for more details.

The results for the Hamming distance distribution for a 2 × 2 system are presented in Fig. 2c, d. The measured distribution is visibly influenced by hardware noise. We observe that this effect can be effectively modeled as a process in which the bit strings undergo random and independent flips with a probability p(t) that depends on the circuit cycle depth t. Using this approximation, we recover the original distributions conditioning the recovery on their mean and variance, which correspond to the classically simulated values of the time correlation and QFI (see Supplementary Information, Section III C for more details). The experiment demonstrates that the Hamming distance rapidly converges to two distinct, non-thermal distributions at odd and even periods. The shapes of these distributions are determined by the localization properties of the system. Unlike in the one-dimensional Ising model40, the distributions are nearly Gaussian, reflecting a higher degree of thermalization and exhibiting only a small degree of skewness.

The quantity defined in Eq. (5) enables us to construct order parameters that distinguish existing regimes of the system as we vary the parameters ϵ and ϕ in the model specified by Eqs. (2) and (4). Recent work30 has explored classical simulations of this problem for the special case ϵ = 0. For non zero values of ϵ, quantum correlations significantly increase the complexity of the problem (see Supplementary Information, Section III B for further details). To characterize the regimes, we define two order parameters:

ΔMBL = 1/NT ∑t=0^T ∑i=1^N |si ⟨Zi(t)⟩|, (7a)

ΔDTC = 1/NT ∑t=0^T ∑i=1^N (−1)^t si ⟨Zi(t)⟩, (7b)

where N is the total number of qubits and T ≫ 1 represents the maximum cycle depth. The parameter ΔMBL reflects the stability of many-body localization. It vanishes in the ergodic regime, whereas both the trivial spin-glass and discrete-time crystalline regimes exhibit nonzero values for any T. In contrast, the parameter ΔDTC is expected to be nonzero only in the discrete-time crystalline regime. The results are presented in Fig. 3a for T = 30. The color plots reveal the presence of a glassy regime for small values of ϕ and a discrete-time crystalline regime for ϕ near π/2, both occurring at nonzero values of ϵ. We provide a simple scaling analysis of the ϕ-dependence of the order parameter in Fig. 3b. As the linear size of the system increases from 2 (in the 2 × 2 system) to 3 heavy hexagons (in the 3 × 3 and 3 × 7 systems), we observe an increase in the order parameter for certain values of ϕ. This provides evidence of a transition between localized and ergodic regimes.

Another key parameter for characterizing potential phase transitions is the quantum Fisher information (QFI). The logarithmic growth of QFI has previously been used as evidence for the existence of a many-body localized phase34,52. The QFI values shown in Fig. 3c, provide evidence of very slow growth of entanglement within the glassy and discrete-time crystal regimes. In those regimes, the limiting value of QFI above the shot-noise limit (FQ = 1) also indicates that the correlations may not have completely thermalized; see Section IV A of the Supplementary Information for more detail. In contrast, for points located in the center of the parameter space, we observe a faster rate of QFI growth, consistent with the behavior expected in an ergodic regime.

As illustrated in Fig. 1a, we expect that the polarized initial state will exhibit behavior similar to that of the Néel state. Trends in the correlator order parameter of Eq. (5) derived from experimental data, for several other random initial states, are presented in the Supplementary Information, Section IV B. A notable exceptions to this expectation are fully polarized states or states close to them by Hamming distance. Polarized states are distinct because they are eigenstates of the spin-flip component of the Hamiltonian in Eq. (3). Consequently, these states exhibit long-time prethermal behavior, even for ϵ ∼ 1. Figure 4 shows the renormalized hardware data for the polarized initial state, emphasizing the pronounced difference in the amplitude of discrete time crystal oscillations between polarized and Néel states. The observed stability of these oscillations closely resembles the behavior of quantum many-body scars42,43 and cat-scar DTCs44,45. Differences in the rate of QFI are also observed for the polarized initial state when ϵ is perturbed away from the Ising limit (ϵ = 0), indicating that the symmetries of the underlying Hamiltonian play a significant role in the correlations dynamics.

The stability of the fully polarized state in the XXZ Hamiltonian under periodic transverse kicks raises questions about weaker ergodicity breaking. The introduction of periodic kick pulses with generic phases (i.e., excluding nπ/2 for integer n) break the U(1) symmetry of the XXZ Hamiltonian, potentially leading to thermalization through eigenstate mixing during Floquet evolution. Remarkably, the fully polarized state shows robustness near ϕ = 0 and π/2, resulting in a phase diagram distinct from earlier findings40, and this stability may stem from stable state evolution trajectories that prevent eigenstate mixing43. Such trajectories, resembling those in cat scar DTCs, suggest an unconventional mechanism preserving magnetic order even in the presence of perturbations44,45.

Methods
Signal recovery from noisy data
The presence of hardware noise requires robust techniques for extracting signals, such as the expectation values of Pauli operators, from noisy measurements. While conventional error mitigation methods offer reasonable improvements for short-time dynamics, they fail to recover signals in circuits deepen enough to explore time-crystalline behavior within our framework. To address this limitation, we introduce a physics-inspired recovery procedure that leverages specific characteristics of the noisy dynamics. In particular, we show that the influence of noise on a collective observable O can be effectively modeled as

⟨O⟩noisy = f(t)⟨O⟩ + c(t), (8)

where ⟨⋅⟩noisy and ⟨⋅⟩ denote noisy and noiseless expectation values, respectively. Here, f(t) represents an amplitude attenuation factor that decays asymptotically, i.e., f(t) → 0 as t → ∞, while c(t) denotes an offset capturing the observable’s long-time behavior. We further assume that the offset exhibits a double-periodic structure, such that c(t) = c(t + 2). This effect can be illustrated by Fig. S3 in the Supplementary Information where the parameters of the amplitude attenuation function and the offset are estimated from the actual data.

Assuming that f(t) exhibits a weak dependence on the model parameters, the observable ⟨O⟩ can be recovered using the following expression:

⟨O⟩ = ⟨O⟩noisy − c(t) / ⟨O⟩0noisy − c0(t) (9)

Here, ⟨O⟩0noisy and ⟨O⟩0 represent noisy and noiseless expectation values of a “reference” observable, measured on the device for the spin-flip interaction strength ϵ and X-gate rotation angle ϕ corresponding to a Clifford point (ϵ = 0 and ϕ = 0 or ϕ = π/2). The four free parameters, i.e., odd-period and even-period values for each offset parameter c(t) and c0(t), are determined via a convex optimization procedure. This procedure fits the proposed ansatz to a classical simulation of the observable on either 2 × 2 or 3 × 3 heavy hexagons. These parameters are subsequently used to reconstruct ⟨O⟩ on 3 × 7 heavy hexagons. Intuitively, this procedure cancels the multiplicative attenuation factor f(t) by normalizing against Clifford point data, while the additive offsets are learned from classical simulations. For further details, as well as extensions to more sophisticated models required to recover the mean square values of two-point correlators and Hamming distances, we refer the reader to Supplementary Information, Sections III B and III C.

Classical simulation
We employ Matrix Product State (MPS) methods by mapping the decorated hexagonal topology onto a one-dimensional chain with long-range interactions. In this approach, the unitary operators in Eq. (3) are represented as Matrix Product Operators (MPOs). We then utilize standard algorithms to contract an MPO with an MPS48,56, enabling us to obtain a representation of the time-evolved wave function. While the bond dimension of the MPS representing the state grows rapidly with each application of UF, the bond dimensions required to express the MPOs for each U(k)XXZ remain relatively modest. These bond dimensions primarily depend on how the two-dimensional topology is mapped onto the one-dimensional chain. A key limitation of the MPS approach is that the bond dimension of the MPO can grow exponentially with the number of overlapping bonds. We outline strategies to minimize this effect in Supplementary Information, Section II A; see also ref. 57. This method is implemented using the TeNPy software package58,59.

We also utilize classical methods that leverage a tensor network architecture matching the topology of our model60. Specifically, we employ an algorithm based on belief propagation51 to determine a gauge transformation that maps the two-dimensional network into an approximate canonical form, referred to as the “Vidal-gauge" in ref. 51, inspired by related work in one dimension61. Single-site observables are computed by applying a rank-one approximation to the environment of the site; see ref. 51 and Section II B of Supplementary Information for further details. Higher-weight observables can also be evaluated using this approach by appending a Clifford circuit to the state and measuring the corresponding single-site operator. This procedure reproduces the expectation value of a two-site operator in the original state, albeit at the cost of increasing the entanglement within the system. For example, we employ this method to compute two-point correlation functions below. These calculations are implemented using the software package ITensorNetworks.jl62, which is built on top of the package ITensors.jl63.

Implementation details
The circuits were designed and parameterized by c(ϕ, ϵ, s, T) for a particular ϕ, ϵ point on the phase diagram, where s is the index of the seed used to calculate the disorder in J, and T ∈ [0, Tmax] is the number of Floquet repetitions. Thus, a total of Tmax + 1 circuits are required to characterize the dynamics of each point on the phase diagram for a given disorder set. The circuits were then transpiled with the preset pass manager optimization level set to 1. This level of transpilation includes single-qubit gate optimization and inverse cancellation. The layout transpilation step was disabled to ensure a manual mapping of physical to virtual qubits. To limit the memory required for the auxiliary classical computing steps for each job after circuit submission, each set of circuits was divided into several circuit splits, for a total of Ncs circuits splits per run. A predetermined fractional number of Ncs, along with all circuits for the two Clifford points ϵ = 0.0, ϕ ∈ {0, 0.5π}, were submitted using Session in Qiskit for a total of NS sessions. All sessions used the SamplerV2 Qiskit primitive. For all sessions, device-level dynamical decoupling, Pauli gate twirling, and Pauli measurement twirling were disabled. Details of each run are given in Supplementary Table II of Supplementary Information.

Determination of statistical uncertainty in experimental data
There are two sources of error relevant to assessing the significance of our results. The first is the statistical measurement error, including shot noise arising from the finite number of shots. The second stems from the variability of random-disorder realizations of the exchange-coupling strength Jij in Eq. (4). For all experimental data, we employ a standard bootstrap procedure applied to the data to estimate the shot noise error. Our procedure generates 200 sampled bitstring sets containing the same number of shots as the original experiment and reports this statistical error in terms of standard deviations (σ) of the resulting spread of the plotted quantity. The error is also propagated through and amplified by the renormalization procedure. To quantify the finite-size error, we report the error in the mean value for the order parameters shown in Eqs. (5) and (6), respectively, as

ΣΔ(t) = √1/N(N−1) ∑i=1^N (si ⟨ψt|Zi|ψt⟩ − Δ(t))^2

Σχ(t) = √1/M(M−1) ∑⟨i,j⟩ (⟨ψt|ZiZj|ψt⟩^2 − χ(t))^2 (10)

These quantities are expected to vanish in the limit N, M → ∞, indicating self-averaging of the signal due to the large number of qubits. Values of ΣΔ(t) are shown in Figs. 2a, 4b, and Σχ(t) are shown in Fig. 2b as shaded regions between Δ(t) ± ΣΔ(t) for the polarization order parameter and between χ(t) ± Σχ(t) for the correlation order parameter. For visual clarity, the even and odd Floquet cycles are treated separately.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors use the IBM Heron r2 quantum processor to simulate discrete time crystals (DTCs) in two dimensions with a\nFloquet operator governed by anisotropic Heisenberg couplings. The authors observe the first signatures (the subharmonic\nresponse) of a two-dimensional DTC in the disordered case. This feature is relevant, as disorder has been hypothesized to\nstabilize the phase via many-body localization (MBL). The authors also map out the corresponding phase diagram as a\nfunction of spin-flip interaction and X-gate angle. These findings are made possible by extensive error mitigation and\nsupplemented by one- and two-dimensional tensor network simulations.\nThis work lies at the intersection of several relevant research fields (quantum computing, dynamical phases of matter, two-\ndimensional many-body localization, quantum many-body scars) and demonstrates the stability of two-dimensional DTCs for\nthe first time. In passing, the authors also note that their work shows the stability of two-dimensional MBL on the time scale of\ntheir quantum simulations. They also observe prethermal quantum many-body scars. Overall, I consider the manuscript a\nhighly relevant piece of work.\nThe manuscript is well-written, accessible, and its main claims are well supported. I only have three minor comments:\n1. The dots in Fig. 1b are extremely small and barely visible. I’d suggest removing them, as they don’t improve the\npresentation.\n2. In the caption of Fig. 2 and in line 128 the authors are trying to refer to a section in the Supplemental Material which is,\nhowever, missing.\n3. The sentence “This behavior is depicted in Fig. 3a.” in line 225 seems out of place, as the previous sentence refers to Fig.\n3c.", "reviewer_2": "The manuscript presents an impressive experimental realization of a 2D discrete time crystal (DTC) using up to 144 qubits\non IBM quantum processors. The authors investigate a 2D system governed by anisotropic Heisenberg interactions,\nextending previous studies that primarily focused on 1D Ising-like couplings. A key strength of this work lies in combining\nstate-of-the-art quantum hardware with sophisticated tensor network methods and noise mitigation strategies. The\ncomprehensive analysis reveals a rich phase diagram encompassing time-crystalline, MBL spin-glass, and ergodic phases,\ndemonstrating the tunability of these phases through various control parameters.\nIt is quite impressive that the authors implement the many-body dynamics with 144 qubits up to 450 layers of two-qubit\ngates. We would recommend publication of this work in Nature Communications, provided that the authors satisfactorily\naddress the following questions/comments.\nMajor questions:\n1. All plots presented in this work, particularly those showing experimental data, do not contain error bars. The authors state\nthat \"using a single disorder realization is sufficient for describing sums of single-qubit observables, given the large number\nof qubits employed.\" We would like to see the authors verify this argument by adding error bars calculated from different\nqubits, or from multiple experimental runs if feasible. This is crucial for assessing the statistical significance of the observed\nphenomena and the reliability of the phase boundaries.\n2. The finite-size scaling in Figure 3b is inspiring. However, its generality is limited by the lack of error bars and the\nconsideration of only two different system sizes (2×2 and 3×7). We strongly suggest the authors add data for another\nintermediate system size, for example, the 3×3 system, to provide more robust evidence for the finite-size scaling.\nWe also have some minor points that we would like the authors to respond to:\n1. As the authors correctly noted, the existence of MBL in 2D systems is a controversial topic and still under debate.\nTherefore, we suggest modifying the sentences, for example, in line 82 \"The disorder in the coupling induces a many-body\nlocalization regime,\" to a more cautious phrasing.\n2. In Figure S14, the authors show the Edwards-Anderson spin glass order parameter for different system sizes. It appears\nthat the order parameter decreases as the system size increases. We would expect that this order parameter would\ndecrease to zero in the thermodynamic limit, and hence not be in the ordering phase. Is this true?\n3. In Figure 3c, the QFI of the DTC phase increases to a larger value than that of the ergodic case. Could the authors provide\nmore intuition or a physical explanation for this?\n4. In a previous paper Nat. Commun. 15, 8963 (2024), the authors have implemented the 2D time-crystalline order with the\nplanar toric code on a superconducting processor. A proper citation to this paper and some other closely related ones should\nbe added.\n5. There are some section numbers missing in the manuscript. For instance, in the Figure 2 caption, it reads, \"Noise\nrecovery is performed using the renormalization methods described in Sections ? and III of Supplementary Information.\"\nSimilarly, in line 128, it states, \"...be found in Supplementary Information, Sections I and ?.\"", "reviewer_3": "The manuscript reports the experimental realization of a two-dimensional discrete time crystal (DTC) with anisotropic\nHeisenberg coupling in a decorated hexagonal lattice system. By using 144-qubit processor experiments and two-\ndimensional tensor network simulation techniques, the authors successfully reveal a rich phase diagram encompassing\nspin-glass, ergodic, and time-crystalline phases in a two-dimensional system governed by anisotropic Heisenberg\ninteractions with periodic driving. Experimental results demonstrate that the synergistic interplay between disorder and\nperiodic driving can maintain stable quantum dynamical behavior in two-dimensional systems.\nThey introduce the spin-flip coupling to construct the anisotropic Heisenberg model, which can be obtain more interesting\nphenomenon compared to Ising model. They observed the fully polarized state demonstrates a stable subharmonic\nresponse that resembles the behavior of quantum many body scars including cat-scar DTCs. These achievements not only\nexpand our understanding of nonequilibrium quantum phase transitions but also lay the foundation for investigating how\ndriven systems establish connections between quantum coherence and emergent nonequilibrium thermodynamics.\nHere I have a few questions/comments below for the authors to address.\n1. Under the weak coupling condition of ε=0.05, how can we demonstrate that the system exhibits characteristics of\nthe anisotropic Heisenberg model rather than behavior of a perturbed Ising model?\n2. The experimental results in the manuscript are primarily based on the Néel state as the initial state. For other initial states,\ncan similar discrete time crystal phase behaviors be observed? How do different initial states influence phase diagram?\n3. The manuscript repeatedly mentions the impact of hardware noise on experimental results, but lacks analysis of its\nspecific sources and quantitative effects on the outcomes. This needs to be supplemented and improved.\n4. The manuscript selects φ=0.45π as a representative DTC parameter (Fig 2a). Why choice this angle? Was this angle\nsystematically optimized to maximize subharmonic response robustness?"}, "author_response": "Response to the First Reviewer\nReviewer: “The authors use the IBM Heron r2 quantum processor to simulate discrete time crystals (DTCs) in\ntwo dimensions with a Floquet operator governed by anisotropic Heisenberg couplings. The authors observe the first\nsignatures (the subharmonic response) of a two-dimensional DTC in the disordered case. This feature is relevant, as\ndisorder has been hypothesized to stabilize the phase via many-body localization (MBL). The authors also map out the\ncorresponding phase diagram as a function of spin-flip interaction and X-gate angle. These findings are made possible\nby extensive error mitigation and supplemented by one- and two-dimensional tensor network simulations.\nThis work lies at the intersection of several relevant research fields (quantum computing, dynamical phases of matter,\ntwo-dimensional many-body localization, quantum many-body scars) and demonstrates the stability of two-dimensional\nDTCs for the first time. In passing, the authors also note that their work shows the stability of two-dimensional MBL\non the time scale of their quantum simulations. They also observe prethermal quantum many-body scars. Overall, I\nconsider the manuscript a highly relevant piece of work.\nThe manuscript is well-written, accessible, and its main claims are well supported.”\nAuthors: We thank the Reviewer for the positive and constructive assessment. Below we address their comments.\nReviewer: “I only have three minor comments:\n1. The dots in Fig. 1b are extremely small and barely visible. I’d suggest removing them, as they don’t improve the\npresentation.”\nAuthors: FollowingtheReviewer’ssuggestion, wehaveremovedthedotsenumeratinglayersinFig.1bandadjusted\nthe caption accordingly for legibility. The caption references the gate order as:\nEdge colors encode the application order of two-qubit gates, as shown in the rightmost panel: red (1st), blue\n(2nd), green (3rd). Single-qubit gates are applied beforehand and their depiction is omitted for clarity.\nReviewer: “2. InthecaptionofFig. 2andinline128theauthorsaretryingtorefertoasectionintheSupplemental\nMaterial which is, however, missing.”\nAuthors: We thank the Reviewer for catching this typo. In the revised version, we have corrected the reference in\nthe caption of Fig. 2 as:\nNoise recovery is performed using the renormalization methods described in Methods and Supplementary Infor-\nmation, Section III.\nReviewer: “ 3. The sentence “This behavior is depicted in Fig. 3a.” in line 225 seems out of place, as the previous\nsentence refers to Fig. 3c.”\nAuthors: We agree with the Reviewer and have removed the sentence in question because it did not properly refer\nto the figure.\nReviewer: “The manuscript presents an impressive experimental realization of a 2D discrete time crystal (DTC)\nusing up to 144 qubits on IBM quantum processors. The authors investigate a 2D system governed by anisotropic\nHeisenberg interactions, extending previous studies that primarily focused on 1D Ising-like couplings. A key strength\nof this work lies in combining state-of-the-art quantum hardware with sophisticated tensor network methods and noise\nmitigation strategies. The comprehensive analysis reveals a rich phase diagram encompassing time-crystalline, MBL\nspin-glass, and ergodic phases, demonstrating the tunability of these phases through various control parameters.\nIt is quite impressive that the authors implement the many-body dynamics with 144 qubits up to 450 layers of two-\nqubit gates. We would recommend publication of this work in Nature Communications, provided that the authors\nsatisfactorily address the following questions/comments.”\nAuthors: We are grateful to the Reviewer for their positive assessment and recommendation to publish our work.\nTheir questions and comments are addressed below.\nReviewer: “Major questions:\n1. All plots presented in this work, particularly those showing experimental data, do not contain error bars. The\nauthors state that ”using a single disorder realization is sufficient for describing sums of single-qubit observables,\ngiven the large number of qubits employed.” We would like to see the authors verify this argument by adding error\nbars calculated from different qubits, or from multiple experimental runs if feasible. This is crucial for assessing the\nstatistical significance of the observed phenomena and the reliability of the phase boundaries.”\nAuthors: We thank the Reviewer for raising this point. Two independent sources of uncertainty are relevant: (i)\nfinite-shot measurement noise and (ii) finite-size disorder variability, as pointed out by the Reviewer. Accordingly,\nwe now include two types of error bars in Figs. 2a,b, 3b,c, and 4b. Statistical error is represented by standard error\nbars for all plots estimated via a nonparametric bootstrap (200 resamples) across shots for each data point. Except\nin Fig. 3c, these error bars show 5σ standard errors. In Fig. 3c (quantum Fisher information), we report 2σ errors\ndue higher error values. To illustrate the finite-size error of the order parameters in Figs. 2a,b and 4b, we also plot\na shaded band indicating the standard error of the mean across individual polarizations/correlators (see our changes\nto the Methods section quoted below). This band reflects finite-size spread in the mean of the order parameters due\nto the finite number of qubits and therefore quantifies the degree of disorder self-averaging.\nWe have modified the text accordingly as follows.\nCaption of Fig. 2:\nWhere present, error bars indicate statistical errors corresponding to five standard deviations (5σ), estimated\nvia a nonparametric bootstrap procedure with 200 resamples. The shaded regions represent the standard error\nof the mean used as the measure of finite-size error (see Methods).\nCaption of Fig. 3:\nErrorbarsandshadedregionsinbandcencompassthestatisticalerrorestimatedbyanonparametricbootstrap\nprocedure (see Methods). Panel b includes a 5σ error, panel c includes a 2σ error.\nCaption of Fig. 4:\nError bars in b encompass the statistical error of 5σ as estimated by a nonparametric bootstrap procedure;\nshaded region shows the error of the mean used as the measure of finite-size error (see Methods).\nWe have moved the modified statement about disorder self-averaging from line 116 to line 175 in the colored\nmanuscript:\nNotably, we find that a single disorder realization suffices to capture the mean of single-qubit observables, owing\ntothelargenumberofqubitsemployed. Toillustratethis,weplottheerrorofthemeanfromEq.(5)asashaded\nband in Fig. 2a,b, which remains relatively small compared to the signal.\nWe have added the following Methods section (see line 524 in the colored manuscript):\nThere are two sources of error relevant to assessing the significance of our results. The first is the statistical\nmeasurement error, including shot noise arising from the finite number of shots. The second stems from the\nvariability of random-disorder realizations of the exchange-coupling strength J in Eq. (4). For all experimental\nij\ndata, we employ a standard bootstrap procedure applied to the data to estimate the shot noise error. Our\nprocedure generates 200 sampled bitstring sets containing the same number of shots as the original experiment\nandreportsthisstatisticalerrorintermsofstandarddeviations(σ)oftheresultingspreadoftheplottedquantity.\nTheerrorisalsopropagatedthroughandamplifiedbytherenormalizationprocedure. Toquantifythefinite-size\nerror, we report the error in the mean value for the order parameters shown in Eqs. (5) and (6), respectively, as\n(cid:118)\n(cid:117) N\n(cid:117) 1 (cid:88)(cid:16) (cid:17)2\nΣ (t)=(cid:116) s ⟨ψ |Z |ψ ⟩−∆(t) ,\n∆ N(N −1) i t i t\ni=1\n(10)\n(cid:118)\n(cid:117) 1 (cid:88)(cid:16) (cid:17)2\nΣ (t)=(cid:117) ⟨ψ |Z Z |ψ ⟩2−χ(t) .\nχ (cid:116)M(M −1) t i j t\n⟨i,j⟩\nThesequantitiesareexpectedtovanishinthelimitN,M →∞,indicatingself-averagingofthesignalduetothe\nlarge number of qubits. Values of Σ (t) are shown in Figs. 2a and 4b, and Σ (t) are shown in Fig. 2b as shaded\n∆ χ\nregions between ∆(t)±Σ (t) for the polarization order parameter and between χ(t)±Σ (t) for the correlation\n∆ χ\norder parameter. For visual clarity, the even and odd Floquet cycles are treated separately.\nReviewer: “2. The finite-size scaling in Figure 3b is inspiring. However, its generality is limited by the lack of\nerror bars and the consideration of only two different system sizes (2×2 and 3×7). We strongly suggest the authors\nadd data for another intermediate system size, for example, the 3×3 system, to provide more robust evidence for the\nfinite-size scaling.”\nAuthors: We agree that the 3×3 system size data is helpful. In the revision of our manuscript, we have added the\n3×3 results to Fig. 3b. We would like to note that these data are not very different from the 3×7 data. This is\nlikelybecausetherelevantscalingparameteristhelinear size ofthesystemwhichisoftendeterminedbytheshortest\ndistancebetweenoppositeboundaries. Sinceboththe3×7and3×3systemssharethesamelinearsize,theresulting\ncurvesaresimilar. Unfortunately,thedevicegeometrydoesnotallowustomeasureasubsetoflinearsizelargerthan\nthree decorated hexagons.\nWe have added the following comment to line 225 in the colored manuscript:\nWe provide a simple scaling analysis of the ϕ-dependence of the order parameter in Fig. 3b. As the linear size\nof the system increases from 2 (in the 2×2 system) to 3 heavy hexagons (in the 3×3 and 3×7 systems), we\nobserveanincreaseintheorderparameterforcertainvaluesofϕ. Thisprovidesevidenceofatransitionbetween\nlocalized and ergodic regimes.\nReviewer: “We also have some minor points that we would like the authors to respond to:\n1. As the authors correctly noted, the existence of MBL in 2D systems is a controversial topic and still under debate.\nTherefore, we suggest modifying the sentences, for example, in line 82 ”The disorder in the coupling induces a many-\nbody localization regime,” to a more cautious phrasing.”\nAuthors: The Reviewer correctly notes that the existence and stability of a long-lived many-body-localized (MBL)\nphase in two dimensions remain under active debate. Consistent with the literature (e.g., [1, 2]), we deliberately use\nthe term “MBL regime” as a cautious alternative to “MBL phase”. This term indicates that the observed dynamics\nareconsistentwithshort-timeMBLphenomenology,whileasymptoticstabilityandabonafidethermodynamicphase\ntransition cannot be asserted. We agree that this clarification should be made explicit, and we have revised the text\naccordingly.\nIn the current version, we have added a footnote (see line 114 in the colored manuscript):\nThe existence of a thermodynamic localized phase in two dimensions is still debated. We follow the literature\n(e.g., [1, 2]) and use the term “regime” instead of “phase,” indicating short-time MBL phenomenology with a\ncrossover to ergodic behavior rather than a confirmed thermodynamic transition.\nThroughout the text, we have also changed the term “phase” to “regime” when describing many-body localization.\nsizes. It appears that the order parameter decreases as the system size increases. We would expect that this order\nparameter would decrease to zero in the thermodynamic limit, and hence not be in the ordering phase. Is this true?”\nAuthors: The Reviewer’s observation is correct regarding the figure displaying raw (noisy) data. The limiting value\nofthedisplayednoisyEdwards-Andersonparameterreacheszero,unlikeinnoiselesssystems. Inparticular,FigureS14\nshows the parameter computed for all spin pairs,\n1 (cid:88)\nχ = ⟨Z Z ⟩2,\nSG N2 i j\ni,j\ncalculated with raw experimental data, reflecting noisy dynamics. In this case, the Edwards-Anderson parameter\ndecays due to vanishing correlators ⟨Z Z ⟩| until it reaches a trivial 1/N value coming from the diagonal terms\ni j i̸=j\n⟨Z Z ⟩ ≡ ⟨I⟩ = 1. Thus, even for large N, a zero value of the Edwards-Anderson parameter from raw experimental\ni i\ndata is not sufficient to draw a conclusion about the actual phase of the system.\nTo assist readers in making the connection with the raw experimental data, we have placed the Edwards-Anderson\nspin glass order parameter section, including Fig. S14 and the discussion, as a subsection of a new section of the\nSupplementary Information, Trends in order parameters derived from raw data. We have added the above equation\ntothissubsectionandincludednewpointsinthediscussionofthefigure(seeline928inthecoloredmanuscript):\nAs described in the main text, if the set of spin pairs M in Eq. 6 is expanded to all qubit pairs, one recovers an\nEdwards-Anderson spin glass order parameter\nχ (t)= 1 (cid:88) ⟨Z Z ⟩2. (S32)\nSG N(N −1) i j t\ni,j\nThe raw values of χ are shown in Fig. S16, and are similar to the raw data shown in Fig. 2b. The parameter\nSG\ndecays due to vanishing correlators ⟨Z Z ⟩| until it reaches a trivial 1/N value coming from the diagonal\ni j i̸=j\nterms ⟨Z Z ⟩ ≡ ⟨I⟩ = 1 in Eq. (S32). We note that because of the noisy data, a zero value of this parameter\ni i\nfrom raw experimental data (even for the large N investigated in this work) is not a sufficient condition to see\nthe actual phase of the system.\nReviewer: “3. In Figure 3c, the QFI of the DTC phase increases to a larger value than that of the ergodic case.\nCould the authors provide more intuition or a physical explanation for this?”\nAuthors: WethanktheReviewerforthisniceobservationandinterestingquestion! Toaddressthis, andinresponse\nto the Reviewer, we provide the following intuition and explanation which has been added as a new section in the\nSupplementary Information (see line 876 in the colored manuscript):\nIn our analysis of the Hamming distributions of Section IIIC, we define the variance of the Hamming distance\ndistribution in Eq. (S.30), which corresponds with the quantum Fisher information order parameter F (t).\nQ\nSeparating out the contributions of the diagonal elements, we obtain\nN\nF (t)=1− 1 (cid:88) ⟨Z ⟩2+ 1 (cid:88) C (t). (S31)\nQ N i t N ij\ni=1 i̸=j\nwhere C (t) = ⟨Z Z ⟩ −⟨Z ⟩ ⟨Z ⟩ is two-point correlation function. The first term corresponds to F = 1,\nij i j t i t j t Q\nindicating the vanishing of local observables, e.g., due to thermalization or noise. The second term is always\nnegative and corresponds to persistent independent spin polarizations in the Z direction, while the third term is\nalways positive and represents the spread of correlations in the system. At early times, the system is polarized\nand uncorrelated, which leads to F = 0. As the system evolves, spins depolarize, making the contribution of\nQ\nthe second term smaller, while simultaneously correlations develop, leading to F > 0 at t > 0. In the ergodic\nQ\nphase, one can expect that both polarization and correlations vanish rapidly, resulting in F → 1 after a few\nQ\ncycles. In contrast, in a DTC regime one may expect that the correlations remain sufficiently strong and the\nsumofthepolarizationandcorrelationtermsremainspositive(assistedbytherelativenumberoftermsbetween\nthe two contributions: N terms in the polarization sum and N(N−1) terms in the correlator sum). The values\nof F > 1 at later times shown in Fig. 3c in the main text appear to be the result of persistent system-wide\nQ\ncorrelations, characteristic of the DTC regime.\nWehavealsoaddedthefollowingpointtothediscussioninthemaintext(seeline234inthecoloredmanuscript):\nQ\nmay not have completely thermalized; see Section IVA of the Supplementary Information for more detail.\nReviewer: “4. In a previous paper Nat. Commun. 15, 8963 (2024), the authors have implemented the 2D time-\ncrystalline order with the planar toric code on a superconducting processor. A proper citation to this paper and some\nother closely related ones should be added.”\nAuthors: WethanktheReviewerforpointingouttheworkofXianget al.,whichexploresadistincttwo-dimensional\nrealization of time-crystalline order via a periodically driven surface code. We have added a citation to this work and\na related topological-ordered DTC citation in our introduction when mentioning other two-dimensional realizations\nof time-crystalline order. We modified the following sentence in line 73 of the colored manuscript:\nBuilding on prior work [3–6], we extend these experimental investigations to a two-dimensional system on a\n“decorated” (heavy) hexagonal lattice and introduce spin-flip coupling.\nReviewer: “5. There are some section numbers missing in the manuscript. For instance, in the Figure 2 caption,\nit reads, ”Noise recovery is performed using the renormalization methods described in Sections ? and III of Supple-\nmentary Information.” Similarly, in line 128, it states, ”...be found in Supplementary Information, Sections I and\n?.”\nAuthors: We thank the Reviewer for spotting these typos. In the updated version, we have fixed these issues.\nReviewer: “The manuscript reports the experimental realization of a two-dimensional discrete time crystal (DTC)\nwith anisotropic Heisenberg coupling in a decorated hexagonal lattice system. By using 144-qubit processor experi-\nments and two-dimensional tensor network simulation techniques, the authors successfully reveal a rich phase diagram\nencompassing spin-glass, ergodic, and time-crystalline phases in a two-dimensional system governed by anisotropic\nHeisenberg interactions with periodic driving. Experimental results demonstrate that the synergistic interplay between\ndisorder and periodic driving can maintain stable quantum dynamical behavior in two-dimensional systems. They\nintroduce the spin-flip coupling to construct the anisotropic Heisenberg model, which can be obtain more interesting\nphenomenon compared to Ising model. They observed the fully polarized state demonstrates a stable subharmonic re-\nsponsethatresemblesthebehaviorofquantummanybodyscarsincludingcat-scarDTCs. Theseachievementsnotonly\nexpand our understanding of nonequilibrium quantum phase transitions but also lay the foundation for investigating\nhow driven systems establish connections between quantum coherence and emergent nonequilibrium thermodynamics.”\nAuthors: We thank the Reviewer for highlighting the significance of our experimental and theoretical work, and we\nrespond to their questions and comments below.\nReviewer: “Here I have a few questions/comments below for the authors to address.\n1. Under the weak coupling condition of ϵ=0.05, how can we demonstrate that the system exhibits characteristics of\nthe anisotropic Heisenberg model rather than the behavior of a perturbed Ising model?\nAuthors: We thank the Reviewer for raising this important question. There exists measurable differences between\nperturbing ϕ from the integrable point ϕ=π/2 (Ising perturbation) and adding a small ϵ (Heisenberg perturbation).\nThe main distinction between these regimes lies in their symmetries. The Ising regime exhibits discrete Z Ising\n2\n(parity) symmetry, whereas the Heisenberg regime possesses continuous U(1) (charge conservation) symmetry.\nThe Ising regime is insensitive to the choice of initial state, a feature also noted in previous works [7, 8]. In contrast,\nthe Heisenberg regime can strongly depend on the initial state. In particular, Fig. 4 demonstrates the difference in\nthe order parameter between the N´eel and polarized states, as well as the clear asymmetry of the phase diagram in\nthe ϵ and ϕ directions for the polarized state (compare to Fig. 3a).\nIn the new version of the manuscript, we have improved the presentation of this difference, see line 85 of the colored\nmanuscript:\nUnlike the Ising model, which is largely insensitive to the choice of initial state [7, 8] and governed by discrete\nZ symmetry, the addition of the Heisenberg terms to the Ising Hamiltonian, even as a perturbation, exhibits\n2\ncontinuous U(1) symmetry. As a result, we find initial states which exhibit strikingly different behavior to the\nprototypical example of a N´eel initial state or a randomly chosen product state. We show that for certain values\nof the spin-flip coupling, oscillations originating from a N´eel state decay rapidly, whereas a fully polarized state\ndisplays a remarkably stable subharmonic response reminiscent of quantum many-body scars [9, 10], including\ncat-scar DTCs [11, 12]. This asymmetry, reflected in the phase diagrams we present in this work, underscores\nthe fundamentally different dynamical behaviors arising from Ising versus Heisenberg interactions.\nTo strengthen our presentation, we have performed a diagnostic of the quantum Fisher information growth between\ndifferentinitialstatesintheIsingandperturbed-Ising(toHeisenberg,ϵ=0.05andϵ=0.1)models. Wepresentthese\nresults in a new Fig. S14 in Supplementary Information (see below).\nWehaveincludedadiscussiontoaccompanythisfigureintheSupplementaryInformation(seeline895inthecolored\nmanuscript):\nForparameterrangeswheretheMBLandDTCregimes,identifiedbythetwo-timecorrelator,aresharedbyboth\nthe N´eel (Fig. 3) and polarized (Fig. 4) initial states, noisy F (t) exhibits logarithmic growth. This behavior\nQ\nis remarkable given that these parameter ranges encompass unperturbed (ϵ = 0) and perturbed (ϵ > 0) kicked\nIsing models. The curvature of F (t) displays a pronounced dependence on the initial state when spin-flip\nQ\ncomponents are introduced (ϵ > 0), as shown in the experimental data of Fig. S14. For Floquet cycles prior to\ncomplete depolarization, the kicked Ising model in Fig. S14a shows little discernible difference between the N´eel\nandpolarizedinitialstatesasinbothcasesitsgrowthoriginatesfromtheunderlyingnoise. Incontrast,atvalues\nofϵthatperturbativelymovethesystemintotheHeisenberg-perturbedregime, Fig.S14b,c, thecurvatureofF\nQ\ndiffers dramatically between these two initial states.\n8\n1.5\n1.0\n0.5\n0.0\n0 2 4 6 8 10\n.ofnI\nrehsiF\nmutnauQ\na b c\nNéel\nPolarized\n0 2 4 6 8 10 0 2 4 6 8 10\nFloquet Cycle\nFigure S14. F (t), given in Eq. (S.30), derived from experimental raw data for the N´eel and polarized initial states, at the\nQ\n(ϕ,ϵ)pointsonthephasediagramlocatedata,(0.0,0.0),b,(0.0,0.05),andc,(0.0,0.1). TheshotnoiselimitatF =1(gray\nQ\nline), is added for clarity.\nWe reference the above discussion in the main text, see line 250 in the colored manuscript:\nDifferences in the rate of of QFI are also observed for the polarized initial state when ϵ is perturbed away from\nthe Ising limit (ϵ = 0), indicating that the symmetries of the underlying Hamiltonian play a significant role\nin the correlations dynamics. Further analysis derived from experimental data is presented in Supplementary\nInformation, Section IVB.\nReviewer: “2. The experimental results in the manuscript are primarily based on the N´eel state as the initial state.\nFor other initial states, can similar discrete time crystal phase behaviors be observed? How do different initial states\ninfluence phase diagram?”\nAuthors: Weexamineinitial-statedependenceinboththepreviousandcurrentversions(seeFig.4b). Inparticular,\nafullypolarizedinitialstate(allspinsalignedalong+z)yieldsaphasediagramthatdiffersmarkedlyfromthatofthe\nN´eel state and displays pronounced robustness to the Heisenberg perturbation. This robustness is expected because\nthe fully polarized product state is an eigenstate of the Heisenberg term, which suppresses its effect on the ensuing\ndynamics. The resulting long-lived response is reminiscent of scar-like behavior, as noted in the previous version of\nthe manuscript.\nBy contrast, other initial states do not exhibit such stability and behave similarly to the N´eel state. To substantiate\nthispoint,weaddedfrequencyspectraforseveralrandomlypolarizedinitialstatestotheSupplementaryInformation.\nHere “randomly polarized” denotes spins with independently sampled orientations (specified in the Supplementary\nInformation, Section IVB). Although these supplementary plots do not use the renormalization scheme detailed in\nour work (and thus represent raw experimental data), the existence of the peak in the frequency plot is relatively\ninsensitive to noise. From these data, we find that the randomly polarized initial states exhibit the same qualitative\nbehavior as the N´eel state for nonzero Heisenberg perturbation ϵ, in contrast to the fully polarized state.\nWe added new Fig. S15 to the Supplementary Information (see below). In addition to the figure, we have added the\nfollowing discussion to Supplementary Information (see line 907 in the colored manuscript):\nThepresenceofsubharmonicresponsesinlocalobservablesunderaperiodicdriveisoneoftheprimaryindicators\nof discrete time-crystalline behavior. A Fourier transform of these observables, such as the two-time correlator\norder parameter defined in Eq. (5), reveals these subharmonic responses through a pronounced peak at characteristic\nfrequency ω =ω /2, where ω is the drive frequency. In the main text, evidence of such a subharmonic\nD D\nresponseisobservedacrossavarietyofparameterregimesandinitialstates,includingboththeN´eel(Fig.3)and\npolarized (Fig. 4) initial states. Notably, the polarized initial state exhibits a robust subharmonic signal over a\nwider range of spin-flip couplings ϵ compared to the N´eel initial state. Other initial states do not display this\nstabilityandbehavenearlysimilarlytotheN´eelinitialstate. InFig.S15,theFouriertransformsofthetwo-time\ncorrelatororderparameter,derivedfromexperimentaldata,areshownforseveralrandominitialstatesalongside\n9\na b c\n0.0 0.5 1.0 0.0 0.5 1.0 0.0 0.5 1.0\n)\n(F\nmrofsnart\nreiruoF\nNéel\nRandom 1\nRandom 2\nRandom 3\nPolarized\nFrequency / D\nFigure S15. The Fourier transform of the order parameter ∆(t), defined in Eq. (5) and derived from experimental raw data\nis shown as a function of the frequency ω in units of the drive frequency ω . The components for ω > ω /2 represent the\nD D\noriginalspectrumfoldedintothisband. Differentinitialstatesareplotted,whichincludetheN´eel(bluesquares)andpolarized\n(orangecircles)initialstates,andthreerealizationsofrandominitialstates(purpleinvertedtriangles,greendiamonds,andred\ntriangles),asdescribedinSectionIVB.Panelscorrespondtopointsonthephasediagram(ϕ,ϵ)locatedata,(0.45π,0.05),b,\n(0.35π,0.15), and c, (0.3π,0.2). Curves for different initial states exhibit a vertical displacement for illustration purposes.\nthosefortheN´eelandpolarizedinitialstates. Eachrandomstateisgeneratedasaproductstateofindependently\nand randomly assigned maximally polarized spins s ∈ {−1,1}. Near the Clifford point in Fig. S15a, the peak\ni\nat ω =ω /2 is pronounced for all initial states, though the polarized initial state shows a larger amplitude than\nD\ntheotherstates. Astheparametersaretunedawayfromtheintegrablepoint,bydecreasingthekick-angleϕand\nincreasing spin-flip strength ϵ, the polarized state retains a half-frequency peak, in contrast to the other initial\nstates, see Fig. S15b,c. In Fig. S15c, as the kick angle approaches the ergodic regime at ϕ = π/4 the original\nhalf-frequency peaks disappears, indicating the cessation of time-crystalline behavior.\nWe reference this expanded discussion in the main text, see line 240 in the colored manuscript:\nWeexpectthatmostinputstateswillexhibitbehaviorsimilartothatoftheN´eelstate. Trendsinthecorrelator\norder parameterof Eq.(5) derived fromexperimental data, for severalother random initialstates, are presented\nin the Supplementary Information, Section IVB. A notable exceptions to this expectation are fully polarized\nstates or states close to them by Hamming distance.\nReviewer: “3. The manuscript repeatedly mentions the impact of hardware noise on experimental results, but lacks\nanalysis of its specific sources and quantitative effects on the outcomes. This needs to be supplemented and improved.”\nAuthors: We thank the Reviewer for this suggestion. We would like to emphasize that the previous version of the\nmanuscript does discuss the role of noise. For example, the Methods section models the effect of dephasing and\namplitude damping in Eq. (8) and suggest an error mitigation scheme. We also have a Supplementary Information\nSection I that provides comprehensive analysis of the device noise.\nWerecognize,however,thatwemaynotadequatelyguidethereadertotherelevantpartsofthepaper. Inthecurrent\nversion, we have improved this aspect. We now provide explicit references in the main text and Methods section to\nthe Supplementary Information where the noise analysis is presented.\nInaddition, wehaveaddedanewfiguretoSupplementaryInformationshowingtheeffectofnoiseontheactualorder\nparameter data, derived from an analysis of Clifford dynamics. We hope that this expanded discussion provides a\nclear and sufficient understanding of the role of noise in our results.\nWehavemadethefollowingchangesinresponsetotheReviewer’scomment. InSupplementaryInformation, wehave\nadded a figure illustrating the quantitative effect of noise:\nWe have accompanied this figure with the comments in Supplementary Information (see line 605 in the colored\nmanuscript):\nAn independent, observable-dependent quantification of noise can be obtained by evaluating the observable in\nthe Clifford regime. At this point, the noiseless values of |∆(t)| and χ(t) in Eqs. (5) and (6) both equal unity,\n10\n1.00\n0.75\n0.50\n0.25\n0.00\n0.25\n0.50\n0.75\n1.00\n0 10 20 30 40 50\nFloquet cycles\nretemarap\nredrO\nexperiment\noffset\ndecay curve\nFigure S3. Quantitative noise analysis at the Clifford point (marked as Clifford 2 in Fig. 3 in the main text), used as a\nreference for error mitigation. In noiseless simulations, the order parameter alternates as ∆(t)=(−1)t, where t is the number\nofFloquetcycles. Withnoise,theobservedsignaltakestheform∆ (t)=∆(t)f(t)+δ,wheref(t)accountsfordecoherence\nnoisy\nnoise (dashed blue curves) and δ is an offset due to amplitude damping noise (red line); see Eq. (S.5). The solid black line\ncorresponds to zero signal. Approximating f(t) ≈ Ae−t/T, with T the characteristic decay time, allows extraction of the\nparameters A = 0.973, T = 14.5, and δ = −0.063 (even cycles) and δ = −0.019 (odd cycles), providing a quantitative\ncharacterization of the noise. These values are subsequently used to mitigate errors at non-Clifford points (see Section III).\nwhich allows us to isolate the impact of noise. An example of this analysis is shown in Fig. S3 for ∆(t) at the\nClifford point 2 in Fig. 3 in the main text (ϵ=0, ϕ=π/2). The data reveals two distinct effects of noise: (i) a\ndecay of the signal due to coherence loss, and (ii) a systematic offset arising from the non-unital component of\nthe noise (amplitude damping). Fig. S3 indicates a decoherence time of approximately T ≈ 15 Floquet cycles,\ncorresponding to D ≈135 two-qubit gate layers. The offset is approximately δ ≈−4×10−2, which requires the\ndedicated error-mitigation strategy described in Section III.\nIn Methods, we have added a reference to this figure (see line 462)\nThiseffectcanbeillustratedbyFig.S3intheSupplementaryInformationwheretheparametersoftheamplitude\nattenuation function and the offset are estimated from the actual data.\nIn the main text, we have added a reference to the noise model in line 123 of the colored manuscript:\nFurther details regarding the device architecture and implementation, including comprehensive description of\nthe device noise, can be found in Methods and Supplementary Information, Section I.\nReviewer: “4. The manuscript selects ϕ = 0.45π as a representative DTC parameter (Fig 2a). Why choice this\nangle? Was this angle systematically optimized to maximize subharmonic response robustness?”\nAuthors: We thank the Reviewer for raising this question. We use ϕ = 0.45π,ϵ = 0.05 as our operating point\nin Fig. 2 because it balances two limiting regimes. It is sufficiently detuned from the suspected phase-transition\npoint—where the DTC order parameter is strongly suppressed—and from the Clifford (integrable) point—where the\ndynamics become trivial. At ϕ = 0.45π,ϵ = 0.05, the subharmonic response remains robust while entanglement\ngrowth is moderate, enabling reliable tensor-network benchmarks and a transparent comparison between experiment\nand classical simulation. For these reasons, we regard ϕ = 0.45π as the optimal working point for resolving DTC\ndynamics in our setup.\nWe have added the following clarification to line 164 of the colored manuscript:\nChosen to lie between the suspected transition to the ergodic regime and the Clifford (integrable) point, this\npoint yields a robust DTC signal while keeping entanglement moderate and enabling reliable tensor-network\nsimulations.\n11"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have satisfactorily addressed my concerns and those of the other reviewers. I therefore maintain my very\npositive view of this work (consistent with the other referee reports).", "reviewer_2": "I went through the revised manuscript, the reports from other referees and the author’s replies to all three referees. In this\nround, the authors have made considerable efforts to improve the manuscript and address all the referees’ comments and\nsuggestions. In particular, they have added the system size 3\\times 3 results to Fig. 3b to obtain better finite-size scaling. As\nfar as I am concerned, they have fully addressed my concerns and their responses are satisfying. The revised manuscript\nhas been improved substantially. I am happy to recommend its acceptance in Nature Communications with the current\nversion.", "reviewer_3": "I thank the authors for their detailed and constructive responses to my previous comments, and for the corresponding\nrevisions made to the manuscript. They have thoroughly addressed all the points raised, providing additional clarifications,\ndata, and analysis where needed. These revisions have strengthened the manuscript and improved its clarity and rigor.\nThis work represents a significant experimental and theoretical achievement in realizing and characterizing a two-\ndimensional discrete time crystal with anisotropic Heisenberg interactions on a noisy quantum processor. The successful\ndemonstration of this non-equilibrium phase, along with the comprehensive phase diagram revealing spin-glass,\nergodic, and time-crystalline regimes, provides valuable insights into the stability of quantum coherence in driven many-body\nsystems. The observed scar-like behavior in the polarized state particularly advances our understanding of weak ergodicity\nbreaking in two dimensions. With the issues I raised now resolved, I find the manuscript complete and interesting for\npublication. I am pleased to recommend its acceptance."}, "author_response": ""}]

Decision letter