Stabilization of Discrete Time-Crystaline Response on a Superconducting Quantum Computer by increasing the Interaction Range

The simulation of complex quantum many-body systems is a promising short-term goal of noisy intermediate-scale quantum (NISQ) devices. However, the limited connectivity of native qubits hinders the implementation of quantum algorithms that require long-range interactions. We present the outcomes of a digital quantum simulation where we overcome the limitations of the qubit connectivity in NISQ devices. Utilizing the universality of quantum processor native gates, we demonstrate how to implement couplings among physically disconnected qubits at the cost of increasing the circuit depth. We apply this method to simulate a Floquet-driven quantum spin chain featuring interactions beyond nearest neighbors. Specifically, we benchmark the prethermal stabilization of the discrete Floquet time-crystalline response as the interaction range increases, a phenomenon never observed experimentally. Our quantum simulation addresses one of the significant limitations of superconducting quantum processors, namely, device connectivity. It reveals that nontrivial physics involving couplings beyond nearest neighbors can be extracted after the impact of noise is properly taken into account in the theoretical model and consequently mitigated from the experimental data.


I. INTRODUCTION
Quantum computers are a revolutionary technology that have the potential to transform our society by solving problems that classical computers cannot [1].However, these machines are still subject to uncontrollable noise and errors that limit their performance, which are far from the threshold required for error correction.Despite these limitations, recent progress in the realm of noisy intermediate-scale quantum (NISQ) devices represents an exciting opportunity for many-body physics by introducing new laboratory platforms with unprecedented control and measurement capabilities [2].Quantum simulation of the dynamics of more and more complex quantum many-body systems is expected to be one of the most promising short-term goals of NISQ quantum computing devices, with intriguing applications in diverse areas ranging from quantum chemistry [3][4][5] and material science [6] to high-energy physics [7].
Various experimental platforms have been tested for quantum computing, among the others we can cite: trapped ions, [8][9][10][11] neutral Rydberg atoms [12][13][14], coherent photons [15,16], nuclear spins in molecule [17,18], NV centers [19,20] and superconducting qubits [21,22].Each of them has its own advantages and drawbacks [23].In this paper, we focus on superconducting quantum processors.Superconducting qubits are relatively easy to fabricate and can be densely packed, allowing for the construction of large-scale quantum computers.This makes them a promising platform for scaling up quantum computing applications [24].Moreover they can be manipulated with a wide range of microwave frequencies, making them versatile and flexible for implementing various quantum gates [21].
Thanks to this flexibility, the number of quantum simulations implemented on noisy superconducting devices has steadily risen in recent years, also thanks to the possibility to easily access these machines remotely, allowing to benchmark a number of phenomena which were not or were very little experimentally corroborated before.Among this plethora of studies we may mention the following: the observation of disorder stabilized discrete time crystal phases, [25,26] the realization of topologically ordered states, dynamical topological phases and topological edge states [27][28][29].One should also cite the observation of Leggett-Garg's inequalities violations [30], the validation of dynamical scalings [31,32] and several studies in the context of quantum thermodynamics [33,34].On the other hand, the performance of superconducting NISQ devices is limited by the presence of various sources of noise and decoherence, whose impact grows with the depth and complexity of the quantum circuit realized, limiting the investigation of non-local effects and complex geometries.
In this context, one main limitation of superconducting NISQ devices is their extremely limited connectivity, since superconducting qubits are typically arranged in a one or two-dimensional grid with nearest-neighbor connectivity, making it challenging to implement quantum algorithms that require long-range interactions [1].In this paper, we aim to advance the field of digital quantum simulation on superconducting quantum hardware by investigating the possibility of reproducing the dynamics of systems with couplings beyond nearest neighbors.To achieve this, we utilize the universality of the quantum processor native gates to implement couplings among physically unconnected qubits.While the depth of the resulting quantum circuit increases with the effective range of the interaction, we show that careful consideration of gate noise, measurement errors, and statistical errors enables the removal of their effects from the raw results.The resulting error-mitigated data closely reproduce the theoretical expectations.
More specifically, we implement the quantum simulation of a Floquet-driven quantum spin chain featuring interactions beyond nearest neighbors on IBM quantum superconducting processors.Indeed, the quantum circuit structure utilized by IBM quantum computers is well suited for implementing discrete Floquet driving protocols [2], making it a natural choice for such applications [26].
Our focus is on the stabilization of discrete Floquet time-crystalline response as the interaction range increases.Discrete Floquet time crystals (DFTCs) are nonequilibrium many-body phases of matter that display a novel form of spatiotemporal order.In particular, in such phases the discrete time translation symmetry of the Floquet driving is broken and an order parameter exhibits persistent oscillations with a period which is an integer multiple of the period of the drive [55][56][57][58][59][60][61][62][63][64].The possibility of generating a DFTC in clean systems has been studied in the context of long-range interacting models, and our quantum simulation on IBM quantum processors constitutes its first experimental benchmark.Our results demonstrate the potential of superconducting quantum computing platforms to simulate quantum systems featuring interaction ranges going beyond the limits imposed by hardware connectivity and offer insights into the fundamental physics of long-range systems.

II. THE MODEL
The kicked Ising spin chain is a prototypical model for the investigation of Floquet driven quantum systems, widely studied from a theoretical point of view [65][66][67][68][69].In this paper, we consider a driven quantum spin chain described by a time-dependent Hamiltonian of the form where the time dependence is generated by a time periodic driving with period T of the transverse magnetic field h(t).The driving takes the form The effect of this impulsive magnetic field applied at integer multiples of the driving period t = nT is to impose a global rotation of every spin by an angle 2ϕ along the x-axis.Accordingly, the Floquet dynamics is obtained by periodically intertwining the evolution generated by the Ising Hamiltonian at zero transverse field with the instantaneous kick operator, The resulting evolution operator for a single step of the Floquet protocol reads The system is initialized at t = 0 in the fully polarized state with positive magnetization along the ẑ direction |ψ(0)⟩ = | . . .↑↑↑ . . .⟩ = | . . .000 . . .⟩, where | ↑⟩ and | ↓⟩ denote the eigenstates of the Z Pauli matrix with eigenvalues +1 and −1 respectively.In our case, these eigenstates correspond to the computational basis of the quantum processor, with the convention | ↑⟩ = |0⟩ and | ↓⟩ = |1⟩.The simplest realization of the time-crystalline spatiotemporal order is obtained by taking the kick operator K ϕ to rotate each spin by an angle π around a transverse axis x.In this case, the kick operator is given by As a result, the time-evolved state after n kicks, |ψ(n)⟩ = U n F |ψ(0)⟩, exhibits a sequence of perfect jumps between the | . . .↑↑↑ . . .⟩ and | . . .↓↓↓ . . .⟩ states, leading to a persistent non-vanishing value of the order parameter in both space and time.The order parameter is given by (7) This is the simplest example of a subharmonic response, where the period of the order parameter evolution is twice the period of the Floquet driving.However, this behavior depends on the finely tuned choice of the kick angle 2ϕ = π.To observe a proper discrete time-crystalline phase of matter, the spatiotemporal order must be stable to sufficiently weak perturbations of the Hamiltonian parameters ϕ = π/2 + ϵ, in the thermodynamic limit N → ∞.This condition is generally not satisfied as the presence of the external driving leads to the exponential decay of the magnetization, ruling out long-lived oscillations.Protecting ordering against relaxation requires a mechanism to control the impact of dynamically generated excitations [64].
In clean systems, the possibility of generating a DFTC has been studied in the context of long-range interacting models [65][66][67][68][69][70], where the interaction between different lattice sites decays as a power law.However, for any finite R and in the absence of disorder, the system magnetization exhibits an exponential decay with the number of Floquet steps: The decay rate γ ϵ,R goes to zero as the perfect kick case is approached, i.e., for ϵ → 0.Moreover, as shown in Ref. [68], γ ϵ,R is deeply affected by the interaction range.In the small ϵ limit, we have that Therefore, increasing the interaction range exponentially enhances the order parameter lifetime.This difference in decay rate should already be apparent when comparing the nearest neighbor R = 1 and the next-to-nearest neighbor R = 2 cases.One of the main results of our digital quantum simulation is to demonstrate this increase in the order parameter lifetime.

III. QUANTUM CIRCUIT IMPLEMENTATION OF THE FLOQUET DYNAMICS
The aim of our quantum simulation is to implement the dynamics of the Floquet-driven quantum spin chain with tunable interaction range, introduced in the previous section, on an IBM quantum processor.Specifically, we utilize the ibmq mumbai 27-qubit processor, whose topology is depicted in Fig. 1a (further technical details can be found in the Supplemental Material (SM) Appendix).Our quantum circuit is optimized using the available connectivity and native gates of the device, including the controlled-NOT gate (CNOT), the identity gate ID, rotations along the z axis R Z , the NOT gate X, and the SX = √ X gate (see Appendix A).
We notice that the Floquet unitary evolution operator at stroboscopic times t = nT can be obtained by applying the unitary operator corresponding to each Floquet step U F n times, i.e., U (nT ) = (U F ) n .Importantly, no Trotter approximation is required, which is a significant advantage of Floquet drivings, making them well-suited for quantum circuit implementation [2].Furthermore, the kicked Floquet protocol of interest can be further decomposed into the successive application of the kick operator K ϕ and the Ising evolution operator V (see Eq. ( 5)).The former can be expressed in terms of singlequbit gates, corresponding to local rotations along the x-axis, and the latter can be written as a product of mutually commuting unitaries that connect pairs of qubits at progressively larger distances as the interaction range is increased, i.e., respectively.Each V r can be implemented by applying the general method to effectively realize r-range interactions introduced in the previous section.
The quantum circuit corresponding to a single Floquet step is shown in Fig. 1(b), where blue gates represent nearest-neighbor Ising interactions V j,j+1 , red gates represent Ising interactions beyond nearest neighbors V j,j+r , and green gates represent the final kick rotation K ϕ applied equally to each qubit.In particular, as shown in Fig. 1(c), the unitary operator associated to nearestneighbor Ising interactions can be decomposed in terms of the elementary gates as On the other hand, the limited processor connectivity does not allow for a simple decomposition of r-range Ising interactions.The idea to overcome this limitation is to exchange the qubit states by applying a sequence of SWAP gates among the couples of physically connected qubits that lie between j and j + r.By doing so, the initial state of qubit q j+r is effectively encoded in qubit q j+1 .Specifically, we achieve this by applying the gate sequence Next, we apply V j,j+1 on the two connected qubits q j+1 and q j .Finally, we need to bring back the state encoded in qubit q j+1 to the r-neighbor qubit q j+r .This is achieved by applying the inverse sequence of SWAP gates S † r .Summarizing, we obtain the quantum circuit identity shown in Fig. 1(e), reading This enables us to realize the desired tunable-range interactions among physically unconnected qubits.However, there is a trade-off involved: the implementation of these interactions requires the insertion of 2(r − 1) additional SWAP gates into the quantum circuit.Since each two-qubit gate typically introduces noise, it becomes imperative to optimize our quantum circuit for each Floquet step of the dynamics (as depicted in Fig. 1(b)).This optimization involves breaking down each operation into the native gates of the quantum device and strategically reordering our quantum circuit to minimize the involvement of two-qubit gates.A comprehensive description of this procedure is provided in Appendix A.
Furthermore, we will need to mitigate effect of the noise as the interaction rate increases.We address this problem in the following sections.

IV. THE ROLE OF NOISE AND NOISE MITIGATION
The analysis of the raw experimental data clearly demonstrates that the decay of magnetization is predominantly influenced by the effect of noise, as illustrated in Fig. 2. The figure depicts the absolute value of the average magnetization |m z | as a function of the stroboscopic time n for nearest neighbor (R = 1, in blue) and next-to-nearest neighbor (R = 2, in red) interactions.The raw experimental data (triangles in Fig. 2) are obtained by running n repetitions of the quantum circuit corresponding to a single Floquet step U F , as depicted in Fig. 1(b), on the ibmq mumbai quantum processor using N = 18 qubits.At the end of each quantum evolution, a projective measurement of each qubit in the Z basis is performed.To collect sufficient statistics, the experiments for each value of n and R are repeated over a sample of size N = 2 13 , allowing us to compute the sample average ⟨Z i ⟩ over the measurement outcomes.Finally, the spatial average of the magnetization over different sites of the processor is computed as where N = 18 in our case.
Estimating the statistical error from multiple instances of each quantum simulation to evaluate the sample mean E(m z ) and standard deviation σ(m z ) is not feasible due to the time taken to produce each magnetization estimate.Instead, we rely on the statistical tool of bootstrapping, which is further detailed in Appendix C, to generate resampled data from the empirical measurement outcomes.Since the bootstrapped data conforms to the central limit theorem, we may assume normality and evaluate E(m z ) and σ(m z ) from these artificially generated samples of data.The results for the statistical averages E(m z (n)) are represented by squares in Fig. 2, while we use two standard deviations, 2σ(m z (n)), as statistical errors, depicted as error bars in the plots.Notably, we observe that the statistical error increases with the number of Floquet steps involved in the dynamics n.This can be understood by a simple statistical argument: we are trying to sample a quantity, the modulus of the magnetization |m z |, which exponentially decreases with n.Consequently, the resolution with which we can estimate this quantity deteriorates as m z approaches the value m z ∼ e −γn ∼ 1/ √ N , i.e., the statistical uncertainty due to the finite size of the sample increases as we approach the stroboscopic time n ∼ (1/2γ) ln N .
The decay of magnetization with stroboscopic time n can be described by an exponential fit |E(m z (n))| = ae −bn , which is obtained using a weighted least squares regression method.This approach accounts for points with high statistical uncertainty, penalizing them in the extrapolation.The resulting exponential decay is depicted as a dashed line in Fig. 2, showing a rapid decline with increasing n.Notably, the decay rate is more pronounced for next-to-nearest neighbor interactions (R = 2) compared to nearest neighbor interactions (R = 1).This discrepancy can be attributed to the fact that the quantum circuit implementing next-to-nearest neighbor interactions involves more gates, resulting in larger noise effects.
In order to effectively simulate desired physical phenomena in a quantum system, it is crucial to account for and mitigate the detrimental effects of noise.Real-world quantum hardware is susceptible to various sources of errors, such as noisy gates, environmental decoherence, and spurious time dependence of circuit parameters [2].To explicitly model these errors, a common approach is to consider one-and two-qubit depolarizing channels that act on the system's state ρ.Specifically, after each singlequbit gate acting on qubit i, the single-qubit channel ϕ 1q i is applied, while after each two-qubit gate on bond (i, j), the two-qubit channel ϕ 2q i,j is applied.These channels are defined as [2,71] where σ 1,i = X i , σ 2,i = Y i , and σ 3,i = Z i are the Pauli matrices for qubit i, and σ α,i and σ β,j are the corresponding matrices for qubits i and j, respectively.By studying the dynamics of the Z i operators under these depolarizing channels, we can estimate the magnetization decay rate induced by the noisy gates.
To isolate the effect of noise, we consider the case of perfect kick dynamics with ϵ = 0.Under this condition, Z i is invariant under the two-qubit Ising interaction gates and simply acquires a minus sign under the π rotation around the x-axis.However, after each two-qubit gate, Z i decays under ϕ 2q i,j as and after each single-qubit gate as Overall, Z i decays to −e −γ dep Z i , over one noisy Floquet step with perfect kicks, with γ dep given by where Q 2q,R and Q 1q,R are the number of two-qubit and single-qubit gates involved in a Floquet step quantum circuit with R-neighbor interactions.A naive estimate of these numbers based on the general method previously introduced would yield Q 2q,R = 2R 3 + 2 and Q 1q,R = 2R + 5, respectively (see Appendix A for more details).However, for the specific case of the kicked Ising model considered in our quantum simulations, we were able to optimize the quantum circuits corresponding to R = 1, 2 Floquet steps, reducing the number of two-qubit native gates, involved in the quantum circuit longest path, to Q 2q,R = 9R 2 − 14R + 7, with a reduction from cubic to quadratic range dependence.In particular for R = 1, 2 we have More details on the quantum circuit optimization strategy are provided in Appendix A. Another source of noise arises from the finite decoherence time T 1 of the qubits, which introduces an additional time scale contributing to the magnetization decay.Taking into account all the contributions, we can estimate the decay rate of magnetization for a Floquet step with imperfect kicks of an angle ϕ = π/2 + ϵ to be approximately given by where τ R represents the time required to practically implement the Floquet step on the quantum hardware.This can be estimated as where τ 1q and τ 2q denote the time needed to execute each single-qubit and two-qubit gate, respectively, while τ m represents the readout time required for measurements.
Estimates of these quantities, as obtained from the engine calibration, are provided in the Supplementary Material.
A third source of errors arises from readout errors, which can be modeled as a stochastic process where the outcome of a qubit-state measurement (in the Z computational basis) is randomly flipped with a probability of p m away from its correct value [2].Specifically, if we define the probability that qubit i points up (down) at time n as Π ± = ⟨(1 ± Z i (n))/2⟩, then the result of the noisy measurement process is Z i = ±1, with a probability of Π± (n) = Π±(1 − p m ) + Π ∓ p m .Accordingly, the estimate for the expectation value of Z i becomes Hence, averaging over positions yields mz = (1 − 2p m )m z , i.e., a damping by a time-independent and range-independent overall prefactor C m = (1 − 2p m ).
The inclusion of noise in our model provides a compelling explanation for the rapid exponential decay of magnetization, as observed in Fig. 2.Moreover, by inserting the estimated values of the parameters p 1 , p 2 , τ , and T 1 , which were extracted from the calibration data provided by IBM and detailed in the SM Appendix, we find that the calculated decay rate is in good agreement with that obtained from fitting the experimental data with a stroboscopic time dependence of the form predicted by our theoretical model, This understanding of the noise effect justifies our exploration of the possibility of mitigating it through a technique called zero noise extrapolation.The basic idea of zero noise extrapolation is to intentionally increase the noise level by amplifying the depth of the quantum circuits by a factor of s through a procedure called circuit folding (see Appendix B for a more detailed description).Subsequently, we perform quantum simulations for different noise scales s, and for each scale, we extract the magnetization decay rate from the measured data.Our noise model then allows us to theoretically estimate the decay rate at noise scale s as Accordingly, a linear fit of the measured decay rates with respect to the parameter s enables us to separate the contribution coming from the noise, γ noise = γ dep,R + τ R /T 1 , from γ ϵ,R , which represents the decay rate due to the internal system thermalization that destroys the time crystalline order at finite ϵ, and should be stabilized by the presence of longer-range interactions.More precisely, γ ϵ,R is obtained as the zero noise extrapolation of the decay rate in (24), γ ϵ,R ≈ Γ s=0,R .The results of this procedure are shown in Fig. 3, where the measured decay rate is plotted as a function of the noise scale s.To estimate Γ s,R and its uncertainty δΓ s,R , we first estimate the magnetization as a function of the stroboscopic time n at different values of s and R from the measured data, along with the corresponding statistical uncertainty from the standard deviation obtained through the statistical bootstrap method, E(m z )±2σ(m z ).Then, the decay rate and its uncertainty are obtained through the exponential fit In particular, the exponential fit is performed using weighted least squares regression, and the last two points with R = 2 and s > 1.5 are excluded from the fitting data (empty points in Fig. 3).This exclusion is justified by the fact that the decay rate for these points falls within the range of 0.2 < Γ s>1.5 < 0.8, and thus, the magnetization can be reliably estimated only for stroboscopic times n < n * ≈ (1/2γ) ln N , where 6 < n * < 22.Therefore, not all the time steps 1 < n < 16 considered in the exponential fit of m z (n) from which we extracted this decay rate are within the reach of our statistical resolution.The difficulty of establishing a reliable bootstrapestimated value confirms this phenomenon, as shown in Fig. 3, where the statistical error bars for these points are significantly larger than those for the other points, indicating the challenge of obtaining a trustworthy value for the magnetization in this regime.Despite the failure of the bootstrap procedure, we include these data as empty points in the plot for completeness, noting that the corresponding error bars are sufficiently large that the resulting fit is still compatible with these unreliable values within ±2σ(m z ).
Remarkably, upon extrapolation to the zero noise limit, the decay rate of the R = 2 case is found to be smaller than that in the nearest neighbors case R = 1.Specifically, we obtain Γ 0,2 ± δΓ 0,2 < Γ 0,1 ± δΓ 0,1 . ( Most significantly, we find that the decay rate analytically predicted from the theoretical model, γ ϵ,R ≈ ϵ 2R+1 , is compatible with the extrapolated values within the estimated uncertainty, i.e., indicating that the extrapolated decay rate is consistent with the theoretical expectations within the statistical uncertainty δΓ 0,R , which has been estimated by extrapolating δΓ s,R to s = 0.

V. CONCLUSION
In this paper we have demonstrated the potential of superconducting quantum hardware for advancing the field of digital quantum simulation by investigating the possibility of implementing quantum dynamics of systems with couplings beyond nearest neighbors.We have utilized the universality of the native gates in quantum processors to implement couplings among physically disconnected qubits and carefully mitigated the effects of gate noise, measurement errors, and statistical errors from the raw results.Our focus has been on the stabilization of discrete Floquet time-crystalline response as the interaction range increases, and we have implemented, on IBM quantum superconducting processors, a quantum simulation of a Floquet-driven quantum spin chain with interactions beyond nearest neighbors.
Our results, as shown in Fig. 2, reveal that the magnetization dynamics under a kicked Floquet driving exhibits a fast exponential decay dominated by noise in the raw data, with a faster decay rate for larger interaction ranges due to the increased depth of the quantum circuit.However, after applying the zero noise extrapolation procedure, we were able to separate the role of noise from the true decay caused by the dynamical generation of excitations in the system during the Floquet driving.The mitigated data for the magnetization decay rate, in Fig. 3, show a clear trend of slower decay in the zero noise limit compared to the raw data, indicating the effectiveness of our error mitigation approach.
Furthermore, we have estimated the statistical uncertainty of the mitigated data, represented by the error bars in Fig. 3, which grows with the noise level as expected.Importantly, we have observed that the regions corresponding to the zero noise extrapolated values of the magnetization decay rate for different interaction ranges do not intersect within the considered values of the parameters, indicating that the effects of noise have been effectively removed from the data.
We have also compared our results with the theoretical expectations, and found that the magnetization decay rate analytically predicted from the theoretical model, γ ϵ,R ≈ ϵ 2R+1 , is compatible with the final results of our quantum simulation within the estimated uncertainty.This agreement between theory and experiment provides evidence for the validity of our approach in simulating quantum systems with interaction ranges beyond the limits imposed by hardware connectivity.
It is important to highlight that,in principle, the quantum circuit we used for implementing beyond nearestneighbor interactions on quantum superconducting computers has the potential to enable tunable interaction ranges.However, there is a fundamental limitation that led us to restrict our quantum simulations to cases where R = 1, 2. This is due to the fact that increasing the effective range of interaction will inevitably require a deeper quantum circuit, which in turn increases the impact of noise.Practically speaking, we anticipate that already at R = 3, 4, the noise level may be significant enough to prevent a reliable estimate of the magnetization decay.Another way of looking at this is that for the result to be within our resolution, the noise level should be at most comparable with the decay rate γ ϵ,R ≈ ϵ 2R+1 .For R ≥ 3, this noise level could already exceed the threshold for fault-tolerant quantum error correction and therefore surpass the capabilities of the available NISQ devices.
Nevertheless, there are several steps that can be taken to overcome this problem.For example, the noise level of quantum devices is expected to decrease as the performance of available quantum computers improves.Additionally, classical algorithms can be used to further optimize the quantum circuit we proposed, reducing the number of two-qubit gates involved.Furthermore, it would be instructive to benchmark our results on different experimental platforms that naturally allow for the implementation of long-range interactions, such as trapped ions or Rydberg atoms devices.These exciting problems are beyond the scope of the present work and hence we leave them for future research.
Another crucial aspect to emphasize is that in our quantum simulation, we utilize 18 functional qubits out of the total 27 nominal qubits of the ibmq mumbai device.This distinction is significant since the number of functional qubits is often smaller in other applications.For instance, in quantum chemistry applications, as described in Ref. [72], noise mitigation strategies are employed to simulate molecules whose Hamiltonian can be encoded in only three qubits in order to achieve the desired chemical accuracy.This exemplifies the suitability of Floquet dynamics, similar to those analyzed in our work, for implementation on NISQ quantum computers.
In conclusion, our quantum simulation on IBM quantum superconducting processors has demonstrated the potential of these platforms for simulating quantum systems with couplings beyond nearest neighbors and has offered insights into the fundamental physics of long-range systems.Our error mitigation approach has been effective in removing the effects of noise and measurement errors from the raw results, and the mitigated data are in good agreement with theoretical expectations.This work opens up new possibilities for studying quantum systems with long-range interactions and paves the way for further advancements in the field of digital quantum simulation on superconducting quantum hardware.[V j,j+l , V j,j+r ] = 0 for all l, r, and [V j,j+1 , SWAP j,j+l ] = 0 to maximize the number of adjacent V j,j+1 SWAP j,j+l , and thereby increase the number of CNOT gates that cancel out.for a circuit implementing a sequence of Ising interactions of ranges r = 1, . . .R, we can cancel up to 2(R − 1) CNOT gates using this trick.The depth of each subcircuit of this form is then given by To realize a kicked Ising model with interaction range up to R in a chain of N qubits, we can divide the N qubits into subsets of size 2R that can be processed in parallel.To compute the circuit depth, which refers to the number of operations in the longest path, we can focus on only one subset at a time.Each subset contains R subcircuits with interaction ranges r = 1, . . ., R, following the form shown in Fig. 4(b).The depth of each of these subcircuits is 6R − 4. The remaining R subcircuits include interactions of range r = 1, . . ., l, where l varies from l = 1 to l = R−1, corresponding to a depth of 6l−4 for each circuit.By summing up all the contributions, we obtain the optimized number of CNOT gates as: Finally, we observe that the last sequence of SWAP gates in the circuit shown in Fig. 4b is only necessary if we need to apply different gates on different qubits after that.If this is not the case, we can simply substitute the SWAP gates with a relabeling of the qubit numbers, which must be taken into account when reading the final measurement outcomes.This fact allows us to eliminate (R − 1) SWAP gates and 3(R − 1) CNOT gates in the last subcircuit of this form.Therefore, we obtain: As a final remark, we note that for large values of maximum range R, and hence large circuit complexity, addi-tional simplifications of the circuit may be possible by using optimized relabelings of the qubit numbers during the evolution, which can increase the number of parallelizable operations.However, such an optimization strategy is circuit and range dependent, and can only be carried out numerically or in an approximate manner.On the other hand, for the case of R = 2 that we considered in our quantum simulation, we can claim that our circuit is optimal with respect to the number of CNOT gates involved.
behavior under different noise conditions.

Appendix C: Statistical Bootstrapping
We utilize the statistical technique of bootstrapping to quantify the uncertainty in our magnetization estimates.In an ideal scenario we would repeat quantum experiments multiple times to obtain a comprehensive understanding of the "true" magnetization distribution.However, this approach is impractical due to the significant time required for each magnetization estimate.Instead, we conduct the experiment once and generate resampled measurement data from the empirical distribution using bootstrapping, a widely used statistical technique.This method makes the statistical analysis very convenient, and it is then becoming a common practice to estimate the statistical errors in digital quantum simulations [72].
Let us assume we perform an N -shot quantum experiment and obtain a collection of N outcomes.Each measurement outcome is represented as a string of 0s and 1s, denoted as Z 1,a . . .Z N,a , where Z i = 0, 1, N is the number of measured qubits, and the index a labels the different outcomes (a = 1, . . ., N ).The magnetization associated with each string can be computed by averaging over the qubits as follows: This gives us the set of magnetization values m z,a with a = 1, . . ., N .We define the empirical magnetization distribution P 1 (m z,a ) as the histogram of the m z,a set.The average over this empirical distribution, denoted as m z , corresponds to the experimentally obtained quantum expectation value on the final state of the system and can be expressed as The bootstrapping approach involves resampling from the empirical measurement distribution P 1 (m z,a ).We sample elements from the set m z,a (or, equivalently, from the set of strings {Z 1,a . . .Z N,a }) N times to create a new set of measurement outcomes, and from this, a new empirical distribution P 2 (m z,a ).We repeat this process as many times as possible given the available computational resources, say M repetitions, to obtain a set of distributions P 1 , P 2 , . . ., P M .From each of these distributions, we can compute the average m α z with α = 1, . . ., M , and from the histogram of the set of averages, we obtain their distribution Π(m α z ).Since each resampling is independent, the distribution of averages should tend to a Gaussian in the large M limit, according to the central limit theorem.Accordingly, we can define our estimator for m z and its statistical error as the average of the

Π(m
and its standard deviation, This method enables us to obtain error bars in Figs. 2 and 3 of the main text as E(m z ) ± 2σ(m z ). Figure 5 shows, as an example, the distributions Π(m α z ) obtained through M = 1000 resamples of the measured data of three quantum simulations with range R = 1, noise scale s = 1.4,and number of Floquet steps n = 0, 5, 8, respectively.These are compared with Gaussian distributions with the same mean and standard deviation, finding good agreement.We notice that, as expected, the mean E(m z ) is smaller for a larger number of Floquet steps n, signaling the magnetization exponential decay.Moreover, distributions at later stroboscopic times become broader, signaling the growth of the statistical error due to the fact that we are trying to sample a quantity which is exponentially decaying with n.

Supplementary Materials
Appendix D: Calibration details of the quantum processor In this section, we provide all the calibration details of the ibmq mumbai backend, which is one of the 27-qubit IBM Falcon processors used in our experiments.Figure 6 shows the topology of the processor and the qubit labeling numbers.We present the corresponding coupling map for two-qubit gates in Figure 7, which includes the CNOT gate error p 2 (panel a) and the CNOT gate length τ 2q (panel b) for each pair of physically connected qubits on the device.Table I reports

Statistical bootstrap data
In this section we provide the complete dataset produced by applying the bootstrap procedure to the measured outputs of each quantum simulation performed on the ibmq mumbai processor in this study.In particular Fig. 8 shows the bootstrap data for nearest neighbor Ising interactions R = 1, different noise scales s = 1, 1.2, 1.4, 2.6 and different stroboscopic times n = 0, 5, 8, 10, 12, while Fig. 9 shows the same the bootstrap data for next to nearest neighbor Ising interactions R = 2, different noise scales s = 1, 1.2, 1.4, 2.6 and different stroboscopic times n = 0, 5, 8.

FIG. 1 .
FIG. 1.Quantum circuit implementation of Floquet dynamics with varying interaction range.(a) Topology of the ibmq mumbai quantum processor.Black links represents the physical connections among the qubits on the quantum hardware, blue and red links represent the nearest neighbors R = 1 and next-to-nearest neighbor R = 2 Ising interactions we effectively implemented among physically unconnected qubits during our quantum simulation.(b) Quantum circuit implementing a single Floquet step for a kicked Ising model with R range interactions.(c) Quantum gate implementation of nearest-neighbor Ising interaction among qubit j and j + 1.(d) Quantum gate implementation of next-to-nearest neighbor Ising interaction among qubit j and j + 2. (e) Quantum gate implementation of r-neighbor Ising interaction among qubit j and j + r.

FIG. 2 .
FIG.2.Modulus of the magnetization |mz| as a function of the stroboscopic times n for nearest neighbor R = 1 (blue points) and next-to-nearest neighbor R = 2 (red points) Ising interactions.Triangles represent the raw experimental data measured on our quantum simulation of the ibmq mumbai quantum processor, which involves N = 18 qubits undergoing a kicked Ising dynamics with kick angle ϕ = π/2 + ϵ with ϵ = 0.2.Square points and the corresponding error bars represent the estimators for the average magnetization and its statistical error E(mz) ± 2σ(mz), obtained through statistical bootstrapping.Dashed lines represent the best fit of the data with an exponential decay.

1 FIG. 3 .
FIG.3.Decay rate of magnetization as a function of noise scale s for R = 1 (blue points) and R = 2 (red points).Error bars represent two standard deviations σ(mz) estimated through statistical bootstrapping.Dashed lines indicate the best linear fit obtained using weighted least squares regression.Empty points are excluded from the fitting data.

FIG. 4 .
FIG. 4. Quantum circuit optimization techniques using circuit identities.(a) Cancellation of CNOT gates in adjacent Vj,j+1 and SWAP gates.(b) Rearrangement of the quantum circuit to implement Ising interactions with r = 1, . . ., R ranges while maximizing the number of adjacent Vj,j+1 and SWAP gates.

8 FIG. 5 .
FIG. 5. Bootstrap distribution of averages obtained by resampling M = 1000 times the measured data of three quantum simulations with range R = 1, noise scale s = 1.4,and different numbers of Floquet steps n = 0, 5, 8.The number of bins considered in each histogram is 100.

FIG. 6 .FIG. 7 .
FIG.6.Topology of the ibmq mumbai quantum processor with numbered qubits.Grey qubits were not involved in our quantum simulations.

TABLE I .
Calibration data for each qubit of the processor of the ibmq mumbai quantum processor at the time our experiments were performed.