Error-Mitigated Simulation of Quantum Many-Body Scars on Quantum Computers with Pulse-Level Control

Quantum many-body scars are an intriguing dynamical regime in which quantum systems exhibit coherent dynamics and long-range correlations when prepared in certain initial states. We use this combination of coherence and many-body correlations to benchmark the performance of present-day quantum computing devices by using them to simulate the dynamics of an antiferromagnetic initial state in mixed-field Ising chains of up to 19 sites. In addition to calculating the dynamics of local observables, we also calculate the Loschmidt echo and a nontrivial connected correlation function that witnesses long-range many-body correlations in the scarred dynamics. We find coherent dynamics to persist over up to 40 Trotter steps even in the presence of various sources of error. To obtain these results, we leverage a variety of error mitigation techniques including noise tailoring, zero-noise extrapolation, dynamical decoupling, and physically motivated postselection of measurement results. Crucially, we also find that using pulse-level control to implement the Ising interaction yields a substantial improvement over the standard CNOT-based compilation of this interaction. Our results demonstrate the power of error mitigation techniques and pulse-level control to probe many-body coherence and correlation effects on present-day quantum hardware.

Quantum simulation is one of the most natural problems in which to expect a quantum computer to outperform a classical one.In particular, simulating the time evolution of a quantum many-body system with local interactions on a quantum computer requires a number of local quantum gates scaling polynomially with the number N of qubits and linearly with the total simulation time T [1].Classical simulation techniques, in contrast, require resources scaling exponentially with N , restricting most studies of quantum many-body dynamics to small systems and/or early times.There is thus a relatively clear path to quantum advantage for the quantum simulation problem, assuming the availability of quantum hardware that can simulate the dynamics of a quantum many-body system with sufficiently large N and T .
The present generation of quantum hardware operates in the so-called noisy intermediate-scale quantum (NISQ) regime [2].NISQ devices have enough qubits (∼ 10 1 -10 3 ) to potentially evade classical simulability, but coherence times and gate fidelities that are too small to simulate dynamics beyond the early-time regime, where tensor-network methods [3] are often still applicable.As hardware continues to improve, a major challenge of the NISQ era is to determine strategies to maximize the utility of such devices despite their imperfections.One approach to this problem is to use variational quantum algorithms [4,5], which reduce the required circuit depth for tasks including time evolution [6][7][8][9][10][11][12][13][14] at the cost of requiring many circuit evaluations.An alternative approach is to apply a growing toolbox of hardware-and noise-aware quantum error mitigation techniques to a relatively simple quantum simulation algorithm, e.g. the first-order Trotter method of Ref. [1].Error mitigation strategies can be applied at the hardware level, e.g. by applying dynamical decoupling pulses to idle qubits [15][16][17] or designing optimized pulse sequences to reduce gate execution times [18,19].Other strategies, such as randomized compilation [6,20] and zero noise extrapolation [6,[21][22][23], run additional circuits that are logically equivalent to the target circuit and perform postprocessing of the data to estimate the hypothetical noiseless result.Recently, Ref. [19] showed that combining these error-mitigation techniques yields first-order Trotter simulations of quantum dynamics whose accuracy as measured by local observables is competitive with tensor-network methods, at least for low bond dimensions.
In this work, we test the ability of present-day quantum hardware to simulate nontrivial many-body dynamics and correlation effects in a prototypical interacting quantum system: the one-dimensional mixed-field Ising model (MFIM).This model and its variants have been simulated on NISQ hardware in several recent works [24][25][26].Our goal here is to model a particular physical phenomenon known as quantum many-body scars (QMBS) [27][28][29][30] (see Refs. [31][32][33] for reviews).This phenomenon was observed experimentally in an analog quantum simulation of the MFIM using Rydberg atoms in optical tweezers [34], where coherent oscillations of the local Pauli expectation values ⟨Z i (t)⟩ were observed after a quantum quench from the Néel state |Z 2 ⟩ = |010 . ..⟩ or its Z 2 conjugate |Z ′ 2 ⟩ = |101 . ..⟩.This came as a surprise, since the MFIM is known to be nonintegrable for generic values of the transverse and longitudinal fields and the Néel states |Z 2 ⟩ and |Z ′ 2 ⟩ have a finite energy density relative to the ground state of the model.In such a case, reasoning based on the eigenstate thermalization hypothesis (ETH) [35,36] leads to the expectation that the dynamics from these initial states should rapidly deco-arXiv:2203.08291v2[quant-ph] 1 Nov 2023 here [37,38].Intriguingly, the oscillatory dynamics also encode coherent oscillations in space: Ref. [39] argued based on finite-size numerics that the dynamics from the Néel state exhibits long-range connected unequal-time correlations at wavenumber π.This finding also defies intuition based on the ETH and suggests that the oscillatory dynamics observed in the experiment [34] are fundamentally many-body in nature, and cannot be explained by the precession of free spins.QMBS are known to occur in a variety of other models [28,[40][41][42][43][44][45][46][47][48][49], and have been observed in several analog quantum simulation experiments [34,50,51], but we focus here on a digital quantum simulation approach.
Motivated by these experimental and theoretical results, we use IBM quantum processing units (QPUs) to perform a first-order Trotter simulation of the MFIM in the regime with QMBS for up to 39 Trotter steps on systems of up to 19 qubits.One goal of the study is to use the coherent dynamics in the scarred regime to benchmark the performance of these QPUs; another is to use QPUs to verify the presence of nontrivial connected unequal-time correlations in the dynamics of the Néel state.Such correlations, which demonstrate the inextricably many-body nature of the oscillatory dynamics, are challenging to measure in analog quantum simulators and have yet to be probed experimentally.Our QPU results demonstrate that they can be accessed using digital quantum simulation.
To optimize the performance of the devices, we make use of a variety of techniques.First, we implement quantum simulation of the Ising interaction using a scaled cross-resonance pulse and show that this implementation outperforms the naive compilation of the Ising evolution operator using two controlled-NOT (CNOT) gates.Second, we apply an arsenal of error mitigation techniques, including zero-noise extrapolation, Pauli twirling, dynamical decoupling, readout error mitigation, and, where appropriate, physically motivated postselection of computational-basis measurement outcomes.Many of these methods were applied to simulate the transversefield Ising model on the heavy hexagon lattice for at most 20 Trotter steps in Ref. [19]; we will comment on areas where our implementation differs from theirs, the most important being our use of postselection and a new approach to Pauli twirling of non-Clifford gates.
We combine these techniques to calculate spatially averaged local observables, as well as more sensitive probes of the dynamics including the Loschmidt echo and a nontrivial finite-wavenumber connected correlation function.We find that combining pulse level control and error mitigation techniques extends by roughly a factor of two the timescales over which nontrivial oscillatory dynamics can be observed.Taken together, our results demonstrate that these nontrivial many-body effects can be probed, at least in the early time regime, on present-day quantum hardware.
The remainder of the paper is organized as follows.In Sec.I we define the mixed-field Ising model and the related observables that we will calculate on the QPU.In Sec.II, we describe the scaled cross-resonance pulse ("scaled-R ZX " for short) implementation of the Ising interaction.We present data benchmarking its performance against the more standard implementation of the interaction using two CNOT gates.In Sec.III, we present data for simulations of a 12-and 19-site chain using both the two-CNOT and scaled-R ZX implementations, and use this to motivate a brief discussion of the error mitigation techniques we use.We show our main results in Sec.IV, which includes error-mitigated results for the staggered magnetization, the Loschmidt echo, and the finite-wavenumber connected unequal-time correlator.Finally, conclusions and outlook are discussed in Sec.V.

A. Model
In this work we simulate the dynamics of a chain of L spin-1/2 degrees of freedom generated by the Hamiltonian where n i = I−Zi 2 , and where Z i and X i are Pauli operators on site i.This model arises when studying chains of trapped Rydberg atoms with rapidly decaying van der Waals interactions, where the operator n i is interpreted as an occupation number for the local atomic Rydberg state and the operator X i induces transitions between the ground and Rydberg state.In this paper, we will label computational basis (CB) states using the eigenstates |0⟩ i and |1⟩ i of the Z i operator for which n i |0⟩ i = 0 and n i |1⟩ i = |1⟩ i .Rewriting Eq. (1.1) in terms of Pauli matrices, we obtain The Hamiltonian (1.2) is an Ising model with transverse and longitudinal fields (the MFIM), but in which the strength of the longitudinal field is tied to the interaction strength V as a consequence of the model's origin in Eq. (1.1).Note that the model is written with open boundary conditions, and that the longitudinal field strength is reduced by a factor of two on the first and last sites of the chain (i = 1, L).
The model (1.2) is nonintegrable for generic values of V and Ω. QMBS emerge in the limit V ≫ Ω, which is known in the Rydberg-atom literature as the "Rydberg blockade" regime [52,53].In this limit, computational basis states with different eigenvalues of The remaining angles are given by θ X i = 2Ω∆t, θ Z i = −4V ∆t, and decouple into sectors separated by an energy scale ∼ V .When the system is initialized in a CB state with ⟨ L−1 i=1 n i n i+1 ⟩ = 0, such as the Néel states |Z 2 ⟩ and |Z ′ 2 ⟩, the probability of finding nearest-neighbor sites in the configuration |11⟩ is heavily suppressed.The subspace of the full Hilbert space in which no two nearestneighbor sites are in the state |11⟩ is known as the Fibonacci Hilbert space, as the number of such states scales as φ L , where φ is the golden ratio.An effective model for the system in this limit is known as the "PXP model," which arises from the projection of H X into the Fibonacci Hilbert space [27,34].
To simulate the system's dynamics under the Hamiltonian (1.2), we employ a first-order Trotter decomposition of the unitary evolution operator over a time ∆t: ZiZi+1/2 is shown in Fig. 1.Evolution over a time T = n∆t is obtained by applying the circuit (1.3) n times.

B. Observables
To probe the dynamics of the model (1.2) in the QMBS regime V ≫ Ω, we use the QPU to measure three dynamical properties.First, to characterize the oscillations, we measure the dynamics of the expectation value of the staggered magnetization operator, This operator takes its extremal values ⟨Z π ⟩ = ±L when the system is prepared in the Néel state |Z ′ 2 ⟩ or |Z 2 ⟩, respectively.Exact simulations of a quantum quench from the |Z 2 ⟩ state show a weakly damped coherent oscillation of ⟨Z π (t)⟩ with frequency ω ≃ 1.33 Ω [27], indicating that the system's state is periodically cycling between the two Néel states.In contrast, the expectation based on the ETH is that ⟨Z π (t)⟩ would decay rapidly, on a timescale ∼ 1/Ω, to its thermal value of 0. These coherent oscillations arise due to the presence of a tower of eigenstates with roughly equal energy spacings in the many-body spectrum [27].These "scar states" have high overlap with the Néel states, so preparing the system in one of these initial states projects the ensuing dynamics strongly onto this set of special eigenstates.⟨Z π (t)⟩ is calculated on the QPU by performing Trotter evolution out to time t and measuring the state in the computational basis.
Second, we measure the Loschmidt echo, where |ψ(0)⟩ is taken to be |Z 2 ⟩.When the system is prepared in this initial state, the scarred eigenstates give rise to sharp periodic revivals of the Loschmidt echo to a value of order 1, with period matching that of the oscillations in ⟨Z π (t)⟩ [29].This behavior is highly atypical-L(t) is expected to decay to zero exponentially fast in quantum quenches of nonintegrable models from typical high-energy-density initial states [54,55].Thus, the Loschmidt echo is a much more sensitive quantity than ⟨Z π (t)⟩, which is built from expectation values of local observables.To measure the Loschmidt echo on the QPU, we perform Trotter evolution of the |Z 2 ⟩ state and measure in the CB to determine the probability to be in the state |Z 2 ⟩ after a time t.Third, we measure the correlation function where with  These long-range correlations arise due to the presence of off-diagonal long-range order in the scarred eigenstates [39,56].This nontrivial correlator has yet to be measured experimentally.Measuring it on a QPU is challenging, but achievable using the ancilla-free protocol of Ref. [57], which we describe further in Sec.IV B.

II. PULSE-LEVEL IMPLEMENTATION OF THE ISING INTERACTION
Implementing the Trotter circuit in Eq. (1.3) requires realizing the two qubit gate R ZZ (θ) on the device.One approach to solving this problem is to decompose R ZZ (θ) into a basis gate set, the most standard of which includes CNOT and arbitrary one-qubit rotations.In this basis, R ZZ (θ) can be realized by applying two CNOT gates on either side of an R Z (θ) gate on the target qubit [58,59], as shown in Fig. 2(a).
An alternative approach is to leverage knowledge of the basic set of pulses used to generate two-qubit gates at the hardware level.On IBM QPUs, well-calibrated CNOT gates are built by adding single qubit gates before and after the R ZX (π/2) gate [60], which is defined via R ZX (θ) = e −iθZcXt/2 , where c and t denote the control and target qubits, respectively.These R ZX (π/2) gates are realized using an echoed cross resonance pulse sequence described in further detail in Appendix A and in Ref. [60].The ability to implement an R ZX (θ) gate opens another route to realize the R ZZ (θ) gate simply by dressing the R ZX (θ) gate with R Y (±π/2) gates on the target qubit, see Fig. 2(b).To realize the R ZX (θ) gate with arbitrary rotation angle on the QPU, we scale the pulse amplitudes and durations used to generate the R ZX (±π/2) gate in the manner described in Refs.[18,61] and summarized in Appendix A. We therefore refer to this pulse-level implementation of the R ZZ gate as the "scaled-R ZX " im-plementation.The pulse sequences are programmed using Qiskit pulse [60,62]; examples of pulse schedules used in our simulations are shown in Appendix A (Fig. 8).
Since the scaled-R ZX approach uses fewer crossresonance pulses than the two-CNOT implementation of R ZZ , we expect the former to yield pulse schedules with shorter overall duration than the latter.The duration of the pulse schedule that realizes the two-CNOT implementation of R ZZ (θ) is independent of the rotation angle θ, since the R Z (θ) gate is simply realized as a phase shift on the pulse schedule [60].In contrast, in the scaled-R ZX implementation of R ZZ , the cross-resonance pulse duration depends roughly linearly on θ for θ above a certain threshold and is constant below that threshold (see Appendix A).The pulse durations in ns of the 2-CNOT and scaled-R ZX implementations of R ZZ (θ) on the IBM QPU Casablanca (ibmq casablanca) are shown as a function of θ in Fig. 2(c).Despite the θ-dependence of the scaled-R ZX pulse duration, there is a wide range of interaction strengths V and Trotter time steps ∆t giving an angle θ = 2V ∆t such that the scaled-R ZX implementation has the shorter pulse duration of the two methods.
The total duration of the pulse schedule that realizes a given quantum gate is positively correlated with the gate's error rate.We therefore expect that the scaled-R ZX implementation of R ZZ should have a lower error rate than that of the two-CNOT implementation of the same gate.To compare the error rates of the R ZZ gates rates realized using the two approaches, we measure the fidelity of R ZZ (θ) at 12 different angles from θ = 0.2 to θ = 2.4 using quantum process tomography (QPT) [63,64], which is built into IBM's Ignis module, with state preparation basis {|0⟩ , |1⟩ , |X + ⟩ , |Y + ⟩} and measurement basis {X, Y, Z} for each qubit.In order to obtain an estimate of the gate error that is decoupled from state preparation and measurement (SPAM) errors, we use a simple scheme relying on gate folding.Letting G = R ZZ (θ), we consider the sequence of logically equivalent gates G, GG † G, and GG † GG † G, which correspond to a "scale factor" of λ = 1, 3, and 5, respectively.For each λ, we use QPT to estimate the average gate fidelity, which includes SPAM errors.This fidelity decreases with λ because gate folding increases the gate noise in the circuit.We then fit the resulting data points to a linear model F 0 − ϵλ, which is justified under the assumptions that ϵ is small and that G and G † have identical error rates.The slope ϵ is an estimate of the error rate that is free of SPAM errors, since these errors do not scale with λ.To obtain the results plotted in Fig. 2(d), we ran the above procedure on the IBM QPU Casablanca using qubits q1 and q3, with 1024 shots for each measurement and using complete readout error mitigation [65] (see also Appendix A).We also repeated QPT four times for each scale factor to collect better statistics; the linear fit to extract ϵ was performed over the full data set.The error bars for each θ value represent the standard deviation of the slope calculated from the covariance matrix of the linear fit for that θ.
Fig. 2(d) shows that the fidelity of the scaled-R ZX implementation of R ZZ (θ) decreases with increasing θ, while the fidelity of the two-CNOT implementation remains almost constant.Because the experiment was implemented across two days, there are also fluctuations in the fidelity due to calibration drifts (see also Ref. [18]).At the smallest value of θ = 0.2, our fidelity results indicate that the scaled-R ZX approach realizes an ∼ 80% reduction in the error rate of the R ZZ (θ) gate relative to the two-CNOT implementation.This decreases to a ∼ 50% error reduction at the largest value of θ = 2.5.This is consistent with the pulse duration results plotted in Fig. 2(c), which show that the pulse durations of the two-CNOT and scaled-R ZX implementations approach one another with increasing θ.
We now comment on our choice of Hamiltonian parameters V and Ω for simulating the system's dynamics in the regime with QMBS.As mentioned in Sec.I, we need V ≫ Ω in order to be in the QMBS regime.However, according to Fig. 2(c), reducing V yields a shorter pulse duration for the scaled-R ZX implementation of R ZZ , thereby reducing the gate error.Moreover, Ω determines the frequency of the oscillations that are characteristic of QMBS, so choosing a larger Ω is desirable in order to manifest more oscillation periods within a fixed time window.The choice of ∆t is essential as well, since the Trotter error scales to leading order as V Ω(∆t) 2 .After trying many sets of parameters, we settled on V = 1, Ω = 0.24, and ∆t = 1 as optimal parameters to simulate the system's dynamics in the QMBS regime.These parameters correspond to an R ZZ rotation angle θ = 2.0, where the data in Fig. 2(d) indicate that the scaled-R ZX approach yields a ∼ 57% error reduction relative to the two-CNOT implementation.Although the sizable value of ∆t incurs substantial Trotter error (see Appendix A for a comparison between Trotter and exact dynamics), the Trotter circuit with ∆t = 1 nevertheless exhibits pronounced co- herent oscillations with period 2π/(1.33Ω)≃ 19.68 for this choice of parameters.For this choice of parameters, the dimensionless quantity V t is simply the number of Trotter steps.
We note in passing that our Trotter circuit (1.3) could alternatively be viewed as a Floquet circuit due to the relatively large value of ∆t.That is, we can view this circuit as simulating not the dynamics under the timeindependent Hamiltonian in Eq. (1.2), but rather the time-dependent Hamiltonian where f ± (t) = [1 ± sgn(sin 2π ∆t t)]/2, whose evolution operator over a time ∆t is precisely the right-hand side of Eq. (1.3).QMBS have been studied in several Floquet variants of the mixed-field Ising and PXP models, see e.g.Refs.[66][67][68][69][70], so it is by now well established that they can exist in this periodically driven setting.Our results demonstrate their existence in the model (2.1).

III. ERROR MITIGATION TECHNIQUES
Based on the discussion in the previous section, we expect that the scaled-R ZX implementation of the R ZZ gate will outperform the two-CNOT implementation when performing Trotter evolution of the Néel state |Z 2 ⟩ under the Hamiltonian (1.2).We test this hypothesis by running Trotter simulations using both approaches on the IBM QPUs Guadalupe (ibmq guadalupe) and Toronto (ibmq toronto) for chains with L = 12 and 19 sites, respectively.For each R ZZ implementation, we execute 39 Trotter steps using the parameters Ω = 0.24, V = 1, and ∆t = 1.At each time t, we measure ⟨Z i (t)⟩ for all sites i using 8192 shots.We repeat the time evolution procedure 20 times and average the results over these trials [71].To quantify the accumulation of error during the simulation, we define where ⟨Z i (t)⟩ SV is the result obtained from Trotter evolution on a noiseless statevector simulator and ⟨Z i (t)⟩ QPU is the result from the QPU.The results of these simulations are shown in Fig. 3. Figs.3(a) and (b) show the dynamics of the staggered magnetization density ⟨Z π (t)⟩ /L [see Eq. (1.4)] for the 12-and 19-qubit simulations, respectively.Note that we did not apply any error mitigation techniques here in order to directly compare the two implementations of the R ZZ gate.For both simulations, both implementations clearly fail to reproduce the expected oscillations.The simulation using the scaled-R ZX approach does slightly outperform the two-CNOT approach-in particular, one weak oscillation over a scar period V t ≃ 20 is barely observable for the scaled-R ZX data.However, both approaches perform poorly compared to the ideal Trotter results from the statevector simulator.For the 12-qubit system, the results begin to deviate strongly from the ideal Trotter curve after V t = 8.For the 19-qubit system, the calculations disagree markedly even at early times.Figs.3(c) and (d) show that the accumulated error grows dramatically before V t = 25, after which it appears to saturate.The saturated value of the error is nearly the same for the 12-and 19-qubit calculations.
There are many error sources that contribute to these results.One source is quantum thermal relaxation, whose effects are quantified by the qubit relaxation time T 1 and the qubit dephasing time T 2 .On IBM QPUs, these timescales are roughly T 1,2 ∼ 100 µs, which limits the total circuit depth that can be executed on the device.Another source of error comes from imperfect operation of the physical gates.On IBM QPUs, single-qubit gates have error rates of ∼ 10 −4 to 10 −3 , while CNOT gates (which involve the use of R ZX gates as described in Sec.II) have error rates of ∼ 5×10 −3 to 1.5×10 −2 .Thus, two-qubit gates provide the dominant source of gate er-ror.Finally, readout errors can occur in which the device misreports the state of a qubit as |0⟩ when it is actually |1⟩ and vice versa.Readout error rates can range from ∼ 2 × 10 −2 to 6 × 10 −2 , which is even larger than the CNOT error rate.However, readout only occurs once per circuit, whereas each circuit can use many two-qubit gates.(Note that readout error likely accounts for the majority of the difference between the QPU and ideal Trotter results for both the 12-and 19-qubit simulations at time t = 0.) The average two-qubit error rates were nearly identical for ibmq guadalupe and ibmq toronto (1.8 × 10 −2 and 1.6 × 10 −2 , respectively) when the data shown in Fig. 3 were obtained.However the average readout error rate for ibmq guadalupe was 3.6 × 10 −2 , while for ibmq toronto it was 5.4 × 10 −2 .In fact, the highest readout error rate among the qubits we used on ibmq guadalupe was 9 × 10 −2 , while on ibmq toronto it was 2.16 × 10 −1 .Therefore, we believe that the poorer agreement with the ideal Trotter simulation at early times that was observed for the results obtained for the 19-qubit chain on ibmq toronto is due to the increased readout error rate on that device at the time the experiments were performed.
To reduce the impact of these various error sources, we implemented an arsenal of error mitigation techniques.The simplest of these techniques is readout error mitigation, which is built into Qiskit Ignis.We also implement dynamical decoupling using an X π − X −π pulse sequence to reduce decoherence errors [15][16][17].To reduce the effect of stochastic gate error in a circuit execution, we use the Mitiq package [72] to implement zero-noise extrapolation (ZNE) using random gate folding.This method scales the gate noise by performing gate folding G → GG † G on randomly chosen two-qubit gates throughout the circuit.This results in a noise scale factor λ that can be noninteger, in contrast to the simpler global gate folding procedure described in Sec.II.ZNE is best justified in the case where gate errors result in a stochastic quantum channel.For this reason, we also implement Pauli twirling [6,21,73], in which two-qubit gates are dressed with random Pauli gates chosen so as not to affect the outcome of a circuit execution in the zero-noise limit.Averaging the results of many of these random circuit instances reduces the gate noise to a stochastic form.The details of our implementation of these techniques are explained in more depth in Appendix A.
Finally, to enhance the signatures of the characteristic oscillatory dynamics, we implement postselection of measurement data to exclude CB measurement outcomes in which nearest-neighbor sites were measured to be in the state |11⟩.As discussed in Sec.I, the probability of such configurations appearing when evolving the Néel state under the Hamiltonian (1.2) is heavily suppressed when V /Ω is large.In Appendix A, we show strong numerical evidence that this is the case for the ratio V /Ω ≃ 4.17 used in this work.There we also show, however, that the Trotter error due to the large step size ∆t = 1 substantially increases the probability of generating these "forbidden" configurations.Thus, postselection mitigates the effect of both Trotter and gate error.In this paper we will always compare QPU results using postselected data with exact results in which the Trotter dynamics generated by Eq. (1.3) are projected into the Fibonacci Hilbert space before calculating the observable of interest.
We note that the set of error mitigation techniques we use for this work (including the scaled-R ZX implementation of the R ZZ gate) is similar to that used in Ref. [19].We highlight here a few important differences upon which we expand in Appendix A. First, Ref. [19] takes a different approach to ZNE wherein the noise is scaled by scaling the duration and amplitude of the crossresonance pulses.In contrast, our ZNE scheme treats the scaled-R ZX pulse schedule for the R ZZ (θ) gate as a custom gate which is folded in the same way as other gates, including CNOT.Ref. [19] also takes a different approach to Pauli twirling of the R ZZ (θ) gate, which is a non-Clifford gate for generic θ.In particular, Ref. [19] performs Pauli twirling using only random Pauli operators from the set {II, XX, Y Y, ZZ}.Our Pauli twirling method, described in Appendix A, uses the full set of two-qubit Pauli operators to perform the twirling, and we prove that it results in a stochastic noise channel.

A. Zπ and Loschmidt Echo
We now test the degree to which the error mitigation strategies outlined in Sec.III improve the results shown in Fig. 3.The results of fully error-mitigated calculations of ⟨Z π (t)⟩ /L and D(t) on ibmq guadalupe and ibmq toronto are shown in Fig. 4. For these simulations, we performed ZNE with random gate folding scale factors λ ∈ {1, 1.5, 2.0} and perform Pauli twirling with 10 random circuit instances.As in Fig. 3, we evolve for 39 Trotter steps and use 8192 shots per circuit execution.Unlike Fig. 3, we use postselected QPU data and compare with the Fibonacci-projected Trotter evolution (black).With the extra circuits needed to perform ZNE and Pauli twirling, the total number of circuits run on each device is now 40 × 10 × 3 = 1200 (including t = 0).ZNE is performed with a linear extrapolation to λ = 0 for each Trotter step, with 10 data points for each scale factor.Each data point in Figs.4(a In particular, for the 12-qubit calculation performed on ibmq guadalupe, the scaled-R ZX approach yields a staggered magnetization density ⟨Z π (t)⟩ /L that is in good agreement with the Trotter simulation until roughly V t = 15.This is roughly a twofold improvement relative to the case without error mitigation.For both the 12-and 19-qubit calculations, the scaled-R ZX approach yields visible oscillations up to the final time V t = 39, which covers about two oscillation periods.In contrast, even with error mitigation the two-CNOT implementation of the Trotter circuit yields visible oscillations only over one period.In Figs. 4 (c) and (d), we see that the accumulated error D(t) is substantially reduced as compared to the unmitigated results shown in Fig. 3.For the scaled-R ZX implementation of the Trotter circuit, error mitigation leads to a reduction of the total accumulated error D(t = 39V −1 ) by 56% (29%) for the 12-qubit (19-qubit) calculation.In contrast, error mitigation of the two-CNOT implementation of the Trotter circuit results in a 37% (22%) error reduction for the 12-qubit (19-qubit) calculation.Compared to results where only postselection is applied (see Appendix A), the results in Fig. 4 show an error reduction of 37% (29%) for the scaled-R ZX implementation and Vt FIG. 5.
Loschmidt echo L(t) of the initial state |Z2⟩ for a chain of 12 qubits calculated on ibmq guadalupe using the same error mitigation techniques described in Fig. 4 and the text.Data obtained from simulations using the two-CNOT implementation of the RZZ gate (red) barely show any tendency toward a revival around V t = 20, whereas data obtained using the scaled-RZX implementation (blue) show a more pronounced revival.Both revivals are nowhere near as pronounced as the one obtained from the Fibonacci-projected ideal Trotter simulation (black), further indicating the effect of errors on the simulation despite the error-mitigation measures used.
23% (20%) for the two-CNOT implementation for the 12-qubit (19-qubit) calculation.As in Fig. 3 (c) and (d), the rate of error accumulation for the scaled-R ZX implementation is always less than it is for the two-CNOT implementation.Interestingly, with error mitigation the final accumulated error is larger for the 19-qubit calculation than for the 12-qubit calculation [cf.Fig. 3].
To further test the extent to which error mitigation improves the accuracy of the results obtained from the QPUs, we calculate the Loschmidt echo L(t) [see Eq. (1.5)] for a 12-qubit system using the same data set from which the results of Fig. 4 (a) and (c) were obtained.The results are shown in Fig. 5.The scaled-R ZX implementation of the Trotter circuit shows a faint but noticeable revival near the first oscillation period V t = 2π/(1.33Ω)≃ 19.68, while the two-CNOT implementation shows hardly any revival at all.Note that the fact that the Loschmidt echo is measured to take a finite value on the device after ∼ 20 Trotter steps is remarkable given the exponential sensitivity of the Loschmidt echo to changes in the state |ψ(t)⟩.The fact that the scaled-R ZX implementation of the Trotter circuit yields a substantially enhanced revival is further evidence of the performance advantage offered by that approach.
To enhance the Loschmidt echo signal, we also consider the effect of counting shots in which the measure-ment outcome differs from the Néel state |Z 2 ⟩ by a single bit flip.This makes the metric L(t) more forgiving by counting instances where the system almost returns to the initial state.The results, shown in Fig. 5, demonstrate that this protocol indeed boosts the amplitude of the first revival in L(t).We observe greater enhancement of the first revival for the scaled-R ZX implementation of R ZZ than for the two-CNOT implementation, consistent with our other results.In both cases, the signal enhancement is localized in time near the first revival time but becomes more diffuse at later times.

B. Connected Correlation Function
We now discuss how we measure the nontrivial correlation function C Y (t) [see Eq. (1.6)] on a QPU.The correlation function C Y (t) is of the form ⟨O(t)O(0)⟩ where O(t) = e iHt Oe −iHt is a Hermitian operator that can be expanded in the basis of Pauli strings.One way to measure such a correlator on a quantum computer is to use a so-called indirect measurement technique based on the Hadamard test [74,75].This approach uses an ancilla qubit and relies on the ability to apply Pauli gates controlled by the ancilla.For an operator like O = Y π , which is a sum of many local Pauli strings, this means that the ancilla qubit must be able to couple to all qubits in the chain.Achieving this on superconducting qubit QPUs with nearest-neighbor connectivity requires substantial gate overhead, making this approach somewhat impractical for our purposes.We therefore opt instead for a direct measurement approach, which avoids the use of ancilla qubits at the cost of running more circuits with fewer gates.We now describe this approach, which is based on the proposal of Ref. [57].
C Y (t) can be calculated as a sum of many local correlators, i.e., ) where the operators (P Y P ) i are defined in Eq. (1.6c).This can be further simplified using our knowledge that the initial state is |Z It remains to evaluate the local correlators In Ref. [57], it is shown that these can be calculated as with Here, UTrotter denotes the Trotter circuit.In the beginning of each circuit, a different initial state is prepared.At the end of each circuit, the expectation value of (P Y P )j is measured.
In the above expressions, U (t) is the evolution operator out to time t, which is approximated on the QPU by the Trotter circuit.Note that both Eqs.(4.3b) and (4.3c) can be formulated as the expectation value of (P Y P ) j in a particular time-evolved state.Quantum circuits to evaluate these expectation values are shown in Fig. 6.
To prepare the initial state needed to evaluate Eq. (4.3b) on the device, we act on the ith qubit with either an identity or an X gate (depending on the choice of + or −, respectively), followed by a Hadamard gate and an S gate.(Note that i is even, so the initial state of the ith qubit is always |1⟩ in the |Z 2 ⟩ state).
We have calculated C Y (t) for chains of L = 5 and L = 12 sites on ibmq casablanca and ibmq guadalupe, respectively, for both the QMBS regime (V = ∆t = 1 and Ω = 0.24) and the chaotic regime, where we use parameters V = 1, Ω = 2 and ∆t = 0.16.The calculation of C Y (t) on the QPU proceeds as follows.For each i and j, we need to evaluate the four circuits shown in Fig. 6 to calculate (P Y P ) j (t) M Y i =±1 and ⟨(P Y P ) j (t)⟩ ±Yi .
Note that [(P Y P ) j , (P Y P ) j+2 ] = 0, so we can measure (P Y P ) j for all even and all odd j in one shot-in fact, this is one of the advantages of the direct measurement scheme.The total number of circuits needed to calculate C Y (t) is therefore 4×2×⌊L/2⌋, where ⌊•⌋ denotes the floor function.For each of these circuits, we employ ZNE with random gate folding for scale factors λ ∈ {1.0, 1.5, 2.0}.For the L = 5 calculations, we employ Pauli twirling with 8 random circuit instances in both the QMBS and chaotic regimes for each scale factor used in ZNE.For the L = 12 calculations, we used 10 random circuit instances for the QMBS regime and 9 for the chaotic regime.For all cases, we evolve the system for 30 Trotter steps and measure the system with 8192 shots at each step [76].For the 5-qubit system, we evaluate a total of 11520 circuits to measure C Y (t); for the 12-qubit system, we evaluate 43200 circuits for the QMBS case and 38880 circuits for the chaotic case.We use the scaled-R ZX implementation of the Trotter circuit, since the results of Sec.IV A indicate that this approach outperforms the two-CNOT implementation.
The results of these calculations are shown in Fig. 7(a,c) for the QMBS regime and in panels (b,d) for the chaotic regime.We show results for |C Y (t)| for ease of visualization, as Eqs.(4.3) demonstrate that C Y (t) is generically complex.[Error bars on individual data points are calculated as described at the beginning of Sec.IV A.] In Fig. 7(a), we see that for L = 5 the calculation of C Y (t) exhibits good quantitative agreement with the ideal Trotter calculation out to V t ≃ 15, and qualitative agreement is maintained throughout the full time window until V t = 30.In particular, oscillations with the expected period π/(1.33Ω)≃ 9.84 are visible throughout the time evolution window (recall that we plot the absolute value of the correlation function).In contrast, in the chaotic regime for L = 5 shown in Fig. 7(b), we find that |C Y (t)| exhibits one approximate revival followed by a rapid decay and incoherent dynamics after time V t = 2.Note that, in the chaotic regime, Trotter dynamics in- volves a smaller R ZZ rotation angle than in the QMBS regime, resulting in shorter cross resonance pulse durations and higher R ZZ gate fidelities.Consequently, the dynamics exhibit excellent quantitative agreement with the ideal Trotter simulation for approximately 22 Trotter steps.The results of the L = 12 calculation are shown in Fig. 7(c) and (d) for the QMBS and chaotic regimes, respectively.While the results for the chaotic regime remain in very good agreement with the ideal Trotter simulation due to the smaller R ZZ rotation angle discussed above, the results in the QMBS regime begin to differ from the ideal Trotter results around V t = 9.This is likely due to the fact that the calculation of C Y (t) involves summing O(L 2 ) terms, each of which suffers from gate and readout errors and is the result of a separate zero noise extrapolation.Despite the more drastic accumulation of error for L = 12, underdamped oscillations close to the correct frequency are clearly visible throughout the simulation time window.These results provide a clear demonstration that the coherence and long-range many-body correlations present in the dynamics of the Néel state in the QMBS regime can be probed on current quantum devices.

V. CONCLUSIONS AND OUTLOOK
In this work, we have used pulse-level control and a variety of quantum error mitigation techniques to simulate the dynamics of a spin chain with QMBS on IBM QPUs for chains of up to 19 qubits.QMBS constitute an intriguing quantum dynamical regime characterized by nontrivial many-body coherence and long-range correlations.We probed this physics by measuring the dynamics of three quantities: the staggered magnetization ⟨Z π (t)⟩, the Loschmidt echo L(t), and the connected unequal-time correlation function C Y (t).We found that ⟨Z π (t)⟩ and C Y (t) exhibit reasonable quantitative agreement with the ideal Trotter simulation at early times, and visible oscillations with the correct frequency over 39 Trotter steps.In contrast, the Loschmidt echo L(t) exhibits only the faintest of revivals, indicating the substantial impact of various noise sources including thermal relaxation and gate and readout errors.Nevertheless, the qualitative features of QMBS are pronounced on time scales beyond which L(t) decays to zero, indicating the presence of coherent many-body dynamics with long-range correlations that cannot be explained by the precession of free spins.To obtain these results, we found it essential to use a pulse-level implementation of the Trotterized Ising interaction relying on amplitude-and duration-scaled crossresonance pulses and to apply a number of error mitigation techniques.
These results provide a physics-based benchmark of QPU performance in the NISQ era.Thus, it would be interesting to repeat these experiments on systems with higher quantum volume (QV ), which is a metric that takes into account the number of qubits as well as the error rate of the device [77].Such devices include ibmq washington (127 qubits, QV = 64) or ibmq kolkata (27 qubits, QV = 128), which was used in Ref. [19].In contrast, the IBM QPUs used for our 12-and 19-qubit simulations, ibmq guadalupe and ibmq toronto, have QV = 32.Furthermore, given the variety of available tools for error mitigation, it will be important to undertake a systematic exploration of the optimal implementations of techniques like ZNE and Pauli twirling for non-Clifford gates defined at the pulse level.As part of this effort, it would also be worthwhile to consider alternatives to ZNE including probabilistic error cancellation [21,[78][79][80][81], virtual distillation [82][83][84], or Clifford data regression [24,85].
As devices with higher QV become available, it will be interesting to pursue further the calculation of nontrivial multi-time correlation functions on quantum devices.One quantity that can be computed using the methods proposed in Ref. [57] and further developed here is the transport of conserved quantities.For example, in the MFIM (1.2), the most natural conserved quantity to consider is the energy.The two-point correlation function of the energy density can be used to measure the timescales associated with energy transport.In the chaotic regime of the model, such transport is expected to be diffusive.
Higher-order correlation functions associated with transport beyond the linear-response regime can also be considered, as well as out-of-time-ordered correlation functions (OTOCs) which can be used to characterize quantum chaos [86][87][88][89][90].
Appendix A: Details on Simulation Methods

Scaled-RZX Implementation of the RZZ Gate
In this section, we discuss how to scale the amplitude and duration of the cross resonance (CR) pulses used to realize the R ZX (π/2) to obtain the more general gate R ZX (θ) with arbitrary rotation angle.The R ZX (π/2) gates used to generate the CNOT gate are realized by echoed cross-resonance pulses composed of CR(±π/4) sandwiching an X-echoed π pulse on the control qubit to eliminate the ZI and IX terms in the cross-resonance Hamiltonian [60].Rotary pulses are applied to the target qubit to suppress the interaction IY and ZZ terms [92].The CR pulse has a square-Gaussian shape consisting of a square pulse with width W ( π 2 ) whose boundaries are smoothed into Gaussians with standard deviation σ.We denote the amplitude of the pulse by A( π 2 ).The total area enclosed by the square-Gaussian pulse is therefore where n σ is the number of standard deviations of the Gaussian tails that are chosen to be included in the pulse shape.An example of a pulse schedule for the two-CNOT implementation of the R ZZ gate on ibmq casablanca is shown in Fig. 8(a).This pulse schedule contains one copy of the echoed CR pulse schedule described above for each CNOT gate.
To change the angle of rotation of the R ZX gate from π/2 to an arbitrary θ, we follow the method outlined in Ref. [18] (see also Ref. [61]).When θ > Pulses labeled "CR(π/4)" and "CR(−π/4)" are the cross-resonance pulses described in the text.Other Gaussian pulses correspond to single-qubit gates.On the D1 channel, the overlapping symbols above two circular arrows read V Z(π).This denotes that a "virtual" rotation by π around Z [91] has been implemented using a pulse delay.The symbols below the two central Gaussian pulses read X(π/2) (dark pulse) and Y (π/2) (light pulse), respectively, indicating rotations by π/2 around the X and Y axes, respectively.On the D3 channel, the overlapping symbols below the two central Gaussian pulses read Y (−π/2) (light pulse) and X(π/2) (dark pulse), respectively.(b) Pulse schedule for the scaled-RZX implementation of RZZ (θ = 2.0) on ibmq casablanca.Pulses labeled "Gaus-sianSquare" are the scaled cross-resonance pulses discussed in the text.The scaled-RZX implementation uses half as many cross-resonance pulses as the two-CNOT implementation.reduce the amplitude to Therefore, just like in the circuit depicted in Fig. 2(b), we can sandwich the R ZX (θ) with R Y (π/2) and R Y (−π/2) pulses on the second channel to create the R ZZ (θ) gate pulse schedule shown in Fig. 8(b).As discussed in the main text, this pulse schedule can have a substantially shorter duration than the pulse schedule that implements R ZZ (θ) using two CNOT gates.

Trotter Evolution
In Fig. 9, we compare exact diagonalization (ED) results for the dynamics of the staggered magnetization under the Hamiltonian (1.2) with the results of a Trotter simulation at L = 12 and ∆t = 1.Both simulations use model parameters V = 1 and Ω = 0.24, which are appropriate for the QMBS regime.Significant differences between the ED and Trotter curves are visible starting around V t = 5.However, note that the Trotter circuit dynamics retains the coherent oscillations visible in the ED result.Thus, despite the large time step and the appreciable Trotter error, the Trotter circuit still exhibits strong signatures of QMBS.
We also plot for reference the dynamics of ⟨Z π ⟩ calculated with respect to the Fibonacci-projected ED and Trotter dynamics.The projection has much less effect on the ED dynamics than on the Trotter dynamics.This indicates that postselection of measurement outcomes (as described below) has a larger effect when a larger Trotter step size is used.

Readout Error Mitigation
To mitigate the effect of readout errors on QPU results, one can use readout error mitigation methods as described in Ref. [65] and built into Qiskit Ignis.Suppose that C ideal is a vector containing the list of measurement counts for each computational basis state in the absence of readout error, and that C noisy is the same quantity with readout error.The relationship between C ideal and C noisy can be characterized by a readout error matrix M defined as To obtain the ideal result from the noisy result, one can invert the readout error matrix: Qiskit Ignis supports several methods to obtain the readout error matrix M .One is to approximate M as the tensor product of readout error matrices for each qubit as follows: where ϵ j and η j are the readout error rates for 0 → 1 and 1 → 0 respectively.This method is attractive because the ϵ j and η j can be estimated by executing two circuits to prepare the states |0 . . .0⟩ and |1 . . .1⟩ and then measuring all qubits in the CB.Another method that we call "complete readout error mitigation" prepares and measure all 2 L CB states from |0 . . .0⟩ to |1 . . .1⟩.This allows the elementwise extraction of M for all computational basis states, but is more costly to perform as it involves measuring exponentially many CB states.In this paper, we apply complete readout error mitigation for smaller systems (L = 5) and use the tensor-product approximation for larger systems (L = 12, 19).

Postselection
In the QMBS regime of the Hamiltonian (1.2), the probability of the system being in a CB state with two consecutive 1s is heavily suppressed by the strong nearest-neighbor interaction, as discussed in the main text.Fig. 10  Since the scarred dynamics we are trying to model take place primarily within the Fibonacci Hilbert space [27,29], we take this as evidence that we can amplify their signatures by calculating observables with measurement outcomes postselected to lie within this space.To calculate the postselected dynamics, we remove CB measurement results containing two or more consecutive 1s from the counting dictionary and calculate the expectation value using the rest of the data.
The effect of postselection on our QPU results is shown in Fig. 11, which plots the same dataset as in Fig. 3 with postselection applied.These results already show an enhancement of the oscillatory signal relative to Fig. 3, even with no error mitigation measures besides postselection applied.The results of Fig. 4 clearly demonstrate the utility of applying further error mitigation techniques beyond postselection.

Zero Noise Extrapolation and Pauli Twirling
In order to reduce the effect of gate noise on the calculation of, e.g., an expectation value on a QPU, one can measure this expectation value at different noise scales and extrapolate to the zero-noise limit to estimate the ideal expectation value [6,[21][22][23].To increase the noise scale we can use unitary folding, which acts on a gate G as Here G † = G −1 so the above operation increases the depth of the circuit without changing its logical action on the input state in the unrealistic case where G is implemented noiselessly on the QPU.In the realistic case where G is noisy, this folding operation increases the effect of noise for that gate by a factor of ∼ 3.In this paper, we use the Mitiq package [72] to randomly fold the gates in our circuit to achieve noise scale factors λ ∈ {1.0, 1.5, 2.0} for the full circuit.Since Mitiq does not support folding of the R ZZ (θ) gate, we first generate folded circuits using CNOT gates as placeholders for R ZZ (θ) gates, and then replace all CNOT gates with R ZZ (θ) and R ZZ (−θ) gates as appropriate.Moreover, we also apply the Pauli twirling technique which converts a two qubit gate G's errors into a stochastic form characterized by the effective noise superoperator [6] where F G is the fidelity, σ α c and σ β t (α, β = 0, 1, 2, 3 correspond to Pauli matrices 1, σ x , σ y , σ z , respectively) are Pauli operators acting on a control qubit c and a target qubit t, and ϵ α,β are error probabilities.The quantity [σ α c σ β t ] is a superoperator that acts on a quantum state with density matrix ρ as [σ α c σ β t ]ρ = σ α c σ β t ρσ α c σ β t .To convert the error into this stochastic form, one sandwiches the two-qubit gate G with Pauli gates σ α c σ β t and σ γ c σ δ t with α, β = 0, 1, 2, 3 and γ, δ chosen such that σ γ c σ δ t = G † σ α c σ β t G.This ensures that the sandwiched operator σ α c σ β t Gσ γ c σ δ t = G.Note that choosing γ and δ in this way is possible only if G preserves the Pauli group, i.e., if G is a Clifford gate.After sandwiching, one averages the true noise superoperator over the twoqubit Pauli group [i.e., over all 16 possible pairs (α, β)] to obtain Eq. (A8).In practice, it is sufficient to generate circuits with randomly chosen (α, β) for each two-qubit gate, compute the output of each circuit, and average the result over as many randomly generated circuits as possible.
In this paper, we use two kinds of two-qubit gates, namely CNOT and R ZZ (θ).For the CNOT gate, we randomly select α, β = 0, 1, 2, 3 and choose (γ, δ) so that the random Pauli gates do not affect the CNOT operation, as described in Fig. 12(a).As discussed above, choosing γ and δ in this way is only possible because CNOT is a Clifford gate.However, the R ZZ (θ) gate is a non-Clifford gate for generic values of θ ̸ = 0, π.To solve this problem, we divide the set of Pauli index pairs {(α, β)} into two sets: S C , containing index pairs corresponding to two-qubit Pauli strings that commute with σ 3 c σ 3 t , and S A , containing index pairs corresponding twoqubit Pauli strings that anticommute with σ 3 c σ 3 t .If the randomly selected pair (α, β) ∈ S C , we replace R ZZ (θ) by σ α c σ β t R ZZ (θ)σ α c σ β t .If the pair (α, β) ∈ S A , we replace R ZZ (θ) by σ α c σ β t R ZZ (−θ)σ α c σ β t .The circuits resulting from this procedure are shown in Fig. 12.Note that this modified Pauli twirling procedure assumes that R ZZ (θ) and R ZZ (−θ) have the same gate error channel, which we denote by the noise superoperator N R ZZ .
To verify that the non-Clifford Pauli twirling procedure described above still results in a stochastic error channel for R ZZ (θ), we assume that the action of the noisy R ZZ (θ) gate can be expressed in superoperator form as N R ZZ U R ZZ .Here, U R ZZ is a superoperator that acts as U R ZZ ρ = R ZZ (θ)ρR ZZ (−θ), and the noise superoperator N R ZZ is expressed in Kraus form as N R ZZ ρ = h E h ρE † h with Kraus operators E h = α,β a h;α,β σ α c σ β t satisfying h E h E † h = 1.The action of the twirled noisy R ZZ (θ) gate on a state ρ can then be written as and NR ZZ is the effective noise of R ZZ (θ) after twirling.Since two-qubit Pauli strings with (α, β) ∈ S C commute with R ZZ (θ), the first term can be written as Since two-qubit Pauli strings with (α, β) ∈ S A anticommute with R zz (θ), the second term becomes Summing up two terms above, we find that the effective noise superoperator of R zz (θ) is given by NR ZZ = 1 16 Using σ α σ β σ α = [2δ α,β − (2δ α,0 − 1)(2δ β,0 − 1)]σ β , one can show that this effective noise superoperator takes the stochastic form (A8).This modified Pauli twirling procedure can also be applied to other two-qubit unitary gates generated by Pauli strings.

Dynamical Decoupling
In a quantum computer, physical two-qubit gates have different execution times.When we stack one-and twoqubit gates into a Trotter circuit, there are some idle qubits suffering from thermal relaxation and white noise dephasing.To reduce the decoherence, one can apply appropriate pulse sequences to stabilize the idle qubits during this waiting period.Here, we utilize the pulse sequence τ iq /4 − X π − τ iq /2 − X −π − τ iq /4 with ±π pulse X ±π = R X (±π) and delay time τ iq = (T idle − 2t x,π ).Here, T idle is the idle time of the qubit and t x,π is the duration of the X ±π pulse [15][16][17].

Appendix B: Site-dependent local magnetization results
In this Appendix we show error-mitigated results for the site-wise local magnetization ⟨Z i (t)⟩ as a function of scaled simulation time V t for a 12-site chain in Fig. 13 and a 19-site chain in Fig. 14.The calculations are carried out on ibmq guadalupe and ibm toronto, respectively.We have use the full set of error mitigation techniques described in the main text and Appendix A. Consistent with other observables discussed in the main text, the local magnetizations measured with the scaled-R ZX implementation on QPU are generally in better agreement with the ideal Trotter simulations than the results from the two-CNOT implementation.As an example, in Fig. 13(e), the oscillatory behavior of ⟨Z 5 (t)⟩ for the 12site model in the second oscillation cycle is still visible with the scaled-R ZX implementation, but is completely washed out by noise with the two-CNOT implementation on ibmq guadalupe.The accuracy of local magnetization measurement also shows clear site-dependence, tied to the heterogeneity of qubit quality and native gate fidelity.For example, with the scaled-R ZX implementation in the 12-site model, an oscillation for two cycles can be clearly observed for ⟨Z 2 (t)⟩, while ⟨Z 6 (t)⟩ shows only a weaker first period of oscillation, as shown in Fig. 13(b) and (f).The site-dependence of local magnetization accuracy becomes more evident for the 19-site model calculations on ibmq toronto.For instance, while the oscillation of ⟨Z 3 (t)⟩ is still well reproduced over two cycles, ⟨Z 13 (t)⟩ is almost entirely dominated by noise as shown in Fig. 14(c) and (m).

FIG. 1 .
FIG. 1. Schematic of the state-preparation circuit and the Trotter circuit for one time step ∆t in an L = 5 site chain.The first part of the circuit prepares the Néel state |01010⟩ from the polarized state |00000⟩ using X gates.The Trotter circuit uses single-qubit rotations RX (θ X i ) = e −iθ X

FIG. 2 .
FIG. 2. Benchmarking the two-CNOT and scaled RZX implementations of the RZZ (θ) gate.All data were taken on the IBM QPU Casablanca (ibmq casablanca).(a) Standard implementation of the RZZ (θ) gate using two CNOTs.(b) Implementation of the RZZ (θ) gate using an RZX (θ) gate and appropriate single-qubit rotations.(c) Duration in ns of the pulse schedule on ibmq casablanca as a function of θ = 2V ∆t ≤ 2.5 for the two-CNOT and scaled-RZX implementations of RZZ (θ).The scaled-RZX implementation has a substantially shorter pulse duration for all θ considered.(d) Fidelity of the two-CNOT and scaled-RZX implementations of RZZ (θ) obtained using quantum process tomography on ibmq casablanca as described in the text.The fidelity of the scaled-RZX implementation decreases with increasing θ while that of the two-CNOT implementation remains roughly constant.Data were accumulated over two days, so fluctuations in the fidelity due to calibration drifts are visible.

FIG. 3 .
FIG. 3. Unmitigated Trotter simulation of the staggered magnetization density ⟨Zπ(t)⟩ /L [(a),(b); see Eq. (1.4)] and accumulated error D(t) [(c),(d); see Eq. (3.1)] from the initial state |Z2⟩.Data for chains of 12 [(a),(c)] and 19 [(b),(d)] qubits were obtained using ibmq guadalupe and ibmq toronto, respectively.Data for both the two-CNOT (red) and scaled-RZX (blue) implementations of the Trotter circuit are shown, with noiseless Trotter simulation results (black) for reference.No error mitigation was used here to directly compare the two RZZ implementations.Error bars for each point in (a) and (b) represent the standard deviation of the data over 20 trials.Error bars in (c) and (d) are calculated from those in (a) and (b) by propagation of errors.Oscillations of ⟨Zπ(t)⟩ /L over roughly one period are observed for the scaled-RZX implementation and are barely discernible for the two-CNOT implementation.The two-CNOT implementation also accumulates more error than the scaled-RZX implementation.

FIG. 4 .
FIG. 4. Error-mitigated Trotter simulation of the staggered magnetization density ⟨Zπ(t)⟩ /L [(a),(b); see Eq. (1.4)] and accumulated error D(t) [(c),(d); see Eq. (3.1)] from the initial state |Z2⟩.Data for chains of 12 [(a),(c)] and 19 [(b),(d)] qubits obtained using ibmq guadalupe and ibmq toronto, respectively.Data for both the two-CNOT (red) and scaled-RZX (blue) implementations of the RZZ gate are shown, with noiseless Fibonacci-projected Trotter simulation results (black) for reference.Each data point in (a) and (b) is the result of linear ZNE with scale factors λ ∈ {1.0, 1.5, 2.0} for 10 random Pauli-twirling circuits.Error bars for all data points represent uncertainty in the ZNE and are calculated as described in the main text.The error mitigation strategies outlined in Sec.III result in a substantial improvement for both implementations of the RZZ gate relative to the results shown in Fig. 3.

FIG. 6 .
FIG. 6. Circuits used to calculate the correlation function CY (t) on the QPU.(a) Circuit to calculate ⟨(P Y P )i⟩ M Y i =±1 [Eq.(4.3b)].(b) Circuit to calculate ⟨(P Y P )i⟩ ±Y i [Eq.(4.3c)].Here, UTrotter denotes the Trotter circuit.In the beginning of each circuit, a different initial state is prepared.At the end of each circuit, the expectation value of (P Y P )j is measured.

FIG. 7 .
FIG. 7. Dynamics of the correlator CY (t) [see Eq. (1.6)] from the initial state |Z2⟩ in (a,c) the QMBS regime (V = ∆t = 1, Ω = 0.24) and (b,d) the chaotic regime (V = 1, Ω = 2, ∆t = 0.16).Panels (a,b) are calculated for a chain of 5 qubits using ibmq casablanca, and panels (c,d) are calculated for a chain of 12 qubits using ibmq guadalupe.The calculation uses ZNE for scale factors λ ∈ {1.0, 1.5, 2.0} with 8 (a,b), 10 (c), and 9 (d) random circuit instances for Pauli twirling.Error bars representing the uncertainty in the ZNE were calculated as in Fig. 4. For L = 5, oscillations with relatively slowly decaying amplitude are clearly visible throughout the simulation time window in the QMBS regime (a).For L = 12, these oscillations remain coherent but exhibit a more rapid decay due to the accumulation of gate and readout errors.In the chaotic regime (b,d), the correlator rapidly decays after a single approximate revival and exhibits good agreement with the ideal Trotter simulation results over the full simulation time window.

FIG. 8 .
FIG.8.Pulse schedules for the two implementations of the RZZ gate.(a) Pulse schedule for the two-CNOT implementation of RZZ (θ = 2.0) on ibmq casablanca.Pulse durations are measured in units of dt = 0.2222 ns.Pulses labeled "CR(π/4)" and "CR(−π/4)" are the cross-resonance pulses described in the text.Other Gaussian pulses correspond to single-qubit gates.On the D1 channel, the overlapping symbols above two circular arrows read V Z(π).This denotes that a "virtual" rotation by π around Z[91] has been implemented using a pulse delay.The symbols below the two central Gaussian pulses read X(π/2) (dark pulse) and Y (π/2) (light pulse), respectively, indicating rotations by π/2 around the X and Y axes, respectively.On the D3 channel, the overlapping symbols below the two central Gaussian pulses read Y (−π/2) (light pulse) and X(π/2) (dark pulse), respectively.(b) Pulse schedule for the scaled-RZX implementation of RZZ (θ = 2.0) on ibmq casablanca.Pulses labeled "Gaus-sianSquare" are the scaled cross-resonance pulses discussed in the text.The scaled-RZX implementation uses half as many cross-resonance pulses as the two-CNOT implementation.
σ ), one can change the area under the square-Gaussian pulse to α(θ) = θ π/2 α by adjusting the width of the square pulse as follows:

FIG. 9 .
FIG. 9. Comparison of the dynamics of the staggered magnetization density ⟨Zπ(t)⟩ /L in the QMBS regime (V = 1, Ω = 0.24) obtained from exact diagonalization (ED) (light green) and ideal Trotter simulations (grey) with ∆t = 1.Although the large time step incurs substantial Trotter error, clear long-lived oscillations of the staggered magnetization are visible for both evolutions.The green and black curves indicate Fibonacci-projected ED and ideal Trotter dynamics, respectively.
(a) shows the dynamics of ⟨ L−1 i=1 n i n i+1 ⟩ starting from the |Z 2 ⟩ state under the Hamiltonian (1.2) with Ω = 0.24V for L = 12 and 19 (dark and light green curves, respectively).This indicates that the Hamiltonian dynamics generated by Eq. (1.2) produces a negligible number of pairs of consecutive 1s over the times we

FIG. 10 .FIG. 11 .
FIG. 10.(a) The dynamics of the number of nearest-neighbor pairs of Rydberg excitations, L−1 i=1 nini+1, obtained from ED and ideal Trotter with ∆t = 1 at L = 12 and L = 19.The number of such nearest-neighbor pairs is subextensive even for the Trotter dynamics.(b) The weight of the time-evolved state |ψ(t)⟩ in the Fibonacci Hilbert space, ⟨ψ(t)|P fib |ψ(t)⟩, plotted using ED and ideal Trotter at L = 12 and L = 19.Even though the Trotter dynamics loses more weight in the subspace over the simulation window than the ED dynamics, the Trotter-evolved state retains a finite weight in the subspace at these system sizes.Taken together, these results justify our use of the postselection technique described in this Appendix.

FIG. 13 .FIG. 14 .
FIG.13.Complete list of error-mitigated local magnetization results ⟨Zi⟩ versus time V t for a 12-site chain measured on ibmq guadalupe using the scaled-RZX and two-CNOT implementations (blue and red, respectively).The ideal Trotter simulation data (black) are also shown for reference.